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 s a t = K d r y + β 2 M , (1) μ s a t = μ d r y = μ , (2) where K s a t and K d r y represent the bulk moduli for the saturated and dry porous rocks, μ s a t and μ d r y denote the shear moduli for the saturated and dry porous rocks, β 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 = β 2 M = 1 − K d r y / K s 2 ϕ / K f + 1 − ϕ / K s − K d r y / K s 2 , (3) where K f and K s denote the bulk moduli of the fluid and saturated rock frame and ϕ represents the porosity. Therefore, formula (2) can be reformulated as follows: K s a t = K d r y + 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 p p θ = 1 − c d r y 2 c s a t 2 sec 2 θ 4 Δ f f + c d r y 2 4 c s a t 2 sec 2 θ − 2 c s a t 2 sin 2 θ Δ μ μ + 1 2 − sec 2 θ 4 Δ ρ ρ , (5) where θ denotes the incidence angle, ρ signifies the density, and c s a t 2 and c d r y 2 represent the velocity ratios of P and S wave in saturated and dry rocks, respectively.

Second, in terms of petrophysics, Sun (2000) proposed a petrophysical model based on Biot theory, introducing the pore structure parameter γ to effectively depict the effect of the pore shape on seismic velocity. Given its inclusion of both the pore structure parameter γ and porosity ϕ , Sun’s model exhibits superior applicability for characterizing reservoirs with a complex pore structure compared to conventional petrophysical models: K s a t = 1 − F k ϕ K s + F k ϕ K f , (6) F k = 1 − 1 − ϕ γ 1 − 1 − ϕ γ K f K s + 1 − K f K s ϕ , (7) μ s a t = μ s 1 − ϕ γ . (8)

It is worth noting that the pore structure parameter γ represents the pore shape, particularly the aspect ratio α , 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: α < 0.08 when γ > 8 and 1 > α > 0.3 when 2 < γ < 3 .

According to the derivation by Han et al. (2004), it can be concluded that K s ≫ K f when ϕ is extremely small ( ϕ < 15 % ). Additionally, in tight sandstone or carbonate with low porosity, γ exhibits a significant influence on the seismic response, which greatly contributes to heterogeneity. Therefore, formula (6) can be simplified as follows: K s a t = 1 − ϕ γ K s + 1 − 1 − ϕ γ K f . (9)

It can be found by comparing formula 9, 4that f = 1 − 1 − ϕ γ K f , (10) K d r y = 1 − ϕ γ 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 = 1 − 1 − γ ϕ + 1 2 γ γ − 1 ϕ 2 − ⋯ 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 ϕ and γ by direct inversion, formula 12 in the first-order approximation is substituted into formula 5 to obtain a novel fluid prediction formula as follows: R p p θ = sec 2 θ 1 4 − c d r y 2 4 c s a t 2 Δ K f K f + sec 2 θ 4 − c d r y 2 2 c s a t 2 sec 2 θ + 2 c s a t 2 sin 2 θ Δ ϕ ϕ + 1 2 − c d r y 2 4 c s a t 2 sec 2 θ − 1 2 Δ γ γ + c d r y 2 4 c s a t 2 sec 2 θ − 2 c s a t 2 sin 2 θ Δ f m f m + 1 2 − 1 4 sec 2 θ Δ f s f s , (13) where f m = ϕ μ represents the dry rock matrix term (Yin et al., 2014) and f s = γ ρ denotes the structural density term.

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 = Δ K f / K f ¯ , R ϕ = Δ ϕ / ϕ ¯ , R γ = Δ γ / γ ¯ , R f m = Δ f m f m ¯ , and R f s = Δ f s / f s ¯ ; consequently, formula 13 can be rewritten as follows: A p p θ 1 A p p θ 2 ⋮ A p p θ M B p p θ 1 B p p θ 2 ⋮ B p p θ M C p p θ 1 C p p θ 2 ⋮ C p p θ M D p p θ 1 D p p θ 2 ⋮ D p p θ M E p p θ 1 E p p θ 2 ⋮ E p p θ M R K f R ϕ R γ R f m R f s = R p p θ 1 R p p θ 2 ⋮ R p p θ M , (14) where A p p θ i = sec 2 θ 1 4 − c d r y 2 4 c s a t 2 ⋅ I N × N , B p p θ i = sec 2 θ 4 − c d r y 2 2 c s a t 2 sec 2 θ + 2 c s a t 2 sin 2 θ ⋅ I N × N , C p p θ i = 1 2 sec 2 θ − 1 2 − c d r y 2 4 c s a t 2 sec 2 θ ⋅ I N × N , D p p θ i = c d r y 2 4 c s a t 2 sec 2 θ − 2 c s a t 2 sin 2 θ ⋅ I N × N , E p p θ i = 1 2 − 1 4 sec 2 θ ⋅ I N × N , and I N × N denotes a unit matrix of N × N .

