ABSTRACT

Seismic modeling in heterogeneous media is accomplished by using either approximate or fully numerical methods. A popular approximate method is ray-Born modeling, which requires the computation of 3D integrals. We have developed an integration technique for accurate and, under certain circumstances, efficient evaluation of the ray-Born integrals in the time domain. The 3D integrals are split into several 2D integrals, each of which gives the wavefield at a certain time, so that the waveform at each time step is computed independently of all other times. We compute seismograms for 3D heterogeneous acoustic media using this technique and compare these seismograms with seismograms computed using two other modeling methods: frequency-domain ray-Born modeling and finite-difference modeling of the acoustic wave equation. Our method can also be applied to elastic ray-Born modeling. Velocity models with smooth scatterers and the SEG/EAGE overthrust model are used for comparison. The ray-Born seismograms computed using the time- and frequency-domain ray-Born modeling methods are identical, as expected. The comparison between the ray-Born modeling and the finite-difference-modeling method indicates that the waveforms are similar for both types of velocity models. We evaluate the discrepancies in terms of multiple scattering and multipathing.

INTRODUCTION

Modeling of seismic waves is very important for studying the Earth’s structure. For example, seismic algorithms such as full-waveform inversion require it. In practice, this modeling is often done using the acoustic wave equation, solved either by using fully numerical methods or by using approximate methods (e.g., Carcione et al., 2002; Virieux and Lambaré [2015], chapter 1.05). Fully numerical methods are popular and very useful, but they are not computationally efficient at higher frequencies, especially in large 3D models. They also have limitations in the case of complicated structural models (e.g., finite-difference seismograms contain artifacts in the case of strongly undulating interfaces). Approximate methods do not necessarily suffer from these disadvantages; for example, ray theory is particularly valid at high frequencies in smooth models. Therefore, it is sometimes beneficial to compute wavefields using approximate methods.

A useful technique for modeling singly scattered waves is the (first-order) ray-Born approximation. Several papers have investigated the theoretical aspects of the Born (e.g., Knopoff and Hudson, 1964; Hudson, 1977; Hudson and Heritage, 1981) and ray-Born approximations (see references listed by Chapman [2004], section 10.3 and Červený [2005], section 2.6.2). In studies focusing on the numerical aspects of ray-Born modeling, Červený and Coppoli (1992) generate synthetic seismograms in the case of a background medium with two curved interfaces and a perturbation consisting of isolated scatterers. For a model with similar characteristics, Gibson et al. (1993) use ray-Born modeling to examine reflections from vertical seismic profile data. Coates and Charrette (1993) compare the 2D isotropic elastic generalized Born approximation (Coates and Chapman, 1991) and the conventional Born approximation to finite-difference modeling. In a study focusing on 3D ray-Born inversion, Thierry et al. (1999) compute synthetic traces from a migrated image and compare the results with real data. In a review paper, Moser (2012) pays special attention to the advantages of acoustic ray-Born modeling for modeling diffractions from specific structural discontinuities, such as edges and tips. Šachl (2013) investigates discretization effects of the elastic ray-Born scattering integral and suggests a cosine window to attenuate boundary diffractions. Tengesdal (2013) and Minakov et al. (2017) discuss 2D acoustic ray-Born and finite-difference modeling and show synthetic- and real-data examples as well as an application to full-waveform inversion. All of the above papers compute the ray-Born integrals, in either the time or frequency domain, by summing over individual volume or area elements (e.g., Coates and Charrette, 1993). An alternative method for the computation of ray-Born seismograms in the time domain was presented by Spencer et al. (1997) and Keers et al. (2002). That method used a specific integration over tetrahedrons.

We present a new method for 3D ray-Born modeling in the time domain using isochron surfaces. Integrals over isochrons are used to accurately compute seismograms for a band-limited signal. We compare the results with seismograms computed using frequency-domain ray-Born and finite-difference modeling. The method can be seen as an extension to 3D integrals of what was done for 1D Wentzel-Kramers-Brillouin-Jeffreys (WKBJ) seismograms by Chapman (1978), where the band limiting is done before the discretization and for 2D Kirchhoff seismograms by Haddon and Buchen (1981), where the integrals are computed over isochrons.

