Abstract
The linear traveltime interpolation has been a routine method to compute first arrivals of seismic waves and trace rays in complex media. The method assumes that traveltimes follow a linear distribution on each boundary of cells. The linearity assumption of traveltimes facilitates the numerical implementation but its violation may result in large computational errors. In this paper, we propose a new way to mitigate the potential shortcoming hidden in the linear traveltime interpolation. We use the vertex traveltimes in a calculated cell to introduce an equivalent homogeneous medium that is specific to the cell boundary from a source. Therefore, we can decompose the traveltime at a point on the cell boundary into two parts: (1) a reference traveltime propagating in the equivalent homogeneous medium and (2) a perturbation traveltime that is defined as the difference between the original and reference traveltimes. We now treat that the traveltime perturbation is linear along each boundary of cells instead of the traveltime. With the new assumption, we carry out the bilinear interpolation over traveltime perturbation to complete traveltime computation in a 3-D heterogeneous model. The numerical experiments show that the new method, the linear traveltime perturbation interpolation, is able to achieve much higher accuracy than that based on the linear traveltime interpolation.
1 INTRODUCTION
Computing seismice traveltimes from a point source to nodes in a 3-D discrete model is an important part in seismic modelling, traveltime tomography, velocity inversion and migration (Leidenfrost et al.1999; Lan & Zhang 2013; Zhang et al.2014; Huang & Schuster 2015). The traveltime interpolation methods are a class of the grid-based methods for seismic traveltime computation. Asakawa & Kawanaka (1993) proposed a 2-D traveltime computation method based on linear traveltime interpolation (LTI). They showed that the LTI method can be regarded as an advanced version of the finite difference solution of the eikonal equation and has high speed and high accuracy. Li & Ulrych (1993), Cardarelli & Cerreto (2002) and Zhang et al. (2004) improved the LTI methods in 2-D models. Recently, Zhang et al. (2011) extended the traveltime interpolation method to a 3-D model. Huang et al. (2011) proposed 3-D traveltime computation in irregular mesh discretization models using the bilinear traveltime interpolation and wave-front expansion. The LTI methods assume that the traveltime is linear along each boundary of cells, and then the traveltime at a point on a cell boundary can be obtained by linearly interpolating the traveltimes at the adjacent discrete points. However, the traveltimes along a cell boundary may not be linear. For example, the wave front from a point source in a 3-D homogenous model is spherical and thus its traveltime is not linear along a cell boundary. Our numerical tests also show that traveltime distribution appears nonlinear on a cell boundary in a constant-velocity gradient model (see Fig. 2). Therefore, the assumption of linear traveltime may result in large errors in the computed traveltimes likely in the region where wave fronts change rapidly.
Geometry of traveltime computation in a rectangular cell. S is a source, P1, P2, P3, and P4 are the vertices at which the traveltimes are assumed to be already computed. Dash-dotted lines denote the true rays from S to Pi, and solid lines denote the rays in the equivalent homogenous medium. R is a vertex at which the traveltime needs to be computed. P is the intersection of the ray trajectory from S to R and its position needs to be determined.
3-D surface chart of traveltime (ms) on a x–y horizontal cell boundary ranging from 500 to 510 m at a depth of z = 600 m in the constant-velocity gradient model. The left-front vertex on the boundary coincides the projection of the source locating on the top surface of the model. (a) Computed wave-front traveltimes; (b) reference traveltimes in the equivalent medium; (c) traveltime perturbations between (a) and (b).
In this paper, we present a new way to overcome the potential errors that can be induced by the LTI. Starting point of our method is to define an equivalent homogeneous medium between a source and a calculated cell boundary in terms of the vertex traveltimes of the cell boundary. Given a defined equivalent homogeneous medium, reference traveltime can be computed from the source to a point on the cell boundary. Correspondingly, the difference between the computed wave-front traveltime and the reference traveltime is termed as traveltime perturbation. In other words, we express the traveltime at a point on a cell boundary as a sum of reference traveltime and its perturbation. To facilitate traveltime computation, we assume that the traveltime perturbation at a point on a cell boundary can be obtained by linearly interpolating the known traveltime perturbations between the adjacent discrete points. Combining the linear traveltime perturbation interpolation (LTPI) with the group marching method (Kim 2002) for wave-front expansion, we derive a new method for traveltime computation in 3-D heterogeneous isotropic models and demonstrate it with numerical experiments.
2 LINEAR TRAVELTIME PERTURBATION INTERPOLATION
In the grid-based method, a 3-D heterogeneous medium is discretized into a series of regular or irregular cells. Each cell is assumed to be homogenous and a ray within it propagates along a straight path. To illustrate the method, let us begin with the traveltime computation in a cell. In Fig. 1, S is a source, and P1, P2, P3 and P4 are vertices on the left boundary of the rectangular cell. The wave-front traveltimes t1, t2, t3 and t4 at P1, P2, P3 and P4 are assumed to be already computed. Based upon the known taveltimes t1, t2, t3 and t4, next we want to find the traveltime from S to R, a vertex on the right boundary of the cell.
For an arbitrarily inhomogeneous model, the traveltime perturbations are generally non-zero since straight-ray-based reference traveltimes are always different from the original one. For a homogeneous model, the traveltime perturbations will be automatically zero since both reference and true rays completely coincide.
In eq. (6), we decompose traveltime tP into the sum of reference traveltime seqlSP and perturbation traveltime ΔtP. The reference traveltime is easily calculated by definition. We only need to determine ΔtP.
The essence of our new method is to decompose an original traveltime into reference and perturbed traveltimes. We observe that the reference traveltime takes up a major proportion of the original one and the perturbed traveltime a very small part in the decomposition, and both vary nonlinearly on a cell boundary (see Fig. 2). Unlike the LTI where the original traveltimes are assumed to be bilinear on a cell boundary we make an approximation that the perturbed traveltimes vary bilinearly on a cell boundary for computational purpose. Because the magnitudes of perturbed traveltimes remain very small, the reference traveltimes plus linearly interpolated traveltime perturbations still hold the basic property of nonlinear variation of the original traveltimes on a cell boundary. Therefore, we anticipate that it is more reasonable to assume the linearity of traveltime perturbations instead of traveltimes.
For a complex heterogeneous model that is discretized into a series of cells, the wave fronts will propagate outwards from a source throughout the entire model and be sampled at all discrete nodes. The fast marching method (FMM) is an effective way to simulate wave-front expansion (Sethian & Popovici 1999), and the group marching method (GMM) is an optimal version of the FFM (Kim 2002). We thus use the GMM to march wave fronts (for details see Zhang et al.2011).
3 NUMERICAL EXAMPLES
In the numerical tests, we computed the wave-front traveltimes using the LTI method and the LTPI, respectively, and compared their accuracies and computational efficiencies.
3.1 Model with linear increase of velocity with depth
The first model we chosen is a 3-D vertically inhomogeneous medium with a velocity increasing linearly with depth. For this velocity model, the analytical solution of traveltimes and rays from a source to any point is available (Slawinski & Slawinski 1999). Therefore, we can unambiguously evaluate the numerical solutions against the exact one.
The domain of a 3-D model covers 1000 × 1000 × 1000 m3. Its velocity v increasing linearly with depth z is given as v(z) = 1000 + vgz m s−1, where vg is a velocity gradient. We set up a velocity gradient vg = 5 s−1 in this test, which is far larger than those used in other publications (e.g. vg = 3 s−1 in Klimeš & Kvasnička 1994 and Bai et al.2007; vg = 0.8 s−1 in Slawinski & Slawinski 1999). A source is located at the centre of the top surface of the model. The model is discretized into different rectangular cells with the grid spacings 40, 20, 10, 5 and 2 m, respectively. The two methods were executed on a PC with CPU of Intel Core2 Quad Processor Q8400 and Base frequency 2.66 GHz, and RAM 4G.
To understand how the computed wave-front traveltimes, the reference traveltimes and the perturbation traveltimes vary in this model, we use a discrete cell size of 10 m and present the three traveltime quantities on a horizontal x–y boundary ranging from 500 to 510 m at depth z = 600 m. Note that the left-front vertex locates just below the source. What we can see from Fig. 2 are as follows: (1) All three traveltimes vary nonlinearly on the cell boundary; (2) Reference traveltimes capture the majority of wave-front traveltimes in terms of amplitude and variations; (3) In contrast to the reference traveltimes, the perturbations are very small and their nonlinearity in space is weak or mild. For this example, using then LTI to determine point P and its traveltime on the boundary would deal with traveltime variation of about 0.1 ms. On the other hand, using the LTPI to determine point P and its traveltime would only involve with a variation of traveltime perturbation less than 0.02 ms. Thus, small traveltime perturbations make feasible to approximately obtain Δtp via linear interpolation. Meanwhile the total traveltime (reference traveltime plus linearly interpolated perturbation) still holds the basic property of nonlinear variation of the original traveltime on a cell boundary. The LTPI attempts to maintain nonlinear traveltime variation but the LTI simply ignores such variation that can exist physically. This validates our key introduction of an equivalent homogeneous medium.
The average L1 absolute errors between the numerical and exact solutions at all grid nodes for the different grid spacing are summarized in Table 1. From Table 1, we see that the LTI method can induce significant computation errors. From the coarse cell to the small one, the LTI has the average absolute errors ranging from high 10.97 ms to low 1.41 ms. In contrast, the average absolute errors in the LTPI ranges from 1.48 to 0.40 ms. Apparently, the CPU times taken by the LTPI are slightly larger than the LTI for the same cell size. However, to achieve the high accuracy (average absolute error of 0.57 ms) as the LTPI with the cell size of 10 m, the LTI method requires a much smaller cell size than 2 m and would consume enormous computational resources. Fig. 3 shows the exact and computed wave fronts at a horizontal slice with the cell size of 10 m. The wave fronts derived from the LTPI overlap well with the exact ones. But the LTI-based wave fronts systematically lag the exact ones and appear severely distorted near the source region. The wave-front distribution and the absolute error distribution clearly highlights that the LTPI is far superior to the LTI.
The traveltimes (s) and their absolute errors (ms) on the horizontal slice at z = 600 m with grid spacing of 10 m. (a) Traveltime contours. Solid curves indicate the theoretical traveltimes, the dashed ones indicate the traveltimes computed using the LTI, and the dotted ones indicate the traveltimes computed using the LTPI. (b) Absolute errors using the LTI. (c) Absolute errors using the LTPI.
The errors of traveltimes computed using LTPI and LTI and the CPU times.
| LTPI | LTI | |||||
|---|---|---|---|---|---|---|
| Grid spacing | Average absolute | Average relative | CPU times | Average absolute | Average relative | CPU times |
| (m) | errors (ms) | errors (%) | (s) | errors (ms) | errors (%) | (s) |
| 40 | 1.48 | 0.51 | 1 | 10.97 | 3.47 | 1 |
| 20 | 0.79 | 0.27 | 5 | 5.48 | 1.77 | 5 |
| 10 | 0.57 | 0.19 | 38 | 3.54 | 1.15 | 31 |
| 5 | 0.48 | 0.15 | 309 | 2.35 | 0.76 | 244 |
| 2 | 0.40 | 0.11 | 4976 | 1.41 | 0.47 | 3989 |
| LTPI | LTI | |||||
|---|---|---|---|---|---|---|
| Grid spacing | Average absolute | Average relative | CPU times | Average absolute | Average relative | CPU times |
| (m) | errors (ms) | errors (%) | (s) | errors (ms) | errors (%) | (s) |
| 40 | 1.48 | 0.51 | 1 | 10.97 | 3.47 | 1 |
| 20 | 0.79 | 0.27 | 5 | 5.48 | 1.77 | 5 |
| 10 | 0.57 | 0.19 | 38 | 3.54 | 1.15 | 31 |
| 5 | 0.48 | 0.15 | 309 | 2.35 | 0.76 | 244 |
| 2 | 0.40 | 0.11 | 4976 | 1.41 | 0.47 | 3989 |
The errors of traveltimes computed using LTPI and LTI and the CPU times.
| LTPI | LTI | |||||
|---|---|---|---|---|---|---|
| Grid spacing | Average absolute | Average relative | CPU times | Average absolute | Average relative | CPU times |
| (m) | errors (ms) | errors (%) | (s) | errors (ms) | errors (%) | (s) |
| 40 | 1.48 | 0.51 | 1 | 10.97 | 3.47 | 1 |
| 20 | 0.79 | 0.27 | 5 | 5.48 | 1.77 | 5 |
| 10 | 0.57 | 0.19 | 38 | 3.54 | 1.15 | 31 |
| 5 | 0.48 | 0.15 | 309 | 2.35 | 0.76 | 244 |
| 2 | 0.40 | 0.11 | 4976 | 1.41 | 0.47 | 3989 |
| LTPI | LTI | |||||
|---|---|---|---|---|---|---|
| Grid spacing | Average absolute | Average relative | CPU times | Average absolute | Average relative | CPU times |
| (m) | errors (ms) | errors (%) | (s) | errors (ms) | errors (%) | (s) |
| 40 | 1.48 | 0.51 | 1 | 10.97 | 3.47 | 1 |
| 20 | 0.79 | 0.27 | 5 | 5.48 | 1.77 | 5 |
| 10 | 0.57 | 0.19 | 38 | 3.54 | 1.15 | 31 |
| 5 | 0.48 | 0.15 | 309 | 2.35 | 0.76 | 244 |
| 2 | 0.40 | 0.11 | 4976 | 1.41 | 0.47 | 3989 |
3.2 Model with velocity varies in three coordinate directions
The second model is a full 3-D velocity varying medium, given as the fig. 6 in Klimeš & Kvasnička (1994). This model was used in Bai et al. (2007) to test their shortest-path method. The model size is 40 × 40 × 20 km3. The velocity at the top surface is 3 km s−1 everywhere. The velocity at the bottom oscillates between 5 and 6 km s−1, with minima and maxima 10 km apart, and the velocity at a point on the bottom between the minima and maxima is determined by 2-D cubic spline interpolation. The velocity along each vertical line varies linearly between the velocity values at the top and the bottom surfaces, shown in Fig. 4. The source is located at the geometric centre of the bottom of the model. Sixty-six receivers are placed on the top surface of the model (see the triangles in Fig. 4a). Klimeš & Kvasnička (1994) computed the traveltimes at the 66 receivers, using the Complete Ray Tracing program package by Červený et al. (1988), as the exact values for the comparison with the traveltimes obtained with their method. We also use the exact values listed in table 1 in Klimeš & Kvasnička (1994) to compute the traveltime errors of our method.
Model with velocity (km s−1) varies in three coordinate directions. (a) Velocity isolines along the bottom of the model. The triangles indicate the 66 receiver positions projected on the bottom. (b) Velocity isolines in the vertical section y = 5 m.
We discretized the 3-D model into different rectangular cells with the grid spacings 10, 5, 2, 1, 0.5, 0.2, and 0.1 km, respectively. The LTI and LTPI methods were used to compute the traveltimes at the 66 receivers distributed irregularly on the top surface of the model. Fig. 5 shows the average of the L1 absolute relative errors between our computed and the exact values at the 66 receivers versus the grid spacings. From coarse to fine cells, the LTPI outperforms the LTI. For example, with a cell size of 0.5 km, the LTPI achieves the relative error of 0.09 per cent; the LTI 0.81 per cent. The results show that the use of smaller sized cells would not help reduce the error levels of the LTI to those of the LTPI.
Average of L1 absolute relative errors at the 66 receivers for the model in Fig. 4 as the grid spacing is decreased from 10 to 0.1 km. The solid line indicates the errors using the LTPI, and the dotted one indicates the LTI.
4 CONCLUSIONS
We have proposed a new method to compute traveltimes in 3-D complex heterogeneous models based on bilinear interpolation of traveltime perturbations (the LTPI). In the new method, we introduce an equivalent homogeneous medium and decompose the original traveltimes into associated reference traveltimes and perturbations specific to each calculated cell boundary. The reference travleltimes take into account the major portions of the original ones and vary nonlinearly on cell boundary, and the traveltime perturbations generally remain small. The linear approximation to traveltime perturbation on a cell boundary makes computation efficient and allows us to maintain the nonlinearity of traveltime on a cell boundary. This is the fundamental difference from the LTI that ignores the existence of nonlinear traveltime variations. The numerical examinations of the computed traveltimes against the exact ones reveal that the LTPI-based method is far superior to the LTI-based method in terms of accuracy and computational time. Next, we generalize the LTPI to the case of irregular cells so that we can effectively simulate more complex models in practice.
This work was supported by the National Natural Science Foundation of China (Grant Nos 41074077 and 41230318), and the Specialized Research Fund for the Doctoral Program of Higher Education (Grant No. 20130132110023). We thank the editor, Prof. Jeannot Trampert, the reviewer, Dr Yunsong Huang, and an anonymous reviewer for their constructive comments that help improve the paper.
REFERENCES
