Bayesian fluid prediction by decoupling both pore structure parameter and porosity

Carbonate reservoirs exhibit complex pore structure, which significantly affects the elastic properties and seismic response, as well as the prediction of physical parameters. As one of the main factors impacting fluid prediction, pore structure parameter directly involves in few inversion methods. In order to directly predict pore structure parameter in inversion, a novel quantitative reflection coefficient formula is proposed, that integrate Russell's poroelasticity theory with Sun's petrophysical model. This formula separates fluid bulk modulus from porosity and pore structure parameter, allowing for accurate determination of pore-fluid distribution through Bayesian framework. Both theoretical model analysis and multi-component digital core experiments of carbonates validate the importance of pore structure parameter on fluid identification. The practical application of carbonate reservoirs in Sichuan Basin demonstrates that the proposed fluid factor, eliminating the prediction illusion caused by heterogeneity in porosity and pore structure parameter within strata, provides more precise and reliable predictions compared to the Russell fluid factor. Furthermore, the similarity between the Russell fluid factor obtained directly from the Russell approximation and the Russell fluid factor calculated indirectly from the proposed method confirms the stability and accuracy of the new reflection coefficient formula.

1 Introduction

With the advancement of exploration and development technology, the targets of exploration and development are becoming increasingly intricate. There has been a decline in reserves and grades of both oil and gas, with a growing proportion of low porosity and low permeability. In terms of carbonate reservoirs, predicting reservoirs is much challenging due to their deep burial depth and weak seismic responses. Moreover, they exhibit strong heterogeneity and significant variations in physical properties and thickness over short distances, as well as complex pore types. All these characteristics contribute to heightened difficulty in identifying oil and gas within them. Therefore, one of the current challenges in geophysics lies in effectively identifying fluids within target reservoirs with complex porous media.

The effect of porosity and pore types on the effective elastic properties of reservoirs is very critical in hydrocarbon prediction. A study conducted by Sun et al. (1997) showed that the variations in the pore aspect ratio can lead to changes in wave velocity exceeding 2000 m/s or more. Currently, significant advancements have been achieved in the research of pore structure parameters. For example, Zimmerman (1986) and Kachanov et al. (1994) examined the influence of two- and three-dimensional pore shape on elastic properties, respectively. Berryman (1999) derived approximate analytical expressions for the elastic parameters of dry/saturated fractured rock based on DEM theory. Jiang et al. (2012) proposed a new petrophysical modeling method using the Gassmann equation and Eshelby–Walsh ellipsoidal fracture theory to extract a parameter characterizing variations in the pore structure, thereby demonstrating its profound influence on the elastic properties of the rock. Deng et al. (2015) provided an analytical expression for the initial minimum aspect ratio of soft pores based on poroelasticity theory and extended the “squirt flow” model based on the characteristic aspect ratio by adding soft pores iteratively to analyze the influence of complex pore distribution on the squirt flow and its possible velocity dispersion characteristics. Wang et al. (2016) established a pore-scale numerical simulation method for elastic wave propagation in porous media based on the extraction of the pore structure parameter from digital core images in view of the complex pore structure and significant heterogeneity of carbonate rocks, verified its reliability with natural core data, and quantitatively analyzed the effect of carbonate pore structures on the propagation velocity and the scattering attenuation of elastic waves. He et al. (2012) proposed a unified expression by consolidating various petrophysical models for porous media, including Pride, Geertsma, and Keys-Xu, aiming to establish a more universal application across different rock types. Subsequently, He et al. (2018) further emphasized the significance of the pore aspect ratio as a fundamental benchmark due to the ambiguous physical definition and inconsistent quantitative representation of the pore structure parameter. These studies collectively indicate that the current research on the pore structure parameter primarily focuses on small-scale cores and logs but lacks applicability for fluid prediction at the seismic scale.

In the past few decades, there has been significant advancement in the fluid identification technology of seismic reservoirs, from a “qualitative” approach based on seismic amplitude anomaly in the 1980s to a “quantitative” method relying on the fluid factor in the present stage. With the development of a prestack seismic inversion, the fluid sensitive terms have evolved from assessing the relative variation in elastic parameters to presently recognizing physical parameters with distinct petrophysical significance. Biot (1941) and Gassmann (1951) both proposed the construction methods of the fluid factor for porous fluid-saturated rocks. Russell et al. (2003) and Russell et al. (2011) used the Biot–Gassmann theory to refine the Aki–Richards approximation under saturated fluid conditions and introduced the fluid factor, which can be directly involved in seismic fluid detection. Yin et al., (2014) derived a seismic reflection coefficient formula that incorporates porosity by combining the Russell approximation with the Nur model. Du et al. (2019) then extended this formula to PP–PS joint inversions, which further improves the fluid prediction in heterogeneous reservoirs but does not consider the pore structure parameter. Zong et al. (2012) and Zong et al. (2015) established a direct relationship among fluid factors and P- and S-wave moduli based on the petrophysical model of porous elastic media to obtain a new Zoeppritz approximation formula based on the P- and S-wave moduli, which circumvents the issue related to accurately determining density fluid factor calculation and has been successfully applied in an exploration area in eastern China. Sun, et al. (2015) proposed a novel fracture fluid factor that can simultaneously detect fracture development and fluid properties by combining P-wave anisotropic fracture prediction and Russell fluid factor into the Cartesian coordinate system in order to address the challenge of fluid identification in anisotropy, which achieved promising application in igneous areas. Sun et al. (2016) employed sequential Gaussian simulation and Metropolis sampling algorithm based on Bayesian’s theoretical framework to directly estimate the Russell fluid factor, which enhanced the accuracy of fluid factor identification. Although substantial developments have been made in recent years regarding the algorithms of fluid prediction, the issue of fluid prediction under non-homogeneous conditions still faces challenges such as strong multi-solution and less precise prediction. This can be attributed to the omission of pore structure parameters in the widely used Gassmann fluid equation for bidirectional media containing fluid, making it complicated and challenging to accurately quantify the impact of pore structure. Fan et al. (2019) and Zong et al. (2019) introduced the squirt flow model into an inversion, along with the relevant parameters of the pore structure, making an initial step toward integrating the pore structure parameter in seismic reservoir prediction. Li et al. (2021) also proposed a new method combining the pore structure with the Bayesian non-linear simultaneous inversion of physical parameters, further improving the importance and involvement of the pore structure in the reservoir prediction. However, the aforementioned methods including pore structures can only be obtained in an indirect way. In other words, traditional physical parameters are first calculated, and then, specific petrophysical theories are employed to derive the pore structure and other fluid parameters. This approach does not directly incorporate the pore structure parameter as the factor in the inversion process of seismic reflection coefficient. Therefore, it is imperative to further study the inversion method of the pore structure from indirect to direct inversions in order to streamline the procedure, while improving the inversion accuracy.

