Effective porosity inversion—a key technology for lithology prediction and fluid identification—plays a significant role in oil and gas exploration. Based on the analytical expression of the vertical reflection coefficient of fast P-wave at the interface of fluid-bearing porous media, a seismic record is described as a function of the porosity by combining the quantitative relationship among rock skeleton modulus, rock matrix modulus, and effective porosity. Considering the nonlinear relationship among them, the simulated annealing (SA) method is used to solve the nonlinear inverse problem, and effective porosity inversion is realized by utilizing the reflection coefficients, porosities, and interpretation results of well logs as prior constraints. A series of numerical analyses shows that reasonable constraints can make the inversion converge rapidly to the optimal solution, and the inversion results of porosity are effective and stable and have high resolution and strong noise immunity. The calculations of well data and sandstone reservoir data further verify the practicability of this method, and the solutions are in good agreement with well log porosity.

Porosity is an important factor that controls the storage and transportation of pore fluids in rock and sedimentary bodies, which is defined as the ratio of pore volume to rock’s total volume. In reservoir prediction, the real storage capacity of rock pore to hydrocarbon is defined as effective porosity. Effective porosity prediction plays a significant role in the identification of lithology and pore fluids and reserves estimation of oil and gas reservoirs. Rock core test and well log data can be used to estimate effective porosity, but these methods are expensive and cannot obtain a two- or three-dimensional porosity data volume.

Generally, the porosity prediction of oil and gas reservoirs can be realized by combining core test data, drilling and logging data, and various seismic attributes. Aminu and Olorunniwo [1] used multivariate linear transforms to realize multiattribute predictions of lithology and porosity from seismic data. Khoshdel and Riahi [2] analyzed acoustic impedance by multiple attribute regression and neural networks to predict porosity. Jalalalhosseini et al. [3] first extracted attribute information from seismic data, then uses well log information to determine the relationship between a series of attribute information and porosity, and finally used this relationship to realize porosity prediction. For further improving the accuracy of porosity prediction, a few scholars have combined multiattribute analysis with geostatistics. Doyen [4] proposed a geostatistical method called “cokriging” to predict the porosity distribution of channel-sand reservoir, using the spatial autocorrelation function to describe the lateral change of porosity, and the correlation function to construct the dependence between porosity and acoustic properties. As P-wave attenuation and acoustic impedance are highly sensitive to porosity and saturation, Pang et al. [5] used the two attribute crossplots to establish rock physics templates and further realized the prediction of porosity, saturation, and permeability. The deep learning (DL) method has a good prediction effect on porosity because it can better capture the nonlinear relationship between seismic data and well data. Leiphart and Hart [6] used the standard linear regression and probabilistic neural network to predict porosity distribution; the latter—according to them—can obtain better porosity distribution. Pramanik et al. [7] and Maurya et al. [8] realized the direct prediction of effective porosity by integrating data including seismic data, P-wave impedance, and well logs, respectively. They used multiattribute transformation, linear regression, neural network (NN), etc., to improve the nonlinear methods and obtained better effective porosity distribution.

Compared with the multiattribute qualitative analysis of porosity, the combination of petrophysical modeling and seismic inversion can realize the semiquantitative or quantitative prediction of rock physical parameters (porosity, clay content, water saturation, etc.). Bachrach [9] realized the prediction of porosity and saturation based on the Bayesian theory, stochastic model, and rock physics, believing that these two parameters needed to be synchronously predicted and that the uncertainty of porosity will lead to an error in saturation prediction, and vice versa. Calderon and Castagna [10] integrated rock physics, seismic inversion, and multiattribute transforms to estimate reservoir rock properties of porosity and lithology. Leite and Vidal [11] proposed using sparse-spike inversion to convert seismic amplitude into acoustic impedance, and thereafter, combining with the NN method to predict effective porosity. Kumar et al. [12] proposed a method to directly predict porosity using post-stack seismic data and well data, which has been verified by considerable well data and porosity data. Figueiredo et al. [13, 14] proposed two Bayesian inversion methods: the single Gaussian prior distribution and Gaussian mixture model (GMM) as the prior distribution, to obtain the lithofacies, porosities, and acoustic impedances from seismic data and well data.

