Second-Order Approximation of the Seismic Reflection Coefficient in Thin Interbeds

As most of the lithostratigraphic reservoirs in China are thin interbeds, the study of seismic responses in thin interbeds is an integral part of lithologic reservoir exploration. However, at present, the research on seismic reflection coefficients of thin interbeds in exploration seismology is still weak, which leads to the lack of theoretical basis for the subsequent interpretation of amplitude variation with offset (AVO) related to thin interbed. To solve this problem, in this paper, we proposed second-order approximate equations of the seismic reflection coefficients in thin-bed and thin-interbed layers. Under the assumption of a small impedance contrast in layered media, we made a second-order approximation with a more evident physical meaning to the reflection coefficient calculation method proposed by Kennett. Then, based on the test of the single thin-layer theoretical model, it was confirmed that the second-order approximation equation of the PP-wave (reflected compressional wave) is accurate at incident angles less than 30°, and that of the PS-wave (converted shear wave) is accurate at wider incident angles. Finally, based on the single-thin-bed equations, the approximate equations of seismic reflection coefficients in thin interbeds were established, the validity of which was verified by the theoretical model. Our equations will be applicable to the calculation of PP- and PS-wave reflection coefficients in thin interbeds where internal multiples are difficult to suppress and transmission loss is hard to accurately compensate. This lays a theoretical foundation for improving the seismic prediction accuracy of lithologic reservoirs.


Introduction
The technology of amplitude versus offset (AVO) or amplitude versus incident angle (AVA) has been used to predict rock properties and fluid in the subsurface reservoirs for decades [1][2][3][4]. Conventional AVO analysis is based on Zoeppritz equations or their approximations; however, these equations are not suitable for thin-bed problems due to the mixing of the reflected waves, converted waves, and internal multiples, which are hard to separate [5][6][7].
The problem of thin-bed layers has been studied for decades. Typically, a single layer with a thickness of less than a quarter of a wavelength is defined as a thin layer in seismology [8,9]. Brekhovskikh [10,11] gave the displacement potential reflections and transmissions of plane waves propagating in layered elastic media. Kennett [12] proposed the recursive algorithm for the calculation of total reflection and transmission coefficients for a stack of layers, in which the unconditional stability for all frequencies and slownesses was improved. The recursive algorithm solved the overflow problems in the calculation of exponential functions for high frequencies and slownesses [13,14]. Liu and Schmitt [15] proposed an exact analytical solution in the acoustic wave domain and simulated the reflection amplitude and AVO responses of a single thin bed for arbitrary incident angles. Pan and Innanen [6] simplified the elastic multi-layered medium equation based on Breshkovskikh's equations to analyze the relationship between the AVO effect and frequency in thin elastic media. Yang et al. [16] derived thin-bed reflection and transmission coefficients defined in terms of displacements and approximated them to be in a quasi-Zoeppritz matrix form under the assumption of a very thin thickness for the thin bed. Yang et al. [17] proposed a low-order series approximation of thin bed for PP-wave, and estimated the P-wave (compressional wave) impedance ratios and thickness. Lu et al. [18] equated the set of thin interlayers with vertical transverse isotropic stratigraphy and used the reflection coefficients of anisotropic media for AVO inversion study.
In general, the existing approximate equations of thin-bed reflection coefficients are based on the assumption of small incident angle and weak impedance contrast. Additionally, these equations are complicated and their physical meaning is not clear, so they do not easily provide theoretical support for the subsequent AVO interpretation of single thin beds or thin interbeds. In this paper, we performed a numerical simulation analysis on the single thin-bed model based on the Kennett equation and found that the high-order multiples contribute less to the overall reflection coefficient of the single thin bed. Thus, we deduced the second-order approximate equation with a more evident physical meaning of the reflection coefficient in a single thin-bed, which is suitable for both PP-and PS-waves. Further, we used a recursive algorithm [12] to extend the second-order approximation equation of the thin-bed reflection coefficient to the thin-interbed model and verified its accuracy by numerical simulation analysis.

