Pre-Stack Seismic Inversion With L1-2-Norm Regularization Via A Proximal DC Algorithm And Adaptive Strategy

Seismic inversion in geophysics is an effective way to obtain underground rock properties from seismic survey data on the Earth’s surface. In particular, we can obtain much more information to characterize subsurface geological structure and lithology via pre-stack seismic inversion, with offset information added to the inversion, than by post-stack seismic inversion. However, pre-stack seismic inversion is usually a nonlinear and complicated process. In this article, we adopt a L1-2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L_{1 - 2}$$\end{document}-norm as a constraint on pre-stack seismic inversion, promoting the generation of a sparse solution. We also propose a novel pre-stack seismic inversion method that reduces the complexity of the solving method by utilizing an objective function decomposition scheme. Comparison of calculation time, accuracy and sparsity of the inversion solutions indicates that the proposed algorithm has better accuracy and robustness. Moreover, considering the difficulty of regularization parameter selection, we develop an adaptive parameter selection strategy based on generalized Stein unbiased risk estimation (G-SURE) and incorporate it into the solving algorithm. The adaptive approach finds an appropriate regularization parameter in each iteration and obtains the optimal solution directly, which is beneficial for improving computational efficiency. A synthetic data test verifies that the adaptive method can converge to the optimal solution iteratively in the case of arbitrary initial regularization parameters. Finally, in application to real field data, we explain why the adaptive method is the better choice even though adaptive and non-adaptive methods can obtain solutions with similar accuracy. A novel pre-stack seismic inversion method is proposed based on a proximal difference- of-convex algorithm (pDCA) A new adaptive regularization parameter selection strategy is proposed based on Generalized Stein unbiased risk estimation (G-SURE) Verification that one of the regularization parameters has a limited effect in L1-2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L_{1 - 2}$$\end{document}-norm A novel pre-stack seismic inversion method is proposed based on a proximal difference- of-convex algorithm (pDCA) A new adaptive regularization parameter selection strategy is proposed based on Generalized Stein unbiased risk estimation (G-SURE) Verification that one of the regularization parameters has a limited effect in L1-2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L_{1 - 2}$$\end{document}-norm


