Fourier coefficients-based stepwise Bayesian inversion for elastic and fracture parameters using azimuthal seismic data

The inversion of elastic and fracture parameters from azimuthal seismic data plays a critical role in charactering naturally fractured reservoirs. We propose a Fourier coefficients-based stepwise Bayesian inversion method to estimate these reservoir parameters with improved accuracy. Utilizing the relationship between the azimuthal PP-wave reflection coefficient and Fourier series expansion, we first derive the Fourier coefficient equation for horizontally transverse isotropic (HTI) media. We then conduct a sensitivity analysis of azimuthal Fourier coefficients with respect to P- and S-wave velocities, density, and fracture weaknesses. The results indicate that the zeroth-order Fourier coefficient exhibits greater sensitivity to P- and S-wave velocities and density compared to normal and tangential fracture weaknesses, whereas the second-order Fourier coefficient is more responsive to fracture weaknesses than the fourth-order Fourier coefficient. Based on these sensitivity observations, we develop a stepwise Bayesian inversion approach that involves (1) computing the cosine and sine components of the Fourier coefficients for the azimuthal seismic data, (2) estimating normal and tangential fracture weaknesses using use the second-order Fourier coefficient within a Bayesian framework, and (3) recovering background P- and S-wave velocities and density using the zeroth-order Fourier coefficient along with the previously estimated fracture weaknesses. Both synthetic and field data applications confirm the robustness and effectiveness of the proposed Fourier coefficients-based stepwise Bayesian inversion method for estimating elastic and fracture parameters in naturally fractured reservoirs.

Introduction

Natural fractures, prevalent geological features within subsurface rock formations, play a crucial role in controlling reservoir permeability and influencing fluid migration pathways during hydrocarbon production (Sayers and Dean, 2001; den Boer and Sayers, 2018; Pan et al., 2021; Pan and Liu, 2024; Cao et al., 2025). In this context, amplitude variation with offset and azimuth (AVOAz) inversion has emerged as a key technique for quantitatively predicting fracture parameters due to its superior vertical resolution. AVOAz inversion establishes a mapping relationship between seismic data acquired at varying azimuths and the anisotropic medium’s response, enabling the extraction of essential fracture parameters such as fracture density, normal weakness, and tangential weakness—critical for effective reservoir characterization (Mallick et al., 1988; Bachrach et al., 2009; Pan et al., 2020; Ma et al., 2022; Liu et al., 2023; Xie et al., 2024; Li et al., 2025).

Fracture development induces seismic anisotropy, which significantly alters wave propagation characteristics (Schoenberg, 1980; Hudson, 1981; Pan et al., 2017). Traditional isotropic reflection coefficient equations fail to accurately model wavefield response in anisotropic media, prompting the development of various anisotropic AVOAz inversion frameworks. Pšenčik and Martins (2001) formulated a linearized AVOAz relationship for complex anisotropic media, though its inherent nonlinearity poses challenges in practical applications. Shaw and Sen (2006) introduced a joint AVOAz inversion approach, simultaneously estimating fracture weakness parameters and fluid indicators, offering new insights into multiparameter coupled inversion. Bachrach (2015) addressed the uncertainties and non-uniqueness issues in orthorhombic AVOAz inversion, laying the foundation for subsequent regularization techniques. Narhari et al. (2015) successfully predicted fracture density in deep carbonate reservoirs using pre-stack orthorhombic anisotropic AVOAz inversion. Pan et al. (2018) proposed a Bayesian AVOAz inversion framework incorporating prior constraints, significantly improving the estimation accuracy of fracture weakness. Luo et al. (2019) demonstrated the applicability of orthorhombic anisotropic parameter inversion in the Longmaxi Formation shale reservoir, though the process remains unstable due to the high number of parameters involved Isotropic elastic parameters contribute more significantly to reflection coefficients than fracture weakness parameters, and strong nonlinear interactions between anisotropic and elastic parameters lead to substantial parameter coupling effects, complicating the inversion process (Downton and Gray, 2006; Bachrach, 2015; Li et al., 2020; Zhang et al., 2020).

To address these challenges, recent research has focused on differential seismic feature extraction. Downton and Russell (2011) introduced an azimuthal P-wave reflection coefficient reconstruction method based on Fourier decomposition, enabling fracture orientation identification by isolating anisotropic gradient features. Downton and Roure (2015) interpreting the azimuthal Fourier coefficients for anisotropic and fracture parameters in HTI and orthorhombic anisotropic media. Barone and Sen (2018) applied a Fourier-based azimuthal amplitude variation model to the joint inversion of fracture density and orientation in the Haynesville shale reservoir. Building on these advancements, this study integrates the Fourier-form P-wave reflection coefficient equation proposed by Downton and Russell (2011) with the HTI medium P-wave reflection coefficient equation developed by Pan et al. (2018). This integration leads to a fracture weakness inversion method based on reflection amplitude Fourier coefficients, improving the accuracy and robustness of parameter estimation by separately inverting isotropic and fracture weakness parameters.

Our findings demonstrate that combining Fourier decomposition with AVOAz inversion enables high-resolution estimation of normal and tangential fracture weaknesses—key parameters for reservoir characterization and resource management. The proposed method enhances understanding of fracture behavior in subsurface formations and provides a practical tool for optimizing exploration and production strategies in fractured reservoirs. Validation with synthetic and field data confirms the method’s effectiveness and feasibility, paving the way for broader applications across diverse geological settings.

Theory and method

Azimuthal Fourier coefficients characterized by fracture weaknesses in HTI media

For a single set of vertical, parallel fractures embedded in an isotropic background rock, the fractured medium exhibits horizontally transverse isotropic (HTI) symmetry, with the horizontal axis of symmetry perpendicular to the fracture orientation (Bakulin et al., 2000). The linear-slip model (Schoenberg and Sayers, 1995) is commonly used to describe the fracture-induced anisotropy, where the fracture properties are characterized by the dimensionless normal and tangential fracture weaknesses:

${\delta }_N=\frac{\left(\lambda +2\mu \right)K_N}{1+\left(\lambda +2\mu \right)K_N}(1)$

${\delta }_T=\frac{\mu K_T}{1+\mu K_T}(2)$

where $\lambda$ and $\mu$ are the first and second Lamé constants of the isotropic background rock, $K_N$ and $K_T$ denote the normal and tangential compliances induced by the fractures, while ${\delta }_N$ and ${\delta }_T$ denote the normal and tangential fracture weaknesses. These fracture weaknesses range between 0 and 1, reflecting the impact of fractures on seismic wave propagation (Schoenberg and Douma, 1988). According to Hudson’s aligned penny-shaped crack model (Hudson, 1981), the normal and tangential fracture weaknesses can be further expressed as:

${\delta }_N=\frac{4{\rho }_f}{3g\left(1-g\right)+3\alpha \left(K^{\prime }+4/3{\mu }^{\prime }\right)/\left(\pi \mu \right)}(3)$

${\delta }_T=\frac{16{\rho }_f}{3\left(3-2g\right)+12\alpha {\mu }^{\prime }/\left(\pi \mu \right)}(4)$

where $g=V_P^2/V_S^2$ represents the ratio of the square of S-to-P wave velocity for the host rocks, $\alpha =b/a\ll 1$ represents the aspect ratio of fractures, in which a and b are the semimajor and semiminor axis of the penny-shaped cracks with the form of oblate spheroid, respectively, and ${\rho }_f=\nu 〈a^3〉$ represents the fracture density, in which $\nu$ and $〈·〉$ are the number of fractures per unit volume and the volume averaging, respectively (Hudson, 1980).

Characterized by the fracture weaknesses, the PP-wave reflection coefficient in an HTI medium can be expressed as (Pan et al., 2018):

$R_{PP}^{HTI}\left(\theta ,\phi \right)=a_{V_P}\left(\theta \right)R_{V_P}+a_{V_S}\left(\theta \right)R_{V_S}+a_{\rho }\left(\theta \right)R_{\rho }+a_{{\delta }_N}\left(\theta ,\phi \right)R_{{\delta }_N}+a_{{\delta }_T}\left(\theta ,\phi \right)R_{{\delta }_T}(5)$

where

$a_{V_P}\left(\theta \right)={\sec }^2\theta (6)$

$a_{V_S}\left(\theta \right)=-8g{\sin }^2\theta (7)$

$a_{\rho }\left(\theta \right)=1-4g{\sin }^2\theta (8)$

$a_{{\delta }_N}\left(\theta ,\phi \right)=-\frac{1}{2}{\sec }^2\theta {\left[1-2g\left({\sin }^2\theta {\sin }^2\phi +{\cos }^2\theta \right)\right]}^2(9)$

$a_{{\delta }_T}\left(\theta ,\phi \right)=2g{\sin }^2\theta {\cos }^2\phi \left(1-{\tan }^2\theta {\sin }^2\phi \right)(10)$

and $R_{V_P}=\Delta V_P/2V_P$, $R_{V_S}=\Delta V_S/2V_S$, $R_{\rho }=\Delta \rho /2\rho$, $R_{{\delta }_N}=\left({\delta }_{N2}-{\delta }_{N1}\right)/2$, and $R_{{\delta }_T}=\left({\delta }_{T2}-{\delta }_{T1}\right)/2$. Here $\theta$ and $\phi ={\phi }_{obs}-{\phi }_f$ denote the angle of incidence and the azimuth angle between the observed azimuth ${\phi }_{obs}$ and the azimuth of fracture symmetry axis ${\phi }_f$, respectively; $R_{{\delta }_N}$ and $R_{{\delta }_T}$ represents the differences of normal and tangential fracture weaknesses between the upper and lower layers, respectively.

Downton and Russell (2011) formulated the azimuthal PP-wave reflection coefficient in terms of a Fourier series expansion so Equation 5 can be rewritten as

$R_{PP}^{HTI}\left(\theta ,\phi \right)=R_0\left(\theta \right)+{\displaystyle \sum _{\xi =1}^{\infty }}{R}_{\xi }\left(\theta \right)\cos \left(\xi \left(\phi -{\phi }_{\xi }\left(\theta \right)\right)\right)(11)$

where $R_{\xi }\left(\theta \right)$ represents Fourier coefficients varying with the angle of incidence, including the magnitude $R_{\xi }\left(\theta \right)$ and the phase ${\phi }_{\xi }\left(\theta \right)$. Therefore, the azimuthal reflectivity at a specified incidence angle $\theta$ can be expressed as a linear superposition of sinusoidal components, and Each nth order component is characterized by three defining parameters: amplitude modulation governed by coefficient $R_{\xi }$, spatial periodicity scaled by integer index $\xi$, and phase retardation controlled by term ${\phi }_{\xi }$. Only the n = 0, 2, and 4 magnitude weighting terms are nonzero (Downton and Roure, 2015), and Equation 11 can be thus simplified to

$R_{PP}^{HTI}\left(\theta ,\phi \right)=R_0\left(\theta \right)+R_2\left(\theta \right)\cos \left(2\left(\phi -{\phi }_2\left(\theta \right)\right)\right)+R_4\left(\theta \right)\cos \left(4\left(\phi -{\phi }_4\left(\theta \right)\right)\right)(12)$

with magnitudes

$R_0\left(\theta \right)=a_{V_P}\left(\theta \right)R_{V_P}+a_{V_S}\left(\theta \right)R_{V_S}+a_{\rho }\left(\theta \right)R_{\rho }+a_{{\delta }_{N0}}\left(\theta \right)R_{{\delta }_N}+a_{{\delta }_{T0}}\left(\theta \right)R_{{\delta }_T}(13)$

$R_2\left(\theta \right)=a_{{\delta }_{N2}}\left(\theta \right)R_{{\delta }_N}+a_{{\delta }_{T2}}\left(\theta \right)R_{{\delta }_T}(14)$

$R_4\left(\theta \right)=a_{{\delta }_{N4}}\left(\theta \right)R_{{\delta }_N}+a_{{\delta }_{T4}}\left(\theta \right)R_{{\delta }_T}(15)$

where

$a_{{\delta }_{N0}}\left(\theta \right)=-\frac{1}{2}{\sec }^2\theta +2g\left[1-g+\frac{1}{2}{\sin }^2\theta +\left(\frac{1}{2}-\frac{3}{8}g\right){\sin }^2\theta {\tan }^2\theta \right](16)$

$a_{{\delta }_{T0}}\left(\theta \right)=g{\sin }^2\theta \left(1-\frac{1}{4}{\tan }^2\theta \right)(17)$

$a_{{\delta }_{N2}}\left(\theta \right)=g{\sin }^2\theta \left[2g-1+\left(g-1\right){\tan }^2\theta \right](18)$

$a_{{\delta }_{T2}}\left(\theta \right)=g{\sin }^2\theta (19)$

$a_{{\delta }_{N4}}\left(\theta \right)=-\frac{1}{4}{g}^2{\sin }^2\theta {\tan }^2\theta (20)$

$a_{{\delta }_{T4}}\left(\theta \right)=\frac{1}{4}g{\sin }^2\theta {\tan }^2\theta (21)$

To calculate the Fourier coefficients in Equation 12, the form of a Fourier series can be written as the sum of trigonometric functions, i.e.,

