Low-Dimensional Multi-Trace Impedance Inversion in Sparse Space with Elastic Half Norm Constraint

Abstract

Recently, multi-trace impedance inversion has attracted great interest in seismic exploration because it improves the horizontal continuity and fidelity of the inversion results by exploiting the lateral structure information of the strata. However, computational inefficiency affects its practical application. Furthermore, in terms of vertical constraints on the model parameters, it only considers smooth features while ignoring sharp discontinuity features. This leads to yielding an over-smooth solution that does not accurately reflect the distribution of underground rock. To deal with the above situations, we first develop a low-dimensional multi-trace impedance inversion (LMII) framework. Inspired by compressed sensing, this framework utilizes low-dimensional measurements in sparse space containing the maximum information of the signal to construct the objective function for multi-trace inversion, which can significantly reduce the size of the inversion problem and improve the inverse efficiency. Then, we introduce the elastic half (EH) norm as a vertical constraint on the model parameters in the LMII framework and formulate a novel constrained LMII model for impedance inversion. Because the introduced EH norm takes into account both the smoothness and blockiness of rock impedance, the constrained LMII model can effectively raise the inversion accuracy of complex strata. Finally, an efficient alternating multiplier iteration algorithm is derived based on the variable splitting technique to optimize the constrained LMII model. The performance of the developed approaches is tested using synthetic and practical data, and the results prove their feasibility and superiority.

1. Introduction

Impedance is a critical parameter to describe the properties of underground rock, and its accurate estimation is essential for mineral zone prediction and deposit characterization [1,2,3]. In mineral exploration, seismic inversion is a major tool to obtain the impedance of subsurface rocks. Affected by the band-limited nature of seismic data, incomplete data coverage, and noise interference, seismic inversion is usually ill-conditioned. Generally, it is difficult to obtain a stable and unique impedance profile [4,5]. To overcome the ill-conditions of impedance inversion and obtain reliable inversion results, regularized inversion techniques that introduce additional prior information have been extensively studied by exploration geophysicists.

Currently, the regularized inversion methods available for impedance prediction can be divided into two main categories: Type one is Tikhonov regularization inversion. This kind of method employs the L2 norm [6,7] or Gaussian prior distribution [8] to constrain the solution and force it to satisfy the smoothness condition. Unfortunately, these methods can blur sharp impedance edges and produce over-smooth solutions [9,10], which is detrimental to fine reservoir characterization. Another type is sparsity regularized inversion. This type of approach constrains the solution by using sparse norms or long-tailed distributions, imposing it to exhibit blocky features. Typical sparse regularization for impedance inversion includes total-variation regularization [11,12,13,14], L0 norm [15], L1 norm [16,17], L1-2 norm [18,19], Lp (0 < p < 1) norm [20], and Cauchy distribution [21,22]. Although sparse regularization preserves important edge information, it is inadequate for describing the prior information of rock impedance. The reason is that subsurface strata usually contain edge discontinuities that separate two smooth regions [9,10], while sparse regularization can only identify a single sharp property and cannot effectively characterize the smooth structures of stratigraphy. Similar to sparse regularization, Tikhonov regularization is also unable to model both smooth and sharp features of subsurface strata. Therefore, the regularized inversion methods mentioned above are ineffective in obtaining unbiased impedance estimation [14,23]. Additionally, these methods perform a trace-by-trace inversion, which ignores the lateral coherence between seismic traces [24,25,26,27,28].

Recently, several multi-trace impedance inversion (MII) methods were developed to eliminate the lateral instability problem caused by trace-by-trace inversion. Hamid and Pidlisecky present a laterally constrained MII algorithm [24], which constructs a constraint term in the horizontal direction using a second-order derivative operator to improve the lateral continuity of the inverted impedance. Karimi discusses the limitations of the MII method in the presence of tilted layers [29]. Hamid and Pidlisecky integrate seismic dip information into the impedance inversion and develop an MII method with structural constraints [26]. This method rotates the derivative operator by employing dip information to force the inversion results to honor the local structure, which enhances the adaptability of the MII method to high-steep layers. Yin et al. propose a cross-correlation-driven MII method that uses cross-correlation to describe the structural characteristics of stratigraphic reflections and raises the stability of impedance inversion [30]. Zhang et al. extract the local structural direction using the structure tensor and design a structure-oriented regularization method for impedance inversion [31]. Although the above MII techniques have made important contributions to improving the spatial continuity and stability of inversion results, there are still two critical issues to be addressed. One is the efficiency issue. The above MII approaches need to expand the forward operator via the Kronecker product to construct the objective function. Since the size of the forward operator is exponentially related to the number of traces and sampling points, the inverse problem has a large scale when both are enormous. Solving the large-scale inverse problem is so time-consuming that the MII method is challenged in practical application. Another one is the accuracy issue. The above MII approaches all construct the constraint terms with the L2 norm, which implies these methods only consider the smoothness and ignore the sharp discontinuity features in terms of vertical constraints on the model parameters. Thus, these methods cannot produce an optimal solution that accurately reflects the impedance distribution.

