3D two-point ray tracing for heterogeneous, weakly transversely-isotropic media

TL;DR

This paper introduces a fast 3D two-point ray-tracing method for heterogeneous, weakly transversely-isotropic media using Chebyshev approximation and conjugate gradient optimization, enabling efficient and accurate traveltime calculations.

Contribution

It presents a novel ray-tracing technique combining Chebyshev approximation with conjugate gradient methods for complex anisotropic media.

Findings

High accuracy in traveltime computation with few polynomial terms

Rapid convergence of the conjugate gradient method

Effective handling of smoothly heterogeneous, anisotropic media

Abstract

A two-point ray-tracing technique for 3D smoothly heterogeneous, weakly transversely-isotropic media is based on Fermat's principle and takes advantage of global Chebyshev approximation of both the model and curved rays. This approximation gives explicit derivatives of traveltimes with respect to ray parameters and allows one to use the rapidly converging conjugate gradient method to compute traveltimes. The method is fast because, for smoothly heterogeneous media, approximation of rays by only a few polynomials and a few conjugate gradient iterations provide excellent precision of traveltime calculation.

Tables

Table 1. Table 1: Computation times and accuracy of the traveltime calculations for the 13 rays in Figure 2 on a Sun-4 workstation, as a function of M 𝑀 M (the number of Chebyshev polynomials) and N 𝑁 N (the number of points used in integration).

		Maximum error	Computation
$M$	$N$	in traveltime (ms)	time (s)
3	4	6.46	2.13
	5	2.80	2.79
	6	0.45	3.06
	7	0.12	3.76
	8	0.12	4.05
	9	0.12	4.90
4	4	4.82	4.04
	5	1.41	4.50
	6	0.37	5.34
	7	0.10	6.01
	8	0.05	6.75
	9	0.02	7.65
5	4	0.98	4.85
	5	0.81	6.75
	6	0.06	7.50
	7	0.04	9.16
	8	0.02	10.20
	9	0.01	11.27

Equations

$$m_1(\mathbf{x})$$

$$m_2(\mathbf{x})$$

$$m_3(\mathbf{x})$$

$$m_4(\mathbf{x})$$

$$m_5(\mathbf{x})$$

$$m_6(\mathbf{x})$$

$$m_7(\mathbf{x})$$

$$\mathbf{c}(\mathbf{x})$$

$$c_3(\mathbf{x})$$

$$\mathbf{x}\in [\mathbf{a},\mathbf{b}].$$

$$\mathbf{x}(s)\equiv x_i(s)=\sum _{k=1}^{M_i}{r}_{i,k}{T}_{k-1}(s),(s\in [0,1];i=1,2,3),$$

$$\mathbf{x}(0)$$

$$\mathbf{x}(1)={\mathbf{x}}^{receiver},$$

$$t_{\mathrm{Q}}={\int }_{{\mathbf{x}}^{source}}^{{\mathbf{x}}^{receiver}}{p}_{\mathrm{Q}}(m_{\eta }(\mathbf{x}),\mathbf{x})d\mathbf{x},(\mathrm{Q}=\mathrm{P},\mathrm{SV}\mathrm{or}\mathrm{SH}),$$

$$t_{\mathrm{Q}}={\int }_0^1{\mathcal{T}}_{\mathrm{Q}}({\mu }_{\eta },\mathbf{r},s)ds,$$

$${\mathcal{T}}_{\mathrm{Q}}({\mu }_{\eta },\mathbf{r},s)=p_{\mathrm{Q}}(m_{\eta }(\mathbf{r},s),\mathbf{r},s)R(\mathbf{r},s),$$

$$R(\mathbf{r},s)={\left[\sum _{i=1}^3{\dot{x}}_i^2(s)\right]}^{1/2}$$

$$\frac{\partial t_{\mathrm{Q}}}{\partial r_{i,l}}={\int }_0^1{\mathcal{D}}_{\mathrm{Q},i,l}({\mu }_{\eta },\mathbf{r},s)ds=0,(\mathrm{Q}=\mathrm{P},\mathrm{SV}\mathrm{or}\mathrm{SH};i=1,2,3;l=3,\dots ,M_i),$$