Introduction
It is necessary to obtain velocity and density information representing the characteristics of an underground medium based on the pre-stack seismic gather with the help of optimization theory. This process is called seismic inversion, which is essentially an optimization problem in mathematics (Wang 2016). Seismic inversion methods generally include linear and nonlinear methods. Nonlinear inversion methods are easier for solving the global optimal solution because they continuously search in the solution space, but linear inversion significantly reduces the computational cost of nonlinear methods and is more suitable for fast inversion of large-scale seismic data (Maurya et al. 2020). However, seismic inversion is generally an underdetermined problem, making the linear inversion algorithm more dependent on the initial model. To make linear inversion converge at the global extreme value and reduce the instability and multiplicity of the inversion solutions, a regularization constraint is generally used to solve the inverse problem in geophysics (Tarantola 2005).
The sparse constraint is common in signal processing; due to its ability in noise suppression, it has been applied in geophysics for a long time. Claerbout and Muir (1973) proved that L 1 -norm could obtain much better results than L 2 -norm in most cases for seismic numerical modeling. Mainly, with compressed sensing (CS) technology in signal processing, more researchers have focused widely on deconvolution, seismic inversion and seismic data reconstruction (Kazemi and Sacchi 2014;She et al. 2019;Ma 2013). In seismic inversion, Zhang and Castagna (2011) applied the basis pursuit (BP) algorithm to post-stack seismic data to solve the L 1 -norm constrained objective function. This practice shows that the inverted impedance has a higher resolution than the traditional method. Total variation (TV) regularization is a special form of sparse constraint, which is not the L 1 -norm of the parameters calculated directly but the L 1 -norm of the calculated parameters after difference calculation. In this way, sparse edges can be controlled, and the parameters obtained by inversion show characteristics of periodic variation. In seismic inversion, TV regularization can highlight vertical variation characteristics of strata, and the elastic parameters of discontinuous variation can be obtained by inversion (Mozayan et al. 2018). Scholars have applied various methods, such as Iteratively Reweighted Least Squares (IRLS), the split-Bregman and the Alternating Direction Method of Multipliers (ADMM), to retrieve elastic parameters from post-stack data (Zhang et al. 2014;Liu and Yin 2015;Pan et al. 2017). Compared with post-stack data, pre-stack data contain more lithology and fluid information, and a relatively stable solution can also be obtained after extending the L 1 -norm regularization method to the pre-stack AVA inversion (Li and Zhang 2017;Zhi et al. 2016).
L 1−2 -norm is a recently proposed sparse constrained regularization term. It was first addressed by Lou et al. (2015) in the context of nonnegative least squares problems and group sparsity with applications to spectroscopic imaging. Some efforts have been made in an application and solving algorithm which has made L 1−2 -norm gradually applied to geophysics. Wang et al. (2018) make use of the sparsity of L 1−2 -norm to compensate for the attenuation of seismic data and effectively improve signal resolution. Wang et al. (2019) carried out pre-stack seismic inversion under the constraint of L 1−2 -norm regularization and improved the lateral continuity of inversion results by using f-x prediction filtering and obtained a high-quality elastic parameter inversion result. Huang et al. (2021) used L 1−2 -norm regularization to carry out AVA joint inversion based on time domain matching of PP-and PS-waves by using the DTW algorithm, and the algorithm showed good stability.
In fact, a more complex iterative algorithm is needed to solve the sub-problems in each iteration for pre-stack seismic inversion, a large-scale problem, by using L 1−2 -norm, because the traditional difference-of-convex algorithm (DCA) may lead to high cost in large-scale inverse problems (Gotoh et al. 2018). Therefore, in this article, we propose a novel pre-stack seismic inversion scheme by using L 1−2 -norm. The objective function is composed of a misfit function and a constraint term of L 1−2 -norm. Based on the proximal DCA, we develop an optimization algorithm by reformulating the objective function as the difference between two convex functions. In each iteration of DCA, we extrapolate the last solution to obtain the starting point of the new iteration and then use the soft thresholding algorithm to calculate the optimal solution of the current iteration. Moreover, we introduce an adaptive regularization parameter selection method into the new algorithm and propose a strategy to solve the problem that the amplitude of the inversion parameters cannot be well recovered in some cases. To verify the effectiveness of the algorithm and the adaptive parameter selection method, we use synthetic data and actual data to test, respectively. The inversion results verify that the new method is effective for the inverse problem constrained by L 1−2 -norm in the pre-stack seismic inversion, and the adaptive parameter selection method is appropriate.

Pre-Stack Seismic Forward Model
According to the seismic convolution model (Robinson 1967), the pre-stack seismic gather can be represented as reflectivity coefficients for different angles convoluted with the seismic wavelet. A noise term should also be added for real field seismic records. Therefore, for an N-trace pre-stack angle gather, the forward modeling can be expressed as where i , i , i and i are the pre-stack seismic trace data, the seismic wavelet source, the reflectivity coefficient series and the noise at the ith angle of incidence i , respectively. ⟨ * ⟩ represents the convolution operation. In Eq. 1, the reflection coefficient controls the amplitude of the seismic reflection wave and directly reflects the difference of impedance between the upper and lower layers of the reflection interface.
Mathematically, the pre-stack seismic forward modeling can be written as, omitting the noise term for simplicity, where d represents the pre-stack seismic gather, G is the forward modeling operator, and m represents the model parameters vector.