Exact Reflection Coefficient Equation in Stratified Media
Our single thin-bed reflection coefficient equation is deduced based on the exact reflection coefficient equation in stratified media proposed by Kennett [13,14]. For the stratified strata, as shown in Figure 1, A and C are the top and bottom interfaces, respectively, and B is an interface inside the model. The down-going wave (denoted by subscript 'D') propagates in area AC (the area between interfaces A and C), generating reflected waves, converted waves, and internal multiples. These wavefields are mixed together, and the up-going waves are denoted by subscript 'U'. Kennett [13,14] derived the following equation using the recursive algorithm to calculate the overall reflection and transmission coefficient matrices of region AC: where r A D and t A D are the single-interface reflection coefficient and transmission coefficient matrices at interface A, respectively: where the superscripts 'PP', 'PS', 'SS', and 'SP' represent the four wave modes (P-wave reflected/ transmitted to P-wave, P-wave reflected/transmitted to shear wave, shear wave reflected/ transmitted to shear wave, shear wave reflected/transmitted to P-wave,) of single interface reflections or transmissions, which can be calculated by Zoeppritz equations [19].R BC D is obtained from the total reflection matrix R BC D of region BC by phase shift to the bottom of interface A: where E is the phase income for downward propagation through the layer. Furthermore, the following exists: where ω is the angular frequency; h indicates the thickness of region AB; q α and q β are the vertical slownesses for P-and S-waves (shear wave): where α and β are the P-and S-wave velocities, respectively, and p is the horizontal slowness.
2 elastic media. Yang et al. [16] derived thin-bed reflection and transmission coefficients defined in 49 terms of displacements and approximated them to be in a quasi-Zoeppritz matrix form under the 50 assumption of a very thin thickness for the thin bed. Yang et al. [17] proposed a low-order series 51 approximation of thin bed for PP-wave, and estimated the P-wave (compressional wave) impedance 52 ratios and thickness. Lu et al. [18] equated the set of thin interlayers with vertical transverse isotropic 53 stratigraphy and used the reflection coefficients of anisotropic media for AVO inversion study. 54 In general, the existing approximate equations of thin-bed reflection coefficients are based on 55 the assumption of small incident angle and weak impedance contrast. Additionally, these equations 56 are complicated and their physical meaning is not clear, so they do not easily provide theoretical 57 support for the subsequent AVO interpretation of single thin beds or thin interbeds. In this paper, we 58 performed a numerical simulation analysis on the single thin-bed model based on the Kennett 59 equation and found that the high-order multiples contribute less to the overall reflection coefficient 60 of the single thin bed. Thus, we deduced the second-order approximate equation with a more evident 61 physical meaning of the reflection coefficient in a single thin-bed, which is suitable for both PP-and 62 PS-waves. Further, we used a recursive algorithm [12] to extend the second-order approximation 63 equation of the thin-bed reflection coefficient to the thin-interbed model and verified its accuracy by 64 numerical simulation analysis.

68
Our single thin-bed reflection coefficient equation is deduced based on the exact reflection 69 coefficient equation in stratified media proposed by Kennett [13,14]. For the stratified strata, as shown 70 in Figure 1, A and C are the top and bottom interfaces, respectively, and B is an interface inside the 71 model. The down-going wave (denoted by subscript 'D') propagates in area AC (the area between 72 interfaces A and C), generating reflected waves, converted waves, and internal multiples. These 73 wavefields are mixed together, and the up-going waves are denoted by subscript 'U'. Kennett [13,14] 74 derived the following equation using the recursive algorithm to calculate the overall reflection and 75 transmission coefficient matrices of region AC： 76

Assumption of Small to Medium Impedance Contrast
We extracted a single thin bed model from the layered medium model shown in Figure 1. As shown in Figure 2, Layers 1 and 3 are the adjacent overlying and underlying strata of the single thin bed Layer 2, respectively. According to Equation (1), the total reflection matrix at top interface 1 can be written as whereR (2) D is the phase-shifted reflection matrix relative to interface 2.

