Data-driven method for an improved linearised AVO inversion

Linearised approximations of the Zoeppritz equation are widely used for amplitude versus offset (AVO) forward modelling and prestack seismic inversion. However, current linearised approximations typically exhibit large errors for strong-contrast media and large incidence angles (even for those less than the critical angle), leading to poor quality estimation of elastic parameters. We therefore develop a data-driven method to improve the linearised AVO approximation, based on the concept that the inaccuracy of the conventional linearised AVO forward operator is solely responsible for generating residuals of the PP reflection coefficients calculated by the conventional linearised approximation compared with using the Zoeppritz equation. Therefore, if true model parameters from well-logging data are known, the residuals of the PP reflection coefficients can be used to estimate linear modifications of the inaccurate linearised AVO forward operators at well sites. We apply estimated linear modifications to original inaccurately linearised forward operators to obtain an improved linear AVO operator (ILAO). Because ILAOs can only be estimated at the well sites by the data-driven method, we construct spatially variant ILAOs using structure-guided interpolation for those at a distance from the well sites. Numerical example results reveal that the improved linear AVO approximation is more accurate than the Aki–Richards approximation for strong-contrast media and large incident angles. Moreover, the accuracy and resolution of the inverted Pand S-wave velocities and density using ILAOs are higher than those using the conventional Aki–Richards operators. Finally, application to the field data successfully demonstrates the feasibility and effectiveness of the improved linearised AVO inversion method.


Introduction
Amplitude versus offset (AVO) prestack seismic inversion, which estimates subsurface elastic parameters from incidence angle gathers, plays an important role in reservoir characterisation. AVO inversion is based on the concept that seismic reflections depend on the incident angles and contrasts of elastic parameters over the reflected interfaces. The fundamental theory of AVO is described by the Zoeppritz equation (Zoeppritz 1919), which presents the energy distribution as functions of elastic parameters and incident angle when a plane P-wave arrives at an interface between two layers of semiinfinite media.
At present, AVO inversion methods are extensively based on linearised approximations of the Zoeppritz equation because of their computational simplicity and low time consumption. Aki & Richards (1980) proposed a linear approximation (hereafter, the Aki-Richards approximation) related to P-and S-wave velocities and density for weak-contrast media and small incident angles (less than 30°). Wiggins et al. (1983) proposed an intercept-gradient-curvature linear approximation dependent on P-and S-wave velocities and density, which is the basis for most AVO analysis. Shuey (1985) then transformed the Wiggins approximation into an equation depending on P-wave velocity, density and Poisson's ratio, which is widely used in the characterisation of gas reservoirs (Castagna & Swan 1997). Fatti et al. (1994) proposed a further approximation related to P-and S-wave impedance and density, while Gray et al. (1999) developed a new approximation expressed by Lamé parameters using fundamental rock properties. Russell et al. (2011) considered the porous fluid effect to propose a linear approximation applied to direct estimation of fluid indicators. All these approximations are reformulations of the Aki-Richards approximation, as Russell (2014) concluded. However, the Aki-Richards approximation is a linear physical relationship that only reflects the first-order approximation of AVO behavior suitable for weak-contrast media and small incidence angles (Buland & Omre 2003;Mallick 2007;Russell et al. 2011;Zong et al. 2013;Mallick & Adhikari 2015). With the development of oil and gas seismic exploration, the targets are moving to complex reservoirs and the large offsets are commonly used in seismic data acquisition. Therefore, it is necessary to develop AVO inversion technology suitable for a large incident angle and strong-contrast media. Many researchers have proposed quasi-linear or nonlinear AVO inversion methods based on the Zoeppritz equation or wave equation. For example, Wang (1999) made quadratic and quartic approximations of the Zoeppritz equation as the function of the ray parameter, and then proposed a ray-impedance concept for AVO inversion (Wang 2003). Chen & Wei (2012) presented a Zoeppritz-based inversion method to estimate the four ratios of elastic parameters in the ray parameter domain. Zhu & McMechan (2012) proposed an AVO inversion method that calculated the Fréchet derivatives of the Zoeppritz equation analytically and then tested the method for a single interface reflection. Zhi et al. (2016) solved the Zoeppritz equation by combining the Levenberg-Marquardt algorithm with the Tikhonov regularisation. Pan et al. (2017) proposed an AVO inversion method based on the Zoeppritz equation by using the Markov chain Monte Carlo process. Although these methods provide better inversion results, they are relatively timeconsuming owing to the complexity and nonlinearity of the inverse problem. Considering this, Mallick (2007Mallick ( , 2015 presented a prestack-waveform inversion method with acceptable time consumption based on the parallelisation of the genetic algorithm. Further, some researchers have provided new ideas for AVO seismic inversion. For example, Causse et al. (2007a,b) proposed a model-based linear AVO approximation with basic functions by applying singular value decomposition to the modeled AVO reference curves. She et al. (2018) presented a data-driven AVO inversion method by pre-learning the dictionaries of elastic parameters from well logs and then applying a sparse representation of subsurface model parameters to constrain the inverse problem for the whole survey. Liu & Wang (2020) adopt a subspace method for AVO inversion, to mitigate the ambiguity among different elastic parameters.
The purpose of this study is to improve the accuracy of the conventional linearised AVO approximation for strongcontrast interfaces and large incident angles (less than the critical angle) while maintaining the advantage of the low computational cost of the linearised AVO inversion. In using the widely used linear Aki-Richards approximation, if the true elastic parameters are known from the well location, the problematic residuals of the PP reflection coefficients calculated using the Aki-Richards approximation compared with using the Zoeppritz equation will be caused only by the inaccurate Aki-Richards operator. Therefore, we estimate the modification of the Aki-Richards operator by solving a linear inverse problem, then, after applying the modification to the operator, we obtain an improved linear AVO operator (ILAO) with higher accuracy for AVO forward modelling. As described above, the improvement to the conventional linearised AVO approximation requires the use of well-log data to obtain corrections. Therefore, this is a data-driven method. Note that because the Zoeppritz equation under the plane-wave assumption is not suitable for super-critical angles (Alhussain et al. 2008), the large incident angles mentioned next are less than the critical angles.
Since the ILAO is a numerical operator without analytical expression, we construct ILAOs away from well sites using a structure-guided interpolation based on the seismic migrated image. Once we obtained the spatially variant ILAOs for the whole survey, a Bayesian linearised AVO inversion (Buland & Omre 2003) is performed to estimate the elastic parameters from seismic angle gathers. Figure 1 shows the technical solution of the AVO inversion based on the improved linear AVO approximation.