To address the above issues as much as possible, we propose a low-dimensional multi-trace impedance inversion (LMII) method with elastic half (EH) norm constraints in this paper. Specifically, we first develop a novel LMII inversion framework. Enlightened by compressed sensing [32,33], this framework uses low-dimensional measurements in a sparse domain containing the maximum information of the signal to construct the objective function for multi-trace inversion, which can effectively reduce the size of the inversion problem and improve the inversion efficiency. Subsequently, we introduce the sum of L1/2 norm and L2 norm, called EH norm [34], as a vertical constraint on the model parameters in the LMII framework to obtain an unbiased impedance estimation. As a combination of sparse regularization and Tikhonov regularization, the EH norm is able to characterize both blocky and smooth features of underground rock. Therefore, the LMII model constrained by the EH norm can effectively raise the inversion accuracy compared with the traditional MII method that only considers smoothness. Then, based on the variable splitting technique [35,36], we derive an efficient alternating multiplier iteration algorithm to optimize the constrained LMII problem. Finally, the effectiveness and superiority of the developed methods are verified by synthetic and field data.

2. Methodology

2.1. MII Method

According to the Robinson convolution model [37], seismic records can be represented as a convolution of source wavelet and reflection coefficient series:

$$\mathbf{s}=\mathbf{w}\ast \mathbf{r}+n,$$

(1)

where $\mathbf{s}\in {ℝ}^M$ denotes the seismic record, $\mathbf{w}\in {ℝ}^M$ denotes the source wavelet, $\mathbf{r}\in {ℝ}^M$ denotes the reflection coefficient series, and $\mathbf{n}\in {ℝ}^M$ denotes the noise. Reformulating Equation (1) as matrix-vector product yields:

$$\mathbf{s}=\mathbf{W}\mathbf{r}+n,$$

(2)

where $\mathbf{W}\in {ℝ}^{M\times M}$ is the circular matrix formed by source wavelet time-shifting. Under the small reflection coefficient hypothesis [37], the following linearized relationship exists between the subsurface impedance and reflection coefficient:

$$\mathbf{r}\left(t\right)=\frac{1}{2}\ln \left(\frac{\overline{\mathbf{z}}\left(t+1\right)}{\overline{\mathbf{z}}\left(t\right)}\right)=\frac{1}{2}\left(\ln \left(\overline{\mathbf{z}}\left(t+1\right)\right)-\ln \left(\overline{\mathbf{z}}\left(t\right)\right)\right),$$

(3)

where $\overline{\mathbf{z}}$ denotes impedance and $t$ represents time. By arranging the reflection coefficients along the time axis, one obtains:

$$\mathbf{r}=\frac{1}{2}\mathbf{D}\mathbf{z},$$

(4)

where $\mathbf{D}$ is the first-order derivative operator in the time direction and $\mathbf{z}$ represents the natural logarithm of the stratigraphic impedance. Combining Equations (2) and (4) produces the following forward problem

$$\mathbf{s}=\frac{1}{2}\mathbf{W}\mathbf{D}\mathbf{z}+\mathbf{n}=\mathbf{G}\mathbf{z}+n.$$

(5)

In the multi-trace case, Equation (5) can be expanded to the following form

$$\underset{\mathbf{S}}{\underbrace{\left[\begin{array}{c}{\mathbf{s}}_1\\ {\mathbf{s}}_2\\ ⋮\\ {\mathbf{s}}_N\end{array}\right]}}=\underset{\mathbf{G}}{\underbrace{\left[\begin{array}{cccc}{\mathbf{G}}_1& \mathbf{0}& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& {\mathbf{G}}_2& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& \mathbf{0}& \ddots & \mathbf{0}\\ \mathbf{0}& \mathbf{0}& \mathbf{0}& {\mathbf{G}}_N\end{array}\right]}}\underset{\mathbf{Z}}{\underbrace{\left[\begin{array}{c}{\mathbf{z}}_1\\ {\mathbf{z}}_2\\ ⋮\\ {\mathbf{z}}_N\end{array}\right]}}+\underset{\mathbf{N}}{\underbrace{\left[\begin{array}{c}{n}_1\\ n_2\\ ⋮\\ n_N\end{array}\right]}}.$$

(6)

To obtain the subsurface impedance, a common method is to formulate a least-square optimization problem:

$${\mathbf{Z}}^{*}=\arg \underset{Z}{\min }\frac{1}{2}{\Vert \mathbf{S}-\mathbf{G}\mathbf{Z}\Vert }_2^2,$$

(7)

where ${\mathbf{Z}}^{*}$ denotes the optimal estimation of natural logarithmic impedance $\mathbf{Z}$. Due to the ill-conditioning feature of $\mathbf{G}$ and the presence of noise, Equation (7) is often ill-posed. Thus, Equation (7) needs to introduce prior information through regularization terms to obtain a stable and reliable solution. In the conventional MII method [24], lateral structure as well as reference model constraints are inserted to improve the inversion accuracy. Specifically, the objective function of the conventional MII method can be expressed as

$$\begin{array}{ll}{\mathbf{Z}}^{*}& =\arg \underset{\mathbf{Z}}{\min }\frac{1}{2}{\Vert S-\mathbf{G}\mathbf{Z}\Vert }_2^2+\lambda {\Vert {\mathbf{Z}}_{ref}-\mathbf{Z}\Vert }_2^2+\mu {\Vert \mathbf{C}\mathbf{Z}\Vert }_2^2,\\ & =\arg \underset{\mathbf{Z}}{\min }\frac{1}{2}{\Vert \overline{\mathbf{S}}-\overline{\mathbf{G}}\mathbf{Z}\Vert }_2^2\end{array}$$

