Travel times and ray paths for acoustic and elastic waves in generally anisotropic media

Wavefield travel time tomography is used for a variety of purposes in acoustics, geophysics and non-destructive testing. Since the problem is non-linear, assessing uncertainty in the results requires many forward evaluations. It is therefore important that the forward evaluation of travel times and ray paths is efficient, which is challenging in anisotropic media. Given a computed travel time field, ray tracing can be performed to obtain the fastest ray path from any point in the medium to the source of the travel time field. These rays can then be used to speed up gradient based inversion methods. We present a forward modeller for calculating travel time fields by localised estimation of wavefronts, and a novel approach to ray tracing through travel time fields. These methods have been tested in a complex anisotropic weld and give travel times comparable to those obtained using finite element modelling while being computationally cheaper.


Introduction
Ultrasonic non-destructive evaluation (NDE) is an umbrella term that describes a range of techniques used to check for flaws during the manufacture and maintenance of safety-critical infrastructure built from materials such as metals and concrete [1]. Usually this involves imaging the medium using algorithms that embed the assumption that the materials are homogeneous and isotropic [2]. However, for many materials of interest this is not the case. For example, in common polycrystalline welds, there exists a locally anisotropic grain structure which causes waves to scatter and refract. Using homogeneous assumptions in this case can dramatically reduce the reliability of the detection and characterisation of material defects. Given correct information about the spatially varying material properties, the heterogeneity and anisotropy can be compensated for to improve the accuracy of estimated defect properties [3].
In geophysical travel time tomography of the solid Earth's structure, orientations and properties of rocks and of the fluid reservoirs that they contain are imaged using travel times of seismic or electromagnetic waves measured between wave energy sources and receivers. Almost all rocks are significantly heterogeneous and anisotropic, and such effects must be accounted for and indeed imaged for purposes ranging from assessing earthquake hazards to monitoring subsurface dynamics [4]. In both NDE and geophysics, it is therefore necessary to model the effects of heterogeneous, anisotropic materials on waves.
Traveltime tomography is a procedure that allows us to use observed travel times (the fastest time for waves to travel between two known points on or in the object or medium) to estimate the spatially varying material properties within that object by comparing observed travel times with modelled time of flight data [5]. Traveltime tomography is often preferred over full waveform techniques because calculating travel times is computationally cheaper than modelling the full scattered wave field [6]. Due to the complexity of the travel time tomography problem, Monte Carlo sampling methods are often deployed to find the solution and assess its uncertainty, but these methods are computationally expensive for large datasets and high dimensional parametrised systems [7,8,9]. Variational Bayesian Inversion (VBI) approaches can offer a computationally efficient approach for these inverse problems by formulating them as an optimisation problem instead of a sampling problem [10,11,12]. Nevertheless, for many tomographic inversion techniques wave field data have to be modelled many times so computational efficiency is key.
We would like to image the interior of anisotropic materials, and in particular to invert travel time data collected on the surface of an object for grain orientation maps. In order to achieve this we need a forward model that calculates travel times of first arriving energy travelling between source and receiver locations. We also want corresponding ray paths along which the energy travelled. These ray paths can be used to increase the efficiency of inverse methods that require gradients of the travel times with respect to parameters that represent properties of the medium. Travel time fields can be used to find the fastest ray path between two given points in some medium. Fast Marching Methods (FMM) [13,14,15,16] and Fast Sweeping Methods [17,18,19,20,21] are two such efficient methods for calculating travel time fields which have similar performance [22]. In this work we employ a Fast Marching Method and while the FMM has been used to calculate travel time fields in isotropic media [13,15], an extension is required for anisotropic media. Typically, previous extensions of the FMM to anisotropic media have limited to weak, elliptical or tilted transversely isotropic media [23,24]. For generally anisotropic media, the Anisotropic Multi-stencil Fast Marching Method (AMSFMM) [14] can be used, where materials are defined using vectors describing the group and phase velocity curves, facilitating the study of any material with smooth velocity curves and 90 • rotational symetry. However the AMSFMM is unsuitable for ray tracing as its travel time errors are directionally biased. Note that although ray tracing can be performed in isotropic media using gradient descent through the travel time field (since the ray path is perpendicular to the wavefront), this is not the case for anisotropic media and so a different ray tracing method is required.
In this paper we present a novel forward modelling algorithm for calculating travel time fields, the Anisotropic Locally Interpolated Fast Marching Method (ALI-FMM), and a novel ray tracing method which uses ALI-FMM to calculate the fastest ray path in anisotropic media. Specifically, we study the arrival times of longitudinal waves in two dimensional orthotropic media. However, the framework retains the flexibility of the AMSFMM to facilitate the study of generally anisotropic media restricted only by smoothness and the presence of 180 • rotational symmetry. The success of the methodology is examined in the context of reconstructing a complex and strongly anisotropic velocity structure derived from a real weld. This represents an example typical of non-destructive testing applications. These methods are thus shown to enable travel time tomography to be applied in highly heterogeneous and anisotropic media.

