Due to the highly nonlinear nature of the exact Zoeppritz equation, it is challenging to apply it to practical seismic inversion, and scholars have developed various forms of approximate equations for the reflection coefficients (Shuey, 1985; Fatti and Smith, 1994; Zong et al., 2012a; Zong et al., 2012b; Zong et al., 2012c; Zhang S. et al., 2017; Li et al., 2022). One of the most commonly used approximate equations is the Zoeppritz approximation equation proposed by Aki and Richards (2001): R P P A ki θ = 1 2 sec 2 ⁡ θ Δ V P V P ¯ − 4 γ ⁡ sin 2 ⁡ θ Δ V S V S ¯ + 1 2 1 − 4 γ ⁡ sin 2 ⁡ θ Δ ρ ρ ¯ (1) where θ is the angle of incidence; Δ V P , Δ V S and Δ ρ denote the differences of the P-wave velocity V P , S-wave velocity, V S and density ρ values across the boundary, and: γ = V S ¯ V P ¯ 2 (2) where V P ¯ and V S ¯ represent the average velocity of the P-wave and S-wave, respectively.

Constructing rational forms of fluid indicators plays a crucial role in reservoir fluid prediction, the most common forms being the Lamé parameters and their combinations with other elastic parameters. Goodway et al. (1997) established a fluid indicator (expressed as the product of the first Lamé constant and the density), which can sensitively characterize the rock feature and is a commonly used elastic parameter for fluid prediction. Therefore, the fluid indicator of the gas-bearing shale reservoir can be expressed as: F = λ ρ (3)

Young’s modulus E, among other commonly used brittleness indices, is a standard parameter for characterizing the brittleness of shale reservoir rocks. Ongoing research has shown that while Young’s modulus can reflect the influence of brittle minerals like quartz in shale, it has been observed that the presence of organic matter, fluids, and porosity in gas-bearing shale reservoirs can significantly weaken the sensitivity to brittle shales. Consequently, the effectiveness of Young’s modulus as an indicator of shale reservoir brittleness is limited (Kumar et al., 2012; Pan et al., 2020; Ma et al., 2023). To address this, Chen et al. (2014) proposed a new brittleness index, BI = E/λ. Through rock physics modeling, this index has been shown to effectively capture the content of minerals, including quartz, and the combined effects of organic content, fluids, and porosity in shale reservoirs, making it more effective in characterizing gas-bearing brittle rocks. To facilitate the derivation of the approximate equation that follows, we express BI in terms of Poisson’s ratio alone: B I = E λ = 1 + σ 1 − 2 σ σ (4)

First, the V P and V S are related to the P-wave modulus M and S-wave modulus μ in an isotropic elastic medium as follows: V P = M ρ , V S = μ ρ (5)

According to the chain rule of multivariable calculus, we can obtain the following equations about Δ V P and Δ V S : Δ V P = Δ M ∂ V P ∂ M + Δ ρ ∂ V P ∂ ρ = 1 2 Δ M 1 ρ M − 1 2 − 1 2 Δ ρ M ρ − 3 2 (6) Δ V S = Δ M ∂ V S ∂ μ + Δ ρ ∂ V S ∂ ρ = 1 2 Δ μ 1 ρ μ − 1 2 − 1 2 Δ ρ μ ρ − 3 2 (7)

Dividing both sides of Equations 6, 7 by V P and V S , respectively. We obtain the following equations: Δ V P V P = 1 2 Δ M M − Δ ρ ρ (8) Δ V S V S = 1 2 Δ μ μ − Δ ρ ρ (9)

In isotropic media, the M and μ have the following relationship: M = E 1 − σ 1 + σ 1 − 2 σ (10) μ = E 2 1 + σ . (11)

Then, we rewrite the fluid indicator F as follows: F = E σ ρ 1 + σ 1 − 2 σ (12)

