Abstract

The seismic attenuation should be considered while accounting for the effect of anisotropy on the seismic wave propagating through a saturated fractured porous medium. Based on the modified linear-slip theory and anisotropic Gassmann’s equation, we derive an analytical expression for a linearized PP-wave reflection coefficient and an azimuthal attenuation elastic impedance (AAEI) equation in terms of fluid/porosity term, shear modulus, density, dry normal and tangential fracture weaknesses, and compressional (P-wave) and shear (S-wave) attenuation parameters in a weak-attenuation isotropic background rock containing one single set of vertical aligned fractures. We then propose an AAEI inversion method to characterize the characteristics of fluids and fractures using two kinds of constrained regularizations in such a fractured porous medium. The proposed approach is finally confirmed by both the synthetic and real data sets acquired over a saturated fractured porous reservoir.

1. Introduction

With the increasing demand for oil and gas around the world, fractured reservoirs have become the focus of the geophysical exploration of hydrocarbons [1]. When the fractures develop to a certain scale, the seismic wave propagation in such fractured reservoirs may result in shear-wave splitting, azimuthal velocity anisotropy, etc. [2, 3]. A single set of vertical aligned fractures along the dominant horizontal direction in a homogeneous isotropic background medium can be regarded as a long-wavelength equivalent horizontal transversely isotropic (HTI) medium [4, 5]. Therefore, the detection and characterization of the underground fractures using seismic data is crucial for the exploration and exploitation of fractured hydrocarbon reservoirs.

The equivalent medium theory (EMT) is used to equalize an inhomogeneous fractured medium to a homogenous anisotropic medium under the long-wavelength assumption, creating a model which connects the seismic reflection characteristics with the fracture parameters. Hudson’s [6, 7] isolated penny-shaped fracture model is the most commonly used theoretical model of a fracture equivalent medium, in which the equivalent elastic property of isolated fractured rock is derived based on the scattering theory analysis of the average seismic wave field of thin penny-shaped elliptical fractures. It is the first model used to describe the effect of fracture density and orientation on the anisotropic characteristics of sparsely distributed fractured media. Schoenberg’s [8, 9] linear-slip theory is another theoretical model of fracture equivalent medium, and it ignores the shape and microstructure of the fractures, in which the fractures are regarded as a thin layer surrounded by a two-dimensional infinite plane. The fractured rock is modelled as a fractured porous medium filled with the fluid which can move from fracture to interconnected background pores and vice versa, during the wave propagation. The elastic properties of such a medium are affected by the characteristics of fluid so they vary between two extremes which are the gas-saturated (dry) and the fluid-saturated (wet) porous fractured media [4].

The propagation of seismic waves in fractured media not only produces the fluid flow between the fully saturated or partially saturated fractures but also compresses the matrix pores to generate a pressure gradient and results in the fluid flow between the fractures and the matrix pores [10]. Thomsen’s [11] equant-porosity model of fractured porous media firstly considers the effect of fluid flow between spherical equant pores and aligned fractures. Under the low-frequency conditions, the fluid pressure of fractured media strikes a new balance in the half cycle of the seismic wave when there is fluid exchange between fractures and matrix pores. At the moment, the characteristics of seismic wave propagation are consistent with the anisotropic Gassmann’s [12] equation. However, Thomsen’s model is based on a special, idealized fracture geometry (so-called penny-shaped fractures), which is limited to small fracture density (i.e., sparse distribution) conditions. To overcome the penny-shaped fracture model assumption, Gurevich [13] combined the anisotropic Gassmann’s [12] equation with the more generally linear-slip (LS) theory to propose a saturated fractured porous model, who derived the exact analytical expressions for the stiffness matrix of a fluid-saturated fractured porous rock, related to the dry elastic properties of isotropic background, matrix moduli, porosity, dry fracture compliances, and saturated fluid modulus. On this basis, Sil et al. [14] and Huang et al. [15] analyzed the effects of fluid substitution on elastic properties and reflection coefficients in saturated fractured porous media with HTI symmetry and orthorhombic symmetry, respectively. Using the anisotropic Gassmann equation and the linear-slip theory model, Pan and Zhang [16] derived the weakly anisotropic approximations of fluid substitutions and reflection coefficients for a set of parallel, aligned vertical fractures (i.e., an equivalent HTI medium) and two sets of orthogonal vertical fractures (i.e., an equivalent orthorhombic medium), respectively.

