Prestack seismic facies-controlled joint inversion of reservoir elastic and petrophysical parameters for sweet spot prediction

Abstract

The elastic and petrophysical parameters of a reservoir can be directly applied to lithology prediction and fluid identification. Existing seismic joint inversion methods for estimating the elastic and petrophysical parameters of a reservoir are primarily based on either the Gassmann equation, with which these parameters are inverted from prestack seismic data through stochastic optimization methods, or Wyllie’s modified equation, with which these parameters are inverted from poststack seismic data using deterministic optimization methods. Regardless of the stochastic or deterministic inversion of reservoir parameters, only the microscopic connection between parameters is considered, without considering the constraints of the macroscopic geological background. The purpose of this work is to develop a strategy for estimating the elastic and petrophysical parameters of a reservoir based on the Gassmann equation and seismic facies using deterministic seismic prestack inversion. We employ the Gassmann equation to construct the relationship between the prestack seismic data and petrophysical parameters and use a low-pass filter matrix to obtain the combination of seismic facies and seismic inversion and improve the robustness of the inversion results. We treat the joint posterior probability of elastic and petrophysical parameters as the objective function under a Bayesian framework by expanding the objective function with the Taylor formula; the joint equations composed of elastic and petrophysical parameters are obtained through differentiation. The conjugate gradient method is subsequently used to find the optimal solutions for the P-wave velocity, S-wave velocity, density, porosity, water saturation, and clay content. Theoretical calculations and actual data inversion results prove the feasibility and applicability of the method.

Introduction

Seismic inversion is an important technique that is employed to obtain the elastic and petrophysical parameters of reservoir lithologies for reservoir prediction and description. With the increasing difficulty of exploration and development, pure elastic parameter inversion cannot fully meet the needs of these predictions and descriptions. As a consequence, using prestack seismic data to carry out the elastic and petrophysical parameter inversion has become a common interest of frontier science.

There are two kinds of methods used for inverting elastic and petrophysical parameters, namely the two-step inversion and joint inversion methods (Bosch et al., 2010). Two-step inversion first employs seismic inversion to obtain the elastic parameters and then estimates the petrophysical parameters, through the deterministic or statistical relationship between the elastic and petrophysical parameters (Avseth et al., 2001; Buland et al., 2008; Grana and Della Rossa, 2010; Lerche and Noeth, 2002; Lessenger and Cross, 1996; Mukerji et al., 2001; Saltzer et al., 2005; Thomsen and Noeth, 2001). Due to the complexity and variability of reservoirs, the elastic parameters of different lithologies often have some range of overlap, and the error caused by the inversion of elastic parameters may lead to further error in estimations of reservoir petrophysical parameters. These shortcomings limit the extensive application of the two-step inversion method in reservoir prediction and description (Bachrach, 2006; Sengupta and Bachrach, 2007).

The joint inversion method combines geological a priori information of reservoir parameters, rock physics models, and seismic forward models and obtains reservoir elastic and petrophysical parameters directly from the seismic data by iterating the applied forward modeling. Joint inversion usually adopts a stochastic optimization method (Bosch et al., 2007; Loures and Moraes, 2006). The model parameters and seismic data are regarded as random variables, and a large number of selectable models are randomly generated and calculated, following which the calculated results are compared with the real seismic data to determine whether to accept the model based on the predetermined information (Eidsvik et al., 2004; Rimstad and Omre, 2010). The advantage of stochastic inversion is that the inversion process is not completely dependent on the initial model and theoretically will not fall into local extrema. The stochastic inversion can be used to predict reservoir elastic and petrophysical parameters, which can either generate a great deal of cognition about the correlation between seismic data and petrophysics or calculate the marginal probability of the petrophysical parameters and quantitatively evaluate the uncertainty of the inversion results (Grana et al., 2012; Gunning and Glinsky, 2007; Kjønsberg et al., 2010; Larsene et al., 2006). However, the shortcomings of stochastic inversion are also obvious, as the computational workload is large and the efficiency is low, which makes it difficult to fully meet the requirements of large-scale inversion.

Joint inversion methods based on deterministic inversion often use the mathematical equation or algebra method and obtain the optimal solution of the inversion problem through repeated iteration (González et al., 2007; Spikes et al., 2007). The advantage of the deterministic inversion is that the convergence speed is fast, and it is thus suitable for large-scale production tasks. Bosch (2004) and Bosch et al. (2009) established the connection between petrophysical parameters and wave impedance using Wyllie’s time average equation in a Bayesian theoretical framework and used Newton’s iterative method to conduct petrophysical parameter inversion by using poststack seismic data. Although the above method achieved good results, Bosch et al. (2009) noted two defects of this approach: first, poststack seismic data do not contain information regarding the amplitude versus incident angle, which renders an estimation of the S-wave velocity and density difficult, and it is accordingly difficult to guarantee a reliable estimate of the reservoir parameters; second, Wyllie’s time average equation is not suitable for application to a mudstone-bearing sandstone and thus cannot be used for porosity calculations of gas-bearing sandstones (Elliott and Wiley, 1975; Saleh and Castagna, 2004).

Currently, prestack inversion primarily focuses on the prior distribution of inversion parameters (Alemie and Sacchi, 2011; Buland and Omre, 2003; Downton, 2005; Karimi et al., 2010; Kidney et al., 1992). In fact, the prior distribution of the inversion parameters is only related to the vertical resolution and sparsity of the inversion results. Meanwhile, another important component of the inversion outcome, the low frequency component, cannot be obtained from the seismic data directly and requires additional constraints. There are relatively few studies on the low frequency component constraint term for stable inversion results. Tian et al. (2013) used a point constrained term and a smoothed regular constraint term to improve the antinoise properties of the inversion and supplement the low frequency component of the inversion results. Niepsuj (2012) tried to utilize different initial models to invert actual data originating from the same work area and demonstrated the importance of the initial model, especially the low frequency component, for inversion results. Numerous methods (Douma and Naeini, 2014; Sams and Saussus, 2013; Wang and Huang, 2009) have been suggested to construct an initial model, and researchers have long sought ways to incorporate log data and seismic facies into initial models. Hunag and Kelkar (1996) and Huang et al. (1995) described acoustic impedance inversion using a heuristic algorithm in which an initial model from log data accelerated the convergence of the final predicted model. Kemper and Gunning (2014) specified initial values with different seismic facies for the joint inversion. Li et al. (2015) demonstrated an achievable higher resolution for inversion of seismic facies constraints by using the conventional tectonic interface and seismic facies interface to invert the actual data. Point constrained and smoothed regular constraint and other constrained methods can stabilize the inversion results of the low frequency components, but inappropriate constraint weighting options may reduce the resolution of the indirect conversion parameters, Therefore, how to accurately invert low-frequency components is a problem to be solved.