Introduce intermediate variables Η A , Η B and Η C such that they are: Η A = 1 − 2 σ σ , Η B = 1 + σ , Η C = 1 − σ 1 − 2 σ (13)

Substituting Equation 13 into Equations 4, 10–12, we get the following equations F = E ρ Η A Η B , B I = Η A Η B , M = E Η C Η B , μ = E Η B (14)

Similarly, same as Equations 8, 9, we get the following equations: Δ F F ¯ = Δ E E ¯ − Δ Η A Η ¯ A − Δ Η B Η ¯ B + Δ ρ ρ ¯ , Δ B I B I ¯ = Δ Η A Η ¯ A + Δ Η B Η ¯ B (15) Δ M M ¯ = Δ E E ¯ + Δ Η C Η ¯ C − Δ Η B Η ¯ B , Δ μ μ ¯ = Δ E E ¯ − Δ Η B Η ¯ B (16)

The term about Δ Η A Η ¯ A has the following definition: Δ Η A Η ¯ A = 2 σ 1 − σ 2 σ 1 + σ 2 − 4 σ 1 σ 2 (17) where σ 1 and σ 2 denote the Poisson’s ratios of the upper layer and lower layer, respectively, with the following expressions: σ 1 = σ − Δ σ 2 , σ 2 = σ + Δ σ 2 (18)

Substituting Equation 18 into Equation 17, assuming relatively weak variations in elastic parameters when crossing the elastic interface, the higher order terms about Δ σ 2 can be approximated as zero neglected, we can obtain as: Δ Η A Η ¯ A = 2 σ 4 σ 2 − 2 σ + Δ σ 2 Δ σ σ ¯ = 1 2 σ − 1 Δ σ σ ¯ (19)

Similarly, we can get results for the terms about Δ Η B Η ¯ B and Δ Η C Η ¯ C : Δ Η B Η ¯ B = σ σ + 1 Δ σ σ ¯ (20) Δ Η C Η ¯ C = σ 2 σ − 1 σ − 1 Δ σ σ ¯ (21)

Next, substituting Equations 19–21 into Equations 15, 16 yields: Δ F F ¯ + Δ B I B I ¯ = Δ E E ¯ + Δ ρ ρ ¯ , Δ B I B I ¯ = 3 σ 2 σ − 1 σ + 1 Δ σ σ ¯ (22) Δ M M ¯ = Δ E E ¯ + − 2 σ 3 + 4 σ 2 2 σ − 1 σ − 1 σ + 1 Δ σ σ ¯ , Δ μ μ ¯ = Δ E E ¯ − σ σ + 1 Δ σ σ ¯ (23)

Meanwhile, a Poisson’s ratio relationship subscript exists as: σ = 1 − 2 γ 2 1 − γ (24)

By observing the above equations, we first substitute Equations 8, 9, 23 into Equation 1, and then through Equation 22, we replace the E term and σ term in it with the F term and the BI term, and replace σ in it with Equation 24. Finally, a new approximate equation that allows to obtain both the fluid indicator and the brittleness index has been established: R P P F B D θ = Κ F θ Δ F F ¯ + Κ B I θ Δ B I B I ¯ + Κ ρ θ Δ ρ ρ ¯ (25) with Κ F θ = 1 4 sec 2 ⁡ θ − 2 γ ⁡ sin 2 ⁡ θ Κ B I θ = 4 γ 2 − 7 γ + 3 6 γ 2 − 8 γ + 3 1 2 γ ⁡ sec 2 ⁡ θ − 2 γ ⁡ sin 2 ⁡ θ Κ ρ θ = 1 2 1 − sec 2 ⁡ θ + 4 γ ⁡ sin 2 ⁡ θ (26)

To verify the accuracy of the newly proposed approximate equation in Equation 25, we tested and analyzed it using the three-layer gas-bearing sandstone and shale model given by Goodway et al. (1997) based on the measured data, and the relevant elastic parameters are shown in Table 1.