The above-mentioned fractured models are independent of frequency. However, seismic waves undergo strong velocity dispersion and attenuation when propagating in a fractured porous medium; thus, their elastic responses are frequency-dependent [10, 14]. Based on the traditional Hudson’s model, Hudson [17] took the effect of fluid flow in fractured media into consideration and analyzed the anisotropy and attenuation characteristics of seismic wave propagation caused by fluid flow in partially saturated fractures. Hudson et al. [18] considered the effects of fluid flow not only between parallel-connected fractures but also between fractures and matrix pores. Using the fluid diffusion equation, Hudson et al. [19] ignored the interaction between fractures and derived the equivalent stiffness matrix based on the effect of fluid flow between the fractures and the matrix pores, which was consistent with Thomsen’s [11] fractured porous model under the high-frequency condition, but inconsistent with the anisotropic Gassmann’s [12] equation under low-frequency conditions. Combining Thomsen’s fractured porous model with the poroelastic BISQ (Biot-Squirt) model, Parra [20] considered the local squirt flow between matrix pores and fractures and studied the dispersion and attenuation characteristics of seismic wave velocity in saturated fractured porous rocks. However, it mainly suited for the sonic and ultrasonic bands and seriously underestimated the velocity dispersion and attenuation in seismic bands. Using the squirt flow mechanism, Chapman et al. [21] derived the frequency-dependent equivalent complex stiffness tensor based on the effects of fluid flow between matrix pores, randomly arranged cracks, and aligned fractures. Jacobsen et al. [22] studied the acoustic characteristics of fluid flow in complex porous media based on the T-matrix method, which can better predicted reservoir parameters such as permeability. Ba [23] derived the wave propagation equation in dual-porosity media, which focused on the mathematical and physical significance of dual-porosity wave theory. Based on Biot’s theory of porous media, Tang [24] introduced the effect of squirt flow in fractures on the elastic media, and proposed the unified theory of elastic wave in fractured porous media. Brajanovski et al. [25] analyzed the variation of P-wave attenuation with normal fracture weakness and the variation of inverse quality factor with frequency in saturated fractured rocks. Chichinina et al. [26, 27] proposed a modified Schoenberg’s linear-slip model to introduce the attenuation anisotropy of fractured media and derived the P- and S-wave attenuation, which was consistent with the attenuation anisotropy factor of plane wave propagation derived by Zhu and Tsvankin [28]. Based on the modified linear-slip model proposed by Chichinina et al. [26], Chen and Innanen [29] considered the seismic attenuation of background media and further improved the modified linear-slip theory. But they did not consider the influence of matrix pores and fluid flow. Based on the traditional linear-slip model, Pan and Zhang [16] considered the effects of matrix porosity and fluid flow but ignored the influence of background attenuation and fracture-induced attenuation.

In this paper, we consider the influence of matrix porosity and fluid flow and integrate the background attenuation and the fracture-induced attenuation to study the characteristics of attenuation anisotropy in fractured porous media. Using the modified Schoenberg’s linear-slip model proposed by Chichinina et al. [26, 27], we first derive the stiffness tensor of a fracture-induced equivalent HTI attenuation medium, which relates to the elastic properties of isotropic background, porosity, matrix modulus, fluid modulus, fracture weaknesses, and attenuation parameters. A linearized PP-wave reflection coefficient and an azimuthal attenuation elastic impedance (AAEI) equation in terms of fluid/porosity term, shear modulus, density, fracture weaknesses, and attenuation parameters are then derived, and an iteratively AAEI inversion approach with regularizations are proposed to characterize the fluids and fractures. In the following sections, we present the methodology, the synthetic tests, and applications of the proposed AAEI inversion approach on a real data set.

2. Theory and Methods

2.1. Derivation of Stiffness Tensor for Fractured Porous Media Based on Modified Linear-Slip Model

To derive the weak-anisotropy linearized approximation of PP-wave reflection coefficient in a fractured porous medium, we first derive the stiffness tensor based on the modified linear-slip (LS) model proposed by Chichinina et al. [26, 27].