L 1−2 -Norm in Seismic Inversion
The inverse problem of seismic inversion is highly ill-posed which leads to multiple solutions, especially for pre-stack seismic. By adding prior knowledge into the inverse problem, the regularization method in mathematics can solve the problem and obtain a consistent solution with the prior information. Because seismic inversion is mainly for the reflection (2) = coefficient sequence in the time domain, and the reflection coefficient of underground media has obvious sparsity, we can use sparsity and add it into the inversion process. That is the sparse regularization in the optimization of inverse problems in mathematics.
At present, there are three commonly used regularization operators: L 0 -norm, L 1 -norm and L 2 -norm, while L 0 -norm and L 1 -norm are referred to as sparse regularization constraints. As an optimal convex approximation of L 0 -norm, L 1 -norm can guarantee sparsity and it has better solving characteristics. However, the optimization problem containing L 0 -norm belongs to an NP-hard problem, which is challenging to solve. Therefore, L 1 -norm is generally used to construct the objective function, which aims to solve the sparse solution in image processing or geophysical inverse problems (Oldenburg et al. 1983;Yin et al. 2015b;Hamid and Pidlisecky 2015). As a recently emerging sparse constraint method, L 1−2 -norm is getting more attention (Wang et al. 2018Huang et al. 2021).
By comparing the four different regularization terms ( L 0 , L 1−2 , L 1 and L 2 ), the advantage of L 1−2 norm in sparse constraints can be illustrated. The three-dimensional surface comparison of the four different norms ( L 0 , L 1−2 , L 1 and L 2 ) function values is shown in Fig. 1. The function values are calculated by different norm terms using a set of twodimensional data, and the projection of different surfaces on the two-dimensional plane is Fig. 1 Comparison of 3D surface and 2D contour of the values for the different norms. a L 0 -norm. b L 1-2norm. c L 1 -norm. d L 2 -norm a contour map. Under the assumption of the misfit function being unchanged, the process of minimizing the objective function can be regarded as finding the solution of minimizing the regularization term in three-dimensional space. The lowest point of the surface is the optimal point to be searched. After the 3D surface is projected to a 2D plane, it can be found that the lowest point is close to the x-axis and y-axis on the 2D contours. An obvious conclusion is that the closer the norm curve is to the x-axis, the more possibility there is of obtaining the sparse solution by inversion. Therefore, the solution of L 0 -norm regularization is the sparsest among these regularization terms. Although the solution of L 2 -norm regularization is not sparse, it will make the solving process fast and stable. Another important conclusion is that L 1−2 -norm is better than the traditional L 1 -norm in sparsity as a regularization penalty term, and it is more likely to get a sparse solution.
According to the pre-stack seismic forward equation, Eq.
(2), the pre-stack seismic inversion estimates the model parameters m by using the pre-stack seismic data d . Moreover, Eq. (2) represents an underdetermined linear equation. Considering the layered distribution of underground media, a sparse regularization constraint is helpful to obtain sparse reflection coefficient sequences. As in the case of the seismic deconvolution method (Oldenburg et al. 1983), a stable and sparse solution is inverted by using the L 2 -norm misfit function with L 1 -norm regularization In Eq. (3), is a trade-off parameter to balance the first error, or misfit, term and the second term.
By adopting L 1−2 -norm to regularize the pre-stack seismic inverse problem, the constructed objective function includes a general misfit function and a sparse constraint regularization term H( ) is the L 1−2 -norm regularization penalty term, and ∈ (0, 1] is a constant, which can promote the generation of a sparse solution (Lou et al. 2015).

Inversion Algorithm
For the objective function of the pre-stack seismic inversion in Eq. (4), we can use the difference-of-convex algorithm (DCA) to solve. In general, we need to set H( ) in the form of the difference between two scalar norms Therefore, the algorithm then transforms the solution of Eq. (4) into alternate iterations of two variables where k is a subgradient of H 2 k and ⟨⋅⟩ denotes the inner product of two vectors. The key to Eq. (6) is solving k+1 in a new convex optimization problem. A common solution is to introduce additional variables and use the Lagrange multiplier method to construct a form, which can be solved with the help of the ADMM algorithm (Yin et al. 2015a;Wang et al. 2019).
The mathematical derivation has shown that the stability and convergence of DCA depend on the concrete decomposition form of H( ) (Tao and An 1998). To get the closed-form solution of each subproblem, Gotoh et al. (2018) designed a new type of DCA, which is called proximal DCA (pDCA) where L is the Lipschitz constant, and it can be obtained by calculating max svd * . Through this algorithm, each component of DC decomposition will be strongly convex. Equation (4) can be written as By taking the derivative with respect to in Eq. (8) and making it equal to zero, we can solve the extremum problem of the objective function J( ) to obtain the inversion solution in each iteration. The solution of Eq. (8) is given by Wen et al. (2018) In order to accelerate the convergence speed, we can extrapolate a point k in the direction of gradient descent on the line between k−1 and k according to the FISTA (Aster et al. 2012) 2 . FISTA is an upgraded version of ISTA, which is often used to solve the L 1 -norm regularized problem. Let = k+1 , We use a variable k to represent the constant term. Equation (11) can be rewritten as The variable k is the algebraic result of the solution k in the last iteration Taking the derivative with respect to in Eq. (12),

