Seismic amplitude inversion for the transversely isotropic media with vertical axis of symmetry

Transverse isotropy with a vertical axis of symmetry is a common form of anisotropy in sedimentary basins, and it has a significant influence on the seismic amplitude variation with offset. Although exact solutions and approximations of the PP‐wave reflection coefficient for the transversely isotropic media with vertical axis of symmetry have been explicitly studied, it is difficult to apply these equations to amplitude inversion, because more than three parameters need to be estimated, and such an inverse problem is highly ill‐posed. In this paper, we propose a seismic amplitude inversion method for the transversely isotropic media with a vertical axis of symmetry based on a modified approximation of the reflection coefficient. This new approximation consists of only three model parameters: attribute A, the impedance (vertical phase velocity multiplied by bulk density); attribute B, shear modulus proportional to an anellipticity parameter (Thomsen's parameter ε−δ); and attribute C, the approximate horizontal P‐wave phase velocity, which can be well estimated by using a Bayesian‐framework‐based inversion method. Using numerical tests we show that the derived approximation has similar accuracy to the existing linear approximation and much higher accuracy than isotropic approximations, especially at large angles of incidence and for strong anisotropy. The new inversion method is validated by using both synthetic data and field seismic data. We show that the inverted attributes are robust for shale‐gas reservoir characterization: the shale formation can be discriminated from surrounding formations by using the crossplot of the attributes A and C, and then the gas‐bearing shale can be identified through the combination of the attributes A and B. We then propose a rock‐physics‐based method and a stepwise‐inversion‐based method to estimate the P‐wave anisotropy parameter (Thomsen's parameter ε). The latter is more suitable when subsurface media are strongly heterogeneous. The stepwise inversion produces a stable and accurate Thomsen's parameter ε, which is proved by using both synthetic and field data.

