Next Article in Journal
Li2O-Based Cathode Additives Enabling Prelithiation of Si Anodes
Next Article in Special Issue
Productivity Prediction of Fractured Horizontal Well in Shale Gas Reservoirs with Machine Learning Algorithms
Previous Article in Journal
Influence of Double-Limb Double-Plate Connection on Stable Bearing Capacity of Quadrilateral Transmission Tower
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Prestack Seismic Inversion via Nonconvex L1-2 Regularization

1
School of Electronic and Information Engineering, Chongqing Three Gorges University, Chongqing 404000, China
2
State Key Laboratory of Oil and Gas Reservoir Geology and Exploitation, Chengdu University of Technology, Chengdu 610059, China
*
Author to whom correspondence should be addressed.
Appl. Sci. 2021, 11(24), 12015; https://doi.org/10.3390/app112412015
Submission received: 25 October 2021 / Revised: 8 December 2021 / Accepted: 14 December 2021 / Published: 17 December 2021
(This article belongs to the Special Issue Geomechanics and Reservoirs: Modeling and Simulation)

Abstract

:
Using seismic data, logging information, geological interpretation data, and petrophysical data, it is possible to estimate the stratigraphic texture and elastic parameters of a study area via a seismic inversion. As such, a seismic inversion is an indispensable tool in the field of oil and gas exploration and development. However, due to unknown natural factors, seismic inversions are often ill-conditioned problems. One way to work around this unknowable information is to determine the solution to the seismic inversion using regularization methods after adding further a priori constraints. In this study, the nonconvex L 1 2 regularization method is innovatively applied to the three-parameter prestack amplitude variation angle (AVA) inversion. A forward model is first derived based on the Fatti approximate formula and then low-frequency models for P impedance, S impedance, and density are established using logging and horizon data. In the Bayesian inversion framework, we derive the objective function of the prestack AVA inversion. To further improve the accuracy and stability of the inversion results, we remove the correlations between the elastic parameters that act as initial constraints in the inversion. Then, the objective function is solved by the nonconvex L 1 2 regularization method. Finally, we validate our inversion method by applying it to synthetic and observational data sets. The results show that our nonconvex L 1 2 regularization seismic inversion method yields results that are highly accurate, laterally continuous, and can be used to identify and locate reservoir formation boundaries. Overall, our method will be a useful tool in future work focused on predicting the location of reservoirs.

1. Introduction

Prestack amplitude variation angle (AVA) inversion methods can be used to estimate multiple elastic parameters in a given study area by including observational logging and petrophysical data [1]. Those methods are of great interest to those in the field of reservoir fluid identification and evaluation. There are two main types of prestack seismic inversions: prestack direct inversion and elastic impedance inversion.
Prestack direct inversion, also known as AVA/AVO inversion, is based on the P-wave approximate reflectivity formulas in the Zoeppritz equation and uses angle gather data as constraints in the prestack elastic parameter inversion [2,3]. For example, Downton improved the accuracy of the AVO inversion results by using the covariance matrix to remove the correlation between the initial elastic parameter constraints [4]. To improve the identification of thin interbeds, Zhang et al. derived a detailed multi-angle prestack forward model by extending the theory of poststack reflectivity decomposition(the dipole decomposition principle)to the inversion of prestack P/S velocities and densities [5,6]. However, due to the lack of far-angle gather seismic data, the density parameters were used as the initial inversion constraints in those specific application examples. To highlight the “ b l o c k y ” inversion effect, Pérez et al. performed a three-parameter prestack inversion using a weighted L 1 , 2 mixed norm [7]. In this inversion method, both the correlation of the inversion parameters and the correlation of the model parameters are considered; using the fast iterative shrinkage-thresholding algorithm (FISTA), they were able to solve the inversion objective function and achieve high-resolution inversion results [8]. In a later study, Li et al. quantified the Young’s modulus, Poisson’s ratio, and density values using a combination of L 1 , L 2 , and total variation (TV) norm regularization methods [9]. Of these three types of methods, L 1 norm regularization reflects the sparsity of the inversion results, L 2 norm regularization forces the inversion results to conform to the initial constraints, and TV norm regularization enhances the lateral continuity of the inversion results. To eliminate the P/S velocity ratio as a parameter in the traditional approximate reflectivity equation, Fu et al. first re-derived the reflectivity equation to define the P/S impedance and then developed the P/S impedance inversion method [10]. Using high-order TV regularization, She et al. ran multiple three-parameter prestack inversions that accounted for various topographical characteristics in multiple study areas [11]. Those methods are sensitive to noise and the amount and quality of the observational data used in the inversion. To improve the stability of the inversion method, Nie et al. proposed a novel P/S impedance inversion used in conjunction with L 1 2 regularization, and they demonstrated the applicability and feasibility of their method, however, which method also faces difficulties in producing accurate density values [12].
The other types of prestack seismic inversion methods are elastic impedance inversion. Those methods are undertaken in two parts. In the first part, the elastic impedance is inverted from partially angularly stacked seismic data. In the second part, the to-be-inverted prestack elastic parameters are extracted from the elastic impedance. The concept of elastic impedance was introduced by Connolly and later extended by Whitcombe to provide a theoretical basis for prestack inversion [13,14]. Zong et al. derived the Zoeppritz approximate equation for P/S-wave moduli and performed an inversion to determine the prestack P/S-wave moduli [15]. In a later study, Zong et al. derived the approximate reflectivity equations for Young’s modulus, Poisson’s ratio, and density based on the assumption of plane wave incidence, and obtained Young’s modulus and Poisson’s ratio, further enriching the prestack inversion method [16]. Liu et al. performed an inversion to quantify the Gassmann fluid term and the shear modulus by applying the basis pursuit decomposition technique and the theory of dipole reflectivity decomposition to a prestack inversion. Using this inversion method, they also conducted a sensitivity analysis to determine how the elastic parameters changed in the presence of various reservoir fluids [17]. Furthermore, they made significant improvements to the prestack AVO inversion method; by refining the Zoeppritz approximate equation, they located and characterized hydrocarbon-bearing subsurface reservoirs using an inversion that was constrained by various elastic parameters such as the Young’s modulus, Poisson’s ratio, Gassmann fluid term, and brittleness index values [18,19,20]. Zhang et al. first inverted the reflectivity of seismic data with different angles and then extracted the brittleness index from those reflectivity values [21]. Wang et al. used nonconvex L 1 2 regularization to extract stable elastic impedance values from seismic data and then employed a Bayesian inversion framework to back out the P/S velocities and the densities from the elastic impedance values [22]. However, this method requires the use of logging data to establish the relationship between the elastic impedance and the inverted parameters. Therefore, the reliability, accuracy, and precision of this inversion method must be determined in future studies.
Meanwhile, the high dimensionality of multi-angle seismic data increases the ill-posedness and decreases the stability of a prestack AVA inversion. Typically, this issue is resolved using regularization methods such as the Tikhonov and many variations of the L 1 regularization methods. However the Tikhonov regularization method often results in over-smoothed solutions, it is not an ideal technique to use with applications related to identifying and characterizing formation boundaries. More importantly, some existing AVO inversion methods that use L 1 regularization techniques fail to generate the optimal sparse solutions because they exacerbate the existing prestack AVA high dimensionality issue by requiring additional model constraints.
With the potential to generate inversion results that are both optimally sparse and highly accurate, nonconvex optimization algorithms have become quite popular in recent years. Lai et al. proposed an L p -norm minimization algorithm that is solved by iteratively reweighting the least-squares solution [23]. Yin et al. proposed L 1 /L 2 and L 1 2 minimization algorithms and confirmed the superiority of nonconvex optimization algorithms in terms of solution accuracy and efficiency using both theoretical examples and actual case studies [24]. Because of its good performance in promoting sparsity, nonconvex optimization algorithms have been widely applied to geophysical problems. Zhong et al. proposed a L 1 / 2 norm regularization to reconstruct highly incomplete seismic data and obtained obviously better results [25]. Recently, Huang et al. used L 1 2 norm regularized logarithmic absolute misfit function to improve the stability and fidelity of the prestack seismic inversion [26]. In the same year, they combined with the logarithmic absolute-criterion-based misfit function with L 1 2 norm-based penalty for P-P and P-SV waves joint inversion and achieved good performance in both resolution and accuracy [27].
With this motivation, we employ nonconvex L 1 2 regularization to simultaneous P/S-impedance and density inversion from pre-stack seismic data. In this study, we derive the multi-angle, multi-reflection interface prestack AVA forward model using the approximate reflectivity equation given by Fatti [2]. We assemble the covariance matrix to remove the correlation between the parameters that act as constraints in our method. We then use low-frequency model to determine the prestack inversion objective function in a Bayesian inversion framework and complete the inversion by applying a nonconvex L 1 2 regularization algorithm. In our validation experiments, we not only quantify the performance of our method in synthetic tests and in field data test, but we also explore how our method addresses existing issues such as prestack AVA inversion instability and the difficulties that arise in inverting for the density.

