According to convolution theory, a prestack seismic gather can be regarded as the result of the convolution between reflection coefficients at different angles and the seismic wavelet. Considering an incident angle, the expression of the prestack convolution model matrix can be written as: D θ = W · R m , θ + n θ , (1) where · in the middle represents multiplication, D θ is the seismic record vector at an incident angle of θ , W is the seismic wavelet matrix, R m , θ represents the reflection coefficient sequence calculated from the model parameters m at angle of θ and n θ is the random noise vector.

In prestack inversion, the Aki-Richards approximate equations are commonly used. The expressions for PP- and PS-wave reflection coefficients are as follows (Aki and Richards, 1980): R PP θ = 1 2 sec 2 ⁡ θ Δ V P V P − 4 V S 2 V P 2 sin 2 ⁡ θ Δ V S V S + 1 2 1 − 4 V S 2 V P 2 sin 2 ⁡ θ Δ ρ ρ , (2) R PS θ , φ = tan ⁡ φ V P 2 V S 4 V S 2 V P 2 sin 2 ⁡ θ − 4 V S V P cos ⁡ θ ⁡ cos ⁡ φ Δ V S V S − 1 − 2 V S 2 V P 2 sin 2 ⁡ θ + 2 V S V P cos ⁡ θ ⁡ cos ⁡ φ Δ ρ ρ , (3) where R PP represents the PP-wave reflection coefficient, and R PS represents the PS-wave reflection coefficient. θ is the average of the incident and transmitted angles for the PP-wave, φ is the average of the reflection and transmission angles for the PS-wave. V P , V S and ρ are the average values of P-wave velocity, S-wave velocity, and density on both sides of the elastic interface, respectively. Δ V P , Δ V S and Δ ρ represent the differences in P-wave velocity, S-wave velocity, and density across the elastic interface.

Equations 2, 3 can be simplified into a linear summation of three parameter reflection coefficients and expressed in matrix form as: R PP θ R PS θ , φ = C 1 θ C 2 θ C 3 θ 0 C 4 θ , φ C 5 θ , φ r P r S r ρ , (4) where r P = Δ V P V P , r S = Δ V S V S and r ρ = Δ ρ ρ . C 1 , C 2 , C 3 , C 4 and C 5 are the coefficients related to the angles θ and φ respectively, and are expressed as: C 1 θ = 1 2 sec 2 ⁡ θ , C 2 θ = − 4 V S 2 V P 2 sin 2 ⁡ θ , C 3 θ = 1 2 1 − 4 V S 2 V P 2 sin 2 ⁡ θ , C 4 θ , φ = tan ⁡ φ v p 2 v s 4 V S 2 V P 2 sin 2 ⁡ θ − 4 V S V P cos ⁡ θ ⁡ cos ⁡ φ , C 5 θ , φ = − tan ⁡ φ V P 2 V S 1 − 2 V S 2 V P 2 sin 2 ⁡ θ + 2 V S V P cos ⁡ θ ⁡ cos ⁡ φ . (5)

In Equation 5, the angle φ can be converted to θ using Snell’s Law. Therefore, for multi-channel seismic data with different incident angles, the PP-wave and PS-wave reflection coefficients defined by the approximate equations can be expressed in matrix form as: R PP θ 1 ⋮ R PP θ i R PS θ 1 ⋮ R PS θ i = C 1 θ 1 C 2 θ 1 C 3 θ 1 ⋮ ⋮ ⋮ C 1 θ i C 2 θ i C 3 θ i 0 C 4 θ 1 C 5 θ 1 ⋮ ⋮ ⋮ 0 C 4 θ i C 5 θ i × r P r S r ρ ,    R         C         r (6) where R PP θ i = R PP , 1 θ i R PP , 2 θ i ⋯ R PP , n θ i T represents the PP-wave reflection coefficient vector at an incident angle of θ i , n is the number of sampling points. Similarly, R PS θ i is the PS-wave reflection coefficient vector. r P = r P , 1 r P , 2 ⋯ r P , n T represents the reflection coefficient vector for P-wave velocity. Similarly, r S and r ρ are the reflection coefficient vectors for S-wave velocity and density, respectively. C 1 θ i , C 2 θ i , C 3 θ i , C 4 θ i and C 5 θ i are all diagonal matrices with similar expressions. For an example of C 1 θ i , and its expression can be given by: C 1 θ i = C 1 , 1 θ i ⋯ 0 0 ⋯ C 1 , 2 θ i ⋯ 0 0 ⋯ ⋯ ⋯ 0 0 ⋯ C 1 , n θ i . (7)