Thus, the solution is given by
Since h k and ∕L are constants in Eq. (15), we can get the solution by the soft thresholding algorithm (Aster et al. 2012). When k > ∕L , if < 0 , then sgn( ) = −1 and = k + ∕L > 0 , and there is a contradiction. But if > 0 , then sgn( ) = 1 and = k − ∕L > 0 , which is the correct solution in accordance with the conditions. In the same way, we can also derive the solution when k < ∕L . In summary, the solution can be expressed as According to the above statements, the procedure for our inversion method is as follows: (1) Set the initial solution 0 , tolerance value and the number of iterations N.

Strata Model
A multilayer geological model ( Fig. 2) is created and used to generate a pre-stack seismic gather dataset to elucidate the robustness, convergence speed and stability of the proposed inversion method. The detailed model parameters are shown in Table 1, including the P-wave velocity V P , S-wave velocity V S and bulk density of each layer. The model consists of 21 layers with different thicknesses, which generates 20 spike reflection coefficients at the interface between two adjacent strata. According to the vertical variation of the model parameters, there will be some small reflection coefficients in the spike sequence (e.g., 2nd, 9th, 11th). The pre-stack reflection coefficient series r(t, ) at different incidence angles are generated by the Aki-Richards approximate equation (Aki and Richards 1980). We assume that the wavelet is a zero-phase Ricker wavelet whose peak frequency is 30 Hz. The synthetic gathers with different incidence angles (5°, 15° and 25°) can be obtained by convolution. Moreover, we add Gaussian random noise with different signal-to-noise ratios (SNR, of 15, 10 and 5, respectively) to verify the robustness of the inversion method in the following section. Noise-free and noisy records are shown in Fig. 3.

Analysis of L 1−2 -Norm Sparse Constraint.
To illustrate the sparsity advantage of the L 1−2 -norm over the L 1 -norm, inversion methods based on different regularization terms are applied to the synthetic data. First, we adopted the FISTA (Pérez et al. 2012) method to do the pre-stack inversion, which utilizes the L 1 -norm regularization term to solve the objective function. Furthermore, our proposed algorithm is implemented, in which L 1−2 -norm regularization is adopted. Both algorithms stop the iteration process under the same convergence condition, that is, Eq. (21). The tolerance value is set to 1 × 10 −5 . The input data are synthetic data with SNR = 5 in the inversion. The inversion results are shown in Fig. 4, including the reflection coefficient inversion results for three angles. The gray bars represent the true reflectivity series, and the black bars are the inversion results.
In FISTA, the regularization parameter is assigned to be 5 × 10 −2 , which was selected by the L-curve. In our proposed method, there are two regularization parameters that need to be determined, of which = 8 × 10 −3 and = 1 × 10 −3 in this example. The iteration numbers of the new algorithm and FISTA are 163 and 59, respectively. For the environment we tested (i7-10700F processor, 32G memory, 64-bit processing system), it took about 3.3520 s and 0.4292 s, respectively. The comparison shows that FISTA has a faster speed when reaching the same convergence condition. At the same time, we carried out error analysis of the different inversion results. The correlation coefficient (CC) is 0.9762, and the normalized root-mean-square error (NRMSe) is 1.46% in FISTA. Correspondingly, our proposed method adopting L 1−2 -norm regularization can significantly reduce the error of the solution, the CC is increased to 0.9971, and the NRMSe is reduced to 0.55% . Since most of the elements in the sparse solution are zero, the difference in error analysis may not be as significant as in Fig. 4. Obviously, the inversion results based on L 1 -norm regularization are not consistent with the amplitude of the real reflection coefficient at some time points. L 1−2 -norm regularization can make up for this deficiency under the condition of sparsity. The analysis of the first numerical experiment shows that the L 1−2 -norm regularization can significantly improve the accuracy of the inversion solution compared with the common L 1 -norm regularization.