In this paper, we employ the Gassmann equation to construct the relationship between prestack seismic data and the corresponding petrophysical parameters, and we treat the joint posterior probability of the elastic and petrophysical parameters as the objective function under the Bayesian framework. We use the Taylor approximation to expand the gradient of the objective function to obtain the joint equations representing the elastic and physical properties of the reservoir, wherein the gradient of the elastic parameters to the physical parameters is represented by the difference form. By constructing a matrix expression of low-pass filtering, a seismic facies constraint is applied to the seismic inversion, which establishes a robust low frequency trend. The conjugate gradient method is then used to invert the P-wave velocity, S-wave velocity, density, porosity, water saturation, and clay content. Compared with smoothed regular constraints, seismic facies constraints are relatively insensitive to constrained weight coefficients and can obtain an accurate low frequency trend. Compared with other petrophysical models, the Gassmann equation can better establish the relationship between the effective elastic parameters of fluid-saturated rocks and pore fluids. Compared with the stochastic optimization method, the deterministic optimization method improves the inversion velocity on the basis of guaranteeing the validity and accuracy of the method. Theoretical calculations in addition to real data inversion result from the Sangtamu oilfield prove the feasibility of the method.

Methodology

Seismic forward modeling

Assuming a small change in the elastic parameters of the reflection interface, Aki and Richard (1980) approximated the reflection coefficient using the relative change of velocity and density. Based on the Aki and Richard (1980) approach, Fatti (1994) constructed an approximate equation expressed as the relative wave impedance

$$R(\theta )=\frac{1}{2}(1+{\tan }^2\theta )\frac{\Delta I_P}{I_P}-4{\gamma }^2{\sin }^2\theta \frac{\Delta I_S}{I_S}-(\frac{1}{2}{\tan }^2\theta -2{\gamma }^2{\sin }^2\theta )\frac{\Delta \rho }{\rho }$$

(1) where $R(\theta )$ is the reflection coefficient of the P-wave; θ is the average of the incident angle and the transmission angle; I_P, I_S, and ρ are the averages of the P-wave impedance, S-wave impedance, and density of the reflection interface, respectively; and $\Delta I_P$ , $\Delta I_S$ , and $\Delta \rho$ are the differences of the P-wave impedance, S-wave impedance, and density of the reflection interface, respectively. $\gamma =V_S/V_P$ represents the velocity ratio of the S-wave and P-wave. In the case of n interfaces and an offset m, the prestack seismic forward model can be expressed as

$$g(x)=G_{\mathrm{mm}\times 3n}{x}_{3n\times 1}+n_{{\mathrm{mm}\times }_1}$$

(2) where $g(x)$ is the prestack seismic forward model, $G$ is the prestack seismic forward model operator, $n$ is the noise vector, and $x$ represents the elastic parameter vector to be estimated which is composed of the logarithm of the P-wave impedance, the logarithm of the S-wave impedance, and the logarithm of the density.

Rock physics model of joint inversion

The elastic properties of rocks form the basis of the propagation of seismic waves in elastic media. When the same lithology contains a different fluid, the elastic properties will be changed correspondingly; this is the physical basis of fluid replacement, which is typically used to predict the velocity of fluid-saturated rocks. In this paper, the bulk modulus and shear modulus of a sandstone skeleton is first calculated by applying the Hertz–Mindlin contact theory and a Hashin–Shtrikman boundary condition. Then, the Gassmann equation is employed to estimate the corresponding moduli for fluid-saturated rocks, as well as establish the relationship between the petrophysical parameters $[\phi ,S_w,\text{clay}{]}^T$ and elastic parameters $[L_p,L_s,L_d{]}^T$ , and is subsequently used to calculate the P- and S-wave velocities and density, where φ, S_w, and $\text{clay}$ are the porosity, water saturation, and clay content of the rock, respectively, and L_p, L_s, and L_d are the logarithms of the P-wave impedance, S-wave impedance, and density, respectively.

The bulk modulus and shear modulus of the sandstone skeleton beyond the grain contact can be estimated using an unconsolidated sandstone model. Assume that the equivalent sphere grains are random packing in close proximity, the porosity of the skeleton of the initial unconsolidated sandstone is ${\phi }_0$ with a value of 0.36, and the average coordination number is 9 (Mavko et al., 1998; Murphy, 1982). The Hertz–Mindlin theory (Mavko et al., 1998; Mindlin, 1949) provides the effective bulk modulus and shear modulus of the equivalent arbitrarily arranged spherical grain under hydrostatic pressure

$$K_{\mathrm{HM}}=[\frac{C^2(1-{\phi }_0{)}^2{\mu }_{\mathrm{mat}}^2}{18{\pi }^2(1-\upsilon {)}^2}P{]}^{\frac{1}{3}}$$

(3)

$${\mu }_{\mathrm{HM}}=\frac{5-4\upsilon }{5(2-\upsilon )}[\frac{3C^2(1-{\phi }_0{)}^2{\mu }_{\mathrm{mat}}^2}{2{\pi }^2(1-\upsilon {)}^2}P{]}^{\frac{1}{3}}$$