Combining linearized petrophysical parameters with pre-stack AVO approximation can realize more accurate inversion of rock physical parameters. Grana [15] used the first-order Taylor series approximations to linearize the petrophysical forward model and thereafter used the Bayesian approach to realize linear AVO inversion of rock and fluid properties. Combining an AVO linear approximation and a linearized petrophysical model, Lang and Grana [16] and Guo et al. [17] proposed a linear seismic-petrophysics forward model to realize the inversion of rock attributes such as porosity, clay volume, and water saturation. Combined with Biot-Gassmann’s poroelasticity and Russell’s linear AVO approximation, Li et al. [18] proposed a multitrace Bayesian AVO petrophysics inversion method to realize the inversion of porosity, fluid bulk modulus, and stiff-pore volume fraction; however, the error of linearized petrophysics was larger than in the original model. With the development of the DL method, a few geophysicists have proposed the seismic-petrophysics nonlinear inversion method. Fang and Yang [19] and Yasin et al. [20] proposed a novel hybrid genetic method, combining it with the partially saturated Biot/Squirt (BISQ) model to predict reservoir parameters, including porosity, water saturation, and permeability. Using the GMM and a deep neural network (DNN), Wang et al. [21] proposed the Gaussian mixture model deep neural network (GMM-DNN) method to predict porosity.

The relationship between porosities and seismic attributes is not simple linear, which is usually solved by nonlinear inversion methods. Simulated annealing (SA), which is a global optimization algorithm to simulate the physical annealing process, can be used to solve nonlinear inversion problems. Kirkpatrick et al. [22] and Cerny and Dreyfus [23] first proposed this method, and Bertsimas and Tsitsiklis [24] further discussed its convergence characteristics and applications. Based on different cooling characteristics and convergence conditions, its deformations were developed [25, 26]—very fast simulated annealing (VFSA) is one of them [27, 28]. Combining with the exact Zoeppritz equation and a fully petrophysical model, Guo et al. [29] and Gholami et al. [30] realized the combined inversion of elastic parameters and rock-physics parameters based on the fast simulated annealing (FSA) method.

In this study, seismic records are first expressed as a function of porosity based on the theories of reflection coefficient approximation of fast P-wave and petrophysics of porous media. A nonlinear objective function for porosity inversion is established, as the nonlinear relationship between seismic records and porosities. The SA method is further used to realize the porosity inversion. In the inversion process, the fast P-wave reflection coefficients, porosities, and interpretation results of well log are introduced into the objective function as the prior constraint information to ensure the validity and stability of the inversion results.

2.1. Forward Model

2.1.1. Fast P-Wave Reflection Coefficients

Following Zhou et al. [31], the analytical expression of the reflection coefficient of fast P-wave perpendicularly incident at the interface of porous media saturated with nonviscous fluids can be expressed as
(1)rIps21Ips11+X2X1+Y2Y1Ips21+Ips11+X2+X1+Y2+Y1,
where r denotes the fast P-wave reflection coefficients. This expression consists of three parts with definite geophysical significance: one part includes the fast P-wave impedance terms, Ips11 and Ips21, generated only by rock skeleton; the other refers to the fast P-wave impedance coupling terms, X1 and X2, produced by rock skeleton and pore fluids; the last presents the fast P-wave impedance terms, Y1 and Y2, generated only by pore fluids. A detailed description of related parameters is shown in the appendix.
A certain relationship exists among the bulk (and shear) modulus of rock skeleton, the bulk (and shear) modulus of rock matrix, and effective porosity, which can be quantitatively described via Equations (2) and (3) [32, 33].
(2)Kd=1.0ϕKs1.0+ckϕ,(3)μd=1.0ϕμs1.0+cμϕ.

Here, ϕ is the effective porosity, and Kd and μd represent the bulk and shear moduli of rock skeleton, respectively. Ks and μs denote the bulk and shear moduli of rock matrix, respectively, which can be computed via Voigt-Reuss-Hill (VRH) averaging [34]. ck and cμ are the adjustment parameters of porous media, which can be obtained from P-wave velocity, S-wave velocity, and density in well data.

2.1.2. Fast P-Wave Seismic Records

