Journal of Earth Science  2018, Vol. 29 Issue (6): 1359-1371 PDF     0
Prestack Multi-Gather Simultaneous Inversion of Elastic Parameters Using Multiple Regularization Constraints
Shu Li1,2,3, Zhenming Peng1,2, Hao Wu1,2
1. School of Optoelectronic Information, University of Electronic Science and Technology of China, Chengdu 610054, China;
2. Center for Information Geoscience, University of Electronic Science and Technology of China, Chengdu 610054, China;
3. School of Information Science and Engineering, Jishou University, Jishou 416000, China
ABSTRACT: Inversion of Young's modulus, Poisson's ratio and density from pre-stack seismic data has been proved to be feasible and effective. However, the existing methods do not take full advantage of the prior information. Without considering the lateral continuity of the inversion results, these methods need to invert the reflectivity first. In this paper, we propose multi-gather simultaneous inversion for pre-stack seismic data. Meanwhile, the total variation (TV) regularization, L1 norm regularization and initial model constraint are used. In order to solve the objective function contains L1 norm, TV norm and L2 norm, we develop an algorithm based on split Bregman iteration. The main advantages of our method are as follows:(1) The elastic parameters are calculated directly from objective function rather than from their reflectivity, therefore the stability and accuracy of the inversion process can be ensured. (2) The inversion results are more in accordance with the prior geological information. (3) The lateral continuity of the inversion results are improved. The proposed method is illustrated by theoretical model data and experimented with a 2-D field data.
KEY WORDS: elastic parameter    pre-stack inversion    multi-gather    regularization

0 INTRODUCTION

In recent years the shale gas prospecting has gradually become a research focus of unconventional oil and gas exploration. Several scholars have done significant works on shale gas seismic exploration. Rickman et al. (2008) pointed out that by using pre-stack seismic inversion we can obtain Poisson's ratio (PR) and Young's modulus (YM) from pre-stack seismic data. These two parameters are very useful in determining the brittleness of rocks and valuable for identifying shale gas reservoirs. In general, shale gas reservoirs have low PR and high YM (Harris et al., 2011; Sena et al., 2011). Inspired by previous studies, Zong et al. (2012) derived the mathematical relationship between the P wave reflection coefficient and YM, PR and density (referred as YPD equation). Based on YPD equation, the relationship between pre-stack seismic data and these three elastic parameters can be established directly. Their work provides a feasible method for direct inversion of these three parameters. In pre-stack inversion, density term often has adverse effects on the inversion results (Zhang et al., 2015). To eliminate the influence of density term in the YPD equation on the inversion results, Zong et al. (2013) derived a mathematical equation describing the relationship between the elastic impedance and the parameters YM and PR. They proposed a method to invert elastic impedance, and then YM and PR are retrieved from elastic impedance. Their experiments showed that this method can get more stable PR and YM.

Seismic data usually have low signal-to-noise ratio. In order to improve inversion results, some scholars have introduced the L1 norm regularization (Baraniuk, 2007; Donoho, 2006) technique into seismic inversion, and have gotten a better result. To improve the resolution ability of thin layer, Zhang et al. (2013) proposed a basis pursuit inversion (BPI) algorithm for pre-stack seismic inversion. BPI algorithm is based on the Basis pursuit (Chen et al., 2001) which is an L1 norm optimization algorithm. In their experiments, the strata boundaries have been delineated very clearly by using BPI algorithm. Because of the advantages of BPI in pre-stack inversion, Yin et al. (2015) proposed a BPI based inversion method to get the "fracability" of shale. Their inversion method consists of two successive steps (referred as two-step inversion): Firstly, imposing the L1 norm sparse constraint on the reflectivity of YM, PR and density, and solving the L1 norm optimization problem to obtain the sparse reflectivity. Secondly, obtaining YM, PR and density from their reflectivity.

In seismic inversion, L1 norm regularization technique is generally used to get the sparse solution in depth/time direction, without considering the lateral continuity. To reinforce lateral continuity, multi-trace/multichannel simultaneous inversion is used in post-stack seismic inversion. Hamid and Pidlisecky (2015) introduced lateral constraint into seismic impedance inversion and presented a multi-trace impedance inversion method. In their experiments the lateral continuity is enhanced by using multi-trace simultaneous inversion. Yuan et al. (2015) showed that multi-trace simultaneous inversion is useful for stabilizing inversion process and utilizing the lateral continuity of underground structures. In addition, obtaining blocky solutions is also a quite meaningful goal for seismic inversion. In blocky image, the boundaries between different strata have been enhanced, so geologists can easily interpret the geological structures and thin layers can be detected more easily (Pérez et al., 2014). To get blocky inversion results, Zhang et al. (2014) took the total variation (TV) (Rudin et al., 1992) regularization into seismic impedance inversion. Gholami (2015) used the TV regularization and multi-trace simultaneous inversion to derive blocky structures, and got better lateral continuity. Their experiments showed that we can get blocky impedance structures by using TV regularization.

From the above introduction, we can discover that obtaining YM, PR and density by employing pre-stack seismic inversion is possible and useful. Zong et al. (2012) have derived the forward equation and proposed the inversion algorithm. To polish up the inversion results, Yin et al. (2015) introduced the BPI technique. Their endeavors lay a foundation for inverting YM and PR from pre-stack data. However, there are still some works that should be done to improve the inversion results. The defects of the previous researches are as following: (1) No particular measures are taken to enhance the lateral continuity of the inversion results, (2) the prior information hasn't been used sufficiently due to fewer regula- rization constraints, (3) the two-step inversion methods will produce cumulative errors. The third disadvantage can use the formula $\mathit{\boldsymbol{E}}(t) = \mathit{\boldsymbol{E}}({t_0})\exp \left[ {2\int_{{t_0}}^t {{\mathit{\boldsymbol{R}}_E}(\eta)d\eta } } \right]$ (Yin et al., 2015) to explain. It's a formula used to obtain YM from its reflectivity. Where E(t) is the YM at time=t, RE is the reflectivity of YM, and E(t0) is the YM at time t=0. We can find that YM depends on its value at time t=0 and the reflectivity inversion result. If E(t0) or RE is inaccurate, E(t) will not be accurate. Actually, E(t0) is usually very difficult to be known precisely. Furthermore, there is a difference between the reflectivity inversion result and the true reflectivity. Meanwhile, E(t) is obtained from its reflectivity by integrating, which can result in accumulating the error of the reflectivity at the present time to E(t) at subsequent times.

Motivated by the multi-trace/multichannel simultaneous inversion used in post-stack seismic inversion, we propose multi-gather simultaneous inversion for pre-stack seismic data. In seismic impedance inversion, blocky impedance structures can be obtained by using TV regularization (Gholami, 2015; Zhang et al., 2014), the anti-noise ability of inversion methods can be improved by exploiting L1 norm regularization, and the non-uniqueness can be reduced through adding initial model constraint. They inspired us to impose these regularization constraints on pre-stack seismic inversion. In order to reduce the cumulative errors caused by two-step inversion, based on the split Bregman iterative algorithm (Goldstein and Osher, 2009), we developed an algorithm to obtain the YM, PR and density directly, while we don't need to invert their reflectivity first. Our method has advantages in stability, accuracy, anti- noise ability and lateral continuity.