(4) where $K_{\mathrm{HM}}$ and ${\mu }_{\mathrm{HM}}$ are the effective bulk modulus and shear modulus of the rock skeleton, respectively; C is the average coordination number; ${\phi }_0$ is the porosity of the unconsolidated sandstone skeleton; P is the hydrostatic pressure; υ and ${\mu }_{\mathrm{mat}}$ are the Poisson’s ratio and shear modulus of the grain, respectively.

Based on the calculation of the effective bulk and shear moduli of the rock skeleton, we can use the modified Hashin-Shtrikman (Mavko et al., 1998) upper boundary theory to solve the effective modulus of the sandstone skeleton with a different porosity

$$K_{\mathrm{dry}}=[\frac{\phi /{\phi }_0}{K_{\mathrm{HM}}+\frac{4}{3}{\mu }_{\mathrm{mat}}}+\frac{1-\phi /{\phi }_0}{K_{\mathrm{mat}}+\frac{4}{3}{\mu }_{\mathrm{mat}}}]-\frac{4}{3}{\mu }_{\mathrm{mat}}$$

(5)

$${\mu }_{\mathrm{dry}}=[\frac{\phi /{\phi }_0}{{\mu }_{\mathrm{HM}}+\frac{1}{6}{\mu }_{\mathrm{mat}}}+\frac{1-\phi /{\phi }_0}{{\mu }_{\mathrm{mat}}+\frac{1}{6}\xi {\mu }_{\mathrm{mat}}}]-\frac{1}{6}\xi {\mu }_{\mathrm{mat}}$$

(6)

$$\xi =\frac{9K_{\mathrm{mat}}+8{\mu }_{\mathrm{mat}}}{K_{\mathrm{mat}}+2{\mu }_{\mathrm{mat}}}$$

(7)

$$K_{\mathrm{mat}}=\frac{1}{2}(\frac{{\mathrm{clayK}}_c+(1-{\mathrm{clayK}}_d)}{1-\phi }+\frac{1-\phi }{\mathrm{clay}/K_c+(1-\mathrm{clay})/K_d})$$

(8)

$${\mu }_{\mathrm{mat}}=\frac{1}{2}(\frac{\mathrm{clay}{\mu }_c+(1-\mathrm{clay}{\mu }_d)}{1-\phi }+\frac{1-\phi }{\mathrm{clay}/{\mu }_c+(1-\mathrm{clay})/{\mu }_d})$$

(9) where φ is the porosity; $K_{\mathrm{dry}}$ and ${\mu }_{\mathrm{dry}}$ are the effective bulk modulus and shear modulus of the sandstone skeleton with a different porosity, respectively; $K_{\mathrm{mat}}$ is the bulk modulus of the solid grain; K_c and K_d are the bulk modulus of mudstone and sandstone, respectively; ${\mu }_c$ and ${\mu }_d$ . are the shear modulus of mudstone and sandstone, respectively.

The Gassmann–Biot theory (Mavko et al., 2009) is used to predict the resultant increase in the modulus of the fluid-saturated rock

$$\frac{K_{\mathrm{sat}}}{K_0-K_{\mathrm{sat}}}=\frac{K_{\mathrm{dry}}}{K_0-K_{\mathrm{dry}}}+\frac{K_{\mathrm{fl}}}{\phi (K_0-K_{\mathrm{fl}})}$$

(10)

$${\mu }_{\mathrm{sat}}={\mu }_{\mathrm{dry}}$$

(11) where φ is the lithology porosity; $K_{\mathrm{sat}}$ and ${\mu }_{\mathrm{sat}}$ are the effective bulk and shear moduli of the fluid-saturated rock, respectively; $K_{\mathrm{dry}}$ and ${\mu }_{\mathrm{dry}}$ are the effective bulk and shear moduli of dry rock, respectively; K₀ is the effective bulk modulus of the rock minerals; and $K_{\mathrm{fl}}$ is the effective bulk modulus of the pore fluid, which can be expressed as

$$K_{\mathrm{fl}}=(\frac{S_w}{K_{\mathrm{brine}}}+\frac{1-S_w}{K_{\mathrm{hc}}}{)}^{-1}$$

(12) where S_w is the water saturation, $K_{\mathrm{brine}}$ is the effective bulk modulus of the formation water, and $K_{\mathrm{hc}}$ is the effective bulk modulus of oil or gas.

Extrapolation calculation controlled by seismic facies

Seismic facies represent sedimentary facies interpreted from seismic data, as variations among sedimentary facies induce variations in seismic waveform characteristics (Huang et al., 2016). In order to accurately describe subsurface sedimentary features, we introduce facies interfaces into the seismic data, which accordingly divide the target strata into several different assumed sedimentary units. Different seismic facies display different ranges of velocity, impedance, and other reservoir parameters. Therefore, seismic facies can control seismic inversion with their variable low frequency tendency. In addition, seismic facies interfaces also control the time window size of the inversion parameters.

The low frequency component of the inversion output plays an important role in reservoir prediction and fluid identification. Incorrectly estimated low frequency components may lead to erroneous estimation results. The correctness of the low frequency component primarily depends on the following two aspects: a relatively accurate initial model and the robustness of the inversion algorithm, which can ensure that the low frequency component of the inversion result does not deviate from the low frequency component of the initial model.

In the proposed work, the low-pass filtering process is expressed as a matrix, and the L₂ norm of the difference between the low frequency components of the inversion result and the logging data is minimized to ensure the accuracy of the low frequency information of the inversion result. Therefore, the mathematic expression of the seismic facies-constrained function (i.e. the low frequency soft constraint term) is