Based on the model mentioned above, we calculated the reflection coefficients and the residuals between the approximate and exact equations for different interfaces at incidence angles from 0° to 40° using the exact Zoeppritz equation, the Aki-Richard equation, and our proposed equation, respectively.

Figures 1A, 2A show the response curves in the PP-wave reflection coefficients of the three equations at the reflection interfaces of overlying sandstone and underlying shale and the reflection interfaces of overlying shale and underlying gas-bearing sandstone, respectively. Figures 1B, 2B show the variation in the difference between the reflection coefficients calculated using the Aki-Richards equation and our proposed equation, compared to those calculated using the exact Zoeppritz equation, for a range of incidence angles from 0° to 40°. From Figures 1, 2, we can find that the new approximate equation of reflection coefficient and the exact Zoeppritz equation are the same within the range of 40° of incidence angle, and their errors are within the acceptable range, which can be used in the seismic inversion of AVO elastic parameters.

In geophysical problems, a noisy seismic signal can be represented by the reflection coefficient Equation 25 convolution with the seismic wavelet: S = w * R P P F B D + N ∼ (27)

Convert the above equation into a matrix multiplication operation: S = W R + N ∼ (28) where S ∈ R m × n denotes seismic data and R ∈ R m × n denotes the reflection coefficient, N ∼ ∈ R m × n denotes random noise; m and n are the number of sampling points and traces of seismic data, respectively. W is the wavelet kernels matrix obtained by time-shifting the seismic wavelet Q: W = Q 1 0 ⋯ 0 ⋮ Q 1 ⋱ ⋮ Q i ⋮ ⋱ 0 0 Q i ⋮ Q 1 ⋮ ⋱ ⋱ ⋮ 0 ⋯ ⋱ Q i m × m − 1 (29)

A linear relationship between R and the P-wave impedance Z can be established (Yilmaz, 2001), as shown in the following equation: R i , j = Z i + 1 , j − Z i , j Z i + 1 , j + Z i , j ≈ 1 2 ln ⁡ Z i + 1 , j − ln ⁡ Z i , j (30) where the symbol “ln” denotes the logarithmic operator; the subscripts “i and j” indicate the sampling point locations. Based on Equation 30, note L = ln ⁡ Z , and introduce the difference matrix D to realize the subtraction operation in Equation 30 to obtain the reflection coefficient expression as follows: R = 1 2 Δ Z Z = 1 2 D L (31) where D represents the difference matrix: D = − 1 1 0 ⋯ 0 0 − 1 1 ⋱ ⋮ ⋮ ⋱ ⋱ ⋱ 0 0 ⋯ ⋯ − 1 1 m − 1 × m (32)

The relevant conversion of the elastic parameters in the reflection coefficient approximation Equation 25 can be further denoted as: Δ F F = D L F , Δ B I B I = D L BI , Δ ρ ρ = D L ρ (33) where L F = ln ⁡ F , L BI = ln ⁡ B I , L ρ = ln ⁡ ρ .

In summary, an approximate equation for the seismic record S can be obtained: S θ = Κ F W D L F + Κ B I W D L BI + Κ ρ W D L ρ (34)

If the number of angles is N, write Equation 34 in matrix form for different incidence angles: S θ 1 S θ 2 ⋮ S θ N = Κ F θ 1 W D Κ B I θ 1 W D Κ ρ θ 1 W D Κ F θ 2 W D Κ B I θ 2 W D Κ ρ θ 2 W D ⋮ ⋮ ⋮ Κ F θ N W D Κ B I θ N W D Κ ρ θ N W D L F L BI L ρ (35)