Aki-Richards approximation
In an isotropic and elastic medium, an approximation of the Zoeppritz equation for a weak-contrast interface is (Aki & Richards 1980;Buland & Omre 2003): where In equations (1-4), is the incident angle; V P , V S and represent P-and S-wave velocities and density, respectively; ΔV P , ΔV S and Δ are the corresponding contrasts over a reflected interface andV P ,V S and̄are the corresponding averages. For a time-series model {V P (t), V S (t), (t)}, the extended form of the reflection coefficient in equation (1) is (Stolt & Weglein 1985;Buland & Omre 2003) where t is the travel time and A ′ V P (t, ), A ′ V S (t, ) and A ′ (t, ) are the coefficients in equations (2-4), respectively, with time-dependent velocitiesV P (t) andV S (t). Practically,V P (t) andV S (t) can be represented by a low-frequency initial model. If , then the matrix expression of equation (5) is where ] T is the Aki-Richards operator.

Improved linear AVO approximation
We know that the Aki-Richards approximation is a first-order approximation of the AVO response suitable for weak-contrast media and small incident angles. However, if the true elastic model is known in advance, the residual of the PP reflection coefficients ΔR calculated using the Aki-Richards approximation compared with using the Zoeppritz equation (see Appendix A) will be caused only by the Aki-Richards approximate operator. Therefore, according to equation (6), we have where ΔA(t, ) = [ΔA V P (t, ), ΔA V S (t, ), ΔA (t, )] T is the modification of the Aki-Richards operator. Once ΔA is estimated from equation (7), then, ILAO is Here, we introduce how to estimate ΔA from equation (7). Because there are three components in ΔA for each incident angle, the solution of equation (7) is an under-determined problem. To ensure the over-determination of the inverse problem, we consider a strategy of sliding time window to estimate ΔA. For each incident angle j with a time window with the length of 2 * L + 1, the extended form of equation (7) is where We extend equation (9) to improve the stability of the inverse problem by including the regularisation and smoothness constraints: where is the regularisation parameter, is the smoothing coefficient and I and D are matrices of 6 × 6 for the identity matrix and the first-order differential matrix, respectively. According to equation (11), ΔA i,j is estimated by using the least-square algorithm of a linear inverse problem. Further, we obtain the ΔA j for each given incident angle j using a sliding time-window strategy.