1 METHODOLOGY 1.1 Forward Equation

YM and PR are key elastic parameters in the seismic exploration of shale gas, which can be used to measure the brittleness of rocks. Due to the distinguished research of Zong et al. (2012), the relationship between the pre-stack seismic data, YM, PR and density can be obtained by exploiting the YPD approximation equation. In our work, we used the YPD equation to establish the forward equation. The YPD equation can be expressed as

 $\begin{array}{l} \mathit{\boldsymbol{R}}(\theta) = (\frac{1}{2}{\sec ^2}\theta - 4k{\sin ^2}\theta)\frac{{\Delta \mathit{\boldsymbol{E}}}}{{2\mathit{\boldsymbol{E}}}} + (\frac{1}{2}{\sec ^2}\theta \frac{{(2k - 3){{(2k - 1)}^2}}}{{k(4k - 3)}} + \\ \quad \quad \quad 4k{\sin ^2}\theta \frac{{1 - 2k}}{{3 - 4k}})\frac{{\Delta \mathit{\boldsymbol{\sigma }}}}{{2\mathit{\boldsymbol{\sigma }}}}\; + (1 - \frac{1}{2}{\sec ^2}\theta)\frac{{\Delta \mathit{\boldsymbol{\rho }}}}{{2\mathit{\boldsymbol{\rho }}}} \end{array}$ (1)

where θ represents the incident angle; R(θ) indicates the reflectivity; △E/2E, △σ/2σ and △ρ/2ρ are the reflectivity of YM, PR and density, respectively; k is the square of the average velocity ratio of S-to-P waves. If we let

 $\begin{array}{l} {\mathit{\boldsymbol{R}}_E} = \Delta \mathit{\boldsymbol{E}}/2\mathit{\boldsymbol{E}}, \;{\mathit{\boldsymbol{R}}_\sigma } = \Delta \sigma /2\sigma, \;{\mathit{\boldsymbol{R}}_\rho } = \Delta \mathit{\boldsymbol{\rho }}/2\mathit{\boldsymbol{\rho }}\\ {\mathit{\boldsymbol{a}}_\theta } = \frac{1}{2}{\sec ^2}\theta - 4k{\sin ^2}\theta, \;{\mathit{\boldsymbol{c}}_\theta } = 1 - \frac{1}{2}{\sec ^2}\theta \\ {\mathit{\boldsymbol{b}}_\theta } = \frac{1}{2}{\sec ^2}\theta \frac{{(2k - 3){{(2k - 1)}^2}}}{{k(4k - 3)}} + 4k{\sin ^2}\theta \frac{{1 - 2k}}{{3 - 4k}} \end{array}$

then (1) can be rewritten as

 $\mathit{\boldsymbol{R}}(\theta) = {\mathit{\boldsymbol{a}}_\theta }{\mathit{\boldsymbol{R}}_E} + {\mathit{\boldsymbol{b}}_\theta }{\mathit{\boldsymbol{R}}_\sigma } + {\mathit{\boldsymbol{c}}_\theta }{\mathit{\boldsymbol{R}}_\rho }$ (2)

In seismic inversion, if we assume the elastic parameter δ(t) changes continuously over time (Russell and Dommico, 1988) and its reflectivity r(t) is less than 0.3, then the relationship between r(t) and δ(t) can be expressed as Eq. (3) (Walker and Ulrych, 1983)

 $\left[ {\begin{array}{*{20}{c}} {r(1)}\\ {r(2)}\\ \vdots \\ {r(n - 1)} \end{array}} \right] = \frac{1}{2}\left[ {\begin{array}{*{20}{c}} { - 1}&1&{}&{}&{}\\ {}&{ - 1}&1&{}&{}\\ {}&{}& \ddots & \ddots &{}\\ {}&{}&{}&{ - 1}&1 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{\rm{ln}}\;(\delta (1))}\\ {{\rm{ln}}\;(\delta (2))}\\ \vdots \\ {{\rm{ln}}\;(\delta (n))} \end{array}} \right]$ (3)

Yilmaz (2001) pointed out that a representative reflectivity is about 0.2 for a strong reflecting layer. In pre-stack inversion, the PR, density and YM are closely related to the physical properties of the rock itself. If the physical properties of the rock only have a slight relative change in the stratigraphic boundaries, these elastic parameters and their reflectivity will meet Eq. (3).

Let D be the matrix

 $\mathit{\boldsymbol{D}} = \frac{1}{2}\left[ {\begin{array}{*{20}{c}} { - 1}&1&{}&{}&{}\\ {}&{ - 1}&1&{}&{}\\ {}&{}& \ddots & \ddots &{}\\ {}&{}&{}&{ - 1}&1 \end{array}} \right]$ (4)

Let

 $\begin{array}{l} \mathit{\boldsymbol{r}} = {\left[ {\begin{array}{*{20}{c}} {r(1)}&{r(2)}& \cdots &{r(n - 1)} \end{array}} \right]^T}\\ {\mathit{\boldsymbol{L}}_\delta } = {\left[ {\begin{array}{*{20}{c}} {{\rm{ln}}\;(\delta (1))}&{{\rm{ln}}\;(\delta (2))}& \cdots &{{\rm{ln}}\;(\delta (n))} \end{array}} \right]^T} \end{array}$

where n denotes the number of sampling points of δ(t), T represents the transpose, then (3) can be written in the matrix equation

 $\mathit{\boldsymbol{r}} = \mathit{\boldsymbol{D}}{\mathit{\boldsymbol{L}}_\delta }$ (5)

Following a similar derivation to the Eq. (5), the following equations can be obtained.

 ${\mathit{\boldsymbol{R}}_E} = \mathit{\boldsymbol{D}}{\mathit{\boldsymbol{L}}_E}, \;{\mathit{\boldsymbol{R}}_\sigma } = \mathit{\boldsymbol{D}}{\mathit{\boldsymbol{L}}_\sigma }, \;{\mathit{\boldsymbol{R}}_\rho } = \mathit{\boldsymbol{D}}{\mathit{\boldsymbol{L}}_\rho }$

where LE = ln E, Lσ = ln σ, and Lρ = ln ρ are the logarithm of YM, PR and density, respectively. Taking these three equations into (2), we can derive the reflectivity equation in terms of LE, Lσ and Lρ, which can be written as

 $\mathit{\boldsymbol{R}}(\theta) = {\mathit{\boldsymbol{a}}_\theta }\mathit{\boldsymbol{D}}{\mathit{\boldsymbol{L}}_E} + {\mathit{\boldsymbol{b}}_\theta }\mathit{\boldsymbol{D}}{\mathit{\boldsymbol{L}}_\sigma } + {\mathit{\boldsymbol{c}}_\theta }\mathit{\boldsymbol{D}}{\mathit{\boldsymbol{L}}_\rho }$ (6)