Fast Marching Method in Anisotropic Media
A travel time field gives the shortest travel time from a fixed source point to all other points in the domain I = X × Y = {(x i , y j )|x i ∈ X, y j ∈ Y } where X and Y are the sets of x and y coordinates in the grid. For isotropic media FMM has often been used to calculate travel time fields, and works in a similar way to Dijkstra's algorithm to solve boundary problems. In this paper the FMM is used to solve the eikonal equation: where τ (x, y) is the shortest travel time from the source to the point (x, y). In isotropic media the phase velocity is independent of the direction of travel, in which case v(x, y, τ /| τ |) is the speed at (x, y). For anisotropic media v(x, y, τ /| τ |) is the phase velocity at location (x, y) in the direction which is normal to the wave front, which in turn is dependent on the source location. FMM works by giving each point one of three states: "known", "close" or "far" relative to the current wave front. When a grid point has a "known" state, the travel time is fixed. A "close" point, has a travel time assigned, however its value is not fixed and can be updated. A "far" point currently has no travel time. To propagate the wave front we take the point with lowest travel time from the set of points in the "close" state (times can be stored in a minimum heap for efficiency) and set it to a "known" point. We then update the travel times for all surrounding "far" and "close" points using a finite difference approximation. This will either involve updating the time for a "close" point or calculating a travel time for a "far" point and setting it to be a "close" point. This is repeated until all points are in the "known" state. This method differs from Dijkstra's algorithm when updating the times as an upwind finite difference approximation is used with multiple surrounding points with "known" or "close" states, while Dijkstra uses a single point and uses the travel time between two points. For locally anisotropic media Anisotropic Multi-Stencil Fast Marching Method (AMSFMM) [14] was developed to calculate travel time fields. This method uses a series of stencils on a 25 point finite difference scheme, where the direction of the stencil arms determines the velocity used for the finite difference approximation. Since those directions are constant across the whole grid, the travel times in those directions are very accurate, but travel time estimates in other directions can have large errors, especially in strongly anisotropic media. This bias towards particular directions makes this method unsuitable for ray tracing: when the geometry is rotated to lie on a different grid, travel times and hence ray paths may vary greatly. The errors can be reduced by using more stencils, but this is computationally more expensive and requires points to be used further from the grid point where the travel time is being estimated. This may in turn require a finer grid to be used in order to reduce errors at interfaces between different materials. AMSFMM also requires 90 • rotational symmetry in the group velocities, while many materials have only 180 • rotational symmetry which limits the types of materials and applications to which it can be applied. For these reasons we seek a different finite difference scheme to update travel times in FMM. The method proposed herein, uses a finite difference scheme similar to that of Eaton 1993 [25], where an estimate of the wave front is calculated by linear interpolation to obtain points in space with equal travel time. Connecting such points produces an estimate of the wave front locally, which can then be used to estimate the travel time at new points on the grid. Since the estimated wave front can be oriented in any direction, we employ the full anisotropic velocity curve rather than a limited set of fixed directions as in AMSFMM. We also use a square grid instead of the hexagonal one used in Eaton 1993 [25]. Use of the hexagonal grid means that the maximum angle between the estimated wavefront normal and the preferential wavefront direction is 30 • , whereas in our framework this value is reduced to 22.5 • (see Section 4.1). Additionally, in applying a hexagonal grid, either the top/bottom or left/right boundaries of the domain do not align with the grid. In this case the maximum angle between the estimated wavefront normal and the preferential wavefront direction can reach 60 • , while with implementation of our triangular stencils this remains at 22.5 • . AMSFMM and ALI-FMM uses group and phase velocity curves which can be obtained by solving the Christoffel equation from a material's stiffness tensor and density. There are at most 21 (independent) stiffness tensor elements which are often written as a 6 x 6 symmetric matrix, but for many materials such as transversely isotropic media there are only a few independent elements.

Calculating velocity curves
In anisotropic media, where the speed at which a wave travels is dependent on its direction of incidence, there are two different relevant and measurable velocities. The group velocity gives the speed at which a packet of energy comprising a band of neighbouring wave vectors travels in a certain direction, while the phase velocity gives the speed of wave fronts with individual wave vectors ( Figure 1). For travel time tomography a material can be defined by group and phase velocity curves. These velocities can be derived from the Christoffel equation given a stiffness tensor and density, and so the group and phase velocity curves contain information about those parameters of the medium [3].
Since the velocity in anisotropic media is dependant on the angle of incidence, we want to determine the longitudinal or so-called P-wave (longitudinal is used for the remainder of this paper) group and phase velocities in a medium with constant stiffness tensor c for different angles. Hooke's law [26] requires that, where σ ij is the stress tensor, ε ql is the elastic strain tensor and C ijql is the elastic stiffness tensor. Using the position vector x and displacement vector u = (u i ), Newton's second law can be rewritten in tensor notation as using the strain-displacement relationship and symmetry of the tensor notation we obtain Cauchy's law of motion, and thus we can look for a three dimensional plane wave solution of the form where A is a constant, p i are the components of the polarisation vector which describes the direction of the phase vector, k i = kn i where k is the wave number and n i is a component of the group vector direction n. The derivatives in equation (3) are then given by, Substituting these into equation (3) gives which rearranges into the Christoffel equation where δ ij is the Kronecker delta function M is the Christoffel matrix defined as M ij = C iqlj k 2 n q n l . This can be solved as an eigenvalue problem to determine the phase velocities of all wave polarisation modes (quasi-longitudinal, quasi-shear horizontal and 4 quasi-shear vertical) using equation (8). The wave front normal n j is given from the eigenvectors and the phase velocity used in equation (1) is given by where λ i are the eigenvalues of M. Since longitudinal waves are the fastest, the largest eigenvalue corresponds to the longitudinal wave. The group velocity is the speed in which a packet of energy goes in one direction and is given by We can multiply the Christoffel equation by p i to obtain which gives Alternatively since the group direction is the same as the wave front normal, we can determine the angle between the group and phase vectors (shown in Figure 1) and use that to determine the group velocity.

Solving the Christoffel equation at run time
In our forward model, materials can be defined by using a table of velocities at each integer angle; since changing the density of a material scales all velocities, the same table can be used to calculate the velocity for any density. However, if we want to use a variety of stiffness tensors, then we are limited to materials enumerated in the table of velocities. If we wish to use different tensors at each grid point this is potentially very computationally and memory intensive, especially if applied in 3D. This is because the Christoffel equation would have to be solved at every grid point, for every angle. Instead we can solve the Christoffel equation at run time: this is computationally cheaper because the equation is solved for a single angle rather than generating the full velocity curve. Phase velocities are calculated using the method described above. Calculating group velocities is more difficult as it is calculated using the group angle rather than the phase angle. To achieve this we determine the phase angle from the group angle; the phase angle is then used to solve the Christoffel equation. We limit ourselves to materials in which wave energy that travels from a source to any other point on a given 2D plane will not go out plane (into the third dimension). We use orthotropic materials which have three orthogonal planes of symmetry; these include transversely isotropic, cubic and isotropic materials (for isotropic materials there is no need to solve the Christoffel equation as materials are defined by a single velocity). For transversely isotropic materials we can use the axis of symmetry to calculate velocities in 3D media. The stiffness matrix for these materials takes the form In transversely isotropic media the stifness matrix is typically written with c 11 = c 22 , c 13 = c 23 and c 55 = c 66 . Since this is isotropic in the XY plane, we are obtaining velocity curves from the YZ plane, therefore our phase vector is   0 cos(θ) sin(θ)   for phase angle θ. We use this to determine the Christoffel matrix for a general phase angle which gives One of the eigenvalues is cos 2 (θ)c 66 + sin 2 (θ)c 55 , however the eigenvector is (1, 0, 0) which corresponds to a shear wave which vibrates out of plane. The other two eigenvalues are calculated by solving a quadratic equation. The solutions are (α + γ) ± (α − γ) 2 + 4β 2 2 where α = cos 2 (θ)c 22 + sin 2 (θ)c 44 , β = cos(θ)sin(θ)(c 23 + c 44 ) and γ = cos 2 (θ)c 44 + sin 2 (θ)c 33 . Taking the largest root (α + γ) + (α − γ) 2 + 4β 2 /2 gives the eigenvector for the longitudinal velocity since it is the fastest. We know that the eigenvector for the group velocity is for group angle φ. The properties of eigenvectors are used to determine the phase angle for a fixed group angle by solving (M − λI 3 )X = 0 for some eigenvalue λ ∈ R. There are two unknowns, θ and λ. Expanding the equation gives By using the double angle formula and rearranging equations 16 and 17 we obtain We divide equations 18 and 19 by cos(φ) and sin(φ) which can not be done for angles which are multiples of 90 • ; however since both velocity curves are smooth functions and are symmetrical about these angles (materials are orthotropic), the group and phase angles must be equal in these directions. Equations 18 and 19 are used to eliminate λ which gives Since the equation is in the form δ cos(2θ) + η sin(2θ) + κ, where δ, η, κ ∈ R, the substitution y = tan( 2θ 2 ), where θ = tan −1 (y) is used to obtain sin(2θ) = 2y 1+y 2 and cos(2θ) = 1−y 2 1+y 2 . This simplifies to a quadratic equation which we can easily solve , and two solutions are obtained. However, one is for a shear wave mode and since we want the longitudinal wave mode the solution which gives the largest eigenvalue is used. For φ ∈ [0, 180] the remainder is used since tan has multiple solutions when inverted (we seek angles between 0 • and 180 • ). For φ ∈ [0, 180] we get the solution If φ ∈ [0, 180] we can use the 180 • rotational symmetries to calculate the corresponding values, after which we obtain the phase velocity The group velocity can also be obtained using equation (14) where θ = θ − φ. For all of these equations, the only stiffness tensors required are c 22 , c 23 , c 33 and c 44 as the remaining stiffness tensors do not affect the velocities (when in plane). When these methods are implemented some of the numbers are very large and can cause overflow errors so we put stiffness tensors in MPa and scale velocities by 1000 to reduce the magnitude of numbers in the calculations. This scaling works for 64 bit numbers but does not for 32 bit numbers as the values are still too large. We also round angles when they are within 0.01 • of a multiple of 90 to bound these numbers which ensures the singularities of tangent function do not cause problems.

Anisotropic Locally Interpolated Fast Marching Method (ALI-FMM)
We now present a new algorithm based on the Fast Marching Method [13] for calculating the travel time from a source to all points on a grid. ALI-FMM works by approximating the wave front using linear interpolation [25] over a series of stencils which are used in an upwind finite difference approximation to propagate the wave front.

Stencils
We use a family of 32 finite difference stencils to approximate the wave front as locally planar, 16 of which form a square, and the other 16 form a triangle with a fourth grid point on the edge of the triangle. Three of the points in each stencil have travel time estimates and the other locates the point at which a travel time is to be estimated. We use the points with travel time estimates to calculate an estimate of the wave front. We begin by assigning labels to the stencil points. Tables 1 and 2 define points A, B, C, D for square and triangular stencils, respectively, and diagrams showing the geometries of each stencil are given in the appendix (Figures 23 and 24). In each case points A, B and C have travel time estimates and point D is The later condition implies that stencils occur in pairs that use the same points for A and D but with positions of B and C swapped. When these conditions are not satisfied the stencil is rejected. We find the point E on line AC that has the same travel time as B using linear interpolation of the travel times between A and C: We approximate the wavefront locally by the straight line joining B and E (Figure 2).

Travel time estimate from wave front
In order to calculate an estimate for the travel time at D we solve the eikonal equation (Equation 1). For anisotropic media, v(x, y, τ /| τ |) is the phase velocity at (x, y) in the direction which is locally normal to Table 2: Labels of the points used in each of the triangular stencils for the finite difference approximation using grid coordinates (i, j). Stencils are shown graphically in figure 24. the wave front. Since we have a locally linear wave front estimate we can use the orientation of its normal vector as an estimate to the direction of τ (x, y), | τ (x, y)| is obtained from the phase velocity curve using that direction. Combining these we get an estimate of τ (x, y) which is used to determine a travel time estimate. We start by finding the point F , which is the closest point to D on line EB (see Figure 2). The location of F is obtained by finding the line perpendicular to EB passing through D, and then finding the point of intersection between this line and the wavefront estimate. The distance and direction from F to D is then obtained. Line EB can be written as for λ 1 ∈ R. The direction of line F D can be found by rotating The normal to the wave front is the phase angle/vector. We therefore use the phase velocity and distance to determine the difference in travel times (τ D − τ F ) between F and D. Since line EB is our approximation of the wave front we assume where v(x F , y F , θ ) is the phase velocity of the material at F at angle θ . For this we use the material orientation at D and the direction of the normal to the wave front to calculate the effective angle θ = (θ D − θ n + 90)mod180, where θ D is the material orientation at D, θ n is the direction of line EB (θ n ± 90 is the direction of the normal vector, perpendicular to the wave front). Since we assume the velocity has 180 • symmetry, we take the remainder. These velocities can either be obtained using linear interpolation with a table of velocities or by solving the Christoffel equation at run time:

Choosing which stencil to use
Square stencils are often more accurate than the triangular stencils as there is less bias toward some directions due to the symmetry of the stencil. Triangular stencils are only used when there are insufficient points with travel time estimates to use a square stencil. This should only occur when there is high curvature on the wave front. Triangular stencils are also used when the grid point is at the edge of the grid (some of stencils 5 to 8 can then not be used as they require a point which is not on the grid). For these stencils we want the group vector from the estimated wave front to D to intersect line EB as this is the region where we have performed a wave front estimate. This is more likely to occur the closer point F is to the middle of line EB. For square stencils, we want to use the stencil which has a wave front which is closest to line CB. This will also have lower error in the linear interpolation, and brings F closer to D making the estimation of τ D − τ F smaller and subsequently the error. Since we get closer to line CB the closer the travel times at C and B are, the stencil with the smallest difference in travel time between C and B is used. For triangular stencils we want the perpendicular bisector of points B and E to pass close to point D, so that F is close to the middle of line EB. This can only happen for the second stencil in each pair and points EDB form an isosceles triangle. This happens when the length of line ED is √ 2 (using grid indexing), the length of line AE is 2 − √ 2 and the line EC is √ 2 − 1. Therefore a in equation (26) would equal 2 − √ 2. We therefore find the stencil where τ B is closest to ( √ 2 − 1)τ A + (2 − √ 2)τ C for the second stencil in a pair (same calculation for the first stencil in each pair except that the positions of B and C are swapped). A stencil is only used when τ D > τ A , so no issues arise when the wave front is travelling in the opposite direction. There are two sizes of the square stencils (stencils 1 to 4 and 5 to 8), and where both are feasible the smaller stencils are more likely to be chosen because the difference in travel time between C and B is more likely to be smaller. These stencils use points closer together (the small square has side lengths of 1 and the larger one has length √ 2), meaning that the small stencils are often more accurate. The same applies when comparing stencils 1 to 4 and 9 to 16 when the grid point is on the edge of the grid. In rare cases in which no stencils can be used, the finite difference approximation from AMSFMM [14] is used. A flowchart for choosing which stencil to use is given in Figure 3.

Propagation around source
When wave propagation is initiated around a localised source, we commonly have too few grid points to represent the curvature of the wave fronts accurately. If the orientation of the medium's anisotropy is constant in space in the vicinity of the source we can use straight rays to calculate travel times out to a certain distance before changing to the finite difference method. Unfortunately we can not use this method for our applications in highly heterogeneous media. We therefore use a finer grid around the source (red and green grids in Figure 4) and assign orientations of the medium on the fine grid to be the orientation of the closest point on the original grid. The reason we do not use interpolation of the materials is because it would make our ray tracing method more computationally expensive since travel times of straight rays can not be calculated on the coarser grid. This would also cause issues around the source as it would not be homogeneous. We also use a very fine grid around the source (green grid in Fig 4) so that there is no change in the orientation map sufficiently close to the source (points with x and y coordinates from 3.5 to 4.5 in Figure 4). This allows us to use straight rays to calculate travel times to points on the latter grid. These points can then be used to begin wave front propagation using the finite difference method. Since Here we use a relative grid density of 9 around the source for a distance of 1(green region), before moving to a relative grid density of 3 at a distance of 3 (red region), and finally onto the original grid (black).
using the finer grid is much more computationally expensive than the original grid, once the wave front reaches the boundary of a defined square region centred around the source, we move to a coarser grid and repeat this process until we are on the original grid (see Figure 4). To create finer grids we increase the number of grid points in each direction by an odd number. An integer value means that we need not use interpolation to determine travel times on the coarser grid, and an odd number means that when assigning material properties to grid points there is a unique closest point on the original grid. We have used a factor of 3 for each grid refinement starting with a relative grid density of 27, and changing to a coarser grid when we calculate a travel time at a point which is at a distance of more than 1, then 6, and finally 13 grid points from the source (on the original grid) either horizontally or vertically. By using a finer grid and/or propagating further from the source before changing grids we can obtain more accurate travel times, but with an increased computational cost. This is because the change in wave front angle between grid points is smaller, so the approximation of the flat wave fronts (see Figure 2) is more accurate. For some applications we may want to keep the finer grid around the source instead of moving all points onto the coarser grid if there is a large difference in velocities spatially. This is because the wave front might not have moved far enough from the source in some directions, which would give inaccurate travel time estimates on the coarser grid due to the high curvature. In FMM all points neighbouring a "known" point are either "known" or "close". This ensures that the wave front continues to propagate. When moving to a coarser grid, "known" and "far" on the finer grid can be next to each other on the coarser grid which does not allow the wave front to propagate. Points which are "known" but have a neighbouring point which is "far" are added to the minimum heap ("close" points are stored in the minimum heap, and when they are removed the surrounding "close" and "far" points are updated). This allows the neighbouring points to be changed to "close" points and obtain travel time estimates, while keeping the travel time at "known" points fixed. Source factorisation methods [21] are an alternative approach for modelling around the source, however using grid refinement allows us to reduce the errors around the source by using finer grids or using them over a larger area (see Table 3).

Results
We tested ALI-FMM using velocities obtained by solving the Christoffel equation for an austenitic steel, a cubic material with stiffness tensor elements c 11 = 203.6GP a, c 12 = 133.5GP a and c 44 = 129.8GP a and density of σ = 7850kg/m 3 [27]. The stiffness matrix for a cubic material takes the form In this example we know the true traveltimes for homogeneous media (anisotropic orientations at every point are equal) which can be calculated using a straight ray from the source to any point on the grid, we can calculate the errors when using AMSFMM [14] and the new method herein. ALI-FMM commits the largest errors where the velocity has higher angular curvature. This is expected since approximating the wavefront locally by a straight line is less accurate for higher curvatures. In AMSFMM we have little to no error in the directions of the stencil arms, however there is error between these directions. The errors are also very dependant upon the orientation of the medium relative to the grid (see Figures 6 and 7). In Figure 7 the curvature of the wavefront is high between the directions of some of the stencil arms in AMSFMM which increases the error, in such orientations the ALI-FMM performs much better in these cases and has very similar error magnitudes to those in Figure  6. Figure 8 shows ALI-FMM has a lower and much more consistent average error, and the maximum error is approximately a factor of 4 smaller than AMSFMM for different orientations of this geometry. Figure 9 shows that the median, upper and lower quartiles are lower in the new method except for a few angles of which we encounter similar errors in the lower quartile. One way to improve the accuracy is to use a finer grid by adding additional grid points, so that the finite difference method is performed on a smaller area. When creating a finer grid we increase the number of grid points in each direction by a factor which we call the subgrid size which is an odd integer. We can assess the error for different sizes by using a homogeneous medium in which we can obtain the true travel times using straight rays. In order to have a fair comparison we only compare the travel times at points on the original grid. The original grid is 21 by 21 with the source at the center. The finer grids around the source remain the same as the previous results for a subgrid size of 1, however for the other sizes we use a relative grid density of 9 for a distance of 2 grid points, and a relative grid density of 3 for a distance of 5 grid points      (distances are measured on the original grid). The grid used around the source is different to the original grid since we have already created a finer grid, hence there is a region around the source which has the same material properties as the source. An anisotropic orientation of 0 • was used as this had the lowest errors in the AMSFMM results. Table 3 shows that the error in ALI-FMM is similar when using a subgrid size of 1 and 3. This is likely caused by the differences in the finer grids around the source. The errors for AMSFMM decrease quickly on the finer grid, however they start to increase after a subgrid size of 15. In ALI-FMM, these errors continue to decrease for finer grids and should indeed continue to decrease for all finer grids. The errors when we solve the Christoffel equation during run time are lower than those obtained using a table of velocities despite using the same methods. This is because the velocity calculations are more accurate as there is no linear interpolation between integer angles. Another way to increase accuracy is by using a finer grid over a larger area for the source which will have a similar affect to the subgrids, but only around the source.
All of the different codes use the same implementation of the Fast Marching Method for the binary tree and for storing the states ("far", "close" and "known") of each point. These codes only differ in the finite difference method used and the finer grid around the source. This ensures that we can make a fair comparison of the compute times for different sizes of grids. A square grid with the source at the center has been used for benchmarking. The methods have a complexity of O(nm log(nm)) in Big O notation where n and m are the number of grid points in each direction. This complexity is because there are nm points and we use a minimum heap structure to keep track of all points in the "close" state, which has logarithmic complexity for binary tree operations. However the operations in the heap are much faster than the finite difference methods making the complexity close to linear. Figure 11 shows that AMSFMM is faster for small grid sizes, which is due to the finer grid around the source. Despite our finite difference method being slower than those used in AMSFMM, ALI-FMM is more efficient for larger grids. This is because when using AMSFMM, the stencil is chosen after the finite difference method is applied to all four stencils. ALI-FMM chooses the stencil before the finite difference method is applied to the stencil. As expected when we solve the Christoffel equation during run time, the time required is higher, but the added flexibility of using stiffness tensors at all grid points makes this method applicable to more complex problems.

Ray Tracing in Anisotropic Media
Finding the fastest ray path in isotropic media can be done by gradient descent through the travel time field from a receiver back to the source. This works because the ray paths are always locally perpendicular to the wave front in isotropic media. In anisotropic media the wave front is rarely perpendicular to the ray path so gradient descent will not produce accurate ray paths. Since the gradient gives the phase direction, the corresponding group direction is the direction of the ray path (see Figure 1), which can be obtained by solving the Christoffel equation (see Section 3). This would requires us to know the relationship between group and phase directions which is unknown when using a table of velocities. In this section we show how travel time fields calculated above can nevertheless be used to find rays.

Finding points along fastest ray paths
In order to find the fastest ray paths in anisotropic media, the grid is divided into sections by a series of planes. The ray path will be able to cross a plane once, and so the choice of these planes is important. The planes are chosen such that the points of intersection with the ray path have similar spacings. Since we use these points to parametrise the ray, we use equispaced horizontal, vertical or diagonal intersection planes passing through grid points, shown in Figure 12.
We use four sets of planes because the choice of planes can affect the accuracy of the path. For example in Figure 12, the horizontal planes are likely to perform poorly because we can only fit 2 planes between the source and receiver, providing a ray path approximated by only four points (source, receiver and two intersection points). In addition, the ray path can not travel above the source or below the receiver as all points must lie on a plane; in heterogeneous media we expect rays to deviate from such light constraints.
The spacing between planes is chosen so that the distance between intersection points is small enough to obtain sufficient detail about the ray path and large enough that the range of angles for lines BC is sufficiently finely sampled on the boundary to get an accurate position for the next point. By putting the geometry onto a finer grid we can get a more accurate ray path which allows for a more detailed path. By using a finer grid we can use all possible planes on the original grid while having enough points for an accurate ray path (see Figure 13). The ray path is calculated using Fermat's principle, starting at the source and adding one intersection point at a time. Assuming that the first n points on the ray path have been Figure 12: Planes (green) used for ray path calculation between source (red star) and receiver (blue triangle) Figure 13: The intersection between all thick black/blue lines are points from the original grid while all intersections of black and blue lines are points in the finer grid (In this case a 3x3 subgrid is used, but a finer grid is often used). The blue line is the plane where the next point on the ray path must lie. The red lines are sections of the ray path, where solid lines between A and B are fixed and the dashed line between B and C is a possible section of ray paths. The green lines are contours of the travel time field with source at D. calculated, the next point is determined by finding the fastest path from the last calculated point on the ray path to the receiver. Let A be the source, B be the last calculated point on the ray path, C be some grid point on the plane where the next point on the path may lie and D is the receiver. Since the next point lies on a predetermined plane, we find the smallest travel time from B to D through C shown in Figure 13. By finding the minimum total travel time of all points along the plane, we find the next point on the ray path. This is true because all other points on the plane result in a higher travel time and so by Fermat's principle the fastest ray path does not pass through these points.
The travel time from B to C is calculated assuming a straight ray, since our path will not have any points between them. The time can be obtained from the group velocity curve and the orientation of the closest point on the orientation map. For this step we can use the original grid as this is more efficient to calculate and since the geometry is identical, the travel time is unchanged. The fastest path from C to D is taken from the travel time field with D as the source which can be obtained using ALI-FMM or any other method. This is the smallest travel time from D to C, which is equal to the smallest travel time from C to D, since the travel time is identical when source and receiver are swapped. The smallest travel time from B to D through C under our constraints is the sum of these two travel times. Since the travel time field was only calculated at grid points, the minimum value on a plane is estimated by fitting a quadratic function to each local minimum and the two neighbouring points. We then find the minimum over these quadratic curves to estimate the global minimum. This allows the ray to travel to any point on the plane, instead of only to grid points giving a more accurate and smoother ray path. When we have multiple local minima there may be an alternate ray path, however this path will be slower or at best equally fast.
The choice of plane can make a difference to the accuracy of the ray path. However, we can calculate a ray path for each of the plane directions and integrate along the path to find the fastest time. The best results are often obtained when the planes are nearly perpendicular to the ray path. For more complex examples and to increase efficiency we can apply the same methods, but switch between plane directions based on the orientation of the last section of the ray. This allows us to obtain an accurate ray path when the ray intersects some planes more than once for all plane directions. This also reduces the computational cost as only one ray path is calculated. The plane direction is chosen by finding the plane for which the angle between the direction of the last section of the ray path and that plane is largest. For the initial point after the source, the direction from the source to the receiver is used. The position of the plane is calculated using the same spacing for the fixed plane directions, but using the distance between the plane and the last calculated point rounded to the nearest integer, to ensure that the plane passes through grid points. We stop when the distance between the last calculated point and the receiver is less than 1.6 grid points (on the Figure 14: Ray path in isotropic media with a velocity gradient showing plane directions used for ray tracing. Source is on the left and receiver is on the right. original grid). This number is chosen arbitrarily, so that we do not pass the receiver and the spacing between the last two grid points is not too large. Since the angle between the direction of the ray and the plane is large, we can also limit the number of intersection points in the plane which are tested at each step. This does limit how much the ray can change direction, however with a sufficient number of points this should not effect the ray path. Figure 14 shows where each of the different plane directions are used in an example with a velocity gradient. The first point is calculated using a vertical plane (black); we then use three different plane directions as the direction of the ray changes, and finally add the last point in the ray when sufficiently close to the receiver (blue).

Results
All results from herein are obtained using grid spacings of 0.001 metres (1 mm) and using the same stiffness tensor as for the earlier results (c 11 = 203.6GP a, c 12 = 133.5GP a and c 44 = 129.8GP a and density of σ = 7850kg/m 3 ). Ray tracing was first applied using travel time fields from AMSFMM [14], but the error in travel-time fields is heavily biased towards low times in certain directions which cause issues in ray tracing. This often causes the ray paths to travel close to these directions and the four different planes give different ray paths due to this bias. These issues should be present for any ray tracing approach using these travel time fields unless the directional bias is compensated for.
In the homogeneous medium in Figure 15 the true solution is a straight ray from source to receiver (the anisotropic orientation is the same at all grid points). For AMSFMM the ray using the horizontal planes has the worst fit, which is to be expected as the angle between the direction of travel and the planes are much closer than for the other planes. This also means that there are fewer points in the ray path. Even though these ray paths are not straight, the travel times along these paths are lower than AMSFMM, except for the ray obtained using horizontal boundaries. Travel times of ray paths are generally lower than those obtained directly from AMSFMM travel time fields despite being calculated from AMSFMM travel time fields, with exceptions when the planes do not fit the path well or the true ray path is aligned along stencil directions (Table 4). This reduction in error is likely due to the ray tracing using the full velocity curve to calculate the travel time along line BC, which makes it more likely to travel in directions between stencil arms where the travel time field predicts the wavefront is slower than the true solution. This means that the ray path has less bias towards the directions of stencil arms than the travel time fields. The ray paths obtained using ALI-FMM are much straighter and all give similar paths for each of the plane directions as the travel times fields have less bias than AMSFMM.
The ray tracing method has also been applied with a circular anomaly using the same stiffness tensor we have used previously but with anisotropic orientations of 60 • inside and 15 • outside of the circle. In Figure  16 the source and receiver have the same x coordinate so there are no vertical planes between them. This means that the ray paths obtained using vertical planes consists of only two points which are the source   and receiver. For the other plane directions, the rays obtained using AMSFMM are biased towards certain directions and give different paths for the different plane directions due to the errors in the travel time fields similar to Figure 15a, while the rays using ALI-FMM have far less bias and the different plane directions give similar paths and are much straighter. We then constructed a ray path by changing boundary orientation along the ray as described above (Figure 17). This method was only tested using ALI-FMM as ray tracing using AMSFMM is heavily biased and gives far less accurate ray paths due to the error in the travel time fields. The corresponding ray path in Figure 17 is similar to those obtained using a fixed plane direction, but since the plane direction may change along the ray, the problem that occurs when using a vertical boundary is unlikely to occur.  Integrating travel times along the ray path will always give a time that is greater than or equal to the true travel time obtained on the discretized grid. The exact optimal path can be computed exactly. Assume that the disk where the material coefficients are modified has radius one, so that its boundary can be parameterized as (cos θ,sin θ), θ ∈ [0, 2π]. Then the optimal path is the concatenation of three straight segments, from (x 0 , y 0 ) the source point, to (cos θ 1 ,sin θ 1 ) a first point on the modified region boundary, to (cos θ 2 ,sin θ 2 ) a second point on the same boundary, to the target point (x 1 , y 1 ). The travel times along each of these segments can be computed using the material properties and the expressions of section 3. Then optimizing over θ 1 , θ 2 ∈ [0, 2π] one can obtain the optimal path. This path (referred to as the 'True Path') is plotted in Figure 17 and produces a travel time of 1.828731e-05s. Note this 'true' travel time is in fact greater than that calculated using our proposed methods. This is because the discretisation of the underlying geometry (specifically the staircasing of the circle boundary) permits the ray to travel in preferential directions for longer. The optimal solution on this discretised geometry (obtained using a brute force approach) produces the shortest traveltime across all methods (1.826223e-05s to 7 significant figures). Since the rays obtained in Figure 16 using horizontal planes gives a straight ray, this is by far the slowest ray path and shows that even though we use the same material throughout, the fastest path is not always a straight line from source to receiver since we have an anisotropic medium with varying anisotropic orientation. Table 4 also shows that a large improvement in the travel time estimates can be obtained when using ray tracing with AMSFMM in comparison to times obtained directly from the travel time field. ALI-FMM also manages to obtain lower travel times (less error) for all estimates of the travel time compared to AMSFMM, except for the vertical planes which have already been discussed.
In anisotropic media there are only a few examples where the true solution is known other than for homogeneous media. One way to construct an example is to deform a homogeneous medium using a change of coordinates [21]. In isotropic media there are some examples where the true solution is known. One of these is when there is a constant velocity gradient in which case the solution is a section of a circle [28]. When the gradient is in the form (g, 0), with source (x 0 , y 0 ) and receiver (x 1 , y 1 ), the circle has equation , v 0 is the velocity at the source and θ 0 is the initial ray direction relative to the velocity gradient. Unfortunately the value of θ 0 is unknown, however the location of the receiver is known and can be used to find the location of the circle. Equation (31) is used to solve for θ 0 : From the equation of the circle, the center is at (x , y ) = − v0 g , v0cotθ0 g , therefore (x, y) = − v0 g + x 0 , v0cotθ0 g + y 0 and the radius of the circle is v0 gsinθ0 . A ray path has been calculated for a velocity gradient with velocity 3000m/s at x = 0 and 7200m/s at x = 200, with a source at (x 0 , y 0 ) = (1, 30) and receiver at (x 1 , y 1 ) = (199, 180) (using grid coordinates, the length of a cell is used for scaling after the ray is calculated), which when using grid coordinates gives the parameters g = 21 and v 0 = 3021. From these parameters the true solution is calculated as a section of a circle with center (-142.857, 425.571) and radius of 420.918. Figure 18 shows that the calculated ray path is very close to the true solution. The maximum distance between these lines is approximately 0.3 grid points away from the true solution and true travel times are calculated by integrating along the ray path (using cell width of 0.001) giving travel times using the discretised grid of 5.0884409e-05 for the estimated ray and 5.0884053e-05 for the true solution, giving an error of 0.0007%. The error obtained is dependent on the error in the travel time field and the geometry used. For anisotropic media the error is likely to be higher. From Table 4, the error for the ray tracing using different plane directions is 0.096% (3 decimal places), on the discretised grid. However if the errors in the travel time field were lower, the errors in ray tracing are expected to be lower.

Ray tracing through an anisotropic weld
Ray tracing was performed using a model based on an anisotropic weld obtained using electron back-scatter diffraction [29] (Figure 19). The model was simplified using k-means clustering to remove some of the high spatial frequency heterogeneity to create a more reasonable boundary between the weld and parent material. Sources/receivers were also changed so as to act more similarly to a point source and were put on the top and bottom so ray tracing could be applied to the model (see Figures 19 and 20). The crystal orientations were also changed from a 3D orienation to 2D orientation since we require 2D velocity curves for this version of our codes. Travel times were estimated using A-scans obtained using finite element (FE) software Onscale with the driving function in time shown in Figure ??. While the driving function is not realistic, it helps to measure more accurate travel time estimates from the modelled waveform due to the rapid onset of energy.
Travel times were estimated by finding the first minimum/maximum above an arbitrarily chosen threshold. A second threshold is then calculated to be the pressure at the min/max divided by 100, and the travel time is estimated as the first time point above this threshold. Due to the method used for estimating these travel times there will be some error in the travel times obtained.
Comparing travel times calculated along ray paths and those from FE in Figure 22, we see that differences are relatively small, however for source 19 the differences are larger. The travel times along many ray paths are lower than those from the FE software, and since ray paths always give travel times that are higher than the true solution on a discretised grid, this is a good result. Figure 19: Model of a weld (from [29]). Background colours represent the material at each point (see legend) with materials pmt3 being transducers, medm and back15 are part of the transducer array, watr is water, stst is an isotropic steel and all other materials are different anisotropic orientations of the same austenitic steel which make up the weld structure.

Discussion
To improve accuracy within the travel time fields a higher order finite difference approximation could be used by adding an additional point in each stencil to get three points with equal travel times for the estimated wavefront. These three points allow the wavefront to be estimated including its curvature. This would require a different method for the finite difference approximation, and there may be numerous challenges using such methods. If there is any bias towards higher curvature in the wavefront, the curvature in the modelled wavefront may increase and cause unstable and inaccurate results. Another way to reduce errors is to change the density of grid points based on the curvature in the wavefront since the highest errors occur in regions of high curvature. This would have a similar effect as the introduction of finer grids around the source which are used due to high curvature in that area. Since we would only increase the grid density when there is high curvature, this would be computationally more efficient than increasing the grid density for the whole grid. Note that shear wave (S-wave) arrival times could also be computed using the ALI-FMM finite difference calculations by substituting in the shear-vertical or shear-horizontal velocity curves and considering the other eigenvalues/solutions in Section 3. However, since this method was developed with traveltime tomography applications in mind, this study restricted examination to longitudinal waves only since these are the fastest and hence first arriving, meaning they are easier to extract from the collected data. ALI-FMM can be used in 3D wave propagation using the same techniques used in this paper, however the wavefront must be approximated using a 2D plane. This would require a new set of stencils to approximate the wavefront. Since three points with equal travel times are then required to estimate the wavefront, an additional point is required in the stencil compared to the 2D stencils used in this paper (hence a total of 5 points). The wavefront can also be estimated using curvature with an additional point using the same method proposed above for the 2D case. Ray tracing can be used in 3D by finding where the ray path crosses a set of 2D planes. Due to the additional plane directions that would be used and the number of points within each plane, ray tracing without adaptively changing plane direction would be computationally expensive.
The model used for ray tracing in this paper is input by giving material properties at all points in a grid. When the geometry is put onto a finer grid the detail in the geometry remains unchanged as grid points are assigned the material properties of the closest point in the original grid. If a different method is used that gives a smoother interface between regions, the accuracy of the ray path should improve. For example, in Figure 17 the geometry in the finer grid can be chosen using the equation of the circle. Modelling the geometry using Voronoi tessellation may be especially useful in NDT applications with polycrystalline materials [30]. Additionally once the ray path has been calculated we can remove points when we know that the ray must travel in a straight line. This would improve accuracy since the fastest ray path between two points in a homogeneous medium is a straight line. Another way to do this would be to look at where the ray path changes material, instead of using an intersection plane during the ray tracing which would be computationally cheaper. For example in Figure 17, the ray path should only change direction at the interface between the outside and inside of the circular region, so the ray path could be reduced to four points. These various improvements are left for further work, should they be required.