Joint PP and PS AVA seismic inversion using exact Zoeppritz equations

The technique of amplitude variation with offset (AVO), or amplitude variation with incidence angle (AVA), has been used to predict lithology and fluid properties for decades (Ostrander, 1982; 1984; Sun et al., 2008). The Zoeppritz equations describe the relationship between reflection and transmission coefficients for a given incidence angle, and elastic media properties (P- and S-wave velocities, as well as density). Because the Zoeppritz equations are nonlinear functions with respect to these properties, many approximations have been made to linearize them to simplify calculations. These approximations have played a key role in the estimation of elastic parameters by AVO/AVA inversions.

P-wave AVO/AVA methods have been developed over decades and remarkable advances have been achieved. Since Ostrander (1982) proposes an AVO method to extract lithology and pore-fluid information from prestack seismic amplitudes, AVO has been used with varying rates of success (Larsen, 1999). Shuey (1985) develops a gradient-intercept method that estimates zero-offset reflectivity and changes in Poisson’s ratio. Smith and Gidlow (1987) estimate rock properties by a weighted stacking method using time- and offset-variant weights to PP data samples before stacking. Roberts et al. (2002) use an AVO method to enhance the reservoir characterization based on the long-offset seismic data. However, a few years of experience in P-mode AVO observation leads to results that are sometimes too ambiguous to interpret (Jin, 1999).

Compared with P-wave approaches alone, joint AVO inversion of PP and PS data can provide better estimates of elastic parameters (Kurt, 2007); however, most AVO studies use two-parameter inversion methods because three-parameter inversions are often plagued by numerical instability (Ursenbach and Stewart, 2008). Stewart (1990) extends the weighted stacking method to invert PP- and PS-waves jointly for rock property estimations. Fatti et al. (1994) improve the weighted stacking method to directly estimate fractional P and S impedance, but their method requires the angle of incidence less than 35° and P- to S-wave velocity ratios between 1.5 and 2.0. Xu and Bancroft (1997) express the P-P and P-S reflection coefficient equations with the elastic parameters (⁠$λ$ and $μ$⁠, or $κ$ and $μ$⁠) and extract these parameters directly without conversion from velocity. Larsen (1999) presents a method to invert PP and PS AVO gathers simultaneously to estimate P and S impedance. Ursenbach and Stewart (2008) derive inversion error expressions for various two-parameter inversion methods, and obtain flexible conversion formulas that can convert the results of any two-parameter method to those of any other two-parameter method. However, two-parameter methods lose important density information, an AVO attribute that is useful for inferring fluid saturation (Downton, 2005).

For this reason, a significant amount of work has been done on elastic three-parameter estimations. Jin (1999) presents a real data example of a joint P- and S-wave AVO analysis that estimates elastic parameters for reservoir characterization, and proposes that S-wave data gave the opportunity to estimate S-wave velocity and density variations. Jin et al. (2000) use singular-value decomposition (SVD) to stabilize the linearized PS system of equations, and they obtain good results for synthetic and field data. Buland and Omre (2003) develop a new linearized AVO inversion method using Bayesian techniques to obtain P- and S-wave velocities, as well as density distributions. Mahmoudian and Margrave (2004) formulate the joint AVO inversion using SVD methods to obtain P- and S-impedance contrasts, as well as density contrast. Three-parameter inversion methods are also investigated by Downton and Lines (2004), Mahmoudian (2006), and Veire and Landrø (2006), but all these studies solve approximations of the Zoeppritz equations. The underlying assumption of these approximations is that the incidence angle is not too large and the contrast in elastic parameters is weak (Aki and Richards, 1980; Zhi et al., 2013); however, these approximations affect the accuracy of calculations, and they are particularly unsuitable in the case of a strong contrast between layers, such as that between a coal seam and unconsolidated sandstone. Therefore, a joint AVO/AVA inversion method using exact Zoeppritz equations would represent a more robust method of estimating elastic parameters using multicomponent seismic data (Zhi et al., 2013).