Combining the reflectivity equation and the convolution model, we can get the forward equation

 $\mathit{\boldsymbol{S}} = \left[ {\begin{array}{*{20}{c}} {{S_{\theta 1}}}\\ {{S_{\theta 2}}}\\ \vdots \\ {{S_{\theta N}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{a_{\theta 1}}{W_{\theta 1}}\mathit{\boldsymbol{D}}}&{{b_{\theta 1}}{W_{\theta 1}}\mathit{\boldsymbol{D}}}&{{c_{\theta 1}}{W_{\theta 1}}\mathit{\boldsymbol{D}}}\\ {{a_{\theta 2}}{W_{\theta 2}}\mathit{\boldsymbol{D}}}&{{b_{\theta 2}}{W_{\theta 2}}\mathit{\boldsymbol{D}}}&{{c_{\theta 2}}{W_{\theta 2}}\mathit{\boldsymbol{D}}}\\ \vdots & \vdots & \vdots \\ {{a_{\theta N}}{W_{\theta N}}\mathit{\boldsymbol{D}}}&{{b_{\theta N}}{W_{\theta N}}\mathit{\boldsymbol{D}}}&{{c_{\theta N}}{W_{\theta N}}\mathit{\boldsymbol{D}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{L}}_E}}\\ {{\mathit{\boldsymbol{L}}_\sigma }}\\ {{\mathit{\boldsymbol{L}}_\rho }} \end{array}} \right]$ (7)

where the subscript θi(i=1, 2…N) represents the ith incident angle, Wθi is the wavelet matrix of the ith incident angle. Let L=(LE, Lσ, Lρ)T, and let A be the form of (8)

 $\mathit{\boldsymbol{A}} = \left[ {\begin{array}{*{20}{c}} {{a_{\theta 1}}{W_{\theta 1}}\mathit{\boldsymbol{D}}}&{{b_{\theta 1}}{W_{\theta 1}}\mathit{\boldsymbol{D}}}&{{c_{\theta 1}}{W_{\theta 1}}\mathit{\boldsymbol{D}}}\\ {{a_{\theta 2}}{W_{\theta 2}}\mathit{\boldsymbol{D}}}&{{b_{\theta 2}}{W_{\theta 2}}\mathit{\boldsymbol{D}}}&{{c_{\theta 2}}{W_{\theta 2}}\mathit{\boldsymbol{D}}}\\ \vdots & \vdots & \vdots \\ {{a_{\theta N}}{W_{\theta N}}\mathit{\boldsymbol{D}}}&{{b_{\theta N}}{W_{\theta N}}\mathit{\boldsymbol{D}}}&{{c_{\theta N}}{W_{\theta N}}\mathit{\boldsymbol{D}}} \end{array}} \right]$ (8)

then Eq. (7) can be rewritten as

 $\mathit{\boldsymbol{S}} = \mathit{\boldsymbol{AL}}$ (9)
1.2 Our Method

Our aim is to obtain accurate and stable inversion results of YM, PR and density, and to improve their lateral continuity. In impedance inversion, multi-trace inversion can be used to explore spatial continuities of strata, which enlightened us to use it in the pre-stack seismic inversion. However, the pre-stack inversion is different from post-stack inversion. In post-stack inversion, each seismic trace corresponds to a trace of inversion result. While, in pre-stack inversion, multiple seismic traces from a gather corresponds to a trace of elastic parameter. In order to improve the lateral continuity of the pre-stack inversion results, we extend the multi-trace inversion in post-stack inversion to the multi-gather simultaneous inversion in the pre-stack inversion.

By using the multi-gather simultaneous inversion, the initial objective function is established by following steps. In our method, Si (i = 1, 2 … N) represents the ith pre-stack gather, which can be seen from Fig. 1.

We arrange each gather to a column vector. ${\mathit{\boldsymbol{\tilde S}}_i}$ is a column vector obtained from the gather Si. This process can be expressed as ${\mathit{\boldsymbol{\tilde S}}_i} = {\rm{vec}}({\mathit{\boldsymbol{S}}_i})$. vec() represents an operation to arrange a pre-stack gather into a column vector. We get the data vector d through arrange ${\mathit{\boldsymbol{\tilde S}}_i}(i = 1, \; 2, \; \cdots N)$ into a long column vector. This can be expressed as $\mathit{\boldsymbol{d}} = {\left[ {\begin{array}{*{20}{c}} {{{\mathit{\boldsymbol{\tilde S}}}_1}} & {{{\mathit{\boldsymbol{\tilde S}}}_2}} & \cdots & {{{\mathit{\boldsymbol{\tilde S}}}_N}} \end{array}} \right]^{\rm{T}}}.$ Similarly, to meet the requirement of multi-gather simultaneous inversion, we should arrange the logarithm of YM, PR and density to a long column vector. Let ${\mathit{\boldsymbol{L}}_i}(i = 1, \; 2, \; \cdots N)$ be the logarithm of YM, PR and density about ith gather, Li can be written as ${\mathit{\boldsymbol{L}}_i} = {\left[ {{{\rm{L}}_{i, \; E}}\; {\mathit{\boldsymbol{L}}_{i, \; \sigma }}\; {\mathit{\boldsymbol{L}}_{i, \; \rho }}} \right]^{\rm{T}}}$. Then the model parameter vector m can be written as $\mathit{\boldsymbol{m}} = {\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{L}}_1}} & {{\mathit{\boldsymbol{L}}_2}} & \cdots & {{\mathit{\boldsymbol{L}}_N}} \end{array}} \right]^{\rm{T}}}$.

We also need to change the linear operator of forward equation accordingly, which is the matrix A in Eq. (9). Note that, if we invert multiple gathers simultaneously, the linear operator needs to be extended. Let G be the new linear operator, which can be expressed as G=kron(I, A) where kron represents the Kronecker product, I is an identity matrix, A is defined in Eq. (8). We also can write G as

 $\mathit{\boldsymbol{G}} = \left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{A}}&{}&{}&{}\\ {}&\mathit{\boldsymbol{A}}&{}&{}\\ {}&{}& \ddots &{}\\ {}&{}&{}&\mathit{\boldsymbol{A}} \end{array}} \right]$ (10)

Now, we can derive the initial objective function as

 $J(\mathit{\boldsymbol{m}}) = \mathop {\min }\limits_{\boldsymbol{m}} \frac{\eta }{2}\left\| {\mathit{\boldsymbol{Gm}} - \mathit{\boldsymbol{d}}} \right\|_2^2$ (11)

where η is the regularization parameter of data fidelity term.

