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
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 prestack seismic inversion we can obtain Poisson's ratio (PR) and Young's modulus (YM) from prestack 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 prestack 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 prestack 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 signaltonoise ratio. In order to improve inversion results, some scholars have introduced the L_{1} 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 prestack seismic inversion. BPI algorithm is based on the Basis pursuit (Chen et al., 2001) which is an L_{1} 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 prestack 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 twostep inversion): Firstly, imposing the L_{1} norm sparse constraint on the reflectivity of YM, PR and density, and solving the L_{1} norm optimization problem to obtain the sparse reflectivity. Secondly, obtaining YM, PR and density from their reflectivity.
In seismic inversion, L_{1} norm regularization technique is generally used to get the sparse solution in depth/time direction, without considering the lateral continuity. To reinforce lateral continuity, multitrace/multichannel simultaneous inversion is used in poststack seismic inversion. Hamid and Pidlisecky (2015) introduced lateral constraint into seismic impedance inversion and presented a multitrace impedance inversion method. In their experiments the lateral continuity is enhanced by using multitrace simultaneous inversion. Yuan et al. (2015) showed that multitrace 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 multitrace 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 prestack 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 prestack 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 twostep inversion methods will produce cumulative errors. The third disadvantage can use the formula
Motivated by the multitrace/multichannel simultaneous inversion used in poststack seismic inversion, we propose multigather simultaneous inversion for prestack seismic data. In seismic impedance inversion, blocky impedance structures can be obtained by using TV regularization (Gholami, 2015; Zhang et al., 2014), the antinoise ability of inversion methods can be improved by exploiting L_{1} norm regularization, and the nonuniqueness can be reduced through adding initial model constraint. They inspired us to impose these regularization constraints on prestack seismic inversion. In order to reduce the cumulative errors caused by twostep 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 EquationYM 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 prestack 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 StoP 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 prestack 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 L_{E} = 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 L_{E}, 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 i^{th} incident angle, W_{θi} is the wavelet matrix of the i^{th} incident angle. Let L=(L_{E}, 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) 
Our aim is to obtain accurate and stable inversion results of YM, PR and density, and to improve their lateral continuity. In impedance inversion, multitrace inversion can be used to explore spatial continuities of strata, which enlightened us to use it in the prestack seismic inversion. However, the prestack inversion is different from poststack inversion. In poststack inversion, each seismic trace corresponds to a trace of inversion result. While, in prestack inversion, multiple seismic traces from a gather corresponds to a trace of elastic parameter. In order to improve the lateral continuity of the prestack inversion results, we extend the multitrace inversion in poststack inversion to the multigather simultaneous inversion in the prestack inversion.
By using the multigather simultaneous inversion, the initial objective function is established by following steps. In our method, S_{i} (i = 1, 2 … N) represents the i^{th} prestack gather, which can be seen from Fig. 1.
Download:


We arrange each gather to a column vector.
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.
L_{1} norm constraint is used to enhance the antinoise performance. Based on the initial objective function, we impose the L_{1} 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
$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
${\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
$\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 nonuniqueness 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
As can be seen from Eq. (18), the objective function has a complicated form with L_{1} norm, L_{2} 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 b_{x}, b_{y} and b_{R} 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, (d_{x}, d_{y}) 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 d_{x} and d_{y}, 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 L_{2} 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 d_{x} and d_{y}. 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 EXPERIMENTSWe 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:


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:


Figure 4 shows the initial elastic parameters model generated from the true model by Gauss lowpass 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:


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 dashdot 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:


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:


Download:


Download:


Our method is a kind of modelbased 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 10^{7}(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:


Download:


Download:


Download:


By comparing Figs. 5–8 and Figs. 9–12, 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 backdiagonal of the synthetic seismic section. Since it's not convenient to display the noisy synthetic prestack 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:


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 multigather 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:


Download:


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. 16a–16c. 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 16d–16f 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 lowpass filtering.
Download:


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 dashdot 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 ladderlike shape, which is helpful for improving the resolution of the inversion results, and highlighting the boundaries of layers.
Download:


Download:


In this article, we introduced the existing methods that can invert YM, PR and density from prestack seismic data, and analyzed the limits of these methods. In order to make improvement, we proposed multigather 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 multigather 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 computational efficiency of our method is the future work.
ACKNOWLEDGMENTSThis 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/s1258301709057.
APPENDIXAlgorithm 1: Prestack multigather simultaneous inversion with multiple regularization constraints
Input: η, γ,
1. Initialize:
2. While
3.
4.
5.
6.
7.
8.
9.
10.
11.
12. End
13.
14.
REFERENCES CITED
Baraniuk, R., 2007. Compressive Sensing. IEEE Signal Processing Magazine, 24(4): 118121. 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): 129159. DOI:10.1137/S003614450037906X 
Donoho, D. L., 2006. Compressed Sensing. IEEE Transactions on Information Theory, 52(4): 12891306. DOI:10.1109/TIT.2006.871582 
Gholami, A., 2015. Nonlinear Multichannel Impedance Inversion by TotalVariation Regularization. Geophysics, 80(5). 
Goldstein, T., Osher, S., 2009. The Split Bregman Method for L1Regularized Problems. SIAM Journal on Imaging Sciences, 2(2): 323343. 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): 284291. DOI:10.1190/1.3567259 
Huang, W., Zhou, H. W., 2015. LeastSquares Seismic Inversion with Stochastic Conjugate Gradient Method. Journal of Earth Science, 26(4): 463470. DOI:10.1007/s1258301505538 
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 MixedNorms. 2014 SEG Annual Meeting, Denver. 31063111. https: //doi.org/10.1190/segam20140801.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. 2124. https://doi.org/10.2118/115258MS

Rudin, L. I., Osher, S., Fatemi, E., 1992. Nonlinear Total Variation Based Noise Removal Algorithms. Physica D:Nonlinear Phenomena, 60(1/2/3/4): 259268. 
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): SC19SC36. DOI:10.1190/INT20130122.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): 758764. DOI:10.1190/1.3609090 
Walker, C., Ulrych, T. J., 1983. Autoregressive Recovery of the Acoustic Impedance. Geophysics, 48(10): 13381350. 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. PreStack Basis Pursuit Seismic Inversion for Brittleness of Shale. Petroleum Science, 12(4): 618627. DOI:10.1007/s1218201500563 
Yuan, S. Y., Wang, S. X., Luo, C. M., et al., 2015. Simultaneous Multitrace Impedance Inversion with TransformDomain Sparsity Promotion. Geophysics, 80(2): R71R80. DOI:10.1190/geo20140065.1 
Zhang, F. C., Dai, R. H., Liu, H. Q., 2014. Seismic Inversion Based on L1Norm Misfit Function and Total Variation Regularization. Journal of Applied Geophysics, 109: 111118. DOI:10.1016/j.jappgeo.2014.07.024 
Zhang, J., Liu, H. S., Tong, S. Y., et al., 2015. Estimation of Elastic Parameters Using TwoTerm Fatti Elastic Impedance Inversion. Journal of Earth Science, 26(4): 556566. DOI:10.1007/s1258301505645 
Zhang, R., Sen, M. K., Srinivasan, S., 2013. A Prestack Basis Pursuit Seismic Inversion. Geophysics, 78(1): R1R11. 
Zong, Z.Y., Yin, X.Y., Zhang, F., et al., 2012. Reflection Coefficient Equation and PreStack Seismic Inversion with Young's Modulus and Poisson Ratio. Chinese Journal of Geophysics, 55(11): 37863794. 
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): N35N42. DOI:10.1190/geo20120529.1 