$$t_{\mathrm{Q}}=\mathbf{P}\cdot \mathbf{C}({\mathcal{T}}_{\mathrm{Q}})$$

$$\frac{\partial t_{\mathrm{Q}}}{\partial r_{i,l}}=\mathbf{P}\cdot \mathbf{C}({\mathcal{D}}_{\mathrm{Q},i,l}),$$

$$V(x_3)=V_0\exp (g{x}_3),$$

$$t(x_1,x_3)=\frac{\sqrt{2(\cosh g{x}_3-\cos g{x}_1)}}{g{V}_0\exp (g{x}_3/2)}$$

$$ϵ$$

$$\delta$$

$$\gamma$$

$$c_1$$

$$c_2$$

$$y_i=1+\frac{x_i-b_i}{b_i-a_i},(i=1,2,3).$$

$$x_i=a_i+(b_i-a_i)y_i,(i=1,2,3).$$

$$T_{\mathbf{k}-1}(\mathbf{y})\equiv T_{k_i-1}(y_i),(i=1,2,3;k_i=1,\dots ,N_i)$$

$$m({\bm{x}})\equiv m(x_{1},\,x_{2},\,x_{3})=\sum_{k_{1}=1}^{N_{1}}\sum_{k_{2}=1}^{N_{2}}\sum_{k_{3}=1}^{N_{3}}\biggl{(}\mu(k_{1},\,k_{2},\,k_{3})\prod_{i=1}^{3}[T_{k_{i}-1}(y_{i})]\biggr{)}.$$

$$\mu({\bm{k}})\equiv\mu(k_{1},\,k_{2},\,k_{3})={\biggl{(}\prod_{i=1}^{3}N_{i}\biggr{)}}^{\!\!-1}\sum_{j_{1}=1}^{N_{1}}\sum_{j_{2}=1}^{N_{2}}\sum_{j_{3}=1}^{N_{3}}\biggl{(}m(\chi_{j_{1}},\,\chi_{j_{2}},\,\chi_{j_{3}})\prod_{i=1}^{3}[T_{k_{i}-1}(\lambda_{j_{i}})]\biggr{)},$$

$${\lambda }_{j_i}=\frac{1}{2}\left(1+\cos \frac{2j_i-1}{2N_i}\pi \right),$$

$${\chi }_{j_i}=a_i+(b_i-a_i){\lambda }_{j_i},(i=1,2,3;j_i=1,\dots ,N_i).$$

$$x_i(s)=\sum _{k=1}^{M_i}{r}_{i,k}{T}_{k-1}(s),(s\in [0,1];i=1,2,3).$$

Full text

I Introduction

Elastic anisotropy is widespread in the earth, and many papers discuss the influence of velocity anisotropy on traveltimes and amplitudes of seismic waves. A review of various anisotropic phenomena and an extensive list of references can be found in Crampin and LovellCrampinLovell1991 . It is difficult to separate anisotropy from heterogeneity because their influences on seismic wavefields resemble each another. Even shear-wave splitting, commonly attributed to anisotropy, can be caused by strong heterogeneityGrechkaMcMechan1995 . One way of measuring seismic anisotropy is cross-well traveltime tomographyPrattChapman1992 ; Michelena1993 ; Michelenaetal1993 that, under certain conditions, allows one to discriminate the effects of anisotropy and heterogeneity on recorded traveltimes. Because anisotropy is essentially a 3D phenomenon, anisotropic 3D ray tracing is an essential tool for its investigation.

Most of the existing ray-tracing approaches fall into one of three groups: methods based on finite-difference solution of the eikonal equationVidale1988 ; Dellinger1991 ; vanTrierSymes1991 , shooting methodsCervenyMolotkovPsencik1977 ; Langanetal1985 ; GajewskiPsencik1987 ; VirieuxFarra1991 ; GuestKendall1993 , and bending methodsUmThurber1987 ; Protheroetal1988 ; Schneideretal1992 ; FariaStoffa1994 , derived, respectively, from Huygens’ principle, Snell’s law, and Fermat’s principle. Methods based on finite-difference solution of the eikonal equation have been applied to 2D anisotropic traveltime calculationDellinger1991 ; QinSchuster1993 . Shooting and bending methods are used for 3D ray tracingUmThurber1987 ; Protheroetal1988 ; VirieuxFarra1991 in isotropic media; only shooting methods have currently been applied for ray tracing in 3D heterogeneous anisotropic mediaCerveny1972 ; GajewskiPsencik1987 ; GuestKendall1993 .