L1 norm constraint is used to enhance the anti-noise performance. Based on the initial objective function, we impose the L1 norm constraint on the reflectivity, and then the objective function becomes the following form

 $J(\mathit{\boldsymbol{m}}) = \mathop {\min }\limits_{\bf{m}} \frac{\eta }{2}\left\| {\mathit{\boldsymbol{Gm}} - \mathit{\boldsymbol{d}}} \right\|_2^2 + {\left\| \mathit{\boldsymbol{R}} \right\|_1}\;\;s.t.\;\;\mathit{\boldsymbol{R}} = \mathit{\boldsymbol{\tilde Dm}}$ (12)

where $\mathit{\boldsymbol{\tilde D}} = {\rm{kron}}(\mathit{\boldsymbol{I}}, \mathit{\boldsymbol{D}})$ is the Kronecker product of identity matrix I and matrix D defined in (4). Reflectivity R also need to meet the requirements of multi-gather simultaneous inversion. Let ${\mathit{\boldsymbol{R}}_i}(i = 1, \; 2, \cdots N)$ be the reflectivity of elastic parameters about ith gather, which is written as ${\mathit{\boldsymbol{R}}_i} = {[{r_{i, E}}\; {r_{i, \sigma }}\; {r_{i, \rho }}]^{\rm{T}}}.$ Then R is gotten by arranging ${\mathit{\boldsymbol{R}}_i}(i = 1, \; 2, \cdots N)$ into a column vector. It is expressed as $\mathit{\boldsymbol{R}} = {\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{R}}_1}} & {{\mathit{\boldsymbol{R}}_2}} & \cdots & {{\mathit{\boldsymbol{R}}_{\rm{N}}}} \end{array}} \right]^{\rm{T}}}.$ Now, we can use Lagrangian multiplier method to convert Eq. (12) into an unconstrained optimization problem.

 $J(\mathit{\boldsymbol{m}}) = \mathop {\min }\limits_\mathit{\boldsymbol{m}} \frac{\eta }{{\rm{2}}}\left\| {\mathit{\boldsymbol{Gm}} - \mathit{\boldsymbol{d}}} \right\|_{\rm{2}}^{\rm{2}} + {\left\| \mathit{\boldsymbol{R}} \right\|_1} + \frac{\gamma }{2}\left\| {\mathit{\boldsymbol{R}} - \mathit{\boldsymbol{\tilde Dm}}} \right\|_2^2$ (13)

where γ is a Lagrangian multiplier.

Another important goal of this article is to obtain the blocky inversion results. When we introduce the TV constraint into the objective function, Eq. (13) becomes the following form

 $\begin{array}{l} J(\mathit{\boldsymbol{m}}) = \mathop {\min }\limits_\mathit{\boldsymbol{m}} \frac{\eta }{2}\left\| {\mathit{\boldsymbol{Gm}}{\bf{ - }}\mathit{\boldsymbol{d}}} \right\|_2^2 + {\left\| \mathit{\boldsymbol{R}} \right\|_1} + \frac{\gamma }{2}\left\| {\mathit{\boldsymbol{R}} - \mathit{\boldsymbol{\tilde Dm}}} \right\|_2^2 + \\ \quad \quad \quad \frac{\varphi }{2}\sum\limits_{i, j} {\sqrt {({\nabla _x}\mathit{\boldsymbol{m}})_{i, j}^2 + ({\nabla _y}\mathit{\boldsymbol{m}})_{i, j}^2} } \end{array}$ (14)

where φ is the regularization parameter of TV term and $\sum\limits_{i, j} {\sqrt {({\nabla _x}\mathit{\boldsymbol{m}})_{i, j}^2 + ({\nabla _y}\mathit{\boldsymbol{m}})_{i, j}^2} }$ is the TV of m. For simplicity, the TV of m can be written as the following short-hand notation

 ${\left\| {({d_{x, }}\;{d_y})} \right\|_2} = \sum\limits_{i, j} {\sqrt {({\nabla _x}\mathit{\boldsymbol{m}})_{i, j}^2 + ({\nabla _y}\mathit{\boldsymbol{m}})_{i, j}^2} }$ (15)

where ${d_x} = {\nabla _x}\mathit{\boldsymbol{m}}$ and ${d_y} = {\nabla _y}\mathit{\boldsymbol{m}}$ represent the difference of m in x and y direction, respectively. Now, (14) can be rewritten as

 $\begin{array}{l} J(\mathit{\boldsymbol{m}}) = \mathop {\min }\limits_{\bf{m}} \frac{\eta }{2}\left\| {\mathit{\boldsymbol{Gm}}{\bf{ - }}\mathit{\boldsymbol{d}}} \right\|_2^2 + {\left\| \mathit{\boldsymbol{R}} \right\|_1} + \frac{\gamma }{2}\left\| {\mathit{\boldsymbol{R}} - \mathit{\boldsymbol{\tilde Dm}}} \right\|_2^2 + \frac{\varphi }{2}{\left\| {({d_x}, \;{d_y})} \right\|_2}\;\\ \;s.t.\;{d_x} = {\nabla _x}\mathit{\boldsymbol{m}}, \;{d_y} = {\nabla _y}\mathit{\boldsymbol{m}} \end{array}$ (16)

Equation (16) is a constrained optimization problem, which can be transformed into an unconstrained optimization problem.

 $\begin{array}{l} J(\mathit{\boldsymbol{m}}) = \mathop {\min }\limits_{\bf{m}} \frac{\eta }{2}\left\| {\mathit{\boldsymbol{Gm}}{\bf{ - }}\mathit{\boldsymbol{d}}} \right\|_2^2 + {\left\| \mathit{\boldsymbol{R}} \right\|_1} + \frac{\gamma }{2}\left\| {\mathit{\boldsymbol{R}} - \mathit{\boldsymbol{\tilde Dm}}} \right\|_2^2 + \\ \quad \quad \quad \frac{\varphi }{2}{\left\| {({d_x}, {d_y})} \right\|_2}\; + \frac{\lambda }{2}\left\| {{d_x} - {\nabla _x}\mathit{\boldsymbol{m}}} \right\|_2^2 + \frac{\lambda }{2}\left\| {{d_y} - {\nabla _y}\mathit{\boldsymbol{m}}} \right\|_2^2 \end{array}$ (17)

where λ is a Lagrangian multiplier.

In order to stabilize the inverse operation and reduce the non-uniqueness of the inversion results, we add initial model constraint to the objective function. Equation (18) shows the final objective function

 $\begin{array}{l} J(\mathit{\boldsymbol{m}}) = \mathop {\min }\limits_{\bf{m}} \frac{\eta }{2}\left\| {\mathit{\boldsymbol{Gm}}{\bf{ - }}\mathit{\boldsymbol{d}}} \right\|_2^2 + {\left\| \mathit{\boldsymbol{R}} \right\|_1} + \frac{\gamma }{2}\left\| {\mathit{\boldsymbol{R}} - \mathit{\boldsymbol{\tilde Dm}}} \right\|_2^2 + \\ \quad \quad \;\;\;\;\frac{\varphi }{2}{\left\| {({d_x}, \;{d_y})} \right\|_2}\; + \frac{\lambda }{2}\left\| {{d_x} - {\nabla _x}\mathit{\boldsymbol{m}}} \right\|_2^2 + \frac{\lambda }{2}\left\| {{d_y} - {\nabla _y}\mathit{\boldsymbol{m}}} \right\|_2^2 + \\ \quad \quad \quad \frac{\alpha }{2}\left\| {\mathit{\boldsymbol{m}} - \mathit{\boldsymbol{\bar m}}} \right\|_2^2 \end{array}$ (18)

