Target-oriented wavefield tomography using demigrated Born data

Target-oriented Born wavefield modeling

Seismic images can be obtained by applying shot-profile migration to the recorded data as follows:

$$m({\bf x},{\bf h}) = \frac{1}{H({\bf x},{\bf h})}\sum_{\omega}\sum_{{\bf x}_s}G_0^{*}({\bf x}-{\bf h},{\bf x}_s,\omega) \sum_{{\bf x}_r}W({\bf x}_r,{\bf x}_s)$$

$$\times G_0^{*}({\bf x}+{\bf h},{\bf x}_r,\omega)d({\bf x}_r,{\bf x}_s,\omega)$$ (1)

where ${}^{*}$ denotes adjoint, $\omega$ is the angular frequency, and $m({\bf x},{\bf h})$ is the migrated image as a function of both image point ${\bf x}=(x,y,z)$ and subsurface half offset ${\bf h}=(h_x,h_y,h_z)$ . The frequency-domain data $d({\bf x}_r,{\bf x}_s,\omega)$ is recorded at receiver position ${\bf x}_r=(x_r,y_r,0)$ due to a source located at ${\bf x}_s=(x_s,y_s,0)$ ; $W({\bf x}_r,{\bf x}_s)$ is the acquisition mask operator, which contains unity values where we record data and zero values where we do not; $G_0({\bf x},{\bf x}_s,\omega)$ and $G_0({\bf x},{\bf x}_r,\omega)$ are the Green's functions connecting the source and receiver, respectively, to image point ${\bf x}$ . The Green's functions are obtained using a starting velocity model ${\bf v}_0$ . Operator $H({\bf x},{\bf h})$ is the diagonal of the subsurface-offset-domain imaging Hessian defined as follows (Tang, 2009; Plessix and Mulder, 2004; Valenciano, 2008):

$$H({\bf x},{\bf h}) = \sum_{\omega}\sum_{{\bf x}_s}\sum_{{\bf x}_r... ...bf h},{\bf x}_s,\omega)\vert^2\vert G({\bf x}+{\bf h},{\bf x}_r,\omega)\vert^2,$$ (2)

the inverse of which partially compensates for the image distortion caused by uneven subsurface illumination and removes the effects of the original acquisition geometry (Tang, 2009; Rickett, 2003; Plessix and Mulder, 2004; Symes, 2008). By doing so, the initial image $m({\bf x},{\bf h})$ contains only the effects of the initial velocity used for migration and is independent from the way the original data is recorded.

We then use Born modeling (Stolt and Benson, 1986) to simulate a new data set based on the initial image obtained using equation 1. The modeling can be performed using arbitrary source functions and arbitrary acquisition geometries. The modeled data set using encoded areal sources reads (Tang, 2008)

$${\widetilde d}({\bf p}_s,{\widetilde {\bf x}}_r,\omega) = \sum_{{... ...\bf x}}+{\bf h},{\widetilde {\bf x}}_r,\omega) m({\widetilde {\bf x}},{\bf h}),$$ (3)

where ${\widetilde {\bf x}}=({\widetilde x},{\widetilde y},{\widetilde z})$ is an image point within the selected target zone; ${\widetilde {\bf x}}_r=({\widetilde x}_r,{\widetilde y}_r,{\widetilde z}_r)$ is the receiver location used for modeling. It can differ substantially from the original receiver location ${\bf x}_r$ , which can be only on the surface, i.e., ${\bf x}_r=(x_r,y_r,0)$ ; $G_0({\widetilde {\bf x}},{\bf p}_s,\omega)$ is the Green's function obtained using the encoded areal source, it can be written as a weighted sum of point-source Green's functions as follows:

$$G_0({\widetilde {\bf x}},{\bf p}_s,\omega) = \sum_{\widetilde {\b... ...,{\widetilde {\bf x}}_s,\omega)\alpha({\widetilde {\bf x}_s},{\bf p}_s,\omega),$$ (4)

where $\alpha({\widetilde {\bf x}}_s,{\bf p}_s,\omega)$ and ${\bf p}_s$ are the encoding function and encoding index, respectively, to be specified later. Similar to the new receiver location ${\widetilde {\bf x}}_r$ , the new source location ${\widetilde {\bf x}}_s=({\widetilde x}_s,{\widetilde y}_s,{\widetilde z}_s)$ can also be substantially different from the original source location ${\bf x}_s$ . Note in particular that the Green's function used for modeling ( $G_0({\widetilde {\bf x}},{\bf p}_s,\omega)$ ) is computed using the same velocity model (${\bf v}_0$ ) as that used for migrating the original data. Therefore, the modeling process undoes the effect of the starting velocity model ${\bf v}_0$ (at least in kinematics), resulting in a starting-velocity-independent data set.

The synthesized areal data can be imaged by using conventional migration of areal-source data as follows:

$${\widehat m}({\widetilde {\bf x}},{\bf h}) = \sum_{\omega} \sum_{{\bf p}_s} G^{*}({\widetilde {\bf x}}-{\bf h}... ...lde {\bf x}}_r,\omega) {\widetilde d}({\bf p}_s,{\widetilde {\bf x}}_r,\omega),$$ (5)

where Green's functions $G$ 's are obtained using velocity model ${\bf v}$ , which can be the same or different from the starting velocity model ${\bf v}_0$ , to be discussed later. Imaging the areal-source data using equation 5, however, generates crosstalk artifacts (Tang, 2008; Romero et al., 2000; Liu et al., 2006). To attenuate the crosstalk, the encoding function $\alpha$ is chosen such that

$$\sum_{{\bf p}_s} \alpha({\widetilde {\bf x}}_s,{\bf p}_s,\omega) ... ...bf p}_s,\omega) \approx \delta({\widetilde {\bf x}}_s-{\widetilde {\bf x}}_s'),$$ (6)

where $\delta(\cdot)$ is the Dirac delta function. For plane-wave-phase encoding, $\alpha({\widetilde {\bf x}}_s,{\bf p}_s,\omega)=A(\omega)e^{i\omega{\bf p}_s \cdot {\widetilde {\bf x}_s}}$ , with $A^2(\omega) =\vert\omega\vert$ in two dimensions and $A^2(\omega) = \vert\omega\vert^2$ in three dimensions, and ${\bf p}_s= (p_{s_x},p_{s_y},0)$ is the ray parameter for the source plane waves at depth level ${\widetilde z}_s$ .

Substituting equations 4 and 3 into equation 5 yields

$${\widehat m}({\widetilde {\bf x}},{\bf h}) = \sum_{{\widetilde {\... ...f x}},{\bf h},{\widetilde {\bf x}}',{\bf h}')m({\widetilde {\bf x}}',{\bf h}'),$$ (7)

where

$$\Delta G({\widetilde {\bf x}},{\bf h},{\widetilde {\bf x}}',{\bf h}') = \sum_{\omega} \sum_{\widetilde {\bf x}_s} \sum_{\widetilde {\bf x... ...x}}_s,\omega) G_0({\widetilde {\bf x}}'-{\bf h}',{\widetilde {\bf x}}_s,\omega)$$

$$\times G^{*}({\widetilde {\bf x}}+{\bf h},{\widetilde {\bf x}}_r,\omega) G_0({\widetilde {\bf x}}'+{\bf h}',{\widetilde {\bf x}}_r,\omega).$$ (8)

