Least-squares wave-equation inversion of time-lapse seismic data sets - A Valhall case study

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}.$$ (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 following regression:

$$\begin{array}{ccc} \left [ \begin{array}{cc} {\bf L}_{0} & {\... ...0} \\ {\bf0} \\ \hline {\bf0} \end{array} \right ] \end{array},$$ (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 (Ayeni and Biondi, 2011)

$$\begin{array}{c} \left [ \begin{array}{ccc} {\bf H_0 } & {\bf... ...\\ \hline {\bf0} \\ {\bf0} \\ \end{array} \right ], \end{array}$$ (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. 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}$ ). Equation 3 can be extended to multiple seismic data sets (Ayeni and Biondi, 2010). Alternatively, we can re-write equation 3 to invert directly for the time-lapse image and a static baseline image (Ayeni and Biondi, 2011). Due to 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, which has been estimated to a high accuracy. However this method can be applied where an accurate monitor velocity has been available. The updated inversion problem is then given by (Ayeni and Biondi, 2011)

$$\begin{array}{c} \left [ \begin{array}{ccc} {\bf H }_{0} & {\... ...\\ \hline {\bf0} \\ {\bf0} \\ \end{array} \right ], \end{array}$$ (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 computed with the monitor geometry but with the baseline velocity. Note that 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},$$ (5)

while 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}.$$ (6)

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}),$$ (7)

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 )$ represent 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 number of off-diagonal elements of the Hessian that are computed -- which represents the size of the point spread function (PSF) at each image point (Lecomte and Gelius, 1998; Valenciano et al., 2006; Chavent and Plessix, 1999). $\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. Note that because of symmetry, only one half of the approximate Hessian is required. In this paper, we follow the phase-encoding approach of Tang (2009) to efficiently compute the target-oriented Hessian. The spatial regularization operators in equation 4 are non-stationary directional Laplacians (Hale, 2007), whereas the temporal constraint is the difference between the aligned images. Further review of the methodology is given by Ayeni and Biondi (2011,2010)

Least-squares wave-equation inversion of time-lapse seismic data sets - A Valhall case study