This is a forward model of seismic records d that can be simplified as follows: d = G r (36) in which d = S θ 1 , S θ 2 , ⋯ S θ N T N × m × n (37) G = Κ F θ 1 W D Κ B I θ 1 W D Κ ρ θ 1 W D Κ F θ 2 W D Κ B I θ 2 W D Κ ρ θ 2 W D ⋮ ⋮ ⋮ Κ F θ N W D Κ B I θ N W D Κ ρ θ N W D N × m × n × 3 (38) r = L F , L BI , L ρ T 3 × m × n (39) where G denotes the forward operator; r is the target elasticity parameter obtained from the inversion.

We can obtain the primary objective function J r of the inversion from the above forward model: J r = min r d − G r 2 2 (40) where 2 denotes the L 2 norm.

Then, to improve the stability of the inversion process and to reduce the multiplicity of solutions in seismic inversion, we added the initial model constraint terms and sparse constraint term into Equation 40: J r = min r d − G r 2 2 + μ 1 L F − l F 2 2 + μ 2 L B I − l B I 2 2 + μ 3 L ρ − l ρ 2 2 + λ D y r 1 (41) where 1 denotes the L 1 norm; l F , l BI , l ρ are the results of logarithmic operations of the initial models for the fluid indicator F, brittleness index BI, and density ρ , μ 1 , μ 2 , μ 3 are the initial model constraint parameters for elasticity parameters, and λ is the constraint parameter for the sparsity term, and D y denotes the vertical difference operator.

In Equation 41, it can be observed that the sparse constraint term of the objective function only includes differential information in the vertical direction. However, anisotropic total variation (ATV) is a sparse constraint that extracts sparse information from seismic data in both vertical and horizontal directions, which has achieved satisfying seismic inversion results (Wu et al., 2020; Zhao et al., 2023).

To fully exploit sparse information, we use the L _P (0 < p < 1) norm instead of the L ₁ norm for sparse constraints. As a result, we can give the expression for anisotropic total variation constrained by the L _P norm (ATpV) as follows: A p T V r = r D x p p + D y r p p (42) where p p denotes the L P norm, r = L F , L BI , L ρ T ∈ R 3 m × n is the parameter to be inverted, D x is the horizontal difference operator, D y 0 is part of the vertical difference operator, D y is the vertical difference operator: D x = − 1 0 ⋯ 0 1 1 − 1 ⋱ ⋮ 0 0 1 ⋱ ⋮ ⋮ ⋮ ⋱ ⋱ − 1 0 0 0 ⋯ 1 − 1 n × n − 1 , D y 0 = − 1 1 0 ⋯ 0 0 − 1 1 ⋱ ⋮ ⋮ ⋱ ⋱ ⋱ 0 0 ⋯ ⋯ − 1 1 0 0 ⋯ 0 − 1 m − 1 × m , D y = D y 0 0 0 0 D y 0 0 0 0 D y 0 3 m − 1 × m (43)

Introducing the anisotropic total variation as a sparse regularization constraint in Equation 41, the inversion objective function can be given as follows: J r = min r d − G r 2 2 + μ 1 L F − l F 2 2 + μ 2 L B I − l B I 2 2 + μ 3 L ρ − l ρ 2 2 + λ r D x p p + D y r p p (44)

Given that Equation 44 involves the L ₂ and L _P norms, conventional methods cannot solve it directly. Therefore, we introduced the Alternating Direction Method of Multipliers (ADMM), which can decompose the inversion objective function into several easier-to-solve sub-objectives. The optimal solution can efficiently be obtained for the original objective function by alternating among subproblems (He et al., 2022; Zhao et al., 2024).

First, the objective function is transformed into a constrained optimization problem by introducing the Lagrange multiplier term in Equation 44, replacing r x and r y with r D x and D y r , respectively, to obtain: J r = min r d − G r 2 2 + μ 1 L F − l F 2 2 + μ 2 L B I − l B I 2 2 + μ 3 L ρ − l ρ 2 2 + λ r x p p + r p p   s . t r x = r D x , r y = D y r (45)