In order to mitigate the impact of complex pore structure on fluid prediction in carbonate reservoirs, a novel inversion method is proposed that integrates Sun’s petrophysical model and Russell’s poroelasticity theory to derive a decoupled reflection coefficient formula for fluid bulk modulus $K_f$, porosity $\varphi$, and pore structure parameter $\gamma$. A two-layer theoretical geological model is used for forward modeling simulation to verify the accuracy of the newly derived formula and assess the impact of porosity $\varphi$ and the pore structure parameter $\gamma$ on the seismic response. Then, an improved decorrelation method is employed to eliminate correlations between multiple parameters, ensuring stability and robustness during the inversion process. In the practical application in a specific exploration area within the Sichuan Basin, this study first demonstrates through digital core experiments that the pore structure parameter $\gamma$ has great influence on elastic moduli. Furthermore, the results based on the actual data show that the fluid bulk modulus $K_f$ as the proposed fluid indicator eliminates the interference of the pore structure ($\varphi$ and $\gamma$) and can provide more accurate and reliable fluid distribution than Russell’s fluid factor $f$.

2 Methods

Establishing a direct mathematical relationship between theoretical or empirical petrophysical models and expressions of seismic reflection coefficients is the key to quantitative characterization of fluid prediction involving the pore structure.

2.1 Poroelasticity theory

First, one of the widely used quantitative formulas of seismic response for fluid identification is Biot–Gassmann theory (Krief et al., 1990), which elucidates the pore/fluid interaction in homogeneous porous media and defines the functional relationship between seismic velocity and fluid. Its corresponding formula is shown as follows:

$K_{sat}=K_{dry}+{\beta }^{2}M,(1)$

${\mu }_{sat}={\mu }_{dry}=\mu ,(2)$

where $K_{sat}$ and $K_{dry}$ represent the bulk moduli for the saturated and dry porous rocks, ${\mu }_{sat}$ and ${\mu }_{dry}$ denote the shear moduli for the saturated and dry porous rocks, $\beta$ represents the Biot coefficient, and $M$ signifies the modulus of pressure that drives water into strata without changing its volume. Russell et al. (2003) derived the Gassmann formula based on the pore-elasticity theory and found that the fluid/porosity term $f$ was related to multiple petrophysical parameters, including porosity, which could be quantitatively characterized by saturated elastic parameters:

$f={\beta }^{2}M=\frac{{\left[1-\left(K_{dry}/K_s\right)\right]}^2}{\left\{\left(\varphi /K_f\right)+\left[\left(1-\varphi \right)/K_s\right]-K_{dry}/{K_s}^2\right\}},(3)$

where $K_f$ and $K_s$ denote the bulk moduli of the fluid and saturated rock frame and $\varphi$ represents the porosity. Therefore, formula (2) can be reformulated as follows:

$K_{sat}=K_{dry}+f.(4)$

In addition, by substituting the fluid factor $f$ into Aki–Richard approximation, Russell proposed a formula for the seismic response that can directly obtain $f$ by inversion (Russell et al., 2011) as follows:

$R_{pp}\left(\theta \right)=\left[\left(1-\frac{{c_{dry}}^2}{{c_{sat}}^2}\right)\frac{{\sec }^2\theta }{4}\right]\frac{\Delta f}{f}+\left(\frac{{c_{dry}}^2}{4{c_{sat}}^2}{\sec }^2\theta -\frac{2}{{c_{sat}}^2}{\sin }^2\theta \right)\frac{\Delta \mu }{\mu }+\left(\frac{1}{2}-\frac{{\sec }^2\theta }{4}\right)\frac{\Delta \rho }{\rho },(5)$

where $\theta$ denotes the incidence angle, $\rho$ signifies the density, and ${c_{sat}}^2$ and ${c_{dry}}^2$ represent the velocity ratios of P and S wave in saturated and dry rocks, respectively.

2.2 Reflection coefficient formula for decoupling porosity and pore structure parameter

Second, in terms of petrophysics, Sun (2000) proposed a petrophysical model based on Biot theory, introducing the pore structure parameter $\gamma$ to effectively depict the effect of the pore shape on seismic velocity. Given its inclusion of both the pore structure parameter $\gamma$ and porosity $\varphi$, Sun’s model exhibits superior applicability for characterizing reservoirs with a complex pore structure compared to conventional petrophysical models:

$K_{sat}=\left(1-F_k\varphi \right)K_s+F_k\varphi K_f,(6)$

$F_k=\frac{1-{\left(1-\varphi \right)}^{\gamma }}{\left[1-{\left(1-\varphi \right)}^{\gamma }\right]\frac{K_f}{K_s}+\left(1-\frac{K_f}{K_s}\right)\varphi },(7)$

${\mu }_{sat}={\mu }_s{\left(1-\varphi \right)}^{\gamma }.(8)$

It is worth noting that the pore structure parameter $\gamma$ represents the pore shape, particularly the aspect ratio $\alpha$, which does not fall within a numerical range from 0 to 1. Through digital core simulation, Zhao et al. (2021a), Zhao et al. (2021b) clarified their non-linear relationship: $\alpha <0.08$ when $\gamma >8$ and $1>\alpha >0.3$ when $2<\gamma <3$.

According to the derivation by Han et al. (2004), it can be concluded that $K_s\gg K_f$ when $\varphi$ is extremely small ($\varphi <15\%$). Additionally, in tight sandstone or carbonate with low porosity, $\gamma$ exhibits a significant influence on the seismic response, which greatly contributes to heterogeneity. Therefore, formula (6) can be simplified as follows:

$K_{sat}={\left(1-\varphi \right)}^{\gamma }{K}_s+\left[1-{\left(1-\varphi \right)}^{\gamma }\right]K_f.(9)$

It can be found by comparing formula 9, 4that

