# Nonlinear AVA Inversion Based on a Novel Quadratic Approximation for Fluid Discrimination

Authors: Lin Zhou; Xingye Liu; Tianchun Yang; Jianping Liao; Mingfeng Zhu; Gan Zhang
Journal: Geofluids
 (2020)

---

## Abstract

Fluid discrimination plays an important role in reservoir exploration and development. At present, the fluid factors used for fluid discrimination are estimated by linear AVA inversion methods based on the linear approximations of the Zoeppritz equations. However, the Zoeppritz equations show that the relationship between prestack AVA reflection coefficients and reservoir parameters is highly nonlinear. Therefore, inversion methods based on linear approximations will seriously influence the nonuniqueness and uncertainty of inversion results. In this paper, a nonlinear inversion based on the quadratic approximation is carried out to reduce the nonuniqueness and uncertainty of fluid factor. Firstly, in order to directly invert the fluid factor, a novel quadratic approximation in terms of the fluid factor (ρf), shear modulus, and density on both sides of the reflection interface is derived based on poroelasticity theory. Then, a nonlinear inversion objective function is constructed using the novel quadratic approximation in a Bayesian framework, and the Gauss-Newton method is adopted to minimize this objective function. The synthetic data example shows that the new method can obtain reasonable fluid factor inversion results even in low SNR (signal-to-noise ratio) case. Finally, the proposed method is also applied to field data which shows that it can effectively discriminate reservoir fluids.


---

## Body

## 1. Introduction

Fluid discrimination is an important step in reservoir exploration and development. In order to characterize the reservoir fluid information by geophysical data, various fluid factors have been proposed and applied to fluid discrimination. Smith and Gidlow [1] first defined the concept of fluid factor as the combination of P- and S-wave velocities. Goodway et al. [2] suggested that Lamé petrophysical parameters (λρ, Lamé modulus × density and μρ, shearmodulus×density) were significantly better than the P- and S-wave velocities in fluid detection and lithology identification. Based on poroelasticity theory [3, 4], Russell et al. [5, 6] defined the Russell fluid factor. Russell fluid factor is derived based on the Biot–Gassmann theory and has become the most commonly used and classic fluid factor. Castagna et al. [7] and Li and Castagna [8] studied the application of AVO intercept and gradient in AVO classification and realized reservoir hydrocarbon discrimination based on the AVO classification research results. Based on a physical modeling study, Wandler et al. [9] demonstrated that the AVO intercept and gradient can be used as traditional hydrocarbon indicators. Feng et al. [10] suggested that the fluid factor ρf (density×Russellfluidfactor) had a better performance in fluid discrimination. Yin and Zhang [11] defined the effective pore-fluid bulk modulus to raise the sensitivity of fluid discrimination. Their researches show that the introduction of effective pore-fluid bulk modulus will lead to an increase in the number of target parameters that need to be inverted, which will seriously affect the stability of nonlinear inversion algorithms. Zong and Yin [12] suggested that Young’s impedance (YI) and Poisson ratio impedance (PI) have great potential in fluid discrimination and lithology identification of unconventional reservoirs. However, the physical meaning of these two parameters is undefined. Considering the effectiveness of the fluid factors and the stability of nonlinear inversion algorithms, we select the fluid factor ρf as one of target parameters of the nonlinear inversion method to directly detect the fluids.Prestack seismic data contain abundant and reliable information about fluids and lithology compared to poststack seismic data. Therefore, various amplitude variation with offset (AVO) or amplitude variation with incidence angle (AVA) inversion methods based on prestack seismic data for fluid discrimination have been developed rapidly in recent years. Russell et al. [6] derived a linear approximate formula in terms of f, μ, and ρ (the Russell fluid factor, shear modulus, and density) and achieved linear inversion using standard least-squares inversion method. Zong et al. [13] derived a linear approximation equation that contains the P- and S-wave moduli and then realized the indirect estimation of the Russell fluid factor based on a petrophysical relationship. In order to avoid the cumulative error caused by indirect calculation, Zong et al. [14] realized the direct inversion of Russell fluid factor by using Russell linear approximation. Du and Yan [15] derived new linear approximations which containing the fluid factor ρf for PP and PS wave reflection coefficients and then realized the joint inversion. Yin and Zhang [11], Zong et al. [16], and Du et al. [17] derived the fluid matrix decoupled linear approximation equation and realized the linear inversion of the fluid modulus. Zong and Yin [12] achieved direct inversion of Young’s and Poisson impedances for fluid discrimination. It should be noted that the forward operators of above inversion methods are linear approximation equations of the exact Zoeppritz equations. Actually, the exact Zoeppritz equations show that the relationship between reservoir parameters and prestack AVA reflection coefficients is complicated and highly nonlinear. If nonlinear equation degenerates into linear equation, it is easy to appear the case that the reflection coefficients calculated by the linear and nonlinear equations are similar, but the input parameter values are quite different [18]. Therefore, inversion methods based on linear approximations will cause strong nonuniqueness and uncertainty in inversion results. Stovas and Ursin [19, 20] derived quadratic approximations for reflection and transmission coefficients in isotropic media. Compared with linear approximations, the relationship between reservoir parameters and AVA reflection coefficients described by these quadratic approximations is closer to the real situation. Based on these quadratic approximations, Rabben et al. [18] and Aune et al. [21] achieved nonlinear inversion using different methods. They showed that the nonlinear inversion methods based on quadratic approximations can effectively reduce the nonuniqueness and uncertainty of inversion results, which provides an important theoretical guarantee for the research content of this paper. In order to improve the estimated accuracy of fluid factor and reduce its uncertainty, a novel quadratic approximation of exact Zoeppritz equations is derived on the basis of the existing quadratic approximation. Since the forward operator we used is a nonlinear function, the objective function is also nonlinear. At present, intelligent algorithm (support vector machine, quantum particle swarm, etc.) and deterministic algorithm (Gauss-Newton method, Generalized Linear Inversion method, etc.) [22–28] are popular to solve the nonlinear inversion problems. In order to solve nonlinear objective function quickly and stably, the classical Gauss-Newton method is chosen in this paper.We proposed a nonlinear inversion for the fluid factorρf based on a quadratic approximation. First, a new quadratic approximation of the exact Zoeppritz equations which containing the fluid factor ρf, shear modulus, and density is derived based on poroelasticity theory. Then, a nonlinear objective function is constructed using this quadratic approximation in a Bayesian framework, and the classical Gauss-Newton method is applied to solve this nonlinear problem. Finally, we apply the proposed method on both synthetic and field data, indicating the feasibility and effectiveness of the proposed method and draw some conclusions.


