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 (), 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.

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 - and -wave velocities. Goodway et al. [2] suggested that Lamé petrophysical parameters (, Lamé modulus × density and , ) were significantly better than the - and -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 () 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 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 , , 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 - and -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 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.) [2228] 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 based on a quadratic approximation. First, a new quadratic approximation of the exact Zoeppritz equations which containing the fluid factor , 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 , , 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 , , and to improve the estimated accuracy of fluid factor.

The existing quadratic approximation of the exact Zoeppritz equations in terms of - and -wave impedances and density is given by [18, 20] where and represent - and -wave impedances, is the incidence angle of wave, is the reflection angle of PS converted wave, and is density. is the background ratio of the saturated rock.

Firstly, we establish the relationship between the wave impedances and the fluid factor where is Russell fluid factor. The expression of Russell fluid factor is shown below where , and are the - and -wave velocities and density of saturated rock, respectively, represents the square of the dry-rock velocity ratio. Russell et al. [5, 6] and Wang [30] had discussed the effect of 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 -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 , 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 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.

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 function and the prior distribution function are assumed following the Gaussian distribution and Cauchy distribution, respectively. Their expressions can be expressed as follows where represents observed data, represents the nonlinear forward operator (the new quadratic approximation), is the target parameters vector, and is the noise variance. is the size of target parameters, and is the mean vector of target parameters. , of which is a covariance matrix that contains the statistical correlations among the three target parameters. Matrix is a matrix defined as

Substituting equations (11) and (12) into Bayes’ formula [31], we have where is the posterior probability distribution. Then, we need to find the maximum value of . This problem can be converted into a problem that solves the minimum value of the following objective function where is the regularization term and 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 and are the first and second partial derivatives of the nonlinear objective function with respect to the parameter vector , respectively. Their expressions are shown below

In equations (17) and (18), the matrix is the partial derivative of the nonlinear forward model based on the new quadratic approximation with respect to (Jacobian matrix). The derivation of this matrix () is given in Appendix A. The matrices and are the first and second partial derivatives of the regularization term with respect to . 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 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: , c: , d: ). 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 . 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.

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

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

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.

Appendix

A. Derivation of the Jacobian Matrix

The nonlinear forward model can be expressed as follows where represents the angle-dependent wavelet matrix and represents the -wave reflection coefficients vector. Then, the Jacobian matrix can be written as

For a certain angle , equation (A.2) can be written as where is the wavelet matrix corresponding to the incidence angle , and where represents the reflection coefficient of th reflection interface when the incident angle is . Next, we need to calculate the partial derivatives of the reflection coefficient with respect to the target parameters on both sides of the interface. Assuming that the target parameters on both sides of the th reflection interface are . Then, we have

Assuming that

Substituting equations (A.8) and (A.9) into equation (A.7), the following expressions can be obtained: where

Finally, extending equation (A.3) to the th incidence angles situation, the Jacobian matrix will become a matrix.

B. Derivation of and

The regularization term is equation (15) is given blow

The first partial derivative of the regularization term with respect to can be written as [32]

Taking the derivative of with respect to where , we have where .

Finally, the first partial derivative of the regularization term with respect to can be expressed in the matrix form shown below

The second partial derivative of the regularization term with respect to can be written as:

Taking the second partial derivative of with respect to where and , we have where and represent the and rows of matrix , 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.

Acknowledgments

This work was financially supported by the National Natural Science Foundation of China (41904116, 41874156), the Natural Science Foundation of Hunan Province (2020JJ5168), and the Doctoral fund of Hunan University of Science and Technology (E518B2).