Construction of spatially variant ILAOs
The ILAO A ′′ can be estimated for every sampling point of the wells with coordinates of x. Using a set of K estimated ILAOs , the interpolation is applied to extend the ILAOs laterally for each seismic angle gather away from the well sites. The interpolation process is to construct a function q(x) : R n → R using the known samples (A ′′k , x k ), such that q(x k ) = A ′′k (k = 0, 1, ..., K). Structureguided interpolation is expressed as (Niu et al. 2018) where D(x) is the tensor field that represents the orientation, coherence and dimensionality of stratigraphic and structural features in the migrated section and is the weight parameter called 'tension' that can control the smoothness and accuracy of the interpolation result.
By solving equation (12) using a preconditioned conjugate gradient algorithm, we obtain the function q(x) to map the coordinate space to the value space of the ILAO. Then, we obtain spatially variant ILAOs for 2D or 3D prestack AVO inversions.

Linearised AVO inversion
A linearised AVO inversion method (Buland & Omre 2003) is performed to invert the elastic parameters from prestack seismic data (see Appendix B). The singular difference between the proposed inversion method improved by the data-driven process and the method proposed by Buland & Omre (2003) is that the ILAO is used instead of the Aki-Richards operator.

AVO forward modelling by ILAO
We test the effectiveness of the improved linear AVO approximation on real well-logging data with a sampling interval of 0.1 m. Figure 2a-c shows the well logs of P-and S-wave velocities and density converted to time domain with a sampling interval of 1 ms. Figure 2d illustrates the velocity reflectivities beyond the two dashed lines corresponding to strong-contrast media, and the strongest contrast interface is indicated by a black arrow. Figure 3a and 3b shows the PP reflection coefficients of all reflected interfaces generated using the Zoeppritz equation and the Aki-Richards approximation, respectively. The maximum incident angle is 56°, which is smaller than the critical angle. As illustrated in figure 3d, the Aki-Richards equation approximates the Zoeppritz equation relatively well at small incident angles and weak-contrast interfaces, while the residuals of the PP reflection coefficients evidently increase with the incident angle, becoming even as large as 0.27 for the incident angle of 56°at the strong-contrast interface. Based on the residuals of the PP reflection coefficients in figure 3d, we can estimate the modification of the Aki-Richards operator by solving equation (11), and then obtain the ILAO. The regularisation parameter and the smoothing coefficients used for estimating the ILAO are 0.005 and 0.01, respectively. The basis of selecting and is explained in the discussion. Figure 3c shows the PP reflection coefficients of all the reflected interfaces generated using the ILAO. Compared with figure 3d, the residuals of the PP reflection coefficients generated using the ILAO in figure 3e are much smaller, especially at large incident angles. By comparing figure 2d with figure 3, we also find a significant improvement in the AVO forward modelling using the ILAO for strong-contrast interfaces. Figure 4 shows a comparison of the Aki-Richards operator and the ILAO for the three components. We find that the components of the Aki-Richards operators are relatively smooth, while the corresponding components of the ILAO show more complex changes over time, especially at large incident angles. Furthermore, it should be noted that according to equations (2-4), the component A V P of the Aki-Richards operator is constant at each incident angle, and the components A V S and A are in perfect positive linear correlation, while, contrastingly, the three components of the ILAO are found to vary with time and incident angle, which is thought to better characterise the complexity of the subsurface media. In this way, the ILAO effectively improves the accuracy of AVO forward modelling.

