First-break traveltime tomography with the double-square-root eikonal equation

Appendix A: Causal discretization of DSR eikonal equation

To simplify the analysis, we consider first the DSR branch as shown in Figure 1 and described by equation 1. We assume a rectangular 2-D velocity model $v (z,x)$ and thus a cubic 3-D prestack volume $T (z,r,s)$ with $r$ and $s$ axes having the same dimension as $x$. After an Eulerian discretization of both $v$ and $T$, we denote the grid spacing in $z$ as $\Delta$, and in $x$, $r$ and $s$ as $\delta$.

update1
Figure 17. An implicit discretization scheme. The arrow indicates a DSR characteristic. Its root is located in the simplex $T^z T^r T^s$.

In Figure A-1, we study the traveltime at grid point $\mathbf{y} = (z,r,s)$ and its relationship with neighboring grid points $\mathbf{y}^z = (z+\Delta,r,s)$, $\mathbf{y}^r = (z,r-\delta,s)$ and $\mathbf{y}^s = (z,r,s+\delta)$ with a semi-Lagrangian scheme. According to the geometry in Figure 1, in the $(z,r,s)$ space the DSR characteristic (Duchkov and de Hoop, 2010) is straddled by $\mathbf{y}^z \mathbf{y}^r \mathbf{y}^s$.

In order to compute $T$, we could continue along this characteristic up until its intersection with the simplex $\mathbf{y}^z \mathbf{y}^r \mathbf{y}^s$. Suppose the intersection point is $\mathbf{\tilde{y}} = (\tilde{z},\tilde{r},\tilde{s})$ and $\alpha_i$'s are its barycentric coordinates, i.e.

$$\alpha_z,\alpha_r,\alpha_s \in [0,1];\,\,\, \alpha_z + \alph... ...mathbf{y}^z + \alpha_r \mathbf{y}^r + \alpha_s \mathbf{y}^s\;.$$ (18)

This leads to the following discretization:

$$T \equiv T (\mathbf{y}) = \smash{\displaystyle\min_{\mathbf{... ...\sqrt{(z-\tilde{z})^2 + (s-\tilde{s})^2}}{v (z,s)} \right\}\;.$$ (19)

Here we further assume that $v (z,r)$ and $v (z,s)$ are locally constant and that ray-segments between $z$ and $z + \Delta$ are well approximated by straight lines. This also means that a linear interpolation in $T$ within the simplex $\mathbf{y}^z \mathbf{y}^r \mathbf{y}^s$ is exact, i.e. $T (\mathbf{\tilde{y}}) = \alpha_z T^z + \alpha_r T^r + \alpha_s T^s$, where $T^i = T (\mathbf{y}^i)$ for $i = z,r,s$. The minimization over all possible intersection points in equation A-2 guarantees a first-arrival traveltime.

Defining the ratio in grid spacing as $\mu \equiv \Delta / \delta$ and denoting $v_r = v (z,r)$ and $v_s = v (z,s)$, equation A-2 can be re-written with the barycentric coordinates in A-1 as

$$T = \smash{\displaystyle\min_{(\alpha_z,\alpha_r,\alpha_s)}}... ...{\delta}{v_s} \sqrt{\alpha_s^2 + \mu^2 \alpha_z^2} \right\}\;.$$ (20)

Figure A-1 is based on a particular direction of the diving wave: rightward from the source and leftward from the receiver, as in Figure 1. This yields the above positions of $\mathbf{y}^s$ and $\mathbf{y}^r$, and the formula A-3 becomes an update from the $\mathbf{y}^z \mathbf{y}^r \mathbf{y}^s$ quadrant. Since in general the direction of a diving wave is not a priori known, we compute one such update from each of the lower quadrants and take the smallest amongst them as a value of $T$.

To explore the causal properties of equation A-3, we first assume that the minimum is attained at some $\mathbf{\tilde{y}}^* = \xi_z \mathbf{y}^z + \xi_r \mathbf{y}^r + \xi_s \mathbf{y}^s$ such that $\xi_i > 0$ for $i = z,r,s$. From the Kuhn-Tucker optimality conditions (Kuhn and Tucker, 1951), there exists a Lagrange multiplier $\lambda$ such that