(8)

where $\overline{S}=\left[\begin{array}{c}S\\ \sqrt{2\lambda }{\mathbf{Z}}_{ref}\\ \mathbf{0}\end{array}\right]\in {ℝ}^{3MN}$, $\overline{\mathbf{G}}=\left[\begin{array}{c}\mathbf{G}\\ \sqrt{2\lambda }\mathbf{I}\\ -\sqrt{2\mu }\mathbf{C}\end{array}\right]\in {ℝ}^{3MN\times MN}$, ${\mathbf{Z}}_{ref}$ is the low-frequency reference model constructed from geological knowledge or well log, $\mathbf{I}$ is the identity matrix, $\mathbf{C}$ is the second-order difference operator in the spatial direction, $\lambda$ and $\mu$ are the weight factors. Equation (8) has an analytical solution, that is:

$$Z^{*}={\left({\overline{G}}^T\overline{G}\right)}^{-1}{\overline{G}}^T\overline{S},$$

(9)

where ${\overline{G}}^T$ represents the transpose of $\overline{G}$. Equation (9) involves the inverse operation of large-scale matrices and using it directly for impedance inversion is very time-consuming. Alternatively, the conjugate gradient method [38], which does not involve an inverse matrix, is usually adopted to solve $Z^{*}$. However, it still suffers from computational inefficiency due to the limitation of matrix scale. In order to address this issue, a novel LMII framework is developed in this paper.

2.2. LMII Framework

Given an orthogonal transform matrix $\mathbf{\Phi }$, whose columns are orthogonal atoms ${\left\{{\phi }_i\right\}}_{i=1}^M$, any M-dimensional discrete signal can be represented as a linear combination of these atoms:

$$\mathbf{x}={\displaystyle \sum _{i=1}^M{\phi }_i{\alpha }_i}=\mathbf{\Phi }\mathit{\alpha },$$

(10)

where ${\alpha }_i=\langle {\phi }_i,\mathbf{x}\rangle$ is the projection coefficient of signal $\mathbf{x}$ on the i-th atom and $\langle •,•\rangle$ denotes the Euclidean inner product operator. Compressive sensing [32,33] proved that, if the projection $\mathit{\alpha }$ is sparse, then K $\left(K\ll M\right)$ non-adaptive measurements of $\mathit{\alpha }$ observed by a measurement matrix satisfying certain conditions are sufficient to accurately recover $\mathbf{x}$. Mathematically, the above K-dimensional measurements can be expressed as:

$$\mathbf{f}=\mathbf{\Psi }\mathit{\alpha },$$

(11)

where $\mathbf{f}\in {ℝ}^K$ is the low-dimensional measurement and $\mathbf{\Psi }\in {ℝ}^{K\times M}$ is the measurement matrix. Theoretically, the condition to be satisfied by the measurement matrix can be described by the restricted isometry (RIP) property [39] as follows:

$$\left(1-{\delta }_K\right){\Vert \mathit{\alpha }\Vert }_2^2\le {\Vert \mathbf{\Psi }\mathit{\alpha }\Vert }_2^2\le \left(1+{\delta }_K\right){\Vert \mathit{\alpha }\Vert }_2^2,$$

(12)

where ${\delta }_k\in \left[0,1\right]$ is the isometry constant of $\mathbf{\Psi }$. Essentially, ${\delta }_k$ quantitatively describes the preservation for the information contained in the original signal by low-dimensional measurements. ${\delta }_K=0$ indicates complete preservation, while ${\delta }_K=1$ indicates no preservation. To recover the signal exactly, ${\delta }_K$ needs to be as small as possible, which means that the subsets of K columns of $\mathbf{\Psi }$ are close to a set of standard orthogonal basis [40].

Inspired by the above idea, we develop a novel LMII framework in sparse space. Concretely, by employing orthogonal transform to convert the multi-trace inversion objective function from the spatio-temporal domain to sparse space, one can obtain:

$${\mathbf{Z}}^{*}=\arg \underset{\mathbf{Z}}{\min }\frac{1}{2}{\Vert S^{\prime }-{\mathbf{F}}_1\mathbf{G}\mathbf{Z}\Vert }_2^2+\lambda {\Vert {{\mathbf{Z}}^{\prime }}_{ref}-{\mathbf{F}}_2\mathbf{Z}\Vert }_2^2+\mu {\Vert \mathbf{C}\mathbf{Z}\Vert }_2^2,$$

(13)