## 2. Theory and Method

### 2.1. Derivation of a Novel Quadratic Approximation in terms ofρf, μ, and ρ

Lithology prediction also plays an important role in reservoir exploration [29]. Shear modulus is an important characteristic parameter for lithology prediction. Therefore, in order to reduce the uncertainty of lithology prediction results, shear modulus μ is also selected as one of the target parameter of the nonlinear inversion. Through detailed comparison, Rabben et al. [18] have proved the effectiveness of the quadratic approximation derived by Stovas and Ursin [19, 20] in reducing the uncertainty of inversion results and improving the accuracy of inversion results. On this basis, we derive a novel quadratic approximation in terms of ρf, μ, and ρ to improve the estimated accuracy of fluid factor.The existing quadratic approximation of the exact Zoeppritz equations in terms ofP- and S-wave impedances and density is given by [18, 20]
(1)rPPθ,φ=12cos2θΔIPIP−4sin2φΔISIS−12tan2θ1−4cos2θγsat2Δρρ+tanθtanφ4γsat21−1+1γsat2sin2θΔISIS2−4γsat21−32+1γsat2sin2θΔISISΔρρ+1γsat21−2+1γsat2sin2θ−14Δρρ2,where IP=ρVP and IS=ρVS represent P- and S-wave impedances, θ is the incidence angle of P wave, φ is the reflection angle of PS converted wave, and ρ is density. γsat=VP/VSsat is the background VP/VSsat ratio of the saturated rock.Firstly, we establish the relationship between the wave impedances and the fluid factorρf where f is Russell fluid factor. The expression of Russell fluid factor f is shown below
(2)f=ρVP2−γdry2ρVS2,where VP, VS, and ρ are the P- and S-wave velocities and density of saturated rock, respectively, γdry2=VP/VSdry2 represents the square of the dry-rock velocity ratio. Russell et al. [5, 6] and Wang [30] had discussed the effect of γdry2 on the reflection coefficients in detail. It is usually estimated from well logging data and is treated as a constant in the inversion algorithm [13, 14]. Then, we have
(3)ρf=ρ2VP2−γdry2ρ2VS2=IP2−γdry2Is2=1−γdry2γsat2IP2.Rearranging equation (3) as
(4)IP2=ρf+γdry2Is2.Taking the derivatives of both sides of equation (4), we can obtain
(5)2IPΔIP=Δρf+2γdry2ISΔIS.Combining equations (3) and (5), the following equation can be obtained
(6)ΔIPIP=121−γdry2γsat2Δρfρf+γdry2γsat2ΔISIS.ForS-wave impedance, we have
(7)IS2=ρμ.Taking the derivatives of both sides of equation (7), we obtain
(8)ΔISIS=12Δμμ+Δρρ.Substituting equations (6) and (8) into equation (1), a new quadratic approximation in terms of fluid factor ρf, shear modulus, and density is derived
(9)rPPθ,φ=14cos2θ1−γdry2γsat2Δρfρf+14cos2θγdry2γsat2−2sin2φΔμμ+14cos2θγdry2γsat2−2sin2φ−12tan2θ+2sin2θγsat2Δρρ+tanθtanφ1γsat21−1+1γsat2sin2θΔμμ2+sin2θγsat2ΔμμΔρρ−14Δρρ2.Based on equation (9), we only estimate the reflectivity of the fluid factor, shear modulus, and density. Then, we need to use these results to perform trace integral operation to obtain the parameter values of the fluid factor, shear modulus, and density. In fact, the trace integral calculation not only depends on the initial values but also brings cumulative errors, thus affecting the final estimation accuracy of the fluid factor and other parameters. Therefore, in order to invert the fluid factor, shear modulus, and density directly, we rewrite equation (9) as the following form
(10)rPPθ,φ=14cos2θ1−γdry2γsat22ρf2−ρf1ρf2+ρf1+14cos2θγdry2γsat2−2sin2φ2μ2−μ1μ2+μ1+14cos2θγdry2γsat2−2sin2φ−12tan2θ+2sin2θγsat22ρ2−ρ1ρ2+ρ1+tanθtanφ1γsat21−1+1γsat2sin2θ2μ2−μ1μ2+μ12+sin2θγsat22μ2−μ1μ2+μ12ρ2−ρ1ρ2+ρ1−142ρ2−ρ1ρ2+ρ12,where subscripts 1 and 2 represent the upper and lower media. Equation (10) is the novel quadratic approximation in terms of the fluid factor, shear modulus, and density on both sides of the interface. In order to test the calculation accuracy of the new quadratic approximation shown in equation (10), two representative models shown in Table 1 are used to compare the reflection coefficients. Using these two models can effectively verify the calculation accuracy of the new quadratic approximation. Figures 1(a) and 2(a) show the comparison of reflection coefficients calculated by the exact Zoeppritz equations, Aki-Richards approximation, and Russell linear approximation, and the new quadratic approximation for the two models’ incidence wave is P wave with incidence from 0°–45°. Figures 1(b) and 2(b) show the difference of the reflection coefficients between the exact Zoeppritz equations and the three approximations, respectively. Note that the novel quadratic approximation of the exact Zoeppritz equations derived is satisfactory in calculation accuracy, which ensures the estimated accuracy of fluid factor.Table 1
Parameters of model 1 and model 2.
LithologyVPm/sVSm/sρg/cm3fGpaμGpaγdry2(-)Model 1Water sand305015952.237.5095.6732.333Gas sand278016652.082.6235.766Model 2Shale289812902.42510.9514.0352.333Gas sand285716662.2753.8386.314Figure 1
Comparison of reflection coefficients calculated by the exact Zoeppritz equations, Aki–Richards approximation, Russell linear approximation, and the new quadratic approximation of model 1. (a) Reflection coefficients calculated by different formulas; (b) difference of the reflection coefficients between the exact Zoeppritz equations and their three approximations, respectively.
(a)(b)Figure 2
Comparison of reflection coefficients calculated by the exact Zoeppritz equations, Aki–Richards approximation, Russell linear approximation, and the new quadratic approximation of model 2. (a) Reflection coefficients calculated by different formulas; (b) difference of the reflection coefficients between the exact Zoeppritz equations and their three approximations, respectively.
(a)(b)


