Riemannian Wavefield Extrapolation

Performing wavefield extrapolation on topographic computational meshes computed through conformal mapping requires parameterizing the acoustic wave-equation by a set of variables that describe the coordinate system. In 2-D, we denote these variables the extrapolation direction, $\t$, (analogous to depth in Cartesian wavefield extrapolation), and the direction orthogonal, $\gamma$ (analogous to horizontal offset in Cartesian wavefield extrapolation). Variables $\t$ and $\gamma$ are related to the topographic coordinate system point set through $(\t,\gamma)=(\Re(z_{topo}^{cs}),\Im(z_{topo}^{cs}))$.Figure  presents a sketch of the topographic coordinate system geometry.

 

topocoord Figure 4 Cartoon illustrating the topography coordinate system. Variable $\t=\t(x,z)$ is the extrapolation direction and parameter $\t^t$ may be considered a topographic ``front''. Variable $\gamma=\gamma(x,z)$ is the coordinate across the extrapolation step at a constant $\t$ step, and parameter $\gamma^g$ may be considered a topographic ``ray''.

The 2-D acoustic wave-equation for wavefield, ${\cal U}$, at frequency, $\omega$, governing propagation in topographic coordinates is Sava and Fomel (2004),  

$$\frac{1}{\alpha J} \left[ \frac{\partial}{\partial\t} \left(... ... U}}{\partial\gamma} \right) \right]= -\omega^2 s^2 {\cal U},$$ (3)

where s is the slowness of the medium, $\alpha$ a distance scaling factor in the extrapolation direction $\t$, and J a Jacobian of transformation of coordinate $\gamma$ (analogous to a geometrical ray spreading factor). Parameters $\alpha$ and J are defined by,

$$\begin{eqnarray} \alpha=& \left[ \frac{\partial x}{\partial\t} \frac{\partial x}... ...\frac{\partial z}{\partial\gamma} \right]^{\frac{1}{2}}, \nonumber\end{eqnarray}$$ (4)

where x and z are the coordinates of the underlying Cartesian basis. Note that parameters $\alpha$ and J are solely components of the coordinate system and are independent of the extrapolated wavefield values.

Analogous to wavefield continuation on a Cartesian mesh, a dispersion relation must be specified that forms the basis for all derived extrapolation operators in a topographic coordinate system. The relation being sought is the wavenumber along the extrapolation direction, $k_\tau$. Following Sava and Fomel (2004), the partial derivative operators in (3) are expanded out to generate a second-order partial differential equation with non-zero cross derivatives. Fourier-domain wavenumbers are then substituted for the partial differential operators acting on wavefield, ${\cal U}$, and the quadratic formula is applied to yield the expression for $k_\tau$,  

$$k_\tau= \frac{{\rm i}\alpha}{2 J} \frac{\partial}{\partial \... ...amma- \frac{\alpha^2}{J^2} k_\gamma^2 \right]^{\frac{1}{2}}.$$ (5)

One relatively straightforward way to apply wavenumber $k_\tau$ in an extrapolation scheme is to develop the topographic coordinate system equivalent to a phase-screen extrapolation operator Sava (2004). In the following example, we treat solely the kinematic, one-way propagation of recorded wavefields. This asymptotic approximation leads us to drop the first order partial differential terms in (5),  

$$k_\tau=\pm \sqrt{a^2 \omega^2 - b^2 k_\gamma^2},$$ (6)

where $a=s\alpha$ and $b=\alpha/J$. The expansion of $k_\tau$ about reference parameters a₀ and b₀ is,

$$k_\tau\approx \kkt0 + \left. \frac{\partial k_\tau}{\partial... ..._\tau}{\partial b} \right\vert _{a_0,b_0} \left( b-b_0\right),$$ (7)

where subscript denotes reference. Partial derivatives with respect to parameters a and b are,

$$\begin{eqnarray} \left.\frac{\partial k_\tau}{\partial a}\right\vert _{a_0,b_0} ... ...ga\frac{b_0}{a_0}\left(\frac{k_\gamma}{\omega}\right)^2, \nonumber\end{eqnarray}$$ (8)

where the square root function in the denominator has been expanded using a Padé approximation. The choice of numerical constants $c_1=\frac{1}{2}$ and c₂=0 yields a 15 finite-difference term. Thus, the phase-screen approximation for extrapolation wavenumber, $k_\tau$, is,

$$\begin{eqnarray} k_\tau\approx \kkt0 + \omega\left(a-a_0\right) + \omega\frac{ \... ...(\frac{b_0}{a_0}\right)^2\left(\frac{k_\gamma}{\omega}\right)^2}. \end{eqnarray}$$ (9)

This expression can be generalized to include multiple reference media through a phase-shift plus interpolation (PSPI) approach Gazdag and Sguazzero (1984) over the two parameters; however, this extension is not treated here. The approximation for wavenumber, $k_\tau$, given in (9) is used in a conventional wavefield extrapolation scheme that extends the recorded wavefield away from the acquisition surface to the required subsurface locations. This involves solving a one-way wave-equation which, in discrete extrapolation steps of $\Delta \t$, requires a recursive computation of the following:

$$\mathcal{U}(\t+\Delta \t,\gamma,\omega)=\mathcal{U}(\t,\gamma,\omega)\;{\rm e}^{{\rm i}k_\tau\Delta \t}.$$ (10)

Our prestack migration example is computed using a shot profile migration code. This involves extrapolating the source and receiver wavefields, $\S$ and ${\cal R}$, independently using,

$$\begin{eqnarray} \S_{\t+\Delta \t} = & \S_{\t} {\rm e}^{-{\rm i}k_\tau\Delta \t}... ...\t} = & {\cal R}_{\t} {\rm e}^{ {\rm i}k_\tau\Delta \t}, \nonumber\end{eqnarray}$$ (11)

and applying an imaging condition at each extrapolation level to generate image, ${\cal I}(\t,\gamma)$,

$${\cal I}(\t,\gamma)= \sum_{i} \sum_{w} \S(\t,\gamma,\omega;{... ...overline{ {\cal R} \left(\t,\gamma,\omega;{\bf s_i}\right) },$$ (12)

where the line over the receiver wavefield indicates complex conjugate. Image ${\cal I}(\t,\gamma)$ is then mapped to a Cartesian coordinate system using sinc-based interpolation operators in the neighborhood of each mapped point to generate the final image, ${\cal I}(x,z)$.