Pre-Stack Nonlinear Direct Exact Inversion of Fracture Parameters in Deep Shale Reservoirs

Abstract

A conventional linear pre-stack inversion method under the conventional stationary convolution model is limited by the assumptions of weak formation contrast change and small angle incidence and fails to take into account the amplitude attenuation of seismic wave propagation. Meanwhile, the resolution and precision of oil and gas evaluation and fracture characterization of shale reservoirs under complex geological conditions are low because the compaction and non-connectivity characteristics of deep shale reservoirs are not fully considered. Therefore, porous rock pores are divided into connected pores and disconnected pores. Combined with the effect of compaction on dry rock skeleton, a petrophysical model considering the compaction and pore dysconnectivity of deep shale reservoir is developed. The quantitative relationship between transverse isotropy with a vertical axis of symmetry (VTI) stiffness matrix, rock physical properties, and fracture parameters is established in this model. It provides a more accurate scheme for the original physical modeling of deep shale. This relationship is incorporated into the exact VTI reflection coefficient equation, and a nonstationary convolution operator is derived by using the attenuation theory of seismic wave propagation. A nonstationary pre-stack nonlinear direct inversion method of fracture parameters of shale reservoirs with horizontal fractures is proposed, which Improves the resolution and accuracy of shale reservoir gas bearing and fracture characteristics prediction. It provides a new way to accurately characterize the fracture development and oil-bearing property of shale reservoirs. A model test and field data test verify the effectiveness of this method.

1. Introduction

Pre-stack seismic inversion, as an effective means to describe subsurface elastic and physical property changes through seismic amplitude information, has gradually become a hot topic in the field of seismic petrophysical reservoir characterization and fluid identification [1,2,3,4,5]. The pre-stack inversion method based on a Bayesian framework can effectively solve the invalidity problem of conventional amplitude variation with the offset (AVO) inversion by reasonably introducing prior information and improve the accuracy and stability of inversion [6,7,8].

With the increasing demand for oil and gas exploration, different types of complex reservoirs have begun to be explored and studied by researchers [9,10]. When the underground medium shows anisotropic characteristics due to the existence of fractures, the seismic petrophysical theory of conventional reservoirs is not enough to give an accurate description of underground oil and gas reservoirs. Anisotropic unconventional petrophysical and seismic inversion theories based on transverse isotropy with a vertical axis of symmetry (VTI) or transverse isotropy with a horizontal axis of symmetry (HTI) theories have gradually developed and provided guidance for reservoir identification of unconventional reservoirs [11,12,13,14,15,16]. In order to make full use of multi-wave data in anisotropic media, Hou et al. [17] proposed a multi-wave joint pre-stack inversion method based on Bayesian theory. This method uses Ruger approximation to describe the reflection coefficient in anisotropic media and uses a least square optimization algorithm to inverse longitudinal and shear wave velocity, density, and Thomsen anisotropy parameters [18]. Based on the Bayesian theory, the method assumes that the likelihood function obeys Gaussian prior distribution and the improved Cauchy distribution of a prior model. The method is tested on the noise-free and noisy synthetic multi-wave data, and the inversion results are obtained with higher accuracy, anti-noise ability, and stability than the conventional p-wave data. In view of the characteristics of azimuth–anisotropy in formations with vertical fractures, Chen et al. [19] introduced a fluid factor index from the anisotropic petrophysical model to analyze the relationship between fluid factor and fluid type, saturation, and fracture aspect ratio and then studied the seismic response characteristics of the medium under different fluid filling. An inversion method of AVAZ (amplitude variation with incident angle and azimuth) is proposed to estimate fluid factors, which can effectively identify fluid types in fractured reservoirs. It has good anti-noise performance. Based on the study of VTI anisotropic media, Wang et al. [20] proposed an improved seismic inversion method for SH-SH waves in VTI media, which deduced a new approximate formula for SH-SH wave reflection coefficient. At the same time, it has higher accuracy and contains only two independent parameters to be retrieved (shear wave impedance and horizontal SH wave velocity). Based on this equation, a VTI medium SH-SH wave inversion method is developed, which has lower uncertainty than conventional PP wave inversion and reduces the need for large angle data in PP wave anisotropy inversion. In order to obtain independent shear wave anisotropy parameters, a new method based on physical rock relationships is proposed. The new method can accurately predict shear wave impedance, horizontal SH wave velocity, and anisotropy parameters, which provides a reliable solution for lithology interpretation and reservoir prediction of anisotropic strata and expands the application potential of shear wave seismic exploration. Based on the exact VTI equation, Lang et al. [21] proposed a PP wave reflection coefficient equation for VTI media that can overcome some of the limitations of the existing linear approximation. Compared with the exact anisotropic reflection coefficient equation, the new reflection coefficient equation is replaced by relative values. After reducing the number of inversion parameters, the inversion results of five reflectance parameters are obtained by using the Bayesian inversion framework of a Markov chain Monte Carlo (McMC) algorithm [22]. This method further improves the stability of inversion and provides an improved scheme for anisotropic seismic interpretation of oil–gas shale reservoirs. Furthermore, the theoretical and methodological breakthroughs in permeability of underground reservoirs and probabilistic analysis in soils and rocks have also contributed valuable results to reservoir characterization [23,24,25].