### 2.2. Nonlinear Inversion Based on the Novel Quadratic Approximation

In order to improve the stability and accuracy of nonlinear inversion algorithm, the nonlinear inversion objective function is constructed in a Bayesian framework. The likelihood functionPd∣m and the prior distribution function Pm are assumed following the Gaussian distribution and Cauchy distribution, respectively. Their expressions can be expressed as follows
(11)Pd∣m∝exp−12σd2d−QmTd−Qm,(12)Pm∝exp−2∑i=1Nln1+m−μTΦim−μ,where d represents observed data, Q represents the nonlinear forward operator (the new quadratic approximation), m=fluidfactor;shearmodulus;density3N×1 is the target parameters vector, and σd2 is the noise variance. N is the size of target parameters, and μ is the mean vector of target parameters. Φi=DiTψ−1Di, of which ψ is a 3×3 covariance matrix that contains the statistical correlations among the three target parameters. Matrix D is a 3×3N matrix defined as
(13)Djki=1,ifj=1andk=i,1,ifj=2andk=i+N,1,ifj=3andk=i+2N,0,other.Substituting equations (11) and (12) into Bayes’ formula [31], we have
(14)Pm∣d∝Pd∣mPm∝exp−12σd2d−QmTd−Qm−2∑i=1Nln1+m−μTΦim−μ,where Pm∣d is the posterior probability distribution. Then, we need to find the maximum value of Pm∣d. This problem can be converted into a problem that solves the minimum value of the following objective function
(15)Om=12d−QmTd−Qm+βRm,where Rm=2∑i=1Nln1+m−μTΦim−μ is the regularization term and β=σd2 controls the weight of the prior information.The classical Gauss-Newton method is adopted to solve the above nonlinear problem. According to the Gauss-Newton method, the iterative formula can be expressed as(16)mj+1=mj−Hmj−1γmj,where γm and Hm are the first and second partial derivatives of the nonlinear objective function with respect to the parameter vector m, respectively. Their expressions are shown below
(17)γm=∂Om∂m=∂Qm∂mTQm−d+β∂Rm∂m,(18)Hm=∂2Om∂m2≈∂Qm∂mT∂Qm∂m+β∂R2m∂m2.In equations (17) and (18), the matrix ∂Qm/∂m is the partial derivative of the nonlinear forward model based on the new quadratic approximation with respect to m (Jacobian matrix). The derivation of this matrix (∂Qm/∂m) is given in Appendix A. The matrices ∂Rm/∂m and ∂R2m/∂m2 are the first and second partial derivatives of the regularization term Rm with respect to m. Their derivations are given in Appendix B.