Some attempts at solving the exact Zoeppritz equations have been made. For example, Tigrek et al. (2005) develop a method to relate the offset-dependent seismic reflections to the local stress distributions, and obtain accurate result by joint AVA analysis based on the full-Zoeppritz equations. Wang et al. (2011) estimate model parameters with the exact Zoeppritz equations using a generalized linear inversion. Zhu and McMechan (2012) propose a new algorithm using the exact Zoeppritz equations to obtain four independent parameters (a density ratio and three velocity ratios) for PP reflections. Zhi et al. (2013) rewrite the exact Zoeppritz equations as a function of four independent parameters to invert PP- and PS-waves jointly; however, these inversion results are elastic parameter contrasts, not the elastic properties themselves.

In this paper, we present a method of joint AVA inversion using the exact Zoeppritz equations for PP and PS seismic data. The inverted parameters are P- and S-wave velocities, as well as density. Using a least-squares approach, we establish a joint objective function that is used to judge whether the difference between simulated and observed records is small enough to end the inversion. To apply the method to field data, some preprocessing is needed, such as true-amplitude processing of the PP and PS data and compression of the PS AVA data set to the PP traveltime domain. The performance of the proposed method is determined by comparison with the synthetic and field seismic data. We also compare the results of the linear and nonlinear joint AVA inversion methods.

Additionally,$i1$ and $i2$ are the P-wave angles of incidence and transmission across the interface, respectively, and $j1$ and $j2$ denote the S-wave angles of incidence and transmission, respectively. The value $p$ is the ray parameter, which is constant for all layers, and $p=sini1/α1=sini2/α2=sinj1/β1=sinj2/β2$ (Aki and Richards, 1980).

To ensure the inversion accuracy, a mean shift method was applied to constrain the inversion result. After some iterations, we derive an elastic parameter model that is close to the prior model by extrapolating from the well model along horizons. If the ratio between the mean value of the inverse elastic parameters and the prior model parameters is large, then this ratio is used to scale the inverse model to improve its proximity to the prior model. Mean shift can be considered as the low-frequency compensation of the model.

To apply the above method to field PP- and PS-waves, special processing techniques are needed for: (1) P- and S-mode separation from vertical and horizontal components, (2) PS to PP time alignment, and (3) wavelet extraction at different incidence angles.

In multicomponent seismology, multiple field shot records are measured (shown as “mode leakage”), allowing for potential cross-contamination of P-wave energy in the horizontal component and S-wave energy in the vertical component. Suppressing mode leakage by P- and S-wave separation is an important step in prestack multicomponent seismic data processing. For this study, we applied the method introduced by Lu et al. (2012) to separate different wave modes.

The inversion output is the elastic parameter vector, rather than the change in an elastic parameter, so wavelets must be defined. Some prestack processing procedures, such as Q-compensation, event stretching, and compression, may lead to the changes of PP- and PS-wavelets at each incidence angle. Wavelets should therefore be extracted from PP and PS data sets separately for the joint inversion, at the total incidence angle. Most other wavelet-extraction methods for PP-waves are suitable for our inversion method. For example, a statistical method of kurtosis maximization by constant-phase rotation developed by Baan (2008) can produce a stationary, constant-phase wavelet (Jonathan and Baan, 2011).

To test our method against synthetic data, we used a five-layer 1D model with high-contrast interfaces (shown in Table 1), and synthetic PP and PS AVA data sets within critical angles, as shown in Figure 1. To correlate the PP and PS events, we change the polarities of PS-waves to match those of PP-waves when they are opposite. These data sets are generated with a 1-ms sample rate through convolution of the reflectivity derived from the exact Zoeppritz equation and Ricker wavelets with dominant frequencies of 40 Hz for PP-waves and 30 Hz for PS-waves in PP time, respectively. Interference, such as noise and multiples, was not taken into account in the forward modeling. Two different inversion models were tested: (1) the nonlinear joint PP and PS AVA inversion based on exact PP and PS reflection coefficients presented above (NJI) and (2) a linear joint PP and PS AVA inversion based on Aki and Richards’ (1980) approximate PP and PS reflection coefficients (LJI). The linear single PP AVA inversion based on Aki and Richards’ (1980) approximate PP reflection coefficients (LSI) is also tested. But the results of LSI are closed to that of LJI, so we just show the results of NJI and LJI for comparison. For joint inversions, the weight factor $w$ was set to 0.5 (i.e., the two data sets were assumed to have equal quality). For nonlinear and linear inversions, the Jacobian matrix $G$ was constructed using equation 8, although the nonlinear inversion used the exact Zoeppritz solutions, whereas the linear inversion used Aki and Richards’ (1980) approximations.