AVO inversion by ILAO
We use the linearised AVO inversion method to test the validity of the improved linear AVO approximation for estimating the elastic parameters. We generate the synthetic noise-free angle gathers by a zero-phase Ricker wavelet with a dominant frequency of 35 Hz, convoluting with the PP reflection coefficients calculated using the Zoeppritz equation. The range of incident angles is from 0 to 56°, while the random noise added in the synthetic data follows a standard normal distribution. Figure 5a-d shows the synthetic seismic angle gathers with signal-to-noise ratios (S/N) of 100, 10, 5, and 2, respectively. We generate the initial models of P-and S-wave velocities and density used for inversion by smoothing the real well logs in figure 2. Figure 6 illustrates the inversion results using the Aki-Richards operator and the ILAO shown in figure 4 for different noisy angle gathers in figure 5. We find that the inverted results using the proposed approximation match the true well logs substantially better than those using the Aki-Richards approximation. Compared with the velocity reflectivities shown in figure 2d, the improved inverted results using the ILAO show an increased correlation with the strong-contrast media. Moreover, the resolution of the inversion results improves significantly using the proposed approximation. We adopt the correlation coefficient between the true well logs and estimated results as a quantitative parameter for evaluating the inversion results. Figure 7 shows that the correlation coefficient of using the ILAO is higher than that of using the Aki-Richards operator, even with S/N being as low as 2.

Field data application
We locate a 2D field seismic data with an S/N of 8 to test the proposed data-driven AVO inversion method, which consists of 1000 angle gathers with a sample interval of 1 ms and a CDP spacing of 25 m in the Tarim Basin, northwest China. The range of incident angles is from 2 to 40°, and two wells from the seismic section are applied to estimate the ILAOs at the well sites. Figure 8 shows the small-, medium-and large-angle stack sections. The vertical black lines in figure 8 indicate the positions of wells X and Y, and the corresponding CDP locations are 459 and 1260, respectively. The upper and lower bounds shown in figure 8 are horizons A and B, respectively. The structure of this region is relatively gentle, except for the faults. Besides, the dominant frequency of the statistical seismic wavelet is 14 Hz. Figure 9 shows the logging curves of P-and S-wave velocities, density, the velocity reflectivities for the P-and S-waves (i.e. ΔV P ∕V P and ΔV S ∕V S ) and the corresponding seismic-well-tie angle gathers at the sites of wells X and Y, respectively. The logging curves are also constrained between horizons A and B, where the velocity reflectivities between the two dashed black lines in figure 9d and 9i correspond to weak-contrast media, while the area beyond the two dashed black lines corresponds to strong-contrast media. Three strong-contrast segments are identified for well X and four for well Y, shown as green rectangles. The variation trends of the elastic parameters of well X and well Y are similar, meaning that the lithology between the two wells varied gradually as opposed to abruptly.
We estimate the ILAOs at the two well sites using the data-driven method. The regularisation parameter and smoothing coefficient are 0.01 and 0.02, respectively. Figure 10 parts a-c and d-f illustrate the difference between the conventional Aki-Richards operator and the ILAO for wells X and Y, respectively. By comparing figures 9 and 10, we find that ILAO is Figure 9. P-and S-wave velocity, density, the velocity reflectivities for P-wave (red line) and S-wave (blue line), and the corresponding seismic-well-tie angle gather for well X (a−e) and well Y (f−j). In d and i, the two dashed black lines constrain the velocity reflectivities between −0.1 and 0.1, while green rectangles indicate the strong-contrast media and the rest are weak-contrast media. nearly identical to the Aki-Richards operator for weak-contrast media, while for strong-contrast media (labeled A, B, C and D in the figures), the differences between the Aki-Richards operator and ILAO increase as the contrast becomes stronger. Figure 11 illustrates the validity of the ILAOs estimated at well X and well Y for AVO forward modelling. We use the PP reflection coefficients calculated by the Zoeppritz equation as references. In figure 11, the accuracy of the PP reflection coefficients calculated by the ILAO is substantially higher than that calculated by the Aki-Richards operator, especially at large incident angles. Note that for the interfaces where the residuals of the PP reflection coefficients are quite large in figure 11a and 11d, there is a significant improvement for AVO forward modelling using the ILAOs ( figure 11b and 11e).
Next, we laterally interpolate the ILAOs estimated at the two wells to the entire target region. The interpolation is constrained by the structural features of the seismic migrated image. Figure 12 shows the conventional Aki-Richards operators and the interpolated ILAOs along the 2D seismic section at three different incident angles. Comparing figure 12a-c 11 Figure 11. Analysis of the AVO forward operator at all reflected interfaces for well X (a−c) and well Y (d−f). Parts (a) and (d) are the differences of the PP reflection coefficients calculated between the conventional Aki-Richards operator and the Zoeppritz equation. Parts (b) and (e) are the differences of the PP reflection coefficients calculated between the ILAO and the Zoeppritz equation. Parts (c) and (f) are the root-mean-square (RMS) error of the PP reflection coefficients calculated using the Aki-Richards operator (blue line) and the ILAO (red line).
with figure 12d-f, we find that the Aki-Richards operators are very smooth, which only include low frequencies, while the interpolated ILAOs contain more high-frequency information, revealing the complexity of the subsurface media, especially at large incident angles and strong-contrast interfaces.
Finally, we estimate P-and S-wave velocities and density by using the linearised AVO inversion method. To compare the effect of the conventional Aki-Richards operators and the interpolated ILAOs on the inversion results, we use the same wavelet, prior information and initial models. Comparisons of the inverted P-and S-wave velocities and densities (figures 13-15) show that the resolution of the inversion results for strong-contrast media using the interpolated ILAOs is higher than that using the conventional Aki-Richards operators, as indicated by the black arrows. In this way, certain thin geological features can be better characterised using the proposed inversion method. For weak-contrast media, the results of inversion with the Aki-Richards operator and the interpolated ILAOs are nearly identical.
Comparisons of the inverted results at well X and well Y are shown in figures 16 and 17, respectively. We find that the inversion results using the interpolated ILAO match better with the true logging data, while the inversion results using the Aki-Richards operators are relatively poor. In addition, there is a good relationship between the strong-contrast media and the improvement of inversion results using the interpolated ILAOs, which demonstrates the advantage of the improved linear AVO approximation over the conventional Aki-Richards approximation for large incident angles and strong-contrast media.