The results in this paper are based on the acoustic wave equation, but the concepts can be extended to the elastic case. We give a brief description of this in the “Discussion” section.

In the “Theory” section, we review the acoustic ray-Born modeling and give the theoretical background of the integration technique. Next, we discuss the numerical implementation of the ray-Born integrals. For the sake of completeness, this section also gives a brief description of the two-point ray-tracing algorithm as well as the frequency-domain ray-Born integrals and finite-difference modeling. In the “Results” section, we present seismograms computed using two types of velocity models with different scattering characteristics. The first type has smooth velocity perturbations, analogous to that used in ray-based tomographic model building (or full-waveform inversion with low-frequency velocity updates). The second type, the 3D SEG/EAGE overthrust model of Aminzadeh et al. (1997), contains rough and rapidly varying perturbations.

THEORY

For a heterogeneous velocity model c(x), the propagation of acoustic waves can be described by a Green’s function G(x,s,t), where x denotes the spatial position of observation, s denotes the source position, and t denotes time. The term G satisfies the wave equation (e.g., Morse and Feshbach, 1953; Ikelle and Amundsen, 2005)
c2(x)Gtt(x,s,t)ΔG(x,s,t)=δ(xs)δ(t),
(1)
where the source is a Dirac delta function in time and space, which acts at time t=0, Gtt=2G/t2 and Δ is the Laplace operator. Expressions for the solution of equation 1 are often more convenient in the frequency domain than those in the time domain. For this reason, we also give the frequency-domain formulation, with ω denoting the angular frequency, and we adopt the sign convention for the Fourier-transform pair from Chapman (2004, pp. 58 and 59). The wave equation in the frequency domain becomes the Helmholtz equation:
ω2c2(x)g(x,s,ω)+Δg(x,s,ω)=δ(xs),
(2)
where g(x,s,ω) denotes the Green’s function in the frequency domain.
In forward modeling, equation 2, for a given c, needs to be solved for g. If perturbation theory is used to solve equation 2, then one decomposes c into a background velocity c0 and perturbation term c1:
c=c0+c1.
(3)
Often, the velocity field c0 is smooth, or it is only a function of depth, and c1 is more rapidly varying. The solution to equation 2 for c=c0 is denoted by g0. The modeling problem then is to find g when g0 is known. Using perturbation theory, an expression for g in terms of the background Green’s function g0 (e.g., Červený, 2005; Ikelle and Amundsen, 2005) can be found:
g(x,s,ω)=g0(x,s,ω)+ω2Dg0(x,x,ω)V(x)g(x,s,ω)dx.
(4)
Here, the “scattering term” V is defined as
V=c02c2,
(5)
and D is the scattering region, i.e., the region where V is nonzero. In compact notation, if we define the operator L as
(Lg)(x)=ω2Dg0(x,x,ω)V(x)g(x,ω)dx,
(6)
then equation 4 can be rewritten as
g=g0+Lg
(7)
or
g=(IL)1g0,
(8)
where I is the identity. The operator (IL)1 can be expanded, assuming that it is invertible, as a power series in terms of L (e.g., Griffel, 2002):
g=(I+L+L2+L3+)g0.
(9)
Here, each term describes scattering of a certain order: L describes the first-order scattering, L2 describes the second-order scattering, and so on. The series is convergent if L<1 for some suitable norm. This is often the case, especially if V is sufficiently small and D is not too large.
In practice, the explicit expression for g as given by the scattering series 9 is approximated using only the first two terms. This gives the Born approximation to the wavefield:
gg0+Lg0=g0+g1.
(10)
In the above expression, the singly scattered waves g1 are given as a volume integral over wavefields propagating through the background medium:
g1(x,s,ω)=ω2Dg0(x,x,ω)V(x)g0(x,s,ω)dx.
(11)
As it is assumed that c1c0, we can rewrite V:
V=1c021(c0+c1)22c1c03.
(12)
If we furthermore use r to indicate receiver position, then the Born-scattering integral 11 becomes
g1(r,s,ω)ω2Dg0(r,x,ω)2c1(x)c03(x)g0(x,s,ω)dx.
(13)
A clear physical interpretation of equation 13 can be given (e.g., Snieder, 2001): The scattered wave g0(x,s,ω) propagates in the background medium from the source located at s to a scatterer at point x, interacts with the scatterer, and then propagates as g0(r,x,ω) from the scatterer to the receiver at r. The total scattered wavefield g1 is then found by summation of the integrand over all scattering points x and multiplication by ω2.