where ${\mathbf{F}}_1=kron\left({\mathbf{I}}_{N\times N},{\mathbf{\Phi }}_1\right)$, ${\mathbf{F}}_2=kron\left({\mathbf{I}}_{N\times N},{\mathbf{\Phi }}_2\right)$, $kron\left(•,•\right)$ denotes the Kronecker product operator, ${\mathbf{\Phi }}_1$ and ${\mathbf{\Phi }}_2$ denote the employed sparse transform, and ${\mathbf{S}}^{\prime }={\mathbf{F}}_1\mathbf{S}$ and ${{\mathbf{Z}}^{\prime }}_{ref}={\mathbf{F}}_2{\mathbf{Z}}_{ref}$ denote projections of the seismic profile and the low-frequency reference model in sparse space, respectively. Afterwards, the measurement matrix is designed according to the RIP property and the low-dimensional measurement is performed in the sparse space to obtain the following LMII framework:

$$\begin{array}{ll}{\mathbf{Z}}^{*}& =\arg \underset{\mathbf{Z}}{\min }\frac{1}{2}{\Vert \hat{\mathbf{S}}-{\mathbf{M}}_1{\mathbf{F}}_1\mathbf{G}\mathbf{Z}\Vert }_2^2+\lambda {\Vert {\hat{\mathbf{Z}}}_{ref}-{\mathbf{M}}_2{\mathbf{F}}_2\mathbf{Z}\Vert }_2^2+\mu {\Vert \mathbf{C}\mathbf{Z}\Vert }_2^2\\ & =\arg \underset{\mathbf{Z}}{\min }\frac{1}{2}{\Vert \tilde{S}-\tilde{\mathbf{G}}\mathbf{Z}\Vert }_2^2\end{array}$$

(14)

where ${\mathbf{{\rm M}}}_1=kron\left({\mathbf{I}}_{N\times N},{\mathbf{\Psi }}_1\right)$, ${\mathbf{{\rm M}}}_2=kron\left({\mathbf{I}}_{N\times N},{\mathbf{\Psi }}_2\right)$, ${\mathbf{\Psi }}_1$, and ${\mathbf{\Psi }}_2$ denote the designed measurement matrix and $\hat{\mathbf{S}}={\mathbf{M}}_1{S}^{\prime }$ and ${\hat{\mathbf{Z}}}_{ref}={\mathbf{M}}_2{{\mathbf{Z}}^{\prime }}_{ref}$ represent the low-dimensional measurements of seismic profile and low-frequency reference model in the sparse space, $\tilde{S}=\left[\begin{array}{c}\hat{\mathbf{S}}\\ \sqrt{2\lambda }{\hat{\mathbf{Z}}}_{ref}\\ \mathbf{0}\end{array}\right]\in {ℝ}^{\left({\kappa }_1+{\kappa }_2+1\right)MN}$, $\tilde{\mathbf{G}}=\left[\begin{array}{c}{\mathbf{{\rm M}}}_1{\mathbf{F}}_1\mathbf{G}\\ \sqrt{2\lambda }{\mathbf{M}}_2{\mathbf{F}}_2\\ -\sqrt{2\mu }\mathbf{C}\end{array}\right]\in {ℝ}^{\left({\kappa }_1+{\kappa }_2+1\right)MN\times MN}$, ${\kappa }_1,{\kappa }_2\in \left(0,1\right]$ denote the measurement ratios of seismic profile and reference model, respectively. Since both ${\kappa }_1$ and ${\kappa }_2$ are less than 1, the scale of the LMII problem is significantly smaller than that of the traditional MII problem. Obviously, this provides a reliable guarantee for the efficient inversion of rock impedance.

2.3. Constrained LMII Model via EH Norm

Although the developed LMII framework effectively improves solution efficiency, it still has limitations in accuracy (the conventional MII method also has this problem). Concretely, in terms of vertical constraints on the model parameters, this method only considers smooth features and ignores sharp discontinuity features, which leads to its solution being too smooth to accurately reflect the distribution of underground rock. To address this issue, we introduce the EH norm [34] in the LMII framework as a vertical constraint on the model parameters and formulate a constrained LMII model for impedance inversion. Since the introduced EH norm takes into account both the smoothness and blockiness of the rock impedance, the constrained LMII model can raise the inversion accuracy of complex stratigraphy. Formulaically, the constrained LMII model can be written as:

$${\mathbf{Z}}^{*}=\arg \underset{\mathbf{Z}}{\min }\frac{1}{2}{\Vert \tilde{S}-\tilde{\mathbf{G}}\mathbf{Z}\Vert }_2^2+\mu {\Vert \mathbf{C}\mathbf{Z}\Vert }_2^2+\gamma {\Vert \mathbf{B}\mathbf{Z}\Vert }_{1/2}^{1/2}+\beta {\Vert \mathbf{B}\mathbf{Z}\Vert }_2^2,$$

(15)

where $\stackrel{⌢}{\mathbf{S}}=\left[\begin{array}{c}\hat{\mathbf{S}}\\ \sqrt{2\lambda }{\hat{\mathbf{Z}}}_{ref}\end{array}\right]$, $\stackrel{⌢}{\mathbf{G}}=\left[\begin{array}{c}{\mathbf{{\rm M}}}_1{\mathbf{F}}_1\mathbf{G}\\ \sqrt{2\lambda }{\mathbf{M}}_2{\mathbf{F}}_2\end{array}\right]$, $\mathbf{B}=kron\left({\mathbf{I}}_{N\times N},\mathbf{D}\right)$, ${\Vert •\Vert }_{1/2}^{1/2}$ denotes the L_1/2 norm and $\gamma$ and $\beta$ are weighting factors. Equation (15) can be equivalently converted into the following optimization problem:

$${\mathbf{Z}}^{*}=\arg \underset{\mathbf{Z}}{\min }\frac{1}{2}{\Vert \stackrel{⌢}{\mathbf{S}}-\stackrel{⌢}{\mathbf{G}}\mathbf{Z}\Vert }_2^2+\mu {\Vert \mathbf{A}\mathbf{Z}\Vert }_2^2+\gamma {\Vert \mathbf{B}\mathbf{Z}\Vert }_{1/2}^{1/2},$$

(16)

where $\mathbf{A}=\mathbf{C}+\sqrt{\beta /\mu }\mathbf{B}$. To optimize Equation (16), we derive an efficient alternating multiplier iteration algorithm based on variable splitting technique. Specifically, with the aid of the variable splitting technique [35,36], the following constrained minimization problem is generated by introducing an auxiliary variable $\mathbf{T}=\mathbf{B}\mathbf{Z}$ into Equation (16)

$$\left({\mathbf{Z}}^{*},{\mathbf{T}}^{*}\right)=\arg \underset{\mathbf{Z},\mathbf{T}}{\min }\frac{1}{2}{\Vert \stackrel{⌢}{\mathbf{S}}-\stackrel{⌢}{\mathbf{G}}\mathbf{Z}\Vert }_2^2+\mu {\Vert \mathbf{A}\mathbf{Z}\Vert }_2^2+\gamma {\Vert \mathbf{T}\Vert }_{1/2}^{1/2} \mathrm{s}.\mathrm{t}. \mathbf{T}=\mathbf{B}\mathbf{Z}.$$

(17)

Equation (17) can be minimized in an unconstrained manner by formulating it as an augmented Lagrangian function. To be specific, its augmented Lagrangian function is:

$$ℒ\left(\mathbf{Z},\mathbf{T},\mathbf{c}\right)=\underset{\mathbf{Z},\mathbf{T},\mathbf{c}}{\min }\frac{1}{2}{\Vert \stackrel{⌢}{\mathbf{S}}-\stackrel{⌢}{\mathbf{G}}\mathbf{Z}\Vert }_2^2+\mu {\Vert \mathbf{A}\mathbf{Z}\Vert }_2^2+\gamma {\Vert \mathbf{T}\Vert }_{1/2}^{1/2}+\langle \mathbf{B}\mathbf{Z}-\mathbf{T},\mathbf{c}\rangle +\frac{\varsigma }{2}{\Vert \mathbf{B}\mathbf{Z}-\mathbf{T}\Vert }_2^2,$$

(18)

where $\mathbf{c}$ is a Lagrangian multiplier and $\varsigma$ is a non-negative parameter that controls the convergence rate of the algorithm. The minimizer of Equation (18) is the saddle point of function $ℒ\left(\mathbf{Z},\mathbf{T},\mathbf{c}\right)$, which can be computed by the following iterative steps:

$${\mathbf{Z}}^{\left(k+1\right)}=\arg \underset{\mathbf{Z}}{\min }ℒ\left(\mathbf{Z},{\mathbf{T}}^{\left(k\right)},{\mathbf{u}}^{\left(k\right)}\right)=\arg \underset{\mathbf{Z}}{\min }\frac{1}{2}{\Vert \stackrel{⌢}{\mathbf{S}}-\stackrel{⌢}{\mathbf{G}}\mathbf{Z}\Vert }_2^2+\mu {\Vert \mathbf{A}\mathbf{Z}\Vert }_2^2+\frac{\varsigma }{2}{\Vert \mathbf{B}\mathbf{Z}-{\mathbf{T}}^{\left(k\right)}+{\mathbf{u}}^{\left(k\right)}\Vert }_2^2,$$

(19)

$${\mathbf{T}}^{\left(k+1\right)}=\arg \underset{\mathbf{T}}{\min }ℒ\left({\mathbf{Z}}^{\left(k+1\right)},\mathbf{T},{\mathbf{u}}^{\left(k\right)}\right)=\arg \underset{\mathbf{T}}{\min }\gamma {\Vert \mathbf{T}\Vert }_{1/2}^{1/2}+\frac{\varsigma }{2}{\Vert \mathbf{B}{\mathbf{Z}}^{\left(k+1\right)}-\mathbf{T}+{\mathbf{u}}^{\left(k\right)}\Vert }_2^2,$$

(20)

$${\mathbf{u}}^{\left(k+1\right)}={\mathbf{u}}^{\left(k\right)}+\mathbf{B}{\mathbf{Z}}^{\left(k+1\right)}-{\mathbf{T}}^{\left(k+1\right)},$$

(21)

where $\mathbf{u}=\mathbf{c}/\varsigma$ is the scaled Lagrangian multiplier and $k$ denotes the iteration index. Equation (19) is a quadratic equation whose solution can be obtained by first order necessary conditions. Specifically, letting its derivative of variable $\mathbf{Z}$ be zero yields

$$-{\stackrel{⌢}{\mathbf{G}}}^T\left(\stackrel{⌢}{\mathbf{S}}-\stackrel{⌢}{\mathbf{G}}\mathbf{Z}\right)+2\mu {\mathbf{A}}^T\mathbf{A}\mathbf{Z}+\varsigma {\mathbf{B}}^T\left(\mathbf{B}\mathbf{Z}-{\mathbf{T}}^{\left(k\right)}+{\mathbf{u}}^{\left(k\right)}\right)=0.$$

(22)

By simplifying Equation (22), the optimal solution of Equation (19) can be obtained:

$${\mathbf{Z}}^{\left(k+1\right)}={\left({\stackrel{⌢}{\mathbf{G}}}^T\stackrel{⌢}{\mathbf{G}}+2\mu {\mathbf{A}}^T\mathbf{A}+\varsigma {\mathbf{B}}^T\mathbf{B}\right)}^{-1}\left({\stackrel{⌢}{\mathbf{G}}}^T\stackrel{⌢}{\mathbf{S}}+\varsigma {\mathbf{B}}^T\left({\mathbf{T}}^{\left(k\right)}-{\mathbf{u}}^{\left(k\right)}\right)\right).$$

(23)

Equation (20) is a typical L1/2 norm optimization problem [41], whose solution can be computed by the following half-threshold function:

$${\mathbf{T}}^{\left(k+1\right)}={\mathscr{H}}_{\lambda }\left(\mathbf{B}{\mathbf{Z}}^{\left(k+1\right)}+{\mathbf{u}}^{\left(k\right)}\right),$$

(24)

where