Selecting the regularisation parameters and smoothing coefficients
To estimate the ILAO, we need to determine the values of the regularisation parameters and smoothing coefficients . Figure 18 illustrates the accuracy of the ILAO varying with different and with the model parameters shown in figure 2.
With increasing values of and , the accuracy of the ILAO decreases. Despite this, the accuracy of ILAO is still higher than that of the Aki-Richards operator, especially at large incident angles, as illustrated in figure 19. Figure 6 demonstrates that 12 Figure 12. Comparison of the three components A V P , A V S , and A of the conventional Aki-Richards operators (a-c), and those of the interpolated ILAOs (d-f) along the 2D seismic section. even with a smooth ILAO estimated by setting and to 0.005 and 0.01, respectively, the inverted results using the improved linear AVO approximation still have better agreement than those using the Aki-Richards approximation.

Constructing spatially variant ILAOs
Once the ILAOs are estimated at the well sites, the spatially variant ILAOs can be constructed for the whole survey by structure-guided interpolation. Relatedly, a good interpolation method generates good spatially variant ILAOs, and the choice of interpolation method depends on the complexity of seismic structural features and the computational efficiency of interpolation. Hale (2009) proposed image-guided blended neighbor interpolation by solving nonlinear equations, while Wu (2017) removed the faulting and folding during building 3D subsurface models to follow discontinuous structures across faults. Some commercial software, such as Ikon's Joint Impedance and Facies Inversion and CGG-Jason's RockMod also provide good well-log interpolation tools. In this study, we apply the interpolation method proposed by Niu et al. (2018), which not only has high computing efficiency in terms of solving linear equations, but also works well for characterising thin layers.
In the case of unique geological conditions, such as karst caves, the estimation of ILAOs should be carefully processed because there will be relatively strong local anomalies in the logging data. Our suggestion is to apply the smoothed logging 13 Figure 13. Inverted results of P-wave velocity. Black arrows indicate improvements of the inversion results using the interpolated ILAOs compared with using the conventional Aki-Richards operators for the strong-contrast media. (a) Background model of P-wave velocity, (b) estimated P-wave velocity by the conventional Aki-Richards operators and (c) estimated P-wave velocity by the interpolated ILAOs. curves instead of the real logging data at these local sites to estimate the ILAOs, which prevent the interpolation from bringing excessive strong local anomalies at the wells into the ILAOs away from the wells.