enable amplitude variation with offset (AVO) analysis and modelling of VTI media. Approximations of reflection coefficients in anisotropic media are extensively investigated to provide expressions in simple and closed forms (Banik 1987;Thomsen 1993;Ursin and Haugen 1996;Rüger 1997;Rüger 1998;Vavryčuk and Pšenčík 1998;Zillmer, Gajewski and Kashtan 1998;Vavryčuk 1999;Stovas and Ursin 2003;Shaw and Sen 2004;Zhang and Li 2013). In spite of that, the inversion of amplitudes of reflected PP waves in VTI media is a difficult issue. Plessix and Bork (2000) studied the feasibility of anisotropic parameters estimation from the reflection coefficients and demonstrated that the recovery of all five parameters is difficult in practice. Stovas, Landrø and Avseth (2006) invert the AVO intercept and gradient to estimate net-to-gross and oil saturation for finely layered reservoirs. Lin and Thomsen (2013) presented a method to extract the anisotropy parameter δ based on the difference between real (measured) and synthetic seismic amplitudes.
In this paper, we propose a modified PP-wave reflection coefficient equation that enables the inversion for VTI media.
The new approximation has similar accuracy to the existing linear approximation and contains only three subsurface parameters to be inverted: attribute A, impedance (vertical phase velocity multiplied by bulk density); attribute B, shear modulus proportional to an anellipticity parameter (Thomsen's parameter ε−δ); and attribute C, approximate horizontal P-wave phase velocity. A simultaneous inversion scheme is applied for the estimation of the three parameters, and the misfit function is regularized by a priori information about the model in a Bayesian framework. The stability and uncertainty analyses are performed on inversion results of synthetic data tests and field-data application to show the feasibility of this approach. The stability analysis is based on the relative error and the cross-correlation coefficient between the inversion result and true model (real log), and the uncertainty is represented by the estimated standard deviation of the inverted parameter. The study area is in the south of Sichuan Basin, China, and the target is a lower Silurian-age shale formation. We show that the shale formation can be effectively discriminated from the surrounding limestone and sandstone by using the crossplot of the impedance and the horizontal P-wave phase velocity; the gas-bearing shale can be easily identified through the combination of the anisotropic shear modulus and the impedance. Then, the P-wave anisotropy parameter ε is recovered from inverted subsurface parameters by two methods based on a rock-physics relation and a stepwise inversion, respectively. The former applies an explicit rock-physics relation between V P0 and impedance (attribute A) to estimate the vertical P-wave phase velocity V P0 , which is then used to estimate the anisotropy parameter ε along with the simultaneously inverted vertical P-wave phase velocity (attribute C), while the latter uses small-angle seismic reflection data to invert the vertical P-wave phase velocity V P0 .

Reflection coefficients in the transversely isotropic media with vertical axis of symmetry
The elastic stiffness tensor of a VTI (transverse isotropy with vertical symmetry axis) medium is be represented as (1) Thomsen (1986) suggests a parameterization that enables to analyse the influence of anisotropy on seismic signatures. Thomsen's anisotropy parameters are defined as Five elastic stiffness coefficients of VTI media in equation (1) can be expressed as where V P0 ,V S0 and ρ are the P-and S-wave vertical phase velocities and the bulk density, respectively. Considering a discrete subsurface model, the linear approximation of the PP-wave reflection coefficient for the ith planar interface separating the ith and the i + 1th VTI layers can be expressed as (Ursin and Haugen 1996;Rüger 1997;Vavryčuk and Pšenčík 1998;Stovas and Ursin 2003;Shaw and Sen 2004;Zhang and Li 2013) where θ is the incidence phase angle, and can be approximated to the average angle θ of the incidence phase angle  Table 1. and the transmission phase angle under the assumption of weak impedance contrast and small angle, ρV P0 refers to average values of the two adjacent layers across the interface, refers to contrasts in the properties of the two adjacent layers across the interface, R ISO (θ ) is the linear approximation of PP-wave reflection coefficient for isotropic media (Bortfeld 1961;Richards and Frasier 1976;Aki and Richards 1980) We aim to propose a modified approximation of the PPwave reflection coefficient that enables a simultaneous inversion for the VTI medium. So the new equation shall contain fewer model parameters to be inverted than the existing approximation. For this purpose, equation (4) is firstly expressed as where K = (2V S0 /V P0 ) 2 . Note that sin 2 θ (1 + tan 2 θ )=tan 2 θ and sin 2 θ tan 2 θ =tan 2 θ −sin 2 θ , the above equation is then rearranged as The first term in the above equation can be expressed by the reflection coefficient of normal incidence: Equation (8) can be expressed as by taking the following approximations: Equation (8) is then written as The difference between when |r i |=0.3. Similar approximations are used by Bortfeld (1961) and Oldenburg et al. (1983). We can also have Substituting equations (11) and (12) into (7), we have the following expression: Tsvankin and Thomsen (1994) defined an effective parameter as σ = (V P0 /V S0 ) 2 (ε − δ), and ε − δ is known as the anellipticity parameter. Letting we have the following expression: From here on, we refer to equation (15) as 'the new equation', and A, B and C are three attributes to be inverted for. Similar to the reflection coefficient functions in Stolt and Weglein (1985) and Buland and Omre (2003), equation (15) contains model parameters in the differential form of natural logarithms. It provides explicit forms of inverted attributes, especially for attributes B and C, so that we can easily understand their physical meanings and capabilities of lithology and pore fluids. Besides, actual values of elastic parameters, instead of their relative perturbations with respect to the background models, can be directly estimated by using a novel simultaneous inversion scheme. This eliminates the error induced by the transformation from the latter to the former. Three attributes have explicit physical meaning. The attribute A represents the impedance. The attribute B is a shear modulus multiplied by e σ 4 in terms of the anellipticity parameter (ε − δ). It can be similar to the isotropic shear modulus at moderate anisotropy, when the value of e σ 4 is close to one (e.g. e σ 4 = 1.088 according to the shale parameters in Table 1, in which ε − δ=0.1 and V P0 /V S0 =1.84). The difference between B = ρV S0 2 e σ 4 and ρV S0 2 becomes significant at strong anisotropy (e.g. e σ 4 = 1.18 when ε − δ = 0.2 and V P0 /V S0 =1.84). The value of e σ 4 is determined by rock properties and has a significant influence on the attribute B. An empirical relationship between V P0 /V S0 and ε − δ has been studied by Ryan-Grigor (1997). It shows that shales have measured V P0 / V S0 between 1.6 and 1.8, and some shales do have ε − δ close to zero. However, this elliptical anisotropy is not observed from the data used in this paper. The attribute C = V P0 e ε can be equivalent to the horizontal P-wave phase velocity V P90 (V P90 ≈ V P0 (1+ε)) when ε < 0.3.

Simultaneous inversion method for the transversely isotropic media with a vertical axis of symmetry
Three model parameters in the new equation can be solved in a Bayesian framework by minimizing a quadratic misfit function J regularized by an a priori model of Gaussian where g is the nonlinear forward modelling operator, d is the n × N-by-1 data vector consisting of seismic amplitude as a function of the incident phase angle at each time sample. The number of angle traces and the number of data samples at each trace are n and N + 1, respectively; m = [A, B, C] T is Figure 4 The synthetic seismic angle gathers with different levels of noise: (a) without noise, (b) signal-to-noise ratio = 2, (c) signal-to-noise ratio = 1 and (d) signal-to-noise ratio = 0.5.

Figure 5
Comparison between the real logs (black solid), initial models (grey) and inversion results (black dashed) for synthetic data with different levels of noise: (a) without noise, (b) signal-to-noise ratio = 2, (c) signal-to-noise ratio = 1 and (d) signal-to-noise ratio = 0.5. A = ρV P0 , B = ρV S0 2 e σ 4 and C = V P0 e ε are shown from left to right.
the 3(N +1)-by-1 model vector consisting of the three model parameters to be solved, m prior is the prior model parameter vector; C d is the data covariance matrix and could be practically calculated by assuming uncorrelated noise as C d = σ 2 n I, the model covariance matrix that can be obtained from real logs or rock-physics relations, and C kl (k and l refer to attributes A, B and C) is the covariance of inverted parameters. Given an initial model m prior , the forward modelling can be approximated as where G is the n × N-by-n × N mapping operator from model to data, and it includes the effect of incident angles and wavelets and is given by (Buland and Omre 2003;Downton and Ursenbach 2006;Zong et al. 2012) where W is the n × N-by-n × N convolution matrix containing wavelets of different angles, F is the n × N-by-3N sensitivity matrix of Fréchet derivatives with respect to the model parameters, and D is the 3N-by-3(N + 1) differential operator Table 2 Relative error (RE) and cross-correlation coefficient (CC) between the true models and the inversion results of attributes A, B and C, respectively, and the standard deviation (STD) of parameter estimates for different signal-to-noise ratios of the model parameters in adjacent layers (Stolt and Weglein 1985). The model parameter can be recovered as is the data residual, μ proportional to σ n is the trade-off parameter of the regularization term and can be estimated using varied methods to guarantee an optimal solution (Wang 1999;Downton 2005; Alemie and Sacchi 2011).

E X P E R I M E N T S W I T H S Y N T H E T I C A N D F I E L D D A T A S E T S Accuracy analysis
The linear reflection coefficient approximations of VTI (transverse isotropy with vertical symmetry axis) media are compared with the exact solution and the approximation of isotropic media to analyse the accuracy of the new equation (Fig. 1). Figure 2 shows the relative errors between different approximations and the exact anisotropic reflection coefficient. The exact PP-wave reflection coefficient of VTI media is calculated using the equation in Schoenberg and Protázio (1992). The isotropic approximation is calculated as in equation (5). Table 1 shows the elastic parameters of four single-interface models belonging to different amplitude variation with offset (AVO) classes according to Rutherford and Williams (1989). All the model parameters are derived from well logs in the study area. The class I and II models represent typical shale/sand interfaces; classes III and IV represent the interfaces of shale/gas-bearing shale and limestone/shale, respectively. Thomsen's anisotropy parameters of the shale formation are ε = 0.15, δ = 0.05 and γ = 0.2, according to the well logs. The limestone, sand and gasbearing shale show very weak anisotropy and thus they are considered as isotropic. The models of classes III-a and IV-a Figure 6 Estimated parameter variance (uncertainty) based on the covariance matrix of the new reflection coefficient equation as a function of the maximum angle for different signal-to-noise ratios: (a-d) No noise, signalto-noise ratio = 2, 1 and 0.5, respectively. The minimum angle is zero, and the data are sampled uniformly over this angle range.

Figure 7
Well logs of P-wave velocity, S-wave velocity, bulk density, Thomsen's ε, δ, γ (real log in black, estimation in grey), porosities of clay-bound water, free water and free gas, and mineralogy from left to right.
have the same elastic parameters as classes III and IV, respectively, but the anisotropy parameters of the shale are set to ε = 0.3 and δ = 0.1, to evaluate the accuracy at strong anisotropy. The approximation of isotropic media (equation (5)) deviates from the exact solution, and such deviation increases with angle and the magnitude of anisotropy.
The new approximation (equation (15)) agrees well with the existing equation (equation (4)), and both of them have high accuracy over a wide range of angles (from 0°to 50°) as shown in Figs 1 and 2. In the following sections, we would discuss the feasibility of the inversion scheme using synthetic and real seismic data sets.

Figure 8
The post-stack seismic section and well location.

Inversion results using a synthetic data set
Real well logs in the study area are used to generate the synthetic data to test the inversion scheme, and the logs are transformed from the depth to the time domain (Fig. 3). Similar to the numerical results in Fig. 1, the reflection coefficients are calculated using the exact solution in Schoenberg and Protázio (1992). The forward modelling uses the mixed-phase wavelets estimated from the corresponding seismic sections of a constant angle. Figure 4 shows the synthetic angle gathers with Gaussian noise of different signal-noise ratios. Each gather consists of five traces from 10°to 50°, which is consistent with the real seismic data set in the study area. Figure 5 shows a comparison between the real logs and inverted attributes. Reliable low-frequency initial models are important for obtaining stable inversion results. The real logs are smoothed by a moving average filter with a 40-point span to generate the initial models. The inversion stability is analysed by using relative error (RE) and cross-correlation coefficient (CC) between the true model and the inversion result, and the uncertainty of inversion result is analysed by using estimated standard deviation (STD) of an inverted parameter (Table 2). RE is calculated as where m real (i) and m inv (i) refer to the true and inverted parameters of the ith time sample, respectively. STD, representing the uncertainty of model parameters, is calculated from the square Figure 9 The well logs are calibrated with the seismic section at CDP400, the cross-correlation coefficient between the synthetic waveform (blue) and real seismic waveform (black and red) is 0.83. root of model covariance matrix diagonal elements (model parameter variance), and the model covariance matrix can be estimated as a function of data noise variance and the linear operator F asCm = σ 2 n [F T F ] −1 by assuming uncorrelated uniform noise (Downton 2005). All the three attributes (A = ρV P0 , B = ρV S0 2 e σ 4 and C = V P0 e ε ) can be satisfactorily inverted from the synthetic data with a high signal-to-noise ratio (signal-to-noise ratio >2). Res and STDs increase with the decrease of signal-to-noise ratio, while CCs decrease with the signal-to-noise ratio. All three inverted attributes tend to be very patchy when using noisy data (signal-to-noise ratio <1), and attribute B has more error and larger uncertainty than those of attributes A and C. The above analysis implies that an inversion scheme can be feasible for a field-data  application using the seismic data set with a signal-to-noise ratio >1.
We now analyse the effect of varying the maximum angle in the uncertainties of three inverted model parameters. The uncertainty is represented by the estimated variance the same as the above analysis. Figure 6 shows the uncertainty as a function of the maximum angle for different signal-to-noise ratios. The minimum angle is zero, and the angle used for this analysis is sampled uniformly over the range. The variance of attribute A is very low even using a small maximum angle. The uncertainties of attributes B and C are very similar and are much higher than that of attribute A when using small maximum angles. So for limited angle ranges, the inversion problem can be ill-conditioned and relies on the model con-straint. The uncertainties of attributes B and C show a rapid decline with the maximum angle used to increase. For high signal-to-noise ratio (>2), uncertainties of B and C decrease to small amounts within small maximum angle (30°); while for low signal-to-noise ratio (<1), reliable estimations of B and C require data with a large maximum angle (>40°). In the following field-data experiment, seismic data sets from 10°to 50°are used to ensure the reliability of all three parameters.

Field seismic data inversion and shale-gas reservoir characterization
A field seismic prestack data set from a shale-gas reservoir survey is used to demonstrate the feasibility of our proposed methods. The study area is in the south of Sichuan Basin, in southwest China. Well 3 penetrates a lower Silurian-age shale formation (S 1 l) that is at depth 2080-2390 m and located below a carbonate formation consisting of limestone and calcareous sandstone and deposited on an Upper Ordovician limestone formation (Fig. 7). Mineralogy log shows that S 1 l shale (2080-2390 m) has moderate clay contents and contains a varied volume of quartz and carbonate minerals. The lower S 1 l formation (2210-2390 m) is organic rich and gas saturated, and it consists of a large volume of quartz and has much higher porosity and kerogen volume than the upper S 1 l formation (2080-2210 m). The log of the anisotropy parameter γ shows that this shale formation has strong to moderate anisotropy. The anisotropy parameter γ is calculated using the C 44 and C 66 measured in the borehole. C 44 is calculated using the measured shear-wave velocity and density, and C 66 is  estimated using Stoneley waveforms. The anisotropy parameters ε and δ are not measured directly in the borehole, so we apply a rock-physics-based method to predict the anisotropy parameters for the shale formation (Qian et al. 2016;Zhang, Qian and Li 2017;Zhang 2017Zhang , 2019. Logs of mineralogy, porosities of clay-bound water, free water and free gas, and TOC are used to model all five effective elastic constants of the VTI shale by applying anisotropic effective-medium theories. Three Thomsen's anisotropy parameters are then calculated using the estimated elastic constants, and the well logs of attributes A, B and C are calculated using the measured elastic parameters and predicted anisotropy parameters. The accuracy of predicted anisotropy parameters is demonstrated by the agreement between the predicted γ and the measured γ . The correlation coefficient between the measured and predicted logs is 0.7514. The mean, minimum, maximum values and standard deviation of the absolute difference between the measured and predicted logs are 0.0098, 0.0001, 0.1705 and 0.0426, respectively. The post-stack seismic section along with well location is shown in Fig. 8. The well logs are calibrated with the seismic section at around CDP400 (Fig. 9). The synthetic waveforms and real seismic waveforms are plotted in blue and red (black), respectively. The cross-correlation coefficient between them is 0.83. The shale formation can be interpreted within the time window from 1280 to 1460 ms at the well location. It has generally lower impedance (much lower velocities and slightly lower density) but the higher magnitude of anisotropy than the upper and lower formations. The organic-rich and gas-saturated layer at the bottom of S 1 l (1445-1460 ms) has lower impedance and weaker magnitude of anisotropy than the upper S 1 l formation. The amplitude variation with offset (AVO) corresponding to the top of the gas-bearing shale and the top of the shale formation belong to classes III and IV, respectively, as shown in Fig. 1(c,d). The calibrated well logs are smoothed and extrapolated along with the interpreted horizons to generate the initial models of the inversion. The well logs are also used to calculate the model covariance for regularization. We now use the well logs to discuss the interpretation capability of the inverted attributes. The crossplot of A-C, representing the impedance and V P0 e ε , respectively, is compared with that of A-V P0 in Fig. 10. Logs of three formations (calcareous sandstone, shale, limestone) within the range from 2000 to 2410 m are used for this analysis. In the crossplot of A and C, it is easy to identify the samples of shales (clay volume > 15%, AI < 14 (km/s) × (g/cc), C < 6.0 km/s) from surrounding sandstones and limestones (Fig. 10a). This implies that attributes A and C are useful in shale formation prediction because shales show a distinct A-C relation from upper and lower formations due to its generally pronounced VTI properties. It is difficult to discriminate the shale from the surrounding formations by using the crossplot of AI and V P0 , because there is a good linear relation between them. Once the shale formation is identified, we need to find the fluid anomalies within the shale formation. The crossplots of A-B and AI-S are compared in Fig. 11. Only the logs within the shale formations (2080-2390 m) are used for this analysis. The gassaturated shale shows higher values of the shear modulus and the attribute B than those of the upper shales. Attribute B has a similarly good interpretation capability as the shear modulus in identifying the gas-saturated shale. Note that an interpretation using inverted parameters from the seismic inversion of anisotropic approximation can be more capable than those from the inversion of isotropic approximation because the latter can be biased by the anisotropy.
As to the field-data inversion, the original prestack seismic gathers are transformed from the offset domain to the angle domain, and seismic data within certain angle intervals are stacked to construct the constant-angle sections (Fig. 12). All of the constant-angle seismic sections from 10°to 50°are used to invert three parameters. The calibrated well logs are smoothed to generate the initial models of the inversion. They are also used to calculate the model covariance for regularization. The inversion results of three parameters (A = ρV P0 , B = ρV S0 2 e σ 4 and C = V P0 e ε ) are shown in Fig. 13. The inversion results at CDP400 are compared with corresponding logs as shown in Fig. 14(a-c). Similar to the synthetic data test, the stability and uncertainty analyses are performed to the inversion results. The relative errors (RE) of inverted attributes A, B and C are 4.72%, 8.61% and 4.2%, respectively; the crosscorrelation coefficients (CCs) of three inverted attributes are 0.8746, 0.8273 and 0.8648, respectively; the standard deviations (STDs) or uncertainty of three inverted attributes are 0.0372, 0.0741 and 0.0351, respectively. Although the results of this field-data application have more errors and uncertainties than those of synthetic tests, there is a generally good match between inverted parameters and real logs. Possible improvements may achieve if applying further data processing procedures on the original gathers to reduce noise in offset domain and flatten reflection events. The synthetic angle gather generated using the inverted parameters is also compared with the real seismic angle gather (Fig. 14d,e). The cross-correlation coefficient between the synthetic angle gather and the real angle gather is 0.8202, and the root-mean-square error of data is 0.0114. It is noted that the logs of ε and δ are estimated using the rock-physics-based method mentioned above, and the logs of attributes B and C are calculated using the measured V P0 , V S0 and ρ in boreholes and logs of ε and δ. The shale formation has lower impedance than the upper sandstone and lower limestone formations thus can be interpreted within the time window from 1280 to 1460 ms at the well location. The shale-gas reservoirs are approximately within the range 1440-1460 ms, in which the inverted attributes show lower value compared with the upper shale layers.

Inversion of the P-wave anisotropy parameter ε
Anisotropy parameters are important in imaging and inversion. In this section, we propose two schemes to recover the anisotropy parameter ε hidden in the attribute C. There is an explicit relationship between the vertical P-wave phase velocity and the bulk density for sedimentary rocks (Gardner, Gardner and Gregory 1974;Castagna et al. 1993;Wang 2000). For the study area, a linear relationship between the vertical P-wave phase velocity and the attribute A (impedance) is found by investigating field well-log data for varied sedimentary formations, including shale, limestone and sandstone (Fig. 15a): where c 1 = 0.3251 and c 2 = 0.5657. In this sense, the anisotropy parameter ε can be estimated by using the two inverted attributes A and C. The procedure is shown in Fig. 15(b), in which the vertical P-wave phase velocity V est P0 is estimated from the inverted attribute A using equation (18), and then the anisotropy parameter ε is calculated as ln( C V est

P0
). The REs between inversion results and real logs of V P0 and ε are 2.2% and 9.09%, respectively; and corresponding CCs are 0.97 and 0.959, respectively.
Stepwise inversion is an alternative to the previous rockphysics-based method and is more applicable when the properties of subsurface media are strongly heterogeneous. The independent anisotropy parameter ε can be recovered from the inverted parameter C by using two steps: (1) simultaneous inversion using the isotropic AVO equation for the vertical Pand S-wave phase velocities and density using small-angle reflection seismic data and isotropic approximation as expressed in equation (4); (2) the anisotropy parameter ε is estimated using the inverted vertical P-wave phase velocity (V P0 ) and the parameter C = V P0 e ε as ln( C V est

P0
). The PP-wave seismic amplitude at small angles is less sensitive to anisotropy than at large angles (Alkhalifah and Tsvankin 1995), and thus it can be assumed to be weakly anisotropic or isotropic. Thus, simultaneous inversion of isotropic PP-wave reflection coefficient approximation using only small-angle data ensures a reliable estimation of V P0 . The synthetic test result of the stepwise inversion is shown in Fig. 16, in which the vertical P-wave phase velocity V P0 is inverted from only small-angle data (10°-30°), and the anisotropy parameter ε is calculated using C and V P0 . The REs between this stepwise inversion results and real logs for V P0 and ε are 0.85% and 8.81%, respectively; corresponding CCs are 0.9885 and 0.9765, respectively; and the STD of the V P0 estimate is 0.000247. So the stepwise inversion method provides more reliable result of V P0 than that of the rock-physics-based method, and the inverted anisotropy parameter ε is also improved accordingly. This stepwise inversion method is then applied to field data. Three parameters (A = ρV P0 , B = ρV S0 2 e σ 4 and C = V P0 e ε ) have been inverted by using the simultaneous inversion method. The vertical Pwave phase velocity (V P0 ) is inverted from small-angle seismic data (10°-30°) and the anisotropy parameter ε is estimated accordingly (Fig. 17a,b). There is a good agreement between the inverted V P0 and the real log, the corresponding RE and CC are 4.27% and 0.9885, respectively, and the STD of inverted V P0 is 0.036. Although the inversion result of attribute C is similarly good (RE = 4.2%, CC = 0.8648 as shown in the above subsection) as V P0 , the inverted anisotropy parameter ε has a relatively large error (RE = 10.37%). Even though there is a generally good match between inverted ε and the real log, the corresponding cross-correlation coefficient is 0.7204.

D I S C U S S I O N A N D C O N C L U S I O N
Amplitude inversion for all five parameters of the VTI (transverse isotropy with vertical symmetry axis) medium is difficult in practice. To enable the inversion of the VTI medium, we first propose a new approximation of the PP-wave reflection coefficient that includes three attributes: the impedance, the shear modulus proportional to an anellipticity parameter (ε − δ) and the approximate horizontal P-wave phase velocity. Inverted attributes show good interpretation capability in shale-gas reservoir characterization: shale formation can be easily interpreted according to the nonlinear relationship in the crossplot of attribute A and C = V P0 e ε , gas-bearing shale is identified through the combination of attribute A and B = ρV S0 2 e σ 4 . The difference between the attribute B and shear modulus can be significant when ε − δ > 0.2, even if V P0 /V S0 is small, and the importance of such a difference could be further discussed in the future study. Numerical analysis shows that the new equation has good accuracy over a very wide range of angles, and inversion of the three attributes is performed by using regularization in a Bayesian framework. The P-wave anisotropy parameter ε can be recovered by using a rock-physics-based method and a stepwise-inversion-based method. Both methods produce stable and high-resolution anisotropy parameter ε, but the stepwise-inversion-based method is more suitable and is proved by using both synthetic data and field data. The proposed methods based on the new reflection coefficient approximation provide an efficient way of dealing with the multipleparameter inversion for anisotropic media, and may also be applicable in more complex anisotropic media, e.g. TTI (transverse anisotropy with tilted axis) or orthorhombic media.
The input field seismic data are normal moveout (NMO) corrected common-image-point gathers in the angle domain. In this sense, appropriate methods, for example, the nonhyperbolic moveout equation and double scanning, have been used for NMO correction (e.g. Tsvankin and Thomsen 1994;Alkhalifah and Tsvankin 1995;Grechka and Tsvankin 2002). The simultaneous inversion is performed by using a priori information on the model as the regularization. The practical procedures of conventional seismic inversion of isotropic media, such as wavelet estimation, well ties and initial-model building can be used in this inversion. Note that if there are no borehole measurements of anisotropy parameters (ε and δ) available, rock-physics-based methods can be used to predict the logs of the anisotropy parameters for generating reliable models. The final results show an encouraging comparison with the borehole measurements.

A C K N O W L E D G E M E N T S
This study is sponsored by the National Science and Technology Major Project (2017ZX05018005), the National Natural Science Foundation of China (No. 41474096), CNPC Science Research and Technology Development Project (2019A-3308) and the Science Foundation of China University of Petroleum, Beijing (2462018BJC001).