We develop a ray-bending technique because two-point ray tracing is more convenient for traveltime tomography than shooting methods. We use global Chebyshev approximation of a heterogeneous anisotropic model and curved rays that make the computations fast.

II Methodology

The main idea of the proposed ray-tracing method is that, in smoothly heterogeneous media, rays are smooth curves that can be approximated by smooth basis functions. We choose Chebyshev polynomials as the basis functions because series of these polynomials usually converge more rapidly to the approximated function than any other polynomial-based seriesLanczos1988 . Fermat’s principle is employed for two-point ray tracing.

To parameterize smoothly heterogeneous media, we use a global 3D approximation of the medium, defined as a sum of Chebyshev polynomials and described by 3D Chebyshev polynomial coefficients or Chebyshev spectral components. The main advantage of this kind of parametrization is that it provides explicit expressions for traveltimes, and, if a ray has already been traced, explicit relations for variations of the traveltime as a function of Chebyshev spectral components of all model and ray parameters. Thus, the traveltime derivatives are obtained at almost no additional computation cost, allowing us to apply the conjugate gradient methodPressetal1987 to bend a ray path and minimize the traveltime along it. The absence of the derivatives would entail the use of slower converging minimization methods, such as the Nelder-Mead searchPressetal1987 ; Protheroetal1988 .

The proposed technique possesses both the advantages and disadvantages of two-point ray tracing. In some models more than one ray (and traveltime) may exist between a source-receiver pair. Finding all valid solutions is generally difficult because they may correspond not only to minima but also to maxima or saddle points of the traveltime. The same pertains to anisotropy; cusps at the group velocity surfaces, as well as heterogeneity, can produce more than one ray path between the two end points. Although all these rays may be found by solving special equationsGrechkaObolentseva1993 or by iteration for different initial guessesObolentsevaGrechka1988 , we shall restrict our analysis to weakly transversely isotropic media to assure the absence of cusps and the presence of explicit relations for the group velocitiesThomsen1986 ; Byunetal1989 .

III Two-point Chebyshev ray-tracing

III.1 Model and ray parametrization

We study kinematics of wave propagation in a heterogeneous transversely isotropic (TI) model occupying a 3D rectangular volume, specified by the Cartesian coordinates of its corners ${\bm{a}}=(a_{1},\,a_{2},\,a_{3})$ and ${\bm{b}}=(b_{1},\,b_{2},\,b_{3})$. We define the functions

[TABLE]

corresponding to Thomsen’sThomsen1986 anisotropy parameters of TI model, and the directional cosines

[TABLE]

defining the orientation of the unit vector of the symmetry axis

[TABLE]

Vector ${\bm{x}}$ in equations (1) and (2) denotes a point ${\bm{x}}\equiv(x_{1},\,x_{2},\,x_{3})$, belonging to the model volume

[TABLE]

Parameters $m_{\eta},~{}(\eta=1,\ldots,7)$ specify a weakly TI medium inside the volume. For our objectives, however, it is more convenient to use the 3D Chebyshev spectral components $\mu_{\eta}$ (Appendix A) instead of the functions $m_{\eta}({\bm{x}})$.

Consider the ray ${\bm{x}}(s)$ connecting a source at ${\bm{x}}^{source}$ and a receiver at ${\bm{x}}^{receiver}$ (Figure 1). The Cartesian coordinates $x_{i}$ of the ray are expanded in the series of Chebyshev polynomials $T_{k}(s)$

[TABLE]

where $s$ is the normalized ray length, defined such that

[TABLE]

$r_{i,k}$ are the Chebyshev spectral components of the ray, and $M_{i}$ are the numbers of the polynomials used to approximate the ray in each direction $i=1,\,2,\,3$ (Appendix B).

Fermat’s principle allows us to compute the Chebyshev spectral components $r_{i,k}$ of a ray.

III.2 Traveltime computation

The traveltime along a curved ray connecting the source ${\bm{x}}^{source}$ and the receiver ${\bm{x}}^{receiver}$ is defined as an integral

[TABLE]

where $p_{\rm Q}$ is the group slowness (reciprocal to the group velocity) of P, SV, or SH waves along the ray (Appendix C), and $m_{\eta}$ are the model parameters [equations (1) and (2)]. Using the ray parametrization given by equations (5) and (6), we rewrite the integral (7) as

[TABLE]

where the integrand

[TABLE]

and

[TABLE]

is the length of an element of the ray arc. The dot over a function denotes a derivative with respect to its argument, and derivatives $\dot{x}_{i}$ are determined by equation (B8) in Appendix B, where the rays are described in terms of Chebyshev polynomials.

Fermat’s principle,

[TABLE]

is used to find the unknown ray spectral coefficients $r_{i,l}$. The index $l$ starts from 3 because the known coordinates of the source and the receiver determine the linear components (for $l=1,\,2$) of the ray (see Appendix B for details). The integrand ${\cal D}$ is given by equation (D1) in Appendix D.

Coefficients $r_{i,l}$ that satisfy equations (11) define rays corresponding to minima, maxima, or saddle points of traveltime. For complicated models, the solution of equations (11) is often nonunique and several rays, connecting the same source-receiver pair, may exist. Tracing all such rays is a complicated computational problem because there are no general methods of solving nonlinear systems like equations (11). However, if the solution of system (11) is unique or we want to find only the fastest ray, we can search for the traveltime minimum instead of solving nonlinear equations (11). The conjugate gradient methodPressetal1987 is applied to find the minimum traveltime. Partial derivatives (11) are the components of the traveltime gradient.

To implement this approach we need to compute integrals (8) and (11). Again, the Chebyshev polynomials are helpful for doing this. Expanding the integrands ${\cal T}$ and ${\cal D}$ as functions of the ray length $s$ in Chebyshev series, we replace the integrals (8) and (11) by the dot products

[TABLE]

and

[TABLE]

where ${\bm{C}}(\cdot)$ denotes the direct Chebyshev transform of its argument, and vector ${\bm{P}}$ is defined by equation (E6) in Appendix E.

III.3 Numerical investigation of Chebyshev ray-tracing

We present two numerical examples to illustrate the features of the proposed technique. The efficiency of the method depends on the choice of the two numbers $M$ and $N$. The value of $M$ [equation (B7)] determines the number of unknown Chebyshev spectral coefficients of the ray to be found with the conjugate gradient method. The value of $N$ [equation (E1)] defines the number of points needed along each ray to compute traveltime with the required precision. These two quantities are to be determined by numerical experiments.

First, we test the method for a model that has a known solution. Consider an isotropic medium with velocity

[TABLE]

increasing exponentially with depth $x_{3}$, for parameters $V_{0}=1.5$ km/s and $g=1.5~{}{\rm km}^{-1}.$ Velocity function (14) was approximated by seven Chebyshev polynomials over the depth interval $x_{3}\in[0,~{}1.2]$, and traveltimes were computed. Figure 2 displays the ray trajectories and calculated traveltimes for $M=4$, $N=7$ compared to the analytic solution

[TABLE]

that exists for this velocity modelGoldin1986 . We also performed computations for a set of values $M$ and $N$ listed in Table 1 to study the accuracy of the traveltime calculations and the time required to compute all 13 rays (Figure 2) on a Sun-4 workstation. Only 1 – 2 iterations of the conjugate gradient method were needed to converge, starting from a straight ray as an initial guess. The number of iterations is found to be independent of $M$ and $N$. Table 1 shows that the computation time depends approximately linearly on $N$ and about quadratically on $M$. Increasing $M$ and $N$ beyond the values presented in Table 1 no longer increases the precision of traveltime calculations because of minor errors associated with the Chebyshev approximation of the velocity function (14). These errors can be made arbitrarily small by using more polynomials to approximate the exponential function in equation (14).

The second example is ray tracing in a 3D heterogeneous TI model described by the anisotropy parametersThomsen1986

[TABLE]

The orientation of the symmetry axis is defined by

[TABLE]

The 3D distribution of the P-wave velocity $\alpha_{0}({\bm{x}})$ is shown in Figure 3, whereas the shear velocity is defined as $\beta_{0}({\bm{x}})=\alpha_{0}({\bm{x}})/2$ [equations (C1) and (C4)]. The other model functions, influencing the traveltimes to a lesser degree, are not shown. The source coordinates (in km) are ${\bm{x}}^{source}=[0.05,\,0.05,\,0.1]$, and 36 receivers are located at the coordinates ${\bm{x}}^{receiver}=[0.2\,i,\,0.2\,j,\,0.0],~{}(i,\,j=0,\,\ldots,\,5)$; these are displayed in Figures 4a and 5a.

We use $M_{1}=M_{2}=M_{3}=5$ polynomials to approximate the rays (i. e., 9 unknown spectral components for each ray are to be found) and $N=9$ points along the rays to perform integration. The conjugate gradient method converges after 2 – 3 iterations. Comparing the results with those in the first example, we expect the accuracy of our traveltime calculations to be around $10^{-4}$ s or better.

Figure 4a presents the P-wave ray trajectories in the specified TI medium and the trajectories in the related isotropic medium, for which $\epsilon=\delta=0$ in equations (16a) and (16b). Corresponding traveltimes are shown in Figure 4b. Although the ray paths differ substantially in the anisotropic and isotropic models, the traveltimes do not.

A similar comparison of ray trajectories and traveltimes for SV and SH waves in the same TI model is given in Figure 5. For this model, the heterogeneity influences traveltimes more than the anisotropy does.

IV Discussion

Two-point ray tracing in 3D heterogeneous anisotropic media may be based on Chebyshev approximation of curved rays. The method is inexpensive in smoothly varying media because only a few polynomials are needed to approximate the ray paths. The global Chebyshev model parametrization provides explicit relations for partial derivatives of traveltimes as functions of the ray parameters. This allows us to use the rapidly converging conjugate gradient method to compute the traveltimes.

The computing time depends directly on the values of $M$ (the number of Chebyshev polynomials approximating a ray) and $N$ (the number of integration points per ray). These quantities were determined in the first example by comparing the traveltimes to the known solution. Generally, when the solution is not known, $M$ and $N$ should be determined by numerical experiments. A few representative rays and traveltimes $t(M,N)$ are computed for a set of values $M$ and $N$. As $M$ and $N$ become greater traveltimes usually gradually decrease, whereas the computation time increases, as illustrated in Table 1. We choose $M$ and $N$ that correspond to the minimum computation time for which traveltimes are about 0.05% greater than the minimum $t(M,N)$ for all tested pairs of $M$ and $N$. The values of $M=5$ and $N=9$ were selected by this procedure for the second example.

Although we applied only the conjugate gradient method to compute the ray trajectories and traveltimes, other methods of solving this problem are available. For example, one could compute the Hessian matrix $\partial^{2}t/\partial\bm{r}^{2}$ and use the second order Newton method, which generally reduces number of iterations.

The proposed technique is fast for weakly TI media, where the group velocity explicitly depends on a ray direction. If anisotropy is not weak, we would have to compute the group velocity for a given ray direction numerically. This can be done iterativelyObolentsevaGrechka1988 and involves searching for one phase angle for TI media and for two phase angles for more generally anisotropic media. This search, done for every evaluation of the group velocity, would significantly increase the computational time. On the other hand, the proposed technique applied to 3D isotropic media is faster by a factor of about five, as becomes clear from the examination of equations (C4) for the group velocities. To find the velocities, some of the quantities $m_{j}~{}(j=1,\,\ldots,\,7)$ have to be computed with equation (A4), which takes most of the computational effort. For an isotropic medium though, only the medium velocity itself needs to be computed.

It is expected that increasing the complexity of the model (especially when functions describing the model parameters are discontinuous) would reduce the computational efficiency. Although Chebyshev ray tracing can be applied to models with discontinuities, the number of polynomials (and perhaps the number of required iterations) would no longer be small, and the technique loses its elegance. Another approach for discontinuous parameters is to separate the medium into blocks along surfaces of discontinuity (interfaces) and to construct rays as sequences of their smooth segments. This idea seems to be feasible but lies beyond the scope of this paper.

V Acknowledgments

The research leading to this paper was supported by a UT-Dallas graduate teaching assistantship (V.G.) and by the sponsors of the UTD Geophysical Consortium (G.A.M.). The computations were performed on a Sun-4 workstation in the Center for Lithospheric Studies at the University of Texas at Dallas. V.G. wishes to thank Sergey Fomel for the idea to perform the Chebyshev ray tracing in 3D rather than in 2D and Ilya Tsvankin for helpful discussions. We thank to Spyros Lazaratos, Reinaldo Michelena and an anonymous reviewer whose comments and suggestions improved the manuscript. This paper is Contribution No. 825 from the Department of Geosciences at the University of Texas at Dallas.

Appendix A 3D Chebyshev transform

The 3D Chebyshev transform is constructed as a generalization of the 1D Chebyshev transformLanczos1988 . A vector ${\bm{x}}\equiv(x_{1},\,x_{2},\,x_{3})$ belonging to the volume ${\bm{x}}\in[{\bm{a}},{\bm{b}}]$ [equation (4)] can be transformed to the unit cube ${\bm{y}}\in[\bm{0},\,\bm{1}]$ by the linear transformation

[TABLE]

The corresponding inverse transform is

[TABLE]

Let us define a set of orthonormal Chebyshev polynomials

[TABLE]

inside the cube ${\bm{y}}\in[\bm{0},\,\bm{1}]$ and express a function $m({\bm{x}})$ as the 3D inverse Chebyshev transform

[TABLE]

The 3D Chebyshev spectral components $\mu({\bm{k}})$ are found by the 3D direct Chebyshev transform

[TABLE]

where $\lambda_{j_{i}}$ are the roots of the Chebyshev polynomials

[TABLE]

and $\chi_{j_{i}}$ relate to $\lambda_{j_{i}}$ through equation (A2), as

[TABLE]

Appendix B Rays and its derivatives

The Chebyshev polynomial approximation of a ray ${\bm{x}}(s)$ is defined by equation (5) as

[TABLE]

The source and receiver conditions (6) give six constraints on components $r_{i,k}$

[TABLE]

rewritten as

[TABLE]

using values of the Chebyshev polynomials at the edges of the segment $[0,\,1]$

[TABLE]

as the boundary conditions. The components $r_{i,1}$ and $r_{i,2}$ are found from equations (B3) and substituted into equation (B1) to give

[TABLE]

where the relations

[TABLE]

were used. Equation (B5) shows that the ray ${\bm{x}}(s)$ is specified by

[TABLE]

quantities $r_{i,k}$, given that the source and receiver conditions (B2) are satisfied.

All the required derivatives can be found by explicit differentiation of equation (B5). These derivatives are

[TABLE]

[TABLE]

and

[TABLE]

The dot over the functions in equations (B8) and (B10) denotes a derivative with respect to the argument.

Appendix C The group slowness in weakly TI and its derivatives

Phase velocities of the P, SV, and SH waves as functions of angle $\theta$ between the wavefront normal and the symmetry axis are given byThomsen1986

[TABLE]

Equations (C1), valid for small $\epsilon$, $\delta$, and $\gamma$ (i.e., for weak anisotropy), provide the linear approximation of the phase velocities as functions of $\epsilon$, $\delta$, and $\gamma$. Also the group velocity equals the phase velocity in the linear approximationThomsen1986 . Thus,

[TABLE]

where $\phi$ is an angle between a ray and the symmetry axis.

In the linear approximation, the relation between the angles $\theta$ and $\phi$ readsThomsen1986

[TABLE]

where $\Delta$, being different for P, SV, and SH waves, is a small quantity that has the same order as $\epsilon$, $\delta$, and $\gamma$. Substituting equations (C1) and (C3) into equation (C2) and neglecting quadratic terms in $\epsilon$, $\delta$, $\gamma$, and $\Delta$, we obtain the group velocities in the linear approximation as

[TABLE]

The same approximation for the group slownesses gives

[TABLE]

where

[TABLE]

The approximation for the square of the group slowness, in a form similar to equation (C5a), was derived earlierByunetal1989 .

The quantity $\psi$ is the square of the dot product of vector ${\bm{c}}$ [equations (2) and (3)] and the normalized tangent to the ray $\dot{\bm{x}}/R$ [equations (B8) and (9)];

[TABLE]

where

[TABLE]

Our next task is to calculate the derivatives $\partial p_{\rm Q}/\partial m_{\eta}$, where Q = P, SV, or SH, and $m_{\eta}$ are defined by equations (1), (2). Differentiating equations (C5), we obtain

[TABLE]

where

[TABLE]

[TABLE]

Because the group slowness (C5) depends directly on $\psi$, which in turn, depends on the Chebyshev spectral components of the ray [equations (10), (B1), (C7) and (C8)], the slowness also depends on these components. The corresponding partial derivatives are

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

derivatives ${\partial p_{\rm Q}/\partial\psi}$ are given by equations (C10), derivatives ${\dot{x}_{i}}$ and ${\partial{\dot{x}_{i}}/\partial r_{i,l}}$ – by equations (B8) and (B10).

Appendix D Computation of integrand ${\cal D}$ in equation (11)

To compute the integrand ${\cal D}_{{\rm Q},i,l}=\partial{\cal T_{\rm Q}}/\partial r_{i,l}$ in equation (11) for anisotropic heterogeneous media, consider the following. First, the group slowness $p_{\rm Q}$ in equation (9) depends on the model parameters $m_{\eta}$ that depend on the spectral components $r_{i,l}$, because the medium is heterogeneous. Second, the group slowness also depends directly on $r_{i,l}$ because of anisotropy (i.e., different $r_{i,l}$ define different ray directions and accordingly different values of the slowness at the same point). Third, the components $r_{i,l}$ are constrained by six equations (B2), determining the components $r_{i,1}$ and $r_{i,2}$. Therefore, we have to find the derivatives $\partial{\cal T_{\rm Q}}/\partial r_{i,l}$ for $l>2.$

Differentiating equation (9) yields

[TABLE]

$({\rm Q=P,SV~{}or~{}SH};~{}i=1,\,2,\,3;~{}l=1,\,\ldots,\,M_{i}),$

where the summation from 1 to 7 and the presence of the second term in the brackets are due to anisotropy. Derivatives ${{\partial p_{\rm Q}}/{\partial m_{\eta}}}$ in (D1) are given by equations (C9), derivatives ${{\partial p_{\rm Q}}/{\partial r_{i,l}}}$ by equation (C12), and derivatives ${{\partial R}/{\partial r_{i,l}}}$ by equation (C14).

The derivatives of the model parameters $m_{\eta}({\bm{x}})$, obtained from equations (A4), are

[TABLE]

$(\eta=1,\,\ldots,\,7;~{}i=1,\,2,\,3;~{}l=1,\ldots,M_{i}),$

where $N_{i}^{(\eta)}$ are the numbers of polynomials approximating the $\eta^{\rm th}$ parameter along $i^{\rm th}$ coordinate axis, and derivatives ${{\partial x_{i}}/{\partial r_{i,l}}}$ are defined by equation (B9).

Appendix E The Chebyshev integration

Consider the 1D direct Chebyshev transform of a function $f(z)$ at the interval $z\in[0,\,1]$

[TABLE]

where ${\bm{C}}(f)$ is the symbol of the direct Chebyshev transform of function $f$ and the roots of the Chebyshev polynomials $\lambda_{j}$ are given by equation (A6).

Integration in the Chebyshev domain reduces to a simple multiplication by the integral matrix ${\bm{J}}$,

[TABLE]

where the matrix ${\bm{J}}$ isVolovodenko1981

[TABLE]

To evaluate the integral in equation (E2), we apply the inverse Chebyshev transform of the vector ${\bm{J}}\,{\bm{C}}(f).$ If integral at only one point $z$ is needed, the inverse transform becomes a single dot product

[TABLE]

where ${\bm{T}}(z)$ is the vector of the Chebyshev polynomials $T_{k}$ at $z$.

Integrals (8) and (11) are to be computed at $z=1$, therefore, we obtain from equation (E4)

[TABLE]

where

[TABLE]

and the vector ${\bm{T}}(1)$ is given by equation (B4b).