R
(2) where E (2) is the phase income for downward propagation through the thin layer with the thickness h (2) , the function of which is to move the reflection coefficient matrix r (2) D to the bottom of interface 1. Using Taylor's expansion approach, Equation 7 can be expanded as where the total reflection matrix R (1) D is in 2 × 2 form of 3 slownesses for P-and S-waves (shear wave): 87 where α and β are the P-and S-wave velocities, respectively, and p is the horizontal slowness.

94
We extracted a single thin bed model from the layered medium model shown in Figure 1. As 95 shown in Figure 2, Layers 1 and 3 are the adjacent overlying and underlying strata of the single thin 96 bed Layer 2, respectively. According to Equation 1, the total reflection matrix at top interface 1 can 97 be written as 98 Since the current multi-component seismic exploration mainly uses reflected PP-and PS-waves, in this paper, we will only study the reflection coefficients of the reflected PP-and PS-waves. According to the different propagation paths, the PP-and PS-waves relating to R  where (2) D R is the phase-shifted reflection matrix relative to interface 2. (2) (2) (2) (2) D D , = R E r E (8) where E (2) is the phase income for downward propagation through the thin layer with the 100 thickness h (2) , 101

99
the function of which is to move the reflection coefficient matrix r to the bottom of interface 1.

102
Using Taylor's expansion approach, Equation 7 can be expanded as 103 where the total reflection matrix Since the current multi-component seismic exploration mainly uses reflected PP-and PS-waves, 105 in this paper, we will only study the reflection coefficients of the reflected PP-and PS-waves. 106 According to the different propagation paths, the PP-and PS-waves relating to