Analysis of Inversion Solving Method
In the methodology section, we briefly introduced common methods for solving the L 1−2 -norm regularization, one of which is DCA-ADMM. For the detailed derivation of DCA-ADMM, please refer to Yin et al. (2015b). However, it is quite different from our new Compared with the DCA-ADMM method, our proposed algorithm has a more straightforward solving process in L 1−2 -norm regularization since there is no need to update multiple variables alternately. We also use synthetic data with SNR = 5 to compare the advantages and disadvantages of two different solving methods. The inversion results are shown in Fig. 5. As before, the gray and black bars represent the real reflection coefficient and the inversion result. In contrast to the previous comparison, the solutions based on the same regularization have no noticeable difference in Fig. 5. Combined with the results of error analysis, there is The main difference between the two algorithms is sparsity. The previous section introduced the model with 20 reflection interfaces, so in the case of three incidence angles, the L 0 -norm of the real reflection coefficient is 60, which is called sparsity. By calculating the norm of the solution obtained by the two methods, we can determine the ability of the two algorithms to obtain sparse solutions. The sparsity of the solution obtained by using the proposed algorithm in this paper is 80. However, the sparsity of the solution obtained There are a lot of tiny values in many places that are supposed to be zero. Even if values below 1.9 × 10 −4 are ignored, which is the calculation result of multiplying the minimum real reflection coefficient amplitude 1.9 × 10 −3 by 0.1 , and the L 0 -norm is reduced to 158, the sparsity is still lower than the solution obtained by the new method. This can be observed in Fig. 5a-c, where there are fewer small spikes than in Fig. 5d-f. Therefore, our proposed new method is more suitable for sparse constraint inversion.

Analysis of Noise Effects
Different amounts of noise are generally included in pre-stack synthetic data. It is a crucial problem to ensure the stability of the inversion method without the influence of noise. The SNRs of the synthetic data used for inversion analysis are 5, 10 and 15. The corresponding results are shown in Fig. 6 by using the new method. Regardless of the noise level, a suitable regularization parameter can be determined to make the inversion results match well with the actual model. The CC of the inversion result reaches 0.9996 when the SNR is 15. Even if the SNR decreases to 5 due to a noise increase in the gather, the CC is still as high as 0.9971. The analysis shows that the proposed algorithm is robust for noisy seismic data.

Analysis of Regularization Parameters
The last problem that needs to be considered is the selection of regularization parameter in inversion processing. In our proposed algorithm, there are two regularization parameters and that need to be determined. The effect of is like the regularization parameter in Tikhonov regularization, which is used to adjust the weight of the misfit function and regularization constraint. As a weight parameter within the L 1−2 norm, the influence of is rarely analyzed. Here, we update the sparsity of the model parameters by adding or combining thin layers based on the original multilayer model. The sparsity of the two new models is 30 and 90. By changing in a wide range, we analyze the influence of weight parameter on the inversion results by using three models, which are different in sparsity. We implemented the algorithm on the synthetic records related to these models and calculated the NRMSe of the solutions. Figure 7 shows the variation of NRMSe with regularization parameter . We find that the influence of parameter on the inversion results is closely related to its value. According to the analysis, the error of the solution will decrease rapidly when is less than 1 and then settles into a stable state. Likewise, the trend of the variation curve is like the L-curve, where there is an inflection point. More importantly, the influence of on the solution tends to be stable when it is less than 0.1. Therefore, to reduce the difficulty of solving the inverse problem, we set the parameter to a fixed value of 0.01 in the following test.