A linear model and a noise-added model (shown in Figure 2) were used to test the sensitivity of our inversion method to the choice of the initial model. The inversion results (without any mean shift corrections) with the linear initial model using the (a) NJI and (b) LJI methods are shown in Figures 3–5 for P- and S-wave velocities, as well as density, respectively. Although the stratigraphic contrast cannot be seen on the plot of the linear initial model (Figure 2a), NJI (Figures 3a, 4a, and 5a) produces a good result, showing a thin layer that exactly conforms to the true model. In contrast, LJI (Figures 3b, 4b, and 5b) produce poor results because the true model has strong contrast parameters and a small stratum thickness. Figures 6–8 display inversion results for the noise-added initial model, and show that NJI (Figures 6a, 7a, and 8a) produces stable results compared with LJI (Figures 6b, 7b, and 8b).

Figures 3–8 show that inversion methods based on exact Zoeppritz solutions pick up the thin low-impedance layer in the true model. The inverted P- and S-wave velocities are closer to the true model values than are the inverted densities because density exhibits little sensitivity to amplitudes and the use of angle-limited AVA data may have affected its stability (Lines, 1998; Alemie and Sacchi, 2011; Du and Yan, 2013). However, for the model described in Table 1, larger incidence angles will cause critical reflectivity. A mean shift method, which is used to constrain the inversion procedure to improve the inversion accuracy, was applied to the NJI model results when the initial model was linear. As shown in Figure 9, the inversion constrained by the mean shift method is virtually identical to the true model.

Any acquired seismic record will contain a certain amount of noise, and inversion accuracy is highly dependent upon this noise level. To assess the sensitivity of our method to noise, we add random noise with different weights to the synthetic PP and PS seismograms separately. As shown in Figure 10, a noise level of 16% was generated randomly in the time domain, with noise values distributed normally around the synthetic signal.

Figure 11 shows the NJI results overlaid with the true values when 1% noise was added to the AVA records. The inversion result is close to the true model, including the thin layer with large elastic contrasts, even though it is noisy. To suppress the noise for the inversion, we smoothed the inversion model values using a piecewise smoothing algorithm to obtain an improved inversion result. First, as shown in Figure 12, a low-pass filter is applied to get the trend line of the inversion results. Second, the turning points are distinguished from the trend line for dividing the trend line into several segments. Finally, the inversion result is smoothed in each segment. The smoothing results at the noise level of 1% are shown in Figure 13, and the results reveal the thin layer more effectively. The smoothing method also produces stable NJI results when a noise level of approximately 16% was added to the synthetic data (Figure 14).

Our inversion method has also been applied to a field multicomponent data. The inversion results are also compared with real well-log data. The field data are 2D with 1 ms sample rate and 118 common depth points (CDPs) in the Sichuan Basin, West China, and the CDP spacing is 10 m. The angle gather at the well location is shown in Figure 15. At the well location, P- and S-waves, as well as density logs were acquired and carefully processed. Before applying the NJI method, we carry out processing procedures, preserving the true amplitude of the seismic data. The first step is a static correction, which can flatten PP and PS events during the conversion of offset gather to incidence angle gather.

The second step is to calibrate PP and PS events reflected from each geologic interface in the PP traveltime domain with the assistance of synthetic seismograms calculated from well logs (Figure 16). During this procedure, PP and PS time-to-depth log conversions are calculated and used for PS data set compression. Compared with Figure 16, it can be seen that the PP and PS seismograms in the PP traveltime domain shown in Figure 15 correlate well with synthetic seismograms after event correlation. Though the PS-field AVA data set is formed from the raw data with a worse signal-to-noise ratio, its correlation coefficient with the synthetic data set is up to 0.65 after careful processing. The correlation coefficient between the actual and synthetic PP AVA data sets is 0.76. So, we set the weight factor $w$ in equation 5 to 0.5 because the two correlation coefficients are close.