Next, the above optimization problem is converted into an unconstrained optimization problem by introducing dual terms C x and C y in the objective function: J r = min r d − G r 2 2 + μ 1 L F − l F 2 2 + μ 2 L B I − l B I 2 2 + μ 3 L ρ − l ρ 2 2 + λ r x p p + r y p p + η r x − r D x − C x 2 2 + η r y − D y r − C y 2 2 (46) where η represents the weight coefficient of the dual term.

The objective function Equation 46 is decomposed into sub-functions related to r , r x , r y , C x , and C y .

The sub-problem related for r is: J 1 r = min d − G r 2 2 + μ 1 L F − l F 2 2 + μ 2 L B I − l B I 2 2 + μ 3 L ρ − l ρ 2 2 + η r x − r D x − C x 2 2 + η r y − D y r − C y 2 2 (47)

Equation 47 is a convex optimization problem that can be solved directly. Compute the gradient concerning r and simplify: G T G + μ E + η D y T D y r + r η D x T D x = G T d + μ r 0 + η D y T r y i − C y i + r x i − C x i D x T (48) where E is the unit matrix, r 0 = l F l B I l ρ is the initial model for the inversion parameters, μ = μ 1 0 0 0 μ 2 0 0 0 μ 3 .

Given that Equation 48 is difficult to solve for r directly, Gholami (2015) found that it can be transformed into a particular form: the 2D Sylvester equation then facilitates solving for r , the equation as follows: A x + x B = C (49)

Then, Equation 48 can be organized into the form of Equation 49 to obtain the inversion parameter r : A = G T G + μ E + η D y T D y , B = η D x T D x , C = G T d + μ r 0 + η D y T r y i − C y i + r x i − C x i D x T (50)

The next sub-problems on r x and r y are: J 2 r x = min r x λ r x p p + η r x − r D x − C x 2 2 (51) J 3 r y = min r y λ r y p p + η r y − r D y − C y 2 2 (52)

Equations 51, 52 contain L _P norms and cannot be solved by traditional convex optimization algorithms. So, Chartrand and Yin (2016) used the L _P -norm iterative shrinkage thresholding algorithm (ISTA) to solve such problems: r x i + 1 = max r i + 1 D x + C x i − η λ p − 2 · r i + 1 D x + C x i p − 1 , 0 · s i g n r i + 1 D x + C x i (53) r y i + 1 = max D y r i + 1 + C y i − η λ p − 2 · D y r i + 1 + C y i p − 1 , 0 · s i g n D y r i + 1 + C y i (54) in which s i g n x = 1   x > 0 0   x = 0 − 1   x < 0 (55)

The sub-problems for C x and C y can be expressed as: J 4 C x i + 1 = min C x η r x i + 1 − r D x − C x i 2 2 (56) J 5 C y i + 1 = min C y η r y i + 1 − D y r − C y i 2 2 (57)

The solutions to Equations 56, 57 can be found by computing the gradient and setting it to zero: C x i + 1 = C x i + r i + 1 D x − r x i + 1 (58) C y i + 1 = C y i + D y r i + 1 − r y i + 1 (59)

Eventually, the framework of the inversion technique is summarized in Algorithm 1. Meanwhile, we show the specific workflow of this paper to obtain the inversion target parameters in Figure 3.

Algorithm 1

Inversion Pseudocode based on ADMM of the ATpV method.

Algorithm: Multitrace pre-stack inversion of elasticity parameters using ADMM and ISTA.

To demonstrate the reliability and practicality of the simultaneous inversion method for fluid indicator and brittleness index presented in this paper, we conducted three sets of tests using synthetic seismic data derived from actual logging data. First, we obtained the multi-angle pre-stack seismic data by convolving the theoretical Ricker wavelet with the reflection coefficients computed from the actual log data based on the Zoeppritz equation, as shown in Figure 4. Then, to assess the performance of the inversion method under different noise conditions, we introduced various levels of random noise into the data. Moreover, Figures 4B, C show noisy datasets with SNR (Signal-to-Noise Ratio) of 5 and 2, respectively.