$${\mathscr{H}}_{\lambda }\left(\mathbf{x}\right)=\{\begin{array}{ll}{f}_{\lambda }\left({\mathbf{x}}_i\right)& , \left|{\mathbf{x}}_i\right|>\frac{3\sqrt[3]{2}}{4}{\lambda }^{\frac{2}{3}}\\ 0,& \mathrm{otherwise}\end{array},$$

(25)

$$f_{\lambda }\left({\mathbf{x}}_i\right)=\frac{2{\mathbf{x}}_i}{3}\left(1+\cos \left(\frac{2\pi }{3}-\frac{2{\phi }_{\lambda }\left({\mathbf{x}}_i\right)}{3}\right)\right),$$

(26)

$${\phi }_{\lambda }\left({\mathbf{x}}_i\right)=\arccos \left(\frac{\lambda }{8}{\left(\frac{\left|{\mathbf{x}}_i\right|}{3}\right)}^{-3/2}\right),$$

(27)

where $\lambda =\gamma /\varsigma$. Updating the variables $\mathbf{Z}$, $\mathbf{T}$, $\mathbf{u}$ by Equations (21), (23), and (24) until the change of objective function is below the threshold level, we attain the optimal natural logarithmic impedance ${\mathbf{Z}}^{*}$. Note that the computational bottleneck of Equation (23) is the matrix inversion. Similar to solving Equation (9), this equation also needs to be solved by the simple and efficient conjugate gradient algorithm. Finally, perform the following exponential operation:

$$\overline{\mathbf{z}}=\exp \left(Z^{*}\right),$$

(28)

and the final predicted rock impedance can be obtained.

3. Examples

In this section, we implement two experiments (one synthetic and one field) to evaluate the feasibility and effectiveness of our presented methods. For our methods, it is a critical task to select appropriate sparse transforms and measurement matrices. In all examples, we employ the discrete Hartley transform (DHT) [42] as a tool for domain conversion because it can concentrate signal energy in the low-frequency part of the spectrum (see Figure 1). Meanwhile, we employ a binary matrix with orthogonal rows as the measurement matrix, which greatly avoids the loss of information contained in the original signal and can effectively maintain the accuracy of impedance inversion. For comparison, the conventional MII method is used as a benchmark. To be fair, the parameters of comparison methods are carefully tuned several times to yield optimal inversion results in each case.

3.1. Synthetic Example

First, we check the effectiveness of the developed methods with the SEG/EAGE overthrust impedance model (Figure 2a). Figure 2b depicts the noise-free seismic profile synthesized by convolving a 30 Hz Ricker wavelet with the reflection coefficients calculated in Figure 2a. This profile contains 801 traces and 373 sampling points with a sampling interval of 2 ms. The conventional MII, LMII, and constrained LMII methods are used for impedance inversion, and their inversion results are presented in Figure 3a–c, respectively. For the above methods, they are set to the same stopping criterion, i.e., the algorithm stops when the number of iterations is greater than 500 or the tolerance error is less than 10⁻⁶. It should be noted that the low-frequency reference model (Figure 2c) utilized in the above three methods is obtained by performing 10 Hz low-pass filtering on Figure 2a. In addition, the measurement ratios of seismic data and reference model are set to 0.45 and 0.1 in the LMII and constrained LMII methods, respectively. As shown in Figure 3a–c, all these methods effectively recover small-scale geological information such as thin layers, lenses, pinch-outs, and channels. Figure 3d–f depicts the difference between the results inverted by the above three methods and the true impedance. As observed from Figure 3d–f, the residuals of the conventional MIII method and the LMII method are basically the same, while the residuals of the constrained LMII method are significantly smaller than the other two methods. This indicates that introducing the EH norm effectively improves inversion accuracy. In other words, the EH norm is more accurate than the conventional Tikhonov regularization in characterizing the subsurface structure. To visually compare the advantages and disadvantages of the three methods, Figure 4 displays the zoomed-in views of the true impedance and the inversion results of Figure 3. Figure 4 demonstrates that the inversion results obtained by the MII and LMII methods are oversmooth, blurring the stratigraphic boundaries and producing additional pseudo-layer interferences (as shown by the arrows). In contrast, the inversion results yielded by the constrained LMII method show stratigraphic boundaries more clearly, while greatly avoiding erroneous pseudo-layer interferences.

In order to quantitatively evaluate the reliability of these inverted results, root mean square error (RMSE) defined by the following equation is used as the evaluation index:

$$RMSE=\sqrt{\frac{1}{M\times N}{\displaystyle \sum _{i=1}^M{\displaystyle \sum _{j=1}^N{\left({\overline{\mathbf{z}}}_{ij}-{\mathbf{y}}_{ij}\right)}^2}}}$$

(29)

where $\mathbf{y}$ denotes the true impedance, $\overline{\mathbf{z}}$ denotes the predicted impedance, $M$ and $N$ represent the number of sampling points and the number of seismic traces, respectively. Note that a smaller RMSE means higher reliability of the predicted impedance. The quantitative comparison of the inverted results of Figure 3 is listed in

Table 1. Quantitative evaluation index and running time of the inversion results of Figure 3.

Method	RMSE (×10−3)	Time (s)
MII	46.10	261.89
LMII	46.11	135.13
Constrained LMII	29.06	161.17

. As seen in Table 1, the constrained LMII method achieves a smaller RMSE value than the other two methods. This further demonstrates that the constrained LMII method can recover the impedance of subsurface rock more precisely. Moreover, the RMSE values obtained by the MII and LMII methods are basically consistent, but the calculation time of the conventional MII method is obviously longer than that of the LMII method. This provides reliable proof that the LMII method can effectively improve inversion efficiency while maintaining inversion accuracy. It should be pointed out that although the constrained LMII method uses an efficient alternating multiplier iteration algorithm to solve its objective function, it is slightly longer than the LMII method in terms of computational time due to the non-convex optimization involved. Notwithstanding, the constrained LMII method still has an advantage in computational efficiency compared to the conventional MII method.

To test the stability of the developed LMII and constrained LMII methods, we add Gaussian noise with a signal-to-noise ratio (SNR) of 3 (whose mean is 0 and standard deviation is one-third of the root mean square of the signal) to Figure 2b and obtain a noisy seismic profile for impedance inversion (Figure 5). In this case, the reference model and the setting of measurement ratios are the same as in the noise-free experiment. Figure 6a–c presents the results inverted by MII, LMII, and constrained LMII methods, respectively. Figure 6d–f presents the residuals between these inversion results and the true impedance, respectively. Note that the maximum number of iterations and tolerance error for the above methods are set to 500 and 10⁻³. As observed from Figure 6, all inversion methods can recover impedance from noise-contaminated seismic data. Nevertheless, the constrained LMII method is significantly better than the other two methods in terms of accuracy since the constrained LMII method achieves a smaller residual. For better comparison, a zoomed window of the true impedance and the inverted results of Figure 6 are shown in Figure 7. At the same time, Figure 8 compares the impedance curves of the 300th trace of the inversion results shown in Figure 6. It can be seen from Figure 7 and Figure 8 that the impedance interfaces of the MII and LMII methods are blurred, while the pseudo-layer interferences are obvious (as shown by arrows). In contrast, the impedance inverted by the constrained LMII method is much better than the other two methods, where the prediction errors are lower, and the layer interfaces are clear. To quantitatively evaluate the inversion quality, we calculate the RMSE values of the three inversion results and show them in

Table 2. Quantitative evaluation index and running time of the inversion results of Figure 6.

Method	RMSE (×10−3)	Time (s)
MII	79.93	269.31
LMII	79.94	138.28
Constrained LMII	57.92	165.80

. Also, the elapsed time of the above three inversion methods is listed in Table 2. From Table 2, it can be found that the inversion accuracy of the LMII method is essentially the same as that of the traditional MII methods, while the elapsed time is reduced by nearly half. This further illustrates the advantage of the LMII method in improving inversion efficiency. Besides, it can also be found from Table 2 that the constrained LMII method has a smaller RMSE value than the other two methods, which again effectively proves that the constrained LMII method is effective in improving the inversion quality. At the same time, this result also validates that the constrained LMII method has better robustness under noise circumstances.

3.2. Field Example

In this subsection, we apply the developed algorithms to the field data to test their practicability for real data inversion. Figure 9a shows 2D field post-stack data, which consists of 901 traces, each containing 300 sampling points with an interval of 1 ms. Therein, a blind well at the 288th trace (shown as the gray line) to judge the reliability of the predicted result. Note that this blind well does not participate in the inversion procedure. Figure 9b displays the low-frequency reference model used for impedance inversion, which is derived by interpolation from the 12–14 Hz high-cut filtered log data within the work area. Figure 10a–c presents the results predicted by the MII, LMII, and constrained LMII methods, respectively. In this case, the maximum number of iterations and the tolerance error of the above methods are set to 500 and 10⁻⁵, respectively. In addition, the measurement ratios of the seismic record and reference model are set to 0.25 and 0.1 in the LMII and constrained LMII methods, respectively. It is evident from Figure 10 that the inversion results obtained by all three methods well reflect the lateral extension and lithological variation patterns of the geological bodies. Nevertheless, one notices that the constrained LMII method performs better than the other two methods in identifying the stratigraphic boundaries. Specifically, the results inverted by the constrained LMII method (Figure 10c) have obvious blocky characteristics and show clear geological boundaries, while the results inverted by the other two are over-smooth, blurring the lithological boundary features and leading to difficulties in layer identification. Additionally, the conventional MII method and the LMII method also produce additional pseudo-layer interferences (as shown in the elliptical zone). To better illustrate this point, Figure 11a–c presents the measured well logging curve at the 288th trace and the predicted impedance by the above three inversion methods, respectively. Herein, the green line represents the low-frequency reference curve, the black line represents the logging data, and the red line represents the predicted impedance curve. It is clear from Figure 11 that these inversion results show similar trends to the well log, but the inverted curve by the constrained LMII approach is more consistent with the well log (as shown by the arrows). In addition, as seen in the logging curve, a low-impedance shale layer exists near 3.23 s. After inversion by the conventional MII method and the LMII method, pseudo-layer interferences are generated in their results, which led to the appearance of false sand-shale interbeds. In contrast, the constrained LMII method effectively avoids the pseudo-layer interferences and correctly identifies this stratum. It reveals that the constrained LMII method outperforms the other two methods in terms of inversion reliability, which is consistent with the synthetic examples. In this case, the running times of the MII method, the LMII method, and the constrained LMII method are 298.74, 138.19, and 152.21 s, respectively. Clearly, both LMII and constrained LMII methods have smaller computing costs compared with the MII method, which further demonstrates that it is effective to compress the size of the inverse problem by low-dimensional measurements in sparse space.

4. Discussion

Investigating the sensitivity of the impedance inversion algorithm to the estimation accuracy of the source wavelet has important guiding significance for its subsequent practical application. Herein, we take Figure 2 as an example to study the sensitivity of the proposed methods to the source wavelet. Figure 12 shows the relationship curves between the quantitative evaluation index (RMSE) of inversion results by the proposed algorithms and wavelet estimation accuracy. As shown in Figure 12, the smallest RMSE value is obtained when the exact wavelet (i.e., frequency of 30 Hz) is used as the input of the LMII and constrained LMII algorithms. This means that both LMII and constrained LMII algorithms achieve the highest inversion accuracy at this time. With the decrease in the wavelet estimation accuracy, the RMSE value of inversion results yielded by the LMII and constrained LMII methods gradually increases. This result demonstrates that the performance of the LMII and constrained LMII methods has a strong dependence on the estimation accuracy of the source wavelet. Therefore, in order to achieve a satisfactory inversion result, the exact source wavelet must be extracted before applying the LMII and constrained LMII methods for impedance inversion.

Although the developed method effectively improves the computational accuracy and efficiency of the multi-trace impedance inversion, it is not flawless and is prone to suffering from lateral smearing when dealing with highly dipping structures. The reason is that the lateral constraint term constructed by the second-order difference operator only forces lateral smoothing and ignores the morphology of underground structures [26]. Recently, structure-oriented multi-trace impedance inversion methods [26,30,31] have effectively solved the lateral smearing problem and can be better adapted to inversion work in regions with large dips. However, similar to the traditional multi-trace impedance method, these methods also suffer from insufficient inversion efficiency and accuracy. Fortunately, the inversion framework presented in this paper has excellent expansibility. Therefore, combining the proposed inversion framework with structure-oriented multi-trace inversion methods to develop a series of novel impedance inversion algorithms that can adapt to highly dipping structures with high computational efficiency and accuracy will be a future research direction.

5. Conclusions

To address the issues existing in the conventional MII method with low efficiency and accuracy, a low-dimensional multi-trace inversion method with elastic half norm constraints is proposed. First, we develop an LMII inversion framework that efficiently reduces the scale of the MII inversion problem and thus achieves efficient inversion. Concretely, this framework converts the seismic data and reference impedance into the sparse space by an orthogonal transformation and constructs the objective function with a smaller size for multi-trace inversion using the low-dimensional measurements in the sparse domain. Subsequently, on the basis of the developed LMII framework, we introduce the EH norm as a vertical constraint on the model parameters and propose a novel constrained LMII model for inverting the subsurface impedance. Since this inversion procedure considers both the smoothness and blockiness of the subsurface strata, it can significantly improve the inversion precision compared with the traditional MII method, which considers smoothness only. Finally, based on the variable splitting technique, we derive an alternating multiplier iteration algorithm to efficiently solve the constrained LMII problem. Two experiments are conducted to investigate the performance of the developed algorithms. The results indicate that the LMII framework can effectively improve the multi-trace inversion efficiency while keeping the inversion accuracy constant. Meanwhile, it is also shown that the constrained LMII model can evidently promote the quality of impedance inversion while inheriting the efficiency of the LMII framework. In addition, it should be pointed out that the proposed methods have good scalability and can be extended to gravity inversion, magnetic inversion, and other fields to provide new insights for solving more geophysical problems.