For simplicity of expression, formula 14 can be reformulated as follows: G p p m = d p p , (15) where G p p = A p p θ 1 A p p θ 2 ⋮ A p p θ M B p p θ 1 B p p θ 2 ⋮ B p p θ M C p p θ 1 C p p θ 2 ⋮ C p p θ M D p p θ 1 D p p θ 2 ⋮ D p p θ M E p p θ 1 E p p θ 2 ⋮ E p p θ M , m = R K f R ϕ R γ R f m R f s T , and d p p = R p p θ 1 R p p θ 2 ⋯ R p p θ M 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 m | d p p is P m | d p p = P m P d p p | m P d p p , (16) where P m signifies the prior distribution of the model parameters, P m P d p p | m represents the likelihood function that characterizes the noise distribution of seismic gathers, and P d p p 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 d p p | m = 2 π LN C n p − 1 2 ⋅ exp − 1 2 d p p − G p p m T C n p − 1 d p p − G p p m , (17) where C n p = σ n p 2 I and σ n p 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 m = ∏ i − 1 N 2 ψ − 1 2 π 2 1 + m T Φ i m 2 , Φ i = D i T ψ − 1 D i , (18) where ψ signifies the correlation matrix, which can be obtained by maximum expectation estimation. D i denotes the matrix of 5 N × 5 N , and its expression is as follows: D i x y = 1 1 1 1 1 0 i f   x = 1   a n d   y = i i f   x = 1   a n d   y = i + N i f   x = 1   a n d   y = i + 2 N i f   x = 1   a n d   y = i + 3 N i f   x = 1   a n d   y = i + 4 N o t h e r w i s e . (19)

By substituting formulas 17, 18 into Formula 16, the posterior probability density of model parameters can be obtained as follows: P m | d p p ∝ ∏ i − 1 N 2 ψ − 1 2 π 2 1 + m T Φ i m 2 ⋅ exp − 1 2 d p p − G p p m T C n p − 1 d p p − G p p m . (20)

After the algebraic transformation of Formula 20, the objective function F m should be as follows: F m = d p p − G p p m T d p p − G p p m + τ ∑ i = 1 N ln 1 + m T Φ i m , (21) where τ determines the sparsity of inversion results, with larger values leading to sparser outcomes. Conversely, if τ is set too small, it may result in a distorted inversion.

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 ϕ , and pore structure parameter γ ) 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: Σ = σ x 2 σ x y σ x z σ x y σ y 2 σ y z σ x z σ y z σ z 2 , (22) where σ x 2 , σ y 2 , and σ z 2 represent the variance of x , y , and z , respectively, and σ x y , σ x z , and σ y z denote the covariance of the three parameters. The singular value decomposition of the covariance matrix is Σ = v L v T , (23) where v signifies the eigenvector matrix and L denotes the eigenvalue matrix. Upon applying v − 1 to d = x , y , z T , the resulting sequences x ′ , y ′ , z ′ are shown in Figure 3B. However, despite this transformation, the newly obtained data d ′ = x ′ , y ′ , z ′ T still have a weak correlation. To achieve further decorrelation, we apply s − 1 = L − 1 to d ′ , resulting in d ″ = s − 1 d ′ . Figure 3C shows that d ″ = x ″ , y ″ , z ″ T can be referred to as “white data,” where x ″ , y ″ , z ″ are almost completely unrelated.

According to the aforementioned method, one of the covariance matrices C r of R K f , R ϕ , R γ , R f m , and R f s can be decomposed as C r = v u u v T . (24)

The covariance matrix can be extended to two 5 N × 5 N 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 m = d p p − G p p ′ m ′ T d p p − G p p ′ m ′ + τ ∑ i = 1 N ln 1 + m ′ T Φ i m ′ , (25) where G ′ = G V U and m ′ = U − 1 V − 1 m . The reflection coefficients of K f , ϕ , γ , 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.

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 ϕ on elastic parameters, such as those by Han et al. (2004) and Yin et al. (2014), the effect of the pore structure parameter γ on elastic parameters within different components in the digital cores is solely focused on being discussed. The relationship between pore structure parameter γ and saturated bulk modulus K or shear modulus μ of the digital core is shown in Figures 7A, B. The K of each mineral component is greatly affected by γ when γ > 4 , with the K of dolomite experiencing a change of over 40% and quartz approximately 32%, in particular. However, the effect of γ on μ is negligible and can be disregarded, as all three mineral components vary by less than 10% when γ > 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 γ . Therefore, if fluid identification does not take the impact of the pore structure parameter γ into account, its prediction results will exhibit significant multi-solution possibilities and uncertainties.

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 ϕ and pore structure parameter γ . Therefore, by analyzing the curves of ϕ and γ , it becomes discernible as to which factor affects the prediction outcomes: The red-boxed reservoir is primarily influenced by ϕ , whereas the blue-boxed one is affected by γ . Notably, the third reservoir represents a thinly interbedded formation controlled by both ϕ and γ , thereby indicating its highly heterogeneous nature.

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 ϕ (Figure 9B), and pore structure parameter γ (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 ϕ and γ , 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 ϕ and γ through inversion, the prediction of K f accurately reflects reservoir distribution.