Amplitude scale factor
For the inversion of real seismic data, although the estimation of the ILAO is independent of the seismic wavelet, the amplitude scale factors matching synthetic and real data affect the accuracy of the inversion results. Figure 20 shows the effect of the amplitude scale factor on the inversion results at well Y, in which the green and red lines represent the inversion results by the Aki-Richards operator and ILAOs with the optimal amplitude scale factor, respectively, and the blue lines represent the inversion results by the ILAOs with an amplitude scale factor error of 20%. The comparison of red and blue lines illustrates that the inversion results with an inaccurate amplitude scale factor have a relatively poorer match to the true logging data. Despite this, since the ILAO improves the accuracy of the reflection coefficients, the accuracy and resolution of the inversion results with an inaccurate amplitude factor are still higher than those obtained using the Aki-Richards operator, as shown by the blue and green lines in figure 20.

Resolution of the inverted results
The Aki-Richards equation is known as an approximate linear physical relationship that only reflects the first-order AVO behavior. Therefore, the smooth (low-frequency) Aki-Richards operator cannot be applied for large incident angles and strong-contrast media (see figures 3 and 11). However, because logging data directly reveals the complex changes in the subsurface, this known information can be used to modify the Aki-Richards operator (data-driven) (equations (7-8)), thereby bringing high-frequency logging information into the estimated ILAOs (figures 4, 10 and 12). In this way, the estimated ILAO includes higher frequency information compared to the Aki-Richards operator.
Here, we want to emphasise that the ILAO rich in high-frequency information is different from the high-frequency initial model for AVO inversion, in that, because the estimation of the high-frequency ILAO is the result of a combination of the Aki-Richards operator and the well-log data, the high-frequency ILAO A ′′ includes not only the Aki-Richards physical relationship A ′ , but also the modification of the Aki-Richards operator ΔA, which is rich in high-frequency information from the logging data (see figures 4 and 10). Contrastingly, the high-frequency initial model only includes high-frequency information from the logging data. From equations (B1) and (B2) (in Appendix B), we find that the resolution of the inverted results depends on not only the prestack seismic data d obs , but also the ILAO A ′′ . Therefore, we apply high-cut filtering to the ILAO and keep the low-frequency components of the high-cut filtered ILAO consistent with the Aki-Richards operator. The inverted results using the high-cut filtered ILAO are nearly identical to those using the Aki-Richards operator with the 15 Figure 15. Inverted results of density. Black arrows indicate improvements of the inversion results using the interpolated ILAOs compared with using the conventional Aki-Richards operators for the strong-contrast media. (a) Background model of density, (b) estimated density by the conventional Aki-Richards operators and (c) estimated density by the interpolated ILAOs. same low-frequency initial model ( figure 21). However, if the high-frequency components of the high-cut filtered ILAO are still higher than those of the Aki-Richards operator, the resolution of the inverted results using the high-cut filtered ILAO is also higher, compared to that using the Aki-Richards operator with the same low-frequency initial model ( figure 22). Therefore, it is the high-frequency ILAO that drives improvement of the resolution of the inversion results by providing the AVO forward operator with more high-frequency information.

Computational efficiency
Compared to the linear AVO inversion, the computational costs of estimating ILAOs at the well sites and interpolating ILAOs is relatively small. For this field data application, the estimation and interpolation of the ILAOs only required 4 minutes, while the elastic parameter inversion required 280 minutes on a PC computer with a 3.60 GHz CPU and 8 GB of RAM. Because the size of the Aki-Richards operator matrix is exactly equal to that of the ILAO matrix, the computational cost of the AVO inversion using the ILAO is also identical to that using the Aki-Richards operator; that is, the computational efficiency of the proposed data-driven inversion method is consistent with that of the linearised AVO inversion. 16 Figure 16. Inverted results at well X. Gray, dashed green, blue and red lines represent the actual well logs, initial model, estimated results by the Aki-Richards operators and estimated results by the ILAOs, respectively. Black arrows indicate improvements of the inversion results using the ILAOs compared with using the conventional Aki-Richards operators. Green rectangles correspond to the strong-contrast segments.