## 2.1. Derivation of a Novel Quadratic Approximation in terms ofρf, μ, and ρ

Lithology prediction also plays an important role in reservoir exploration [29]. Shear modulus is an important characteristic parameter for lithology prediction. Therefore, in order to reduce the uncertainty of lithology prediction results, shear modulus μ is also selected as one of the target parameter of the nonlinear inversion. Through detailed comparison, Rabben et al. [18] have proved the effectiveness of the quadratic approximation derived by Stovas and Ursin [19, 20] in reducing the uncertainty of inversion results and improving the accuracy of inversion results. On this basis, we derive a novel quadratic approximation in terms of ρf, μ, and ρ to improve the estimated accuracy of fluid factor.The existing quadratic approximation of the exact Zoeppritz equations in terms ofP- and S-wave impedances and density is given by [18, 20]
(1)rPPθ,φ=12cos2θΔIPIP−4sin2φΔISIS−12tan2θ1−4cos2θγsat2Δρρ+tanθtanφ4γsat21−1+1γsat2sin2θΔISIS2−4γsat21−32+1γsat2sin2θΔISISΔρρ+1γsat21−2+1γsat2sin2θ−14Δρρ2,where IP=ρVP and IS=ρVS represent P- and S-wave impedances, θ is the incidence angle of P wave, φ is the reflection angle of PS converted wave, and ρ is density. γsat=VP/VSsat is the background VP/VSsat ratio of the saturated rock.Firstly, we establish the relationship between the wave impedances and the fluid factorρf where f is Russell fluid factor. The expression of Russell fluid factor f is shown below
(2)f=ρVP2−γdry2ρVS2,where VP, VS, and ρ are the P- and S-wave velocities and density of saturated rock, respectively, γdry2=VP/VSdry2 represents the square of the dry-rock velocity ratio. Russell et al. [5, 6] and Wang [30] had discussed the effect of γdry2 on the reflection coefficients in detail. It is usually estimated from well logging data and is treated as a constant in the inversion algorithm [13, 14]. Then, we have
(3)ρf=ρ2VP2−γdry2ρ2VS2=IP2−γdry2Is2=1−γdry2γsat2IP2.Rearranging equation (3) as
(4)IP2=ρf+γdry2Is2.Taking the derivatives of both sides of equation (4), we can obtain
(5)2IPΔIP=Δρf+2γdry2ISΔIS.Combining equations (3) and (5), the following equation can be obtained
(6)ΔIPIP=121−γdry2γsat2Δρfρf+γdry2γsat2ΔISIS.ForS-wave impedance, we have
(7)IS2=ρμ.Taking the derivatives of both sides of equation (7), we obtain
(8)ΔISIS=12Δμμ+Δρρ.Substituting equations (6) and (8) into equation (1), a new quadratic approximation in terms of fluid factor ρf, shear modulus, and density is derived
(9)rPPθ,φ=14cos2θ1−γdry2γsat2Δρfρf+14cos2θγdry2γsat2−2sin2φΔμμ+14cos2θγdry2γsat2−2sin2φ−12tan2θ+2sin2θγsat2Δρρ+tanθtanφ1γsat21−1+1γsat2sin2θΔμμ2+sin2θγsat2ΔμμΔρρ−14Δρρ2.Based on equation (9), we only estimate the reflectivity of the fluid factor, shear modulus, and density. Then, we need to use these results to perform trace integral operation to obtain the parameter values of the fluid factor, shear modulus, and density. In fact, the trace integral calculation not only depends on the initial values but also brings cumulative errors, thus affecting the final estimation accuracy of the fluid factor and other parameters. Therefore, in order to invert the fluid factor, shear modulus, and density directly, we rewrite equation (9) as the following form
(10)rPPθ,φ=14cos2θ1−γdry2γsat22ρf2−ρf1ρf2+ρf1+14cos2θγdry2γsat2−2sin2φ2μ2−μ1μ2+μ1+14cos2θγdry2γsat2−2sin2φ−12tan2θ+2sin2θγsat22ρ2−ρ1ρ2+ρ1+tanθtanφ1γsat21−1+1γsat2sin2θ2μ2−μ1μ2+μ12+sin2θγsat22μ2−μ1μ2+μ12ρ2−ρ1ρ2+ρ1−142ρ2−ρ1ρ2+ρ12,where subscripts 1 and 2 represent the upper and lower media. Equation (10) is the novel quadratic approximation in terms of the fluid factor, shear modulus, and density on both sides of the interface. In order to test the calculation accuracy of the new quadratic approximation shown in equation (10), two representative models shown in Table 1 are used to compare the reflection coefficients. Using these two models can effectively verify the calculation accuracy of the new quadratic approximation. Figures 1(a) and 2(a) show the comparison of reflection coefficients calculated by the exact Zoeppritz equations, Aki-Richards approximation, and Russell linear approximation, and the new quadratic approximation for the two models’ incidence wave is P wave with incidence from 0°–45°. Figures 1(b) and 2(b) show the difference of the reflection coefficients between the exact Zoeppritz equations and the three approximations, respectively. Note that the novel quadratic approximation of the exact Zoeppritz equations derived is satisfactory in calculation accuracy, which ensures the estimated accuracy of fluid factor.Table 1
Parameters of model 1 and model 2.
LithologyVPm/sVSm/sρg/cm3fGpaμGpaγdry2(-)Model 1Water sand305015952.237.5095.6732.333Gas sand278016652.082.6235.766Model 2Shale289812902.42510.9514.0352.333Gas sand285716662.2753.8386.314Figure 1
Comparison of reflection coefficients calculated by the exact Zoeppritz equations, Aki–Richards approximation, Russell linear approximation, and the new quadratic approximation of model 1. (a) Reflection coefficients calculated by different formulas; (b) difference of the reflection coefficients between the exact Zoeppritz equations and their three approximations, respectively.
(a)(b)Figure 2
Comparison of reflection coefficients calculated by the exact Zoeppritz equations, Aki–Richards approximation, Russell linear approximation, and the new quadratic approximation of model 2. (a) Reflection coefficients calculated by different formulas; (b) difference of the reflection coefficients between the exact Zoeppritz equations and their three approximations, respectively.
(a)(b)