where α is the regularization parameter of initial model constraint term, and m represents the initial model vector. The structure of m is similar to that of m.m is composed of ${\mathit{\boldsymbol{\bar m}}_i}\; (i = 1, \; 2, \cdots N),$ which is the initial model vector of ith pre-stack gather. It can be written as $\mathit{\boldsymbol{\bar m}} = {\left[ {\begin{array}{*{20}{c}} {{{\mathit{\boldsymbol{\bar m}}}_1}} & {{{\mathit{\boldsymbol{\bar m}}}_2}} & \cdots & {{{\mathit{\boldsymbol{\bar m}}}_N}} \end{array}} \right]^T}$, and ${\mathit{\boldsymbol{\bar m}}_i}$ can be expressed as ${\mathit{\boldsymbol{\bar m}}_i} = {\left[ {\begin{array}{*{20}{c}} {{{\mathit{\boldsymbol{\bar L}}}_{i, E}}} & {{{\mathit{\boldsymbol{\bar L}}}_{i, \sigma }}} & {{{\mathit{\boldsymbol{\bar L}}}_{i, \rho }}} \end{array}} \right]^T}$. ${\mathit{\boldsymbol{\bar L}}_{i, x}}(x = E, \sigma, \rho)$ is the logarithm of the initial YM, PR and density model.

As can be seen from Eq. (18), the objective function has a complicated form with L1 norm, L2 norm and TV norm terms. Thanks to the work of Goldstein and Osher (2009), we can use the split Bregman iteration algorithm to solve the optimization problem with multiple regularization constraints. The main idea of this algorithm in solving this kind of problem is to "split" it into several simple problems. According to this algorithm, when we solve Eq. (18), we should take the auxiliary variable bx, by and bR into it. Then Eq. (18) can be expressed as

 $\begin{array}{l} J(\mathit{\boldsymbol{m}}) = \mathop {\min }\limits_{\bf{m}} \frac{\eta }{2}\left\| {\mathit{\boldsymbol{Gm}}{\bf{ - }}\mathit{\boldsymbol{d}}} \right\|_2^2 + {\left\| \mathit{\boldsymbol{R}} \right\|_1} + \frac{\gamma }{2}\left\| {\mathit{\boldsymbol{R}} - \mathit{\boldsymbol{\tilde Dm}} - {b_R}} \right\|_2^2 + \\ \quad \quad \quad \frac{\varphi }{2}{\left\| {({d_x}, \;{d_y})} \right\|_2} + \frac{\lambda }{2}\left\| {{d_x} - {\nabla _x}\mathit{\boldsymbol{m}} - {b_x}} \right\|_2^2 + \\ \quad \quad \quad \frac{\lambda }{2}\left\| {{d_y} - {\nabla _y}\mathit{\boldsymbol{m}} - {b_y}} \right\|_2^2 + \frac{\alpha }{2}\left\| {\mathit{\boldsymbol{m}} - \mathit{\boldsymbol{\bar m}}} \right\|_2^2 \end{array}$ (19)

It can be "decoupled" into the simpler sub problems about m, (dx, dy) and R. These sub problems can be written as following equations

 $\begin{array}{l} {\mathit{\boldsymbol{m}}^{k + 1}} = \mathop {\min }\limits_{\bf{m}} \frac{\eta }{2}\left\| {\mathit{\boldsymbol{Gm}} - \mathit{\boldsymbol{d}}} \right\|_2^2 + \frac{\gamma }{2}\left\| {{\mathit{\boldsymbol{R}}^k} - \mathit{\boldsymbol{\tilde Dm}} - b_R^k} \right\|_2^2 + \\ \quad \quad \;\;\frac{\lambda }{2}\left\| {d_x^k - {\nabla _x}\mathit{\boldsymbol{m}} - b_x^k} \right\|_2^2 + \frac{\lambda }{2}\left\| {d_y^k - {\nabla _y}\mathit{\boldsymbol{m}} - b_y^k} \right\|_2^2 + \\ \quad \;\;\;\;\;\frac{\alpha }{2}\left\| {\mathit{\boldsymbol{m}} - \mathit{\boldsymbol{\bar m}}} \right\|_2^2 \end{array}$ (20)
 $\begin{array}{l} (d_x^{k + 1}, \;d_y^{k + 1}) = \mathop {\min }\limits_{{d_x}, {d_y}} \frac{\varphi }{2}{\left\| {({d_x}, \;{d_y})} \right\|_2} + \\ \quad \quad \;\;\;\;\;\;\;\;\;\;\frac{\lambda }{2}\left\| {d_x^k - {\nabla _x}\mathit{\boldsymbol{m}} - b_x^k} \right\|_2^2 + \frac{\lambda }{2}\left\| {d_y^k - {\nabla _y}\mathit{\boldsymbol{m}} - b_y^k} \right\|_2^2 \end{array}$ (21)
 ${\mathit{\boldsymbol{R}}^{k + 1}} = {\rm{shrink}}(\mathit{\boldsymbol{\tilde D}}{\mathit{\boldsymbol{m}}^{k + 1}} + b_R^k, \frac{1}{\gamma })$ (22)

where k represents the kth iteration, and shrink() in Eq. (22) is a shrinkage function defined as

 ${\rm{shrink}}(x, \;\gamma) = \frac{x}{{\left| x \right|}}*{\rm{max}}(\left| x \right| - \gamma, \;0)$ (23)

Equation (21) can be further "decoupled" into the sub problems about dx and dy, which are expressed as

 $d_x^{k + 1} = \max ({h^{k + 1}} - \frac{\varphi }{{{\rm{2}}\lambda }}, 0)\frac{{{\nabla _x}{\mathit{\boldsymbol{m}}^{k + 1}} + b_x^k}}{{{h^{k + 1}}}}$ (24)
 $d_y^{k + 1} = \max ({h^{k + 1}} - \frac{\varphi }{{{\rm{2}}\lambda }}, 0)\frac{{{\nabla _y}{\mathit{\boldsymbol{m}}^{k + 1}} + b_y^k}}{{{h^{k + 1}}}}$ (25)

where

 ${h^{k + 1}} = \sqrt {{{\left| {{\nabla _x}{\mathit{\boldsymbol{m}}^{k + 1}} + b_x^k} \right|}^2} + {{\left| {{\nabla _y}{\mathit{\boldsymbol{m}}^{k + 1}} + b_y^k} \right|}^2}}$ (26)

We didn't introduce the details of the split Bregman iteration algorithm. The interested readers can refer to the work of Goldstein and Osher (2009).