Background Green’s function

The Born integral 13 requires knowledge of g0. In the case of a homogeneous background medium, g0 is known explicitly. However, this situation is geologically not realistic. If we assume a slowly varying c0, then it is natural to use ray theory (Červený, 2005) to compute g0. In such a case, we look for solutions of the form
g0(x,s)=A(x,s)eiωT(x,s),
(14)
and we solve for traveltime T and amplitude A. Equations for T and A are found by inserting 14 into the wave equation and applying the high-frequency approximation (ω). This gives the eikonal equation for T
|T|2=c02,
(15)
and the transport equation for A
AΔT+2A·T=0.
(16)
Characteristics of the eikonal equation 15 are raypaths, described by the kinematic ray equations (Červený, 2005), and they can be used to find T. A solution to the transport equation 16 can be found by solving the dynamic ray-tracing equations (Červený, 2005). Note that it is assumed that there is no multipathing.

The ray-Born approximation

If the background Green’s function g0 is obtained from ray theory, then the Born approximation is called the ray-Born approximation. In this case, the ray-Born integral (e.g., Beylkin, 1985) becomes
g1(r,s,ω)=ω2DArs(x)eiωTrs(x)V(x)dx,
(17)
which in the time domain is
G1(r,s,t)=d2dt2DArs(x)δ(tTrs(x))V(x)dx.
(18)
Here, the traveltime function Trs and amplitude function Ars are defined as
Trs(x)=T(r,x)+T(x,s)
(19)
and
Ars(x)=A(r,x)A(x,s).
(20)
The integral in equation 18 can be computed using the composition of a function and the Dirac delta. This gives
G1(r,s,t)=d2dt2StArs(x)|Trs(x)|V(x)dx,
(21)
with St={x|t=Trs(x)}. For a slowly varying c0, St is a closed 2D surface, which envelopes the source s and receiver r. We call this surface an isochron. From equation 21, we see that at time t the ray-Born seismogram is determined only by the isochron St and the values of Ars, Trs, and V on St.

Later on in this paper, we briefly discuss multiple scattering and multipathing, and it is important to distinguish the two effects. Multipathing occurs when, given a velocity model c and source s, there are two or more rays between s and the receiver r. Multiply scattered waves form that part of the wavefield that is modeled using the higher order terms (L2g0, L3g0, etc.) of the scattering series 9. Therefore, multipathing depends on the source s, the velocity model c, but also the background velocity c0. Singly scattered waves are the waves modeled by the term g1=Lg0. Note that multipathing and multiple scattering are not mutually exclusive, so that a part of the wavefield can be interpreted in terms of multipathing and in terms of multiple scattering.

Discrete ray-Born seismograms in the time domain