According to Equation 1 of convolution model, the forward modelling of prestack seismic gathers at different angles can be expressed in matrix form as: D PP θ 1 ⋮ D PP θ i D PS θ 1 ⋮ D PS θ i = W PP θ 1 ⋱ W PP θ i W PS θ 1 ⋱ W PS θ i × R PP θ 1 ⋮ R PP θ i R PS θ 1 ⋮ R PS θ i D W R (8)

where D PP θ i = D PP , 1 θ i D PP , 2 θ i ⋯ D PP , n θ i T represents the PP-wave seismic data vector at an incident angle of θ i , and similarly, D PS θ i represents the PS-wave seismic data vector. W PP θ i is the PP-wave wavelet matrix corresponding to the angle θ i , and similarly, W PS θ i is the PS-wave wavelet matrix.

According to Equations 6, 8, the noise-contaminated prestack gather forward convolution model can be simplified as: D = W R + n = W C r + n (9)

where D represents the prestack seismic data vector for both PP-wave and PS-wave. W is the large diagonal matrix composed of the wavelet matrices for PP-wave and PS-wave at different angles. C is the coefficient matrix composed of the diagonal matrices C 1 θ i , C 2 θ i , C 3 θ i , C 4 θ i and C 5 θ i . r = r P r S r ρ T is the overall elastic parameter reflection coefficient vector composed of the P-wave velocity, S-wave velocity, and density reflection coefficient vectors. n represents the random noise vector.

The construction of the inversion objective function includes a misfit term and a regularization term. The misfit term describes the residuals between the observed data and the modelled data, and it typically uses the least-squares misfit function due to its computational efficiency. The role of regularization term is to introduce prior constraint information into the inversion process and solve the problem of multiple solutions in inversion. In seismic inversion, the presence of noise and inaccuracies in the wavelet forward operator often make the problem highly ill-posed with significant non-uniqueness. Therefore, selecting an appropriate regularization constraint is crucial for seismic inversion.

Initially, the L ₂-norm constraint was commonly used, which is a type of smoothness constraint. However, in seismic inversion, the primary target is the reflection coefficient series, which tends to be sparse due to the layered nature of the subsurface. To better reflect real conditions, most modern seismic inversion methods introduce sparse constraints to enhance the sparsity of the inversion results. Among different sparse norm regularizations, the L ₀-norm is the best measure of sparsity. However, the optimization problem involving the L ₀-norm is NP-Hard and difficult to solve. As a convex approximation of the L ₀-norm, the L ₁-norm is widely used in seismic inversion. Here, we introduce a sparse norm constraint that has been gradually applied to geophysical inversion in recent years—the L _1-2-norm constraint. It is defined as the difference between the L ₁-norm and the L ₂-norm, and it provides better sparsity than the L ₁-norm.

To illustrate the advantages of the L _1-2-norm constraint, we computed the 3D solution space distributions for different norms, as well as their contour projections on a 2D plane, comparing their sparsity (Wang and Chen, 2022), as shown in Figure 1. The figure not only compares the conventional L ₀, L ₂, and L ₁-norm, but also compares two other norm constraints proposed by scholars to optimize the inversion results, namely, the L _cauchy and L _p-norm. In the 3D solution space, the process of minimizing the objective function can be seen as finding the solution with the smallest regularization term while keeping the misfit function constant, which corresponds to searching for the lowest point on the surface. Therefore, the lowest point of the 3D surface is generally located near the x-axis and y-axis on the 2D contour map. Accordingly, for the distribution of contour lines, the closer they are to the x-axis and y-axis, the sparser the corresponding inversion solution. It can be observed that the L ₀ norm is the sparsest among these constraints, with contour lines distributed along the x-axis and y-axis. Compared to the other norms, the L _1-2-norm’s minimum contour lines are closer to the x-axis and y-axis, making it more similar to the L ₀ norm. Therefore, the sparsity of the L _1-2-norm constraint is superior to the other norm constraints.

