Wave-equation inversion of time-lapse seismic data sets

Methodology

Given a linearized modeling operator ${\bf L}$ , the seismic data ${\bf d}$ for survey ${i}$ due to a reflectivity model ${\bf m}$ is

$${\bf d}_{i}={\bf L}_{i}{\bf m}_{i}.$$ (A-1)

Assuming we have two data sets (baseline ${\bf d}_{0}$ and monitor ${\bf d}_{1}$ ) acquired at different times over an evolving reservoir, joint least-squares migration/inversion involves solving the regression

$$\begin{array}{ccc} \left [ \begin{array}{cc} {\bf L}_{0} & {\... ...0} \\ {\bf0} \\ \hline {\bf0} \end{array} \right ] \end{array},$$ (A-2)

where ${\bf R}_{i}$ and ${{\Lambda}}_i$ are the spatial and temporal regularization operators respectively, and $\epsilon_{i}$ and ${\bf {\zeta }}_i$ are the corresponding regularization parameters. Although we can directly solve equation 2 by minimizing the quadratic-norm of the regression (Ajo-Franklin et al., 2005), we choose to transform it to an image space problem of the form (Appendix A)

$$\begin{array}{c} \left [ \begin{array}{ccc} {\bf H_0 } & {\bf... ...\\ \hline {\bf0} \\ {\bf0} \\ \end{array} \right ], \end{array}$$ (A-3)

where ${\bf H}_{ i}={\bf L}^{^T}_{ i}{\bf L}_{ i}$ is the wave-equation Hessian, and ${\bf R}_{ {ii}}=\epsilon^{2}_{ i}{\bf R}_{ i}^{^T}{\bf R}_{ i}$ and ${{\Lambda}_{ij} = {{\zeta_i {\Lambda}_{ i}}}^{\it ^T} {\bf\zeta}_j {\Lambda}_j}$ are the spatial and temporal constraints. Note that ${\bf L}^{^T}$ , the migration operator, is the adjoint of the modeling operator ${\bf L}$ . The inverted time-lapse image $\Delta \hat {{\bf m}}$ is then the difference between the inverted baseline and monitor images ( ${\hat{{\bf m}}}_{0}$ and ${\hat{{\bf m}}}_{1}$ ). Alternatively, we can re-write equation 3 to invert directly for the time-lapse image and a static baseline image (Appendix A). Furthermore, equation 3 can be extended to multiple seismic data sets (Ayeni and Biondi, 2010).

As shown in Appendix A, if there are physical movements of reflectors and velocity changes (e.g., due to reservoir depletion and compaction) between surveys, the baseline and monitor images will not be aligned. Such misalignments must be accounted for before or during inversion. As is the case in many practical time-lapse monitoring problems, we assume that the monitor data are migrated with the baseline velocity. The updated inversion problem is then given by (Appendix A)

$$\begin{array}{c} \left [ \begin{array}{ccc} {\bf H }_{0} & {\... ...\\ \hline {\bf0} \\ {\bf0} \\ \end{array} \right ], \end{array}$$ (A-4)

where ${\tilde {\bf m}^{b}_{1}}$ and ${\bf\hat m}^{b}_{1}$ are respectively the migrated and inverted monitor images repositioned (warped) to the baseline image. The superscript $^{b}$ on the operators denotes that they are referenced to the baseline image. For example, ${\bf H}^{b}_{1}$ is the Hessian matrix with the monitor geometry but with the baseline velocity. Note that whereas the conventional time-lapse image $\Delta \tilde{\bf m}^{b}$ estimated at the baseline position is given by

$$\Delta \tilde{\bf m}^{b}= \tilde{\bf m}^{b}_{1} - \tilde{{\bf m}}_{0},$$ (A-5)

the inverted time-lapse image $\Delta \hat{\bf m}^{b}$ is given by

$$\Delta \hat{\bf m}^{b}= \hat{\bf m}^{b}_{1} - \hat{{\bf m}}_{0}.$$ (A-6)

The wave-equation Hessian at image point ${\bf x}$ is defined as (Valenciano et al., 2006; Plessix and Mulder, 2004)

$${H} \left ({\bf x}, {\bf x'} \right ) =\sum_{\omega}\omega^{4} \s... ..._{s}, {\bf x},{\omega}) \bar G ({\bf x}_{s}, {\bf x'},{\omega}) \hspace{-0.8cm}$$

$$\sum_{{\bf x}_{r}} G ({\bf x}, {\bf x}_{r},{\omega}) \bar G ({\bf x'}, {\bf x}_{r}, {\omega}) ,$$ (A-7)

where ${\bf x'}$ denotes all image points, $\bar G$ is the complex conjugate of Green's function $G$ at frequency ${\omega}$ , ${\it f_{s}}$ is the source function, and ${\bf x}_{s}$ and ${\bf x}_{r}$ are the source and receiver positions, respectively. For any practical application, it is infeasible (and unnecessary) to compute the full Hessian matrix. Because the problem is posed in the image space, we only need to compute the Hessian for a target region of interest around the reservoir. In addition, we only compute off-diagonal elements sufficient to capture the dominant structure of the Hessian. This target-oriented approximation of the Hessian is given by (Valenciano et al., 2006)

$${H} \left ({\bf x_{T}}, {\bf x_{T+a_{x}}} \right ) = \sum_{\omega... ...bf x_{s}, x_{T},\omega}) \bar G ({\bf x_{s}, x_{T+a_x},\omega}) \hspace{-0.8cm}$$

$$\sum_{\bf x_{r}} G ({\bf x_{T}, x_{r}, \omega}) \bar G ({\bf x_{T+a_x}, x_{r}, \omega}),$$ (A-8)

where ${\bf x_{T}}$ is an image point within the target area, and ${\bf x_{T+a_x}}$ represents points within a small region around ${\bf x_{T}}$ . For any image point, elements of ${\bf H} \left ({\bf x_{T}}, {\bf x_{T+a_{x}}} \right )$ represents a row of a sparse Hessian matrix ${\bf H}$ whose non-zero components are defined by ${\bf a_{x}}$ . Therefore, ${\bf a_{x}}$ defines the how many off-diagonal elements of the Hessian are computed -- which represents the size of the point spread function (PSF) around each image point. Note that because of symmetry, only one half of the approximate Hessian is required. The target-oriented Hessian and computational savings are discussed in further detail by Valenciano et al. (2006) and Tang (2008).

In this paper, the spatial regularization operators are non-stationary dip filters. First, we estimate the local dips on the migrated baseline image using a plane-wave destruction method (Fomel, 2002). Then we construct the operator based on factorization of directional Laplacian representations of the local dip filters (Hale, 2007). The temporal constraint is implemented as a difference between the aligned images. To attenuate multiples and other unwanted artifacts in the data, we perform a high-resolution Radon demultiple on the data. In order to align the baseline and monitor images prior to inversion, we perform pre-stack warping using a cyclic local cross-correlation method (Ayeni, 2010).

In the next section, we apply the proposed method to a field time-lapse data set. We show applications to complete and incomplete monitor data sets.

Wave-equation inversion of time-lapse seismic data sets