Chichinina et al. [26, 27] assumed that the background attenuation is much smaller than the fracture-induced anisotropic attenuation at the LS interfaces, so they only considered the fracture-induced attenuation to characterize the stiffness matrix of fractured rocks and ignored the attenuation induced by the background rocks to simplify the research process of attenuation anisotropy. In this paper, we comprehensively consider the attenuation of background rocks and the fracture-induced anisotropic attenuation , expressing the background isotropic moduli and fracture parameters as the complex forms, i.e., where and and represent the complex P-wave modulus, the first and the second Lamé constants of background isotropic attenuation rocks, respectively, and ; () and () represent the complex normal and tangential fracture weaknesses, respectively, which can be expressed using the fracture-induced P- and S-wave attenuation and as [27]

Pointer et al. [30] demonstrated that the complex tangential fracture weakness parameter does not change with frequency, which is a real number in the range of seismic frequencies, i.e., .

Under the assumption of an isotropic viscoelastic background rock, the complex model parameters can be written as [31, 32] where and represent the background P- and S-wave inverse quality factors.

To describe the effect of saturated fluids on the viscoelastic properties of fractured porous media and consider the effect of fluid flow between the matrix pores and fractures, we express the anisotropic Gassmann’s fluid substitution equation as a complex form, i.e., where is the complex stiffness matrix of saturated fractured porous rocks, is the complex stiffness matrix of dry fractured porous rocks, and can be written as where represents the effective bulk modulus of rock solid particles. In equation (4), is a direct analog of complex Gassmann’s pore space modulus, which is written as where represents an anisotropic analog of complex bulk modulus of dry rocks, is the porosity of fractured rocks, and is the effective bulk modulus of pore fluids.

Chen et al. [33] demonstrated that the imaginary part of the complex reflection coefficient is much smaller than the real part. Therefore, we only emphasize the derivation of the real part of the complex reflection coefficient and ignore the imaginary part. The real part of complex stiffness equation (1) can be thus expressed as where and in which and represent the comprehensive attenuation factors of the background P- and S-wave attenuation and fracture-induced P-wave attenuation.

Substituting equations (5), (6), and (7) into equation (4), and based on the weak anisotropy and weak background attenuation, we express the saturated stiffness matrix of fractured porous rocks in the seismic range as where and in which is the Biot coefficient.

Equation (9) provides the description of a weak-anisotropy approximation for elastic properties of fractured porous rocks.

2.2. Linearized PP-Wave Reflection Coefficient and Elastic Impedance of Fracture-Induced Equivalent HTI Attenuation Media

Based on the scattering theory [34], the PP-wave reflection coefficient of an HTI medium can be expressed as where is the angle of incidence, is the density term of homogenous isotropic background, represents the perturbations in stiffness matrix of saturated fractured porous rocks, and is related to the slowness vector and polarization vector [35].

Therefore, we can combine equations (9) and (11) to derive the PP-wave reflection coefficient of fracture-induced equivalent HTI attenuation media, which is given by where

In equation (12), is the angle of azimuth, and the fracture weaknesses can be estimated using the well log data and rock-physics model (see [36, 37]). Following Connolly [38] and Martins [39], the relationship between the PP-wave reflection coefficient and the azimuthal attenuation elastic impedance (AAEI) can be written as where represents the azimuthal attenuation elastic impedance.

Under the assumption of weak contrasts of elastic parameters (i.e., , , and ), small fracture weaknesses (i.e., and ), and weak attenuation (i.e., and ), the relative contrasts of background elastic moduli and azimuthal attenuation EIs in equation (14) can be substituted as , , , and . Under the assumption of continuity variation in elastic parameters, fracture weaknesses, attenuation parameters, and azimuthal attenuation EIs, they can be further written as , , , , , , , and . Therefore, equation (14) can be written as

By integrating equation (15) and performing corresponding operations, an AAEI equation in the logarithmic domain can be obtained as

For the logarithm operation of equation (16), the final AAEI equation can be expressed as

2.3. Azimuthal Attenuation Elastic Impedance (AAEI) Inversion for Fluid and Fracture Characterization