Based on the prestack gather forward model in Equation 9 and the least-squares misfit function, the objective function with the L ₁-norm constraint can be expressed as: J r = argmin 1 2 W C r − D 2 2 + λ r 1 , (10)

where λ is the trade-off factor. Here, we introduce the L _1-2-norm constraint, and the objective function is then expressed as: J r = argmin 1 2 W C r − D 2 2 + λ r 1 − α r 2 , (11)

where α ∈ 0 , 1 is a constant that enhances the sparsity of the inversion results. r represents the reflection coefficients of the subsurface elastic parameters. To directly invert for the elastic parameters, a new inversion objective function is constructed based on Equation 11, by introducing a first-order difference matrix and incorporating initial model constraints: J m = argmin 1 2 W C L m − D 2 2 + λ L m 1 − α L m 2 + μ 2 m − m 0 2 2 , (12)

where L is the first-order difference matrix, m is the logarithm of the elastic parameters, m 0 is the initial model.

To solve the optimization problem, a two-step iterative strategy is adopted. In the first step, the Difference of Convex Algorithm (DCA) is used to decompose the objective function. In the second step, the Alternating Direction Method of Multipliers (ADMM) is applied to solve the decomposed problem (Aster et al., 2019). In the first step, the DCA algorithm is used to reformulate the objective function as follows: J m = F m − H m , (13)

where F m = 1 2 W C L m − D 2 2 + λ L m 1 + μ 2 m − m 0 2 2 , H m = λ α L m 2 . After decomposition, the objective function takes the form of the difference between two convex functions, which means that the overall function is also convex. In this case, any locally optimal solution is guaranteed to be a globally optimal solution. For Equation 13, the solution is obtained through an alternating iterative process: y k ∈ α H m k , m k + 1 = argmin F m − H m k + y k , m − m k , (14) where k is the number of iterations, y k represents the gradient of H m at m k , and is expressed as: y k = 0 , if  m k = 0 , λ α L T L m k m k 2 , otherwise . (15)

In addition, in Equation 14, , represents the inner product symbol, and m k + 1 denotes the obtained solution. To solve for m k + 1 , this is equivalent to solving a new optimization problem: m k + 1 = argmin 1 2 W C L m − D 2 2 + λ L m 1 + μ 2 m − m 0 2 2 + m , y k . (16)

For this new problem, the second-step strategy is applied, using the ADMM algorithm for solving it. In the ADMM algorithm, a new auxiliary vector x is introduced. The optimization problem then becomes: m k + 1 = argmin 1 2 W C L m − D 2 2 + λ x 1 + μ 2 m − m 0 2 2 + m , y k , s u b j e c t   t o   L m = x . (17)

Using Lagrange multipliers, Equation 17 can be rewritten in the form of an augmented Lagrangian: L m , x , u = argmin 1 2 W C L m − D 2 2 + λ x 1 + μ 2 m − m 0 2 2 + m , y k , +   u T L m − x + ω 2 L m − x 2 2 , (18)

where u is the Lagrange multiplier vector, and ω > 0 is the penalty parameter that controls the convergence speed of the iterations. For Equation 18, the solution can be obtained through alternating iterations among the three variables (Aster et al., 2019): m k + 1 = L T C T W T W C L + μ I + ω L T L − 1 L T C T W T D − y k + μ m 0 + ω L T x k − u k , x k + 1 = s t h r e s h L m k + 1 + u k , λ / ω , u k + 1 = u k + L m k + 1 − x n + 1 , (19)

Where T represents the transpose symbol of the matrix, s t h r e s h represents the soft-thresholding algorithm, and is expressed as: x k + 1 = L m k + 1 + u k − λ / ω , L m k + 1 + u k > λ / ω , L m k + 1 + u k + λ / ω , L m k + 1 + u k < − λ / ω , 0 , otherwise . (20)

The iterative calculation terminates when the maximum number of iterations is reached or when the iteration error is smaller than a predefined threshold, which can be expressed as: m k + 1 − m k 2 1 + m k + 1 ≤ ε ∪ k > M (21)