$\begin{array}{l}{R}_{PP}^{HTI}\left(\theta ,\phi \right)={\displaystyle \sum _{\xi =0}^{\infty }}\left[{\alpha }_{\xi }\left(\theta \right)\cos \xi \phi +{\beta }_{\xi }\left(\theta \right)\sin \xi \phi \right]\\ ={\alpha }_0\left(\theta \right)+{\alpha }_2\left(\theta \right)\cos 2\phi +{\beta }_2\left(\theta \right)\sin 2\phi \\ +\hspace{0.17em}{\alpha }_4\left(\theta \right)\cos 4\phi +{\beta }_4\left(\theta \right)\sin 4\phi \end{array}(22)$

with coefficients

${\alpha }_0\left(\theta \right)=a_{V_P}\left(\theta \right)R_{V_P}+a_{V_S}\left(\theta \right)R_{V_S}+a_{\rho }\left(\theta \right)R_{\rho }+a_{{\delta }_{N0}}\left(\theta \right)R_{{\delta }_N}+a_{{\delta }_{T0}}\left(\theta \right)R_{{\delta }_T}(23)$

${\alpha }_2\left(\theta \right)=a_{{\delta }_{N2}}\left(\theta \right)\cos 2{\phi }_f{R}_{{\delta }_N}+a_{{\delta }_{T2}}\left(\theta \right)\cos 2{\phi }_f{R}_{{\delta }_T}(24)$

${\beta }_2\left(\theta \right)=a_{{\delta }_{N2}}\left(\theta \right)\sin 2{\phi }_f{R}_{{\delta }_N}+a_{{\delta }_{T2}}\left(\theta \right)\sin 2{\phi }_f{R}_{{\delta }_T}(25)$

${\alpha }_4\left(\theta \right)=a_{{\delta }_{N4}}\left(\theta \right)\cos 4{\phi }_f{R}_{{\delta }_N}+a_{{\delta }_{T4}}\left(\theta \right)\cos 4{\phi }_f{R}_{{\delta }_T}(26)$

${\beta }_4\left(\theta \right)=a_{{\delta }_{N4}}\left(\theta \right)\sin 4{\phi }_f{R}_{{\delta }_N}+a_{{\delta }_{T4}}\left(\theta \right)\sin 4{\phi }_f{R}_{{\delta }_T}(27)$

and the coefficients ${\alpha }_{\xi }\left(\theta \right)$ and ${\beta }_{\xi }\left(\theta \right)$ in front of the cosine and sine functions with K-samples azimuthal seismic data are defined as follows:

${\alpha }_{\xi }\left(\theta \right)=\frac{1}{\pi }{\displaystyle \sum _{k=1}^K}{R}_k\left(\theta ,\phi \right)\cos \xi \phi d\phi (28)$

${\beta }_{\xi }\left(\theta \right)=\frac{1}{\pi }{\displaystyle \sum _{k=1}^K}{R}_{\xi }\left(\theta ,\phi \right)\sin \xi \phi d\phi (29)$

The above parameters can be transformed to magnitude and phase as

$R_{\xi }\left(\theta \right)=\sqrt{{\alpha }_{\xi }^2\left(\theta \right)+{\beta }_{\xi }^2\left(\theta \right)}(30)$

${\phi }_{\xi }\left(\theta \right)=\frac{1}{\xi }\arctan \left(\frac{{\beta }_{\xi }\left(\theta \right)}{{\alpha }_{\xi }\left(\theta \right)}\right)(31)$

We then analyze the sensitivity of azimuthal Fourier coefficients to variations in P- and S-wave velocities, density, and fracture weaknesses. Figure 1 illustrates the effect of P- and S-wave velocities, density, and fracture weaknesses on the zeroth-order Fourier coefficient, while Figures 2, 3 depict the effects of normal and tangential fracture weaknesses on the second- and fourth-order Fourier coefficients, respectively. The sensitivity analysis is conducted by varying five key parameters $R_{V_P}$, $R_{V_S}$, $R_{\rho }$, $R_{{\delta }_N}$, and $R_{{\delta }_T}$ within the range of −0.3 to 0.3 in increments of 0.1. Results demonstrate that the zeroth-order Fourier coefficient exhibits higher sensitivity to P- and S-wave velocities and density than to normal and tangential fracture weaknesses. Conversely, the second-order Fourier coefficient demonstrates greater sensitivity to fracture weaknesses compared to the fourth-order coefficient. These findings suggest that the zeroth-order Fourier coefficient is primarily useful for estimating P-wave and S-wave velocities and density, while both the second- and fourth-order Fourier coefficients provide valuable information for characterizing normal and tangential fracture weaknesses. Of course, the second-order Fourier coefficient is more effective than the fourth-order coefficient for inverting fracture weaknesses.

Figure 1

Figure 1. Sensitivity analysis of zeroth-order Fourier coefficients, where (a) only changes $R_{V_P}$, (b) only changes $R_{V_S}$, (c) only changes $R_{\rho }$, (d) only changes $R_{{\delta }_N}$ and (e) only changes $R_{{\delta }_T}$.

Figure 2

Figure 2. Sensitivity analysis of second-order Fourier coefficients, where (a) only changes $R_{{\delta }_N}$ and (b) only changes $R_{{\delta }_T}$.

Figure 3

Figure 3. Sensitivity analysis of fourth-order Fourier coefficients, where (a) only changes $R_{{\delta }_N}$ and (b) only changes $R_{{\delta }_T}$.

Stepwise Bayesian inversion for elastic and fracture parameters using azimuthal seismic data

From the sensitivity analysis of azimuthal Fourier coefficients, we propose a stepwise Bayesian inversion for elastic and fracture parameters using the azimuthal seismic data.

Firstly, we calculate the cosine and sine components of the Fourier coefficients for the azimuthal seismic data. The azimuthal seismic data vector $\mathbf{d}\left(\theta ,\phi \right)$ can be obtained via the convolution model (Pan and Zhang, 2019) and the estimated seismic wavelets as

$\mathbf{d}\left(\theta ,\phi \right)={\mathbf{R}}_{PP}^{HTI}\left(\theta ,\phi \right)·\mathbf{W}\left(\theta ,\phi \right)(32)$

where $\mathbf{W}\left(\theta ,\phi \right)$ represents the azimuthal seismic wavelet matrix. For seismic data with five azimuths and three angles of incidence, the cosine and sine components of the Fourier coefficients of each order including the influence of the wavelet can be obtained according to the discrete Fourier transform:

$A_{\xi }\left(\theta \right)=\frac{1}{\pi }{\displaystyle \sum _{k=1}^5}{d}_k\left(\theta ,\phi \right)\cos \xi \phi d\phi ,\xi =0,2,4(33)$

$B_{\xi }\left(\theta \right)=\frac{1}{\pi }{\displaystyle \sum _{k=1}^5}{d}_k\left(\theta ,\phi \right)\sin \xi \phi d\phi ,\xi =0,2,4(34)$