## 2.2. Nonlinear Inversion Based on the Novel Quadratic Approximation

In order to improve the stability and accuracy of nonlinear inversion algorithm, the nonlinear inversion objective function is constructed in a Bayesian framework. The likelihood functionPd∣m and the prior distribution function Pm are assumed following the Gaussian distribution and Cauchy distribution, respectively. Their expressions can be expressed as follows
(11)Pd∣m∝exp−12σd2d−QmTd−Qm,(12)Pm∝exp−2∑i=1Nln1+m−μTΦim−μ,where d represents observed data, Q represents the nonlinear forward operator (the new quadratic approximation), m=fluidfactor;shearmodulus;density3N×1 is the target parameters vector, and σd2 is the noise variance. N is the size of target parameters, and μ is the mean vector of target parameters. Φi=DiTψ−1Di, of which ψ is a 3×3 covariance matrix that contains the statistical correlations among the three target parameters. Matrix D is a 3×3N matrix defined as
(13)Djki=1,ifj=1andk=i,1,ifj=2andk=i+N,1,ifj=3andk=i+2N,0,other.Substituting equations (11) and (12) into Bayes’ formula [31], we have
(14)Pm∣d∝Pd∣mPm∝exp−12σd2d−QmTd−Qm−2∑i=1Nln1+m−μTΦim−μ,where Pm∣d is the posterior probability distribution. Then, we need to find the maximum value of Pm∣d. This problem can be converted into a problem that solves the minimum value of the following objective function
(15)Om=12d−QmTd−Qm+βRm,where Rm=2∑i=1Nln1+m−μTΦim−μ is the regularization term and β=σd2 controls the weight of the prior information.The classical Gauss-Newton method is adopted to solve the above nonlinear problem. According to the Gauss-Newton method, the iterative formula can be expressed as(16)mj+1=mj−Hmj−1γmj,where γm and Hm are the first and second partial derivatives of the nonlinear objective function with respect to the parameter vector m, respectively. Their expressions are shown below
(17)γm=∂Om∂m=∂Qm∂mTQm−d+β∂Rm∂m,(18)Hm=∂2Om∂m2≈∂Qm∂mT∂Qm∂m+β∂R2m∂m2.In equations (17) and (18), the matrix ∂Qm/∂m is the partial derivative of the nonlinear forward model based on the new quadratic approximation with respect to m (Jacobian matrix). The derivation of this matrix (∂Qm/∂m) is given in Appendix A. The matrices ∂Rm/∂m and ∂R2m/∂m2 are the first and second partial derivatives of the regularization term Rm with respect to m. Their derivations are given in Appendix B.