2. Method

2.1. The Forward Problem

Fatti’s three-term formula, which approximates the relationship between the seismic reflectivity of the P/S impedance and the density for a given data set, is expressed as [2]:
r ( θ ) = A ( θ ) Δ I p I ¯ p + B ( θ ) Δ I s I ¯ s + C ( θ ) Δ ρ ρ ¯
where, A ( θ ) = 1 2 ( 1 + tan 2 θ ) , B ( θ ) = 4 V s 2 V p 2 sin 2 θ , C ( θ ) = 1 2 ( tan 2 θ 4 V s 2 V p 2 sin 2 θ ) , I ¯ p = I p t + 1 + I p t 2 , I ¯ s = I s t + 1 + I s t 2 , ρ ¯ = ρ t + 1 + ρ t 2 , Δ I p = I p t + 1 I p t , Δ I s = I s t + 1 I s t , Δ ρ = ρ t + 1 ρ t , where t is sampling time point, θ is the angle of incidence, I p t , I s t , ρ t are the P impedance, S impedance, and density at t, respectively. The matrix form of Equation (1) is
r ( θ ) = [ A ( θ ) B ( θ ) C ( θ ) ] [ Δ I p I ¯ p Δ I s I ¯ s Δ ρ ρ ¯ ]
In the case of multiple sampling points, Equation (2) is extended as follows:
[ r 1 ( θ ) r 2 ( θ ) r t ( θ ) ] = [ A ( θ ) B ( θ ) C ( θ ) A ( θ ) B ( θ ) C ( θ ) A ( θ ) B ( θ ) C ( θ ) ] × [ r p r s r ρ ]
where, r p = Δ I p I ¯ p , r s = Δ I s I ¯ s , r ρ = Δ ρ ρ ¯ . With multiple sampling points and N incident angles, the matrix form of the approximate reflectivity equation is defined as:
[ R ( θ 1 ) R ( θ 2 ) R ( θ N ) ] R = [ C p ( θ 1 ) C s ( θ 1 ) C ρ ( θ 1 ) C p ( θ 2 ) C s ( θ 2 ) C ρ ( θ 2 ) C p ( θ N ) C s ( θ N ) C ρ ( θ N ) ] C × [ r p r s r ρ ] r
where R θ N = r 1 θ N r 2 θ N r t θ N T , r p , r s , and r ρ are the reflectivity series of the P impedance, S impedance, and density, respectively. C p θ N , C s θ N , and C ρ θ N are the diagonal coefficient matrices composed of A ( θ ) , B ( θ ) , and C ( θ ) , respectively; for brevity, we have not included the forms of these matrices.
According to the seismic convolution model, the relationship between the prestack angle gather seismic data and the reflectivity is expressed as
s θ 1 s θ 2 s θ N s w θ 1 w θ 2 w θ N W × R θ 1 R θ 2 R θ N R
where w θ N is the wavelet matrix corresponding to the N th angle seismic data. Equations (4) and (5) can be written concisely in matrix form as
s = WCr
where s = s θ 1 s θ 2 s θ N T is a multi-angle seismic data series, W is a multi-angle wavelet matrix with the various w θ N matrices on its diagonal, and C is a coefficient matrix comprised of C p ( θ ) , C s ( θ ) , and C ρ ( θ ) .

