The Slope-Attribute-Regularized High-Resolution Prestack Seismic Inversion

Abstract

Prestack seismic inversion can be regarded as an optimization problem, which minimizes the error between the observed and synthetic data under the premise of certain geological/geophysical a priori information constraints. It has been proved to be a powerful approach for reconstructing the subsurface properties and building the elastic parameter models (e.g., P- and S-wave velocity, and density). With respect to the specific expressions of a priori information, the starting model and regularization are expected to be the most widely used and indispensable constraints to reconstruct structural features and subsurface properties. The conventional prestack inversion (trace-by-trace) methods perform well when the geological structure of the target area is not too complex. However, due to the lack of lateral constraint, such trace-independent methods are inevitably limited by their capability of characterization (including accuracy, resolution, and robustness) in the case of geologically complex structures, such as tilted stratum and steep faults. The geological structure-guided constraint, herein referred to as the seismic slope attribute, can be exploited as a lateral constraint integrated into the prestack inversion algorithm. In this work, the seismic slope attribute is introduced to the amplitude variation with offset/angle inversion from two aspects, i.e., starting model building and regularization penalty. Firstly, using the seismic slope attribute, instead of the traditional manual interpreted geological horizons, as a constraint, the well-log data are interpolated to build the initial model. The interpolation algorithm is formulated as solving the inverse problem by using the shaping regularization method rather than the kriging-based algorithm. Secondly, by rotating the coordinate system according to the seismic slope attribute, the directional total variation regularization is used as a constraint to improve the resolution (in both vertical and horizontal directions) and lateral continuity of the inversion results. Finally, the proposed methods are applied to synthetic and real seismic data. Synthetic tests and field data applications demonstrate that the proposed method is capable of revealing complex structural features and achieving stabilized inversion of multi-parameters with less uncertainty.

Appendices

Appendix A: The Exact Zoeppritz Equation

Without loss of generality, the forward problem of seismic wave propagation can be expressed as a nonlinear equation as follows:

$$\begin{aligned} \mathbf{d }=\mathbf{G }(\mathbf{m }). \end{aligned}$$

(A.1)

According to the convolution theory, the seismic data can be considered as the convolution between the stationary wavelet and reflectivity: Then, the P–P seismic data can be simulated by convolving the P–P reflectivity coefficient with the stationary wavelet as:

$$\begin{aligned} \mathbf{d }=\mathbf{W }\mathbf{R }(\mathbf{m }). \end{aligned}$$

(A.2)

With respect to the prestack seismic inversion, the reflectivity coefficients can be obtained by the Zoeppritz’s equation. When a plane-wave propagates onto a surface, according to the Zoeppritz’s equation, the reflection and transmission coefficients can be expressed as:

$$\begin{aligned} \mathbf{AR} =\mathbf{B }, \end{aligned}$$

(A.3)

where

$$\begin{aligned} \mathbf{A }= & {} \begin{bmatrix} \sin \theta _{1} &{} \cos \phi _{1} &{} -\sin \theta _{2} &{} \cos \phi _{2}\\ \cos \theta _{1} &{} -\sin \phi _{1} &{} \cos \theta _{2} &{} \sin \phi _{2}\\ \sin 2\phi _{1} &{} -\frac{V_{s1}}{V_{p1}}\sin 2\phi _{1} &{} -\frac{\rho _{2}V_{p2}}{\rho _{1}V_{p1}}\cos 2\phi _{2} &{} -\frac{\rho _{2}V_{p2}}{\rho _{1}V_{p1}}\sin 2\phi _{2}\\ \sin 2\theta {1} &{} \frac{V_{p1}}{V_{s1}}\sin \phi _{1} &{} \frac{\rho _{2}V_{s2}^{2}V_{p2}}{\rho _{1}V_{s2}^{2}V_{p1}}\cos 2\phi _{2} &{} -\frac{\rho _{2}V_{s2}V_{p1}}{\rho _{1}V_{s1}^{2}}\sin 2\phi _{2}\\ \end{bmatrix},\quad \mathbf{R }= \begin{bmatrix} R_{pp}\\ R_{ps} \\ T_{pp} \\ T_{ps} \end{bmatrix},\nonumber \\ \mathbf{B }= & {} \begin{bmatrix} -\sin 2\theta \\ \cos 2\theta \\ -\cos 2\phi \\ \sin 2\phi \end{bmatrix} \end{aligned}$$

(A.4)

where ${\textit{v}_{p1}}$, ${\textit{v}_{s1}}$, and ${\rho _{1}}$ denote the elastic parameters of the upper layers, and ${\textit{v}_{p1}}$, ${\textit{v}_{s1}}$, and ${\rho _{2}}$ correspond to the counterpart of the lower layers, $\theta _{1}$ and $\phi _{1}$ are the angles of the P- and S-wave reflections, and $\theta _{1}$ and $\phi _{1}$ stand for the angles of P- and S-wave transmissions.

The partial derivative with respect to the parameter m can be expressed as:

$$\begin{aligned}&\frac{\partial \mathbf{R }}{\partial \mathbf{m }}=\mathbf{A }^{-1}\left( \frac{\partial \mathbf{B }}{\partial \mathbf{m }}-\frac{\partial \mathbf{A }}{\partial \mathbf{m }}\mathbf{R }\right) , \end{aligned}$$

(A.5)

$$\begin{aligned}&\begin{aligned} \frac{\partial \mathbf{A }}{\partial \textit{V}_{p1}}&=\frac{1}{V_{p1}}\cdot \\&\begin{bmatrix} 0&{}\tan \phi \sin \phi _{1} &{}\sin \theta _{2} &{}\tan \phi _{2}\sin \phi _{2}\\ 0&{}\sin \phi _{1} &{}\sin \theta _{2}\tan \theta _{2}&{}-\sin \phi _{2}\\ 2(1-\cos 2\phi _{1})&{}(2-\tan ^{2}\phi _{1})\frac{V_{s1}}{V_{p1}}\sin 2\phi _{1} &{}\frac{\rho _{2}}{\rho _{1}}\frac{V_{p2}}{V_{p1}}(3\cos 2\phi _{2}-2) &{}(2-\tan ^{2}\phi _{2})\frac{\rho _{2}}{\rho _{1}}\frac{V_{s2}}{V_{p1}\sin 2\phi _{2}}\\ 0&{}\frac{V_{p1}}{V_{s1}}(2-\cos 2\phi _{1})&{} \frac{\rho _{2}}{\rho _{1}}\frac{V^{2}_{s2}}{V^{2}_{s1}}\frac{V_{p1}}{V_{p2}} \sin 2\theta _{2}\tan ^{2}\theta _{2}&{}\frac{\rho _{2}V_{p1}V_{s2}}{\rho _{1} V^{2}_{s1}}(\cos 2\phi _{2}-2) \end{bmatrix}, \end{aligned}\nonumber \\ \end{aligned}$$

(A.6)

$$\begin{aligned}&\frac{\partial \mathbf{A }}{\partial \textit{V}_{s1}}=\frac{1}{V_{s1}} \begin{bmatrix} 0&{}-\tan \phi _{1}\sin \phi _{1}&{}0&{}0\\ 0&{}-\sin \phi _{1}&{}0&{}0\\ 4\sin ^{2}\phi _{1}&{}(\tan ^{2}\phi _{1}-2)\frac{V_{s1}}{V_{p1}}\sin 2\phi _{1}&{}0&{}0\\ 0&{}\frac{V_{p1}}{V_{s1}}(\cos 2\phi _{1}-2)&{}-2\frac{\rho _{2}V^{2}_{s2} V_{p1}}{\rho _{1}V^{2}_{s1}V_{p2}}\sin 2\theta _{2}&{}2\frac{\rho _{2}}{\rho _{1}} \frac{V_{s2}V_{p1}}{V^{2}_{s1}}\cos 2\phi _{2} \end{bmatrix} \end{aligned}$$

(A.7)

$$\begin{aligned}&\frac{\partial \mathbf{A }}{\partial {\rho }_{1}}=\frac{1}{\rho _{1}} \begin{bmatrix} 0&{}0&{}0&{}0\\ 0&{}0&{}0&{}0\\ 0&{}0&{}\frac{\rho _{2}}{\rho _{1}}\frac{V_{p2}}{V_{p1}}\cos 2\phi _{2}&{} \frac{\rho _{2}}{\rho _{1}}\frac{V_{s2}}{V_{p1}}\sin 2\phi _{2}\\ 0&{}0&{}-\frac{\rho _{2}}{\rho _{1}}\frac{V^{2}_{s2}}{V^{2}_{s1}} \frac{V_{p2}}{V_{p1}}\sin 2\theta _{2}&{}\frac{\rho _{2}}{\rho _{1}} \frac{V_{s2}V_{p1}}{V^{2}_{s1}}\cos 2\phi _{2}\\ \end{bmatrix}, \end{aligned}$$

(A.8)

$$\begin{aligned}&\frac{\partial \mathbf{A }}{\partial \textit{V}_{p2}}=\frac{1}{V_{p2}} \begin{bmatrix} 0&{}0&{}-\sin \theta _{2}&{}0\\ 0&{}0&{}-\sin \theta _{2}\tan \theta _{2}&{}0\\ 0&{}0&{}-\frac{\rho _{2}}{\rho _{1}}\frac{V_{p2}}{V_{p1}}\cos 2\phi _{2}&{}0\\ 0&{}0&{}-\frac{\rho _{2}}{\rho _{1}}\frac{V^{2}_{s2}}{V^{2}_{s1}} \frac{V_{p1}}{V_{p2}}\sin 2\theta _{2}\tan ^{2}\theta _{2}&{}0\\ \end{bmatrix}, \end{aligned}$$

(A.9)