When a fast P-wave incident is vertically at an interface of porous media saturated with pore fluids, two kinds of P-waves, fast and slow, will be reflected and transmitted. The slow P-wave attenuates and disappears rapidly as it propagates in porous media; so, the seismic records received by detectors are, in fact, the reflected fast P-wave. The effective porosity ϕ is a function of time t as ϕt, and the fast P-wave reflection coefficient is a function of the effective porosity as rϕt according to Equations (1), (2), and (3). Following the convolution theory, the fast P-wave seismic records dϕt can be expressed as
(4)dϕt=wtrϕt+nt.
Equation (4) can be discretized and written as the vector forms
(5)dN×1=GN×NrN×1+nN×1.

Here, wt is a seismic wavelet, and nt denotes the noise. The subscript N×1 denotes an N-dimensional vector, and N×N represents an N-dimensional matrix. The sign N is the sampling point of a seismic record, and “” denotes the convolution operator.

2.2. Simulated Annealing Inversion

The SA method is based on the principle of solid annealing, which simulates the annealing process of crystals in physics. First, the solid is heated to a high temperature and then allowed to slowly cool. When heating, solid particles become disordered because of the rising temperature and their internal energy increases. However, when cooling, they tend to be ordered, and the temperature of each particle reaches equilibrium state. Finally, the particles reach the ground state at normal temperature, and their internal energy decreases to a minimum value. This method is a global optimization algorithm with asymptotic convergence.

2.2.1. Objective Function

Following the forward process in Equation (5), the objective function for porosity inversion F is
(6)F=dddata22+βRR022+γff022.

Here, f is predicted porosities, R denotes predicted reflection coefficients, and d represents synthetic seismic records. f0 and R0 are the porosities and reflection coefficients calculated from well log data, respectively, and ddata denotes real seismic records. β and γ are constraint parameters of reflection coefficients and porosities, respectively. All vectors in function are N dimensional, and the symbol of 22 denotes the L2 norm.

In Equation (6), the first term refers to the fitting degree between synthetic seismic records and observed seismic data; the second term denotes the similarity between the reflection coefficients calculated from predicted porosities and those calculated from well data; and the third term represents the similarity between predicted porosities and the well log porosities. The latter two terms are well data constraints on inversion results, which provide more prior information for objective function and improve the stability and antinoise of inversion results. The β and γ control the contribution of well log constraints to the inversion results, and adjusting them can improve inversion accuracy.

2.2.2. Metropolis Criterion

Iteratively updating porosities, the inversion method obtains the optimal solution. In the ith iteration, the objective function value Fi corresponding to porosities fi is first calculated. Thereafter, the porosities are updated to fi+1, and the corresponding value of the objective function Fi+1 is calculated. The difference between the objective function values ΔF is
(7)ΔF=Fi+1Fi,i1.
If the condition of ΔF<0 is satisfied, which implies that the porosity renewal direction decreases the objective function value, the porosities will be modified; otherwise, we will further determine whether ΔF satisfies Equation (8). If the condition is satisfied, the porosities modification is still accepted; otherwise, the modification is refused.
(8)0<p<K,(9)p=expΔFkbT,
where p denotes the probability density, whose value is described using the Metropolis criterion [35] in Equation (9). K is a random number, whose range is 0<K<1. kb is a Boltzmann constant, and exp represents an exponential function. T represents the current annealing temperature and is often expressed as an exponential form
(10)Ti=T00.95i1,i1.

Here, T0 is the initial temperature, and Ti is the ith annealing temperature.

2.2.3. Iteration and Update of Solution

The solution of this method is independent of the porosity initial state, but our research experience shows that the iteration number can be reduced and the running speed can be improved when the initial solution is given properly. In this study, the initial porosities f1 can be obtained by well log porosities f0 after multiple smoothing. If an initial solution is given, the renewal porosities Δfi can be generated by using the temperature-dependent Cauchy-like distribution as
(11)fi+1=fi+ξΔfi,i1,(12)Δfi=Tsignq0.51+1T2q11.

Here, ξ is a renewal coefficient, fi denotes the current predicted porosities, fi+1 represents the updated porosities, and its range is fi+10,0.3. q is an N-dimensional random vector, and the symbol “sign” denotes a positive or negative sign.