$${\mathrm{Const}}_{\mathrm{LowFre}}=\overline{x}-f(x)={\alpha }_{\text{IP}}(S_{\text{IP}}{H}_{\text{IP}}x-\ln L_{\text{IP}}{)}^{\text{T}}(S_{\text{IP}}{H}_{\text{IP}}x-\ln L_{\text{IP}})+{\alpha }_{\text{IS}}(S_{\text{IS}}{H}_{\text{IS}}x-\ln L_{\text{IS}}{)}^{\text{T}}(S_{\text{IS}}{H}_{\text{IS}}x-\ln L_{\text{IS}})+{\alpha }_{\text{D}}(S_{\text{D}}{H}_{\text{D}}x-\ln L_{\text{D}}{)}^{\text{T}}(S_{\text{D}}{H}_{\text{D}}x-\ln L_{\text{D}})$$

(13) where $\overline{x}$ consists of $L_{\text{IP}}$ , $L_{\text{IS}}$ , and $L_{\text{D}}$ ,which is the low frequency component of the model, $L_{\text{IP}}$ , $L_{\text{IS}}$ , and $L_{\text{D}}$ represent the low frequency component of P-wave impedance, shear wave impedance, and density; $f(x)$ is the low-pass filter function; ${\alpha }_{\text{IP}}$ , ${\alpha }_{\text{IS}}$ , and ${\alpha }_{\text{D}}$ are the constraint weight coefficients for the P-wave impedance, S-wave impedance, and density terms, respectively; and $S_{\text{IP}}$ , $S_{\text{IS}}$ , and $S_{\text{D}}$ are the low-pass filter matrices of the P-wave impedance, shear-wave impedance, and density terms, respectively, which perform the low-pass filtering of the inversion parameters. $H_{\text{IP}}$ , $H_{\text{IS}}$ , and $H_{\text{D}}$ are the integral matrices of the P-wave impedance, S-wave impedance, and density, respectively. Take the P-wave impedance for example

$$H_{\text{IP}}=[10\dots 00\dots 011\dots 00\dots 0::\dots ::\dots 011\dots 10\dots 0{]}_{3N\times N}$$

(14)

$$\ln L_{\text{IP}}=[\ln \frac{I_{\text{P}2}}{I_{\text{P}1}}\ln \frac{I_{\text{P}3}}{I_{\text{P}1}}\dots \ln \frac{I_{\text{P}(\text{N}+1)}}{I_{\text{P}1}}{]}_{1\times N}^T$$

(15) wherein $S_{\text{IP}}{H}_{\text{IP}}x$ , $S_{\text{IS}}{H}_{\text{IS}}x$ , and $S_{\text{D}}{H}_{\text{D}}x$ represent the low frequency components of the natural logarithms of the inverted results of the P-wave impedance, the S-wave impedance, and the density, respectively. The mathematical meaning of the above equation can be regarded as the minimum L₂ norm of the difference between the low frequency components of the natural logarithms of the inversion result and the model. The physical meaning is to employ the low frequency component of the inversion result to approach the low frequency component of the model to ensure that the inversion results possess the correct low frequency component to achieve the desired seismic facies-controlled constraint. This constraint does not affect the middle and high frequency components of the inversion result. Therefore, the low frequency soft constraint is essentially different from the smooth regular constraint.

The construction of a low-pass filter matrix is the key to a seismic facies-controlled constraint. Taking the P-wave impedance as an example, its mathematical expression is

$$S_{\text{IP}}=E^{\text{T}}{F}^{-1}\mathrm{QFE}$$

(16) where $F$ and $F^{-1}$ are the positive and inverse Fourier transform matrices, $Q$ is the diagonal matrix composed of the Hanning window functions (which, when the window of the inversion objective interval is too short, will produce an insufficient resolution in the frequency domain or which, due to the truncation effect on the data, will cause the low-pass filtering result to be inaccurate), and the matrix $E$ extends the model parameters three times. Although equation (14) is complex, its physical meaning is clear.

Objective function and algorithm optimization of inversion

The prestack seismic data are also known as the observed data, and the unknown model parameters are the elastic parameters and petrophysical parameters. Taking the observed data and model parameters as random variables, the inversion objective function is the conditional probability of the model parameters on the condition of the observation data, and the inversion result is the corresponding model parameter when the conditional probability is a maximum, which can be expressed as

$$\{y,x\}=\arg {\text{Max}}_{y,x}\{p(y,x|d)\}$$

(17) where $d$ is the prestack seismic data; $x$ is the elastic parameters, which a vector consisting of the logarithms of the P-wave impedance, S-wave impedance, and density; and $y$ is the petrophysical parameters, which is a vector consisting of the porosity, water saturation, and clay content.

Based on Bayesian theory, equation (17) can be expressed as the product of the seismic likelihood function and the a priori function of the model

$$\{y,x\}=\arg {\text{Max}}_{y,x}\{p(d|y,x)•p(y,x)\}$$

(18) where $p(d|y,x)$ is the likelihood function and $p(y,x)$ is the a priori function of the model.

The seismic likelihood function is only directly associated with the elastic parameters of the reservoir, and therefore, the seismic likelihood function can be expressed as

$$p(d|y,x)=p(d|x)$$

(19)

Substituting equation (19) into equation (18) yields

$$\{y,x\}=\arg {\text{Max}}_{y,x}\{p(d|x)•p(y,x)\}$$

(20)

Under the condition of knowing the reservoir petrophysical parameters in combination with the rock physical model, the joint distribution of elastic parameters and physical parameters can be expressed as

$$p(y,x)=p(x|y)•p(y)$$

(21) where $P(y,x)$ is the joint distribution of elastic parameters and petrophysical parameters. Substituting equation (21) into equation (20) yields the following