## 3. Examples

### 3.1. Synthetic Data Example

First, a synthetic profile with single-well logs which are obtained from a sandstone layer in an actual oil field in China is used to verify the feasibility and stability of the new method. The value ofγdry2 is set to 2.333. The parameter curves of the single-well are shown in Figure 3, where the solid lines represent the real fluid factor, shear modulus, and density curves, and the dashed lines represent the initial model data that are obtained by smoothing the true curves. Then, the seismic forward modeling is implemented by convoluting a Ricker wavelet (the main frequency is 30 Hz) with the reflection coefficients generated by the exact Zoeppritz equations. The corresponding noise-free synthetic angle gathers are shown in Figure 4(a). Figure 5(a) shows the inversion results of the proposed method using noise-free data. The new method can converge to an optimal solution in only 4 to 5 iterations. We can see that in the absence of noise, all of the fluid factor, shear modulus, and density can be inverted with a high accuracy. In order to test the stability of the new method, random noises with signal-to-noise ratio (SNR) of 2, 1, and 0.5 are added to the synthetic data. These noise data with different SNR are shown in Figure 4 (b: SNR=2, c: SNR=1, d: SNR=0.5). Figures 5(b)–5(d) show the corresponding inversion results of noise data with corresponding SNRs. We can see that the fluid factor and shear modulus can be estimated stably with a satisfactory accuracy, even when SNR=0.5. However, the estimated density results show more bias when noise exists, suggesting that the estimation of density parameters is seriously affected by noise.Figure 3
True model data (solid line) and initial model data (dashed line).Figure 4
Synthetic prestack angle gathers: (a) noise-free, (b)SNR=2, (c) SNR=1, and (d) SNR=0.5.Figure 5
Inversion results of synthetic data with different SNR. The black line indicates the initial model, the blue line indicates the inversion result, and the red line indicates the true data. (a) Noise-free, (b)SNR=2, (c) SNR=1, and (d) SNR=0.5.
(a)(b)(c)(d)In order to further demonstrate the advantage of the proposed method, the inversion based on Russell linear approximation is implemented on the synthetic data shown in Figure4(a), and the inversion results are shown in Figure 6. From the comparison between Figures 5(a) and 6, it is obvious that the inverted density has been improved, and the inversion accuracy of fluid factor, and shear modulus in the range of 0 to 0.2 s also has been improved, which indicates that the proposed method can effectively improve the accuracy of the inversion results.Figure 6
Inversion results of the traditional method based on Russell linear approximation.


