Summary
Wave front construction (WFC) methods are a useful tool for tracking wave fronts and are a natural extension to standard ray shooting methods. Here we describe and implement a simple WFC method that is used to interpolate wavefield properties throughout a 3-D heterogeneous medium. Our approach differs from previous 3-D WFC procedures primarily in the use of a ray interpolation scheme, based on approximating the wave front as a ‘locally spherical’ surface and a ‘first arrival mode’, which reduces computation times, where only first arrivals are required. Both of these features have previously been included in 2-D WFC algorithms; however, until now they have not been extended to 3-D systems. The wave front interpolation scheme allows for rays to be traced from a nearly arbitrary distribution of take-off angles, and the calculation of derivatives with respect to take-off angles is not required for wave front interpolation. However, in regions of steep velocity gradient, the locally spherical approximation is not valid, and it is necessary to backpropagate rays to a sufficiently homogenous region before interpolation of the new ray. Our WFC technique is illustrated using a realistic velocity model, based on a North Sea oil reservoir. We examine wavefield quantities such as traveltimes, ray angles, source take-off angles and geometrical spreading factors, all of which are interpolated on to a regular grid. We compare geometrical spreading factors calculated using two methods: using the ray Jacobian and by taking the ratio of a triangular area of wave front to the corresponding solid angle at the source. The results show that care must be taken when using ray Jacobians to calculate geometrical spreading factors, as the poles of the source coordinate system produce unreliable values, which can be spread over a large area, as only a few initial rays are traced in WFC. We also show that the use of the first arrival mode can reduce computation time by ∼65 per cent, with the accuracy of the interpolated traveltimes, ray angles and source take-off angles largely unchanged. However, the first arrival mode does lead to inaccuracies in interpolated angles near caustic surfaces, as well as small variations in geometrical spreading factors for ray tubes that have passed through caustic surfaces.
1 Introduction
Ray theoretical methods have long been used for the modelling of seismic wave propagation and are an important element in many seismic imaging procedures. In particular, dynamic ray tracing has proven to be a useful tool in cases where velocity models display lateral heterogeneity. However, standard implementations of dynamic ray tracing have the lateral resolution along the wave front limited by the initial choice of ray directions (Cerveny et al. 1977). This can lead to aliasing of the wave front in the forward problem and uneven coverage when wavefield properties are interpolated for the construction of traveltime tables used in seismic imaging, for example.
The problem of constructing traveltime tables that contain sufficient coverage, amounts to a solution of the two-point ray tracing problem. That is, given a background model, how does one derive the ray(s) or quantities integrated along the ray(s), which connect two points. A solution to the two-point ray tracing problem is an important requirement for forward modelling techniques such as scattering integrals, imaging procedures like seismic migration or tomography, as well as locating seismic sources.
As dynamically traced rays can behave unpredictably even in isotropic velocity models, past efforts to solve the two-point ray tracing problem have turned to alternative methods for solving the eikonal equation. Ray bending algorithms (RB, e.g. Jullian & Gubbins 1977; Um & Thurber 1987) iteratively perturb an initial path estimate to produce a stable traveltime. The final result from RB is dependent on the choice of initial path, and there is no guarantee the technique will converge to a particular arrival path. In contrast shortest path ray tracing (SPRT, e.g. Mosser 1991; Fischer & Lees 1993) is stable and should always converge to the first arrival. SPRTs compute all potential paths in a network of nodes, a search routine is then used to isolate the path with the minimum traveltime between two points. The resolution of SPRTs is limited by the network density, so there is a trade-off between accuracy and computation time which may become significant in 3-D media.
An alternative approach for finding wavefield quantities throughout a medium stems from finite difference and fast-marching finite difference solutions to the eikonal equation (FMES, Vidale Vidale 1988, 1990; Sethian & Poppovici 1999). These methods produce a map of first arrival traveltime throughout the medium. Given a traveltime field, rays can then be traced from points of interest backwards along the traveltime gradient to the source (Rawlinson & Sambridge 2004; de Kool et al. 2006). However, errors in the traveltime gradient near the source can produce substantial errors in ray angles in this region; information that is required in many standard processing applications. Alternatively, geometrical spreading factors and source take-off angles can be calculated without knowledge of the ray paths by initiating the eikonal solver from closely spaced sources (Vidale & Houston 1990). However, the requirement for multiple simulations can potentially lead to a significant increase in computation time. For example, up to 36 traveltime calculations are required to determine amplitudes at a receiver in a 3-D medium. As both SPRTs and FMES deal with first arrival traveltime fields, they cannot be used to model retrograde traveltime branches, although layered approaches can be used to track late stage arrivals such as reflections and head waves (Rawlinson & Sambridge 2004; de Kool et al. 2006).
Wave front Construction (WFC) provides an alternative solution to the two point ray tracing problem, and is the natural extension to standard dynamic ray tracing. Instead of considering rays in isolation, WFC tracks a network of rays as it passes through the medium. Wave front aliasing is prevented by controlling the spatial resolution of the ray network and interpolating new rays where necessary. The idea is not new, one of the first (2-D) implementations was Vinje et al. (1993), and various authors have expanded on the idea with different ray interpolation methods (Ettrich & Gajewski 1996; Lambare et al. 1996; Hildebrand et al. 1999), and extended solutions for 3-D media (Lucio et al. 1996; Vinje et al. 1996, 1999; Gibson et al. 2005; Lee & Gibson 2007).
Our purpose here is to describe a simple, but effective, WFC procedure, which can be used to solve the two-point ray tracing problem in isotropic 3-D media. We introduce two new aspects to 3-D wave front construction: a simple ray interpolation scheme based on approximating the wave front as a ‘locally spherical’ surface and a ‘first arrival mode’, which reduces computation times, where only first arrivals are required. Both the locally spherical approximation and the tracking of first arrivals have previously been implemented in 2-D WFC algorithms (see Ettrich & Gajewski 1996; Vinje et al. 1993, respectively,); however, until now, they have not been extended to 3-D systems. Recent approaches to 3-D WFC require the calculation of derivatives with respect to ray take-off angles at the source (e.g. Vinje et al. 1996), or constrain the distribution of source angles to corners of a quadrilateral cells on a cubed sphere mesh (Lee & Gibson 2007). In our method, a near arbitrary distribution of initial rays can be used, and the calculation of geometrical spreading derivatives along the rays is only required if they are to be interpolated at receivers (i.e. if ray Jacobians are used to calculate amplitudes). We show that our WFC approach can be used to rapidly calculate quantities such as traveltime, ray angles, source angles and amplitudes at points throughout the medium.
2 The Wave Front Construction Method
The basic procedure behind our implementation of WFC follows that of several previous studies (Lucio et al. 1996; Vinje et al. 1996). So here we limit ourselves to a brief discussion of the WFC method, with emphasis on aspects which differ from previous work, such as the ray interpolation method, tracking the first arrival wave front and amplitude estimation.
WFC begins with a set of rays with initial conditions and a triangular network which links adjacent rays on the wave front. Typically, we start with rays which radiate from a point source, but this is not a requirement; in fact, almost any initial distribution of rays could be used, as long as the ray directions are not coplanar. Once the ray field has been initialized, WFC consists of a loop thatthat runs for as long as rays exist in the medium. Each iteration of the loop advances the wavefield by a time-step, Δt, and monitors the ray network. The basic tasks within the propagation loop are as follows:
- 1
Each ray is propagated by a single time-step. We use a fourth-order Runga–Kutta method to time-step differential equations for ray position, slowness and geometrical spreading derivatives (i.e. dynamic ray tracing, Cerveny et al. 1977). Both the current and future ray points are stored for the duration of the loop.
- 2
Rays which exit the model are deleted from the network. WFC is finished once the last three rays are deleted.
- 3
The ray network is monitored for folds in the wave front (caustics). If the WFC is running in first arrival mode then the secondary arriving wave fronts are removed, and the resulting gaps in the wave front are closed by making new links in the ray network. Alternatively, for rays which have passed through caustics, the index of the ray trajectory (or KMAH index, Cerveny 2001) is incremented.
- 4
The distance between all pairs of adjacent rays is calculated. Where this exceeds a pre-defined maximum, the distance between points on the wave front is reduced by a diagonal swap with the opposing points of the surrounding quadrilateral. If a diagonal swap is not possible or it will not reduce the interray distance beneath the threshold, a new ray is interpolated.
- 5
The desired wavefield quantities are interpolated at gridpoints, which lie between current and future wave fronts.
- 6
The wavefield is advanced to the next step, and we return to step 1 for the next propagation step.
The only adjustable parameters are the ray time-step, the maximum distance between rays and the initial ray distribution, making the procedure very straightforward to implement.
2.1 The Interpolation of New Rays
A critical component of WFC is the interpolation of new rays where the existing wave front has insufficient spatial resolution. Previous studies have used a variety of techniques for the interpolation of new rays. Vinje et al. (1993, 1996) use a polynomial interpolation, whereas Hildebrand et al. (1999) used radii of wave front curvature to approximate the wave front. Most recently, Gibson et al. (2005) have used a second order Taylor expansion to interpolate the traveltime at an initial new ray position behind the wave front. The difference between the interpolated traveltime and the wave front time was then used step the new ray to it's correct position on the wave front.
These methods require the use of derivatives of the ray coordinates with respect to source take-off angles to evaluate the curvature of the wave front. However, special care needs to be taken when calculating ray derivatives for rays with take-off angles near the poles of the coordinate system. In standard ray tracing, where the range of take-off angles does not cover an entire hemisphere at the source, this is not an issue; the axis of symmetry in the source angles can always be changed so that there are no rays with take-off angles at the poles. However, in general, this is not practical for WFC methods, which should allow for rays with any take-off angle to be traced.
This problem was recognised by Lee & Gibson (2007), who suggested the use of a cubed sphere mesh at the source and the interpolation of rays within quadrilateral cells, so that ray derivatives are stable over the whole range of source angles. However, this approach limits the allowable distribution of rays to quadrilateral cells on the cubed sphere mesh.
Our aim is to implement a WFC method which can have a near arbitrary distribution of ray take-off angles. So, rather than trying to estimate wave front curvature from the spatial derivatives of the ray parameters, we extend the 2-D approach of Ettrich & Gajewski (1996) to 3-D media.
In the method of Ettrich & Gajewski (1996), radial paths are projected from the two basis rays that we wish to interpolate between to locate a virtual source (see Fig. 1A). The circular arc with the basis rays at either end and the virtual source at the centre is used to approximate the wave front. The extension of this approach to 3-D media is straightforward, the two basis rays can still be considered to have come from a virtual point source, and the wave front can be approximated as being ‘locally spherical’ (as opposed to ‘locally circular’ in the 2-D case). The arc between the two basis rays, with its focus at the virtual point source, then traces a great circle along the locally spherical wave front.
Illustration of the ray interpolation method. (A) Two rays heading into a high velocity region with the same velocity at either end of the interpolation arc, (r1−r2). The new ray position is derived by finding the focusing point, f, and positioning half-way along the circular arc r1−r2. (B) Two rays entering a low velocity region. In this case, it is necessary to backpropagate the rays from r′1 and r′2 to r1 and r2, so that the interpolation can be performed in an area without a significant velocity gradient.
The procedure is illustrated in Fig. 1(A); if a propagation step produces a distance between two rays, d = ǀr′2−r′1ǀ, which exceeds the pre-defined maximum, we attempt interpolation on the previous wave front, (i.e. time-step r1, r2). We locate a virtual source or focusing point f by finding the closest point of approach between the two rays, using their current positions and ray velocities. The new ray position is then placed half way along the spherical arc r1−r2.
The ray unit vector for the new ray is interpolated using a weighted average of the slowness vectors of the two basis rays. The Gaussian weights are determined by the difference between the velocity at the interpolation position and velocity at the position of the basis rays.
If the two points/rays we are interpolating between have the same velocity (as in the illustration in Fig. 1A), then the wave front will in fact be spherical between them, and the interpolation will be exact. However, if there is a velocity gradient between the two rays, the wave front will depart from sphericity, and accordingly the locally spherical wave front approximation is not valid. Such a situation is shown in Fig. 1(B); here there is a significant difference in the velocities at the point where the critical ray distance is exceeded (ǀr′2−r′1ǀ). In these cases, the basis rays are propagated backwards until there is a sufficiently small difference in velocity (<2.5 per cent) at the two ray points. In this case, the interpolation is actually performed at, r1−r2, three time-steps behind the wave front, where the velocity is constant between the two rays and the spherical wave front approximation is valid. The new ray is then propagated back to the original wave front.
Previous authors (e.g. Vinje et al. 1996) restrict the interpolation of new rays to the wave front prior to that where the interray distance threshold was exceeded. In these cases, an additional criterion for the interpolation of new rays, based on a maximum allowable divergence angle between rays, is required to prevent interpolation of new rays in regions of rapidly varying velocity. However, the backpropagation of rays into a sufficiently homogenous region prior to interpolation means that monitoring the divergence angle between rays is not required in our approach.
2.2 The First Arrival Wave Front
There are many cases in seismic imaging where only first arrivals are needed, and one potential draw back of WFC is that, in heterogeneous models, wave front folding (caustics) is common. This produces large numbers of rays and slows down the modelling procedure. Accordingly, we can reduce the computational requirements for WFC by only following the first arrival wave front.
The construction of a first arrival wave fronts in 2-D was described in one earliest WFC publications (Vinje et al. 1993); however, until now, 3-D implementations have not expanded on this. Our approach for constructing a first arrival wave front in 3-D media is similar in principle to that used for the 2-D case. That is, the wave front is monitored for wave front folds (caustics) and their associated secondary arrival rays. The secondary arrivals are removed from the network and new links are made joining the rays either side to prevent gaps forming in the wave front. In the 2-D case there will be only two rays bounding the region of removed wave front, making updating the network a relatively simple matter. However, for a 3-D medium, an additional complication arises as many rays can surround a region of removed wave front.
Our procedure for producing a first arrival wave front begins by checking the ray field for the secondary arrivals, which are created when rays cross the boundary formed by their neighbours over a propagation step (ray x in Fig. 2A). The ray is deleted from the network using the same process used to remove rays which exit the model space (Fig. 2B). The resulting wave front contains a gap, which is patched over to produce a first arrival wave front in Fig. 2(C). Filling in the gap in the wave front is equivalent to re-triangulating the deleted ray's bounding polygon (the blue hexagon a, b, c, d, e, f in Fig. 2). In principle it may be possible to use the bounding polygon at the wave front; however, in practice it is more convenient to project the bounding polygon on to a 2-D surface. The natural 2-D surface is the sphere associated with the ray(s) take-off angles at the source, as the triangular network is always defined and non-overlapping in this coordinate system. After re-triangulation, the updated network is projected back on to the wave front to produce an approximate first arrival surface. In cases where adjacent rays have both formed secondary arrivals, the ray removal and re-triangulation can be applied successively, as the ray distribution always forms a network of regular polygons at the source.
Illustration of the technique for generating a first arrival wave front. The left-hand side shows the arrangement of the the rays in the source coordinate system (i.e. the take-off position for each rays on the unit sphere). The right-hand side shows the corresponding positions on the wave front. (A) The ray segment x–x′ passes through a caustic forming a secondary arrival. (B) The ray is removed from the triangular network. (C) The gap in the wave front created by the removal of the ray is patched by re-triangulating the ray's bounding polygon (a, b, c, d, e, f).
The stability of the re-triangulation is limited by the numerical precision of the source angles. In the examples presented here, source coordinates were stored as 64-bit floats. However, more complex wavefields may require increased precision.
2.3 Interpolation of Wavefield Properties
While the wave front is propagating through the medium, the WFC procedure stores points on the current and future wave fronts, as well as the triangular network. Each triangle in the network corresponds to a three sided ray tube (or volume element) with two triangular faces defined by each of the wave fronts. We evaluate wavefield properties at points which lie within the volume element, by dividing each cell into three tetra-hedra and performing a linear interpolation from the surrounding ray points. Examples of wavefield properties which can be interpolated this way include, traveltimes, ray angles, source take-off angles and ray Jacobians.
2.4 Geometrical Spreading Factors
Alternatively, the triangular network used for WFC provides a natural system for evaluating the derivative in eq. (1) directly from the ray geometry. Each element of the triangular network provides a mapping between three points on the wave front and their corresponding take-off angles at the source. Accordingly, we can approximate the derivative in eq. (1) by taking the ratio of the solid angle at the source and the interpolated wave front area. We refer to this method of determining geometrical spreading factors as the ratio method.
The Jacobian method has the advantage that the ray Jacobian can be evaluated along a ray, and so, the quantity can be interpolated in a similar manner to other ray characteristics. However, a potential issue arises where the ray take-off angle is near the pole of the coordinate system. In this case, the calculation of ǀJǀ will become numerically unstable. The ratio method is not sensitive to variations in take-off angle. However, the pseudo ray tubes made during the re-triangulation of the wave front while using the first arrival mode can cause instabilities, especially in the area near caustic surfaces.
3 Results
We begin by examining the accuracy of our WFC technique using a model with a known traveltime function. The model extends for 5000 m laterally and 2500 m vertically and is sampled by points at intervals of 100 m in the horizontal dimensions and 25 m in depth. Fig. 3(A) shows ray points for wave front construction in the model and part (C) shows a vertical profile through the model. As the model is a linear velocity gradient, there is an analytical solution for the traveltime as a function of horizontal offset from the source (for example see Aki & Richards 1980). Accordingly, Fig. 3(B) shows difference in traveltime between the WFC interpolated value and the analytic (or true) value, as a function of horizontal offset from the source for different WFC run parameters. The most accurate results are obtained where the smallest maximum interray distance is used (100 m), with the magnitude of the time-step having a smaller effect on the accuracy. All the traveltime errors are within the magnitude of respective time-step, and in most cases much smaller. Table 1 shows the number of ray steps taken and the CPU time spent propagating the wavefield for the same model runs. Taken together with the traveltime errors, the CPU times provide an indication of the trade-off between accuracy and computation time.
Results of the wave front construction in a velocity model with a linear gradient. (A) Ray points and wave fronts (intervals of 0.2 s) for wave front construction initiated using rays eminating from 5° spherical caps, with a maximum interray distance (Amax) of 200 m and time-step (dt) of 0.02 s. (B) Traveltime error, WFC interpolated traveltime (Twfc) minus the analytically determined value (Ttrue), for a line of receivers at 500 m depth in the model and spaced at 50 m intervals. Results are shown of four different values of Amax and dt. (C) Velocity profile through the model.
Number of rays steps, and CPU time spent propagating the wavefield for the linear gradient velocity model (see Fig. 3).
We also demonstrate our WFC method by showing P-wave simulations for a isotropic 3-D velocity model derived from a North Sea oil reservoir. The model measures 11 × 18 × 4 km and is sampled by points at 200 m intervals laterally and 50 m intervals in depth. The model is 3-D but has a constant velocity in the first (horizontal) dimension; this eliminates out of plane wave propagation effects and makes the interpretation of results easier. As with the previous model (see above), all calculations are carried out for the full 3-D model, but only 2-D cross sections are shown for clarity.
We initiate the model run using 408 rays distributed on 10° spherical caps, emanating from a source 1000 m from the model edge and 500 m from the top. This source is designed to mimic a microseismic event in the subsurface. The geometry is also similar to the procedure for using WFC to generate traveltime tables for use in seismic migration with surface sources and receivers. In the following simulations, we use a ray propagation step of 0.04 s and a maximum distance between adjacent rays of 300 m. Ray quantities are interpolated onto a 17 × 3 km grid, with lateral resolution of 200 m and vertical resolution of 50 m.
We begin by examining time savings produced by tracking only the first arrival wave front (first arrival mode). Example CPU statistics for various parts of the WFC process are shown in Table 2. It can be seen that the use of the first arrival mode results in a factor of 3 reduction in the number of ray steps. This scales directly into savings in the time spent on ray propagation, addition and deletion. The computational overhead associated with the re-triangulation of the wave front is more than offset by not having to calculate KMAH indices for the first arrival wave front. The result, in this case, is an overall saving of some 65 per cent in the amount of time spent managing the wavefield.
Run time statistics for various parts of the wave front construction in the laterally heterogeneous model shown in Fig. 4.
Fig. 4(A) shows the wave fronts generated using the WFC in all arrivals mode, while ray points for the WFC run in first arrival mode are shown in Fig. 4(B). First arrival rays are stopped where the wave front folds back on itself, forming a caustic surface across the middle of the model. It is also worth noting that occasionally new rays are interpolated near the caustic surface, where re-triangulation has taken place, producing non-geometric ray paths. However, as this only occurs in areas of collapsing wave fronts, these features are short lived.
(A) Wave fronts (intervals of 0.4 s) for the WFC run in all arrival mode. (B) Ray points from WFC run in first arrival mode projected from either side of the cross-section line. (C) Contours of the first arrival traveltime, from WFC (black) and an eikonal solver (grey). Contours are spaced by 0.2 s. Colour backgrounds show the profile through the P-wave velocity model. Note the border around the contoured region in (C) is due to the interpolation grid not covering the entire velocity model.
Fig. 4(C) shows traveltime interpolation using the wave front construction method. For comparison, we also show the results from an eikonal solver (FMES). The eikonal solver used the same grid spacing as the wave front construction method, except near the source, where the eikonal solver uses a finer grid (see de Kool et al. 2006, for details). The CPU times for the eikonal solver were 5.57 and 2.27 s for propagation on the refined source grid and the main grid, respectively. In general, the WFC and eikonal solver traveltimes are in good agreement throughout most of the model, and dissimilarities mainly occur near the source or caustic surface.
The interpolated ray angles at grid nodes are shown in Figs 5(A) and (B). In Fig. 5(A), the ray angles trace a clear path to/from the source; however, the rays angles are scattered and appear unstable in the area near the caustic surface, where re-triangulation of the wave front has taken place.
(A) Ray angles for the first arrival interpolated at grid nodes (WFC run in first arrival mode). (B) As before, except using WFC with all arrivals followed. (C) Source take-off angles determined from WFC with first arrival mode, plotted at grid node positions.
Fig. 5(C) shows the variation in take-off angles at the source interpolated on the grid. Alternatively, by reciprocity, these angles give the incidence angle for the wave front from each point in the grid at a receiver located at the WFC source position. Such information is useful for the construction of traveltime tables for seismic migration. Despite the use of the first arrival mode, the source angles do not suffer from the errors seen in the ray angles around the caustic surface. Note that the caustic surface is marked by a discontinuity in take off angle, and the source angles are uniform over a large area, suggesting that first arrivals for a large part of the model are generated by a relatively small range of take-off angles at the source.
In Fig. 6, we show examples of the geometrical spreading factor. Part (A) shows the results derived using the Jacbian method; a major feature is the increased amplitudes beneath the source position. This is in an area where the rays are in fact diverging. In contrast, this is not seen where the geometrical spreading is calculated directly from the ray geometry (ratio method, Figs 6B and C). Instead these plots show focusing occurs down and to the side of the source location.
Geometrical spreading factors interpolated at grid nodes during WFC. (A) Jacobian method with WFC in first arrival mode. (B) Ratio method with WFC in first arrival mode. (C) As (B) but using WFC with all arrivals being followed. The interpolated ray angles as per Fig. 5 are also shown for reference.
Where the ratio method has been used to determine the geometrical spreading factor (Fig. 6B), we see high amplitudes following the caustic surface where the first arrival wave front is patched together. This is not seen with the Jacobian method when used in conjunction with the first arrival mode (part A). The most stable amplitudes are seen where the ratio method was used to calculate geometrical spreading with all arrivals tracked using the WFC (Fig. 6C); neither increased amplitudes beneath the source nor unstable values around the caustic surface are observed in this case. Comparison of parts B and C shows that the use of the first arrival mode creates small differences in the geometrical spreading factor at the far end of the model, away from the source. Comparison with Fig. 4(B) shows that all of the rays in this region have passed through the caustic surface at some stage.
4 Discussion
Our results show that the use of the first arrival mode can produce significant savings in computing resources required for WFC. Although the CPU times mentioned here may seem small, savings will be significant where the WFC is applied to large models or run for large numbers of sources/receivers. The amount of time saved by the use of the first arrival mode depends on how many fewer rays will be traced. Accordingly, the computational saving will vary depending upon the application of WFC and the available hardware. There is a computational overhead associated with the re-triangulation of the wave front; however, this is not much more than what is required for tracking KMAH indices if all arrivals are being followed, so, there will nearly always be some time saving by running in first arrival mode. The extra speed does, however, come at the cost of limiting results to first arrivals and potential inaccuracies in interpolated ray angles and geometrical spreading factors.
The traveltime fields produced agree well with those derived from an independent eikonal solver. Differences near the source are most likely due to limitations in the finite difference eikonal solver, but discrepancies in the traveltime near the caustic surface are probably due to accuracies in the interpolation in the WFC. The use of the first arrival mode means that grid cells in this region can be interpolated between different sections of wave front; hence eikonal solver traveltimes should be considered more reliable in this region. In the case presented here the CPU times for the WFC with first arrival mode were faster than the corresponding values for the finite difference eikonal solver. Although variations in algorithms, models and run parameters for both the WFC and eikonal solver might produce different results, our results suggest that the determination of traveltimes using WFC in first arrival mode is of, at least, comparable speed to finite difference eikonal solvers.
The agreement between the traveltime fields from the WFC and eikonal solver also shows that the ray interpolation method employed here is sufficient for representing the wave front. The locally spherical interpolation method used in this study represents an improvement over the polynomial based methods used in previous 3-D WFC implementations in two main respects. First, the method is physically intuitive and becomes more exact in clearly defined situations, such as the interray distance being reduced or the medium tending towards homogeneity. Second, the method does not require the calculation of geometrical spreading partial derivatives for ray interpolation, which could further reduce the computational cost of ray tracing.
Our results show that use of the first arrival mode produces errors in the interpolated ray angles and geometrical spreading factors near the caustic surface. The effect is not as pronounced in the interpolation of source angles, although it is presumably still present.
Unlike the ray angles, the effect of the first arrival mode on geometrical spreading factors is not entirely localized around the caustic surface. There are some small differences between the geometrical spreading factors derived with and without first arrival mode (Figs 6B and C), at the far end of the model from the source. These can probably also be attributed to the wave front re-triangulation, which took place much earlier on in the ray tube's history.
The comparison of methods for calculating geometrical spreading factors makes it clear that care must be taken when using the ray Jacobian method. The singularity in the source coordinate system produces unreliable results in the region beneath the source. Here, the amplitudes are clearly in error, as we get high values in an area where the ray field is diverging. As only a few initial rays are used, the error in the ray derivatives can be spread over a wide area, resulting in a large area of biased amplitudes. In our results, this masks out the true focusing pattern down and to the side from the source, which is seen when the ratio of wave front area to solid angle at the source is used.
5 Conclusions
In this study, we have described a simple, but effective 3-D WFC procedure, which can be used to interpolate ray quantities throughout the propagation medium. We have used a simple ray interpolation method where the wave front is assumed to be locally spherical. In regions of steep velocity gradient, where this assumption is not valid, rays are backpropagated to a point where the interpolation can be performed in a sufficiently homogenous region. This means that the method does not pose any limitations on the initial distribution of rays, and the calculation of derivatives with respect to take-off angles is not required for wave front interpolation.
The results of our WFC simulations show that the technique can be used to rapidly determine wavefield quantities such as traveltimes, ray angles, source take-off angles and geometrical spreading factors throughout a medium. However, care must be taken when using ray Jacobians to calculate geometrical spreading factors, as the poles of the source coordinate system produce unreliable values, which can be spread over a wide area as only a few initial rays are traced in WFC.
We have also shown how computation time can be reduced by only following the first arrival wave front. For the most part, the use of the first arrival mode leaves the accuracy of the interpolated traveltimes, ray angles and source take-off angles largely unchanged. However, it does have a significant effect on the values interpolated at gridpoints near caustic surfaces, as well as producing small variations in geometrical spreading factors for ray tubes that have passed through caustic surfaces.
Acknowledgments
The authors would like to thank Nick Rawlinson for the use of his Eikonal solver. This work was funded by NERC/BP-Norway partnership grant NE/E006329/1. The authors also benefited from conversations with James Wookey, Doug Angus and Glenn Jones in relation to this work.
References