where $A_{\xi }\left(\theta \right)$ and $B_{\xi }\left(\theta \right)$ represent the cosine and sine components of the Fourier coefficients for the kth-azimuth seismic data $d_k\left(\theta ,\phi \right)$.

Next, we use the second-order Fourier coefficient to estimate the normal and tangential fracture weaknesses in a Bayesian framework.

Due to the sensitivity of second-order Fourier coefficient to the fracture weaknesses characterized by Equation 14, the coefficients ${\alpha }_{\xi }\left(\theta \right)$ and ${\beta }_{\xi }\left(\theta \right)$ in front of the cosine and sine functions are also sensitive to the fracture weaknesses expressed by Equations 24 and 25. Incorporating Equations 33 and 34, the forward modeling can be expressed in the matrix form as

${\mathbf{d}}_1={\mathbf{G}}_1{\mathbf{m}}_1(35)$

where

${\mathbf{d}}_1={\left[\begin{array}{ccccc}{\mathbf{A}}_2\left({\theta }_1\right)& {\mathbf{B}}_2\left({\theta }_1\right)& \mathrm{L}& {\mathbf{A}}_2\left({\theta }_3\right)& {\mathbf{B}}_2\left({\theta }_3\right)\end{array}\right]}^T(36)$

${\mathbf{G}}_1={\left[\begin{array}{l}\begin{array}{cccc}\mathbf{W}\left({\theta }_1\right){\mathbf{a}}_{{\delta }_{N2}}\left({\theta }_1\right)\cos 2{\phi }_f& \mathbf{W}\left({\theta }_1\right){\mathbf{a}}_{{\delta }_{N2}}\left({\theta }_1\right)\sin 2{\phi }_f& \mathrm{L}& \mathbf{W}\left({\theta }_3\right){\mathbf{a}}_{{\delta }_{N2}}\left({\theta }_3\right)\cos 2{\phi }_f\end{array}\\ \begin{array}{c}\mathbf{W}\left({\theta }_3\right){\mathbf{a}}_{{\delta }_{N2}}\left({\theta }_3\right)\sin 2{\phi }_f\end{array}\\ \begin{array}{cccc}\mathbf{W}\left({\theta }_1\right){\mathbf{a}}_{{\delta }_{T2}}\left({\theta }_1\right)\cos 2{\phi }_f& \mathbf{W}\left({\theta }_1\right){\mathbf{a}}_{{\delta }_{T2}}\left({\theta }_1\right)\sin 2{\phi }_f& \mathrm{L}& \mathbf{W}\left({\theta }_3\right){\mathbf{a}}_{{\delta }_{T2}}\left({\theta }_3\right)\cos 2{\phi }_f\end{array}\\ \begin{array}{c}\mathbf{W}\left({\theta }_3\right){\mathbf{a}}_{{\delta }_{T2}}\left({\theta }_3\right)\sin 2{\phi }_f\end{array}\end{array}\right]}^T(37)$

${\mathbf{m}}_1={\left[\begin{array}{cc}{\mathbf{R}}_{{\delta }_N}& {\mathbf{R}}_{{\delta }_T}\end{array}\right]}^T(38)$

and

${\mathbf{A}}_2\left({\theta }_j\right)={\left[\begin{array}{ccc}{A}_2\left({\theta }_j,t_1\right)& \mathrm{L}& A_2\left({\theta }_j,t_M\right)\end{array}\right]}^T,j=1,2,3(39)$

${\mathbf{B}}_2\left({\theta }_j\right)={\left[\begin{array}{ccc}{B}_2\left({\theta }_j,t_1\right)& \mathrm{L}& B_2\left({\theta }_j,t_M\right)\end{array}\right]}^T,j=1,2,3(40)$

${\mathbf{a}}_{{\delta }_{N2}}\left({\theta }_j\right)=diag\left[\begin{array}{ccc}{a}_{{\delta }_{N2}}\left({\theta }_j,t_1\right)& \mathrm{L}& a_{{\delta }_{N2}}\left({\theta }_j,t_M\right)\end{array}\right](41)$

${\mathbf{a}}_{{\delta }_{T2}}\left({\theta }_j\right)=diag\left[\begin{array}{ccc}{a}_{{\delta }_{T2}}\left({\theta }_j,t_1\right)& \mathrm{L}& a_{{\delta }_{T2}}\left({\theta }_j,t_M\right)\end{array}\right](42)$

${\mathbf{R}}_{{\delta }_N}={\left[\begin{array}{ccc}{R}_{{\delta }_N}\left(t_1\right)& \mathrm{L}& R_{{\delta }_N}\left(t_M\right)\end{array}\right]}^T(43)$

${\mathbf{R}}_{{\delta }_T}={\left[\begin{array}{ccc}{R}_{{\delta }_T}\left(t_1\right)& \mathrm{L}& R_{{\delta }_T}\left(t_M\right)\end{array}\right]}^T(44)$

in which M represents the number of samples for the azimuthal seismic data, and the symbol T represents the transposition of matrix.

We use the Bayesian theory to estimate the normal and tangential fracture weaknesses, in which the posterior probability distribution function (PDF) of the fracture weaknesses $p\left({\mathbf{m}}_1\left|{\mathbf{d}}_1\right.\right)$ can be written as a joint PDF of the prior PDF $p\left({\mathbf{m}}_1\right)$ and the likelihood function $p\left({\mathbf{d}}_1\left|{\mathbf{m}}_1\right.\right)$:

$p\left({\mathbf{m}}_1\left|{\mathbf{d}}_1\right.\right)\propto p\left({\mathbf{m}}_1\right)·p\left({\mathbf{d}}_1\left|{\mathbf{m}}_1\right.\right)(45)$

where $p\left(·\right)$ represents the PDF. Following Pan et al. (2018), the Cauchy PDF can be used as the prior PDF (Alemie and Sacchi, 2011), and the Gaussian PDF used as the likelihood function, thus the posterior PDF can be also expressed as

$p\left({\mathbf{m}}_1\left|{\mathbf{d}}_1\right.\right)\propto {\displaystyle \prod _{i=1}^{2M}}\left(\frac{1}{1+{\mathbf{m}}_1^2/{\sigma }_{{\mathbf{m}}_1}^2}\right)·\exp \left(-\frac{{‖{\mathbf{G}}_1{\mathbf{m}}_1-{\mathbf{d}}_1‖}^2}{2{\sigma }_{{\mathbf{d}}_1}^2}\right)(46)$

where ${\sigma }_{{\mathbf{d}}_1}^2$ and ${\sigma }_{{\mathbf{m}}_1}^2$ represent the variances of seismic data and model parameters, respectively, and the symbol ${‖·‖}^2$ represents the square of the 2-norm.