We start with the expression for the ray-Born approximation
G1(r,s,t)=d2dt2Dars(x)δ(tTrs(x))dx,
(22)
with
ars(x)=Ars(x)2c1(x)c03(x),
(23)
and Trs as given in equation 19. In practice, seismograms are only recorded and computed at discrete equally spaced time points. It is therefore natural to compute seismograms at these times only. The discretization can be done in various ways. One way is to convolve G1 with the Boxcar function B(t) given by the Heaviside step function H(t) as
B(t)=H(t+Δt/2)H(tΔt/2).
(24)
This gives the band-limited integral, where the support is bounded
G^1(r,s,t)=d2dt2Vt(x)ars(x)dx,
(25)
with
Vt(x)={x|tΔt/2<Trs(x)<t+Δt/2}.
(26)
We leave the differential d2/dt2 to be evaluated numerically later. The integral over Vt(x) is the volume between two isochron surfaces StΔt/2 and St+Δt/2 with St±Δt/2={x|t±Δt/2=Trs(x)} (see Figure 1). Note also that Vt=Sτ with τ[tΔt/2,t+Δt/2], i.e., Vt is the union for all isochron surfaces and Sτ1Sτ2=Ø for τ1τ2. Moreover, because there is no multipathing, all surfaces Sτ are well-defined and smooth. For a homogeneous 3D background medium, the isochrons are ellipsoidal (see Figure 2). The ellipsoids become spheres in the special case when s=r. If the background medium is slowly varying, then the surfaces Sτ are deformed, but they are still smooth and ellipsoidal with Sτ+Δt completely enveloping Sτ.
A natural parameterization of Vt is given in terms of coordinates (ξ1,ξ2), which parameterize the surface Sτ and coordinate y. For fixed (ξ1,ξ2), the isochron raypath x(ξ1,ξ2,y) is approximated by a line perpendicular to Sτ. Because of this, the map (x1,x2,x3)(ξ1,ξ2,y) is well-defined and smooth, and so is its Jacobian
J=|(x1,x2,x3)(ξ1,ξ2,y)|.
(27)
Therefore, we have
G^1(r,s,t)=d2dt2Vt(x)ars(x(ξ1,ξ2,y))×|(x1,x2,x3)(ξ1,ξ2,y)|dξ1dξ2dy.
(28)
Because
|(x1,x2,x3)(ξ1,ξ2,y)|=|xy·(xξ1×xξ2)|
(29)
and (x/y)(x/ξi), with i=1 or 2, this can be rewritten as
G^1(r,s,t)=d2dt2Vt(x)ars(x(ξ1,ξ2,y))×|xy||xξ1×xξ2|dξ1dξ2dy.
(30)
We now assume that, for fixed (ξ1,ξ2), ars is constant and also |(x/ξ1)×(x/ξ2)|. We choose these constants, such that ars(x(ξ1,ξ2,y))=ars(x(ξ1,ξ2,y˜)), where y˜ is the value for St, and similarly for |(x/ξ1)×(x/ξ2)|. We denote these constants for y by ars(x(ξ1,ξ2)) and |(x/ξ1)×(x/ξ2)|. It should be stressed that ars varies over each isochron (as a function of (ξ1,ξ2)). In the integration over y in equation 30, it is only assumed that ars is constant as a function of y when going from one isochron to the next isochron. We now have
G^1(r,s,t)=d2dt2Stars(x(ξ1,ξ2))|xξ1×xξ2|×[|xy|dy]dξ1dξ2.
(31)
The integral |x/y|dy is found using the length of the line segment 2Δy (see Figure 1) and from the definition of St and St+Δt/2. This gives
Δy=Δt2|pr+ps|.
(32)
Here, the slowness vectors from the source and receiver at x are denoted by ps=T(s,x) and pr=T(r,x), respectively. Furthermore, using
|pr+ps|=2cosφc0,
(33)
with φ being the half the angle between pr and ps (e.g., Miller et al., 1987), 2Δy is, to first order
2Δy=c0Δt2cosφ.
(34)
Inserting this in equation 31 and using the definition of ars gives
G^1(r,s,t)=d2dt2StArs(x(ξ1,ξ2))Δt|cosφ(x(ξ1,ξ2))|×c1(x(ξ1,ξ2))c02(x(ξ1,ξ2))|xξ1×xξ2|dξ1dξ2
(35)
or
G^1(r,s,t)=d2dt2StArs(x)Δt|cosφ(x)|c1(x)c02(x)dx.
(36)
It should be emphasized that equation 36 is band-limited and therefore scales when Vt changes as a function of Δt. Note that when φ0, for instance, near zero-offset, the integral can be approximated by
G^1(r,s,t)=d2dt2StArs(x)Δtc1(x)c02(x)dx.
(37)
Similarly, for forward scattering, we have φπ, so that equation 37 also holds in this case.

For t=t1, it is not possible to integrate along an isochron because the isochron for t=t1Δt does not exist. In this case, there are two different cases that must be distinguished: (1) t1Δt<Trs<t1Δt/2 and (2) t1Δt/2<Trs<t1. In both cases, the ray-Born integral is given in terms of a volume integral as in equation 25. In our experience, these contributions in practice are not significant. However, this should be included here for the sake of completeness. The numerical implementation of equation 36 and the volume integrals are discussed in the next section.

NUMERICAL IMPLEMENTATION

Ray tracing