$f=\left[1-{\left(1-\varphi \right)}^{\gamma }\right]K_f,(10)$

$K_{dry}={\left(1-\varphi \right)}^{\gamma }{K}_s.(11)$

The exponential term in formula 10 poses challenges for numerical analysis, and therefore, it can be expanded by the Taylor series as follows:

$f=\left\{1-\left[1-\gamma \varphi +\frac{1}{2}\gamma \left(\gamma -1\right){\varphi }^2-\cdots \right]\right\}{K}_f.(12)$

The selection of different orders in approximation (12) can be tailored to meet different accuracy requirements. To achieve the decoupling of $K_f$ from $\varphi$ and $\gamma$ by direct inversion, formula 12 in the first-order approximation is substituted into formula 5 to obtain a novel fluid prediction formula as follows:

$\begin{array}{l}{R}_{pp}\left(\theta \right)={\sec }^2\theta \left(\frac{1}{4}-\frac{{c_{dry}}^2}{4{c_{sat}}^2}\right)\frac{\Delta K_f}{K_f}+\left(\frac{{\sec }^2\theta }{4}-\frac{{c_{dry}}^2}{2{c_{sat}}^2}{\sec }^2\theta +\frac{2}{{c_{sat}}^2}{\sin }^2\theta \right)\frac{\Delta \varphi }{\varphi }\\ +\left[\left(\frac{1}{2}-\frac{{c_{dry}}^2}{4{c_{sat}}^2}\right){\sec }^2\theta -\frac{1}{2}\right]\frac{\Delta \gamma }{\gamma }+\left(\frac{{c_{dry}}^2}{4{c_{sat}}^2}{\sec }^2\theta -\frac{2}{{c_{sat}}^2}{\sin }^2\theta \right)\frac{\Delta f_m}{f_m}+\left(\frac{1}{2}-\frac{1}{4}{\sec }^2\theta \right)\frac{\Delta f_s}{f_s},\end{array}(13)$

where $f_m=\varphi \mu$ represents the dry rock matrix term (Yin et al., 2014) and $f_s=\gamma \rho$ denotes the structural density term.

3 Accuracy analysis

For verifying the rationality of the proposed reflection coefficient in formula 13 for decoupling the pore structure parameter $\gamma$ and porosity $\varphi$ simultaneously, a two-layer geological model based on the logging data from carbonate reservoirs is designed to analyze its accuracy. The model parameters are shown in Table 1. The exact Zoeppritz formula, Aki–Richards approximation, and the proposed formula are used to calculate the reflection coefficient, and their comparison is shown in Figure 1. The reflection coefficient obtained from the new formula exhibits good agreement with that from the exact Zoeppritz formula for incidence angle below 35°, while maintaining similar accuracy to the A&K approximation for incidence angle above 35°. Therefore, Figure 1 shows that there is no loss of computational accuracy in the proposed formula of separating the pore structure despite the increase in input parameters.

TABLE 1

TABLE 1. Parameters of the two-layer strata model.

FIGURE 1

FIGURE 1. Accuracy comparison of reflection coefficients among the Zoeppritz formula (×), Aki–Richard formula (+), and proposed formula of decoupling the pore structure (○).

4 Importance of the pore structure

It should be noted that there are five terms to be solved in the proposed approximation (13), including porosity $\varphi$ and the pore structure parameter $\gamma$. Since these two parameters are not commonly used in a seismic inversion, it is necessary to confirm their considerable influence on the AVO reflection coefficient.

A two-layer model with the same basic elastic parameters ($V_p=4000\mathrm{m}/\mathrm{s}$; $V_S=2350\mathrm{m}/\mathrm{s}$; $\rho =2450{\text{kg}/\mathrm{m}}^3$) is built to reveal the effect of the pore structure on the reflection coefficient by varying only $\varphi$ or $\gamma$, as shown in Figure 2. Figure 2A shows that when $\gamma$ is fixed at 5 and $\varphi$ changes from 1% to 10%, the reflection coefficient exhibits significant changes. Similarly, in Figure 2B, when $\varphi$ is fixed at 6% and $\gamma$ changes from 2 to 8, corresponding to the transition from micro-fractures to vugular pores in pore types, noticeable variations in the reflection coefficient are observed. The forward simulation results indicate that both $\varphi$ and $\gamma$ have a considerable impact on the seismic response, akin to conventional elastic parameters. Upon comparing Figure 2A with (b), it is observed that at low values of $\varphi$, the seismic response of $\gamma$ shows similarities with that of $\varphi$. The consistent variation trend of reflection coefficient in both figures suggests that an increase in the complexity of the pore structure leads to a more pronounced seismic response. Therefore, the influence of $\gamma$ on the variation of seismic reflection should not be neglected, necessitating the consideration of both $\varphi$ and $\gamma$ as pivotal factors in fluid prediction.

FIGURE 2

FIGURE 2. Contribution of porosity $\varphi$ and pore structure parameter $\gamma$ to seismic reflection coefficients. (A) Effect of variation in $\varphi$ on reflection coefficients when $\gamma =5$; (B) Effect of variation in $\gamma$ on reflection coefficients when $\varphi =6\%$.

5 Inversion for decoupling the pore structure

5.1 Bayesian framework

The solution of the proposed method, as indicated by formula 13, involves five parameters to be solved for. This places higher demands on the stability and accuracy of the inversion compared to the traditional AVO method, which typically require only three parameters. The proposed formula requires the association of $M$ (where $M$ takes an integer greater than or equal to five) equations, with each equation corresponding to an incidence angle. Assuming that each seismic trace has $M$ incidence angles and $N$ sampling points, let $R_{K_f}=\Delta K_f/\overline{K_f}$, $R_{\varphi }=\Delta \varphi /\overline{\varphi }$, $R_{\gamma }=\Delta \gamma /\overline{\gamma }$, $R_{f_m}=\frac{\Delta f_m}{\overline{f_m}}$, and $R_{f_s}=\Delta f_s/\overline{f_s}$; consequently, formula 13 can be rewritten as follows:

$\left[\begin{array}{cccc}\begin{array}{c}{A}_{pp}\left({\theta }_1\right)\\ A_{pp}\left({\theta }_2\right)\\ ⋮\\ A_{pp}\left({\theta }_M\right)\end{array}& \begin{array}{c}{B}_{pp}\left({\theta }_1\right)\\ B_{pp}\left({\theta }_2\right)\\ ⋮\\ B_{pp}\left({\theta }_M\right)\end{array}& \begin{array}{c}{C}_{pp}\left({\theta }_1\right)\\ C_{pp}\left({\theta }_2\right)\\ ⋮\\ C_{pp}\left({\theta }_M\right)\end{array}& \begin{array}{cc}\begin{array}{c}{D}_{pp}\left({\theta }_1\right)\\ D_{pp}\left({\theta }_2\right)\\ ⋮\\ D_{pp}\left({\theta }_M\right)\end{array}& \begin{array}{c}{E}_{pp}\left({\theta }_1\right)\\ E_{pp}\left({\theta }_2\right)\\ ⋮\\ E_{pp}\left({\theta }_M\right)\end{array}\end{array}\end{array}\right]\left[\begin{array}{c}{R}_{K_f}\\ R_{\varphi }\\ R_{\gamma }\\ \begin{array}{c}{R}_{f_m}\\ R_{f_s}\end{array}\end{array}\right]=\left[\begin{array}{c}{R}_{pp}\left({\theta }_1\right)\\ R_{pp}\left({\theta }_2\right)\\ ⋮\\ R_{pp}\left({\theta }_M\right)\end{array}\right],(14)$

where $A_{pp}\left({\theta }_i\right)={\sec }^2\theta \left(\frac{1}{4}-\frac{{c_{dry}}^2}{4{c_{sat}}^2}\right)\cdot I_{N\times N}$, $B_{pp}\left({\theta }_i\right)=\left(\frac{{\sec }^2\theta }{4}-\frac{{c_{dry}}^2}{2{c_{sat}}^2}{\sec }^2\theta +\frac{2}{{c_{sat}}^2}{\sin }^2\theta \right)\cdot I_{N\times N}$, $C_{pp}\left({\theta }_i\right)=\left(\frac{1}{2}{\sec }^2\theta -\frac{1}{2}-\frac{{c_{dry}}^2}{4{c_{sat}}^2}{\sec }^2\theta \right)\cdot I_{N\times N}$, $D_{pp}\left({\theta }_i\right)=\left(\frac{{c_{dry}}^2}{4{c_{sat}}^2}{\sec }^2\theta -\frac{2}{{c_{sat}}^2}{\sin }^2\theta \right)\cdot I_{N\times N}$, $E_{pp}\left({\theta }_i\right)=\left(\frac{1}{2}-\frac{1}{4}{\sec }^2\theta \right)\cdot I_{N\times N}$, and $I_{N\times N}$ denotes a unit matrix of $N\times N$.

For simplicity of expression, formula 14 can be reformulated as follows:

$G_{pp}m=d_{pp},(15)$

where $G_{pp}=\left[\begin{array}{cccc}\begin{array}{c}{A}_{pp}\left({\theta }_1\right)\\ A_{pp}\left({\theta }_2\right)\\ ⋮\\ A_{pp}\left({\theta }_M\right)\end{array}& \begin{array}{c}{B}_{pp}\left({\theta }_1\right)\\ B_{pp}\left({\theta }_2\right)\\ ⋮\\ B_{pp}\left({\theta }_M\right)\end{array}& \begin{array}{c}{C}_{pp}\left({\theta }_1\right)\\ C_{pp}\left({\theta }_2\right)\\ ⋮\\ C_{pp}\left({\theta }_M\right)\end{array}& \begin{array}{cc}\begin{array}{c}{D}_{pp}\left({\theta }_1\right)\\ D_{pp}\left({\theta }_2\right)\\ ⋮\\ D_{pp}\left({\theta }_M\right)\end{array}& \begin{array}{c}{E}_{pp}\left({\theta }_1\right)\\ E_{pp}\left({\theta }_2\right)\\ ⋮\\ E_{pp}\left({\theta }_M\right)\end{array}\end{array}\end{array}\right]$, $m={\left[\begin{array}{cccc}\begin{array}{cc}{R}_{K_f}& R_{\varphi }\end{array}& R_{\gamma }& R_{f_m}& R_{f_s}\end{array}\right]}^{\mathrm{T}}$, and $d_{pp}={\left[\begin{array}{cccc}{R}_{pp}\left({\theta }_1\right)& R_{pp}\left({\theta }_2\right)& \cdots & R_{pp}\left({\theta }_M\right)\end{array}\right]}^{\mathrm{T}}$.

The solution of formula 15 provides the five parameters; however, it suffers from severe ill-posedness in practical applications. Therefore, based on the Bayesian framework, the issue of well-posedness can be effectively improved by introducing the prior information of the model parameters into the regularization of inversion. Thus, the posterior probability density distribution of model parameters $P\left(m|d_{pp}\right)$ is

$P\left(m|d_{pp}\right)=\frac{P\left(m\right)P\left(d_{pp}|m\right)}{P\left(d_{pp}\right)},(16)$

where $P\left(m\right)$ signifies the prior distribution of the model parameters, $P\left(m\right)P\left(d_{pp}|m\right)$ represents the likelihood function that characterizes the noise distribution of seismic gathers, and $P\left(d_{pp}\right)$ denotes a constant that can be ignored if the posterior probability distribution function remains unchanged.

The prestack seismic gathers in actual data contain a certain level of noise. Assuming that the noise is uncorrelated with a mean value of 0, following the likelihood function of the Gaussian distribution, thus, the likelihood function of noise can be expressed as follows:

$P\left(d_{pp}|m\right)={\left[{\left(2\pi \right)}^{\mathit{LN}}\left|C_{np}\right|\right]}^{-\frac{1}{2}}\cdot \exp \left[-\frac{1}{2}{\left(d_{pp}-G_{pp}m\right)}^{\mathrm{T}}{C}_{np}^{-1}\left(d_{pp}-G_{pp}m\right)\right],(17)$

where $C_{np}={\sigma }_{np}^{2}I$ and ${\sigma }_{np}^2$ represent the covariance and the mean square error of P-wave noise, respectively.

The prior distribution of the model parameters can be either univariate or multivariate. Considering that five involved parameters, which have certain correlation, need to be obtained by the proposed method, the selection of multivariate distribution can effectively mitigate the ill-conditioning issue. While Gaussian distribution only provides consistent weighting coefficients, which may impact the sparsity of inversion results, Cauchy distribution can yield non-consistent weighting coefficients with sparsity effects and greater geological significance. Therefore, the multivariate Cauchy distribution is adopted as the prior distribution for the model parameters, with its specific formula given as follows:

$P\left(m\right)={\displaystyle \prod _{i-1}^N}\frac{2{\left|\psi \right|}^{-\frac{1}{2}}}{{\pi }^2{\left(1+m^{\mathrm{T}}{\mathrm{\Phi }}_{i}m\right)}^2},{\mathrm{\Phi }}_i=D_i^{\mathrm{T}}{\psi }^{-1}{D}_i,(18)$

where $\psi$ signifies the correlation matrix, which can be obtained by maximum expectation estimation. $D_i$ denotes the matrix of $5N\times 5N$, and its expression is as follows:

${\left[D_i\right]}_{xy}=\left\{\begin{array}{cc}\begin{array}{c}1\\ 1\\ 1\\ \begin{array}{c}1\\ 1\\ 0\end{array}\end{array}& \begin{array}{c}ifx=1andy=i\\ ifx=1andy=i+N\\ ifx=1andy=i+2N\\ \begin{array}{c}ifx=1andy=i+3N\\ ifx=1andy=i+4N\\ otherwise\end{array}\end{array}\end{array}\right..(19)$

By substituting formulas 17, 18 into Formula 16, the posterior probability density of model parameters can be obtained as follows:

$P\left(m|d_{pp}\right)\propto {\displaystyle \prod _{i-1}^N}\frac{2{\left|\psi \right|}^{-\frac{1}{2}}}{{\pi }^2{\left(1+m^{\mathrm{T}}{\mathrm{\Phi }}_{i}m\right)}^2}\cdot \exp \left[-\frac{1}{2}{\left(d_{pp}-G_{pp}m\right)}^{\mathrm{T}}{C}_{np}^{-1}\left(d_{pp}-G_{pp}m\right)\right].(20)$

After the algebraic transformation of Formula 20, the objective function $F\left(m\right)$ should be as follows:

$F\left(m\right)={\left(d_{pp}-G_{pp}m\right)}^{\mathrm{T}}\left(d_{pp}-G_{pp}m\right)+\tau {\displaystyle \sum _{i=1}^N}\ln \left(1+m^{\mathrm{T}}{\mathrm{\Phi }}_{i}m\right),(21)$

where $\tau$ determines the sparsity of inversion results, with larger values leading to sparser outcomes. Conversely, if $\tau$ is set too small, it may result in a distorted inversion.

5.2 Decorrelation method

The statistical correlation between elastic parameters is widely recognized as a major factor contributing to the instability of the prestack inversion. The proposed method involves five parameters, out of which three parameters (fluid bulk modulus $K_f$, porosity $\varphi$, and pore structure parameter $\gamma$) are further decomposed by the fluid factor $f$, resulting in a closer relationship among them, thus leading to a substantial increase in inversion instability. Evidently, relying solely on the prior probability distribution scheme is insufficient for eliminating the parameter correlations. To address this issue, an improved decorrelation method (Wang et al., 2017) is adopted to obtain parameters with the lowest possible correlation, which uses the variance matrix and two linear transformations based on Chen’s previous work (Chen et al., 2007) to convert relevant data into “white data.”

The decorrelation method is illustrated using a three-parameter dataset as an example, since it involves five parameters in the proposed inversion, which is not convenient for figuring. As shown in Figure 3A, the sequences x, y, and z exhibit significant correlation, and their covariance matrix is expressed as follows:

$\mathrm{\Sigma }=\left[\begin{array}{ccc}{\sigma }_x^2& {\sigma }_{xy}& {\sigma }_{xz}\\ {\sigma }_{xy}& {\sigma }_y^2& {\sigma }_{yz}\\ {\sigma }_{xz}& {\sigma }_{yz}& {\sigma }_z^2\end{array}\right],(22)$

where ${\sigma }_x^2$, ${\sigma }_y^2$, and ${\sigma }_z^2$ represent the variance of $x$, $y$, and $z$, respectively, and ${\sigma }_{xy}$, ${\sigma }_{xz}$, and ${\sigma }_{yz}$ denote the covariance of the three parameters. The singular value decomposition of the covariance matrix is

$\mathrm{\Sigma }=vL{v}^T,(23)$

where $v$ signifies the eigenvector matrix and $L$ denotes the eigenvalue matrix. Upon applying $v^{-1}$ to $d={\left[x,y,z\right]}^T$, the resulting sequences $x^{\prime },y^{\prime },z^{\prime }$ are shown in Figure 3B. However, despite this transformation, the newly obtained data $d^{\prime }={\left[x^{\prime },y^{\prime },z^{\prime }\right]}^T$ still have a weak correlation. To achieve further decorrelation, we apply $s^{-1}={\sqrt{L}}^{-1}$ to $d^{\prime }$, resulting in $d^{″}=s^{-1}{d}^{\prime }$. Figure 3C shows that $d^{″}={\left[x^{″},y^{″},z^{″}\right]}^T$ can be referred to as “white data,” where $x^{″},y^{″},z^{″}$ are almost completely unrelated.

FIGURE 3

FIGURE 3. Decorrelation example: (A) original data with correlation, (B) data with weak correlation, and (C) white data.

According to the aforementioned method, one of the covariance matrices $Cr$ of $R_{K_f}$, $R_{\varphi }$, $R_{\gamma }$, $R_{f_m}$, and $R_{f_s}$ can be decomposed as

$Cr=vuu{v}^T.(24)$

The covariance matrix can be extended to two $5N\times 5N$ sparse eigenvector matrix $V$ and eigenvalues matrix $U$ by considering $N$ time samples. Therefore, the transformation observed in Eq. 21 can be represented as follows:

$F\left(m\right)={\left(d_{pp}-G_{pp}^{\prime }{m}^{\prime }\right)}^{\mathrm{T}}\left(d_{pp}-G_{pp}^{\prime }{m}^{\prime }\right)+\tau {\displaystyle \sum _{i=1}^N}\ln \left(1+m{\prime }^{\mathrm{T}}{\mathrm{\Phi }}_i{m}^{\prime }\right),(25)$

where $G^{\prime }=GVU$ and $m^{\prime }=U^{-1}{V}^{-1}m$. The reflection coefficients of $K_f$, $\varphi$, $\gamma$, $f_m$, and $f_s$ can be obtained by using the iterative reweighed least-squares to solve the objective function (25), followed by obtaining the five parameters through the trace integral.

6 Synthetic data test

The proposed inversion of decoupling porosity $\varphi$ and pore structure parameter $\gamma$ involves a total of five parameters, and thus, it is imperative to verify its feasibility and anti-noise capabilities. A set of well data from carbonates in the actual working area of the Sichuan Basin was selected for testing purposes. First, prior to inversion, the measured data underwent the Backus averaging process (Backus et al., 1968) which transformed the data from log scale to seismic scale, followed by time–depth conversion to convert it from depth domain to time domain. Second, the reflection coefficients of the well data at various sampling times and angles (5°–60°) are obtained by forward simulation based on the exact Zoeppritz equation, subsequently convolved with a 30 Hz Ricker wavelet to generate synthetic seismic data, as shown in Figure 4A. Third, Gaussian random noise with the signal-to-noise ratio (S/N) of 2 was added to the synthetic record, as shown in Figure 4B.

FIGURE 4

FIGURE 4. Synthetic seismic gathers without noise (A) and with noise of S/N=2 (B).

Figures 5A, B show the inversion results without and with noise, respectively, displaying the three direct inversion parameters of fluid bulk modulus $K_f$, which characterize reservoirs, porosity $\varphi$, and pore structure parameter $\gamma$. Additionally, the Russell fluid factor $f$ is calculated indirectly using formula 10 or (12) following the proposed inversion. In the absence of noise, both the direct and indirect inversions yield consistent results with well curves. However, due to a weaker seismic response than the fluid bulk modulus $K_f$, porosity $\varphi$, and pore structure parameter $\gamma$ exhibit slightly lower resolution. At S/N = 2, the direct inversion results ($K_f$, $\varphi$, and $\gamma$) show a slight decrease in accuracy compared to noise-free conditions, while errors became more apparent in indirect inversion’s $f$ due to error accumulation from indirect calculation. Nevertheless, even with added noise during the inversion process, the inversion results still maintain a similar trend as observed in the well data, indicating a good stability of the proposed method for practical applications.

FIGURE 5

FIGURE 5. Inversion results from the synthetic data without noise (A) and with noise of S/N=2 (B) (the green line indicates the low-frequency data, the blue line indicates the model data, and the red line indicates the inversion data).

7 Actual data application

The well data from an exploration area in the Sichuan Basin is utilized to validate the feasibility and applicability of the proposed prestack inversion technique for decoupling both pore structure parameter $\gamma$ and porosity $\varphi$ simultaneously. The target reservoir comprises dolomitic gas-bearing carbonate with low porosity but an intricate pore structure, which constitutes one of the influential factors affecting fluid prediction.

7.1 Digital core analysis

Before fluid detection, component analysis and simulation of different pore shapes were conducted on a total of seven cores from three wells in various formations to construct multi-component digital core models for evaluating the influence of the pore structure in this working area. Figure 6A shows a gray-scale slice obtained through a CT scan depicting the actual core. Based on the analysis of the actual core, the multi-component digital core was constructed, as shown in Figure 6B. The digital core comprised two dominant minerals: dolomite (light-colored) and quartz (dark-colored), and the dark-colored component also containing 7% kaolinite. As numerous previous studies have demonstrated the effect of porosity $\varphi$ on elastic parameters, such as those by Han et al. (2004) and Yin et al. (2014), the effect of the pore structure parameter $\gamma$ on elastic parameters within different components in the digital cores is solely focused on being discussed. The relationship between pore structure parameter $\gamma$ and saturated bulk modulus $K$ or shear modulus $\mu$ of the digital core is shown in Figures 7A, B. The $K$ of each mineral component is greatly affected by $\gamma$ when $\gamma >4$, with the $K$ of dolomite experiencing a change of over 40% and quartz approximately 32%, in particular. However, the effect of $\gamma$ on $\mu$ is negligible and can be disregarded, as all three mineral components vary by less than 10% when $\gamma >7$. However, in the simulation of filling oil or gas, the measured fluid bulk modulus $K_f$ remains constant regardless of variations in the three mineral components or changes in the pore structure parameter $\gamma$. Therefore, if fluid identification does not take the impact of the pore structure parameter $\gamma$ into account, its prediction results will exhibit significant multi-solution possibilities and uncertainties.

FIGURE 6

FIGURE 6. Core analysis. (A) Gray-scale slice of core using a CT scan; (B) digital core image (he light-colored component is dolomite, and the dark-colored component is 93% quartz and 7% kaolinite).

FIGURE 7

FIGURE 7. Impact of the pore structure parameter $\gamma$ on bulk modulus $K$ (A) and shear modulus $\mu$ (B) of different mineral components.

7.2 Fluid prediction

Reservoir analysis is initially performed on the well data, using the same well in the anti-noise experiment, as shown in Figure 8. Among the three logged carbonate gas-bearing reservoirs (highlighted by red and blue boxes in Figure 8), Russell fluid factor $f$ exhibited conspicuous manifestation solely within the second reservoir, and fluid bulk modulus $K_f$ demonstrated consistent characteristics in all the three reservoirs. According to Eqs 3, 10, Russell fluid factor $f$ serves as an integrated indicator of fluid within the saturated porous rock, whereas fluid bulk modulus $K_f$ acts as a fluid indicator independent of porosity $\varphi$ and pore structure parameter $\gamma$. Therefore, by analyzing the curves of $\varphi$ and $\gamma$, it becomes discernible as to which factor affects the prediction outcomes: The red-boxed reservoir is primarily influenced by $\varphi$, whereas the blue-boxed one is affected by $\gamma$. Notably, the third reservoir represents a thinly interbedded formation controlled by both $\varphi$ and $\gamma$, thereby indicating its highly heterogeneous nature.

FIGURE 8

FIGURE 8. Factor analysis of affecting reservoir identification (reservoirs in the red boxes are controlled by porosity $\varphi$, and those in blues ones are controlled by the pore structure parameter $\gamma$).

Second, both the proposed inversion and Russell’s $f$ inversion are calculated for fluid prediction. Additionally, $f$ indirectly obtained from the results of the proposed inversion by formula 10 is also obtained for comparison, as shown in Figure 9. The outcomes of fluid bulk modulus $K_f$ (Figure 9A), porosity $\varphi$ (Figure 9B), and pore structure parameter $\gamma$ (Figure 9C) closely match the data from the test well. Moreover, the profile of direct $f$ (Figure 9E) coincides with that of indirect $f$ (Figure 9D), indicating the stability and reliability of the proposed inversion. The result of $K_f$ (Figure 9A), which is decoupled from $\varphi$ and $\gamma$, exhibits a more pronounced distribution pattern in reservoirs compared to that of $f$. Notably, between them, there exists a significant disparity in fluid prediction for the third reservoir—the thin interbedded gas-bearing formation. By magnifying the red box near the well in Figures 9A, E,as shown in Figure 10, it can be observed that $K_f$ (Figure 10A) displays better consistency with verified gas production tests within the thin interbedded gas reservoirs than $f$ does (Figure 10B), which shows that after eliminating the influence of $\varphi$ and $\gamma$ through inversion, the prediction of $K_f$ accurately reflects reservoir distribution.

FIGURE 9

FIGURE 9. Inversion profile of fluid bulk modulus $K_f$ (A), porosity $\varphi$ (B), pore structure parameter $\gamma$ (C), and Russell fluid factor $f$ indirectly (D) and directly (E).

FIGURE 10

FIGURE 10. Magnified partial profiles of fluid bulk modulus $K_f$ (A) from Figure 9A and Russell fluid factor $f$ (B) from Figure 9E.

8 Conclusion

The pore structure parameter $\gamma$, like porosity $\varphi$, is one of the dominant factors affecting fluid prediction in heterogeneous reservoirs. By combining Russell’s poroelasticity theory and Sun’s petrophysical model, a new reflection coefficient formula is proposed, which decouples fluid bulk modulus $K_f$ from porosity $\varphi$ and pore structure parameter $\gamma$, thus eliminating the interference of the pore-related factors on reservoir prediction. Theoretical experiments demonstrate that the proposed inversion retains comparable computational accuracy to the conventional A&K approximation and has good anti-noise ability under the condition of five involved parameters. Both theoretical model and digital core analyses reveal that the pore structure parameter $\gamma$ can have a significant impact on the prediction of heterogeneous reservoirs with low porosity. In the example from the Sichuan Basin, compared with Russell fluid factor $f$, fluid bulk modulus $K_f$ from the proposed inversion method provides more accurate and distinct reservoir distribution in carbonate gas-bearing strata. However, there are still two prerequisites for the application of the proposed inversion. First, as a less commonly used parameter, the pore structure parameter $\gamma$ needs to be obtained from sufficient logging data in the working area under generally industrial process. The method used to obtain the pore structure parameter $\gamma$ in this paper is the empirical method proposed by Zhang et al. (2018a), Zhang et al. (2018b), which is still inadequate in terms of efficiency and verifiability. Second, in order to ensure the robustness of the proposed inversion method with five parameters, high-quality seismic data are also indispensable.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material; further inquiries can be directed to the corresponding author.

Author contributions

RH: methodology and writing-original draft. EW: conceptualization and project administration. WW: methodology and visualization. GY: data curation and formal analysis. XZ: validation and writing-review and editing. WZ: investigation and formal analysis. DH: conceptualization and investigation.

Funding

This work was funded by the Key Technology Research Project of CNPC, grant number 2021ZG03, and the Science and Technology Project of Petrochina, grant number 2022DJ8004.

Conflict of interest

Authors RH, EW, WW, GY, XZ, WZ, and DH were employed by PetroChina.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Backus, G., and Gilbert, F. (1968). The resolving power of gross Earth data. Geophys. J. Int. 16 (2), 169–205. doi:10.1111/j.1365-246x.1968.tb00216.x

CrossRef Full Text | Google Scholar

Berryman, J. G. (1999). Origin of Gassmann’s equations. GEOPHYSICS 64 (5), 1627–1629. doi:10.1190/1.1444667

CrossRef Full Text | Google Scholar

Biot, M. A. (1941). General theory of three-dimensional consolidation. J. Appl. Phys. 12 (2), 155–164. doi:10.1063/1.1712886

CrossRef Full Text | Google Scholar

Chen, J. J., and Yin, X. Y. (2007). Three-parameter AVO waveform inversion based on Bayesian theorem. Chin. J. Geophys. 50 (4), 1251–1260. doi:10.1007/s11442-007-0020-2

CrossRef Full Text | Google Scholar

Deng, J., Zhou, H., Wang, H., Zhao, J., and Wang, S. (2015). The influence of pore structure in reservoir sandstone on dispersion properties of elastic waves. Chin. J. Geophys. 58 (9), 3389–3400. doi:10.6038/cjg20150931

CrossRef Full Text | Google Scholar

Du, B., Yang, W., Zhang, J., Yong, X., Gao, J., and Li, H. (2019). Matrix-fluid decoupling-based joint PP-PS wave seismic inversion for fluid identification. GEOPHYSICS 84 (3), R477–R487. doi:10.1190/geo2017-0376.1

CrossRef Full Text | Google Scholar

Fan, X., Zhang, G., and Zhang, J. (2019). Prediction method of pore structure parameters of tight sandstone. Seg. Tech. Program Expand. Abstr. 89, 3678–3682. doi:10.1190/segam2019-3199837.1

CrossRef Full Text | Google Scholar

Gassmann, F. (1951). Über die elastizitat poroser medien: vierteljahrschrift der naturforschenden. Gesellshaft 96, 1–23.

Google Scholar

Han, D., and Batzle, M. L. (2004). Gassmann’s equation and fluid-saturation effects on seismic velocities. GEOPHYSICS 69 (2), 398–405. doi:10.1190/1.1707059

CrossRef Full Text | Google Scholar

He, X., He, Z., Wang, X., Xiong, X., and Jiang, L. (2012). Rock skeleton models and seismic porosity inversion. Appl. Geophys. 9 (3), 349–358. doi:10.1007/s11770-012-0345-1

CrossRef Full Text | Google Scholar

He, X., Lin, K., Zhang, Z., and He, Z. (2018). Quantitative comparison of the influence of porosity and pore structure on reservoir characteristics. Geophys. Prospect. Petroleum 57 (2), 179–185.

Google Scholar

Jiang, L., Wen, X., Zhou, D., He, Z., and He, X. (2012). The constructing of pore structure factor in carbonate rocks and the inversion of reservoir parameters. Appl. Geophys. 9 (2), 223–232. doi:10.1007/s11770-012-0333-5

CrossRef Full Text | Google Scholar

Kachanov, M., Tsukrov, I., and Shafiro, B. (1994). Effective moduli of solids with cavities of various shapes. Appl. Mech. Rev. 47 (1S), S151–S174. doi:10.1115/1.3122810

CrossRef Full Text | Google Scholar

Krief, M., Garat, J., Stellingwerff, J., and Ventre, J. (1990). A petrophysical interpretation using the velocities of P and S waves. Log. Anal. 31, 355–369.

Google Scholar

Li, H., Zhang, J., Pan, H., and Gao, Q. (2021). Nonlinear simultaneous inversion of pore structure and physical parameters based on elastic impedance. Sci. China Earth Sci. 64 (6), 977–991. doi:10.1007/s11430-020-9773-8

CrossRef Full Text | Google Scholar

Russell, B. H., Gray, D., and Hampson, D. P. (2011). Linearized AVO and poroelasticity. GEOPHYSICS 76 (3), C19–C29. doi:10.1190/1.3555082

CrossRef Full Text | Google Scholar

Russell, B. H., Hedlin, K., Hilterman, F. J., and Lines, L. R. (2003). Fluid-property discrimination with AVO: A biot-gassmann perspective. GEOPHYSICS 68 (1), 29–39. doi:10.1190/1.1543192

CrossRef Full Text | Google Scholar

Sun, R., Yin, X., Wang, B., and Zhang, G. (2016). A direct estimation method for the Russell fluid factor based on stochastic seismic inversion. Chin. J. Geophys. 59 (3), 1143–1150. doi:10.6038/cjg20160334

CrossRef Full Text | Google Scholar

Sun, W., He, Z., Li, Y., Zhang, F., and Zhou, Y. (2015). A new factor of fluid-filled fractures and its application. Chin. J. Geophys. 58 (7), 2536–2545. doi:10.6038/cjg20150728

CrossRef Full Text | Google Scholar

Sun, Y. (2000). Core-log-seismic integration in hemipelagic marine sediments on the eastern flank of the Juan de Fuca Ridge. Proc. Ocean Drill. Program, Sci. Results 168, 21–23. doi:10.2973/odp.proc.sr.168.009.2000

CrossRef Full Text | Google Scholar

Sun, Y., and Goldberg, D. (1997). Effects of aspect ratio on wave velocities in fractured rocks. Seg. Tech. Program Expand. Abstr. 67, 925–928. doi:10.1190/1.1886170

CrossRef Full Text | Google Scholar

Wang, F., Bian, H., Zhang, Y., Duan, C., and Chen, G. (2016). Hilbert-Huang transform combined with smoothed pseudo Wigner-Ville time-frequency distribution to identify reservoir fluid properties. Geophys. Prospect. Petroleum 55 (6), 851–860. doi:10.3969/j.issn.1000-1441.2016.06.010

CrossRef Full Text | Google Scholar

Wang, W., Li, G., Yu, X., Yang, W., Wang, W., Feng, L., et al. (2017). A meta-analysis of adenosine A2A receptor antagonists on levodopa-induced dyskinesia in vivo. Seg. Tech. Program Expand. Abstr. 87, 702–706. doi:10.3389/fneur.2017.00702

PubMed Abstract | CrossRef Full Text | Google Scholar

Yin, X., and Zhang, S. (2014). Bayesian inversion for effective pore-fluid bulk modulus based on fluid-matrix decoupled amplitude variation with offset approximation. GEOPHYSICS 79 (5), R221–R232. doi:10.1190/geo2013-0372.1

CrossRef Full Text | Google Scholar

Zhang, T., Sun, Y., Tian, J., Lu, L., Qin, F., Zhao, X., et al. (2018b). Two-parameter prestack seismic inversion of porosity and pore-structure parameter of fractured carbonate reservoirs: part 2 - applications. Interpretation 6 (4), SM9–SM17. doi:10.1190/int-2018-0019.1

CrossRef Full Text | Google Scholar

Zhang, T., and Sun, Y. (2018a). Two-parameter prestack seismic inversion of porosity and pore structure parameter of fractured carbonate reservoirs: part 1 — methods. Interpretation 6, 1–SM8. doi:10.1190/int-2017-0219.1

CrossRef Full Text | Google Scholar

Zhao, J., Pan, J., Hu, Y., Li, J., Liu, X., Li, C., et al. (2021a). Digital rock physics-based studies on effect of pore types on elastic properties of carbonate reservoir Part 1: imaging processing and elastic modelling. Chin. J. Geophys. 64 (2), 656–669. doi:10.6038/cjg2021O0228

CrossRef Full Text | Google Scholar

Zhao, J., Pan, J., Hu, Y., Li, J., Liu, X., Li, C., et al. (2021b). Digital rock physics-based studies on effect of pore types on elastic properties of carbonate reservoir Part 2: pore structure factor characterization and inversion of reservoir. Chin. J. Geophys. 64 (2), 670–671. doi:10.23736/S0026-4806.20.06558-1

CrossRef Full Text | Google Scholar

Zimmerman, R. W. (1986). Compressibility of two-dimensional cavities of various shapes. J. Appl. Mech. 53 (3), 500–504. doi:10.1115/1.3171802

CrossRef Full Text | Google Scholar

Zong, Z., and Wang, Y. (2019). Amplitude variation with angle inversion for fluid discrimination with the consideration of squirt flow. Seg. Tech. Program Expand. Abstr. 89, 734–738. doi:10.1190/segam2019-3216431.1

CrossRef Full Text | Google Scholar

Zong, Z., Yin, X., and Wu, G. (2012). AVO inversion and poroelasticity with P- and S-wave moduli. GEOPHYSICS 77 (6), N17–N24. doi:10.1190/geo2011-0214.1

CrossRef Full Text | Google Scholar

Zong, Z., Yin, X., Wu, G., and Wu, Z. (2015). Elastic inverse scattering for fluid variation with time-lapse seismic data. GEOPHYSICS 80 (2), WA61–WA67. doi:10.1190/geo2014-0183.1

CrossRef Full Text | Google Scholar