120
We grade the energy of the internal multiples and retain the multiples that contribute most to 121 the single thin-bed reflection. 122 Let the absolute amplitude of the incident wave be A0, and the absolute amplitudes of the first-123 order multiples and second-order multiples in the single thin layer ( Layer 2 ) be A1 and A2, 124 respectively. According to the definition of reflection and transmission coefficients [20], the 125 relationships between A1, A2, and A0 are as follows: 126 Similarly, the absolute amplitude An of the n th -order multiples is 127 Where 128 (1) (2) U D r r .
If n = 0, Equation 14 represents the absolute amplitude of the primary reflected wave. 129 Because the seismic wave energy E is proportional to the square of amplitude and to the primary 130 wave energy EP, we obtain 131 Thus, for the energy of all internal multiples EM, we obtain 132 Because the range of the reflection coefficients  We grade the energy of the internal multiples and retain the multiples that contribute most to the single thin-bed reflection.
Let the absolute amplitude of the incident wave be A 0 , and the absolute amplitudes of the first-order multiples and second-order multiples in the single thin layer ( Layer 2 ) be A 1 and A 2 , respectively. According to the definition of reflection and transmission coefficients [20], the relationships between A 1 , A 2 , and A 0 are as follows: Similarly, the absolute amplitude A n of the n th -order multiples is where If n = 0, Equation (14) represents the absolute amplitude of the primary reflected wave. Because the seismic wave energy E is proportional to the square of amplitude and to the primary wave energy E P , we obtain Thus, for the energy of all internal multiples E M , we obtain Because the range of the reflection coefficients r Energies 2020, 13, 1465 Then, the ratio between the energy of all internal multiples E M and the energy of the primary wave E P is The ratio between the energy of first-order multiples E M1 and the energy of all multiples E M is The ratio of the sum energy of the first-order and second-order multiples to the energy of all multiples E M is Furthermore, the ratio of the sum energy of the first-order, second-order, and third-order multiples to the energy of all multiples E M is In the derivation of Equations (12)- (22), wave types are not distinguished, so they are applicable to both PP-and PS-waves. Let the reflection coefficient be in the range of (−1, 1); the calculated energy ratios are shown in Figure 5, from which we can conclude the following.
(1) The smaller the absolute value of the reflection coefficient, the smaller the ratio of E M /E P . Figure 5a shows that, as the reflection coefficient decreases, E M /E P also decreases. When the absolute value of the reflection coefficient is less than 0.3, the ratio of E M /E P will be less than 0.1, and the energy ratio corresponds to the blue color in Figure 5a.
(2) The sum of the first-and second-order multiple energies occupies the major part of the internal multiples energy. Figure 5b shows that when the absolute value of the reflection coefficient is less than about 0.3, the energy of the first-order multiples occupies more than 95% of the multiple energy. Besides, it can be seen from Figure 5c that when the absolute value of the reflection coefficient is less than about 0.6, the sum energy of the first-and second-order multiples occupy more than 95% of the multiple energy. In addition, Figure 5d shows that when the absolute values of the reflection coefficients are less than about 0.75, the sum of the energy of the first-, second-, and third-order multiples occupies more than 95% of the multiple energy. In actual seismic data, the absolute values of the reflection coefficients of most strata are less than 0.6. Therefore, for most single thin layers, the contribution of multiple waves in the overall reflection coefficient only needs to take the first-and second-order multiples into account.
According to the energy ratio analysis, if the wave impedance contrast is very small-that is, the reflection coefficient is also extremely small-the energy of internal multiples will be very weak. In this case, internal multiples may not be considered, and then Equation (10) can be approximated as Under the condition of moderate impedance contrast (i.e., the absolute value of reflection coefficient is less than 0.6), according to the energy ratio analysis, let n = 2 in Equation (10), and the second-order approximation equation of the total reflection coefficient of the single thin-bed is the expansion of which is According to the energy ratio analysis, if the wave impedance contrast is very small-that is, the 160 reflection coefficient is also extremely small-the energy of internal multiples will be very weak. In 161 this case, internal multiples may not be considered, and then Equation 10 can be approximated as 162 Under the condition of moderate impedance contrast (i.e., the absolute value of reflection 163 coefficient is less than 0.6), according to the energy ratio analysis, let n = 2 in Equation 10, and the 164 second-order approximation equation of the total reflection coefficient of the single thin-bed is 165 the expansion of which is 166 As shown in Figure 6, in Equation 25, each item has a clear physical meaning. The details are as 167 follows: 168 (1) (1) r : the reflection coefficient at the top interface (Interface 1) of the thin bed.
As shown in Figure 6, in Equation (25), each item has a clear physical meaning. The details are as follows: (1) r (1) D : the reflection coefficient at the top interface (Interface 1) of the thin bed.
D : the first-order multiple reflection coefficient of the thin-bed. An additional interaction in the thin-bed is introduced based on the primary reflection.
D : the second-order multiple reflection coefficient of the thin-bed. An additional interaction in the thin-bed is added based on the first-order multiple reflections.

179
It should be noted that the reflection coefficient calculated by Equation 24 is in the slowness-180 frequency domain, and can be written as 181 However, the reflection coefficient used in seismic inversion is mainly in the angle-time domain, 182 and thus a domain conversion is necessary. The entire domain transformation process can be divided 183 into the following two steps.

184
(1) Slowness to incident angle 185 According to Snell's law, there is 186 where θ is the incident angle. Given the incident angle and the P-wave velocity of the overlying 187 formation of the thin layer, the corresponding horizontal slowness can be calculated according to 188 Using an inverse Fourier transform, we can get the reflection coefficient in the intercept time and 191 angle (τ-θ) domain [14,21]: 192 Then, the AVA synthetic seismogram can be obtained by the convolution of the reflection 193 coefficients and wavelets: 194 where W denotes the wavelets. 195

Accuracy Analysis of Second-Order Approximation 196
We used a theoretical model to test the accuracy of our second-order approximate equation  It should be noted that the reflection coefficient calculated by Equation (24) is in the slowness-frequency domain, and can be written as However, the reflection coefficient used in seismic inversion is mainly in the angle-time domain, and thus a domain conversion is necessary. The entire domain transformation process can be divided into the following two steps.
(1) Slowness to incident angle According to Snell's law, there is where θ is the incident angle. Given the incident angle and the P-wave velocity of the overlying formation of the thin layer, the corresponding horizontal slowness can be calculated according to (2) Frequency domain to time domain Using an inverse Fourier transform, we can get the reflection coefficient in the intercept time and angle (τ-θ) domain [14,21]: Then, the AVA synthetic seismogram can be obtained by the convolution of the reflection coefficients and wavelets: where W denotes the wavelets.

Accuracy Analysis of Second-Order Approximation
We used a theoretical model to test the accuracy of our second-order approximate equation (Equation (24)) of the single thin-bed reflection coefficient. In the model, the elastic parameters of the first and third layers are fixed and those in the middle layer change with varying impedance contrast. The parameters of the first and third layers are set as α 1 = 3.094 km/s, β 1 = 1.515 km/s, ρ 1 = 2.4 g/cm 3 , h 1 = 150 m. The thickness of the thin bed is 8 m, and its elastic parameters are shown in Table 1 with different P-wave impedance differences between the thin-bed layer and the top layer. Besides, the S-wave velocities and densities of the thin-bed layer are converted from the P-wave velocity according to Castagna equation [22] and Gardner equation [23], respectively. In the numerical simulation, the Ricker wavelets with a dominant frequency of 40 Hz and 30 Hz are used for PP-and PS-waves, respectively [24]. Taking ∆I = −1.5 as an example, Figures 7 and 8 display the angle gathers of PP-and PS-waves, respectively. It can be seen from Figures 7e and 8e that the difference between the angle gather including all the reverberations and the second-order approximation is very small. Besides, Figures 7d  and 8d show that the second-order approximate angle gather contains multiples, which are aliased together with the primary waves. It is difficult for the existing technology to eliminate these internal multiples without affecting the energy of primary reflections.
To further quantify the accuracy of the second-order approximation equation, we extract the amplitude along with the top interface of the thin bed, where the amplitude of the second-order approximation is A 2 , and the amplitude of the data containing all reverberations is A all . We define the amplitude change ratio A ra as (30) Figure 9 shows the change of A ra with the incident angles of PP-and PS-wave corresponding to the thin-bed model shown in Table 1, from which we can obtain the following rules: (1) With the increase of impedance contrasts, A ra also tends to increase.
(2) If the incident angle exceeds 30 • , the error will rapidly increase when using the second-order approximation equation to calculate the PP-wave reflection coefficient. For the PP-wave, initially, A ra is small and increases very slowly with the increasing angle, but after the angle exceeds 30 • , A ra increases rapidly. (3) The A ra of the PS-wave decreases as the incident angle increases, which shows that the second-order approximation equation is suitable for calculating PS reflections at large angles. However, if the PP-and PS-wave are jointly used for seismic inversion, the applicable angle range of the second-order approximate equation needs to take the intersection of PP-and PS-waves. respectively. It can be seen from Figure 7e and Figure 8e that the difference between the angle gather 210 including all the reverberations and the second-order approximation is very small. Besides, Figures  211  7d and 8d show that the second-order approximate angle gather contains multiples, which are aliased 212 together with the primary waves. It is difficult for the existing technology to eliminate these internal 213 multiples without affecting the energy of primary reflections. 214     approximation is A2, and the amplitude of the data containing all reverberations is Aall. We define the 217 amplitude change ratio Ara as 218  Ara is small and increases very slowly with the increasing angle, but after the angle exceeds 30°, 224 Ara increases rapidly.

225
(3) The Ara of the PS-wave decreases as the incident angle increases, which shows that the second-226 order approximation equation is suitable for calculating PS reflections at large angles. However, 227 if the PP-and PS-wave are jointly used for seismic inversion, the applicable angle range of the 228 second-order approximate equation needs to take the intersection of PP-and PS-waves. In order to test the influence of frequency on the second-order approximate equation, we used waves with different wavelengths (frequencies) to test the single thin-bed model. As shown in Figures 10 and 11, taking ∆I = -1.5 as an example, we calculated the modulus and phase angles of the second-order approximate reflection coefficients of PP-and PS-waves in the frequency domain based on Equations 26 and 27, and compared them with the exact reflection coefficients calculated by the Kennett equation (Equation 7). The wavelengths of the incident waves are set as λ = h, λ = 2h,  Figure 10c,d and Figure 11c,d. The critical angle of the model is about 59.4 • . It can be seen that the error is relatively small when the incident angle is within the critical angle. Especially when the thickness is less than half a wavelength, the relative error is almost zero. If the angle of incidence is greater than the critical angle, the error of the second-order approximation equation is relatively large. shown in Figures 10c, 10d, 11c, and 11d. The critical angle of the model is about 59.4°. It can be seen 245 that the error is relatively small when the incident angle is within the critical angle. Especially when 246 the thickness is less than half a wavelength, the relative error is almost zero. If the angle of incidence 247 is greater than the critical angle, the error of the second-order approximation equation is relatively 248 large. 249 Figure 10. Accuracy analysis of the second-order approximate reflection coefficient of the PPwave in the frequency domain: the modulus (a) and phase angle (b) of the second-order approximate reflection coefficient, the relative modulus error (c), and the relative phase error (d). 12 that the error is relatively small when the incident angle is within the critical angle. Especially when 246 the thickness is less than half a wavelength, the relative error is almost zero. If the angle of incidence 247 is greater than the critical angle, the error of the second-order approximation equation is relatively 248 large. 249 Figure 10. Accuracy analysis of the second-order approximate reflection coefficient of the PPwave in the frequency domain: the modulus (a) and phase angle (b) of the second-order approximate reflection coefficient, the relative modulus error (c), and the relative phase error (d). Figure 11. Accuracy analysis of the second-order approximate reflection coefficient of the PSwave in the frequency domain: the modulus (a) and phase angle (b) of the second-order approximate reflection coefficient, the relative modulus error (c), and the relative phase error (d).

Assumption of Small to Medium Impedance Contrast 251
Since thin-interbed layers are widespread in the Earth, such as interbeds formed by thin sand 252 and thin mud, it is necessary to expand the approximate equation from the thin bed to thin interbed 253 Figure 11. Accuracy analysis of the second-order approximate reflection coefficient of the PS-wave in the frequency domain: the modulus (a) and phase angle (b) of the second-order approximate reflection coefficient, the relative modulus error (c), and the relative phase error (d).

Assumption of Small to Medium Impedance Contrast
Since thin-interbed layers are widespread in the Earth, such as interbeds formed by thin sand and thin mud, it is necessary to expand the approximate equation from the thin bed to thin interbed signal. If several thin layers with small impedance contrasts are superposed together, a set of thin interbeds can be formed. We chose a recursive algorithm [12][13][14] to generate the response of a thin-interbed medium recursively by adding one layer at a time ( Figure 12). Since there is no reflection under the bottom interface, the total reflection coefficient of the bottom R (n+1) D is the same as the reflection coefficient of the single interface r According to Equation (24), the second-order approximate total reflection coefficient at the top of the nth layer is . The second-order approximate total reflection coefficient at any interface K in the thin interbed can be derived by the recursive algorithm:

Accuracy Analysis of Second-Order Approximation 266
Compared with the single thin-bed model, the thin-interbed model will produce more complex 267 internal multiples, which cannot be discriminated by the reflection coefficient curve. Therefore, we 268 will directly use the angle gather to analyze the accuracy of the second-order approximate equation 269 of the thin-interbed model. The thickness of a single thin bed in the thin-interbed model is set as 6 m, 270

Accuracy Analysis of Second-Order Approximation
Compared with the single thin-bed model, the thin-interbed model will produce more complex internal multiples, which cannot be discriminated by the reflection coefficient curve. Therefore, we will directly use the angle gather to analyze the accuracy of the second-order approximate equation of the thin-interbed model. The thickness of a single thin bed in the thin-interbed model is set as 6 m, and the elastic parameters of the thin-interbed model are shown in Table 2. In the numerical simulation, we used the Ricker wavelet with dominant frequencies of 40 Hz and 30 Hz for PP-and PS-waves, respectively. The angle gathers of the PP-and PS-waves of the thin-interbed model are shown in Figures 13 and 14. In order to avoid wide-angle reflection, the incident angle range is set to 0 • -65 • . Figures 13 and 14 reveal that: (1) The simulated angle gather of the second-order approximate equation is very similar to the angle gather including all reverberations (Figures 13e and 14e). (2) The second-order approximation equation effectively calculates the internal multiples (Figures 13d  and 14d), and these multiples are aliased together with the primary wave. It is difficult for existing technology to eliminate these internal multiples without affecting the energy of primary reflections. (3) As can be seen from Figures 13d and 14d, for PP-waves, the internal multiples have a more significant effect on the reflection coefficient at small to medium angles, and the effect decreases as the incident angle increases. For PS-waves, the internal multiples have a large effect on the reflection coefficient at medium to large angles, and the effect becomes more significant as the incident angle increases.

Comparison With Field Multicomponent Data 288
In order to verify the validity of the second-order approximate equation, we simulated the angl 289 gather of field multicomponent data in the Sichuan Basin, West China [24]. Figure 15a and Figure 16  290 are the angle gather of field data at the well location for PP-wave and PS-wave, respectively. Besides 291 in order to compare the simulation accuracy, we simulated the angle gathers using exact Zoepprit 292 equations and second-order approximate equation. The seismic wavelets used for simulation ar 293 statistical wavelets extracted from the corresponding filed data [24]. Moreover, the simulation angl 294 gathers are shown in Figures 15 and 16. In summary, the comparison with field data reveals that: 295 (1) The simulated angle gathers based on the second-order approximate equation will be mor 296 similar to the field data. For example, the field angle gather in the green box in Figure 15a has 297 phase inversion phenomenon, but this phenomenon only appears in the simulated angle gathe 298 based on the second-order approximate equation.

299
(2) For thin interbeds, the simulation results of the second-order approximate equation contain 300 richer information. As shown in Figure 16, the strata in the green box is thin interbeds. Moreover 301 the simulation results based on the second-order approximate equation contain interna 302 multiples and are more similar to the field data. 303

Comparison With Field Multicomponent Data
In order to verify the validity of the second-order approximate equation, we simulated the angle gather of field multicomponent data in the Sichuan Basin, West China [24]. Figures 15a and 16a are the angle gather of field data at the well location for PP-wave and PS-wave, respectively. Besides, in order to compare the simulation accuracy, we simulated the angle gathers using exact Zoeppritz equations and second-order approximate equation. The seismic wavelets used for simulation are statistical wavelets extracted from the corresponding filed data [24]. Moreover, the simulation angle gathers are shown in Figures 15 and 16. In summary, the comparison with field data reveals that: (1) The simulated angle gathers based on the second-order approximate equation will be more similar to the field data. For example, the field angle gather in the green box in Figure 15a has a phase inversion phenomenon, but this phenomenon only appears in the simulated angle gather based on the second-order approximate equation.
(2) For thin interbeds, the simulation results of the second-order approximate equation contain richer information. As shown in Figure 16, the strata in the green box is thin interbeds. Moreover, the simulation results based on the second-order approximate equation contain internal multiples and are more similar to the field data. 16 are the angle gather of field data at the well location for PP-wave and PS-wave, respectively. Besides, 291 in order to compare the simulation accuracy, we simulated the angle gathers using exact Zoeppritz 292 equations and second-order approximate equation. The seismic wavelets used for simulation are 293 statistical wavelets extracted from the corresponding filed data [24]. Moreover, the simulation angle 294 gathers are shown in Figures 15 and 16. In summary, the comparison with field data reveals that: 295 (1) The simulated angle gathers based on the second-order approximate equation will be more 296 similar to the field data. For example, the field angle gather in the green box in Figure 15a has a 297 phase inversion phenomenon, but this phenomenon only appears in the simulated angle gather 298 based on the second-order approximate equation.

299
(2) For thin interbeds, the simulation results of the second-order approximate equation contain 300 richer information. As shown in Figure 16, the strata in the green box is thin interbeds. Moreover, 301 the simulation results based on the second-order approximate equation contain internal 302 multiples and are more similar to the field data.

Discussion
The reflection coefficient equations of the thin-bed and thin-interbed media discussed in this paper are approximated based on the Kennett equation. Our second-order approximate equations have no layer thickness limitation since there is no layer thickness limitation in the Kennett equation. For thick layers, multiples can be suppressed by the predictive deconvolution method. However, for thin-bed and thin-interbed layers, internal multiples are aliased together with primary waves, and it is difficult for the existing technologies to eliminate these internal multiples without affecting the energy of primary reflections. Therefore, for the inversion of a single thin bed and thin interbed, internal multiples must be considered, since this has a significant impact on the overall reflection coefficient.
During the analysis of the single thin-bed model, we found that if the impedance contrast is too large, the error of the approximate equation will increase. Therefore, the approximate equation of the thin-interbed is not suitable for strata with a large impedance contrast. Besides, although there are only two sets of interbed elements in thin-interbed model testing, the approximate equation can be applied to multiple sets of interbed elements.
We did not consider the influence of the incident angle during the derivation of the approximate equation. However, we found that the approximate equation of the PP-wave will have higher accuracy in the range of small and medium angles in the model tests. On the contrary, the approximate equation of the PS-wave can be applied to a broader range of incident angles. However, in order to cooperate with the PP-wave, we limited the application scope of the PS-wave approximation equation to the range of the PP-wave approximation equation. Most of the actual seismic data are acquired within a medium incident angle, because a large incidence angle tends to produce wide-angle reflections. In addition, the distortion caused by normal moveout correction is severe at large incident angles. Thus, reflections at large angles are often removed during data processing. Therefore, the second-order approximate equation is suitable for most seismic data.
The single thin-bed models and thin-interbed models can be divided into multiple types according to the strata parameters and thicknesses. The conclusions and phenomena we discussed above are only for the models involved in this paper. Therefore, there may be different phenomena in other models, which we will test in the future work.
Forward modeling using our equation is implemented in the frequency domain, and then the seismic record can be obtained by the conversion to the time domain. Thus, we should consider the time-frequency conversion when using our equations for AVO inversion, and the time consumption of this process is relatively large. In the next study, we will focus on designing a fast algorithm for our equation to improve the efficiency of AVO inversion.
In the numerical simulation analysis of the thin-interbed model, the thickness of each thin bed is set to a constant, and the elastic parameters of thin beds with the same lithology are also set to be consistent. Therefore, this model is idealized by us, but it does not affect the accuracy analysis of the second-order approximate reflection coefficient equations. The thickness of each thin bed in the actual thin-interbed strata is different, and the elastic parameters of thin beds with the same lithology may vary significantly. The anomalous seismic responses caused by the changes of physical properties of a thin bed in thin interbeds will be further studied in the future.

Conclusions
We analyzed the relationship between the energy of internal multiples and primary waves based on the single thin-bed model. Then, we proposed the second-order approximate reflection coefficient equations with a clear physical meaning for the single thin-bed model. For PP-waves, the approximate equation has a high accuracy at small to medium angles (≤30 • ), while for PS-waves, the equation can be applied to a larger range of angles. Under the assumption of weak impedance contrast, we extended the second-order approximate equation of the single thin-bed model to the thin-interbed model and tested its accuracy through numerical simulation analysis. At present, it is challenging to remove internal multiples in thin interbeds without affecting primary waves. In this case, the second-order approximate reflection coefficient equations can be used to calculate the total reflection coefficients of thin interbeds more conveniently and efficiently, which is conducive to improving the accuracy of seismic inversion.
Author Contributions: J.L. conceived the idea of this research. Z.Y. and J.L. deduced the equations. Z.Y. programmed the codes and performed the simulation tests. The paper was written by all the authors. All authors have read and agreed to the published version of the manuscript.