The final step is to extract wavelets of PP and PS AVA data sets separately. We extract the statistical wavelets over the entire incidence angle range separately after event correlation. Figure 17 shows that both wavelets are zero phase and have similar main frequencies but differing side lobes.

The inverted sections are shown in Figures 18 and 19, where the green and red colors represent the shale and the sand contents of the layer, separately, which are calibrated by well logs. It shows that we can distinguish more thin layers from the sections inverted by the NJI method than by the LJI method. Compared with the other parameters, the velocity-ratio section (Figure 19c) inverted by NJI method shows obviously higher resolution. The near-well inversion results are extracted from the sections and shown in Figure 20. We can see that good calibration at the well position made the near-well inversion results closely match the well logs. Though the inversion results have lower frequency content than the well logs because of the limited frequency of the wavelet, NJI can produce a high-resolution result, which is beneficial to the recognition of thin layer.

Our NJI inversion method requires careful data processing. Some true-amplitude processing techniques are mentioned in the above. Incidence angle is an important factor, which will affect the inversion result seriously. In cases of complicated structures, we can use a ray-tracing method to calculate the raypath at the target interface and derive the accurate incidence angle, which is out of the scope of this paper.

Noise attenuation is also quite important. Noise in the AVA data set will increase the uncertainty in the inversion result, but we can decrease the weight factor of the uncertain component to reduce the noise effect. However, to distinguish thin layer information from the inversion result, we do not want the data to be filtered too hard before inversion.

Our inversion method does not consider anisotropy, so the negative effects of elastic parameter anisotropy on multicomponent seismic data must be eliminated. For vertical transverse isotropy (VTI) anisotropy, we can consider using a nonhyperbolic moveout correction technique to obtain large-angle data. Furthermore, in the case of horizontal transverse isotropy (HTI) anisotropy, phase error, induced by S-wave splitting, must be compensated in prestack processing. Many fast and slow S-wave separation methods can be used to compensate such phase error.

For thin layer models, it is hard to get the interface reflections only because the reflections from the top and bottom interfaces of the thin layer will interfere together and cannot be separated. Moreover, interbed multiples and converted waves are hard to suppress, also they will influence the amplitude and phase of the target reflection. Therefore, to pick up the thin layer from the inversion result of field data, it is necessary to deduce the reflection coefficient formula of the thin layer.

Event matching is crucial to joint PP and PS AVA inversion. The PS wavelet will change with time when compressed to PP time because of the time-variant compression ratio. Therefore, the PS wavelet should be extracted from the compressed data, and the time window chosen for inversion should be short enough to retain stability of the wavelet. It is necessary for us to explore a more effective method to improve event matching and PS-wave extraction in future work.

We have developed a joint AVA inversion method (NJI) using the exact Zoeppritz equations for PP and PS seismic data, and we have tested the model for synthetic and field seismic data. When dealing with high-contrast interfaces, inversion methods based on approximations to the Zoeppritz equations are insufficient, but our NJI method can achieve significantly better results.

The synthetic data example shows that it is possible to reveal a thin layer even with a linear initial model or with small amount of noise. The inverted elastic properties show a high level of random perturbation when a large amount of noise added to the synthetic data; however, the thin layer can still be revealed after a piecewise smoothing operation applied. Our NJI method is suitable for larger incidence angles, so velocity analysis and NMO corrections for large offset data must be applied before AVA gather formation; however, even in small incidence angle cases, the NJI method offers an advantage over approximation methods because it uses fewer theoretical assumptions. Also, the methodology requires noise suppression in AVA gathers, even though stable results can be achieved if the noise level is small.

For field data, our NJI method proved helpful in geologic layer interpretation, even though the PP- and PS-wavelets were narrow. We anticipate that our NJI method would be especially useful in predicting lithology and fluid properties in fields with thin and low-impedance layers.

We greatly appreciate the support of PetroChina Innovation Foundation (2011D-5006-0303), the Natural Science Foundation of China (41104084), and the Fundamental Research Funds for the Central Universities (2-9-2013-093). We are very grateful to Z. Li for sharing the field data. Thanks to T. Chen and an anonymous reviewer, whose comments and suggestions significantly improved the manuscript.