Adaptive Selection Method
Since the algorithm only updates the solution in a critical step and does not need to update multiple variables alternately like ADMM, it is possible to use other methods to obtain regularization parameters adaptively. By reasonably changing the regularization parameters to adjust the weight of the regularization term, L 1−2 -norm regularization can effectively solve the ill-posedness of the inverse problem and help us obtain a sparse solution, as shown in the previous section. Generally, the larger the regularization parameter, the greater the dependence of the inversion solution on the initial model, and the less sensitive to small perturbations (Thore 2015). The strategy to select an appropriate regularization parameter quickly and effectively is an important problem to be solved. Common methods include L-curve (Hansen 1992) and GCV (Golub et al. 1979) and have been widely used in seismic inversion (Gholami 2016;Huang et al. 2017). To adjust the regularization parameters adaptively during inversion iteration, we introduce a parameter selection method based on generalized Stein unbiased risk estimation (G-SURE) (Eldar 2008).
In order to find a function F( ) =̂ with noisy data i to minimize the mean square error (MSE), where ̂ represents an arbitrary estimate of , Stein (1981) proposed Stein's unbiased risk estimate (SURE), which is proved to be better than the common maximum likelihood estimation. To further extend the applicability of the SURE to a broad class of problems, Eldar (2008) proposed a generalized SURE by adding a penalty term to the expression. . (23) E We suppose that there are p components greater than ∕L and q components less than − ∕L in the kth iteration. By arranging them together, we can form two vectors with dimensions p and q , marked as 1 and 2 . In the same iteration, the components of ( ) −1 corresponding to the indexes 1 and 2 in v(F) are recorded as 1 and 2 , respectively. Therefore, we can get the expression of each term in Eq. (27) and We can get that the specific form of Eq. (27) under the framework of the new algorithm To obtain the optimal regularization parameter k in k th iteration, we have that We can conclude that It is worth noting that the solution is unstable if (G T G) −1 u is calculated directly in Eq. (31); we can solve the problem by calculating the generalized inverse G † (Infante-Pacheco et al. 2020). Both singular value decomposition (SVD) and truncated singular value decomposition (TSVD) are alternative methods. If the latter method is selected, the GCV method can determine the truncation parameter.

Applicability of L 1−2 -Norm Regularization
To study the applicability of this adaptive strategy in L 1−2 -norm regularization, we use the noisy synthetic data (SNR = 5) in numerical examples as the input data for verification. Before we apply the adaptive strategy to seismic data, a significant problem is the determination of iterative convergence conditions. Once this condition is met, the inversion process is suspended to avoid meaningless computation. In the numerical experiments, we found that the sparsity of the solution will be continuously reduced if Eq. (21) is taken as the convergence condition in the adaptive process, which will cause the components of the solution to disappear. To illustrate this process, we tested the adaptive parameter selection method with three different regularization parameters as initial values independently and recorded the regularization parameters in the adaptive iterative process. Figure 8 shows the recorded results in this process and depicts the variation of solution and regularization parameters. The change of the solution is represented by the L 1−2 norm, which reflects sparsity. The three initial parameters are 1 × 10 −3 , 1 × 10 −5 and 1 × 10 −7 , corresponding to the light gray, black and dark gray lines, respectively, in Fig. 8. During the initial iteration, there is a significant increase in the regularization term. The reason for this change is that the initial solution of the inversion is zero, and the solution is rapidly searching in the direction of descending gradient, while the solution is not sparse. After several times of jitter reduction, the value of the regularization term will decrease sharply, and the regularization parameter will move away from the low-value range and increase significantly in a few successive iterations. In this case, the solution is optimal in the whole iterative process. If the iterative process does not stop, then a following adaptive iterative process will continue to increase the sparsity of the solution until almost all the elements in the inverse solution are zero, which will obviously obtain a solution that is different from the actual model. Equation (21) is not satisfied at the optimal points, which shows that it is not suitable as an iterative convergence condition. Combined with the above analysis results, we choose the following conditions as the convergence condition during iteration The establishment of a new convergence condition combines the rules of parametric variation in Fig. 8. Obviously,Eq. (35) shows that the slope of the regularization term curve is positive, which is opposite to that of the regularization parameter variation curve, and the difference of the regularization parameters between the k and k + 1 iteration is greater than the tolerance .
In the synthetic data test, we set the tolerance to 0.03 and the parameter to 0.01. Figure 9a-c shows the inversion results using these three regularization parameters directly, and Fig. 9d-f shows the inversion results using the adaptive adjustment strategy, where they are only used as initial values. It is worth mentioning that we only compare the situation when the incidence angle is 5° which is like that of other incidence angles. Obviously, a wrong solution can be obtained by directly applying the inversion method in the case of inappropriate regularization parameters. After applying the adaptive strategy, the sparsity of the solution is greatly improved compared with Fig. 9a-c. Nevertheless, the solution is still not optimal because the amplitude of the reflection coefficient has not been well recovered, even if the reflection coefficient appears at the correct time. Moreover, we can adjust the amplitude of the spike to the appropriate value with the help of a hybrid FISTA leastsquares strategy (Pérez et al. 2013). The main idea of the hybrid strategy is to obtain the optimized solution hybrid by rewriting the forward matrix to , which is a dimension reduction form. The optimized solution hybrid is computed by the least squares method where m spike is the rearranged solution by removing the zero element. The results are shown in Fig. 9g-i. After applying the hybrid strategy, the CC of the inversion results is 0.9977, 0.9973 and 0.9967, respectively. Therefore, the application of the hybrid strategy can effectively improve the quality of inversion. Considering that the selection of fixed regularization parameters in the actual data may lead to uncertainty at the non-well location, the adaptive method can be applied to the inversion of real data as an optional method.