### 3.2. Field Data Example

A small 2D field seismic angle gather data (effective angle range is 3°-34°) which is extracted from a work area in China is used to further demonstrate the feasibility and availability of the proposed method in field data. The stack section of the field data is shown in Figure7. In this figure, the black line denotes the location of Well A, and the black arrow indicates the location of the target reservoir. The target reservoir is a shale layer, located between two sandstone layers. A series of conventional processing procedures, such as amplitude compensation and correction, deconvolution, noise suppression, Normal Moveout (NMO), interbedded multiple suppression processing, and prestack time migration, were performed on the field data to make sure the final prestack angle gathers meet the requirements of prestack AVA inversion.Figure 7
Stacked section from the test area. The black line indicates the location of Well A, and the arrow indicates the target reservoir.Before the inversion algorithm is implemented,γdry2 was estimated based on the statistical analysis of logging data. First, the logging curve is divided into oil layers and nonoil layers. Then, the corresponding fluid factor curve, which is calculated with specified different γdry2, is used to distinguish the divided layers. Finally, the fluid factor curve with the relatively high resolution is retained, and the corresponding value of γdry2 is determined as 2.16. From the black curve shown in Figure 8(a), it can be seen that the fluid factor can well identify the oil-bearing layer, indicating that the estimated value of γdry2 is reasonable. Next, the proposed method is applied to the 2D field data. Figure 8 shows the field inversion results for fluid factor (Figure 8(a)), shear modulus (Figure 8(b)), and density (Figure 8(c)) inverted by using the new method. Black curves indicate real logging data of Well A. Figure 9 shows the comparison of inverted results at the Well A location and real logging data in the time domain. From these figures, we note that the inversion results of fluid factor, shear modulus, and density show good agreements with the logging data and the fluid factor inversion results estimated by the proposed method can well characterize oil layer, which indicating that the proposed method is feasible and valid in application to field data.Figure 8
Field data inversion results: (a) fluid factor, (b) shear modulus, and (c) density. The black curve indicates the real logging data of Well A, and the arrow indicates the target reservoir.
(a)(b)(c)Figure 9
Comparison of inverted results at location of Well A. The black line indicates the initial model data, the red line indicates the real seismic well logging curve, and the blue line indicates the inversion results.


## 3.1. Synthetic Data Example

First, a synthetic profile with single-well logs which are obtained from a sandstone layer in an actual oil field in China is used to verify the feasibility and stability of the new method. The value ofγdry2 is set to 2.333. The parameter curves of the single-well are shown in Figure 3, where the solid lines represent the real fluid factor, shear modulus, and density curves, and the dashed lines represent the initial model data that are obtained by smoothing the true curves. Then, the seismic forward modeling is implemented by convoluting a Ricker wavelet (the main frequency is 30 Hz) with the reflection coefficients generated by the exact Zoeppritz equations. The corresponding noise-free synthetic angle gathers are shown in Figure 4(a). Figure 5(a) shows the inversion results of the proposed method using noise-free data. The new method can converge to an optimal solution in only 4 to 5 iterations. We can see that in the absence of noise, all of the fluid factor, shear modulus, and density can be inverted with a high accuracy. In order to test the stability of the new method, random noises with signal-to-noise ratio (SNR) of 2, 1, and 0.5 are added to the synthetic data. These noise data with different SNR are shown in Figure 4 (b: SNR=2, c: SNR=1, d: SNR=0.5). Figures 5(b)–5(d) show the corresponding inversion results of noise data with corresponding SNRs. We can see that the fluid factor and shear modulus can be estimated stably with a satisfactory accuracy, even when SNR=0.5. However, the estimated density results show more bias when noise exists, suggesting that the estimation of density parameters is seriously affected by noise.Figure 3
True model data (solid line) and initial model data (dashed line).Figure 4
Synthetic prestack angle gathers: (a) noise-free, (b)SNR=2, (c) SNR=1, and (d) SNR=0.5.Figure 5
Inversion results of synthetic data with different SNR. The black line indicates the initial model, the blue line indicates the inversion result, and the red line indicates the true data. (a) Noise-free, (b)SNR=2, (c) SNR=1, and (d) SNR=0.5.
(a)(b)(c)(d)In order to further demonstrate the advantage of the proposed method, the inversion based on Russell linear approximation is implemented on the synthetic data shown in Figure4(a), and the inversion results are shown in Figure 6. From the comparison between Figures 5(a) and 6, it is obvious that the inverted density has been improved, and the inversion accuracy of fluid factor, and shear modulus in the range of 0 to 0.2 s also has been improved, which indicates that the proposed method can effectively improve the accuracy of the inversion results.Figure 6
Inversion results of the traditional method based on Russell linear approximation.