In preparation for ray-Born modeling using equation 36, the quantities Trs, Trs, and Ars are computed using ray tracing. They are found by first computing T(s,x), T(s,x), and A(s,x) from the source s to all points x in the scattering region D. Likewise, T(r,x), T(r,x), and A(r,x) are computed from the receiver r to the scattering points x.

If the velocity model is slowly varying, raypaths, traveltimes T, and wavefront normals T are found by solving the kinematic ray equations given the initial conditions (position and take-off angles). In a similar manner, the amplitude is computed by solving the dynamic ray equations. Both sets of equations are solved numerically with a fourth-order Runge-Kutta method (Press et al., 2007) with Snell’s law applied when the ray encounters an interface. The result is two 3D fans of rays emanating from the source and the receiver.

For a given source s and scattering point x, the two-point ray-tracing problem is to find a ray from s to x. This is done in two steps. First, all points on the ray fan are used to construct a 3D Delaunay tessellation (de Berg et al., 2000), which covers the scattering region D. For each x in D, we then interpolate values for T, T, and A using the tessellation and linear interpolation on a tetrahedron (e.g., Farin, 2002). Likewise, the same procedure is done for a receiver r and scattering point x.

With the two-point ray tracing completed for the source-receiver pair, we compute the traveltime function Trs using equation 19 and the amplitude function Ars using equation 20.

Ray-Born seismograms: Time domain

The isochron can be computed in various ways. In a heterogeneous background with the traveltime function Trs given on a discrete grid, one of several isosurface extraction techniques can be used (e.g., Sutton et al., 2000). We have computed Trs using the ray-tracing procedure described above, and we have used the marching cubes algorithm (Lorensen and Cline, 1987) to extract surfaces. In this case, Sτ is defined on a triangular grid with vertex positions x1=x1(n), x2=x2(n), and x3=x3(n), where n=1,,ntr and ntr is the number of triangles.

The numerical integration on the triangular grid is carried out by computing an average value for the integrand at each triangle times the area of each triangle. We evaluate the integrand 36 at vertex positions xm(n) with m=1,2,3 and compute an average value In as
In=Δt3m=13Ars(xm)|cosφ(xm)|c1(xm)c02(xm).
(38)
Furthermore, we denote the area element of a triangle ΔAn and calculate this as half the area of the parallelogram spanned by x2x1 and x3x1:
ΔAn=12|(x2(n)x1(n))×(x3(n)x1(n))|.
(39)
Having at our disposal the area ΔAn and integrand In on the triangular grid, we compute the integral for Trs=ti as
G^1(r,s,ti)d2dt2n=1ntrInΔAn.
(40)
The second derivative is computed once for each source-receiver pair using numerical differentiation. In the following sections, we show seismograms computed using equation 40.

The volume integral 25 for the two special cases mentioned in the previous section can also be evaluated numerically. For this, the vertices of the first isochron as well as the points on the background ray, which is completely enclosed by the isochron, can be used. In this case, the integral 25 becomes a sum of the volume contributions of the tetrahedrons n=1ntetI^nΔVn, where ntet is the total number of tetrahedrons, ΔVn is the volume of tetrahedron n, and I^n is the integrand of tetrahedron n (which may be computed as the average of the integrand on the vertices of each tetrahedron).

Ray-Born seismograms: Frequency-domain comparison

The computation of the ray-Born approximation in the frequency domain is done using equation 17. The discretization of the integral is accomplished using a regular grid over the scattering region D. By approximating the integrand in each grid element by a constant, the computation of g1(r,s,ω) becomes a sum over these integrand values, multiplied by the volume element.

The comparison with the time-domain ray-Born modeling method and finite-difference modeling is done in the time domain. Therefore, after computation of g1(r,s,ω) for all frequencies, an inverse fast Fourier transform is used.

Finite-difference comparison

The finite-difference method (e.g., Etgen and O’Brien, 2007) is used to solve the wave equation 1. For our implementation, we do not include any absorbing boundary conditions. To avoid edge effects within the time window of interest, all velocity models are extended beyond the area of interest.

RESULTS