Real-Data Application
Finally, we apply the inversion method to real field data to verify the applicability of our proposed new algorithm and adaptive parameters selection strategy. To enhance the SNR of pre-stack gathers, we utilize partially stacked gathers to do the pre-stack inversion from the original gathers. The near-, mid-and far-angles partially stacked gathers are obtained, and the corresponding incidence angles are 8 • , 16 • and 24 • , respectively. The time interval of the target layer ranges from 3100 to 4400 ms on the seismic profile ( Fig. 10a-c), in which there are apparent strong amplitude anomalies and fault development, as it is a typical sandstone reservoir. The frequency range of seismic data is about 4 ~ 80 Hz, and the location of one drilled oil well is at CDP number 292.
Before applying the inversion method, it is necessary to accurately obtain the time-depth relationship (TDR) by well-seismic calibration. The corresponding range of well-logging data on the seismic profile is 3660 to 4100 ms. Three angle-dependent wavelets are estimated from the partially stacked seismic sections. To show the accuracy of the TDR and estimated wavelet, we use the velocity and density logging data in the well to make synthetic records, which were compared with the field data at the well location. The result is shown in Fig. 11, where the x-axis represents incidence angle, and the y-axis represents time. The synthetic data are shown as the black curve and the field data as the red curve. We analyzed the similarity of seismic traces at different incidence angles. The CC is 0.8813, 0.8751 and 0.7300, which represent incidence angles of 8 • , 16 • and 24 • , respectively. The analysis results show that the synthetic records are in good agreement with the field data in the time domain. And the SNR of the far-angle stack data is not as high as the near-angle stack.
In the case of applying the general inversion method to the actual data, we first need to determine the appropriate regularization parameters according to the seismic data at the well location and then extend the parameters to the seismic data volume. The regularization parameter selection method proposed in this paper can adaptively select the appropriate regularization parameters according to the inversion data and the extracted wavelet. It is no longer necessary to extend the regularization parameters determined by synthetic data to the whole dataset.
As an effective parameter to identify reservoir properties, seismic impedance can be obtained via the following equation according to the convolution model where EI is elastic impedance in time domain t, C is an integral matrix, and exp[⋅] denotes exponential operation (Zhang et al. 2014). Here, we adopt Eq. (37) to calculate the seismic impedance from the inverted reflectivity, which is relative impedance. Moreover, for real field data, we can compare it with the result calculated by the logging data. Considering the attenuation in seismic wave propagation, the logging data should be processed by a low-pass filter before calculating relative impedance. At the same time, in the process of calculating seismic impedance by using Eq. (37), if there is an error in the reflection coefficient at a specific time, there will be a low-frequency cumulative error in the whole integration process. To avoid cumulative errors, the usual method is to filter the calculation results, which will cause the results to be band-limited, so the calculation results of logging data need to be processed by the same filter . Figure 12 compares the actual well-logging relative impedance trace and the inverted relative impedance traces by using two inversion methods. The black solid line represents the relative impedance by filtering the well-logging curves. The gray dashed line and gray dotted line represent the inversion results of the adaptive and non-adaptive method, with CC of 0.7585 and 0.7521, respectively. Although the CC between two inversion methods is equivalent, significant errors may occur if the optimal regularization parameters obtained from the trace at the well location are extended to other traces. Through the comparison of inversion results at the well location, the accuracy and effectiveness of our proposed inversion method are verified. Therefore, we inverted the whole seismic section data (shown in Fig. 10). Three angle reflectivity series sections are obtained by using the proposed inversion method and the adaptive parameters selection strategy for real field data processing. Moreover, we can calculate the relative seismic impedance using Eq. 37, which is based on the convolution model of reflectivity coefficient. The inversion results for three-angle stack seismic sections are shown in Fig. 13. Comparing the near-and mid-angles relative impedance sections, at 3700 ms, there is a bright spot, and its impedance is significantly higher than that of the surrounding rocks, so it can be considered as a favorable exploration target. The relative impedance section of the far-angle stacked seismic section (Fig. 13c) is not consistent with the inversion results of the near-and midangles stacked seismic sections, but it is consistent with the pre-stack seismic gathers as shown in Fig. 11 and the stacked seismic section as shown in Fig. 10c.