$$\lambda = T^z + \delta \left( \frac{\mu^2 \xi_z}{v_r \sqrt{... ...ac{\mu^2 \xi_z}{v_s \sqrt{\xi_s^2 + \mu^2 \xi_z^2}} \right)\;;$$ (21)

$$\lambda = T^i + \delta \left( \frac{\xi_i}{v_i \sqrt{\xi_i^2 + \mu^2 \xi_z^2}} \right)\;;\,\,\,i = r,s\;.$$ (22)

Taking a linear combination of A-4 and A-5 to match the right-hand side of A-3, we find that $\lambda = T$ and thus

$$T - T^z = \delta \left( \frac{\mu^2 \xi_z}{v_r \sqrt{\xi_r^... ...mu^2 \xi_z}{v_s \sqrt{\xi_s^2 + \mu^2 \xi_z^2}} \right) > 0\;;$$ (23)

$$T - T^i = \delta \left( \frac{\xi_i}{v_i \sqrt{\xi_i^2 + \mu^2 \xi_z^2}} \right) > 0\;;\,\,\,i = r,s\;.$$ (24)

This means that if $T$ defined by equation A-3 depends on $T^i$ then $T > T^i$ for $i = z,r,s$, or

$$T > \max (T^z, T^r, T^s)\;$$ (25)

and a Dijkstra-like method (Dijkstra, 1959) is applicable to solve the discretized system.

A direct substitution from equations A-6 and A-7 results in

$$\frac{T-T^z}{\Delta} = \sqrt{\frac{1}{v_r^2} - \left( \frac... ...rt{\frac{1}{v_s^2} - \left( \frac{T-T^s}{\delta} \right)^2}\;.$$ (26)

If $T^z$, $T^r$ and $T^s$ are known, then $T$ can be recovered by solving the 4th degree polynomial equation A-9 and choosing the smallest real root that satisfies condition A-8. This gives a three-sided update at $T$. The discretization is implicitly causal and provides unconditional consistency and convergence.

If there is no real root or none of the real roots satisfy A-8, the minimizer $\mathbf{\tilde{y}}$ can not lie in the interior of simplex $\mathbf{y}^z \mathbf{y}^r \mathbf{y}^s$ and at least one of the $\xi_i$s must be zero. If $\xi_z = 0$, it is easy to show that one of the other barycentric coordinates is also zero and equation A-3 simplifies to

$$T = \min \left( T^r + \frac{\delta}{v_r}, T^s + \frac{\delta}{v_s} \right)\;,$$ (27)

which is a causal discretization of the DSR singularity in equation 3. On the other hand, if $\xi_z \not= 0$ and $\xi_r \not= 0$ but $\xi_s = 0$, i.e. the slowness vector at $s$ is vertical, then

$$T = (\xi_z T^z + \xi_r T^r) + \frac{\delta}{v_r} \sqrt{\xi_r^2 + \mu^2 \xi_z^2} + \frac{\Delta}{v_s} \xi_z\;.$$ (28)

A similar Kuhn-Tucker-type argument shows that equation A-11 is also causal: if $\xi_z,\xi_r > 0$, then $T > \max(T^z,T^r)$. In this case, $T$ can be computed by solving

$$\frac{T-T^z}{\Delta} = \sqrt{\frac{1}{v_r^2} - \left( \frac{T-T^r}{\delta} \right)^2} + \frac{1}{v_s}\;.$$ (29)

Equation A-12 is equivalent to setting $\partial T / \partial s = 0$ in equation 1. Analogously, when $\xi_z \not= 0$ and $\xi_s \not= 0$ but $\xi_r = 0$, we have $\partial T / \partial r = 0$ and

$$\frac{T-T^z}{\Delta} = \frac{1}{v_r} + \sqrt{\frac{1}{v_s^2} - \left( \frac{T-T^s}{\delta} \right)^2}\;,$$ (30)

with the causal solution satisfying $T > \max(T^z,T^s)$. Equations A-12 and A-13 provide a two-sided update at $T$. Finally, if $\xi_z \not= 0$ but $\xi_s = 0$ and $\xi_r = 0$, i.e. $\partial T / \partial s = 0$ and $\partial T / \partial r = 0$, we obtain the simplest one-sided update:

$$\frac{T-T^z}{\Delta} = \frac{1}{v_r} + \frac{1}{v_s}\;.$$ (31)

We note that the one-sided update A-14 could be considered a special case of two-sided updates: if $T= T^r$ (or $T^s$), then A-14 becomes equivalent to A-12 (or, respectively, A-13). Similarly, the two-sided updates can be viewed as special versions of the three-sided one: e.g., if $T = T^r = \max(T^z,T^r,T^s)$, then A-9 becomes equivalent to A-13. This means that the causal criteria for formulas A-9, A-12 and A-13 can be relaxed (the inequalities do not have to be strict). This relaxation is used to streamline the update strategy in the Numerical Implementation section.

In Figure A-1 and the corresponding semi-Lagrangian discretization A-2, the ray-path is linearly approximated up to its intersection with the simplex $\mathbf{y}^z \mathbf{y}^r \mathbf{y}^s$ at a priori unknown depth $z + \xi_z \Delta$. An alternative explicit semi-Lagrangian discretization can be obtained in the spirit of Figure 1 by tracing the ray up to the pre-specified depth $z + \Delta$. In Figure A-2, we consider the DSR characteristic being straddled by $\mathbf{y}_*^z \mathbf{y}_*^r \mathbf{y}_*^s$, where $\mathbf{y}_*^z = (z+\Delta,r,s)$, $\mathbf{y}_*^r = (z+\Delta,r-\delta,s)$ and $\mathbf{y}_*^s = (z+\Delta,r,s+\delta)$. Denoting $\mathbf{\tilde{y}}_* = (z+\Delta,\tilde{r}_*,\tilde{s}_*)$ for the intersection point between DSR characteristic and simplex $\mathbf{y}_*^z \mathbf{y}_*^r \mathbf{y}_*^s$, we obtain the following discretization:

$$T = \smash{\displaystyle\min_{\mathbf{\tilde{y}_*} \in \math... ...frac{\sqrt{\Delta^2 + (s-\tilde{s}_*)^2}}{v (z,s)} \right\}\;.$$ (32)

One could perform the same analysis of A-3 through A-9 to equation A-15. For the sake of brevity, we omit the derivation and show the resulting explicit discretization scheme:

$$\frac{T-T_*^z}{\Delta} = \sqrt{\frac{1}{v_r^2} - \left( \fr... ...ac{1}{v_s^2} - \left( \frac{T_*^z-T_*^s}{\delta} \right)^2}\;,$$ (33)

where $T_*^i = T (\mathbf{y}_*^i)$ for $i = z,r,s$. More generally, to account for various possible directions of the diving wave, we can set $T_*^r = \min \left(T (z+\Delta,r-\delta,s), T (z+\Delta,r+\delta,s)\right)$ and $T_*^s = \min \left(T (z+\Delta,r,s-\delta), T (z+\Delta,r,s+\delta)\right)$.

update2
Figure 18. An explicit discretization scheme. Compare with Figure A-1. The arrow again depicts a DSR characteristic with its root confined in the simplex $T_*^z T_*^r T_*^s$.

Compared with equation A-9, equation A-16 does not require solving a polynomial equation. Moreover, $T$ depends only on values in lower z-slices, which means that the system of equations can be solved in a single sweep in the $-z$ direction. Unfortunately, despite this efficiency on a fixed grid, the explicit discretization has a major disadvantage stemming from the requirement that the characteristic should be straddled by $y^z_* y^r_* y^s_*$. This imposes an upper bound on $\mu$ based on the slope of the diving wave. Moreover, since every diving ray is horizontal at its lowest point, the convergence is possible only if $\mu \to 0$ under mesh refinement. In practice, this means that the results are meaningful only if $\Delta$ is significantly smaller than $\delta$. We note that restrictive stability conditions also arise for time-dependent Hamilton-Jacobi equations of optimal control, where sufficiently strong inhomogenieties can make nonlinear/implicit schemes preferable to the usual linear/explicit approach (Vladimirsky and Zheng, 2013).

The above analysis also applies to the first branch of the DSR eikonal equation in Figure 2. However, in the discretized $(z,r,s)$ domain, the slowness vectors at $s$ and $r$ are always aligned in the $z$ direction, either upward or downward. For this reason, there is no DSR characteristic that accounts for the second and third scenarios. We will refer to the first and last branches in Figure 2 as causal branches of DSR eikonal equation, and the left-over two as non-causal branches.

search
Figure 19. When slowness vectors at $s$ and $r$ are pointing in the opposite directions, the ray-path must intersect with line $s-r$ at certain point $q$.

Note that when the slowness vectors at $s$ and $r$ are pointing in opposite directions, there must be at least one intersection of the ray-path with the $z$ depth level in-between. As shown in Figure A-3, ray segments between these intersections fall into the category of causal branches. Thus a search process for the intersections is sufficient in recovering the non-causal branches during forward modeling. Moreover, because we are interested in first-breaks only, the minimum traveltime requirement allows us to search for only one intersection, such as $q$ denoted in Figure A-3:

$$T (z,r,s) = \smash{\displaystyle\min_{q \in (s,r)}} \left\{ T (z,q,s) + T (z,r,q) \right\}\;.$$ (34)

Other possible intersections in intervals $(s,q)$ and $(q,r)$ have already been recovered when computing $T (z,q,s)$ and $T (z,r,q)$, as long as we enable the intersection searching from the beginning of forward modeling. The traveltime of non-causal branches from equation A-17 is then compared with that from causal branches, and the smaller one should be kept.

Unfortunately, this search routine induces considerable computational cost. Moreover, we note that, under a dominant diving waves assumption, the first DSR branch, despite being causal, becomes useless if A-17 is turned-off.

First-break traveltime tomography with the double-square-root eikonal equation