Robinson convolution theory [26] summarizes the propagation process of seismic waves as the convolution process of seismic wavelet and reflection coefficient, which is an ideal model to describe the change in the seismic wave field with time and space. On this basis, Margrave [27] proposed a generalized linear theory to extend the convolution method to non-stationary processes. The theory can be applied to any linear, non-stationary filter with arbitrary time–frequency variation in the time domain, frequency domain, or mixed domain. The method provides an intuitive explanation for the analysis of non-stationary convolution and a new theoretical basis and implementation scheme for the processing of non-stationary signals. Later, Margrave et al. [28] proposed a Gabor deconvolution method for processing non-stationary seismic data, which uses Gabor transform to approximately decompose the non-stationary convolution model into the product of the Gabor transform of source waveform, attenuation function, and reflectivity. At the same time, the source waveform and attenuation function can be estimated directly from the Gabor spectrum of the seismic waveform. The proposed method extends the stationary pulse deconvolution to the case of non-stationary signals. At the same time, the effect of attenuation on the waveform is considered, so as to provide help for better recovery of formation reflectivity. In order to solve the problem of non-stationary propagation of seismic waves in seismic inversion, Chai et al. [29] proposed a method called non-stationary sparse reflectance inversion (NSRI), which can directly extract reflectance sequences from non-stationary seismic data and avoid the use of inverse Q filtering method with instability. It lays a foundation for the non-stationary convolution model to be used in pre-stack seismic inversion.

However, the above existing theories lack sufficient consideration of the compaction and dysconnectivity characteristics of deep shale reservoirs, and the anisotropy inversion methods of the conventional stationary convolution model based on linear approximation are limited by a variety of assumptions and fail to fully consider the attenuation effect of deep seismic wave propagation, which is difficult to meet the accurate description of deep complex geological conditions. Therefore, starting from petrophysical modeling, this study established a petrophysical model of a shale reservoir, considering deep compaction, established the quantitative relationship between VTI medium stiffness matrix and fracture parameters, and then substituted this relationship into the exact VTI reflection coefficient equation, based on the non-stationary convolution model. A non-stationary pre-stack nonlinear direct inversion method of fracture parameters is established to describe the fracture development of underground fractured shale reservoir more accurately and to lay a foundation for reservoir characterization.

2. Amplitude Variation with Angles for VTI Media

First, we aim to obtain nonlinear relationships between P-wave velocity, S-wave velocity, and density and anisotropy parameters in deep tight shale reservoirs through petrophysical modeling. The V-R-H model [30] can be used to obtain the expression of rock matrix containing quartz, clay, and kerogen mixed minerals:

(1)$K_m=\left[C{K}_c+\left(1-C-C_{\ker }\right)K_q+C_{\mathrm{k}}{K}_{\mathrm{k}}+{\displaystyle \frac{1}{C/K_c+\left(1-C-C_{\mathrm{k}}\right)/K_q+C_{\mathrm{k}}/K_{\mathrm{k}}}}\right]/2$

(2)${\mu }_m=\left[C{\mu }_c+\left(1-C-C_{\ker }\right){\mu }_q+C_{\mathrm{k}}{\mu }_{\mathrm{k}}+{\displaystyle \frac{1}{C/{\mu }_c+\left(1-C\right)/{\mu }_q+C_{\mathrm{k}}/{\mu }_{\mathrm{k}}}}\right]/2$

where $K_q\left(GPa\right)$, $K_c\left(GPa\right)$, $K_{\mathrm{k}}\left(GPa\right)$ and ${\mu }_q\left(GPa\right)$, ${\mu }_c\left(GPa\right)$, and ${\mu }_{\mathrm{k}}\left(GPa\right)$ represent the bulk modulus and shear modulus of quartz, clay, and kerogen, respectively.$C\left(GPa\right)$ and $C_{\ker }\left(GPa\right)$ represent the clay content and kerogen content, respectively. Considering the disconnected properties of deep pores due to compaction, we divide the total porosity $\varphi \left(-\right)$ into two parts. If ${\varphi }_c\left(-\right)$ is used to represent the connected porosity, then the disconnected porosity can be represented as $\varphi -{\varphi }_c$. Assuming that disconnected pores are a part of the rock matrix, the equivalent modulus of the rock matrix containing disconnected pores is obtained by using the spherical inclusion model as follows:

(3)$K_{am}={\displaystyle \frac{4K_m{\mu }_m{\varphi }_c}{3K_m\left(1-{\varphi }_c\right)+4{\mu }_m}}$

(4)${\mu }_{am}={\displaystyle \frac{{\mu }_m\left(9K_m+8{\mu }_m\right){\varphi }_c}{9K_m+8{\mu }_m+6\left(K_m+2{\mu }_m\right)\left(1-{\varphi }_c\right)}}$

Pride [31] gave the quantitative relationship between the skeleton modulus of dry rock under diagenetic consolidation and the matrix modulus of rock. Referring to Equations (3) and (4), we can establish the bulk modulus $K_{dry}\left(GPa\right)$ and shear modulus ${\mu }_{dry}\left(GPa\right)$ of dry rock skeleton under an isotropic background as follows:

(5)$K_{dry}={\displaystyle \frac{K_{am}\left(1-{\varphi }_c\right)}{1+e{\varphi }_c}}$

(6)${\mu }_{dry}={\displaystyle \frac{{\mu }_{am}\left(1-{\varphi }_c\right)}{1+1.5e{\varphi }_c}}$

where $e\left(-\right)$ represents the consolidation parameter, which is used to reflect the degree of compaction between rock particles, usually about 2–20. Considering that most shale reservoirs in the study area contain horizontally stratified fractures, we can simplify them into VTI medium. We will apply the Schoenberg linear sliding model [32] to describe the effect of fractures on the skeleton modulus of dry rock by defining two fracture parameters, ${\delta }_N$ and ${\delta }_T$, which reflect the normal and tangential weakness of fractures. The following equations are used to describe the stiffness matrix of rock skeleton in VTI medium containing horizontal fractures:

(7)$C^{dry}=\left[\begin{array}{cccccc}{M}_{dry}\left(1-{\displaystyle \frac{{{\lambda }_{dry}}^2}{{M_{dry}}^2}}{\delta }_N\right)& {\lambda }_{dry}\left(1-{\displaystyle \frac{{\lambda }_{dry}}{M_{dry}}}{\delta }_N\right)& {\lambda }_{dry}\left(1-{\delta }_N\right)& 0& 0& 0\\ {\lambda }_{dry}\left(1-{\displaystyle \frac{{\lambda }_{dry}}{M_{dry}}}{\delta }_N\right)& M_{dry}\left(1-{\displaystyle \frac{{{\lambda }_{dry}}^2}{{M_{dry}}^2}}{\delta }_N\right)& {\lambda }_{dry}\left(1-{\delta }_N\right)& 0& 0& 0\\ {\lambda }_{dry}\left(1-{\delta }_N\right)& {\lambda }_{dry}\left(1-{\delta }_N\right)& M_{dry}\left(1-{\delta }_N\right)& 0& 0& 0\\ 0& 0& 0& {\mu }_{dry}\left(1-{\delta }_T\right)& 0& 0\\ 0& 0& 0& 0& {\mu }_{dry}\left(1-{\delta }_T\right)& 0\\ 0& 0& 0& 0& 0& {\mu }_{dry}\end{array}\right]$

(8)$M_{dry}=K_{dry}+{\displaystyle \frac{4}{3}}{\mu }_{dry},$

(9)${\lambda }_{dry}=K_{dry}-{\displaystyle \frac{2}{3}}{\mu }_{dry},$

where $M_{dry}\left(GPa\right)$ and ${\lambda }_{dry}\left(GPa\right)$ represent the P-wave modulus and Lamme parameters of the isotropic dry rock background, respectively. Mavko and Bandyopadhyay [33] gave the expression of anisotropic stiffness matrix of fluid replacement equation. As a generalized form of the Gassmannn equation, this equation was associated with the stiffness matrix of VTI medium before and after fluid replacement. Based on Equations (3)–(8), it can be obtained:

(10)$C_{ijkl}^{sat}=C_{ijkl}^{dry}+{\displaystyle \frac{\left(K_{am}{\delta }_{ij}-C_{ij\alpha \alpha }^{dry}/3\right)\left(K_{am}{\delta }_{ij}-C_{\beta \beta kl}^{dry}/3\right)}{\left(K_{am}/K_f\right){\varphi }_c\left(K_{am}-K_f\right)+\left(K_{am}-C_{ppqq}^{dry}/9\right)}}$