The advantage of using Cauchy-like distribution over generate renewal vectors is that the search of the optimal solution is conducted in a wide range at high temperature while remaining as close to the current predicted value as possible at low temperature. The Cauchy-like distribution has a flat “tail” that can easily jump out of the local minimum.

2.3. Inversion Process

Assume that there are seismic data ddata, reflection coefficients R0, and porosities f0 obtained from well log. The effective porosity inversion process using the SA method is

  • (1)

    Calculate the adjustment parameters ck and cμ from well log information, and set the constrained parameters β and γ in the objective function; thereafter, give the maximum iteration L, the coefficient ξ, and the initial porosities f1

  • (2)

    Using the porosities of the ith iteration fi, the new porosities fi+1 are calculated using Equation (11), and the objective function values are obtained using Equation (6) before Fi and after Fi+1 for the renewal porosities

  • (3)

    Calculate the difference between the two objective function values ΔF using Equation (7), and judge whether the condition of ΔF<0 is satisfied or not; if satisfied, accept the renewal porosities. Otherwise, continue to judge whether to meet Equation (8); if satisfied, accept the update of porosities. Otherwise, refuse to update

  • (4)

    Determine whether the maximum iteration number L is reached. If reached, then withdraw the cycle, and obtain the optimal porosities. Otherwise, update the current annealing temperature using Equation (10), and return to step 2 to continue updating porosities f

When the method accepts the new solution based on the Metropolis criterion, it accepts not only the better solution but also the worse in a certain range, so as to avoid falling into a local minimum. As the number of iterations increases, the annealing temperature T exponentially decreases, and the method can only accept better deterioration solutions. Finally, the value of T tends to be 0, and the method no longer accepts any deterioration solution.

3.1. Model Test

For testing this method, a single-trace model of sandstone porous media saturated with brine is established. By analyzing the initial solution, constrained parameters, and signal-to-noise ratios (SNRs), the validity, stability, and antinoise of the method are verified.

Figure 1 shows the sandstone porous media model, (a) denotes the porosities of sandstone media, (b) denotes the fast P-wave reflection coefficients, (c) represents the Ricker wavelet with dominant frequency of 30 Hz, (d) refers to synthetic seismic records, and (e) represents initial porosities. The sampling rate and model length are 2 ms and 0.6 s, respectively. Other parameters required by the model include the bulk modulus, shear modulus, and density of rock matrix, which are 38GPa, 44GPa, and 2650kg/m3 [36], respectively; the permeability and porosity factor of porous media are 0.2darcy and 2.0, respectively; the adjustment parameters, ck and cμ, are constant at 6.0; and the bulk modulus and density of brine are 3.0GPa and 1050kg/m3 [37], respectively.

3.1.1. Choose Initial Porosities

For analyzing the effect of the initial solution on inversion results, a low-frequency porosity model and a random vector are selected, respectively, as the initial solution for inversion. The convergence of inversion results is described by the energy ratio of the difference between predicted porosities and well log porosities to the well log porosities, that is , hereinafter referred to as porosity energy ratio.

Figure 2 shows the porosity energy ratio for different initial solutions, and the red and blue solid lines represent the selections of a low-frequency porosity model and a random vector as the initial solutions, respectively. Figures 35 show inverted porosities, inverted reflection coefficients, and reconstructed seismic records with different initial solutions, respectively. Figures (a), (b), (c), and (d) represent the 100th, 200th, 300th, and 400th iterations, respectively. Green solid lines denote the original porosities, and the red and blue dotted lines represent the inversion or reconstruction results with a low-frequency porosity model and a random vector as the initial solution, respectively.

Figures 35 show that our method can inverse ideal porosities with a random vector as the initial solution. The inverted reflection coefficients (Figure 4) and reconstructed seismic records (Figure 5) are, respectively, in good agreement with the original models, which shows that our proposed method is applicable to porosity inversion without the initial solution, and the results are effective and stable. However, when a random vector is chosen as the initial solution, even if the Metropolis criterion is adopted to avoid inversion results falling into a local minimum solution, the method still requires more iterations and the inverted speed is relatively slow. When a low-frequency porosity model is selected as the initial solution, our inversion results can be quickly converged to the optimal solution, as shown in Figure 2.