The AAEI equation in the logarithmic domain can be expressed as the sum of the background isotropic EI and the fracture-induced azimuthal attenuation EI. Based on the logarithmic EI difference between different azimuths, we can realize the inversion for fracture weaknesses and fracture-induced attenuation parameters, while the background elastic moduli can be separately inverted by using the azimuth seismic data with the same observation azimuth as the fracture orientation. Therefore, we can perform the AAEI inversion of fluid parameters, fracture parameters and attenuation parameters for three steps: (1) the inversion of logarithmic attenuation EIs in different azimuths, (2) the inversion of fluid parameters and background elastic moduli using single seismic data with fracture azimuth, and (3) the inversion of fracture parameters and attenuation parameters based on the difference between different azimuthal EIs.

Firstly, we invert the logarithmic attenuation EIs using seismic data with different azimuths. According to equation (14), the linear relationship between logarithmic EIs and azimuthal data can be written as where is the number of incident angles, is the number of azimuthal angles, and is the number of reflection interface. represents the seismic data with () azimuth angles and incident angles, is the difference matrix, is the wavelet matrix, and is the azimuthal attenuation EIs in logarithmic domain, which are given by

Therefore, the logarithmic attenuation EIs can be estimated using the model-based least-squares inversion method [33, 36].

Secondly, we estimate the background fluid/porosity term, the shear modulus, and the density. The least-squares ellipse fitting (LSEF) method is used to estimate the approximate fracture orientation [36]. Then, we use the seismic data with single azimuth of fracture orientation to estimate the background elastic moduli. This is a classic inverse problem of prestack seismic inversion for three parameters and can be resolved (and not explained here) using the inversion method proposed by Downton (2005) and Russell et al. [40].

Using matrix parameterization, the linear expression for the azimuthal differences in the logarithmic attenuation EI can be expressed as a matrix form where is the number of azimuthal difference angles, represents the differences of AAEI in the logarithmic domain, represents the estimated fracture weaknesses and fracture-induced attenuation parameters, and represents the forwarding matrix related to the weighing coefficient matrix in different azimuths, which are given by

Considering the decorrelation between model parameters during the AAEI inversion [36], the decorrelated kernel matrix becomes and the model parameter vector becomes . Therefore, equation (20) can be expressed as

Using the Cauchy probability distribution as a prior probability density function (PDF), and the Gaussian distribution as the likelihood function, the posterior PDF can be solved using the joint PDF of the prior PDF and the likelihood function, i.e., where is the noise variances of seismic data and is the variances of model parameters. Maximizing the posterior equation (23) and combining the initial-model constrained low-frequency regularization term, the objective function can be expressed as where represents the regularization coefficients of model parameters, , and , in which represents the initial values of model parameters. Solving equation (24) leads to where and .

The iteratively reweighted least-squares (IRLS) optimization algorithm [36, 41] is used to solve equation (25), and we can finally get the model parameters .

3. Examples

3.1. Synthetic Examples

The reliability of the AAEI inversion method proposed in this paper is verified by using a synthetic azimuthal prestack seismic data set. Based on the convolution model and a 35 Hz Ricker wavelet, four sets of azimuthal gathers in angle domain are synthesized, as shown in Figure 1(a), where the well log data of fluid modulus, shear modulus, density, fracture weaknesses, and attenuation parameters are displayed in Figure 2(a). It should be noted that the well log data of fluid modulus, shear modulus, and density can be calculated using conventional P- and S-wave velocity and density logging data [40, 42], and the fracture weakness parameters can be estimated using the rock-physics model [36], while the well log data of attenuation parameters can be calculated using the empirical formulas (Haase and Stewart, 2004). Adding appropriate random noises to the noise-free gathers (shown in Figure 1(a)), we can synthesize the noisy gathers with signal-to-noise ratios (SNRs) of 5 and 2, as shown in Figures 1(b) and 1(c), respectively. Figures 2(a) and 2(b) show the inverted model parameters using the AAEI inversion without noise, while Figures 3 and 4 show the corresponding inversion results of model parameters with SNRs of 5 and 2, respectively. From the inversion results, it can be seen that when the synthetic gather does not contain noises or contains appropriate noises, the inversion results of model parameters agree well with the actual data. But as the noises increase, the inversion result of the attenuation parameters becomes not very stable. Figures 5(a)5(c) shows the comparisons between the azimuthal gathers synthesized by the inverted model parameters and the initial gathers. The comparisons show that the difference between the two is small, which further demonstrates the feasibility of the proposed AAEI inversion approach for the model parameters.

