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: δ N = λ + 2 μ K N 1 + λ + 2 μ K N (1) δ T = μ K T 1 + μ K T (2) where λ and μ 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 δ N and δ 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: δ N = 4 ρ f 3 g 1 − g + 3 α K ′ + 4 / 3 μ ′ / π μ (3) δ T = 16 ρ f 3 3 − 2 g + 12 α μ ′ / π μ (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, α = b / a ≪ 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 ρ f = ν a 3 represents the fracture density, in which ν 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 P P H T I θ , φ = a V P θ R V P + a V S θ R V S + a ρ θ R ρ + a δ N θ , φ R δ N + a δ T θ , φ R δ T (5) where a V P θ = sec 2 ⁡ θ (6) a V S θ = − 8 g ⁡ sin 2 ⁡ θ (7) a ρ θ = 1 − 4 g ⁡ sin 2 ⁡ θ (8) a δ N θ , φ = − 1 2 sec 2 ⁡ θ 1 − 2 g sin 2 ⁡ θ ⁡ sin 2 ⁡ φ + cos 2 ⁡ θ 2 (9) a δ T θ , φ = 2 g ⁡ sin 2 ⁡ θ ⁡ cos 2 ⁡ φ 1 − tan 2 ⁡ θ ⁡ sin 2 ⁡ φ (10) and R V P = Δ V P / 2 V P , R V S = Δ V S / 2 V S , R ρ = Δ ρ / 2 ρ , R δ N = δ N 2 − δ N 1 / 2 , and R δ T = δ T 2 − δ T 1 / 2 . Here θ and φ = φ o b s − φ f denote the angle of incidence and the azimuth angle between the observed azimuth φ o b s and the azimuth of fracture symmetry axis φ f , respectively; R δ N and R δ 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 P P H T I θ , φ = R 0 θ + ∑ ξ = 1 ∞ R ξ θ cos ξ φ − φ ξ θ (11) where R ξ θ represents Fourier coefficients varying with the angle of incidence, including the magnitude R ξ θ and the phase φ ξ θ . Therefore, the azimuthal reflectivity at a specified incidence angle θ 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 ξ , spatial periodicity scaled by integer index ξ , and phase retardation controlled by term φ ξ . 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 P P H T I θ , φ = R 0 θ + R 2 θ cos 2 φ − φ 2 θ + R 4 θ cos 4 φ − φ 4 θ (12) with magnitudes R 0 θ = a V P θ R V P + a V S θ R V S + a ρ θ R ρ + a δ N 0 θ R δ N + a δ T 0 θ R δ T (13) R 2 θ = a δ N 2 θ R δ N + a δ T 2 θ R δ T (14) R 4 θ = a δ N 4 θ R δ N + a δ T 4 θ R δ T (15) where a δ N 0 θ = − 1 2 sec 2 ⁡ θ + 2 g 1 − g + 1 2 sin 2 ⁡ θ + 1 2 − 3 8 g sin 2 ⁡ θ ⁡ tan 2 ⁡ θ (16) a δ T 0 θ = g ⁡ sin 2 ⁡ θ 1 − 1 4 tan 2 ⁡ θ (17) a δ N 2 θ = g ⁡ sin 2 ⁡ θ 2 g − 1 + g − 1 tan 2 ⁡ θ (18) a δ T 2 θ = g ⁡ sin 2 ⁡ θ (19) a δ N 4 θ = − 1 4 g 2 ⁡ sin 2 ⁡ θ ⁡ tan 2 ⁡ θ (20) a δ T 4 θ = 1 4 g ⁡ sin 2 ⁡ θ ⁡ tan 2 ⁡ θ (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., R P P H T I θ , φ = ∑ ξ = 0 ∞ α ξ θ cos ⁡ ξ φ + β ξ θ sin ⁡ ξ φ = α 0 θ + α 2 θ cos ⁡ 2 φ + β 2 θ sin ⁡ 2 φ +   α 4 θ cos ⁡ 4 φ + β 4 θ sin ⁡ 4 φ (22) with coefficients α 0 θ = a V P θ R V P + a V S θ R V S + a ρ θ R ρ + a δ N 0 θ R δ N + a δ T 0 θ R δ T (23) α 2 θ = a δ N 2 θ cos ⁡ 2 φ f R δ N + a δ T 2 θ cos ⁡ 2 φ f R δ T (24) β 2 θ = a δ N 2 θ sin ⁡ 2 φ f R δ N + a δ T 2 θ sin ⁡ 2 φ f R δ T (25) α 4 θ = a δ N 4 θ cos ⁡ 4 φ f R δ N + a δ T 4 θ cos ⁡ 4 φ f R δ T (26) β 4 θ = a δ N 4 θ sin ⁡ 4 φ f R δ N + a δ T 4 θ sin ⁡ 4 φ f R δ T (27) and the coefficients α ξ θ and β ξ θ in front of the cosine and sine functions with K-samples azimuthal seismic data are defined as follows: α ξ θ = 1 π ∑ k = 1 K R k θ , φ cos ⁡ ξ φ d φ (28) β ξ θ = 1 π ∑ k = 1 K R ξ θ , φ sin ⁡ ξ φ d φ (29)

The above parameters can be transformed to magnitude and phase as R ξ θ = α ξ 2 θ + β ξ 2 θ (30) φ ξ θ = 1 ξ arctan β ξ θ α ξ θ (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 ρ , R δ N , and R δ 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.

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 d θ , φ can be obtained via the convolution model (Pan and Zhang, 2019) and the estimated seismic wavelets as d θ , φ = R P P H T I θ , φ · W θ , φ (32) where W θ , φ 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 ξ θ = 1 π ∑ k = 1 5 d k θ , φ cos ⁡ ξ φ d φ , ξ = 0 , 2 , 4 (33) B ξ θ = 1 π ∑ k = 1 5 d k θ , φ sin ⁡ ξ φ d φ , ξ = 0 , 2 , 4 (34) where A ξ θ and B ξ θ represent the cosine and sine components of the Fourier coefficients for the kth-azimuth seismic data d k θ , φ .

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 α ξ θ and β ξ θ 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 d 1 = G 1 m 1 (35) where d 1 = A 2 θ 1 B 2 θ 1 L A 2 θ 3 B 2 θ 3 T (36) G 1 = W θ 1 a δ N 2 θ 1 cos ⁡ 2 φ f W θ 1 a δ N 2 θ 1 sin ⁡ 2 φ f L W θ 3 a δ N 2 θ 3 cos ⁡ 2 φ f W θ 3 a δ N 2 θ 3 sin ⁡ 2 φ f W θ 1 a δ T 2 θ 1 cos ⁡ 2 φ f W θ 1 a δ T 2 θ 1 sin ⁡ 2 φ f L W θ 3 a δ T 2 θ 3 cos ⁡ 2 φ f W θ 3 a δ T 2 θ 3 sin ⁡ 2 φ f T (37) m 1 = R δ N R δ T T (38) and A 2 θ j = A 2 θ j , t 1 L A 2 θ j , t M T , j = 1 , 2 , 3 (39) B 2 θ j = B 2 θ j , t 1 L B 2 θ j , t M T , j = 1 , 2 , 3 (40) a δ N 2 θ j = d i a g a δ N 2 θ j , t 1 L a δ N 2 θ j , t M (41) a δ T 2 θ j = d i a g a δ T 2 θ j , t 1 L a δ T 2 θ j , t M (42) R δ N = R δ N t 1 L R δ N t M T (43) R δ T = R δ T t 1 L R δ T t M 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 m 1 d 1 can be written as a joint PDF of the prior PDF p m 1 and the likelihood function p d 1 m 1 : p m 1 d 1 ∝ p m 1 · p d 1 m 1 (45) where p · 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 m 1 d 1 ∝ ∏ i = 1 2 M 1 1 + m 1 2 / σ m 1 2 · ⁡ exp − G 1 m 1 − d 1 2 2 σ d 1 2 (46) where σ d 1 2 and σ 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 J m 1 = G 1 m 1 − d 1 2 + 2 σ d 1 2 ∑ i = 1 2 M ln 1 + m 1 2 / σ m 1 2 + ∑ i = 1 2 τ i χ i − P i m 1 i 2 (47) where τ i represents the regularization coefficients of normal and tangential fracture weaknesses, P i = ∫ t 1 t i · d t represents the integral operator, and χ i = 1 2 ln m 1 i m 1 i 0 , in which m 1 i 0 represents the initial smoothing model of fracture weaknesses.

Minimizing the objective function of Equation 47 yields ψ = Γ m 1 , (48) where Γ = G 1 T G 1 + 2 σ d 1 2 σ m 1 2 C C a u c h y + ∑ i = 1 2 τ i P i T P i (49) ψ = G 1 T d 1 + ∑ i = 1 2 τ i P i T χ i (50) where C C a u c h y 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: δ N t i = δ N t 1 + ∑ t 1 t i − 1 R δ N τ (51) δ T t i = δ T t 1 + ∑ t 1 t i − 1 R δ T τ (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 d 2 = G 2 m 2 (53) where d 2 = A 0 θ 1 L A 0 θ 3 T (54) G 2 = W θ 1 a V P θ 1 L W θ 3 a V P θ 3 W θ 1 a V S θ 1 L W θ 3 a V S θ 3 W θ 1 a ρ θ 1 L W θ 3 a ρ θ 3 W θ 1 a δ N 0 θ 1 L W θ 3 a δ N 0 θ 3 W θ 1 a δ T 0 θ 1 L W θ 3 a δ T 0 θ 3 T (55) m 2 = R V P R V S R ρ R δ N R δ T T (56) and A 0 θ j = A 0 θ j , t 1 L A 0 θ j , t M T , j = 1 , 2 , 3 (57) a V P θ j = d i a g a V P θ j , t 1 L a V P θ j , t M (58) a V S θ j = d i a g a V S θ j , t 1 L a V S θ j , t M (59) a ρ θ j = d i a g a ρ θ j , t 1 L a ρ θ j , t M (60) a δ N 0 θ j = d i a g a δ N 0 θ j , t 1 L a δ N 0 θ j , t M (61) a δ T 0 θ j = d i a g a δ T 0 θ j , t 1 L a δ T 0 θ j , t M (62) R V P = R V P t 1 L R V P t M T (63) R V S = R V S t 1 L R V S t M T (64) R ρ = R ρ t 1 L R ρ t M 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 t i = V P t 1 exp 2 ∑ t 1 t i R V P τ (66) V S t i = V S t 1 exp 2 ∑ t 1 t i R V S τ (67) ρ t i = ρ t 1 exp 2 ∑ t 1 t i R ρ τ (68)

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.

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.