2.2. Establishment of Low-Frequency Model Constraints

To enhance the lateral continuity and stability of the inversion results, our low-frequency P impedance model must satisfy the following equation:
LBr p = ξ p l o w
where ξ p l o w is the low-frequency P impedance model determined using the logging data, and L is the low-frequency filter matrix. B is the integral operator. Similarly, the low-frequency S impedance and density models must satisfy Equations (8) and (9), respectively:
LBr s = ξ s l o w
LBr ρ = ξ ρ l o w
where ξ s l o w and ξ ρ l o w are the low-frequency S impedance and density models established from the logging data, respectively.
In augmented matrix form, these three constraints are expressed as follows:
L B 0 0 0 L B 0 0 0 L B Φ r p r s r ρ r = ξ p l o w ξ s l o w ξ ρ l o w ξ
where the simplified matrix form of this expression is
Φ r = ξ

2.3. Construction of the Inversion Objective Function

According to Bayesian inversion theory and the definition of the L 1 2 regularized prior distribution function, the inversion objective function that simultaneously satisfies the P impedance, S impedance, and density constraints is expressed as (see Appendix A for details):
m ^ = arg min 1 2 W C r s 2 2 + β 2 Φ r ξ 2 2 + λ r 1 r 2
where β is the weight coefficient for the low-frequency P impedance, S impedance, and density models and λ is the regularization tuning parameter that controls the sparsity of the parameters that will be used as input values in the inversion.
If we define G = W C β Φ and b = S β ξ , then Equation (12) can be further simplified into a nonconvex L 1 2 regularization criterion:
r ^ = arg min 1 2 G r b 2 2 + λ r 1 r 2

2.4. Decorrelation Processing of Inversion Parameters

According to prior work in this field and statistical analyses of logging curve data, it has been determined that the P impedance, S impedance, and density values that are input into the above forward model are not independent of one another. An additional problem arises because these three parameters have different dimensions. To improve the accuracy of our inversion results, we decorrelate the P impedance, S impedance, and density parameters using singular value decomposition (SVD).
The relative changes in the parameters to be inverted are expressed in matrix form as
X = e p 1 e p 2 e p N e s 1 e s 2 e s N e ρ 1 e ρ 2 e ρ N
where N is the number of sampling points and e p , e s , and e ρ are the natural logarithmic series of the P impedance, S impedance, and density values obtained from the logging curve, respectively. Based on Equation (18), we define the symmetric covariance matrix as:
C x = X X T N 1 = c 11 c 12 c 13 c 21 c 22 c 23 c 31 c 32 c 33
The SVD of this matrix is
C x = vuv 1 = vuv T
Substituting Equation (15) into Equation (16) and rearranging the variables yields:
u = v T C x v = σ p 2 σ s 2 σ ρ 2 = v T X X T v N 1
where v = v 11 v 12 v 13 v 21 v 22 v 23 v 31 v 32 v 33
For a single-layer interface, u = σ p 2 σ s 2 σ ρ 2 can be used as a priori constraints. However, with a multi-layer interface, we first obtain r via decorrelation during the actual inversion process. Next, we convert r into r :
r = V ˜ r
where, V ˜ = V 11 V 12 V 13 V 21 V 22 V 23 V 31 V 32 V 33 and V ij = diag v i j 1 v i j 2 , v i j N , for i,j = 1,2,3.
Once the decorrelation is completed, we can redefine the inversion objective function shown in Equation (19):
x = arg min 1 2 G x b 2 2 + λ x 1 x 2
where G = G V ˜ , x = r .

2.5. Resolving the Inversion Objective Function