In this section, seismograms computed using the three modeling methods (time-domain ray-Born using isochrons, frequency-domain ray-Born, and finite-difference) discussed in the previous section, are presented for two different velocity models. The first of these velocity models is a smooth model with random reflectivity. The second one is the 3D EAGE/SEG overthrust model (Aminzadeh et al., 1997). This is a complicated velocity model with a rapidly varying velocity. Small differences between the computed seismograms are assessed using single-trace displays. Larger scale differences, with significance for structural interpretation, are assessed using shot records and zero-offset sections.

3D random Gaussian model

First, we present the modeling results for the velocity models with random variations. These random velocities are generated with a Gaussian autocorrelation function characterized by a perturbation strength ϵ (root-mean-square slowness variation) and correlation distance a (Baig and Dahlen, 2004). We use this to create a random realization with weak, medium-strength, and strong perturbations (details are given in Table 1). The correlation distance in all cases is a=50  m and the velocity model covers a volume of 1.8×2.4×2  km. A constant background velocity of 2000  m/s is used. Figure 3 shows two cross sections of the particular realization that was used.

The model grid spacing is 10 m, and this is used for the ray-Born-modeling methods as well as the finite-difference-modeling method. This way, amplitudes, traveltimes, and velocities are all sampled on the same grid. The models are extended in all directions to avoid interference from artificial boundary reflections within the recording length of 2 s. The sampling interval used was 4 ms. Receivers are positioned at a 50 m interval from 0 to 2350 m in the y-direction and at depth of 100 m (Figure 3). Synthetics are generated for a shot located at s=(800,1200,200)  m, using a source time function in the form of a causal Ricker wavelet, with a center frequency of 15 Hz.

For three discrete sample points in time, ti=(0.65,0.85,1.05)  s and a source-receiver distance of 800 m, Figure 4 shows the corresponding isochron surfaces Trs=ti. The values on the surfaces are the values of the integrand of equation 36 for the weak perturbation model. As the background medium is homogeneous, the isochron surfaces are part of a spheroid.

Figure 5a (0.5% perturbation) shows a trace offset by 800 m relative to the source position (see the dashed blue line in Figure 6 for the location of the receiver) for the weak perturbation model (Table 1). We show the finite-difference seismogram and ray-Born seismograms computed in the time and frequency domains. Next, we show traces for the same random realization, but with medium-strength and strong perturbation models in Figure 5b (10% perturbation) and 5c (15% perturbation), respectively. For the models with 0.5%, 10%, and 15% perturbation, the ray-Born seismograms all appear with the same phase, but with a relative difference in the amplitudes. The results computed in the time and frequency domain are basically identical, as expected. In the comparison with finite-difference modeling, the seismograms appear with changes in the waveforms for the medium-strength (Figure 5b) and more so for the strong perturbation model (Figure 5c). This is because finite-difference modeling accounts for complicated wave-propagation effects such as multiple scattering and multipathing that become more predominant for stronger perturbations.

Figure 6 shows synthetic shot gathers for the medium-strength perturbation model. The synthetics generated by the ray-Born method and finite differences are shown in Figure 6a and 6b, respectively. Across the offset positions, we can recognize events from the singly scattered wavefield in both seismograms.

3D SEG/EAGE overthrust model

In this subsection, time-domain ray-Born, frequency-domain ray-Born, and finite-difference synthetics are computed using the acoustic 3D SEG/EAGE overthrust model (Aminzadeh et al., 1997; Operto et al., 2003). The overthrust model represents a succession of sedimentary deposits in a compressional tectonic setting, resulting in strong vertical and lateral velocity variations. It was created to study 3D wave propagation and to test seismic processing algorithms.

Two minor modifications are made to the original model: Edges are extrapolated to avoid artificial boundary reflections, and a seafloor interface is added at a depth of 500 m using a water layer velocity of 1480  m/s. A subvolume of the larger model is displayed as a vertical cross section in Figure 7a and at a depth of 2 km in Figure 7b.

We create the background model by 3D boxcar averaging of the model without the water layer. The background model is shown in Figure 7c and 7d. The difference cc0 then gives the perturbation model, which, in Figure 7e and 7f, is displayed in terms of percent relative to the background. Perturbation strengths have a maximum of 41.1% and a minimum of 21.7%. On average, the positive velocity perturbation is 6.8% and the average value of the negative velocity perturbation is 5.1%. In contrast to the velocity models of the previous section, the perturbation model contains sharp velocity variations, and, therefore, it is quite rough.