Combining the smoothing-model regularization (Pan and Zhang, 2018), the maximum posterior solution can be derived by maximizing the posterior PDF as

$\mathbf{J}\left({\mathbf{m}}_1\right)={‖{\mathbf{G}}_1{\mathbf{m}}_1-{\mathbf{d}}_1‖}^2+2{\sigma }_{{\mathbf{d}}_1}^2{\displaystyle \sum _{i=1}^{2M}}\ln \left(1+{\mathbf{m}}_1^2/{\sigma }_{{\mathbf{m}}_1}^2\right)+{\displaystyle \sum _{i=1}^2}{\tau }_i{‖{\chi }_i-{\mathbf{P}}_i{\mathbf{m}}_{1i}‖}^2(47)$

where ${\tau }_i$ represents the regularization coefficients of normal and tangential fracture weaknesses, ${\mathbf{P}}_i={\displaystyle \underset{t_1}{\overset{t_i}{\int }}}·dt$ represents the integral operator, and ${\chi }_i=\frac{1}{2}\ln \left(\frac{{\mathbf{m}}_{1i}}{{\mathbf{m}}_{1i0}}\right)$, in which ${\mathbf{m}}_{1i0}$ represents the initial smoothing model of fracture weaknesses.

Minimizing the objective function of Equation 47 yields

$\mathbf{\psi }=\mathbf{\Gamma }{\mathbf{m}}_1,(48)$

where

$\mathbf{\Gamma }={\mathbf{G}}_1^T{\mathbf{G}}_1+\frac{2{\sigma }_{{\mathbf{d}}_1}^2}{{\sigma }_{{\mathbf{m}}_1}^2}{\mathbf{C}}_{Cauchy}+{\displaystyle \sum _{i=1}^2}{\tau }_i{\mathbf{P}}_i^T{\mathbf{P}}_i(49)$

$\mathbf{\psi }={\mathbf{G}}_1^T{\mathbf{d}}_1+{\displaystyle \sum _{i=1}^2}{\tau }_i{\mathbf{P}}_i^T{\chi }_i(50)$

where ${\mathbf{C}}_{Cauchy}$ represents the Cauchy sparse matrix (Sacchi and Ulrych, 1995).

We solve Equation 48 by using the iteratively re-weighted least-square (IRLS) algorithm (Sacchi and Ulrych, 1995; Pan et al., 2018) obtain the normal and tangential fracture weaknesses via the second-order Fourier coefficient. We then integrate the reflectivities of normal and tangential fracture weaknesses to recover the fracture parameters based on the smoothing-model data as follows:

${\delta }_N\left(t_i\right)={\delta }_N\left(t_1\right)+{\displaystyle \sum _{t_1}^{t_{i-1}}}{R}_{{\delta }_N}\left(\tau \right)(51)$

${\delta }_T\left(t_i\right)={\delta }_T\left(t_1\right)+{\displaystyle \sum _{t_1}^{t_{i-1}}}{R}_{{\delta }_T}\left(\tau \right)(52)$

We finally estimate the background P- and S-wave velocities and density using the zeroth-order Fourier coefficient and the estimated fracture weaknesses in a Bayesian framework.

Due to the sensitivity of zeroth-order Fourier coefficient to the P- and S-wave velocities and density characterized by Equations 13 and 23, we use Equation 23 to express the forward modeling in the matrix form as

${\mathbf{d}}_2={\mathbf{G}}_2{\mathbf{m}}_2(53)$

where

${\mathbf{d}}_2={\left[\begin{array}{ccc}{\mathbf{A}}_0\left({\theta }_1\right)& \mathrm{L}& {\mathbf{A}}_0\left({\theta }_3\right)\end{array}\right]}^T(54)$

${\mathbf{G}}_2={\left[\begin{array}{l}\begin{array}{ccc}\mathbf{W}\left({\theta }_1\right){\mathbf{a}}_{V_P}\left({\theta }_1\right)& \mathrm{L}& \mathbf{W}\left({\theta }_3\right){\mathbf{a}}_{V_P}\left({\theta }_3\right)\end{array}\\ \begin{array}{ccc}\mathbf{W}\left({\theta }_1\right){\mathbf{a}}_{V_S}\left({\theta }_1\right)& \mathrm{L}& \mathbf{W}\left({\theta }_3\right){\mathbf{a}}_{V_S}\left({\theta }_3\right)\end{array}\\ \begin{array}{ccc}\mathbf{W}\left({\theta }_1\right){\mathbf{a}}_{\rho }\left({\theta }_1\right)& \mathrm{L}& \mathbf{W}\left({\theta }_3\right){\mathbf{a}}_{\rho }\left({\theta }_3\right)\end{array}\\ \begin{array}{ccc}\mathbf{W}\left({\theta }_1\right){\mathbf{a}}_{{\delta }_{N0}}\left({\theta }_1\right)& \mathrm{L}& \mathbf{W}\left({\theta }_3\right){\mathbf{a}}_{{\delta }_{N0}}\left({\theta }_3\right)\end{array}\\ \begin{array}{ccc}\mathbf{W}\left({\theta }_1\right){\mathbf{a}}_{{\delta }_{T0}}\left({\theta }_1\right)& \mathrm{L}& \mathbf{W}\left({\theta }_3\right){\mathbf{a}}_{{\delta }_{T0}}\left({\theta }_3\right)\end{array}\end{array}\right]}^T(55)$

${\mathbf{m}}_2={\left[\begin{array}{ccccc}{\mathbf{R}}_{V_P}& {\mathbf{R}}_{V_S}& {\mathbf{R}}_{\rho }& {\mathbf{R}}_{{\delta }_N}& {\mathbf{R}}_{{\delta }_T}\end{array}\right]}^T(56)$

and

${\mathbf{A}}_0\left({\theta }_j\right)={\left[\begin{array}{ccc}{A}_0\left({\theta }_j,t_1\right)& \mathrm{L}& A_0\left({\theta }_j,t_M\right)\end{array}\right]}^T,j=1,2,3(57)$

${\mathbf{a}}_{V_P}\left({\theta }_j\right)=diag\left[\begin{array}{ccc}{a}_{V_P}\left({\theta }_j,t_1\right)& \mathrm{L}& a_{V_P}\left({\theta }_j,t_M\right)\end{array}\right](58)$

${\mathbf{a}}_{V_S}\left({\theta }_j\right)=diag\left[\begin{array}{ccc}{a}_{V_S}\left({\theta }_j,t_1\right)& \mathrm{L}& a_{V_S}\left({\theta }_j,t_M\right)\end{array}\right](59)$

${\mathbf{a}}_{\rho }\left({\theta }_j\right)=diag\left[\begin{array}{ccc}{a}_{\rho }\left({\theta }_j,t_1\right)& \mathrm{L}& a_{\rho }\left({\theta }_j,t_M\right)\end{array}\right](60)$