3.1.2. Determine Constraint Parameters

Two constraint parameters, β and γ, exist in the objective function, which determine the relative contributions among the seismic records, well log reflection coefficients, and well log porosities to the inversion results. Here, we analyze the influence of different constraint parameter values on porosity inversion results after the 600th iteration, and the specific contents are as follows.

  • (1)

    When the parameter γ is 0 and parameter β is 0.005, 0.01, 0.05, and 0.5, the porosity inversion results are shown in Figure 6; the corresponding porosity energy ratios are shown in Figure 7. The red, green, blue, and black solid lines in Figure 7 represent the porosity energy ratios with setting parameter β as 0.005, 0.01, 0.05, and 0.5, respectively; and (a), (b), (c), and (d) in Figure 6 denote, respectively, the corresponding porosity inversion results. In Figure 6, the blue solid line is the original porosity model, and the red dotted lines represent porosity inversion results

  • (2)

    When parameter β is 0 and parameter γ is 0.005, 0.01, 0.1, and 0.5, the porosity inversion results are shown in Figure 8; the corresponding porosity energy ratios are shown in Figure 9. The red, green, blue, and black solid lines in Figure 9 represent the porosity energy ratios with setting parameter γ as 0.005, 0.01, 0.1, and 0.5, respectively; (a), (b), (c), and (d) in Figure 8 denote, respectively, the corresponding porosity inversion results. In Figure 8, the blue solid line denotes the original porosity model, and the red dotted lines represent porosity inversion results

  • (3)

    When parameter β is 0.3 and parameter γ is the same as in Figure 8, the porosity inversion results are shown in Figure 10; the corresponding porosity energy ratios are shown in Figure 11. The meanings of all solid and dashed lines in Figures 10 and 11 are the same as those in Figures 8 and 9, respectively

According to Figures 7, 9, and 11, the porosity energy ratios decrease as the number of iterations increase under different constraint parameters, and the porosity inversion results converge gradually. Our research shows that if the constraint term of well reflection coefficients exist only in the objective function (parameter γ is 0), the inversion results are poor because of lack of low-frequency information (Figure 6); whereas, if the well porosity constraint term is only retained (the parameter β is 0), the solutions are unstable because of losing the constraint of intermediate variable (Figure 8); good solutions can be obtained only under both constraints of well porosities and well reflection coefficients (as shown in Figure 10, the low-frequency components of the solutions are gradually compensated with the increase of well porosity constraint parameter γ, and finally we obtain better inversion results).

3.1.3. Noise Suppression Analysis

As real seismic data contains noise, Gaussian white noise is added to the model to verify the antinoise of our proposed method. Figure 12 shows the results of inverted porosities and reconstructed seismic records with different signal-to-noise ratios (SNRs); the corresponding porosity energy ratios are shown in Figure 13. The red and blue solid lines in Figure 13 represent the porosity energy ratios corresponding to the SNRs being 5 and 2, respectively. In Figure 12, (a) and (b) are porosity comparisons between the original model (blue solid lines) and inversion results (red dotted lines) corresponding to the SNRs being 5 and 2, respectively, and (c) and (d) are the corresponding comparisons between the original records (blue solid lines) and reconstructed records (red dotted lines), respectively.

According to Figures 12 and 13, the solutions are still in good agreement with the original model at the condition of low SNRs after the 600th iteration, which shows that our method is effective, stable, and robust.

3.2. Real Data Inversion

3.2.1. Inversion Results

Sandstone reservoir data are intercepted to verify the correctness and practicability of our method. Figure 14 shows the real seismic data and inversion results; (a) denotes the real seismic profile through A well with time-domain depth of 3.2 s-4.1 s, the sampling interval of 4 ms, and CDP of 400; (b) represents a porosity low-frequency model, which is constructed using well log porosity, horizon, and seismic data; and (c) is the porosity inversion results after the 2,100th iterations with constrained parameters of reflection coefficient and porosity being 0.3 and 0.6, respectively. Figure 15 shows the porosity energy ratios with the variation of seismic traces; (a), (b), and (c) are the ratios after the 700th, 1,400th, and 2,100th iterations, respectively.