It can be seen from Eqs. (20–26) that the sub problems are easy to deal with. Eq. (20) is the form of the sum of several L2 norm terms. It can be solved effectively by using conjugate gradient algorithm or stochastic conjugate gradient algorithm (Huang and Zhou, 2015). Equation (22) is a standard shrinkage formula used to obtain the sparse solution of the reflectivity R. Equations (24) and (25) are the generalized shrinkage formula used to get the optimal value of dx and dy. The pseudo code of our method is shown in algorithm 1 in the part of appendix, with the implementation details.

In algorithm 1, reshape() is a function used to reshape a vector into matrix. Step 14 extracts the inversion results from I. It is expressed as the form of the MATLAB statement, where N represents the sampling number of elastic parameter in time direction. It can be seen from algorithm 1 that, we do not need to know and also do not need to use the value of the elastic parameters in time t=0. The elastic parameters (E, σ, ρ) are obtained from the iterative process directly, and we do not need to invert their reflectivity first.

2 THEORETICAL MODEL EXPERIMENTS

We use a theoretical model to test the effectiveness of our method. The model contains horizontal, elliptic, wedge layers and thin layers between horizontal and elliptic layers as shown in Fig. 2. The sampling interval in time direction is 1 ms.

 Download: larger image Figure 2. True elastic parameters model. (a) YM, (b) PR, (c) density.

The true model in Fig. 2 and a Ricker wavelet with dominant frequency 35 Hz are used to generate the synthetic seismic data. This process is based on the forward Eq. (7). We have done the inversion experiment in the case of synthetic seismic data with 10% Gaussian random noise and 20% Gaussian random noise. The noisy synthetic seismic section is used as the original seismic record. Three gathers of the synthetic seismic data, one without noise, the other with 10% Gaussian random noise and the third with 20% Gaussian random noise are shown in Fig. 3. In each gather, the incident angles range from 0° to 35°.

 Download: larger image Figure 3. Gathers of the synthetic seismic data (a) without noise, (b) with 10% Gaussian random noise, (c) with 20% Gaussian random noise.

Figure 4 shows the initial elastic parameters model generated from the true model by Gauss low-pass filtering. It is considered as the initial guess of the inversion results. It can be seen that the initial model only retains the low frequency trend of the true model.

 Download: larger image Figure 4. Initial model. (a) YM, (b) PR, (c) density.

After adding 10% Gaussian random noise into the synthetic seismic data, a gather of the noisy seismic data is shown in Fig. 3b. In this case, some weak reflection layers become fuzzy. The inversion results are shown in Fig. 5. We find that the three elastic parameters are obtained clearly. The boundaries between different layers are very clear. Each of the strata can be accurately sketched out. Although the synthetic data is corrupted by the random noise, the effect of noise on the inversion results has been reduced significantly. In addition, the continuity in the horizontal direction is well maintained. Figure 6 shows the inversion results of a gather which extracted from Fig. 5. Blue dash-dot line, black dotted line and red solid line represent the true model, initial model and the inversion results, respectively. It can be seen that the inversion result of PR is very close to its true value. The inversion results of YM and density are deviated from their true value in some layers, but the error is relatively small. Moreover, in the layer boundaries, the inversion result curves are quite steep, which is advantageous in distinguishing different strata and improving the resolution ability.

 Download: larger image Figure 5. Inversion results when adding 10% Gaussian random noise. (a) YM, (b) PR, (c) density.

When we add 20% Gaussian random noise into the synthetic seismic data, some weak reflection layers become more difficult to identify. Figure 3c shows a gather of the noisy synthetic seismic data. In this case, the inversion results are shown in Fig. 7. It can be seen that the inversion results are more affected by noise due to an increase in noise amplitude. Although the inversion results have become worse, the layer boundaries are still clear, and we can distinguish the layers easily. In order to observe the inversion results of a single gather, we extract the inversion results of a gather from Fig. 7. They are shown in Fig. 8. We can find that the inversion results are worse than the situation of adding 10% Gaussian random noise. The result curves at the layer boundaries are not as steep as that in Fig. 6, and the line segments which represent layers are not very straight. However, the inversion result of PR is still very close to its true value. As for YM and density, the deviation between their inversion results and their true value is increased, but the relative size of the parameter values at different time is still correct. Therefore we can draw a conclusion that the inversion results are still satisfactory, even if the seismic data is corrupted by such a strong noise.

 Download: larger image Figure 6. Inversion results of a gather when adding 10% Gaussian random noise. (a) YM, (b) PR, (c) density.
 Download: larger image Figure 7. Inversion results when adding 20% Gaussian random noise. (a) YM, (b) PR, (c) density.
 Download: larger image Figure 8. Inversion results of a gather when adding 20% Gaussian random noise. (a) YM, (b) PR, (c) density.

Our method is a kind of model-based inversion method. In order to test the influence of the initial model on the inversion results, the initial model is set to be constant parameter and repeat the above experiment. The value of the initial YM model, PR model and density model are set to be 2.11 107(m/s)2×(g/cc), 0.29 and 2.1 g/cc, respectively. They are equal to the mean value of their true model. When we use the constant initial model and add 10% and 20% Gaussian random noise into the synthetic seismic data, the inversion results are obtained and displayed in Figs. 9 and 11. Inversion results of a gather from Figs. 9 and 11 are shown in Figs. 10 and 12.

 Download: larger image Figure 9. Inversion results when adding 10% Gaussian random noise and the initial model is set to be constant parameter. (a) YM, (b) PR, (c) density.
 Download: larger image Figure 10. Inversion results of a gather when adding 10% Gaussian random noise and the initial model is set to be constant parameter. (a) YM, (b) PR, (c) density.
 Download: larger image Figure 11. Inversion results when adding 20% Gaussian random noise and the initial model is set to be constant parameter. (a) YM, (b) PR, (c) density.
 Download: larger image Figure 12. Inversion results of a gather when adding 20% Gaussian random noise and the initial model is set to be constant parameter. (a) YM, (b) PR, (c) density.

By comparing Figs. 58 and Figs. 912, we can see that the difference between the inversion results obtained under different initial model is not significant. In other words, in our method the initial model has little effect on the inversion results. In fact, the initial model is used as the initial guess of the inversion result. In addition to the initial model, the inversion results also depends on other constraints (such as the data fidelity term) and inversion methods.

In the above experiments, we add the Gaussian random noise to the synthetic data. However, the noise in field data often has a strong lateral continuity. In order to test the effect of this kind of noise on the inversion results, we add the noise into the back-diagonal of the synthetic seismic section. Since it's not convenient to display the noisy synthetic pre-stack data, we extract a seismic trace from each gather and arrange them into sections, which are shown in Fig. 13. The amplitude of the noise added is equal to 30% and 50% of the maximum seismic amplitude.

 Download: larger image Figure 13. Seismic sections with noise in the back-diagonal. The sections are obtained from the prestack gathers by extracting a trace from each gather. The amplitude of the noise added is equal to 30% (a) and 50% (b) of the maximum seismic amplitude.