Data-driven analysis by real seismic records
This study expected to obtain an improved linear AVO approximation that is closer to the Zoeppritz equation than the Aki-Richards approximation through data-driven approach. However, for super-critical angles, the Zoeppritz equation under plane-wave assumptions is typically not sufficient to describe the spherical-wave reflections (Alhussain et al. 2008). Therefore, in further research, we may consider the pore-elastic, anisotropic or dispersive wave equations instead of Zoeppritz equation modelling, or possibly using real seismic records as the basis of the AVO behaviors at well sites in order to estimate ILAOs in more complex media.

Conclusions
We propose an improved linear AVO approximation using a data-driven method. This method is based on the concept that the residuals of the PP reflection coefficients calculated using the linear Aki-Richards approximation compared with using the Zoeppritz equation can be applied in order to modify the Aki-Richards operator. Because ILAOs can only be estimated at well sites, we construct spatially variant ILAOs for the whole survey based on a structure-guided interpolation constrained by the seismic migrated data.
Synthetic example results show that the improved linear AVO approximation is substantially more accurate than the conventional Aki-Richards approximation for AVO forward modelling for large incident angles and strong-contrast media. Additionally, the quality and resolution of the inverted P-and S-wave velocities and density using the proposed approximation are also substantially higher than those using the Aki-Richards approximation. The 2D field data application confirms the validity and effectiveness of the proposed method for estimating elastic parameters from prestack seismic angle gathers. Figure 17. Inverted results at well Y. Gray, dashed green, blue and red lines represent the actual well logs, initial model, estimated results by the Aki-Richards operators and estimated results by the ILAOs, respectively. Black arrows indicate improvements of the inversion results using the ILAOs compared with using the conventional Aki-Richards operators. Green rectangles correspond to the strong-contrast segments.  Analysis of the effect of the amplitude scale factor on inversion results at well Y. Gray lines represent the actual well logs, green and red lines represent the inversion results by the Aki-Richards operator and the ILAOs with the optimal amplitude scale factor, respectively, and blue lines represent the inversion results by the ILAOs with an amplitude scale factor error of 20%. Black arrows indicate improvements of the inversion results using the ILAOs compared with using the conventional Aki-Richards operators. Green rectangles correspond to the strong-contrast segments.
19 Figure 21. Comparison of the inverted results using the Aki-Richards operator and the high-cut filtered ILAO, while the frequency band of the high-cut filtered ILAO is consistent with that of the Aki-Richards operator. Gray, dashed green, blue and red lines represent the actual well logs, initial model, estimated results by the Aki-Richards operator and estimated results by the high-cut filtered ILAO, respectively. Figure 22. Comparison of the inverted results using the Aki-Richards operator and the high-cut filtered ILAO, while the high-frequency components of the high-cut filtered ILAO are still higher than those of the Aki-Richards operator. Gray, dashed green, blue and red lines represent the actual well logs, initial model, estimated results by the Aki-Richards operator and estimated results by the high-cut filtered ILAO, respectively.

The Zoeppritz equation is
cos 1 − sin 2 cos 2 cos 1 − sin 1 cos 2 sin 2 cos 2 1 − V S 1 V P 1 sin 2 1 − 2 V P 2 1 V P 1 cos 2 2 − 2 V S 2 1 V P 1 sin 2 2 sin 2 1 where V P , V S and represent P-and S-wave velocities and density, and subscript numbers 1 and 2 represent upper and lower layers, respectively; 1 and 1 are reflected angles of the P-and SV-waves, and 2 and 2 are transmitted angles of the P-and SV-waves, respectively; R PP and R PS are reflected coefficients of the P-and SV-waves and T PP and T PS are transmitted coefficients of the P-and SV-waves, respectively.