where ε > 0 is a given tolerance value, M is the maximum number of iterations.

In summary, Figure 2 presents the workflow of the proposed PP- and PS waves joint inversion method. Based on the previous description, the workflow can be divided into three main parts. The first part involves the combination of PP- and PS-wave seismic data. The second part is the construction of the joint inversion objective function. The third part is the optimization inversion solution of three parameters.

In the 1D multilayer model test, the reflection coefficients of the model exhibit significant sparsity, making it suitable for demonstrating the advantages of the L _1-2-norm constraint. A comparison with the conventional L ₁-norm constraint is conducted to verify that the L _1-2-norm offers superior performance. The model curves are shown in Figure 3a, where the left panel represents P-wave velocity, the middle panel represents S-wave velocity, and the right panel represents density, with a time of 720 ms and a time sampling interval of 2 ms. Using this model, reflection coefficients at different angles are calculated using the approximate equations (Equations 2, 3), with a maximum angle of 40°. These are then convolved with a 40 Hz Ricker wavelet to generate prestack angle gathers for both PP- and PS-wave data. Different noise levels are added to the synthetic gathers, as shown in Figures 3b–d, representing noise-free, SNR = 10, and SNR = 5 conditions for the PP- and PS-wave prestack angle gathers.

The synthetic angle gathers from the multilayer model are used directly for joint inversion tests of PP- and PS-waves data to verify the feasibility of the proposed method. Figure 4 shows the comparison of elastic parameter results obtained using L ₁-norm constraint and L _1-2-norm constraint for joint inversion under different signal-to-noise ratios. The blue dashed line represents the initial inversion model, the black solid line represents the true model, and the red solid line represents the inversion results. Figures 4a–c show the results of the joint inversion using the L ₁-norm constraint. In the noise-free case, the inversion results match the true model closely, but as noise increases, the inversion quality declines. Figures 4d–f show the results of the joint inversion using the L _1-2-norm constraint. In the noise-free case, the inversion results coincide with the true model, and although the inversion quality decreases as noise increases, the results show better overall agreement with the true model compared to the L ₁-norm constraint.

Table 1 provides the correlation coefficients and normalized root mean square errors (NRMSE) between the inversion results from Figure 4 and the true model. The table shows that the L _1-2-norm constraint yields higher correlation coefficients and lower NRMSE than the L ₁-norm constraint, indicating better performance. Figure 5 compares the elastic parameter reflection coefficients calculated from the inversion results. The blue thick line and red thin line represent the true model reflection coefficients and the calculated reflection coefficients from the inversion results, respectively. As seen in the figure, with increasing noise, the L _1-2-norm constraint better suppresses small values while preserving the amplitude of larger values, confirming that the L _1-2-norm constraint provides superior sparsity.

In the well log model test, the model curves are shown in Figure 6a. The left panel shows the P-wave velocity, the middle panel shows the S-wave velocity, and the right panel shows the density, with a time of 600 ms and a time sampling interval of 2 ms. Using this model and applying approximate equations, reflection coefficients at different angles are calculated, with the maximum angle being 40°. These are then convolved with a Ricker wavelet of 40 Hz dominant frequency to generate the prestack angle gathers for both PP-wave and PS-wave data. Different levels of noise are added to the synthetic gathers, as shown in Figures 6b–d, representing noise-free, SNR = 10, and SNR = 5 conditions for the PP-wave and PS-wave prestack angle gathers, respectively.

The synthetic angle gathers of PP- and PS-waves data from the well-log model are used directly for inversion tests. The results of only PP-wave inversion and PP-PS joint inversion are compared. Figure 7 shows the elastic parameter results under different signal-to-noise ratio, comparing the outcomes of only PP-wave inversion and PP-PS joint inversion. In the noise-free case, both only PP-wave and PP-PS joint inversion results match the true model. As the noise level increases, the inversion quality decreases, but the joint inversion consistently outperforms the only PP-wave inversion. Table 2 provides the correlation coefficients and normalized root mean square errors between the inversion results from Figure 7 and the true model. The table shows that the PP-PS joint inversion achieves higher correlation coefficients and lower root mean square errors, indicating that the accuracy of the PP-PS joint inversion is superior to that of the only PP-wave inversion.