When the amplitude of the noise added in the back- diagonal of the seismic section is equal to 30% of the maximum seismic amplitude, the inversion results of YM, PR and density are shown in Fig. 14. It can be seen that the influence of the noise on the inversion results is very small. This is mainly due to the use of multi-gather simultaneous inversion, which can reduce the impact of the noise that has a strong lateral continuity. Figure 15 shows the inversion results under the condition of the amplitude of the noise added being equal to 50% of the maximum seismic amplitude. In this situation, the effect of the noise on the inversion results becomes more obvious. Therefore, we can find that our method can reduce the influence of the noise which has a strong lateral continuity as long as its amplitude is not very large (no more than 30% of the maximum seismic amplitude).

 Download: larger image Figure 14. Inversion results when the amplitude of the noise added in the back-diagonal is equal to 30% of the maximum seismic amplitude.
 Download: larger image Figure 15. Inversion results when the amplitude of the noise added in the back-diagonal is equal to 50% of the maximum seismic amplitude.
3 FIELD DATA EXPERIMENTS

The theoretical model inversion experiments prove that our method is feasible in theory. But we still need to test whether it is applicable to field data. The field data used in our experiment are acquired in China, which consist of 300 CDP gathers, and sampling interval in time domain is 1 ms. We extract a wavelet from each seismic trace, and the average of these wavelets is used for the inversion. The partial angle stack seismic sections used are shown in Figs. 16a16c. The mean angles of the small incident angle range, medium incident angle range and large incident angle range are 7°, 20° and 33°, respectively. Figures 16d16f show the initial model of YM, PR and density, respectively. The initial model is obtained by the interpolation and extrapolation of the well data under the guidance of the horizons, and then through Gaussian low-pass filtering.

 Download: larger image Figure 16. Partial angle stack seismic sections and the initial model. (a) Stack section of small incident angle, the mean angle is 7°; (b) stack section of medium incident angle, the mean angle is 20°; (c) stack section of large incident angle, the mean angle is 33°; (d) initial YM model; (e) initial PR model; (f) initial density model.

Figure 17 displays the inversion results by using our method. It can be observed that the inversion results have good lateral continuity, and the reflection horizons are clearly visible. There is a well near CDP 172. In order to observe the details of the inversion results, we compare the inversion results of the seismic traces near borehole with the well logs. Figure 18 shows these curves. Blue solid line and red dash-dot line represent the inversion result and well log, respectively. The well logs of YM and PR are calculated from the P wave velocity, S wave velocity and density data of the well. The calculation details can be referred to the work of Mavko et al. (2009). We can observe that the inversion results of YM and PR have a good fit with the well logs, however, the inversion result of density does not fit very well since the density term is very sensitive to noise and the limited angle range (Russell, 2014). Moreover, the obtained curves have a ladder-like shape, which is helpful for improving the resolution of the inversion results, and highlighting the boundaries of layers.

 Download: larger image Figure 17. Inversion results of YM, PR, and density obtained by using our method. (a) YM, (b) PR, (c) density.
 Download: larger image Figure 18. Well logs and the inversion results of the gather near borehole. (a) YM, (b) PR, (c) density.
4 CONCLUSIONS

In this article, we introduced the existing methods that can invert YM, PR and density from pre-stack seismic data, and analyzed the limits of these methods. In order to make improvement, we proposed multi-gather simultaneous inver- sion for prestack seismic data, introduced multiple regulari- zation constraints into prestack seismic inversion, and presented an algorithm which can obtain the elastic parameters directly and don't need to invert their reflectivity first. We used this method in the theoretical model and field data inversion experiments. The experimental results show that our method is not only applicable to the seismic data contaminated by strong Gaussian random noise, but also can reduce the influence of the noise that has a strong lateral continuity. Meanwhile, the results of the experiment confirm that the impact of the initial model is very little. Moreover, our method can improve the lateral continuity of the inversion results, and can effectively solve the problem of cumulative errors in existing methods.

Since the multi-gather simultaneous inversion is used in our method, the amount of data and the amount of computation in each iteration become larger. How to enhance the compu-tational efficiency of our method is the future work.

ACKNOWLEDGMENTS

This work was supported by the National Natural Science Foundation of China (Nos. 61775030, 61571096, 41301460, 61362018, and 41274127), and the key projects of Hunan Provincial Department of Education (No. 16A174). The authors thank the referees for their valuable suggestions. The authors also thank Chengdu Jingshi petroleum Technology Co., Ltd for providing us with the field data. The final publication is available at Springer via https://doi.org/10.1007/s12583-017-0905-7.

APPENDIX

Algorithm 1: Prestack multi-gather simultaneous inversion with multiple regularization constraints

Input: η, γ, $\mathit{\varphi }$, λ, α, tol

1. Initialize:$d_x^0 = d_y^0 = b_x^0 = b_y^0 = b_R^0 = {\mathit{\boldsymbol{R}}^0} = {\rm{0}}, \; \bar m{\rm{, }}$

2. While ${\left\| {{\mathit{\boldsymbol{m}}^k} - {\mathit{\boldsymbol{m}}^{k - 1}}} \right\|_2}/{\left\| {{\mathit{\boldsymbol{m}}^k}} \right\|_2} > {\rm{tol}}$

3. $\begin{array}{l} {\mathit{\boldsymbol{m}}^{k + 1}} = \mathop {\min }\limits_m \frac{\eta }{2}\left\| {\mathit{\boldsymbol{Gm}} - \mathit{\boldsymbol{d}}} \right\|_2^2 + \frac{\gamma }{2}\left\| {{\mathit{\boldsymbol{R}}^k} - \mathit{\boldsymbol{\tilde Dm}} - b_R^k} \right\|_2^2 + \\ \quad \quad \; \; \frac{\lambda }{2}\left\| {d_x^k - {\nabla _x}\mathit{\boldsymbol{m}} - b_x^k} \right\|_2^2 + \frac{\lambda }{2}\left\| {d_y^k - {\nabla _y}\mathit{\boldsymbol{m}} - b_y^k} \right\|_2^2\\ \; \; \; \; \; \; \; \; \; \; \; + \frac{\alpha }{2}\left\| {\mathit{\boldsymbol{m}} - \mathit{\boldsymbol{\bar m}}} \right\|_2^2 \end{array}$

4. ${h^{k + 1}} = \sqrt {{{\left| {{\nabla _x}{\mathit{\boldsymbol{m}}^{k + 1}} + b_x^k} \right|}^2} + {{\left| {{\nabla _y}{\mathit{\boldsymbol{m}}^{k + 1}} + b_y^k} \right|}^2}}$

5. $d_x^{k + 1} = \max ({h^{k + 1}} - \frac{\varphi }{{{\rm{2}}\lambda }}, 0)\frac{{{\nabla _x}{\mathit{\boldsymbol{m}}^{k + 1}} + b_x^k}}{{{h^{k + 1}}}}$