To find the solution of the above objective function, we apply the convex function difference algorithm in order to decompose Equation (19) into F ( x ) = G ( x ) H ( x ) , where
G ( x ) = 1 2 G x b 2 2 + λ x 1 H ( x ) = λ x 2
When x 0 , x 2 can be expressed as x x 2 ; when x = 0 , 0 x 2 . We then expand Equation (24) into its iterative solution form:
x n + 1 = arg min 1 2 G x b 2 2 + λ x 1 , x n = 0 arg min 1 2 G x b 2 2 + λ x 1 < x , λ x n x n 2 > , x n 0
According to the DCA iterative formula [28], the nonconvex L 1 2 regularized inversion objective function can be simplified as
x ^ = arg min 1 2 G x b 2 2 + λ x 1 + v , x
where v , which is a constant vector defined as v = λ x n x n 2 , x 0 , is the gradient of H ( x ) at X n .
After introducing an intermediate dummy variable z , we can restate the inversion objective function in Equation (23) as
x n + 1 = arg min x R n 1 2 G x b 2 2 + v , x + λ z 1 s . t . x z = 0
Using the augmented Lagrange multiplier method, the final inversion objective function is
L δ ( x , z , y ) = 1 2 G x b 2 2 + v , x + λ z 1 + y T ( x z ) + δ 2 x z 2 2
where y is the Lagrange multiplier and σ is the penalty parameter. Equation (24) is solved by employing the alternating direction method of multipliers [23]; this method has an iterative recursive form that is expressed as
z k + 1 = arg min z L δ x k , z , y k x k + 1 = arg min x L δ x , z k + 1 , y k y k + 1 = y k + δ x k + 1 z k + 1
where k is the number of iterations. The iterative formulas for x and z are
x k + 1 = G T G + δ I 1 G T b v + δ z k y k z k + 1 = shrin k x k + y k / δ , λ / δ
where shrink ( · ) is the soft threshold operator. To ensure fast convergence of the algorithm, we apply an adaptive penalty factor δ k + 1 :
δ k + 1 = min δ max , τ δ k ,
where min ( · ) is the minimum value function and δ max is the maximum value of δ k . We define the constraint weight τ as
τ = τ 0 , δ z k + 1 z k 2 x k 2 < ε , 1 , δ z k + 1 z k 2 x k 2 ε .
where τ 0 and ε are the initial values for these parameters that are set to 0.001.
Once we’ve obtained r from r in Equation (18), we can calculate the corresponding P impedance, S impedance, and density values:
I p = ξ p l o w ( 1 ) exp 2 i = 1 N r p ( i ) , I s = ξ s l o w ( 1 ) exp 2 i = 1 N r s ( i ) , ρ = ξ ρ l o w ( 1 ) exp 2 i = 1 N r ρ ( i ) .

3. Synthetic Examples

To verify the effectiveness of our proposed inversion methodology (hereafter, L 1 2 ), we compare the results of our technique against those of the conventional basis pursuit method (BP) and the prestack AVA inversion method [6]. For single-channel experimental analysis, our model test data consists of the observed logging curve for a given work area, the related anti-noise experimental analysis, and the sensitivity of those inversion results to the angle gathers. First, the logging data are resampled and segmented at a sampling interval of 2 ms to obtain the P velocity curve (red), the S velocity curve (blue), and the density curve (green) (Figure 1); the black curve represents the low-frequency models of those same three parameters. By inputting a zero-phase Ricker wavelet with a dominant frequency of 35 Hz into the Zoeppritz equation forward model, we generate a total of 15 prestack angle gathers with an angle range of 3 to 45 at an interval of 3 (Figure 1d) [29]. The synthetic angle gathers record has significant AVO characteristics.

3.1. Noise-Free Test

The results of the noise-free analysis of the synthetic record shown in Figure 1d are shown in Figure 2. In addition to generating more accurate elastic parameter values than the BP method, the L 1 2 inversion method also produces blocky inversion results that facilitate the identification and characterization of the relevant formation boundaries. As shown in Figure 2b, the BP inversion method produces an overly smoothed solution that is not sparse. In contrast, the improved low-frequency model constraint in our proposed method uses the low-frequency component of the inversion results to approximate the initial low-frequency model. This step not only improves the accuracy of the inversion results but also reduces the sensitivity of the inversion results to the initial low-frequency model.

3.2. Noise Test

Downton pointed out that the prestack AVA inversion method is highly susceptible to noise [4]. As such, it is very important to determine how our L 1 2 inversion method performs in an anti-noise experiment. Multiple synthetic records with signal-to-noise ratios (SNRs) of 10 dB, 5 dB, and 1 dB are shown in Figure 3; the inversion results of these three sets of synthetic records are shown in Figure 4, Figure 5 and Figure 6, respectively. While the results of the L 1 2 inversion method are comparable to those shown in Figure 2, the BP inversion results have become increasingly unstable in the presence of noise, leading to a large error in the density model and distinct inaccuracies in the P/S impedance model. While the overall accuracy of the L 1 2 inversion method results decreases as the SNR of the synthetic record decreases (Figure 5), the L 1 2 inversion results are still more accurate than those of the BP inversion method in terms of the P/S impedance values. At even lower SNR values (Figure 6), the density results from neither the L 1 2 or BP inversion methods are credible; however, the L 1 2 method P/S impedance results are still valid. In summary, while these inversion methods are sensitive to the SNR of the seismic records, the L 1 2 method still produces accurate P/S impedance values, even when using seismic data with a low SNR value.

3.3. Analysis of Factors Influencing the Inversion Results

Downton and Zong et al. both argued that high-quality far-angle gathers are the key to generating accurate density results when performing a prestack AVO inversion [4,15,16]. To study this problem in-depth, we divided the angle gathers ( 3 to 45 ) into five groups of small, medium, and large angles (Table 1). The inversion results corresponding to the five sets of angle gathers are shown in Figure 7, Figure 8, Figure 9, Figure 10 and Figure 11.
Comparison of the five groups of experimental results shows that the P/S impedance results for the L 1 2 and BP inversion methods are not sensitive to the angle gathers. However, per the Fatti approximation formula, errors in the density calculations decrease as the degree of the angle gather increases. Even in the absence of far-angle gather seismic data, the P/S impedance and S impedance results generated by the L 1 2 inversion method are more accurate than those produced by the BP inversion method.

4. Application