$$\{y,x\}=\arg {\mathrm{Max}}_{y,x}\{p(d|x)•p(x|y)•p(y\}$$

(22)

Assume a priori geological information and that the error of the rock physics and the seismic noise obey a Gaussian distribution. Then, equation (22) can be further expressed as

$$\{y,x\}=\arg {\text{Max}}_{y,x}\{a_1\exp [-\frac{1}{2}(g(x)-d{)}^{\text{T}}{C}_d^{-1}(g(x)-d)]\times a_2\exp [-\frac{1}{2}(x-f(y){)}^{\text{T}}{C}_{x|y}^{-1}(x-f(y))]\times a_3\exp [-\frac{1}{2}(y-\overline{y}{)}^{\text{T}}{C}_y^{-1}(y-\overline{y})]\}$$

(23) where $x$ is the vector of the elasticity parameters; $g(x)$ is the prestack seismic forward model; $y$ is the vector of the petrophysical parameters; $f(y)$ is the petrophysical equation with which to establish the relationship between the petrophysical and elastic parameters; $d$ is the vector of the observed seismic data; $C_d$ is the covariance matrix of the observed seismic data; $C_{x|y}$ is the covariance matrix of the elastic parameters with the petrophysical parameters as the condition; $\overline{y}$ is the mean of the petrophysical parameters; $C_y$ is the covariance matrix of the petrophysical parameters; and a₁, a₂, and a₃ are constant terms.

Appending the seismic phase-controlled constraint term to the objective function, we obtain a new objective function, and equation (23) is correspondingly equivalent to

$$\{y,x\}=\arg {\text{Min}}_{y,x}\{\frac{1}{2}(g(x)-d{)}^{\text{T}}{C}_d^{-1}(g(x)-d)+\frac{1}{2}(x-f(y){)}^{\text{T}}{C}_{x|y}^{-1}(x-f(y))+\frac{1}{2}(y-\overline{y}{)}^{\text{T}}{C}_y^{-1}(y-\overline{y})+\lambda (\overline{x}-f(x){)}^{\text{T}}{C}_x^{-1}(\overline{x}-f(x))\}$$

(24)

Using the Taylor formula to expand the gradient of equation (24) yields the following

$$\nabla S(y_{\text{n}+1},x_{\text{n}+1})-\nabla S(y_{\text{n}},x_{\text{n}})={\nabla }^{2}S(y_{\text{n}},x_{\text{n}})•\Delta (y,x)$$

(25) where ${\nabla }^{2}S(y_{\text{n}},x_{\text{n}})$ is the Hesse matrix of $S(y_{\text{n}},x_{\text{n}})$ , $\nabla S(y_{\text{n}},x_{\text{n}})$ , and $\nabla S(y_{\text{n}+1},x_{\text{n}+1})$ are the gradients of $S(y_{\text{n}},x_{\text{n}})$ and $S(y_{\text{n}+1},x_{\text{n}+1})$ , respectively, and $\Delta (y,x)=(y_{\text{n}+1},x_{\text{n}+1})-S(y_{\text{n}},x_{\text{n}})$ .

Let the $\nabla S(y_{\text{n}+1},x_{\text{n}+1})$ term in equation (25) be equal to zero, which yields the following

$$-\nabla S(y_{\text{n}},x_{\text{n}})={\nabla }^{2}S(y_{\text{n}},x_{\text{n}})•\Delta (y,x)$$

(26) where

$${\nabla }^{2}S(y_{\text{n}},x_{\text{n}})=(C_y^{-1}+F^{\text{T}}{C}_{{x|}_y}^{-1}F-C_{{x|}_y}^{-1}{F}^{\text{T}}-{\mathrm{FC}}_{{x|}_y}^{-1}{C}_{{x|}_y}^{-1})+(0002\lambda E^{\text{T}}{C}_x^{-1}E+G^{\text{T}}{C}_d^{-1}G)$$

(27)

$$-\nabla S(y_{\text{n}},x_{\text{n}})=(C_y^{-1}(\overline{y}-y)-F^{\text{T}}{C}_{x|y}^{-1}(f(y)-x)C_{x|y}^{-1}(f(y)-x)+\lambda E^{\text{T}}{C}_x^{-1}(\overline{x}-f(x))+G^{\text{T}}{C}_d^{-1}(d-g(x)))$$

(28) where $G^{\text{T}}$ is the gradient of the prestack forward model to the elastic parameters, which can be expressed as $G^{\text{T}}=\partial (\mathrm{Gx})/\partial x$ ; $E^{\text{T}}$ is the gradient of the seismic facies-controlling model to the elastic parameter, which can be expressed as $E^{\text{T}}=\partial f(x)/\partial x$ ; $F^{\text{T}}$ is the gradient of the petrophysical parameters to the elastic parameters, which can be expressed as $F^{\text{T}}=\partial f(y)/\partial y$ ; the gradient of the reservoir elasticity parameter to porosity is denoted as $F{\text{'}}_{\phi }$ if the petrophysical equation is linear; and $F{\text{'}}_{\phi }$ is a constant. The petrophysical equation is nonlinear, and its mathematical expression is complex; thus, its direct derivation is difficult. Take the gradient of the reservoir elasticity parameter to porosity as an example; $F{\text{'}}_{\phi }$ can be computed in the form of the difference

$$F{\text{'}}_{\phi }=\frac{f(\phi +\Delta \phi ,S_w,\mathrm{clay})-f(\phi ,S_w,\mathrm{clay})}{\Delta \phi }$$

(29)

Similarly, the gradient of the elastic parameters to the water saturation and clay content can be computed in the same form.

To simplify the ${\nabla }^{2}S(y_{\text{n}},x_{\text{n}})$ term in equation (26) to the upper triangular matrix, both sides of equation (26) are left multiplied by $C_{y,x}=(C_{x|y}{\mathrm{FC}}_{x|y}{C}_{x|y}{F}^{\text{T}}{C}_{x|y}+F^{\text{T}}{C}_{y}F)$ , which yields

$$(C_{x|y}[C_y^{-1}(\overline{y}-y)+F[\lambda E^{\text{T}}{C}_x^{-1}(\overline{x}-f(x))+G^{\text{T}}{C}_d^{-1}(d-g(x))]C_{x|y}{\mathrm{FC}}_y^{-1}(\overline{y}-y)+({\mathrm{FF}}^{\text{T}}+(C_{x|y}+{\mathrm{FC}}_y{F}^{\text{T}})C_{x|y}^{-1})(f(y)-x)+(C_{x|y}+F^{\text{T}}{C}_{y}F)\lambda E^{\text{T}}{C}_x^{-1}(\overline{x}-f(x))+G^{\text{T}}{C}_d^{-1}(d-g(x)))=(I{\mathrm{FC}}_{x|y}(2\lambda E^{\text{T}}{C}_x^{-1}E+G^{\text{T}}{C}_d^{-1}G)0I+(C_{x|y}+{\mathrm{FC}}_y{F}^{\text{T}})(2\lambda E^{\text{T}}{C}_x^{-1}E+G^{\text{T}}{C}_d^{-1}G))(\Delta y\Delta x)$$

(30) then we obtain the decoupling form of the elastic parameter and the physical parameter correction.

Equation (30) is thus equivalent to

$$[I+(C_{x|y}+{\mathrm{FC}}_y{F}^{\text{T}})(2\lambda E^{\text{T}}{C}_x^{-1}E+G^{\text{T}}{C}_d^{-1}G)]•\Delta x=C_{x|y}{\mathrm{FC}}_y^{-1}(\overline{y}-y)+({\mathrm{FF}}^{\text{T}}+(C_{x|y}+{\mathrm{FC}}_y{F}^{\text{T}})C_{x|y}^{-1})(f(y)-x)+(C_{x|y}+F^{\text{T}}{C}_{y}F)[{\mathrm{FC}}_y{F}^{\text{T}}+C_{x|y}]G^{\text{T}}{C}_d^{-1}[d-g(x)]$$

(31)

$$\Delta y=C_{x|y}[C_y^{-1}(\overline{y}-y)+F[\lambda E^{\text{T}}{C}_x^{-1}(\overline{x}-f(x)-2\Delta x)+G^{\text{T}}{C}_d^{-1}(d-g(x)-\Delta x)]$$

(32) equations (31) and (32) are the basic iterative equation of elastic parameters and petrophysical parameters inversion.

Let

$$N=C_{x|y}{\mathrm{FC}}_y^{-1}(\overline{y}-y)+({\mathrm{FF}}^{\text{T}}+(C_{x|y}+{\mathrm{FC}}_y{F}^{\text{T}})C_{x|y}^{-1})(f(y)-x)+(C_{x|y}+F^{\text{T}}{C}_{y}F)[{\mathrm{FC}}_y{F}^{\text{T}}+C_{x|y}]G^{\text{T}}{C}_d^{-1}[d-g(x)]$$

and $M=[I+(C_{x|y}+{\mathrm{FC}}_y{F}^{\text{T}})(2\lambda E^{\text{T}}{C}_x^{-1}E+G^{\text{T}}{C}_d^{-1}G)]$

Then equation (31) can be rewritten as

$$M\Delta x=N$$

(33)

L₂ norm has good ability to suppress Gauss noise and super Gauss noise, and can improve the robustness of the algorithm in noisy environment (Yuan et al., 2015). Define the target function for $\Delta x$ as follows

$$J(\Delta x)=\Vert N-M\Delta x{\Vert }^2$$

(34) then the first derivative of the objective function is

$$J\text{'}(\Delta x)=N-M\Delta x$$

(35)

Calculating the first derivative at the initial model $\Delta x_0$ and construct the orthogonal vector a_i constructing conjugate vectors b_i based on orthogonal vectors a_i

$$\{{\alpha }_j=-\frac{a_j^T{a}_j}{a_{j-1}^T{a}_{j-1}}{b}_0=a_0{b}_i=a_i+\sum _{j=1}^{i-1}{\alpha }_j{b}_j$$

(37) where ${\alpha }_j$ is the correction factor of conjugate vector.

Based on equation (37), establishing the iterative updating model for $\Delta x$

$$\{\Delta x_{i+1}=\Delta x_i+{\beta }_{i}b{\beta }_i=-\frac{[N-M\Delta x{]}^T{b}_i}{b_i^{T}M{b}_i}$$

(38) where ${\beta }_i$ is the optimization step length. By using the conjugate gradient method to solve equation (31), the valve of $\Delta x$ can be obtained, and thereafter, by substituting $\Delta x$ into equation (32), $\Delta y$ can be calculated.

Synthetic data example and real data application

To test the feasibility of the prestack seismic facies-controlled joint inversion of reservoir elastic and petrophysical parameters, we utilize synthetic seismic data to test and apply the technique within a real data example. The logging and seismic data used for the synthetic data example and practical application are from the Sangtamu Oilfield located in western China.

Synthetic data example

For the synthetic seismic data test, the logging curve interpretation is first used to obtain the reservoir petrophysical parameter curve (Figure 1(d) to (f)). Subsequently, the petrophysical parameters are converted into elastic parameters by the petrophysical equation, and the Gaussian noise with mean 0 and variance $\varepsilon 1$ will be added to the elastic parameters (Figure 1(a) to (c)). Finally, angle gather data are synthesized using the elastic parameters (Figure 2(a)), and Gaussian noise with mean 0 and variance $\varepsilon 2$ will be added to the angle gather data (Figure 2(b)). During the construction of the synthetic seismic record, multiples and transmission loss are not taken into account, and a Ricker wavelet with a main frequency of 30 Hz is used.

Figure 1.

Reservoir parameters derived from well-log data. (a) P-wave velocity, (b) S-wave velocity, (c) density, (d) porosity, (e) water saturation, and (f) clay content.

Figure 2.

(a) The angle gather generated by well-log data and (b) the angle gather affected by Gaussian noise.

Figure 3 illustrates the inversion result of the two-step inversion, wherein the blue lines represent the well-log data, and the red lines represent the inverted results. Initially, the elastic parameters are obtained through inversion of the angle gather data combined with Gaussian noise, and then, the porosity, water saturation, and clay content are estimated according to the relationship between the elastic parameters and petrophysical parameters. As shown in Figure 3, there is a certain error between the elastic parameters obtained from the two-step inversion method and the real data in the Gaussian noise environment; these errors caused by the inversion of the elastic parameters further lead to errors in the estimation of the reservoir physical parameters.

Figure 3.

The inversion result obtained by applying the two-step forward inversion method (the blue lines are the well-log data, the black lines represent the initial model, and the red lines represent the estimated values). (a) P-wave velocity, (b) S-wave velocity, (c) density, (d) porosity, (e) water saturation, and (f) clay content.

Figures 4 and 5 demonstrate the inversion results of the proposed method with seismic facies weight values of 10³ and 10⁶, respectively. Compared with the two-step forward inversion, the proposed method has relatively high resolution; specifically, the resolution of the inverted petrophysical parameters is improved significantly. Comparing Figures 4 and 5, it is apparent that increasing the constraint weight does not affect the inversion results. This suggests that the algorithm is more robust and is relatively insensitive to the constraint weight.

Figure 4.

The inversion result obtained by applying the proposed method with a seismic facies weight value of 103 (the blue lines are the well-log data, the black lines represent the initial model, and the red lines represent the estimated values). (a) P-wave velocity, (b) S-wave velocity, (c) density, (d) porosity, (e) water saturation, and (f) clay content.

Figure 5.

The inversion result obtained by applying the proposed method with a seismic facies weight value of 106 (the blue lines are the well-log data, the black lines represent the initial model, and the red lines represent the estimated values). (a) P-wave velocity, (b) S-wave velocity, (c) density, (d) porosity, (e) water saturation, and (f) clay content.

Figure 6(a) demonstrates the seismic residual of the proposed method and the two-step inversion evaluation with iterations, we can find out that the proposed method has good steady-state performance and faster convergence speed, Figure 6(b) demonstrates the correlation coefficient of real seismic record and seismic record generated by the inverted parameters evaluation with iterations, we can clearly find that the correlation coefficient of the two methods increases rapidly with the increase of the number of iterations, and the proposed method has a higher speed and correlation coefficient than the two-step inversion, which reflects the superiority of the proposed method on the other hand.

Figure 6.

Seismic residual (a) and correlation coefficient (b) evaluation with iterations (the blue line is the joint inversion and the red line represents the two-step inversion).

Real data application

The Sangtamu region, which is located in the central Lunnan Low Uplift of the Tabei Uplift in Tarim Basin, is an important area where petroleum is enriched and produced. During the Triassic, under the influence of the activity of the southern margin, it experienced a multiphase period of uplift and subsidence, and therefore, the sedimentary characteristics are complicated and variable. The natural gas accumulation conditions are complex, and the structure does not exert a primary influence on gas accumulation. The reservoirs are characterized by strong diagenesis, and the effective pores are mainly secondary pores, which have the characteristics of low porosity, low permeability, and spatial heterogeneity, which is directly expressed as a certain range of overlap in the P-wave velocity, S-wave velocity, and density between the reservoir and the surrounding rock. It is difficult to conduct reliable reservoir prediction and description via traditional reservoir prediction methods. The reservoir petrophysical parameters can be directly applied to reservoir fluid detection and lithology identification, which provides a great advantage in solving associated geological problems. Based on logging and seismic data, we can perform a feasibility analysis of the proposed joint inversion technique for reservoir elastic and petrophysical parameters via prestack seismic facies-controlled joint inversion for hydrocarbon “sweet spot” prediction.

Real data are used to validate the application of the proposed inversion method. True-amplitude processing has been implemented before inversion of the real data. The seismic data were processed by a contractor and the processing sequence was defined to ensure that the final prestack amplitudes should image the reflection strength of the subsurface interfaces as correctly as possible. Figure 7 illustrates the seismic section taken from the Sangtamu region with a fluid interpretation of a well inserted into the section based on drilling results. There are two interpretation layers located at times of 3210 and 3240 ms within the inserted well, wherein red indicates a gas layer and blue indicates a water layer. The water and gas layers all clearly show a positive amplitude; thus, from the amplitude interpretation, we cannot distinguish the different types of reservoirs.

Figure 7.

Seismic section originating from the Sangtamu region.

Figure 8 displays the inverted section of the elastic parameters. The resolution of the inverted result is relatively high, and as such, the reservoir and surrounding rock can be distinguished clearly; however, the elastic parameter properties of the water-bearing layer and the gas-bearing layer are similar, and it is thus difficult to differentiate water and gas directly using elastic parameters. The inverted results show that the variation range for the P-wave velocity of the gas-bearing layer is 4200–4300 m/s, the range of the S-wave velocity is 2100–2300 m/s, and the density range is 2.40–2.55 g/cm⁻³. The ranges for the same parameters for the water-bearing layer are essentially consistent with those of the gas-bearing layer due to the existence of overlaps in the P- and S-wave velocities and the density between the water- and gas-bearing layers, and thus, only a portion of the reservoir can be identified.

Figure 8.

The inversion results of the elastic parameters. (a) P-wave velocity, (b) S-wave velocity, and (c) density.

The ratio of the P-wave and S-wave velocities is one of the most reliable transformation parameters of elastic parameter inversion. Figure 9 illustrates the ratio of the P-wave and S-wave velocities within the inverted section of the reservoir. As demonstrated, the differences in the velocity ratio for the gas- and water-bearing sands are very small, ranging from 2.0 to 2.1. This means that it is difficult to identify the fluid type through a simple application of the ratio of P-wave and S-wave velocities, and in some cases, this method does not even work. During the calculation of the lithological physical parameters using the elastic parameters, the error in the inversion of the elastic parameters will further lead to an error in the petrophysical estimation results.

Figure 9.

The ratio of P-wave and S-wave velocities in the inversion section.

Figure 10 demonstrates the inverted section of the petrophysical parameters. As illustrated, the porosities of the gas- and water-bearing layers are similar, and the porosity is obviously higher than that of the surrounding rock, which is mainly composed of mudstone. However, there are obvious differences in the water saturation of the layers, wherein the gas-bearing layer has a lower value of water saturation relative to that of the water-bearing layer. This is consistent with the basic characteristics of gas and water. From the inversion results for the clay content, we can observe that both the gas- and water-bearing layer possess relatively low clay contents, which is consistent with the drilling data. The significance of the porosity, water saturation, and clay content in the interpretation of the inversion results is readily apparent.

Figure 10.

The inversion section illustrating the petrophysical parameters. (a) Porosity, (b) water saturation, and (c) clay content.

For intuitive delineation about capability and applicability of the proposed method, water- and gas-bearing reservoirs in Figure 9 are chosen in small-scale 3D experimental inversion at target layers of 10 ms top and bottom extended, respectively. Figure 11 displays the horizon slice of water saturation of the gas- and water-bearing layer, within which the dashed box indicates the distribution area of the favorable reservoir based on sedimentary facies division and hydrocarbon accumulation analysis. We can clearly find that the gas-bearing layer has a lower value of water saturation compared with the water-bearing layer at well JF 128, which is the well located in the aforementioned inverted sections. Figure 12 shows the horizon slice of the velocity ratio of the P- and S-waves of the gas- and water-bearing layers, from which it is easily observed that both layers have similar trends, and it is difficult to distinguish the fluid properties through the velocity ratio attribute (i.e. there exist multiple solutions and limitations). We therefore combine the inverted results of elastic parameters with the inverted results of petrophysical parameters, which can effectively distinguish the properties of different fluids.

Figure 11.

The horizon slice of water saturation inversion results. (a) Gas layer and (b) water layer.

Figure 12.

The horizon slice of the velocity ratio of the P- and S-waves. (a) Gas layer and (b) water layer.

Discussion

To solve the problem of cumulative error in the stepwise inversion of reservoir physical parameters, this paper proposes a joint method employing physical and elastic parameters to eliminate the cumulative error. As can be observed from the theoretical data, the resolution of the joint inversion is obviously higher than that of two-step inversion. In the practical data application, the accuracy of the elastic parameter inversion is guaranteed. However, the accuracy of the inversion of the physical parameters is low and necessitates further improvement. The main reasons for this are as follows: First, the physical conditions of the petrological equation are harsh, but the subsurface media are complex and cannot strictly satisfy the assumptions of the physical equation of the lithologies. Second, the limitations of the observed seismic data and the interference of noise cause the seismic signal to contain harmful signals independent of the reservoir characteristics. It is difficult to obtain the elastic parameters such as the P- and S-wave velocities and the density, and it is difficult to guarantee the accuracy of the reservoir physical parameters since seismic data are less sensitive to reservoir physical parameters. Third, the deterministic inversion can only utilize the information within a limited seismic bandwidth, and the deterministic inversion result is the optimal solution with a mean property, which theoretically determines that the resolution is not high. Fourth, the derivative of the petrophysical equation uses a differential form, which loses a certain degree of precision.

Various widely applied reservoir lateral prediction technologies greatly simplify the complexity of the petrophysical problems, which is suitable for reservoirs with simple structures and relatively uniform lateral changes. Wyllie’s time average equation (Wyllie et al., 2002) is broadly employed for its simplicity and can only describe the relationship between velocity and porosity in a medium porosity range. The Raymer equation (Nur et al., 1995) is a modified time averaging equation that applies to different porosity ranges, although both models are suitable for pure sandstone but not for argillaceous sandstone. Compared with other petrophysical models, the Gassmann equation can better establish the relationship between the effective elastic parameters of the fluid-saturated rock and the pore fluid, and it is helpful for understanding the effect of reservoir petrophysical parameters on the elastic properties of the reservoir.

The single sparse constraint may make the inversion unstable, so other constraints should be added for stability. Seismic facies constraints use the low frequency of the geological layer and provide soft constraints, minimizing the L₂ norm of the low frequency difference between the inversion results and the model and improving the robustness of the inversion compared with smooth constraints. However, although the smoothness constraint can stabilize the low frequency component of the inversion result, an improper choice of the constraint weight may lead to ambiguity of the inversion results, whereas the proposed approach can obtain stable low frequency constraints on the inversion result. Since the seismic facies constraints are relatively insensitive to constrained weight coefficients, this technique can obtain an accurate low frequency trend.

It should be noted that the petrophysical properties of the reservoir obtained from inversion have the ability to characterize the gas reservoir, but there are some multiple solutions, we can utilize log curves to build an a priori information model to constrain the inversion and thereby make the inversion results more in-line with the actual situation. In future studies, we can try to improve the applicability of the method by optimizing the algorithm on the basis of a no-differential operation, following which the resolution of the physical parameters inversion will be improved. At the same time, we can attempt to use a parallel algorithm to accelerate the speed of the stochastic inversion and break through the limitations of the seismic frequency bandwidth and the large magnitude of the stochastic inversion.

Conclusions

We have proposed a new prestack seismic simultaneous deterministic inversion of reservoir elastic and petrophysical parameters based on seismic facies constraints. Through the construction of a matrix expression for low-pass filtering, seismic sequence division is applied to the seismic inversion, and the L₂ norm of the difference between the low frequency components of the inversion results and the well logs is minimized. The inversion result is relatively insensitive to the change of the weight coefficient of the low frequency constraint, which reduces the risk of inaccuracy of the low frequency information caused by artificial factors and more effectively guarantees the accuracy of the inversion of the low frequency components. We use the Gassmann equation to construct the relationship between the prestack seismic data and petrophysical parameters, and the elastic and petrophysical parameters are inverted from the prestack seismic data using deterministic optimization methods. Compared with the stochastic optimization method, the proposed method utilizes gradient information of the objective function, wherein the optimal solution can be obtained more quickly in favor of a large-scale production mission. Compared with Bosch’s poststack deterministic optimization, prestack seismic amplitudes have more information that vary with the offset. Finally, compared with Wyllie’s time average equation, the Gassmann equation is better for establishing the relationship between the reservoir elastic parameters and petrophysical parameters, which helps with understanding the impact of reservoir parameters on the elastic properties and reducing the uncertainty of reservoir prediction and fluid identification.