$$\begin{aligned}&\frac{\partial \mathbf{A }}{\partial \textit{V}_{s2}}=\frac{1}{V_{s2}} \begin{bmatrix} 0&{}0&{}0&{}-\tan \phi _{2}\sin \phi _{2}\\ 0&{}0&{}0&{}\sin \phi _{2}\\ 0&{}0&{}2\frac{\rho _{2}}{\rho _{1}}\frac{V_{p2}}{V_{p1}}(1-cos2\phi _{2}) &{}\frac{\rho _{2}}{\rho _{1}}\frac{V_{s2}}{V_{p1}}\sin 2\phi _{2}(2\tan ^{2}\phi _{2}-2)\\ 0&{}0&{}\frac{\rho _{2}}{\rho _{1}}\frac{V^{2}_{s2}}{V^{2}_{s1}}\frac{V_{p1}}{V_{p2}} \sin 2\theta _{2}&{}-\frac{\rho _{2}}{\rho _{1}}\frac{V_{p1}V_{s2}}{V^{2}_{s1}}(2+3\cos \phi _{2}) \end{bmatrix}, \end{aligned}$$

(A.10)

$$\begin{aligned}&\frac{\partial \mathbf{A }}{\partial {\rho }_{1}}=\frac{1}{\rho _{2}} \begin{bmatrix} 0&{}0&{}0&{}0\\ 0&{}0&{}0&{}0\\ 0&{}0&{}-\frac{\rho _{2}}{\rho _{1}}\frac{V_{p2}}{V_{p1}}\cos 2\phi _{2} &{}-\frac{\rho _{2}}{\rho _{1}}\frac{V_{s2}}{V_{p1}}\sin 2\phi _{2}\\ 0&{}0&{}\frac{\rho _{2}}{\rho _{1}}\frac{V^{2}_{s2}}{V^{2}_{s1}}\frac{V_{p2}}{V_{p1}} \sin 2\theta _{2}&{}-\frac{\rho _{2}}{\rho _{1}}\frac{V_{s2}V_{p1}}{V^{2}_{s1}}\cos 2\phi _{2}\\ \end{bmatrix}, \end{aligned}$$

(A.11)

$$\begin{aligned}&\frac{\partial \mathbf{B }}{\partial \textit{V}_{p1}}=\frac{1}{V_{p1}} \begin{bmatrix} 0&0&-4sin^{2}\phi _{1}&0 \end{bmatrix}^{T},\qquad \frac{\partial \mathbf{B }}{\partial \textit{V}_{p2}}=\frac{\partial \mathbf{B }}{\partial \textit{V}_{s1}}=\frac{\partial \mathbf{B }}{\partial \textit{V}_{s2}}=\frac{\partial \mathbf{B }}{\partial \rho _{1}}=\frac{\partial \mathbf{B }}{\partial {\rho }_{2}}=\begin{bmatrix} 0&0&0&0 \end{bmatrix}^{T}, \end{aligned}$$

(A.12)

Appendix B: The Hertz–Mindlin Model

The Hertz–Mindlin contact model (Mindlin 1949) calculates the bulk and shear modulus of two spherical grains in contact. It appears to be the most commonly used contact model to describe seismic parameter changes caused by the pressure changes (Dadashpour et al. 2007). Although the Hertz–Mindlin contact model is proved to be only applicable to perfect elastic contacts of spherical bodies, it works fairly well for sandstones (Avseth et al. 2005). According to he Hertz–Mindlin theory, the effective bulk modulus and shear modulus of a dry random identical sphere packing can be expressed as

$$\begin{aligned}&\textit{K}_{H\text {-}M}=\left[ \frac{c^2(1-\phi _{c})^2G^2P_{\text {eff}}}{18\pi ^2(1-\nu )^2}\right] ^{\frac{1}{n}}, \end{aligned}$$

(B.1)

$$\begin{aligned}&\textit{G}_{H\hbox {-}M}=\frac{5-4\nu }{5(2-\nu )} \left[ \frac{3c^2(1-\phi _{c})^2G^2P_{\text {eff}}}{2\pi ^2(1-\nu )^2}\right] ^{\frac{1}{n}}, \end{aligned}$$

(B.2)

where $\textit{K}_{H\text {-}M}$ and $\textit{G}_{H\text {-}M}$ indicate the bulk and shear modulus calculated by the Hertz–Mindlin model, $\phi _{c}$ denotes critical porosity, $\textit{P}_{\text {eff}}$ represents the effective pressure, and $\textit{G}$ and $\nu$ are the shear modulus and Poisson’s ratio of the solid grains, respectively. n is the coordination number and c denotes the average number of contacts per sphere. In the original Hertz–Mindlin theory, n is equal to 3, which indicates that the variation in velocity is proportional to the 1/6 power of $\textit{P}_{\text {eff}}$.

Some laboratory measurements of samples gave a larger number for n. Vidal et al. (2000) found $\textit{n}=5.6$ for P-wave and $\textit{n}=3.8$ for S-wave for gas-saturated sands, while Landro et al. (2001) used $\textit{n}=5$ for oil-saturated sands.