Figures 5–7 present the results obtained from the synthetic data, the inverted data: fluid indicator F, brittleness index BI, and density (red solid curves), the initial models (green dashed curves), and the original well data (black solid curves), and the distribution of error between them. From Figure 5, in the noise-free case, all elastic parameter inversion results are well inverted and agree very well with the original well curve tracks. Figures 6, 7 show the inversion results when the SNR of the synthetic data is 5 and 2, respectively. We found that when the low-frequency curve is smooth and the model resolution is low, the gap between the obtained results and the actual logging data is more minor. The resolution of the inversion results can be guaranteed. Even at an SNR of 2, the inversion results obtained in this paper can also be matched well with the original well curves and have good noise resistance.

Then, to quantitatively evaluate our results, we use Equations 60, 61 to calculate the inversion result errors ε and RMSE (root-mean-square-errors). ε = X i n v − X w e l l X w e l l × 100 % (60) R M S E = ∑ i = 1 N X i i n v − X i w e l l 2 N (61) where X _inv and X _well represent the inversion result and original well curve, respectively. N is the number of sampling points.

We find that the inversion error for the brittleness index is more significant than that for the indicator index and density in the noise-free case. Still, the inversion error for all results is less than 5%, giving high accuracy. With the increase of noise, the inversion errors and RMSE of the brittleness index significantly increase. However, the errors and RMSE of fluid indicator and density vary within normal limits. Overall, the synthetic data test results can demonstrate that the methodology of this paper is theoretically feasible.

The validity and practicability of the methodology proposed in this paper are verified through a field shale gas reservoir work area. This paper’s field study area (as shown in Figure 8) is located in the Sichuan Basin, China, and the primary target formation is the Silurian Longmaxi Formation. Well A in the work area drilled a high-quality shale gas reservoir, but the reservoir thickness is relatively small. The main target reservoir and lithology of well A are shown in Figure 9. The upper part of T1 is a limy siltstone with interbedded gray and sand. The section from T1 to T2 is a shale Section 1 with locally thinly interbedded argillaceous siltstone. The reservoir section from T2 to T3 is characterized by high porosity and substantial organic matter content, making it a high-quality, gas-bearing shale reservoir (Section 2) with significant hydrocarbon production potential. The lower part of T3 is the underlying limestone layer. In addition, based on actual logging curve values, such as elasticity curves for velocity and density, we can calculate the fluid indicator curve and brittleness index curve using Equations 3, 4. We can find that the fluid indicator curve and the brittleness index curve show apparent low-value and high-value anomalies in the shale 2 section; it is a favorable location for predicting shale reservoirs.

Figure 10 shows the stacked seismic data with center incidence angles of 5° (0°–10°) in Figure 10A, 15° (10°–20°) in Figure 10B, and 25° (20°–30°) in Figure 10C, respectively. The profile has 285 traces of seismic data, and the red vertical line represents the path of well A, which is located on trace 112. In addition, the low-frequency components often lacking in the actual seismic data may lead to unstable inversion results. Therefore, a common approach is introducing initial model constraints in the inversion constraint term to make the results more reliable. Utilizing the seismic stratigraphic information, the fluid indicator, brittleness index, and density data from well A were extrapolated, interpolated, and then low-pass filtered to generate the corresponding low-frequency models, as illustrated in Figures 11A–C.

According to the inversion Algorithm 1 and combined with the approximate equation proposed in Equation 25, we inverted the 2D seismic data. Figures 12A–C shows the final 2D inversion results for the target parameters.

