Potential function ray tracing

The next processing step involves tracing geometric coordinate system rays from the generated PF. The goal here is to develop an orthogonal, ray-based coordinate system related to an underlying Cartesian mesh through one-to-one mappings,

$$\begin{eqnarray} \tau =& \tau (z,x,y) & \; \nonumber \\ \gamma=& \gamma (z,x... ...partial (z,x,y) } \neq 0,\\ \eta =& \eta (z,x,y) & \; \nonumber \end{eqnarray}$$ (7)

where $\tau$ is the wavefield extrapolation direction (equivalent to z in Cartesian), $\gamma$ and $\eta$ are the two orthogonal directions (equivalent to x and y in Cartesian), and J is the Jacobian of the coordinate system transformation. Recorded wavefield, $U(\tau_0,\gamma,\eta)$, is extrapolated from the acquisition surface defined by $\tau_0$ into the subsurface along the rays coordinate system defined by triplet $\left[\tau,\gamma,\eta\right]$.

Geometric rays are traced by solving a first-order ordinary differential equation through integrating the PF gradient field along the gradient direction,  

$$\delta \phi = \int_{a(z_0,x_0,y_0)}^{b(z_1,x_1,y_1)} \vert... ...tial x}+ {\rm d}y \frac{\partial \phi} {\partial y} \right],$$ (8)

where a(z₀,x₀,y₀) is a known lower integration bound at equipotential $\phi(a)$, and b(z₁,x₁,y₁) is an unknown upper integration bound located on equipotential surface, $\phi(b)$, and $\vert\vert\nabla \phi\vert\vert$ is the L² norm of the gradient function. The only unknown parameter is b(z₁,x₁,y₁); hence, Equation (8) is an integral equation with an unknown integration bound. This approach is similar to phase-ray tracing method described in Shragge and Sava (2004b); however, in this case the integration step lengths are now unknown. Note also that the equipotentials of the upper and lower bounding surfaces in Equation (4) require PF steps of $\delta \phi=1/N$.

The following approach locates unknown integration bound, b, on the next equipotential:

1.

Numerically integrate Equation (8) on the interval [$a,a+\delta a$] where $\delta a$ is smaller than the expected step size, and test to see whether $\phi(a)-\phi(a+\delta a) \gt \delta \phi$; if yes goto step 3;

2.

Numerically integrate Equation (8) on next interval $[a+\delta a,a+2\delta a]$ and test whether $\phi(a)-\phi(a+2\delta a) \gt \delta \phi$; if yes goto step 3; if no, repeat step 2 n times until true;

3.

Linearly interpolate between points $a+(n-1)\delta a$ and $a+\delta a$ to find the b corresponding to $\phi(b) = \phi(a)-\delta \phi$.

A geometric ray is initiated at a particular $[\gamma_0, \eta_0]$ on acquisition surface defined by $\tau_0$, and computed by integrating through each successive $\delta \phi$ step until the lower bounding surface $\tau_N$ is reached. This procedure is repeated for all $\gamma$ and $\eta$ acquisition points.