${\mathbf{a}}_{{\delta }_{N0}}\left({\theta }_j\right)=diag\left[\begin{array}{ccc}{a}_{{\delta }_{N0}}\left({\theta }_j,t_1\right)& \mathrm{L}& a_{{\delta }_{N0}}\left({\theta }_j,t_M\right)\end{array}\right](61)$

${\mathbf{a}}_{{\delta }_{T0}}\left({\theta }_j\right)=diag\left[\begin{array}{ccc}{a}_{{\delta }_{T0}}\left({\theta }_j,t_1\right)& \mathrm{L}& a_{{\delta }_{T0}}\left({\theta }_j,t_M\right)\end{array}\right](62)$

${\mathbf{R}}_{V_P}={\left[\begin{array}{ccc}{R}_{V_P}\left(t_1\right)& \mathrm{L}& R_{V_P}\left(t_M\right)\end{array}\right]}^T(63)$

${\mathbf{R}}_{V_S}={\left[\begin{array}{ccc}{R}_{V_S}\left(t_1\right)& \mathrm{L}& R_{V_S}\left(t_M\right)\end{array}\right]}^T(64)$

${\mathbf{R}}_{\rho }={\left[\begin{array}{ccc}{R}_{\rho }\left(t_1\right)& \mathrm{L}& R_{\rho }\left(t_M\right)\end{array}\right]}^T(65)$

We then use the Bayesian theory to estimate the P- and S-wave velocities and density using the estimated normal and tangential fracture weaknesses, and details about the Bayesian inversion for P- and S-wave velocities and density are the same as the inversion method of normal and tangential fracture weaknesses mentioned above. The final estimated P- and S-wave velocities and density can be expressed as

$V_P\left(t_i\right)=V_P\left(t_1\right)\exp \left[2{\displaystyle \sum _{t_1}^{t_i}}{R}_{V_P}\left(\tau \right)\right](66)$

$V_S\left(t_i\right)=V_S\left(t_1\right)\exp \left[2{\displaystyle \sum _{t_1}^{t_i}}{R}_{V_S}\left(\tau \right)\right](67)$

$\rho \left(t_i\right)=\rho \left(t_1\right)\exp \left[2{\displaystyle \sum _{t_1}^{t_i}}{R}_{\rho }\left(\tau \right)\right](68)$

Results

Synthetic examples

The feasibility of the proposed method was validated through numerical experiments using well log data from a single borehole, incorporating P-wave velocity, S-wave velocity, density, normal and tangential fracture weakness parameters derived from rock physics modeling. Synthetic seismic records were generated for five azimuths (10°, 50°, 90°,130°, and 170°) using a 35 Hz Ricker wavelet and Equation 22, as shown in Figure 4. To evaluate the method’s robustness in realistic noise conditions, Gaussian noise with signal-to-noise ratios (SNRs) of 10 and 5 was introduced to the seismic records, as demonstrated in Figures 5, 6. The inversion results for noise-free, SNR = 10, and SNR = 5 cases are illustrated in Figures 7–9, respectively. Comparative analysis indicates strong consistency between the inversion results of the noisy and noise-free cases. As summarized in Table 1, the correlation coefficients for P- and S-wave velocity inversion results remain stable across different noise levels. Although the correlation of density and normal/tangential fracture weaknesses slightly decreases under SNR = 10 and SNR = 5 conditions, it still maintains a good agreement with the reference, demonstrating the noise robustness and algorithmic stability of the proposed inversion approach. Figure 10 demonstrates the noise susceptibility of inverted normal weakness (upper panel) and tangential weakness (lower panel) through box-and-whisker plots comparing error distributions at SNR = 10 and SNR = 5. For the normal weakness parameter, median errors increase by 66.7% (1.5→2.5 × 10⁻⁵) with SNR reduction, accompanied by interquartile range expansion from 0.8 × 10⁻⁵ to 1.6 × 10⁻⁵, indicating doubled dispersion. The tangential weakness parameter exhibits more pronounced sensitivity: median errors escalate 150% (2→5 × 10⁻⁵) while interquartile range triples (1→3 × 10⁻⁵), revealing non-linear instability growth. The error growth rates quantify parameter-specific noise amplification, establishing the tangential weakness parameter as the critical constraint for system robustness. This differential response necessitates SNR-dependent weighting strategies in joint inversion frameworks, particularly below the identified reliability thresholds.

Figure 4

Figure 4. Synthetic seismic record without noise, where (a) azimuth angle is 10°, (b) azimuth angle is 50°, (c) azimuth angle is 90°, (d) azimuth angle is 130°, and (e) azimuth angle is 170°.

Figure 5

Figure 5. Synthetic seismic record with SNR of 10, where (a) azimuth angle is 10°, (b) azimuth angle is 50°, (c) azimuth angle is 90°, (d) azimuth angle is 130°, and (e) azimuth angle is 170°.

Figure 6

Figure 6. Synthetic seismic record with SNR of 5, where (a) azimuth angle is 10°, (b) azimuth angle is 50°, (c) azimuth angle is 90°, (d) azimuth angle is 130°, and (e) azimuth angle is 170°.

Figure 7

Figure 7. Inversion results without noise, where (a) shows the inverted P- and S-wave velocities and density, (b) shows the inverted normal and tangential fracture weaknesses.

Figure 8

Figure 8. Inversion results with SNR of 10, where (a) shows the inverted P- and S-wave velocities and density, (b) shows the inverted normal and tangential fracture weaknesses.

Figure 9

Figure 9. Inversion results with SNR of 5, where (a) shows the inverted P- and S-wave velocities and density, (b) shows the inverted normal and tangential fracture weaknesses.

Table 1

Table 1. Correlation coefficients of inversion results under different noise conditions.

Figure 10

Figure 10. Error analysis for inversion results of normal and tangential fracture weaknesses in different SNRs.

Field data example

To further validate the feasibility of the proposed method, we applied it to actual data from a research area in southwestern China. This study area is characterized by deep gas reservoirs, high formation pressures, and fractures primarily associated with fault zones, making it suitable for modeling as an HTI medium. Initially, we conducted denoising and amplitude-preserving processing on the prestack seismic data. The processed data were then stacked by angle, resulting in seismic profiles for azimuth angles of 10°, 50°, 90°, 130°, and 170°, with three incident angles of 18°, 22°, and 26°. This produced a total of 15 azimuth-specific seismic profiles, with the corresponding shear-wave velocity logs included in Figure 11. We then extracted the zeroth- and second-order Fourier coefficients (both cosine and sine components) from the processed seismic data. The second-order Fourier coefficients were used to estimate the normal and tangential fracture weaknesses. Using the zeroth-order Fourier coefficients, we further inverted for the normal and tangential fracture weaknesses, as well as the P- and S-wave velocities and density. The inversion results are shown in Figures 12, 13, where the profiles for P-wave and S-wave velocities, density, and normal and tangential fracture intensities demonstrate good continuity and alignment with the logging curves. This strong correlation between the inversion results and well log data confirms the feasibility and stability of the proposed method in real-world applications.