(11)${\delta }_{ij}=\left\{\begin{array}{cc}1,& i=j\\ 0,& i\ne j\end{array}\right.$

(12)$C_{ppqq}^{dry}={\displaystyle \sum _{p=1}^3{\displaystyle \sum _{q=1}^3{C}_{ppqq}^{dry}}}$

(13)$C_{ij\alpha \alpha }^{dry}={\displaystyle \sum _{\alpha =1}^3{C}_{ij\alpha \alpha }^{dry}}$

In Equation (10), $K_f\left(GPa\right)$ represents the pore fluid bulk modulus, which can be calculated by Wood’s formula [34] as follows:

(14)$K_f={\left({\displaystyle \frac{S_w}{K_w}}+{\displaystyle \frac{1-S_w}{K_h}}\right)}^{-1}$

where $S_w\left(-\right)$ is the water saturation, and $K_w\left(GPa\right)$ and $K_h\left(GPa\right)$ are the bulk modulus of water and natural gas, respectively. Based on Equations (10)–(13), we can derive the independent stiffness parameters contained in the exact VTI equation as follows:

(15)$C_{11}^{sat}=C_{11}^{dry}+{\displaystyle \frac{{\left(K_{am}-\left(C_{11}^{dry}+C_{12}^{dry}+C_{13}^{dry}\right)/3\right)}^2}{\left(K_{am}/K_f\right){\varphi }_c\left(K_{am}-K_f\right)+\left(K_{am}-\left(2C_{11}^{dry}+C_{33}^{dry}+2C_{12}^{dry}+4C_{13}^{dry}\right)/9\right)}}$

(16)$C_{13}^{sat}=C_{13}^{dry}+{\displaystyle \frac{\left(K_{am}-\left(C_{11}^{dry}+C_{12}^{dry}+C_{13}^{dry}\right)/3\right)\left(K_{am}-\left(2C_{13}^{dry}+C_{33}^{dry}\right)/3\right)}{\left(K_{am}/K_f\right){\varphi }_c\left(K_{am}-K_f\right)+\left(K_{am}-\left(2C_{11}^{dry}+C_{33}^{dry}+2C_{12}^{dry}+4C_{13}^{dry}\right)/9\right)}}$

(17)$C_{33}^{sat}=C_{33}^{dry}+{\displaystyle \frac{{\left(K_{am}-\left(2C_{13}^{dry}+C_{33}^{dry}\right)/3\right)}^2}{\left(K_{am}/K_f\right){\varphi }_c\left(K_{am}-K_f\right)+\left(K_{am}-\left(2C_{11}^{dry}+C_{33}^{dry}+2C_{12}^{dry}+4C_{13}^{dry}\right)/9\right)}}$

(18)$C_{55}^{sat}=C_{55}^{dry}$

The saturated rock’s density ${\rho }_{sat}\left(kg/m^3\right)$ can be estimated by volumetric averaging as follows:

(19)${\rho }_{sat}=\left(C{\rho }_c+\left(1-C-C_{\ker }\right){\rho }_q+C_{\mathrm{k}}{\rho }_{\mathrm{k}}\right)\left(1-\varphi \right)+\left(S_w{\rho }_w+\left(1-S_w\right){\rho }_h\right)\varphi$

where ${\rho }_c\left(\mathrm{kg}/{\mathrm{m}}^3\right)$, ${\rho }_q\left(\mathrm{kg}/{\mathrm{m}}^3\right)$, ${\rho }_{\mathrm{k}}\left(\mathrm{kg}/{\mathrm{m}}^3\right)$, ${\rho }_w\left(\mathrm{kg}/{\mathrm{m}}^3\right)$, and ${\rho }_h\left(\mathrm{kg}/{\mathrm{m}}^3\right)$ represent the density of clay, quartz, kerogen, water, and gas, respectively.

Studying the energy distribution of waves at the interface is an important topic in the field of rock seismology, which provides a theoretical basis for understanding the structure and properties of underground rock layers and for oil and gas exploration and development. Ruger [35] gave the reflection coefficient equation in the case that the incidence angle of P-wave in VTI medium is $\theta \left(°\right)$, the expression is as follows: the expression will be shown by Equations (S1)–(S21) in the Supplementary Material.

Based on the weak anisotropy hypothesis [18], we will make an appropriate simplification of the exact VTI equation. In order to reduce the complexity of the inversion process and improve the stability of the inversion, the phase velocity formula can be simplified and expressed as follows:

(20)$V_p\left(\theta \right)=\sqrt{a_{33}}\left(1+\delta {\sin }^2\theta {\cos }^2\theta +\epsilon {\sin }^4\theta \right)$

(21)$V_s\left(\theta \right)=\sqrt{a_{55}}\left(1+{\displaystyle \frac{a_{55}}{a_{33}}}\left(\epsilon -\delta \right){\sin }^2\theta {\cos }^2\theta \right)$

(22)$\epsilon ={\displaystyle \frac{C_{11}^{sat}-C{\mathrm{sat}}_{33}}{2C_{33}^{sat}}}$

(23)$\delta ={\displaystyle \frac{{\left(C_{13}^{sat}+C_{55}^{sat}\right)}^2-{\left(C_{33}^{sat}-C_{55}^{sat}\right)}^2}{2C_{33}^{sat}\left(C_{33}^{sat}-C_{55}^{sat}\right)}}$

Among them, $\epsilon \left(-\right)$ and $\delta \left(-\right)$ are other Thomsen anisotropy parameters that characterize crack development states. Substitute Equations (1)–(19) into the exact reflection coefficient equation of VTI shown in Equations (S1)–(S21) and Equations (20)–(23) for consideration. $C\left(-\right)$, $\varphi \left(-\right)$, $S_w\left(-\right)$, ${\delta }_N\left(-\right)$, ${\delta }_T\left(-\right)$ are variables that vary with the stratum, while the other parameters do not change with the stratum, and ${\varphi }_c/\varphi$ are constant; we yield the following:

(24)$R_{pp}\left(\theta \right)=f\left(\begin{array}{cccccccccc}{C}^{\left(1\right)}& C^{\left(2\right)}& {\varphi }^{\left(1\right)}& {\varphi }^{\left(2\right)}& S_w^{\left(1\right)}& S_w^{\left(2\right)}& {\delta }_N^{\left(1\right)}& {\delta }_N^{\left(2\right)}& {\delta }_T^{\left(1\right)}& {\delta }_T^{\left(2\right)}\end{array}\right)$

Equation (24) is the derived nonlinear VTI reflection coefficient equation. As the core equation of inversion, it describes the nonlinear relationship between the reflection coefficient $R_{pp}\left(\theta \right)$ and five independent petrophysical parameters ($C\left(-\right)$, $\varphi \left(-\right)$, $S_w\left(-\right)$, ${\delta }_N\left(-\right)$, and ${\delta }_T\left(-\right)$). The acquisition of $R_{pp}^{VTI}$ can be obtained by Klem’s rule or matrix inversion. Scholars have given detailed expressions (such as: Graebner’s research achievement in 1992 [36]), which will not be repeated here.

In order to verify the accuracy of the proposed equation, we designed the parameters of the double-layer shale model as shown in Table 1 to test the accuracy of Equation (24). In Table 1, $C_{33}^{sat}\left(GP\mathrm{a}\right)$ and $C_{55}^{sat}\left(GP\mathrm{a}\right)$ are the two stiffness matrix components of saturated rock, respectively. Figure 1a,b, respectively, show the accuracy comparison of multiple reflection coefficient equations and the absolute values of errors compared to the exact VTI reflection coefficient equations. Compared with linear approximation, the accuracy of Equation (24) is obviously improved, which lays a foundation for the new method to invert deep fracture parameters more accurately.

3. Inversion Theory

Considering that deep seismic wave propagation has strong attenuation characteristics, based on the nonstationary convolution theory, we construct a nonstationary pre-stack nonlinear direct inversion method by deducing the nonlinear pre-stack seismic forward operator under the nonstationary convolution model.

Following the derivation process shown in Equations (S22) and (S33) in the Supplementary Data, the objective function of the final pre-stack direct inversion is as follows:

(25)$F(m)={\displaystyle \sum _{i=1}^p{\left({\mathbf{S}}_p\prime -{\mathbf{G}}_p\prime \left(\mathbf{m}\prime \right)\right)}^{{\rm T}}\left({\mathbf{S}}_p\prime -{\mathbf{G}}_p\prime \left(\mathbf{m}\prime \right)\right)}+2{\sigma }_n^2{\displaystyle \sum _{i=1}^N\ln (1+m{\prime }_i^2/{\sigma }_m^2)}$

(26)$\Lambda ={({\mathbf{\eta }}_i-l_i{\mathbf{m}}_i)}^{{\rm T}}({\eta }_i-l_i{m}_i)$

(27)${\eta }_i=1/2\times \ln (m_i/m_{i0})$

(28)$l_i={\displaystyle {\int }_{t_0}^{t_i}d\tau }$

The McMC nonlinear optimization algorithm [22] is used to solve the nonlinear objective function shown in Equation (25), and the inversion results of shale reservoir fracture parameters can be obtained while the minimum value of the objective function is obtained. The nonlinear optimization solution process based on the McMC algorithm is described in Figure S1 of the Supplementary Materials.

4. Synthetic and Field Data Examples

The synthetic examples are used to evaluate the rationality and feasibility of the proposed method. Based on the nonstationary convolution model, we use a real well interval in a working area for forward synthesis recording. Figure 2 shows the curves of the petrophysical parameter models used for the synthetic record test. The vertical coordinate in Figure 2 represents the time with a sampling interval of 2 ms, and the horizontal coordinate represents the range of each petrophysical parameter. Figure 3a,b show the synthetic seismogram with noise free and signal-to-noise ratio (SNR) [37] of five, respectively. The SNR is the ratio of the noise power to the original noise-free power, which is used to characterize the intensity of the noise in the signal, and the noise intensity in the signal increases with the reduction in the SNR. The ordinate in Figure 3 represents the time, the sampling interval is 2 ms, and the abscissa represents the incidence angle. Based on the noiseless synthesis record shown in Figure 3a, Figure 4 shows the inversion results obtained by using the novel nonstationary pre-stack direct nonlinear inversion method designed by us and the conventional stationary deconvolution nonlinear inversion method, respectively. It can be seen that the inversion results of the new method have a high agreement with the real petrophysical models, which proves that the proposed method can be applied to accurately estimate the fracture parameters of underground shale reservoirs. Compared with the conventional nonlinear inversion method from the stationary convolution model, the new method can effectively improve the accuracy of deep inversion results and protect the reflection characteristics of the formation by considering the attenuation of seismic waves, which provides a theoretical basis for the effective application of the method in deep layers. When the synthetic record contains noise with an SNR of 5, the inversion results of the new method and the conventional nonlinear method are, respectively, shown in Figure 5. It can be seen that the new method can provide more accurate prediction results for the five petrophysical parameters to be retrieved. Meanwhile, due to the use of the non-stationary inversion method, the new method is still superior to the conventional inversion method in terms of the coincidence of inversion results, although the inversion accuracy is somewhat lower than that of the noiseless inversion results. However, it can still provide reliable prediction results for fracture parameters of shale reservoirs. The synthetic seismogram test results fully illustrate the rationality and reliability of the proposed method and provide theoretical support for its application in the practical work area.

Field data examples will be further used to evaluate the application ability of the proposed method in the actual working area. The figures of the field data examples are shown in the Supplementary Materials. Figure S2 shows our fully stacked cross-well seismic profile with actual seismic data from an actual shale gas field in southwest China. The horizontal coordinate in Figure S2 represents the number of CDP sampling points in the profile, including 230 CDP sampling points; the vertical coordinate represents the time range of 2.1–2.45 s, with a time sampling interval of 2 ms. Several sets of fractured gas shale reservoirs have developed in the actual section, and the reservoirs are buried at the travel time about 2.1–2.3 s, located at the positions indicated by the white arrows in the figures. Well A at the position of the black vertical line is a blind well used to verify the inversion results. The seismic data are processed by the contractor, and the partial superimposed seismic section involved in the final inversion is the processed pure wave section.

Figure S3 shows the results of the petrophysical parameters predicted by nonstationary pre-stack nonlinear inversion. It can be seen that the inversion results have a high coincidence with the logging curve, the pre-stack physical parameters can better reflect the property changes in formation porosity, shale content, and water saturation. The two fracture parameters ${\delta }_N$ and ${\delta }_T$ can effectively indicate the fracture weakness property of the reservoir. It can clearly describe the fracture development degree of the shale gas reservoir (the greater the absolute value of fracture parameters, the higher the fracture development degree) and the distribution state of fractured formation, which provides a reliable method for evaluating the fracture and gas bearing of shale gas reservoir. Figure S4 shows the inversion results of the nonlinear pre-stack method using the conventional stationary convolution model. Compared with the results obtained via the new method, it can be seen that the formation where the gas-bearing fractured reservoir develops, indicated by the white arrow in the figures, although it can predict the petrophysical parameters more reliably than the logging results, has low accuracy in predicting the physical properties of shale gas reservoir, fracture parameters, and formation distribution characterization. In addition, at the position circled by the white dotted line in the figures, the nonlinear method of the conventional convolution model is not ideal for the characterization of the thin layer with weak reflection characteristics, and the resolution of the inversion results is low. The comparison between the methods highlights the advantages of the new method in the accuracy and resolution of petrophysical parameters prediction of deep shale gas reservoirs and lays a foundation for its application in practical production to improve the accuracy and resolution of reservoir characterization.

5. Discussion

Aiming at deep shale reservoirs with horizontal fractures, a nonlinear inversion method is proposed to estimate fracture parameters and oil- and gas-bearing indicator parameters of deep shale reservoirs. The method provides a reference for fracture development, distribution characterization, and oil- and gas-bearing characterization of shale gas reservoirs. Therefore, it is necessary to discuss the innovation, shortcomings, and application scope of this method.

5.1. Innovation

The innovation of this method is that it targets deep shale reservoirs, taking into account deep pore compaction and the pore’s partial connectivity [31,38]. Through petrophysical modeling, the quantitative relationship between the stiffness matrix of VTI medium and fracture parameters used to evaluate the degree of reservoir fracture development is obtained, and this relationship is then incorporated into the exact VTI nonlinear reflection coefficient equation. The nonlinear reflection coefficient equations of VTI media with fracture parameters are obtained. Compared with the linear approximation of the conventional VTI medium under weak anisotropy hypothesis [39], this equation not only improves the accuracy of the VTI reflection coefficient under the weak anisotropy hypothesis but also avoids the assumptions of formation reflection coefficient change and small angle incidence inherent in the linear inversion method [40] in the actual inversion, and its applicability in complex geological conditions is stronger. It is helpful to meet the need of high-precision reservoir parameter prediction. On the other hand, we consider the attenuation effect of deep seismic wave propagation in the actual inversion; combined with the stationary convolution model [26], we design a set of non-stationary pre-stack nonlinear direct inversion method for direct inversion prediction of fracture parameters in VTI media. The nonlinear inversion method based on the non-stationary convolution model can effectively protect the weak reflection characteristics of the deep shale reservoir and enhance the resolution of petrophysical parameter prediction in the deep shale reservoir compared with the nonlinear inversion method based on the conventional stationary convolution model. In addition, the petrophysical modeling part of this study complements the existing anisotropic petrophysical theory of shale under deep compaction and provides theoretical guidance for the prediction of deep shale modulus and velocity. The VTI dielectric reflection coefficient equation in this study also provides a clear explanation for the study of the deep consolidation effect and the variation in the reflection coefficient with physical rock parameters and fracture parameters of shale reservoirs with low pore connectivity and lays a foundation for seismic forward response analysis of deep shale oil and gas reservoirs.

5.2. Application Scope and Shortcomings

The method proposed in this study also has disadvantages and limitations. Its disadvantage is that the nonlinear equation has a more complex mathematical form, which increases the inversion time cost compared with the simple linear inversion. One limitation is that due to the simplification of the phase velocity formula of the nonlinear equation, it is still limited by the weak anisotropy hypothesis [18], which makes its applicability in the strong anisotropy working area weak. The nonlinear inversion method proposed in this study is mainly targeted at shale-working areas with horizontal fractures. Since pore disconnectedness and compaction are considered in the petrophysical modeling process, the attenuation of seismic wave propagation is considered in the inversion process, and the nonlinear equations with fewer assumptions are adopted. This method can be applied to deep shale reservoirs or shale reservoirs with horizontal fractures in complex geological environments and has a wide range of applications.

6. Conclusions

In this paper, a novel nonlinear direct inversion method for deep shale fracture parameters is proposed, which solves the main limitations of traditional linear inversion techniques. By integrating petrophysical models that accurately account for pressure effects and disconnected pores, the method improves seismic wave modeling and better characterizes fracture development under complex geological conditions. Then, we establish a quantitative relationship between VTI stiffness matrix and fracture parameters, which is incorporated into the VTI reflection coefficient equation, and the inversion results are more accurate. Compared with traditional methods, this method has higher resolution and accuracy in predicting gas-bearing properties and fracture characteristics. The model and field data tests verify the effectiveness of the method in shale reservoirs with horizontal fractures, which provides a promising tool for shale oil and gas exploration and development. But it is difficult to have good applicability in shale reservoirs with high angle fractures, so improving the application ability of petrophysical model and inversion method in shale reservoirs with high angle fractures is a potential research direction in the future. Furthermore, the implementation of the nonlinear inversion algorithm is more complicated, and the time cost is high. How to improve the inversion efficiency is also the main aim to the next research.

Author Contributions

Program writing, M.W.; writing—original draft preparation, L.Y.; writing—review and editing, T.G.; method applications research, X.S.; implementation of petrophysical modeling, X.Z.; reflection coefficient equation derivation, Y.R. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

Data is unavailable due to privacy or ethical restrictions.

Acknowledgments

This research was supported by Study on metallogenic mechanism of Zichang-Yanchuan salt bed and its role in natural gas accumulation.

Conflicts of Interest

Authors Meng Wang, Liang Yu and Tianchao Guo were employed by the company Gasfield Company of Shaanxi Yanchang Petroleum Group Co., Ltd. Author Xiuyan Song was employed by the company Beijing Precise Energy Technology (PET). The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Table

Figure 1. Precision analysis of reflection coefficient equations: (a) precision comparison; (b) absolute error between the equations and the exact VTI equation.

Figure 2. Petrophysical parameter models for synthetic examples. (a)[Forumla omitted. See PDF.]; (b) [Forumla omitted. See PDF.]; (c) [Forumla omitted. See PDF.]; (d) [Forumla omitted. See PDF.]; (e) [Forumla omitted. See PDF.].

Figure 3. Synthetic seismogram: (a) noise-free; (b) SNR = 5.

Figure 4. Pre-stack nonlinear inversion results under different methods with synthetic seismogram at noise-free: (a)[Forumla omitted. See PDF.]; (b) [Forumla omitted. See PDF.]; (c) [Forumla omitted. See PDF.]; (d) [Forumla omitted. See PDF.]; (e) [Forumla omitted. See PDF.]

Figure 5. Pre-stack nonlinear inversion results under different methods with synthetic seismogram at SNR = 5 (a)[Forumla omitted. See PDF.]; (b) [Forumla omitted. See PDF.]; (c) [Forumla omitted. See PDF.]; (d) [Forumla omitted. See PDF.]; (e) [Forumla omitted. See PDF.].

Table 1

The double-layer shale model.

	$\mathit{C}$ (-)	$\mathit{\varphi }$ (-)	${\mathit{\varphi }}_{\mathit{c}}$(-)	${\mathit{S}}_{\mathit{w}}$(-)	${\mathit{C}}_{\mathbf{k}}$(-)	${\mathit{\delta }}_{\mathit{N}}$(-)	${\mathit{\delta }}_{\mathit{T}}$(-)	${\mathit{C}}_{33}^{\mathit{s}\mathit{a}\mathit{t}}$(GPa)	${\mathit{C}}_{55}^{\mathit{s}\mathit{a}\mathit{t}}$(GPa)	${\mathit{\rho }}_{\mathit{s}\mathit{a}\mathit{t}}$ (kg/m3)
Gas shale	0.8	0.1	0.08	0.1	0.02	0.2	0.4	21.53	4.61	2.26
Water shale	0.8	0.05	0.04	0.9	0.01	0.1	0.2	29.51	7.45	2.41

Supplementary Materials

The following supporting information can be downloaded at https://www.mdpi.com/article/10.3390/pr13020426/s1: Equations (S1)–(S21); Equations (S22)–(S33); Figure S1. Solution flow of McMC algorithm nonlinear optimization; Figure S2: Post-stack seismic profiles for field data examples.; Figure S3. Pre-stack nonlinear direct inversion results of field data under novel nonstationary convolution model. (a) $C$; (b)$\varphi$; (c) $S_w$; (d) ${\delta }_N$; (e) ${\delta }_T$; Figure S4. Pre-stack nonlinear direct inversion results of field data under conventional stationary convolution model. (a) $C$; (b) $\varphi$; (c) $S_w$; (d) ${\delta }_N$; (e) ${\delta }_T$. Table S1. Abbreviations and full names.

References

1. Lu, J.; Yang, Z.; Wang, Y.; Shi, Y. Joint PP and PS AVA seismic inversion using exact Zoeppritz equations. Geophysics; 2015; 80, pp. R239-R250. [DOI: https://dx.doi.org/10.1190/geo2014-0490.1]

2. Yin, X.; Cheng, G.; Zong, Z. Non-linear AVO inversion based on a novel exact PP reflection coefficient. J. Appl. Geophys.; 2018; 159, pp. 408-417. [DOI: https://dx.doi.org/10.1016/j.jappgeo.2018.09.019]

3. Li, K.; Yin, X. Seismic AVO statistical inversion incorporating poroelasticity. Pet. Sci.; 2020; pp. 250-254.

4. Yu, H.; Wang, J.; Lei, Z.; Kuang, Z.; Deng, W.; Lu, L. Characterization of Natural Gas Hydrate Constrained by Well and Seismic Data in Qiongdongnan Basin. Energies; 2024; 17, 5411. [DOI: https://dx.doi.org/10.3390/en17215411]

5. Tao, S.; Guo, Y.; Huang, H.; Li, J.; Chen, L.; Gui, J.; Zhao, G. Deep-Learning-Based Amplitude Variation with Angle Inversion with Multi-Input Neural Networks. Processes; 2024; 12, 2259. [DOI: https://dx.doi.org/10.3390/pr12102259]

6. Buland, A.; More, H. Bayesian linearized AVO inversion. Geophysics; 2003; 68, pp. 184-198. [DOI: https://dx.doi.org/10.1190/1.1543206]

7. Zhang, S.; Yin, X.; Zhang, F. Prestack three term inversion method based on Trivariate Cauchy distribution prior constraint. Oil Geophys. Prospect.; 2011; 46, pp. 737-743.

8. Zhang, F.; Wei, F.; Wang, Y.; Wang, W.; Li, Y. Generalized linear AVO inversion with the priori constraint of trivariate cauchy distribution based on Zoeppritz equation. Chin. J. Geophys.; 2013; 56, pp. 2098-2115.

9. Han, S.; Bai, S.; Tang, Z.; Rui, Y.; Gong, D.; Zhang, J. Nitrogen-rich gas shale logging evaluation and differential gas-bearing characterization of lower Cambrian formation in northern Guizhou, south China. Mar. Pet. Geol.; 2020; 115, 104270. [DOI: https://dx.doi.org/10.1016/j.marpetgeo.2020.104270]

10. Han, S.; Xie, L.; Du, X.; Xiang, C.; Huang, J.; Tang, Z.; Wang, C.; Horsfield, B.; Mahlstedt, N. Insights into organic metagenesis using Raman spectroscopy and high resolution mass spectrometry: A case study of the Shahezi formation, deep Songliao basin, China. Int. J. Coal Geol.; 2023; 265, 104153. [DOI: https://dx.doi.org/10.1016/j.coal.2022.104153]

11. Zeng, L. Fissure and its seepace characteristics in low-permeable sandstone resovior. Chin. J. Geol.; 2004; 39, pp. 11-17.

12. Sa, L.; Yao, F.; Di, B.; Yao, Y. Theory and Method of Seismic Identification of Fracture-Vugg-Type Reservoirs; Petroleum Industry Press: Beijing, China, 2010.

13. Liu, R.; Zhao, J.; Yin, X. Strategy of anisotropic parameter tomography inversion in VTI medium. Oil Geophys. Prospect.; 2017; 52, pp. 484-490.

14. Yang, Y.; Yin, X.; Zhang, B.; Cao, D.; Gao, G. Linearized Frequency-Dependent Reflection Coefficient and Attenuated Anisotropic Characteristics of Q-VTI Model. Energies; 2021; 14, 8506. [DOI: https://dx.doi.org/10.3390/en14248506]

15. Pan, F.; Li, S.; Qin, D.; Chen, K.; Gao, J.; Zhang, G. Direct inversion method for fluid factors and anisotropy parameters in VTI media. Oil Geophys. Prospect.; 2024; 59, pp. 875-886.

16. Xie, Q.; Wu, Y.; Huang, Q.; Hu, Y.; Hu, X.; Guo, X.; Jia, D.; Wu, B. Prediction of Marine Thin Shale Gas Reservoir with Seismic Phase-Controlled Nonlinear Stochastic Inversion. Processes; 2023; 11, 2301. [DOI: https://dx.doi.org/10.3390/pr11082301]

17. Hou, D.; Liu, Y.; Ren, Z.; Wei, X.C.; Chen, T. Multi-wave prestack joint inversion in VTI media based on Bayesian theory. Geophys. Prospect. Pet.; 2014; 53, pp. 294-303.

18. Thomsen, L. Weak elastic anisotropy. Geophysics; 1986; 51, pp. 1954-1966. [DOI: https://dx.doi.org/10.1190/1.1442051]

19. Chen, H.; Yin, X.; Gao, C. AVAZ inversion for fluid factor based on fracture anisotropic rock physics theory. Chin. J. Geophys.; 2014; 57, pp. 968-978.

20. Wang, B.; Zhang, F.; Dai, F.; Li, X. Study on seismic inversion method of SH-SH wave in VTI media. Chin. J. Geophys.; 2023; 66, pp. 2112-2122.

21. Lang, K.; Yin, X.; Zong, Z.; Qin, D. Anisotropic Nonlinear Inversion Based on a Novel PP Wave Reflection Coefficient for VTI Media. IEEE Trans. Geosci. Remote Sens.; 2023; 61, 5913613. [DOI: https://dx.doi.org/10.1109/TGRS.2023.3295800]

22. Sambridge, M.; Mosegaard, K. Monte Carlo methods in geophysical inverse problems. Rev. Geophys.; 2002; 40, 1009. [DOI: https://dx.doi.org/10.1029/2000RG000089]

23. Khormali, A.; Ahmadi, S.; Kazemzadeh, Y. Inhibition of Barium Sulfate Precipitation During Water Injection into Oil Reservoirs Using Various Scale Inhibitors. Arab. J. Sci. Eng.; 2023; 48, pp. 9383-9399. [DOI: https://dx.doi.org/10.1007/s13369-022-07503-z]

24. Ziegler, M.O. Rock Properties and Modelled Stress State Uncertainties: A Study of Variability and Dependence. Rock Mech Rock Eng; 2022; 55, pp. 4549-4564. [DOI: https://dx.doi.org/10.1007/s00603-022-02879-8]

25. Savvides, A.A.; Antoniou, A.A.; Papadopoulos, L.; Monia, A.; Kofina, K. An estimation of clayey-oriented rock mass material properties, sited in Koropi, Athens, Greece, through feed-forward neural networks. Geotechnics; 2023; 3, pp. 975-988. [DOI: https://dx.doi.org/10.3390/geotechnics3040052]

26. Robinson, E.A. Predictive decomposition of time series with application to seismic exploration. Geophysics; 1967; 32, pp. 418-484. [DOI: https://dx.doi.org/10.1190/1.1439873]

27. Margrave, G. Theory of nonstationary linear filtering in the Fourier domain with application to time-variant filtering. Geophysics; 1998; 63, pp. 244-259. [DOI: https://dx.doi.org/10.1190/1.1444318]

28. Margrave, G.; Lamoureux, M.P.; Henley, D.C. Gabor deconvolution: Estimating reflectivity by nonstationary deconvolution of seismic data. Geophysics; 2011; 76, pp. 15-30. [DOI: https://dx.doi.org/10.1190/1.3560167]

29. Chai, X.; Wang, S.; Wei, J.; Li, J.; Yin, H. Reflectivity inversion for attenuated seismic data: Physical modeling and field data experiments. Geophysics; 2015; 81, pp. T11-T24. [DOI: https://dx.doi.org/10.1190/geo2015-0250.1]

30. Mavko, G.; Mukerji, T.; Dvorkin, J. The Rock Physics Handbook; Cambridge University Press: Cambridge, UK, 2009.

31. Pride, S.; Berryman, J.; Harris, J. Seismic attenuation due to wave-induced flow. J. Geophys. Res.; 2004; 109, B01201. [DOI: https://dx.doi.org/10.1029/2003JB002639]

32. Schoenberg, M.; Sayers, C. Seismic anisotropy of fractured rock. Geophysics; 1995; 60, pp. 204-211. [DOI: https://dx.doi.org/10.1190/1.1443748]

33. Mavko, G.; Bandyopadhyay, K. Approximate fluid substitution for vertical velocities in weakly anisotropic VTI rocks. Geophysics; 2009; 74, pp. D1-D6. [DOI: https://dx.doi.org/10.1190/1.3026552]

34. Wood, A. A Textbook of Sound; McMillan, Co.: New York, NY, USA, 1955.

35. Rüger, A. Reflection Coefficients and Azimuthal AVO Analysis in Anisotropic Media; Society of Exploration Geophysicists: Tulsa, OK, USA, 2002.

36. Graebner, M. Plane-wave reflection and transmission coefficients for a transversely isotropic solid. Geophysics; 1992; 57, pp. 1512-1519. [DOI: https://dx.doi.org/10.1190/1.1443219]

37. Bennett, W.R. Spectra of Quantized Signals. Bell Syst. Tech. J.; 1948; 27, pp. 446-471. [DOI: https://dx.doi.org/10.1002/j.1538-7305.1948.tb01340.x]

38. Liu, Q.; Yin, X.; Li, C. Rock elastic modulus estimation for tight sandstone reservoirs with disconnected pores. Geophys. Prospect. Pet.; 2015; 54, pp. 635-642.

39. Lu, J.; Wang, Y.; Chen, J.; An, Y. Joint anisotropic amplitude variation with offset inversion of PP and PS seismic data. Geophysics; 2018; 83, pp. N31-N50. [DOI: https://dx.doi.org/10.1190/geo2016-0516.1]

40. Yin, X.; Wang, H.; Cao, D.; Zhou, Q.; Guo, S. Three term AVO approximation of Kf-Km-ρ and prestack seismic inversion for deep reservoirs. Oil Geophys. Prospect.; 2018; 53, pp. 129-135.