The wavefield is typically simulated for seismic exploration applications through solving the wave equation for a specific seismic source location. The direct relation between the form (or shape) of the wavefield and the source location can provide insights useful for velocity estimation and interpolation. As a result, I derive partial differential equations that relate changes in the wavefield shape to perturbations in the source location, especially along the Earth's surface. These partial differential equations have the same structure as the wave equation with a source function that depends on the background (original source) wavefield. The similarity in form implies that we can use familiar numerical methods to solve the perturbation equations, including finite difference and downward continuation. In fact, we can use the same Green's function to solve the wave equation and its source perturbations by simply incorporating source functions derived from the background field. The solutions of the perturbation equations represent the coefficients of a Taylor's series type expansion of the wavefield as a function of source location. As a result, we can speed up the wavefield calculation as we approximate the wavefield shape for sources in the vicinity of the original source. The new formula introduces changes to the background wavefield only in the presence of lateral velocity variation or in general terms velocity variations in the perturbation direction. The approach is demonstrated on the smoothed Marmousi model.
The wave equation is the central ingredient in simulating and constraining wave propagation in a given medium. No other formulation, including the eikonal equation or ray tracing, can be as conclusive and elaborate (includes most traveltime and amplitude aspects) as does the full elastic wave equation (Aki & Richards 1980). Thus, the wave equation has the near complete far-field description of the wave behaviour for a given velocity function. To reduce the computational cost, P-wave propagation in the Earth's subsurface is approximately simulated by numerically solving the acoustic wave equation. The wavefield in acoustic media is described by a scalar quantity. Kinematically, for P waves in the far field of the source, the acoustic and elastic wave equations are similar; they both yield the same eikonal equation in isotropic media.
The first-order dependence of wavefields, and specifically the acoustic wave equation, on media parameters is described by the Born approximation. It is a single scattering approximation, used in seismic applications to approximate the perturbed wavefield due to a small perturbation of the reference medium. It is also, in its inverse form, used to help infer medium parameters from observed wavefields (Cohen & Bleistein 1977; Panning et al. 2009). In the spirit of the Born approximation, the partial differential equations, introduced here, describe the wavefield shape first- and second-order dependence on a source location perturbation. Data evolution as a function of changes in acquisition parameters goes back to the development of normal moveout and the transformation of common-offset seismic gathers from one constant offset to another (Bolondi et al. 1982). Bagaini & Spagnolini (1996) identified offset continuation (OC) with a whole family of pre-stack continuation operators, such as shot continuation (Bagaini & Spagnolini 1993), dip moveout as a continuation to zero offset (Hale 1991; Alkhalifah 1996), and 3-D azimuth moveout (Biondi et al. 1998). Even residual operators between different medium parameters or acquisition configurations are described by Alkhalifah & Biondi (2004) and Alkhalifah (2005). All these methods are based on a geometrical optics development using constant velocity approximations.
Alkhalifah & Fomel (2010) suggested a plane wave source perturbation expansion for the eikonal equation. Their approach rendered linearized forms of the eikonal equation capable of predicting the traveltime field for a shift in the source location represented by a shift in the velocity field in the opposite direction. The development here follows the same approach applied now to the wave equation. The major result in the eikonal application is the linearized forms of the new perturbation equations extracted from the conventionally nonlinear eikonal partial differential equation.
The source perturbation introduced here is based on a plane wave expansion around the source. For homogeneous or vertically inhomogeneous media, the wavefield calculated for a source in any location along the horizontal surface is valid in all other locations, and as a result no modifications are needed. For laterally inhomogeneous media, this statement is no longer true and the difference in the wavefield depends on the complexity of the lateral velocity variation. In this paper, I develop partial differential equations that approximately predict such changes, and thus, can be used to simulate wavefields for sources at other positions in the vicinity of the original source. Fig. 1 illustrates the concept of wavefield evolution as a function of source perturbation in which the wavefield evolves as a function of source location perturbation l granted that the velocity field changes laterally. In Fig. 1, the wavefield wave front schematic reaction depicts a general velocity decrease with l. In this paper, we consider source perturbation in the lateral direction, which adheres to our familiarity with surface seismic exploration applications. However, equivalent perturbations may be applied in the depth direction as well, as Alkhalifah & Fomel (2010) showed for the eikonal equation, which may have applications in datuming.
In this section, I derive partial differential equations that relate changes in the wavefield shape to perturbations in the source location. We will start by taking derivatives of the wave equation with respect to lateral perturbation and then use the Taylor's series expansion to predict the wavefield form at another source.eq. (1) given by f(x, y, z, t), defined usually at a point, or represented by the wavefield u(x, y, z, t) around time t= 0 as an initial condition. A change in the source location along the surface, while keeping its source function stationary, is equivalently represented, in the far field, by shifting the velocity field laterally by the same amount in the opposite direction and thus can be represented by the following wave equation form: Fig. 2 shows the relation of a single source and image point combination with the velocity shift. To evaluate the wavefield response to lateral perturbations, we take the derivative of eq. (2) with respect to l, where the wavefield is dependent on the source location as well [u(x, y, z, l, t)], which yields eq. (3), and setting l= 0, the location in which we evaluate the equation for the Taylor's series expansion, yields eq. (3) again with respect to l, which yields
Again, by substituting , as well as, Dx into eq. (6), and setting l= 0, the second-order perturbation equation is given by
The Green's Function
The Green's function represents the response and behaviour of the wavefield if the source was a point impulse given theoretically by the Dirac delta function. The function allows us to solve the wave equation using an integral formulation as we convolve the Green's function with a source function. The most complicated part of this type of solution of the wave equation is the construction of the Green's function. Kirchhoff modeling and migration is a special case of this type of integral solution and provides incredible speed upgrades over finite difference implementations with some loss in quality, because of the limitations in the calculation of the phase and amplitude components of the Green's function.eq. (1) with a source function f(x0, t0) is given by eq. (4) eq. (4) has the same form as the wave equation with the same velocity function, eq. (7) with a slightly more complicated source function given by
As we have seen earlier, the equations associated with the wavefield perturbation has a form similar to the wave equation. As a result, any of the methods typically used to solve the wave equation suffices for the perturbation partial differential equations. The most general and straightforward of these methods is the finite-difference approach. With proper space and time grid distribution this method provides acceptable solutions regardless of the complexity of the velocity model. Its only limitation is the relative slow execution speed.
To apply the source perturbation, I first solve the original wave equation for a particular point source as a background field for the perturbation step. In the process, we store the Laplacian evaluation as they are needed for the perturbation calculation. Using an initial condition of Dx(x, z, t= 0) = 0, we can solve for Dx using eq. (4) and add the solution to the original wavefield using eq. (8), and thus, obtain a new wavefield shape approximating that for another source location. For higher-order accuracy, we also solve for Dxx using eq. (7) with a similar initial condition, Dxx(x, z, t= 0)= 0, and include it in the Taylor's series expansion with terms from the solutions of the original background source and Dx.
To avoid potential problems with the storage requirement especially in 3-D, we can solve the wave equation and the corresponding perturbation equations simultaneously, and thus, use the already evaluated derivatives (Laplacian) directly. In this case, the cost of the perturbation finite difference application is similar to the cost of solving the wave equation. Thus, the cost of the first- and second-order expansions to obtain the wavefield for other sources is equivalent to two and three times, respectively, of the cost of solving the wave equation for a single source. However, the information obtained approximates the wavefield for infinite source possibilities in the vicinity of the original source.
Because the perturbation equations have the same form as the wave equation they adhere to the same Courant-Friedrichs–Lewy (CFL) condition (Courant et al. 1928). Thus, the time step is constrained by the grid spacing, which in turn relies on velocity, based on the following formula:
I test the developed partial differential equations on two examples: a simple lens velocity model and the complex Marmousi model. In both cases, I compare the wavefields obtained from a direct finite difference solution of the wave equation for a particular source to that obtained by perturbing the solution from a nearby source using the first- and second-order approximations. In the Marmousi case, I also compare the resulting surface recorded synthetic data for all options.
Because the differential equation depends on velocity changes in the direction of the perturbation, we test the methodology on a model that contains a lens anomaly in an otherwise homogeneous medium with a velocity of 2 km s−1 (Fig. 3). The lens apex is located at 600 m laterally and 500 m in depth with a velocity perturbation of +250 m s−1 (or 12.5 per cent). The lens has a diameter of 300 m. Using this model we will test the accuracy of the first- and second-order perturbation equations.
For a source located at the surface 0.3 km from the origin, I apply a second-order in time and fourth-order in space finite difference approximation to the wave equation as well as the perturbation equations to simulate a source 50 m away from the original source position along the surface. A separate direct finite difference calculation using the wave equation is done for a source at 0.35 km location for comparison. Fig. 4 shows a snap shot at time 0.5 s of the wavefield generated for the source at 0.35 km (left), as well as the snap shot at the same time for the perturbed wavefields using the first-order approximation (middle) and the second-order approximation (right). All three wavefields look similar.
However, if we subtract the actual wavefield for the 0.3 km source from that of the 0.35 km one after we superpose the sources we obtain the difference between the wavefields. This difference occurs only if there is lateral velocity variation. Because there are no lateral variations in the velocity field in Fig. 3 prior to wave front from the source crossing the lens we expect that the difference snapshot plots are zero. However, at time 0.3 s the difference, where the wave front starts to interact with the lens as shown in Fig. 5, appears as expected largest for the unperturbed case (left), while the differences for the perturbed case are much smaller, especially in the case of the second-order expansion. All three plots in Fig. 5 are displayed using the same range, for easy comparison, and this range is maintained for all Figures in this section.
Fig. 6 shows a snap shot at 0.5 s (at the same time as wavefields shown in Fig. 4) of the difference. Again, the second-order approximation shows less difference, and thus, a better match than the first-order approximation and definitely the unperturbed wavefield. In fact, in the unperturbed wavefield a clear polarity reversal at the anomaly apex is evident.
Clearly, the perturbation formulas help reduce the difference between the actual wavefield and the perturbed one. An even closer look suggests that most of the difference is amplitude related.
The Marmousi model
The geometry of the Marmousi is based, somewhat, on a profile through the Cuanza basin Versteeg (1993). The target zone is a reservoir located at a depth of about 2500 m. The model contains many reflectors, steep dips and strong velocity variations in both the lateral and the vertical directions (with a minimum velocity of 1500 m s−1 and a maximum of 5500 m s−1). However, the Marmousi model includes complex discontinuities that pose problems to the perturbation formulation. As a result, we smooth the velocity model to obtain the model in Fig. 7 (right). The point source considered here is a Ricker wavelet with a 15 Hz peak frequency.
Again using the fourth-order finite difference approximation in space and second order in time we solve the wave equation for a source located at the surface at lateral position 5 km, A snap shot of the resulting wavefield at time 1.2 s is shown in Fig. 8 (left). Solving the wave equation for a source located 25 m away results in the snap shot of the wavefield at 1.2 s shown in Fig. 8 (middle). Superimposing the sources for the two fields and subtracting them yields the difference shown in Fig. 8 (right). All three snap shots are plotted at the same scale (and this scale is maintained for all Figures in this section) and thus the difference, which is totally due to lateral inhomogeneity, is relatively large. It is especially large for the parts of the wave front that were exposed to large lateral variations in the smoothed Marmousi model. This difference represents the wavefield we anticipate from the solution of new perturbation equation.
Using the new perturbation partial differential equations, I predict this difference from the original wavefield with source at location 5 km. We then add this difference to that original wavefield using eq. (8), which provides an approximate to the wavefield for a source located at 5.025 km. Fig. 9 shows a 1.2 s snap shot of the wavefield computed directly from a source at 5.025 km (left) and that obtained from the first-order perturbation expansion (middle), as well as, the difference between the two wavefields (right). Clearly, the difference is now less than that in Fig. 8 in which perturbation was not used.
In fact, if we display the differences side by side along with that associated with the second-order perturbation approximation, Fig. 10 demonstrates that the difference decreases considerably for the higher-order perturbation approximation, shown on the right.
One of the main objectives of solving the wave equation is to simulate the behaviour of the wavefield at the surface (the measurement plane) as a function of time. Fig. 11 shows the difference between the 5.025 km source and the 5 km source common-shot gathers after superimposing the sources (left) and compares it with difference between the 5.025 km source and first-order (middle) and the second-order (right) perturbed versions. Clearly, the source gather extracted from using the perturbation equations better resemble the directly evaluated one than the source gather that does not include the perturbation. Specifically, most of the primary reflections in the section are seemingly well modelled by the perturbation approximation, as evidence by the small difference between the directly extracted gather and the perturbed one.
The transformation of the wavefield's shape as a function of source location provides valuable information for many applications, and in particular interpolation, velocity estimation and imaging. All three applications rely in one way or another on the relation between the wavefields and the source location. To predict the content of a missing shot (or receiver) gather, we usually rely on reflection slope information of nearby common shot gathers (Fomel 2002) to extend the information to the missing locations. For homogeneous and vertically inhomogeneous media, the process of interpolation is trivial as the common shot gathers are the same regardless of surface position. The complication occurs when the velocity varies laterally and differences, as we have seen above, can be large. Using these perturbation partial differential equations we can estimate the changes needed to fill in the gaps. This can be done as part of a finite-difference modeling or a reverse time migration process. It also can be done using a point source to generate the wavefield in a forward finite-difference approach or using a boundary condition, as typically the case for reverse time extrapolation of the receiver wavefields recorded at the surface for the purpose of imaging. The equations shown above have no source restrictions and their development are not based on a particular source.
A major drawback of using conventional methods to solve the wave equation is that typically the velocity information and complexity have no baring on the efficiency of obtaining such solutions. In the development here, the perturbed wavefields are only excited by lateral velocity variation and in the absence of such variations we do not need to evaluate the perturbations. This allows us to implement selective computations that depends on the wavefield complexity and isolate areas of contribution based on the velocity field.
Nevertheless, the accuracy of the first- and second-order expansion approximations introduced here depends on the size of the source (or velocity) shift. Unlike, the traveltime version (Alkhalifah & Fomel 2010) the wavefield is highly oscillatory (sinusoidal components)and thus their Taylor's series approximation accuracy is dependent on the wavelength of the perturbed wavefield within the context of the lateral velocity complexity. The accuracy here is synonymous to what we encounter using the Born approximation. However, unlike Born approximation, the source functions in eqs (4) and (7) depend on the lateral velocity variation, not the source perturbation. Specifically, if the lateral velocity change induces perturbations in the wavefield that exceeds a half wavelength, we will encounter aliasing in the construction. This issue effects, more frequently, large dips and large perturbations with respect to the wavelength. However, unlike conventional source or velocity perturbation developments, this wavefield shape perturbation approach is far more stable and explicitly depends on the complexity of the lateral velocity variation as for the case of lateral homogeneity the approach is exact independent of amount of perturbation. Fig. 12 shows a 1.2 s snap shot of the differences between the wavefields obtained directly from a source and that obtained from nearby sources and perturbed a distance of 50 m (left), 75 m (middle) and 100 m (right) for the smoothed Marmousi model using the second-order expansion. As expected the difference (error) is larger for the bigger perturbation. Also, we can observe that the dipping parts of the wavefield have larger errors as the effective change is bigger. Of course, we have to remind our selves that we are dealing with the Marmousi model, which is highly complex, and we can expect better results for smoother models. Also, we can observe that the difference is mainly manifested in the amplitude, where the kinematics (phase) show little difference.
The transformation in the wavefield shape as a function of source location is directly related to the lateral velocity variation. Such transformation is described by partial differential equations that have forms similar to the conventional acoustic wave equation in which their solutions provide the coefficients needed for a Taylor's series type expansion. The source function, for the perturbation equation, depends on the background wavefield of the original source as well as lateral derivatives of the velocity of the medium. As a result, while the second order expansion, which requires solving two PDEs, provide the best approximation of the perturbation in generally smooth velocity models, as expected, and similar to the Born approximation, the accuracy of the approximation here reduces with the size of the source perturbation. However, unlike the Born approximation, the accuracy here depends only on the amount of lateral velocity variation, not on the velocity perturbation acting as a secondary source.
I thank Sergey Fomel for many useful discussions. I thank KAUST for its support. I also thank the editors and the anonymous reviewers for their critical and helpful review of the paper.