Image-space wave-equation tomography in the generalized source domain

Image-space encoded wavefields

The image-space encoded gathers are obtained using the prestack exploding-reflector modeling method introduced by Biondi (2006) and Biondi (2007). The general idea of this method is to model the data and the corresponding source function that are related to only one event in the subsurface, where a single unfocused SODCIG (obtained with an inaccurate velocity model) is used as the initial condition for the recursive upward continuation with the following one-way wave equations:

$$\left\{ \begin{array}{l} \left( \frac{\partial}{\partial z}-i\sqr... ... h};x_m,y_m) \\ Q_D(x,y,z=z_{\rm max},\omega;x_m,y_m) = 0 \end{array} \right.,$$ (A-28)

and

$$\left\{ \begin{array}{l} \left( \frac{\partial}{\partial z}-i\sqr... ... h};x_m,y_m) \\ Q_U(x,y,z=z_{\rm max},\omega;x_m,y_m) = 0 \end{array} \right.,$$ (A-29)

where $I_D({\bf x},{\bf h};x_m,y_m)$ and $I_U({\bf x},{\bf h};x_m,y_m)$ are the isolated SODCIGs at the horizontal location $(x_m,y_m)$ for a single reflector, and are suitable for the initial conditions for the source and receiver wavefields, respectively. They are obtained by rotating the original unfocused SODCIGs according to the apparent geological dip of the reflector. This rotation maintains the velocity information needed for migration velocity analysis, especially for dipping reflectors (Biondi, 2007). By collecting the wavefields at the surface, we obtain the areal source data $Q_D(x,y,z=0,\omega;x_m,y_m)$ and the areal receiver data $Q_U(x,y,z=0,\omega;x_m,y_m)$ for a single reflector and a single SODCIG located at $(x_m,y_m)$.

Since the size of the migrated image volume can be very big in practice and there are usually many reflectors in the subsurface, modeling each reflector and each SODCIG one by one may generate a data set even bigger than the original data set. One strategy to reduce the cost is to model several reflectors and several SODCIGs simultaneously (Biondi, 2006); however, this process generates unwanted crosstalk. As discussed by Guerra and Biondi (2008a,b), random phase encoding could be used to attenuate the crosstalk. The randomly encoded areal source and areal receiver wavefields can be computed as follows:

$$\left\{ \begin{array}{l} \left( \frac{\partial}{\partial z}-i\sqr... ..._m,\omega) \\ Q_D(x,y,z=z_{\rm max},{\bf p}_m,\omega) = 0 \end{array} \right.,$$ (A-30)

and

$$\left\{ \begin{array}{l} \left( \frac{\partial}{\partial z}-i\sqr... ..._m,\omega) \\ Q_U(x,y,z=z_{\rm max},{\bf p}_m,\omega) = 0 \end{array} \right.,$$ (A-31)

where $\widetilde{I}_D({\bf x},{\bf h},{\bf p}_m,\omega)$ and $\widetilde{I}_U({\bf x},{\bf h},{\bf p}_m,\omega)$ are the encoded SODCIGs after rotations. They are defined as follows:

$$\widetilde{I}_D({\bf x},{\bf h},{\bf p}_m,\omega) = \sum_{x_m}\sum_{y_m}I_D({\bf x},{\bf h},x_m,y_m)\beta({\bf x},x_m,y_m,{\bf p}_m,\omega),$$ (A-32)

$$\widetilde{I}_U({\bf x},{\bf h},{\bf p}_m,\omega) = \sum_{x_m}\sum_{y_m}I_U({\bf x},{\bf h},x_m,y_m)\beta({\bf x},x_m,y_m,{\bf p}_m,\omega),$$ (A-33)

where $\beta({\bf x},x_m,y_m,{\bf p}_m,\omega)=e^{i\gamma({\bf x},x_m,y_m,{\bf p}_m,\omega)}$ is chosen to be the random phase-encoding function, with $\gamma({\bf x},x_m,y_m,{\bf p}_m,\omega)$ being a uniformly distributed random sequence in ${\bf x}$, $x_m$, $y_m$ and $\omega$; the variable ${\bf p}_m$ is the index of different realizations of the random sequence. Recursively solving Equations 30 and 31 gives us the encoded areal source data $Q_D(x,y,z=0,{\bf p}_m,\omega)$ and areal receiver data $Q_U(x,y,z=0,{\bf p}_m,\omega)$, which can be collected on the surface.

The synthesized new data sets are downward continued using the same one-way wave equation defined by Equations 13 and 14 (with different boundary conditions) as follows:

$$\left\{ \begin{array}{l} \left( \frac{\partial}{\partial z}+i\sqr... ...bf p}_m,\omega) = \overline{Q_D(x,y,z=0,{\bf p}_m,\omega)} \end{array} \right.,$$ (A-34)

and

$$\left\{ \begin{array}{l} \left( \frac{\partial}{\partial z}+i\sqr... ...(x,y,z=0,{\bf x}_s,\omega) = Q_U(x,y,z=0,{\bf p}_m,\omega) \end{array} \right.,$$ (A-35)

where $\widetilde{D}({\bf x},{\bf p}_m,\omega)$ and $\widetilde{U}({\bf x},{\bf p}_m,\omega)$ are the downward continued areal source and areal receiver wavefields for realization ${\bf p}_m$. The image is produced by cross-correlating the two wavefields and summing images for all realization ${\bf p}_m$ as follows:

$$I_{\rm me}({\bf x},{\bf h}) = \sum_{{\bf p}_m}\sum_{\omega} \widetilde{D}({\bf x},{\bf p}_m,\omega) \widetilde{U}({\bf x},{\bf p}_m,\omega).$$ (A-36)

The crosstalk artifacts can be further attenuated if the number of ${\bf p}_m$ is large; therefore, approximately, the image obtained by migrating the image-space encoded gathers is kinematically equivalent to the image obtained in the shot-profile domain.

From Equation 36, the perturbed image is easily obtained as follows:

$$\Delta I_{\rm me}({\bf x},{\bf h}) = \sum_{{\bf p}_m}\sum_{\omega} \left( \Delta \widetilde{D}({\bf x}... ...}_m,\omega) \widehat{\widetilde{U}}({\bf x}+{\bf h},{\bf p}_m,\omega) + \right.$$

$$\left. \widehat{\widetilde{D}}({\bf x}-{\bf h},{\bf p}_m,\omega) \Delta \widetilde{U}({\bf x}+{\bf h},{\bf p}_m,\omega) \right),$$ (A-37)

where $\widehat{\widetilde{D}}({\bf x},{\bf p}_m,\omega)$ and $\widehat{\widetilde{U}}({\bf x},{\bf p}_m,\omega)$ are the image-space encoded background source and receiver wavefields; $\Delta \widetilde{D}({\bf x},{\bf p}_m,\omega)$ and $\Delta \widetilde{U}({\bf x},{\bf p}_m,\omega)$ are the perturbed source and receiver wavefields in the image-space phase-encoding domain, which satisfy the perturbed one-way wave equations defined by Equations 17 and 18. The tomographic operator ${\bf T}$ and its adjoint ${\bf T}^{\prime }$ can be implemented in a manner similar to that discussed in Appendices B and C, by replacing the original wavefields with the image-space phase-encoded wavefields.

Image-space wave-equation tomography in the generalized source domain