When the same velocity model is used for Born modeling and migration ( ${\bf v}_0={\bf v}$ ), $\Delta G({\widetilde {\bf x}},{\bf h},{\widetilde {\bf x}}',{\bf h}')$ becomes the local Hessian operator or resolution function (Tang, 2009; Lecomte, 2008; Valenciano, 2008) under the new acquisition geometry. It has zero phase and is centered at ${\widetilde {\bf x}}'={\widetilde {\bf x}}$ and ${\bf h}' = {\bf h}$ . Therefore, ${\widehat m}({\widetilde {\bf x}},{\bf h})$ is a filtered version of the original image $m({\widetilde {\bf x}},{\bf h})$ . They have exactly the same kinematics. When the migration velocity is different from the modeling velocity ( ${\bf v}_0 \neq {\bf v}$ ), the two images may substantially differ. Because we want to use migration results to estimate velocity, it is important to demonstrate that the velocity information contained in the prestack image obtained from the data modeled using the proposed procedure is consistent with the velocity information extracted from the prestack image obtained from migrating the originally recorded data set.

Throughout this paper, we perform numerical examples using Green’s functions computed by means of one-way wavefield extrapolation (Biondi, 2002; Stoffa et al., 1990; Claerbout, 1985; Ristow and Rühl, 1994). Although not tested here, Green’s functions obtained using other methods, such as by solving the two-way wave equation, can also be used under this framework. Since the one-way wavefield extrapolator has limited accuracy for large-angle propagation, we only compute horizontal half subsurface offset, i.e., ${\bf h}=(h_x,h_y,0)$ .

Figure 1(a) shows a modified Sigsbee2A velocity model, which contains two square anomalies below the salt: one $10\%$ lower and the other $10\%$ higher than the sediment velocity. We model $268$ shots using the two-way wave-equation with a time-domain finite-difference scheme. The data is recorded with a marine acquisition geometry, the maximum offset for each shot is about $26000$ ft. We refer this two-way data set as ``original data'' hereafter. The migrated prestack image and gathers using the original data and a starting velocity model (Figure 1(b)) are shown in Figure 2. The amplitude of the background image has been normalized by the diagonal of the phase-encoded Hessian (equation 3) to partially compensate for uneven illumination (Tang, 2009). Note the unfocused subsurface-offset-domain common-image gathers (SODCIGs) due to velocity errors. Then we use the target image (Figure 2) and the Born modeling described above to generate $41$ plane-wave-source gathers from depth level ${\widetilde z}_s=15500$ ft, where the take-off angle is from $-45^{\circ}$ to $45^{\circ}$ . The same starting velocity model that was used for migration (Figure 1(b)) has been used for modeling, and the new data set is collected just above the target region, i.e., ${\widetilde z}_r={\widetilde z}_s=15500$ ft. We refer to the Born-modeled data as ``new data'' hereafter. The plane-wave migration result of the new data set using the starting velocity model is shown in Figure 4. Note the same kinematics shown in Figures 2 and 4. This suggests that the velocity information has been successfully preserved using the new data set, which is substantially smaller compared to the original one. Also note that Figure 4 is more blurry than Figure 2 due to the filtering effect of the resolution function $\Delta G$ (equation 8).

bwi-sigsb2c-vmod,bwi-sigsb2c-bvel
Figure 1. (a) The modified Sigsbee2A velocity model and (b) the starting velocity model. [ER]

bwi-sigsb2c-bimg-target-cpst
Figure 2. Initial image obtained using the original two-way finite-difference modeled data. The image has been normalized using the diagonal of the phase-encoded Hessian. Top: zero-subsurface offset image (stacked image), bottom: SODCIGs at surface location $38575$ ,$40000$ ,$41425$ ,$49300$ ,$50000$ and $50700$ ft, respectively. [CR]

bwi-sigsb2c-bhes-target
Figure 3. The diagonal of the phase-encoded Hessian obtained using the background velocity model. View descriptions are the same as in Figure 2. [CR]

bwi-sigsb2c-bimg-born-planes-target
Figure 4. The migrated image and gathers using the new data set. View descriptions are the same as in Figure 2. [CR]

Target-oriented wavefield tomography using demigrated Born data