In the fluid indicator inversion results, the red color indicates that the fluid identification factor is low, which means that the formation in the red part of the profile is likely to have fluids. In the brittleness index inversion results, the red color indicates high brittleness index values, which means that the formation in the red part of the profile is prone to fracture. Figure 12 shows that the brittleness inversion results have apparent high anomalies in the target layer section (T2 to T3). The fluid indicator inversion results have low anomalies in the target layer section (T2 to T3). The combination of the two inversion results indicates that the formation at this location is brittle and contains fluid, which is in agreement with the recognition in Figure 9. The combination of the inversion results is beneficial to characterize the gas-bearing shale reservoir better and helps us identify the target shale reservoir section.

To demonstrate the advantages and reliability of the direct inversion method for fluid indicator F and brittleness index BI proposed in this paper, we followed the inversion strategy outlined in Figure 3. Initially, according to the method proposed in this paper, we invert three elastic parameters V P , V S , and ρ, using the Aki-Richards approximation equation. Subsequently, we indirectly calculated the F and BI at the well location (trace 112) using Equations 3, 4, as illustrated by the green curves in Figure 13. Eventually, we extracted the inversion results of F and BI for the location of well A from Figure 12 and represented them as red curves in Figure 13.

From Figures 13A, B, we found that the indirect inversion results accumulate errors in the calculation, which can lead to inaccurate results at the reservoir location (black dotted box). In addition, in several places in Figure 13 (black arrows), we can also find that the indirect prediction results (green solid curves) significantly exceed the actual well data (black solid curves). In contrast, the simultaneous inversion results for the fluid indicator and brittleness index (red solid curves) proposed in this paper strongly align with the original well curves. Furthermore, the residuals and RMSE analysis in Figure 14 indicate that the simultaneous inversion method yields much smaller residuals and lower RMSE than the indirect calculation method. Therefore, the approximation equation proposed in this paper enables more accurate and reliable simultaneous inversion of fluid indicator and brittleness index.

The brittleness index BI, as an engineering sweet spot, is employed to characterize the fracturing ability of rock, which plays an instrumental guidance role in shale reservoir hydraulic fracturing construction, and the shale gas fluid indicator F, as a geological sweet spot has been revealed to be used to characterize hydrocarbon-bearing enriched zones. The combination of engineering and geological sweet spots is dedicated to the well deployment, horizontal well trajectory design, and effective shale gas production and extraction promotion. Based on the simultaneous seismic inversion method of shale gas fluid factor and brittleness index realized in this article, 3D pre-stacked seismic inversion at the fracturing platform has been utilized to perform the spatial distribution estimation with the above double-sweet spot parameters. The 3D seismic data was acquired by 230 Inlines and 270 Xlines, where wells H1, H3, and H5 are drilled as horizontal fractured wells in the shale reservoir, detecting numerous micro-seismic events at different periods. The inversion result slices of the fluid indicator and brittleness index along the destination layer are displayed in Figures 15A, B, respectively. The dashed circle closures pointed by white arrows reveal the hydrocarbon-enriched favorable fracturing zones, where the low fluid indicator spreading in Figure 15A represents enrichment of shale gas resources, and the high brittleness index distribution in Figure 15B represents satisfactory rock brittleness, which is beneficial to fracture fragmentation and facilitates hydraulic fracturing construction. The black dashed line Ft1 represents the projection of the seismic wave artificial fault interpretation conclusion on the slice. It can be seen that the fluid indicator exhibits uninteresting high values because the fault structure leads to increased stratum permeability, which is unfavorable to shale gas storage. The brittleness index presents low values, indicating that the rock brittleness at the fault location may be inferior because the fracture zone has been filled with more nonbrittle minerals, such as clay. In addition, in Figure 15C, numerous micro-seismic events in Figure 15D detected along three horizontal wells H1, H3, and H5 demonstrate the satisfying positive correlation with the high brittleness index in the fracturing favorable zone, which indicates that the rock is brittle enough and can be easily fractured to induce micro-seismic events. The practical application demonstrates that this article’s simultaneous seismic inversion method of fluid indicator and brittleness index can effectively provide reliable guidance information for unconventional shale gas production.