6. $d_y^{k + 1} = \max ({h^{k + 1}} - \frac{\varphi }{{2\lambda }}, 0)\frac{{{\nabla _y}{\mathit{\boldsymbol{m}}^{k + 1}} + b_y^k}}{{{h^{k + 1}}}}$

7. ${\mathit{\boldsymbol{R}}^{k + 1}} = {\rm{shrink}}(\mathit{\boldsymbol{\tilde D}}{\mathit{\boldsymbol{m}}^{k + 1}} + b_R^k, \frac{1}{\gamma })$

8. $b_x^{k + 1} = b_x^k + ({\nabla _x}{\mathit{\boldsymbol{m}}^{k + 1}} - d_x^{k + 1})$

9. $b_y^{k + 1} = b_y^k + ({\nabla _y}{\mathit{\boldsymbol{m}}^{k + 1}} - d_y^{k + 1})$

10. $b_R^{k + 1} = b_R^k + (\mathit{\boldsymbol{\tilde D}}{\mathit{\boldsymbol{m}}^{k + 1}} - {\mathit{\boldsymbol{R}}^{k + 1}})$

11. $k = k + 1$

12. End

13. $\mathit{\boldsymbol{I}} = {\rm{reshape}}(\exp {\rm{ }}(\mathit{\boldsymbol{m}}))$

14. $\begin{array}{l} \mathit{\boldsymbol{E}} = \mathit{\boldsymbol{I}}(1:N, \; :); \; \; \mathit{\boldsymbol{\sigma }} = \mathit{\boldsymbol{I}}(N + 1:2N, \; :); \; \\ \mathit{\boldsymbol{\rho }} = \mathit{\boldsymbol{I}}(2N + 1:3N, \; :) \end{array}$

REFERENCES CITED
 Baraniuk, R., 2007. Compressive Sensing. IEEE Signal Processing Magazine, 24(4): 118-121. DOI:10.1109/MSP.2007.4286571 Chen, S. S., Donoho, D. L., Saunders, M. A., 2001. Atomic Decomposition by Basis Pursuit. SIAM Review, 43(1): 129-159. DOI:10.1137/S003614450037906X Donoho, D. L., 2006. Compressed Sensing. IEEE Transactions on Information Theory, 52(4): 1289-1306. DOI:10.1109/TIT.2006.871582 Gholami, A., 2015. Nonlinear Multichannel Impedance Inversion by Total-Variation Regularization. Geophysics, 80(5). Goldstein, T., Osher, S., 2009. The Split Bregman Method for L1-Regularized Problems. SIAM Journal on Imaging Sciences, 2(2): 323-343. DOI:10.1137/080725891 Hamid, H., Pidlisecky, A., 2015. Multitrace Impedance Inversion with Lateral Constraints. Geophysics, 80(6). Harris, N. B., Miskimins, J. L., Mnich, C. A., 2011. Mechanical Anisotropy in the Woodford Shale, Permian Basin:Origin, Magnitude, and Scale. The Leading Edge, 30(3): 284-291. DOI:10.1190/1.3567259 Huang, W., Zhou, H. W., 2015. Least-Squares Seismic Inversion with Stochastic Conjugate Gradient Method. Journal of Earth Science, 26(4): 463-470. DOI:10.1007/s12583-015-0553-8 Mavko, G., Mukerji, T., Dvorkin, J., 2009. The Rock Physics Handbook:Tools for Seismic Analysis of Porous Media. Cambridge University Press, New York. 23 Pérez, D. O., Velis, D. R., Sacchi, M. D., 2014. Blocky Inversion of Prestack Seismic Data Using Mixed-Norms. 2014 SEG Annual Meeting, Denver. 3106-3111. https: //doi.org/10.1190/segam2014-0801.1 Rickman, R., Mullen, M. J., Petre, J. E., et al., 2008. A Practical Use of Shale Petrophysics for Stimulation Design Optimization: All Shale Plays are not Clones of the Barnett Shale. SPE Annual Technical Conference and Exhibition, Halliburton. 21-24. https://doi.org/10.2118/115258-MS Rudin, L. I., Osher, S., Fatemi, E., 1992. Nonlinear Total Variation Based Noise Removal Algorithms. Physica D:Nonlinear Phenomena, 60(1/2/3/4): 259-268. Russell, B. H., Dommico, S., 1988. Introduction to Seismic Inversion Methods. Society of Exploration Geophysicists, Tulsa. https://doi.org/10.1190/1.9781560802303 Russell, B. H., 2014. Prestack Seismic Amplitude Analysis:An Integrated Overview. Interpretation, 2(2): SC19-SC36. DOI:10.1190/INT-2013-0122.1 Sena, A., Castillo, G., Chesser, K., et al., 2011. Seismic Reservoir Characterization in Resource Shale Plays:Stress Analysis and Sweet Spot Discrimination. The Leading Edge, 30(7): 758-764. DOI:10.1190/1.3609090 Walker, C., Ulrych, T. J., 1983. Autoregressive Recovery of the Acoustic Impedance. Geophysics, 48(10): 1338-1350. DOI:10.1190/1.1441414 Yilmaz, Ö., 2001. Seismic Data Analysis. Society of Exploration Geophysicists, Tulsa. 1809. https://doi.org/10.1190/1.9781560801580 Yin, X. Y., Liu, X. J., Zong, Z. Y., 2015. Pre-Stack Basis Pursuit Seismic Inversion for Brittleness of Shale. Petroleum Science, 12(4): 618-627. DOI:10.1007/s12182-015-0056-3 Yuan, S. Y., Wang, S. X., Luo, C. M., et al., 2015. Simultaneous Multitrace Impedance Inversion with Transform-Domain Sparsity Promotion. Geophysics, 80(2): R71-R80. DOI:10.1190/geo2014-0065.1 Zhang, F. C., Dai, R. H., Liu, H. Q., 2014. Seismic Inversion Based on L1-Norm Misfit Function and Total Variation Regularization. Journal of Applied Geophysics, 109: 111-118. DOI:10.1016/j.jappgeo.2014.07.024 Zhang, J., Liu, H. S., Tong, S. Y., et al., 2015. Estimation of Elastic Parameters Using Two-Term Fatti Elastic Impedance Inversion. Journal of Earth Science, 26(4): 556-566. DOI:10.1007/s12583-015-0564-5 Zhang, R., Sen, M. K., Srinivasan, S., 2013. A Prestack Basis Pursuit Seismic Inversion. Geophysics, 78(1): R1-R11. Zong, Z.-Y., Yin, X.-Y., Zhang, F., et al., 2012. Reflection Coefficient Equation and Pre-Stack Seismic Inversion with Young's Modulus and Poisson Ratio. Chinese Journal of Geophysics, 55(11): 3786-3794. Zong, Z.-Y., Yin, X.-Y., Wu, G. C., 2013. Elastic Impedance Parameterization and Inversion with Young's Modulus and Poisson's Ratio. Geophysics, 78(6): N35-N42. DOI:10.1190/geo2012-0529.1