In this case study, we used actual seismic data fromsampling the Pearl River Mouth Basin; this data consists of high-quality far-angle gathers suitable for a prestack AVA inversion study. The reservoir formation is comprised of sandstone, while the upper and lower layers surrounding the reservoir are mudstone. Previous work in this area indicates that this reservoir formation, known as the Pearl River Formation, is characterized by shallow deltaic deposits. This formation is deeply buried, structurally complex, and exhibits significant lateral variation [12]. Based on the results of our model experiments, we separated the seismic data into three angle ranges: 3 15 , 18 30 , and 33 45 . We then used this data to generate seismic data sets with center angles of 9 , 24 , and 39 , respectively (Figure 12a–c, respectively).
We then perform prestack AVA inversions with these three data volumes. To enhance the lateral continuity of the inversion results, we use the logging data and horizon information to establish the low-frequency model required for the inversion (Figure 13a–c).
Employing our L 1 2 inversion method, we performed the prestack AVA inversion and obtain three elastic parameters: P impedance, S impedance, and density (Figure 14a–c, respectively). In these figures, the black curves represent the logging curves corresponding to the prestack elastic parameters. The results in Figure 14 indicate that there are low P impedance and density anomalies and high S impedance anomalies. Because shear waves do not propagate in fluid media, the S-wave propagation velocity is comparable to that of the sandstone skeleton, which is in line with the actual geological interpretation of this study area. As shown in Figure 14d, our calculated P/S velocity ratio results not only allows us to pinpoint the location of our target formation, but they are also highly consistent with the P/S velocity ratios determined using the logging data (black curve in Figure 14d). As we have demonstrated, it is feasible to apply our L 1 2 prestack AVA inversion method to the task of reservoir prediction.

5. Discussion

To further explore our inversion results and to identify the factor(s) affecting the density inversion results, we plugged the above five groups of inversion results (Figure 7, Figure 8, Figure 9, Figure 10 and Figure 11) into Equation (10) to generate synthetic seismic records. We then compared these five sets of synthetic angle gathers with the original angle gathers and calculated the difference between the synthetic and observed data gathers (Figure 15, Figure 16, Figure 17, Figure 18 and Figure 19). Based on this analysis, we found that the synthetic angle gathers are highly consistent with the original data. Furthermore, we have determined that accurate synthetic records can only be generated with the P/S impedance inversion results; as such, it is difficult to calculate the density values from an AVA inversion.

6. Conclusions

In this study, we developed a new method for performing a prestack AVA inversion that relies on nonconvex L 1 2 regularization. In our method, both the construction of the inversion objective function and the decorrelation of the input parameters improve the stability of the inversion results. Additionally, we verified the effectiveness of our new method via synthetic tests and conducting a specific case study with real data. Our L 1 2 prestack AVA inversion method produces P/S impedance results that are more accurate and stable than those generated using the BP inversion method. From our P/S impedance results, we calculated P/S velocity ratios that allowed us to easily identify a target reservoir formation. Although prestack AVA inversion methods are sensitive to noise and the amount and quality of the observational data used in the inversion, our proposed method demonstrates better performance and has a high potential in the detection and identification of fluids.

Author Contributions

Project administration, X.W.; writing—original draft preparation, W.N.; writing—review and editing, F.X. and B.L.; Supervision X.W. and X.N. All authors have read and agreed to the published version of the manuscript.

Funding

This study was financially supported by Chongqing Key Laboratory of Geological Environment Monitoring and Disaster Early-warning in Three Gorges Reservoir Area (No. ZD2020A0501), Open Fund (PLC20211105) of State Key Laboratory of Oil and Gas Reservoir Geology and Exploitation (Chengdu University of Technology), youth project of science and technology research program of Chongqing Education Commission of China(Grant No.KJQN202101226),the National Science Foundation of China (Grant No. 41774142), and the Major National Science and Technology Projects of China (Grant No. 2016ZX05002-004-013).

Institutional Review Board Statement

Not Applicable.

Informed Consent Statement

Not Applicable.

Data Availability Statement

The data presented in this study are available on request from the corresponding author. The data are not publicly available due to privacy.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A. Bayesian Derivation of the Objective Function

According to Bayes’ theorem, the posterior probability distribution (PDF) of the estimated parameter r may be expressed as
p ( r | s ) = p ( s | r ) p ( r ) p ( s )
where p ( r | s ) is the PDF of r . The likelihood function, p ( r | s ) , summarizes the strength of the fit between the observed data and r . p ( r ) represents the prior probability distribution of r in the absence of observed data s . p ( s ) is the probability of the entire sample space, which is typically a constant. Therefore, Equation (A1) is equivalent to
p ( r | s ) p ( s | r ) p ( r )
In seismic exploration, the observed data are seismic data. Prior to performing the inversion, the seismic data must undergo preprocessing procedures such as dynamic correction, multiple wave suppression, and attenuation compensation. It is typically assumed that the likelihood function has a Gaussian distribution with a mean of zero and a variance of c d 2 . Therefore, p ( s | r ) may be expressed as
p ( s | r ) exp 1 2 ( W C r s ) T Σ 1 ( W C r s )
where Σ is the covariance matrix that represents noise in the seismic data ( Σ = c d 2 · I ) and I is the identity matrix.
Certain types of a priori information, such as structural and well-log data, significantly improve the stability of the inversion by providing useful constraints for the model results. According to [12,30], there are two types of useful a priori information: (1) low-frequency data that adhere to a (0, c l o w 2 ) Gaussian distribution and (2) constraints on the sparsity of the estimated parameter. The low-frequency information, which is independent of the seismic data and is typically obtained by interpolating and extrapolating the horizon and well-log data, is used to enhance the lateral continuity of the inversion. The sparsity information is used to characterize the strata boundaries. As such, p ( r ) may be defined as
p ( r ) = p L ( r ) p V ( r )
where p L ( r ) is the lateral constraint for the low-frequency information and p V ( r ) is the vertical sparsity constraint for the estimated parameter. The definition of p L ( r ) is
p L ( r ) exp 1 2 L B r ξ l o w T Ω 1 LBr ξ l o w
where L is the low-pass matrix, ξ l o w is the low-frequency prior information ( ξ l o w ( t ) = 1 2 ln I l o w ( t ) I l o w t 0 ), I l o w is the low-frequency P impedance trend that was built from the well-log data, Ω is the covariance matrix of the model parameters ( Ω = c l o w 2 · I ).
The L 1 2 regularized prior distribution function is defined as
p ( r ) = 1 2 π 1 2 exp 1 2 r 1 r 2
which is equivalent to
p ( r ) exp 1 2 r 1 r 2
By substituting the likelihood function (Equation (A3)) and prior information equations (Equations (A5) and (A7)) into the Bayes theorem equation (Equation (A2)), we can obtain the final expression for p ( r | s ) :
p ( r | s ) exp 1 2 ( WCr s ) T Σ 1 ( WCr s ) × exp 1 2 ( LBr ξ l o w ) T Ω 1 ( LBr ξ l o w ) × exp 1 2 r 1 r 2
Simplifying Equation (A8) and taking the natural logarithm on both sides yields
J ( r ) = ln [ p ( r s ) ] 1 2 W C r s 2 2 + β 2 L B r ξ l o w 2 2 + λ r 1 r 2
where λ and β are regularization parameters, which are used to tune the inversion sparsity and low-frequency model constraints. These parameters are defined as λ = c d 2 c y 2 and β = c d 2 c l o w 2 . Maximizing the posterior probability results in the objective function for AVA inversion:
r ^ = arg min 1 2 W C r s 2 2 + λ r 1 r 2 + β 2 LBr ξ low 2 2