Conclusions
We have combined the sparse regularization of L 1−2 -norm and the proximal difference-ofconvex algorithm (pDCA) to implement a pre-stack seismic inversion technique to obtain reflectivity coefficients from pre-stack seismic gather data. The analysis of synthetic data inversion results verified that L 1−2 -norm is better than the traditional L 1 -norm in sparsity as a regularization penalty term. At the same time, the analysis of the inversion solving method and noise effects indicated that our implemented inversion method by using the (37) EI(t) = EI t 0 * exp [2 * C * r(t)], Fig. 9 Inversion results of three strategies for regularization parameter. a-c Fixed regularization parameter strategy, d-f adaptive method and g-i hybrid strategy. The results are estimated with different initial regularization parameters: (a, d, g) 1 × 10 −3 , (b, e, h) 1 × 10 −5 , and (c, f, i) 1 × 10 −7 . The gray and black bars correspond to the actual reflectivity series and the inversion results, respectively 1 3 pDCA is also more suitable for sparse constraint inversion than the DCA-ADMM, which is a valuable inversion algorithm used in seismic inversion.
We presented an adaptive hybrid strategy, discussing the influence of regularization parameter in L 1−2 -norm and deriving the adaptive selection method of regularization parameter based on G-SURE, which could reduce the computation time to obtain the appropriate parameter. The analysis of different sparsity inversions of synthetic data inversion indicated that our adaptive parameter section method could effectively improve the quality of inversion results. Considering that the selection of fixed regularization parameters in real data may lead to uncertainty at the non-well location, the adaptive method can be applied to the inversion of real data as an optional method.
If only one optimization parameter is considered, the computing times of adaptive and non-adaptive methods are the same. However, the limitation of the non-adaptive method is that it needs many tests to obtain a stable regularization parameter. Therefore, Fig. 10 Three partial stacked seismic sections. a Near-angle stacked section, b mid-angle stacked section, and c far-angle stacked section the non-adaptive method is more time-consuming in obtaining a solution with the same precision. The adaptive method simplifies the process of parameter selection, which only requires an initial value to update the regularization parameters in the iterative process. Nevertheless, the adaptive method based on G-SURE needs to use the generalized inverse operator G † to obtain a low-precision initial solution (G T G) −1 u . Each component of the initial solution will be tested to determine whether it is appropriate in the iterative process. There is no doubt that the choice of the initial solution will influence the effect of the adaptive method. In this paper, we use truncated singular value decomposition (TSVD) to obtain a generalized inverse, which has an obvious effect and avoids the introduction of more regularization parameters.
We presented an application to real field data by using the proposed pre-stack seismic inversion and adaptive hybrid strategy for selection of regularization parameters. The results showed that our inversion method is effective and stable for real field seismic Comparison of field data (red) and synthetic data (black) at the well location in the time domain Fig. 12 Comparison of the actual well-logging relative impedance trace and the inverted relative impedance traces by using two inversion methods in the time domain. The black line is relative impedance by filtering from the well-logging curves, the gray dashed line represents the result by the adaptive method, and the gray dotted line represents the result by the non-adaptive method data, which is generally contaminated by ambient noise. In this paper, we calculated the relative impedance from the inverted reflectivity coefficient series in the time domain according to the convolution model. The relative impedance is band-limited because the seismic data are band-limited, and filtering is required to eliminate the low frequency accumulated error in the integration process. More importantly, the elastic parameters cannot be directly inverted from the relative impedance; this requires wide-band impedance data based on a low-frequency model, but this is not the focus of this paper. Considering the wider application of sparse constraints in geophysics as reviewed in the introduction, our proposed inversion method based on L 1−2 -norm can be applied to   where Tr[⋅] denotes the trace of a matrix.