## 3.2. Field Data Example

A small 2D field seismic angle gather data (effective angle range is 3°-34°) which is extracted from a work area in China is used to further demonstrate the feasibility and availability of the proposed method in field data. The stack section of the field data is shown in Figure7. In this figure, the black line denotes the location of Well A, and the black arrow indicates the location of the target reservoir. The target reservoir is a shale layer, located between two sandstone layers. A series of conventional processing procedures, such as amplitude compensation and correction, deconvolution, noise suppression, Normal Moveout (NMO), interbedded multiple suppression processing, and prestack time migration, were performed on the field data to make sure the final prestack angle gathers meet the requirements of prestack AVA inversion.Figure 7
Stacked section from the test area. The black line indicates the location of Well A, and the arrow indicates the target reservoir.Before the inversion algorithm is implemented,γdry2 was estimated based on the statistical analysis of logging data. First, the logging curve is divided into oil layers and nonoil layers. Then, the corresponding fluid factor curve, which is calculated with specified different γdry2, is used to distinguish the divided layers. Finally, the fluid factor curve with the relatively high resolution is retained, and the corresponding value of γdry2 is determined as 2.16. From the black curve shown in Figure 8(a), it can be seen that the fluid factor can well identify the oil-bearing layer, indicating that the estimated value of γdry2 is reasonable. Next, the proposed method is applied to the 2D field data. Figure 8 shows the field inversion results for fluid factor (Figure 8(a)), shear modulus (Figure 8(b)), and density (Figure 8(c)) inverted by using the new method. Black curves indicate real logging data of Well A. Figure 9 shows the comparison of inverted results at the Well A location and real logging data in the time domain. From these figures, we note that the inversion results of fluid factor, shear modulus, and density show good agreements with the logging data and the fluid factor inversion results estimated by the proposed method can well characterize oil layer, which indicating that the proposed method is feasible and valid in application to field data.Figure 8
Field data inversion results: (a) fluid factor, (b) shear modulus, and (c) density. The black curve indicates the real logging data of Well A, and the arrow indicates the target reservoir.
(a)(b)(c)Figure 9
Comparison of inverted results at location of Well A. The black line indicates the initial model data, the red line indicates the real seismic well logging curve, and the blue line indicates the inversion results.


## 4. Conclusions

The proposed method aims to reduce the nonuniqueness and uncertainty of the fluid factor inversion results. Therefore, we derive a novel quadratic approximation of the Zoeppritz equations with the chosen constants of fluid factor, shear modulus, and density and proposed a new nonlinear inversion method to achieve this purpose. In order to directly invert the target parameters instead of indirectly calculating them from the inverted reflectivity, we rewrite the expressions of reflectivity in terms of target parameters on both sides of the interface. Numerical experiments show that the novel quadratic approximation has satisfactory calculation accuracy. Then, the nonlinear inversion objective function is constructed in a Bayesian framework, and the Gauss-Newton method is used to solve this nonlinear problem. Finally, we obtain the iterative updating formula of the fluid factor, shear modulus, and density; achieve the nonlinear direct inversion of these parameters; and avoid the cumulative errors caused by trace integral operation. The synthetic data test shows that the proposed method can estimate the fluid factor information stably and reasonably, even with low SNR. The field data test shows that the proposed method can effectively identify reservoir fluids. Both of them confirm the feasibility and availability of the proposed method.Compared with the state-of-the-art approaches, the advantage of the proposed method is that it can effectively reduce the nonuniqueness and uncertainty of factor inversion results and improve the estimation accuracy of fluid factors. Since the proposed approach uses a second-order nonlinear forward equation, it will increase the instability of the inversion algorithm. Therefore, the stability of the proposed method could by enhanced, which will be studied in the further.



---

*Source: 8860119-2020-12-01.xml*