3.2. Field Data Example

The target reservoir in the area is a fractured carbonate gas-bearing reservoir, and its main lithology is thick gray dolomite. The target reservoir has the characteristics of low porosity and low permeability and has the characteristics of multistage accumulation, which belongs to a structural gas reservoir. The seismic data of the target area is relatively high in SNR, and it has good lateral continuity and high vertical resolution. Before the seismic inversion process, the seismic amplitudes need to be preserved, including fine wave-front diffusion compensation, inverse Q filtering, surface consistency processing, prestack denoising, and multiple wave removal. In addition, when the seismic data is divided into the azimuth, it is necessary to ensure that each of the azimuthal seismic data divided has sufficient coverage to ensure the SNR of the seismic data and that the coverage times of each azimuthal seismic data are as uniform as possible. Figures 6(a)6(d) is the through-well (A and B) azimuthal partial angle-stacked gathers, respectively, in which four azimuths are 22.5° (overlay range 0°-45°), 67.5° (overlay range 45°-90°), 112.5° (overlay range 90°-135°), and 157.5° (overlay range 135°-180°), and three incident angles are 5° (overlay range 0°-10°), 15° (overlay range 10°-20°), and 25° (overlay range 20°-30°). As shown in Figure 6(a), the two black lines are the positions of well A and well B and the red ellipse indicates the target gas-bearing reservoir, in which obviously strong amplitude anomalies can be found. It should be pointed out that the fracture parameters estimated by the rock-physics model and the attenuation parameters calculated by empirical relationship in well A are used to construct the initial model parameters, while well B is used to verify the prediction result.

To estimate the background fluid/porosity term, the shear modulus, and the density, the approximate fracture orientation is first estimated. The least-squares ellipse fitting (LSEF) method is used to estimate the properties of fracture orientation and fracture density. Figure 7 shows the estimated fracture orientation and fracture density parameters. Figure 8 shows the estimated fracture normal, and we find that the fracture normal is about 100° (or 280°). Therefore, the background fluid/porosity term, shear modulus, and density parameters can be estimated using the three angle-stacked seismic data in the azimuth shown in Figure 6(c). On this basis, we then perform the AAEI inversion method to estimate the fracture weaknesses and the attenuation parameters.

Figures 9(a)9(d) is the estimated azimuthal attenuation EI profiles of the four azimuths and three angles. It can be seen from these figures that the inverted EI values of the target reservoir reveal a low-value anomaly. Figures 10, 11, and 12 illustrate the estimated fluid/porosity term, shear modulus, and density parameter profiles; normal and tangential fracture weakness parameter profiles; and attenuation parameter profiles. In these figures, the red ellipses at well A denote the target gas-bearing fractured porous reservoir, while the red ellipses at well B denote the high-quality gas producing area that is drilled. From the inversion results, it can be seen that the inverted fluid/porosity term at well A shows a low value, while the inverted fracture parameters and attenuation parameters show high values. Correspondingly, the area of low-value fluid/porosity term, high-value fracture parameters, and high-value attenuation parameters at well B is a high-yield gas target layer, which is consistent with the fracture development position and the drilling results. So it further confirms the AAEI inversion approach proposed in this paper to estimate the fracture parameters and attenuation parameters.

4. Discussion

According the derivation process of Pan et al. [43], it is obvious that the derived linearized PP-wave reflection coefficient with real background elastic coefficients and complex fracture characteristics can be written in terms of a real part and an imaginary part of the complex reflection coefficient. The real part of the PP-wave complex reflection coefficient is the same as the PP-wave reflection coefficient with real background elastic coefficients and real fracture characteristics. However, the imaginary part of the PP-wave complex reflection coefficient with real background elastic coefficients and complex fracture characteristics is related to the real background elastic coefficients and the attenuation parameters due to the fractures. Following the numerical analysis of Chen et al. [33], we found that the real part of the complex reflection coefficient is much more than the imaginary part. The imaginary part of the complex reflection coefficient can be thus neglected if we invert the background elastic properties, fracture parameters, and attenuation parameters using the complex reflection amplitudes. Therefore, the attenuation parameters caused by fractures cannot be reasonably inverted by only considering the real background elastic coefficients.