3.2.2. Prior Information

Unlike other inversion methods based on the single-phase medium theory, our method based on the porous medium theory requires more lithologic and pore fluid parameters. The results of well interpretation and the rock physics test are added to well log as prior information, and the parameter models related to rock matrix and pore fluids, respectively, are established in this study. Figure 16 shows the parameter models of rock matrix and pore fluids; (a) is the bulk modulus of rock matrix, (b) represents bulk modulus of pore fluids, (c) denotes the density of rock matrix, (d) is the density of pore fluids, (e) represents the shear modulus of rock matrix, and (f) denotes the viscous coefficient of pore fluids. For calculating the adjustment parameters in Equations (2) and (3), the low-frequency models are established using well log data. Figure 17 shows well log parameter models; (a), (b), (c), and (d) represent the low-frequency models of P-wave velocity, SV-wave velocity, average density, and P-wave impedance, respectively.

For real data inversion, although the inverted resolution will be reduced by properly decreasing the constrained parameters of well data, it can increase the contribution of seismic data to the solutions and make the results more dependent on real data and more reliable. Our research shows that the inversion results converge to the optimal solution after the 2,100th iterations as shown in Figure 15(c), and the ideal porosity inversion results which are in good agreement with the well log porosity can be obtained as shown in Figure 14(c). The inversion results provide guidance for the prediction of lithology and pore fluids in sandstone reservoirs and give theory and method support for the exploration and development of oil and gas fields.

Based on the analytical expression of the fast P-wave reflection coefficient with a fast P-wave perpendicularly incident at an interface of porous media saturated with a nonviscous fluid, the effective porosity inversion is realized using the SA method. The closed-form expression can be described by lithological parameters and pore fluid parameters. As the bulk and shear moduli of rock skeleton are the functions of porosity and rock matrix moduli, a nonlinear quantitative relationship between seismic records and porosity is established. The curves and interpretation result of well data, the parameters of pore fluids, and other information are introduced into nonlinear porosity inversion as prior constraints. In sum, our method can realize the direct quantitative inversion of porosity, which is more suitable for the characterization, description, and prediction of sandstone reservoirs. As the inversion process depends on initial models to some extent, the inversion effect at the fault is poor, which is a problem that we need to solve in the next research.

Appendix

Relevant Parameter Expressions Involved in the Equation (1)

The analytical expression of the fast P-wave reflection coefficient is given in equation (1), and the specific parameters and their physical meanings are described as
(A.1)X1=α2Ipf22γ22Ipf12+Ipf22Ips11α1Ipf11γ22Ipf12+Ipf22Ips22,(A.2)X2=α1Ipf12γ12Ipf12+Ipf22Ips21α2Ipf21γ12Ipf12+Ipf22Ips12,(A.3)Y1=α1α1α21+α2γ22Ipf22Ipf12+Ipf22Ipf11,(A.4)Y2=α2α2α11+α1γ12Ipf12Ipf12+Ipf22Ipf21,(A.5)Ipslj=Kl+4/3μlvplj,l=1,2,j=1,2,(A.6)Ipflj=Mlvplj,l=1,2,j=1,2.

Here, Kl,μl,l=1,2 denote the bulk and shear moduli of rock skeleton, respectively. vplj,l=1,2,j=1,2 represent the fast or slow P-wave velocity, and α1,Ml,l=1,2 are the Biot parameters. γlj,l=1,2,j=1,2 denote the ratios of displacement potential amplitude of pore fluids relative to rock skeleton to that of rock skeleton. The subscript l=1,2 represents, respectively, the upper and lower porous media, and j=1,2 denotes the fast and slow P-waves, respectively.

Data associated with this research are confidential and cannot be released.

The authors declare no conflicts of interest.

This study was financially supported by the Natural Science Foundation of Sichuan Province of China (Grant No. 2022NSFSC1140), Open Fund (PLC20211101) of State Key Laboratory of Oil and Gas Reservoir Geology and Exploitation (Chengdu University of Technology). The authors are also grateful to the anonymous reviewers for their constructive comments.

Exclusive Licensee GeoScienceWorld. Distributed under a Creative Commons Attribution License (CC BY 4.0).