The grid spacing is 12.5 m along all spatial axes, and it is the same for frequency-domain ray-Born-modeling and finite differences. To achieve an exact match between time-domain ray-Born modeling and frequency-domain ray-Born modeling, we reduced the vertical grid spacing to a 6.25 m space for the time-domain implementation. This allows for a more precise isochron surface extraction. The source-time function used for the computation of seismograms is defined by a causal Ricker wavelet with a center frequency of 10 Hz. In the case of the ray-Born modeling methods, we include the direct wave traveling in the water layer, but we do not calculate the seafloor reflection.

With a water layer, take-off angles for ray tracing are limited by the critical angle at the water-sediment interface. Therefore, we do not have any real valued raypaths beyond these angles. This limited ray coverage creates an artificial boundary in the scattering volume, with an associated end-point signal (Bender and Orszag, 1999). To avoid these artifacts, traveltime and amplitude values are interpolated by including points just above the seafloor from ray tracing in the water layer, and adding these to the collection of spatial locations from one-point ray tracing before the Delaunay tessellation process. Figure 8 shows traveltimes and amplitudes for a source at s=(2500,5000,200)m after interpolation on the uniform grid. Figure 8 also shows the interpolated ray shadow zone, the triangular wedges between the dashed line and the seafloor (the solid line).

For a source positioned at s=(2500,5000,200)  m and a receiver positioned at r=(2500,5000,100)  m, Figure 9 shows isochron surfaces in the overthrust background model for times Trs=1.8  s (Figure 9a) and Trs=2.36  s (Figure 9b). Values on the surfaces are computed with the integrand in equation 36. Because the background model is heterogeneous, the isochron surfaces look like deformed spheroids as expected.

Ray-Born and finite-difference zero-offset sections

We show an example of zero-offset synthetics computed across a region with complicated faulting as well as a structurally simpler area with near-horizontal layers. Sources (at a depth of 100 m) and receivers (at a depth of 200 m) are positioned from 500 to 5000 m in the x-direction and at 5000 m in the y-direction. This is indicated by the dotted black line in Figure 7b. With an x, y aperture of 4000 m, the wavefield is computed using the two ray-Born-modeling methods and the finite-difference-modeling method.

Seismograms, computed using finite-difference and ray-Born modeling, are shown in Figure 10, for three positions in the overthrust model. The three positions are at 1000, 2500, and 4000 m, as indicated with the dashed blue lines in Figure 11. The seismograms in Figure 10a (1000 m) are extracted from a region with less complicated wave-propagation effects (the velocity model has horizontal and dipping layers). Here, the comparison with the finite differences shows that the phase and amplitude appear to be well-matched, indicating the relative importance of singly scattered waves. We also compare the seismograms shown in Figure 10b (2500 m). In this case, the source and receiver are located in a region where the velocity model has steeply dipping layers and faulted structures. These structures cause complicated wave-propagation effects (such as multiply scattered events) that appear as differences in the modeling results after traveltimes of approximately 1.4 s. In Figure 10c (4000 m), the early-time waveform mismatch in the time interval 1−1.4 s is caused by a triplication in the wavefield that is formed by the shallow syncline. Perturbations at the syncline are relatively strong (up to 30%) and are therefore not accurately captured by the ray-Born approximation.

A zero-offset section for the ray-Born synthetics is shown in Figure 11a and for finite differences in Figure 11b. Strong reflections from the horizontal and dipping layers are indicated in the blue rectangle labeled by I, confirming what was discussed above for Figure 10a. The wavefield becomes more complicated below the anticline, indicated by the blue rectangle labeled II. In this area, ray-Born synthetics exhibit stronger reflections than the finite-difference synthetics. This is likely an effect of multiple scattering and multipathing: Internal reflections in the finite-difference synthetic interfere with the main reflections, and these internal reflections do not exist in the case of ray-Born modeling. In the blue rectangle labeled III, we have indicated the triplication effect that causes the mismatch at early times in Figure 10c. This particular effect is present in the finite-difference synthetic, but poorly modeled in the ray-Born synthetic.