However, the attenuation parameters due to the presence of fractures can be reasonably inverted using the real part of the complex reflection amplitudes when we consider the complex background elastic coefficients and complex fracture characteristics. This is because the attenuation parameters can affect the real part of the complex reflection amplitudes (see equation (12)). As a result, we consider the complex background elastic coefficients while accounting for the effect of the complex fracture characteristics and derive an analytical expression for a linearized PP-wave reflection coefficient and an azimuthal attenuation elastic impedance (AAEI) equation to estimate the fluid/porosity term, the shear modulus, the density, the dry normal and tangential fracture weaknesses, and the compressional (P-wave) and shear (S-wave) attenuation parameters in a weak-attenuation background media permeated by aligned fractures.

Three effects can cause the attenuation of seismic waves, including the geometric spreading, the scattering attenuation, and the intrinsic attenuation. They can be divided into the elastic (geometric dispersion, scattering attenuation) and inelastic (intrinsic attenuation) processes [44]. In exploration geophysics, the interesting intrinsic attenuation is the energy loss to heat and internal friction during the wave propagation, which can result from the fluid flow between aligned fractures and/or randomly distributed pores [18, 27, 30]. The linear-slip (LS) model proposed by Schoenberg [8, 9] describes the elastic seismic wave propagating in such media with welded interfaces, which can be treated in two ways, including bedding planes, joints, cracks, and fractures, or as thin soft layers. When it refers to the imperfectly bonded interfaces of fractures, the attenuation of elastic waves may result from the asperities at the rough fracture surfaces [45, 46] or undulations at the even contracts [26, 47]. Moreover, for the limit of the long wave length, the LS model of Schoenberg [8, 9] is equivalent to the penny-shaped crack model of Hudson [6, 7], and the attenuation mechanisms are also the same [4, 26]. This paper focuses on the prestack seismic inversion for the attenuation parameters, so the attenuation mechanism has not been discussed much.

5. Conclusions

Based on the modified Schoenberg’s linear-slip theory, we integrate the attenuation of background rock and fracture-induced attenuation in this paper and propose an AAEI inversion approach to characterize the saturated fractures in a fracture-induced equivalent HTI attenuation medium. We first derive an analytical expression for the stiffness tensor of a weak-attenuation isotropic background containing one single set of vertical aligned fractures, related to the elastic properties of dry isotropic background, porosity, matrix modulus, saturated fluid modulus, dry normal and tangential fracture weaknesses, and compressional (P-wave) and shear (S-wave) attenuation parameters. We then derive a linearized PP-wave reflection coefficient and an AAEI equation in terms of fluid/porosity term, shear modulus, density, two fracture weaknesses, and two attenuation parameters. We finally propose an iteratively AAEI inversion method with Cauchy-sparse constrained regularization and low-frequency constrained regularization and implement the fluid identification and fracture characterization based on the estimates of sensitive fluid indicator, fracture weaknesses, and attenuation parameters, which are confirmed by both the synthetic and real data acquired over a saturated fractured porous reservoir with HTI symmetry.

Since the background fluid parameters and the fracture parameters are more suitable for the seismic identification of gas and water, the identification for oil and water is not very sensitive. On the contrary, the estimated attenuation parameters are more conducive to the seismic identification for oil and water. Therefore, the prestack seismic inversion for background fluid parameters, fracture parameters, and attenuation parameters is more sensitive to the identification of fluid types filling in the fractures.

Data Availability

The well and seismic data used to support the findings of this study have not been made available because the data used in this paper has not been approved to be used publicly, but the data used to support the findings of this study are available from the corresponding authors upon request.

Conflicts of Interest

The authors declare that there is no conflict of interest regarding the publication of this paper.

Acknowledgments

We would like to acknowledge the sponsorship of the National Natural Science Foundation of China (41674130, 41874145) for funding this research. H. Chen is thanked for giving us very valuable suggestions.