References

  1. Russell, B.H. Prestack seismic amplitude analysis: An integrated overview. Interpretation 2014, 2, SC19–SC36. [Google Scholar] [CrossRef] [Green Version]
  2. Fatti, J.L.; Smith, G.C.; Vail, P.J.; Strauss, P.J.; Levitt, P.R. Detection of gas in sandstone reservoirs using AVO analysis: A 3-D seismic case history using the Geostack technique. Geophysics 1994, 59, 1362–1376. [Google Scholar] [CrossRef]
  3. Gidlow, P.M.; Smith, G.C.; Vail, P. Hydrocarbon detection using fluid factor traces: A case history. In 3rd SAGA Biennial Conference and Exhibition, Expanded Abstracts; European Association of Geoscientists and Engineers: Houten, The Netherlands, 1993; pp. 78–89. [Google Scholar]
  4. Downton, J.E. Seismic Parameter Estimation from AVO Inversion. Ph.D. Thesis, The University of Calgary, Calgary, AB, Canada, 2005. [Google Scholar]
  5. Zhang, R.; Castagna, J. Seismic sparse-layer reflectivity inversion using basis pursuit decomposition. Geophysics 2011, 76, R147–R158. [Google Scholar] [CrossRef] [Green Version]
  6. Zhang, R.; Sen, M.K.; Srinivasan, S. Multi-trace basis pursuit inversion with spatial regularization. J. Geophys. Eng. 2013, 10, 035012. [Google Scholar] [CrossRef]
  7. Pérez, D.O.; Velis, D.R.; Sacchi, M.D. Three-term inversion of prestack seismic data using a weighted L2,1 mixed norm. Geophys. Prospect. 2017, 65, 1477–1495. [Google Scholar] [CrossRef]
  8. Beck, A.; Teboulle, M. A fast iterative shrinkage-thresholding algorithm for linear inverse problems. SIAM J. Imaging Sci. 2009, 2, 183–202. [Google Scholar] [CrossRef] [Green Version]
  9. Li, S.; Peng, Z.; Wu, H. Prestack Multi-Gather Simultaneous Inversion of Elastic Parameters Using Multiple Regularization Constraints. J. Earth Sci. 2018, 29, 1359–1371. [Google Scholar] [CrossRef]
  10. Fu, X.; Zhang, F.; Li, X.; Qian, Z.; Chen, H.; Mei, L. Simultaneous inversion of P-and S-wave impedances based on an improved approximation equation of reflection coefficients. Chin. J. Geophys. 2019, 62, 276–288. (In Chinese) [Google Scholar]
  11. She, B.; Wang, Y.; Zhang, J.; Wang, J.; Hu, G. AVO inversion with high-order total variation regularization. J. Appl. Geophys. 2019, 161, 167–181. [Google Scholar] [CrossRef]
  12. Nie, W.; Wen, X.; Yang, J.; He, J.; Lin, K.; Yang, L. L1-2 minimization for P- and S-impedance inversion. Interpret.-A J. Subsurf. Charact. 2020, 8, T379–T390. [Google Scholar] [CrossRef]
  13. Connolly, P. Elastic impedance. Lead. Edge 1999, 18, 438–452. [Google Scholar] [CrossRef]
  14. Whitcombe, D.N.; Connolly, P.A.; Reagan, R.L.; Redshaw, T.C. Extended elastic impedance for fluid and lithology prediction. Geophysics 2002, 67, 63–67. [Google Scholar] [CrossRef]
  15. Zong, Z.; Yin, X.; Wu, G. AVO inversion and poroelasticity with P- and S-wave moduli. Geophysics 2012, 77, N17–N24. [Google Scholar] [CrossRef]
  16. Zong, Z.; Yin, X.; Wu, G. Elastic impedance parameterization and inversion with Young’s modulus and Poisson’s ratio. Geophysics 2013, 78, N35–N42. [Google Scholar] [CrossRef]
  17. Liu, X.J.; Yin, X.Y.; Wu, G.C.; Zong, Z.Y. Identification of deep reservoir fluids based on basis pursuit inversion for elastic impedance. Chin. J. Geophys. 2016, 59, 277–286. [Google Scholar]
  18. Zhang, B.; Chang, D.; Lin, T.; Marfurt, K.J. Improving the quality of prestack inversion by prestack data conditioning. Interpretation 2014, 3, T5–T12. [Google Scholar] [CrossRef]
  19. Zhang, B.; Zhao, T.; Jin, X.; Marfurt, K.J. Brittleness evaluation of resource plays by integrating petrophysical and seismic data analysis. Interpretation 2015, 3, T81–T92. [Google Scholar] [CrossRef] [Green Version]
  20. Yin, X.; Liu, X.; Zong, Z. Pre-stack basis pursuit seismic inversion for brittleness of shale. Commun. Inf. Syst. 2015, 12, 618–627. [Google Scholar] [CrossRef] [Green Version]
  21. Zhang, F.Q.; Jin, Z.J.; Sheng, X.J.; Li, X.S.; Shi, J.Y.; Liu, X.Q. A direct inversion for brittleness index based on GLI with basic-pursuit decomposition. Chin. J. Geophys. 2017, 60, 3954–3968. [Google Scholar]
  22. Wang, L.; Zhou, H.; Wang, Y.; Yu, B.; Zhang, Y.; Liu, W.; Chen, Y. Three parameters prestack seismic inversion based on L1-2 minimization. Geophys. J. Int. 2019, 84, R753–R766. [Google Scholar] [CrossRef]
  23. Lai, M.; Xu, Y.; Yin, W. Improved iteratively reweighted least squares for unconstrained smoothed Lq minimization. SIAM J. Numer. Anal. 2003, 51, 927–957. [Google Scholar] [CrossRef]
  24. Yin, P.; Esser, E.; Xin, J. Ratio and difference of L1 and L2 norms and sparse representation with coherent dictionaries. Commun. Inf. Syst. 2015, 14, 87–109. [Google Scholar] [CrossRef]
  25. Zhong, W.; Chen, Y.; Gan, S. L1/2 norm regularization for 3D seismic data interprolation. J. Seism. Explor. 2016, 25, 257. [Google Scholar]
  26. Huang, G.; Chen, X.; Luo, C.; Chen, Y. Pre-stack seismic inversion based on-norm regularized logarithmic absolute misfit function. Geophys. Prospect. 2020, 68, 2419–2443. [Google Scholar] [CrossRef]
  27. Huang, G.; Chen, X.; Chen, Y. PP and Dynamic Time Warped P-SV Wave AVA Joint-Inversion With L1-2 Regularization. IEEE Trans. Geosci. Remot Sens. 2020, 59, 5535–5548. [Google Scholar] [CrossRef]
  28. Lou, Y.; Yin, P.; He, Q.; Xin, J. Computing Sparse Representation in a Highly Coherent Dictionary Based on Difference of L1 and L2. J. Sci. Comput. 2015, 64, 178–196. [Google Scholar] [CrossRef]
  29. Zoeppritz, K. On the reflection and penetration of seismic waves through unstable layers. J. Geophys. Eng. 1919, 1, 66–84. [Google Scholar]
  30. Theune, U.; Jensås, I.Ø.; Eidsvik, J. Analysis of prior models for a blocky inversion of seismic AVA data. Geophysics 2010, 75, C25–C35. [Google Scholar] [CrossRef]