Figure 11

Figure 11. Seismic records for different azimuth angles, where (a) shows the azimuth angle of 10°, (b) shows the azimuth angle of 50°, (c) shows the azimuth angle of 90°, (d) shows the azimuth angle of 130°, and (e) shows the azimuth angle of 170° (the black curve represents the P-wave velocity log).

Figure 12

Figure 12. Inversion results of parameter profiles, where (a) shows the inverted P-wave velocity profile, (b) shows the inverted S-wave velocity profile, and (c) shows the inverted density profile (the red curves represent the corresponding logging curves).

Figure 13

Figure 13. Inversion results of parameter profiles, where (a) shows the inverted normal fracture weakness profile, and (b) shows the inverted tangential fracture weakness profile (the red curves represent the corresponding logging curves).

Conclusion

This study introduced a Fourier coefficients-based stepwise Bayesian inversion method for estimating elastic and fracture parameters in fractured reservoirs using azimuthal seismic data. By analyzing the sensitivity of Fourier coefficients to key subsurface properties, we established that the zeroth-order Fourier coefficient is predominantly sensitive to P- and S-wave velocities and density, while the second-order Fourier coefficient is more responsive to fracture weaknesses. Building upon this sensitivity analysis, we developed a stepwise Bayesian inversion approach, where fracture weaknesses are first estimated using the second-order Fourier coefficient, followed by the estimation of background elastic parameters using the zeroth-order Fourier coefficient. Synthetic experiments demonstrated the robustness of the proposed method, with inversion results maintaining strong consistency with noise-free reference data, even under SNR of 10 and 5. Furthermore, application to field data from a research area in southwestern China confirmed the method’s practical feasibility. The inversion results for P-wave and S-wave velocities, density, and fracture intensities exhibited excellent continuity and alignment with well log data, highlighting the effectiveness and stability of the approach in fractured reservoir characterization.

Data availability statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author contributions

JS: Project administration, Supervision, Funding acquisition, Writing – review and editing, Investigation. HX: Data curation, Software, Visualization, Formal Analysis, Writing – review and editing. SC: Project administration, Supervision, Funding acquisition, Investigation, Writing – review and editing. XZ: Investigation, Supervision, Funding acquisition, Writing – review and editing, Project administration. XP: Software, Methodology, Visualization, Resources, Writing – original draft, Data curation, Funding acquisition, Conceptualization, Writing – review and editing, Project administration, Supervision.

Funding

The author(s) declare that financial support was received for the research and/or publication of this article. We would like to express our gratitude to the sponsorship of the program of China National Petroleum Corporation (2023ZZ05-05 and 2023ZZ14YJ05), the National Natural Science Foundation of China (42474172), and the Natural Science Foundation of Hunan Province (2025JJ20036) for their fundings in this research. The authors declare that this study received funding from China National Petroleum Corporation. The funder was not involved in the study design, collection, analysis, interpretation of data, the writing of this article, or the decision to submit it for publication.

Acknowledgments

We also thank the reviewers for their constructive suggestions.

Conflict of interest

Authors JS, SC, and XZ were employed by PetroChina.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The author(s) declare that no Generative AI was used in the creation of this manuscript.

Publisher’s note

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

References

Alemie, W., and Sacchi, M. D. (2011). High-resolution three-term AVO inversion by means of a Trivariate Cauchy probability distribution. Geophysics 76, R43–R55. doi:10.1190/1.3554627

CrossRef Full Text | Google Scholar

Bachrach, R. (2015). Uncertainty and nonuniqueness in linearized AVAZ for orthorhombic media. Lead. Edge 34, 1048–1056. doi:10.1190/tle34091048.1

CrossRef Full Text | Google Scholar

Bachrach, R., Sengupta, M., Salama, A., and Miller, P. (2009). Reconstruction of the layer anisotropic elastic parameters and high-resolution fracture characterization from P-wave data: a case study using seismic inversion and Bayesian rock physics parameter estimation. Geophys. Prospect. 57, 253–262. doi:10.1111/j.1365-2478.2008.00768.x

CrossRef Full Text | Google Scholar

Bakulin, A., Grechka, V., and Tsvankin, I. (2000). Estimation of fracture parameters from reflection seismic data—part I: HTI model due to a single fracture set. Geophysics 65, 1788–1802. doi:10.1190/1.1444863

CrossRef Full Text | Google Scholar

Barone, A., and Sen, M. K. (2018). A new Fourier azimuthal amplitude variation fracture characterization method: case study in the Haynesville Shale. Geophysics 83, WA101–WA120. doi:10.1190/geo2017-0030.1

CrossRef Full Text | Google Scholar

Cao, Z., Zhang, S., Xue, Y., Wang, Z., Du, F., Li, Z., et al. (2025). Disaster-causing mechanism of spalling rock burst based on folding catastrophe model in coal mine. Rock Mech. Rock Eng. doi:10.1007/s00603-025-04497-6

CrossRef Full Text | Google Scholar

den Boer, L. D., and Sayers, C. M. (2018). Constructing a discrete fracture network constrained by seismic inversion data. Geophys. Prospect. 66, 124–140. doi:10.1111/1365-2478.12527

CrossRef Full Text | Google Scholar

Downton, J., and Russell, H. (2011). Azimuthal fourier coefficients: a simple method to estimate fracture parameters. CSEG Rec. Houston, TX: SEG.

Google Scholar

Downton, J. E., and Gray, D. (2006). AVAZ parameter uncertainty estimation. Houston, TX: SEG Technical Program Expanded Abstracts.

Google Scholar

Downton, J. E., and Roure, B. (2015). Interpreting azimuthal Fourier coefficients for anisotropic and fracture parameters. Interpretation 3, ST9–ST27. doi:10.1190/int-2014-0235.1

CrossRef Full Text | Google Scholar

Hudson, J. A. (1980). Overall properties of a cracked solid. Math. Proc. Camb. Philos. Soc. 88, 371–384. doi:10.1017/s0305004100057674

CrossRef Full Text | Google Scholar

Hudson, J. A. (1981). Wave speeds and attenuation of elastic waves in material containing cracks. Geophys. J. Int. 64, 133–150. doi:10.1111/j.1365-246x.1981.tb02662.x

CrossRef Full Text | Google Scholar

Li, L., Zhang, G., Wang, Y., and Hu, G. (2025). Bayesian Fourier coefficient inversion for gravity-induced stress parameters in tilted transversely isotropic gas-bearing reservoirs. Geophysics 90, M59–M73. doi:10.1190/geo2024-0355.1

CrossRef Full Text | Google Scholar

Li, L., Zhang, J., Pan, X., and Zhang, G. (2020). Azimuthal elastic impedance-based Fourier coefficient variation with angle inversion for fracture weakness. Petro. Sci. 17, 86–104. doi:10.1007/s12182-019-00405-0

CrossRef Full Text | Google Scholar

Liu, X., Zhou, H., Guo, K., Li, C., Zu, S., and Wu, L. (2023). Quantitative characterization of shale gas reservoir properties based on BiLSTM with attention mechanism. Geosci. Front. 14, 101567. doi:10.1016/j.gsf.2023.101567

CrossRef Full Text | Google Scholar

Luo, T., Feng, X., Guo, Z., Liu, C., and Liu, X. (2019). Seismic AVAZ inversion for orthorhombic shale reservoirs in the Longmaxi area, Sichuan. Appl. Geophys. 16, 185–198. doi:10.1007/s11770-019-0759-0

CrossRef Full Text | Google Scholar

Ma, Z., Yin, X., Li, K., and Tan, Y. (2022). Fourier coefficients variation with angle for fracture detection and fluid discrimination in tilted transversely isotropic media. Surv. Geophys. 43, 775–813. doi:10.1007/s10712-022-09704-5

CrossRef Full Text | Google Scholar

Mallick, S., Craft, K. L., Meister, L. J., and Chambers, R. E. (1998). Determination of the principal directions of azimuthal anisotropy from P-wave seismic data. Geophysics 63, 692–706. doi:10.1190/1.1444369

CrossRef Full Text | Google Scholar

Narhari, S. R., Al-Qadeeri, B., Dashti, J., Sliva, J., Dasgupta, S., Hannan, A., et al. (2015). Application of prestack orthotropic AVAz inversion for fracture characterization of a deep carbonate reservoir in northern Kuwait. Lead. Edge 34, 1488–1493. doi:10.1190/tle34121488.1

CrossRef Full Text | Google Scholar

Pan, X., Du, X., Zhao, X., Ge, Z., Li, L., Zhang, D., et al. (2021). Seismic characterization of decoupled orthorhombic fractures based on observed surface azimuthal amplitude data. IEEE Trans. Geosci. Remote Sens. 60, 1–12. doi:10.1109/tgrs.2021.3105724

CrossRef Full Text | Google Scholar

Pan, X., and Liu, J. (2024). Stress-dependent PP-wave reflection coefficient for fourier-coefficients-based seismic inversion in horizontally stressed vertical transversely isotropic media. Surv. Geophys. 45, 1143–1176. doi:10.1007/s10712-024-09841-z

CrossRef Full Text | Google Scholar

Pan, X., and Zhang, G. (2018). Model parameterization and PP-wave Amplitude versus Angle and Azimuth (AVAZ) direct inversion for fracture quasi-weaknesses in weakly anisotropic elastic media. Surv. Geophys. 39, 937–964. doi:10.1007/s10712-018-9481-3

CrossRef Full Text | Google Scholar

Pan, X., and Zhang, G. (2019). Fracture detection and fluid identification based on anisotropic Gassmann equation and linear-slip model. Geophysics 84 (1), R85–R98. doi:10.1190/geo2018-0255.1

CrossRef Full Text | Google Scholar

Pan, X., Zhang, G., and Cui, Y. (2020). Matrix-fluid-fracture decoupled-based elastic impedance variation with angle and azimuth inversion for fluid modulus and fracture weaknesses. J.Petro. Sci. Eng. 189, 106974. doi:10.1016/j.petrol.2020.106974

CrossRef Full Text | Google Scholar

Pan, X., Zhang, G., and Yin, X. (2017). Azimuthally anisotropic elastic impedance inversion for fluid indicator driven by rock physics. Geophysics 82, C211–C227. doi:10.1190/geo2017-0191.1

CrossRef Full Text | Google Scholar

Pan, X., Zhang, G., and Yin, X. (2018). Azimuthal seismic amplitude variation with offset and azimuth inversion in weakly anisotropic media with orthorhombic symmetry. Surv. Geophys. 39, 99–123. doi:10.1007/s10712-017-9434-2

CrossRef Full Text | Google Scholar

Pšenčik, I., and Martins, J. L. (2001). Properties of weak contrast PP reflection/transmission coefficients for weakly anisotropic elastic media. Stud. Geophys. Geod. 45, 176–199. doi:10.1023/A:1021868328668

CrossRef Full Text | Google Scholar

Sacchi, M. D., and Ulrych, T. J. (1995). High-resolution velocity gathers and offset space reconstruction. Geophysics 60, 1169–1177. doi:10.1190/1.1443845

CrossRef Full Text | Google Scholar

Sayers, C. M., and Dean, S. (2001). Azimuth-dependent AVO in reservoirs containing non-orthogonal fracture sets. Geophys. Prospect. 49, 100–106. doi:10.1046/j.1365-2478.2001.00236.x

CrossRef Full Text | Google Scholar

Schoenberg, M. (1980). Elastic wave behavior across linear slip interfaces. J. Acoust. Soc. Am. 68, 1516–1521. doi:10.1121/1.385077

CrossRef Full Text | Google Scholar

Schoenberg, M., and Douma, J. (1988). Elastic wave propagation in media with parallel fractures and aligned cracks. Geophys. Prospect. 36, 571–590. doi:10.1111/j.1365-2478.1988.tb02181.x

CrossRef Full Text | Google Scholar

Schoenberg, M., and Sayers, C. M. (1995). Seismic anisotropy of fractured rock. Geophysics 60, 204–211. doi:10.1190/1.1443748

CrossRef Full Text | Google Scholar

Shaw, R. K., and Sen, M. K. (2006). Use of AVOA data to estimate fluid indicator in a vertically fractured medium. Geophysics 71, C15–C24. doi:10.1190/1.2194896

CrossRef Full Text | Google Scholar

Xie, L., Huang, L., Wang, W., Lin, L., and Pan, X. (2024). Fourier-coefficients-based multiscale seismic inversion for elastic and fracture parameters in frequency domain. IEEE Trans. Geosci. Remote Sens. 62, 1–7. doi:10.1109/tgrs.2024.3361654

CrossRef Full Text | Google Scholar

Zhang, G., Li, L., Pan, X., Li, H., Liu, J., and Han, L. (2020). Azimuthal Fourier coefficient inversion for horizontal and vertical fracture characterization in weakly orthorhombic media. Geophysics 85, C199–C214. doi:10.1190/geo2019-0797.1

CrossRef Full Text | Google Scholar