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

Hunan Provincial Key Laboratory of Shale Gas Resource Utilization, Hunan University of Science and Technology, Xiangtan 411201, China College of Geology and Environment, Xi’an University of Science and Technology, Xi’an 710054, China PetroChina Pipeline Compressor-Set Maintenance, Repair&Overhaul Center, Langfang 065000, China Sichuan Water Resources and Hydroelectric Investigation & Design Institute, Chengdu 610072, China


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 μρ, shear modulus × 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 × Russell fluid factor) 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][23][24][25][26][27][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.

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 of P-and S-wave impedances and density is given by [18,20] where I P = ρV P and I S = ρV S 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 = ðV P /V S Þ sat is the background ðV P /V S Þ sat ratio of the saturated rock.

Geofluids
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 where V P , V S , and ρ are the P-and S-wave velocities and density of saturated rock, respectively, γ 2 dry = ðV P /V S Þ 2 dry represents the square of the dry-rock velocity ratio. Russell et al. [5,6] and Wang [30] had discussed the effect of γ 2 dry 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 Rearranging equation (3) as Taking the derivatives of both sides of equation (4), we can obtain Combining equations (3) and (5), the following equation can be obtained For S-wave impedance, we have Taking the derivatives of both sides of equation (7), we obtain Substituting equations (6) and (8) into equation (1), a new quadratic approximation in terms of fluid factor ρf , shear modulus, and density is derived 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 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.

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 function Pðd | mÞ and the prior distribution function PðmÞ are assumed following the Gaussian distribution and Cauchy distribution, respectively. Their expressions can be expressed as follows where d represents observed data, Q represents the nonlinear forward operator (the new quadratic approximation), m = ½fluid factor ; shear modulus ; density 3N×1 is the target parameters vector, and σ d 2 is the noise variance. N is the size of target parameters, and μ is the mean vector of target parameters. Φ i = ðD i Þ T ψ −1 D i , 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 Substituting equations (11) and (12) into Bayes' formula [31], we have where Pðm | dÞ is the posterior probability distribution. Then, we need to find the maximum value of Pðm | dÞ. This problem can be converted into a problem that solves the minimum value of the following objective function where RðmÞ = 2∑ N i=1 ln ð1 + ðm − μÞ T Φ i ðm − μÞÞ is the regularization term and β = σ d 2 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 where γðmÞ and HðmÞ 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 In equations (17) and (18), the matrix ∂QðmÞ/∂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 (∂QðmÞ/∂m) is given in Appendix A. The matrices ∂RðmÞ/∂m and ∂R 2 ðmÞ/∂m 2 are the first and second partial derivatives of the regularization term RðmÞ with respect to m. Their derivations are given in Appendix B.

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 γ 2 dry 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.

Geofluids
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   9 Geofluids 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.
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 Figure 4(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.  10 Geofluids

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 Figure 7. 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.

Geofluids
Before the inversion algorithm is implemented, γ 2 dry 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 γ 2 dry , is used to distinguish the divided layers. Finally, the fluid factor curve with the relatively high resolution is retained, and the corresponding value of γ 2 dry 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 γ 2 dry 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.

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. where W represents the angle-dependent wavelet matrix and rðmÞ represents the P-wave reflection coefficients vector. Then, the Jacobian matrix can be written as where Wðθ i Þ is the wavelet matrix corresponding to the incidence angle θ i , and

ðB:6Þ
where Φ i j and Φ i k represent the j and k rows of matrix Φ i , respectively.

Data Availability
All model data can be obtained by contacting the corresponding author; researchers can replicate the analysis. But the real seismic data in the application section are not available because they involve business secrets.

Conflicts of Interest
The authors declare that they have no conflicts of interest.