DISCUSSION

The main results of this paper are the derivation and computation of ray-Born integrals as integrals over isochron surfaces. First, we give the continuous ray-Born integral 21 and then derive the band-limited ray-Born integral 36, which is used in the numerical implementation. In the derivation of the band-limited integral, we use the normal distance from the isochron St to St+Δt/2 to reduce a volume integral to a surface integral. Then, we described a numerical integration over this surface in equation 40, where the integration is carried out over triangles.

The concepts discussed above are derived for the acoustic case, but they can be applied to the elastic case and also to the generalized ray-Born integrals (Chapman, 2004). For an elastic isotropic model, a natural starting point is to carry out the derivations from equations 2136 using the elastic ray-Born integral (e.g., Dahlen et al. [2000], equation 56) as a starting point.

Waveforms computed using this new method are compared with ray-Born integrals computed on a regular grid in the frequency domain. This latter method is accurate as long as sampling intervals for the integration grid are kept small; however, relative to the new time domain method, the frequency-domain method is not as computationally efficient. This is because a summation is carried out over all volume elements in the scattering region and the summation is repeated for every frequency. The numerical efficiency of the time-domain isochron ray-Born method largely depends on the efficiency of the method used to compute the isochrons. If the isochrons are readily available, as they could be in a 1D background model, the integration over triangles is very cheap and the computational cost is insignificant compared with the frequency-domain method.

In the ray-Born integration method by Spencer et al. (1997) and Keers et al. (2002) (see also Chapman, 2004), the computations are made accurate in the time domain by having an integration grid that closely follows isochrons. This is important for modeling at higher frequencies because the isochron surface tends to have irregular curvature. However, compared with the ray-Born integrals using the isochrons explicitly, as discussed in this paper, the integration for the tetrahedron method is carried out over volumes and with a more involved implementation.

It should be noted that isochrons play an important role in many seismic methods because they associate the traveltime of a point in the subsurface with a source-receiver pair. Some examples are Kirchhoff imaging (Bleistein, 1987), the computation of Fresnel volumes (Červený and Soares, 1992), and finite-frequency kernels (Dahlen et al., 2000). Isochrons have also been studied by Iversen (2004), who presents a system of differential equations to compute isochron rays in a similar manner to how conventional rays are computed using the ray equations.

We also compare the new ray-Born modeling method with finite-difference modeling to assess the singly scattered wavefield as an approximation to the complete wavefield, where multiple scattering and multipathing are present. We broadly conclude that, for weak perturbations, or for regions where the wavefield to a large extent consists of singly scattered waves, the ray-Born approximation is the most successful. Based on its ability to produce multiple-free data, it is natural to consider the ray-Born approximation as a modeling method for full-waveform inversion (e.g., during early iterations when the background model is slowly varying) or for the purpose of testing various inversion and imaging algorithms. Furthermore, there are applications in target-oriented studies, for example, in the modeling of the seismic responses caused by small and localized scatterers.

CONCLUSION

We have presented a physically intuitive method for computing 3D ray-Born integrals over isochrons in the time domain. The integrals are simplified because volume integrals are split in a natural way to several surface integrals. Each surface integral gives the waveform in the ray-Born integral for a certain time. The integration technique was used to compute ray-Born seismograms for two velocity models with different scattering characteristics: slowly and rapidly varying velocities. For both model types, the seismograms computed over surface integrals in the time domain are identical to ray-Born seismograms computed over volumes in the frequency domain.

ACKNOWLEDGMENTS

We thank our colleagues at Schlumberger for their constructive comments that led to improvements of the manuscript. In particular, we thank J. Ø. H. Bakke, M. Nickel, and L. Sønneland. The work was funded by the Norwegian Research Council through the Petromaks2 program, project number #245587. This research was partly done while H. Keers was on sabbatical at the University of Kiel and the University of Washington, and support for this work (by the L. Meltzer Høyskolefond, the Faculty of Mathematics and Natural Sciences of the University of Bergen, the Institute of Geoscience at the Christian Albrechts University of Kiel, the Applied Physics Laboratory at the University of Washington, and especially hosts T. Meier [CAU] and B. Odom [UW]) is gratefully acknowledged.

Freely available online through the SEG open-access option.