Figure 1. Model parameters from real well logs.
Figure 1. Model parameters from real well logs.
Applsci 11 12015 g001
Figure 2. The inversion results generated with noise-free synthetic records using the (a) L 1 2 and (b) BP inversion methods. Red curves represent the inversion results, blue curves represent the observed values, and the black curves were generated using the low-frequency model.
Figure 2. The inversion results generated with noise-free synthetic records using the (a) L 1 2 and (b) BP inversion methods. Red curves represent the inversion results, blue curves represent the observed values, and the black curves were generated using the low-frequency model.
Applsci 11 12015 g002
Figure 3. The synthetic seismic angle gathers with a SNR value of (a) 10 dB, (b) 5 dB, and (c) 1 dB.
Figure 3. The synthetic seismic angle gathers with a SNR value of (a) 10 dB, (b) 5 dB, and (c) 1 dB.
Applsci 11 12015 g003
Figure 4. The inversion results for synthetic seismic data with a SNR value of 10 dB using (a) our L 1 2 regularization inversion method and (b) the BP inversion method. The red, blue, and black curves represent the inverted, model, and initial values, respectively.
Figure 4. The inversion results for synthetic seismic data with a SNR value of 10 dB using (a) our L 1 2 regularization inversion method and (b) the BP inversion method. The red, blue, and black curves represent the inverted, model, and initial values, respectively.
Applsci 11 12015 g004
Figure 5. The inversion results for synthetic seismic data with a SNR value of 5 dB using (a) our L 1 2 regularization inversion method and (b) the BP inversion method. The red, blue, and black curves represent the inverted, model, and initial values, respectively.
Figure 5. The inversion results for synthetic seismic data with a SNR value of 5 dB using (a) our L 1 2 regularization inversion method and (b) the BP inversion method. The red, blue, and black curves represent the inverted, model, and initial values, respectively.
Applsci 11 12015 g005
Figure 6. The inversion results for synthetic seismic data with a SNR value of 1 dB using (a) our L 1 2 regularization inversion method and (b) the BP inversion method. The red, blue, and black curves represent the inverted, model, and initial values, respectively.
Figure 6. The inversion results for synthetic seismic data with a SNR value of 1 dB using (a) our L 1 2 regularization inversion method and (b) the BP inversion method. The red, blue, and black curves represent the inverted, model, and initial values, respectively.
Applsci 11 12015 g006
Figure 7. The inversion results for angle gathers of 3 , 18 , and 33 using (a) our L 1 2 regularization inversion method and (b) the BP inversion method. The red, blue, and black curves represent the inverted, model, and initial values, respectively.
Figure 7. The inversion results for angle gathers of 3 , 18 , and 33 using (a) our L 1 2 regularization inversion method and (b) the BP inversion method. The red, blue, and black curves represent the inverted, model, and initial values, respectively.
Applsci 11 12015 g007
Figure 8. The inversion results for angle gathers of 6 , 21 , and 36 using (a) our L 1 2 regularization inversion method and (b) the BP inversion method. The red, blue, and black curves represent the inverted, model, and initial values, respectively.
Figure 8. The inversion results for angle gathers of 6 , 21 , and 36 using (a) our L 1 2 regularization inversion method and (b) the BP inversion method. The red, blue, and black curves represent the inverted, model, and initial values, respectively.
Applsci 11 12015 g008
Figure 9. The inversion results for angle gathers of 9 , 24 , and 39 using (a) our L 1 2 regularization inversion method and (b) the BP inversion method. The red, blue, and black curves represent the inverted, model, and initial values, respectively.
Figure 9. The inversion results for angle gathers of 9 , 24 , and 39 using (a) our L 1 2 regularization inversion method and (b) the BP inversion method. The red, blue, and black curves represent the inverted, model, and initial values, respectively.
Applsci 11 12015 g009
Figure 10. The inversion results for angle gathers of 12 , 27 , and 42 using (a) our L 1 2 regularization inversion method and (b) the BP inversion method. The red, blue, and black curves represent the inverted, model, and initial values, respectively.
Figure 10. The inversion results for angle gathers of 12 , 27 , and 42 using (a) our L 1 2 regularization inversion method and (b) the BP inversion method. The red, blue, and black curves represent the inverted, model, and initial values, respectively.
Applsci 11 12015 g010
Figure 11. The inversion results for angle gathers of 15 , 30 , and 45 using (a) our L 1 2 regularization inversion method and (b) the BP inversion method. The red, blue, and black curves represent the inverted, model, and initial values, respectively.
Figure 11. The inversion results for angle gathers of 15 , 30 , and 45 using (a) our L 1 2 regularization inversion method and (b) the BP inversion method. The red, blue, and black curves represent the inverted, model, and initial values, respectively.
Applsci 11 12015 g011
Figure 12. Real partial angle stacks for angle gathers of (a) 12 , (b) 24 , and (c) 39 .
Figure 12. Real partial angle stacks for angle gathers of (a) 12 , (b) 24 , and (c) 39 .
Applsci 11 12015 g012
Figure 13. (a) P-impedance, (b) S-impedance, and (c) density low-frequency models.
Figure 13. (a) P-impedance, (b) S-impedance, and (c) density low-frequency models.
Applsci 11 12015 g013
Figure 14. (a) P-impedance, (b) S-impedance, (c) density, and (d) VP/VS inversion results. The black circles denote the positions of gas-bearing sandstone units.
Figure 14. (a) P-impedance, (b) S-impedance, (c) density, and (d) VP/VS inversion results. The black circles denote the positions of gas-bearing sandstone units.
Applsci 11 12015 g014
Figure 15. Synthetic vs. actual records for angle gathers of 3 , 18 , and 33 using (a) the L 1 2 inversion method and (b) the BP inversion method.
Figure 15. Synthetic vs. actual records for angle gathers of 3 , 18 , and 33 using (a) the L 1 2 inversion method and (b) the BP inversion method.
Applsci 11 12015 g015
Figure 16. Synthetic vs. actual records for angle gathers of 6 , 21 , and 36 using (a) the L 1 2 inversion method and (b) the BP inversion method.
Figure 16. Synthetic vs. actual records for angle gathers of 6 , 21 , and 36 using (a) the L 1 2 inversion method and (b) the BP inversion method.
Applsci 11 12015 g016
Figure 17. Synthetic vs. actual records for angle gathers of 9 , 24 , and 39 using (a) the L 1 2 inversion method and (b) the BP inversion method.
Figure 17. Synthetic vs. actual records for angle gathers of 9 , 24 , and 39 using (a) the L 1 2 inversion method and (b) the BP inversion method.
Applsci 11 12015 g017
Figure 18. Synthetic vs. actual records for angle gathers of 12 , 27 , and 42 using (a) the L 1 2 inversion method and (b) the BP inversion method.
Figure 18. Synthetic vs. actual records for angle gathers of 12 , 27 , and 42 using (a) the L 1 2 inversion method and (b) the BP inversion method.
Applsci 11 12015 g018
Figure 19. Synthetic vs. actual records for angle gathers of 15 , 30 , and 45 using (a) the L 1 2 inversion method and (b) the BP inversion method.
Figure 19. Synthetic vs. actual records for angle gathers of 15 , 30 , and 45 using (a) the L 1 2 inversion method and (b) the BP inversion method.
Applsci 11 12015 g019
Table 1. Groups of angle gathers used in inversion experiments.
Table 1. Groups of angle gathers used in inversion experiments.
Small AngleMedium AngleLarge Angle
Group 1 3 18 33
Group 2 6 21 36
Group 3 9 24 39
Group 4 12 27 42
Group 5 15 30 45
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Nie, W.; Xiang, F.; Li, B.; Wen, X.; Nie, X. Prestack Seismic Inversion via Nonconvex L1-2 Regularization. Appl. Sci. 2021, 11, 12015. https://doi.org/10.3390/app112412015

AMA Style

Nie W, Xiang F, Li B, Wen X, Nie X. Prestack Seismic Inversion via Nonconvex L1-2 Regularization. Applied Sciences. 2021; 11(24):12015. https://doi.org/10.3390/app112412015

Chicago/Turabian Style

Nie, Wenliang, Fei Xiang, Bo Li, Xiaotao Wen, and Xiangfei Nie. 2021. "Prestack Seismic Inversion via Nonconvex L1-2 Regularization" Applied Sciences 11, no. 24: 12015. https://doi.org/10.3390/app112412015

APA Style

Nie, W., Xiang, F., Li, B., Wen, X., & Nie, X. (2021). Prestack Seismic Inversion via Nonconvex L1-2 Regularization. Applied Sciences, 11(24), 12015. https://doi.org/10.3390/app112412015

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop