Hybrid Kinematic-Dynamic Approach to Seismic Wave-Equation Modeling, Imaging, and Tomography

Estimation of the structure response to seismicmotion is an important part of structural analysis related tomitigation of seismic risk caused by earthquakes. Many methods of computing structure response require knowledge of mechanical properties of the ground which could be derived from near-surface seismic studies. In this paper we address computationally efficient implementation of the wave-equation tomography. This method allows inverting first-arrival seismic waveforms for updating seismic velocity model which can be further used for estimating mechanical properties. We present computationally efficient hybrid kinematic-dynamic method for finite-difference (FD) modeling of the first-arrival seismic waveforms. At every time step the FD computations are performed only in a moving narrowband following the first-arrival wavefront. In terms of computations we get two advantages from this approach: computation speedup and memory savings when storing computed first-arrival waveforms (it is not necessary tomake calculations or store the complete numerical grid). Proposed approach appears to be specifically useful for constructing the so-called sensitivity kernels widely used for tomographic velocity update from seismic data. We then apply the proposed approach for efficient implementation of the wave-equation tomography of the first-arrival seismic waveforms.


Introduction
Mechanical properties of the ground are crucial information for safe construction.In particular, seismic analysis is an important stage of structural analysis in areas with high earthquake risk.This includes estimation of a structure response to seismic motion.Specific topics here include performance-based seismic design [1], dynamic stress concentration in tunnels and underground structures during earthquakes, slope stability analysis, and design of power transmission tower-lines.Many methods of computing structure response require knowledge of mechanical properties of the ground which could be derived from near-surface seismic studies.It is known that subsurface geologic irregularities may have serious effect on actual site response to ground motion which can be modeled for complicated subsurface structures [2].Thus it is important to construct a detailed seismic model from near-surface seismic studies.
In this paper we address computationally efficient implementation of the wave-equation tomography for constructing seismic velocity models.This method requires multiple numerical computation of seismic wavefield propagation (specifically first-arrival waveforms) which is a computationally challenging problem, especially in 3D.
Numerical methods of seismic wavefield modeling (finite-difference (FD), finite-element, etc.) are widely used in seismic data imaging and inversion.They can be used for computing wavefield propagation in complex subsurface structures but usually are computationally expensive.For this reason a lot of work is done for speeding up these methods.
As mentioned by Boore [3] there is no need to compute the wavefield in some region of the computational grid before the first-arrival front passes this region (it is uniformly zero in this case).This approach known as expanding computational domain was reintroduced by Karrenbach [4].Furthermore, if only first-arrival waveforms are of interest, then numerical computations are required only for the grid points in a narrow region following the first-arrival wavefront.Mentioned idea was introduced by [5] and further developed by Kvasnička and Zahradnřandék [6] and Hansen and Jacobsen [7].
The first-arrival front can be computed by numerically solving the eikonal equation on a regular grid using finitedifference method [8,9] or the fast marching (FM) method [10].Note that the cost of computing first-arrival traveltimes is negligible compared to numerically solving the wave equation.
Similar approach to addressing first arrivals is the socalled SWEET algorithm for computing first-arrival traveltimes and amplitudes using the damped wave solution in the Laplace domain as suggested in [11].This algorithm can be implemented for large velocity models by reducing the number of grid points per wavelength with an out-of-core multifrontal algorithm [12].The advantage of this algorithm is that it computes most energetic wave arrivals.However, it is different from the previously mentioned methods because the Laplace damping is distorting the first-arrival waveforms.Thus it provides stable estimate of finite-frequency traveltimes and amplitudes to be used in Kirchhoff type imaging rather than waveforms which can be used in wave-equation migration/inversion type algorithms.
Efficient method of computing the first-arrival waveforms can be beneficial in the reverse-time migration for constructing seismic images [13,14] and the wave-equation tomography (WT) for velocity model building of Luo and Schuster [15] (finite-frequency generalization of the conventional ray-based traveltime tomography).Both methods consider each shot gather independently and imply computation of two wavefields: forward computation of the so-called source field and adjoint (time-reversed) computation of the so-called receiver field.Construction of seismic images or velocity update kernels requires time integration of a product of these two wavefields.Thus there is a need to store the entire history of forward field in the computer memory which is usually unfeasible.Alternatively, one can store the history of the forward field at the boundaries and then recalculate it backward in time simultaneously with the adjoint computation [16].This can reduce memory requirements for storing wavefields but require extra computations to be made.Computation of the source wavefield only in a narrowband around first arrivals seems to be optimal for mentioned applications as it requires considerably less memory.
In this paper we first describe our approach for computing first-arrival waveforms in a narrow computational band following the moving front of the first arrivals.Then we briefly recall the method of the wave-equation tomography and the reverse-time migration which can take advantage of our modeling approach.Finally we present examples showing the speedup that we get in seismic imaging and velocity model update.In this paper we consider only 2D isotropic models to illustrate the concept.

Hybrid Forward Modeling Technique
We further develop this idea of computing waveforms in a running window and apply it in tomographic velocity model building which can be done in different ways.The most straightforward approach is ray tomography which requires only traveltime picks of the first arrivals.
The key step is regrading the computational arrays in the time increasing order that naturally comes out from the fast marching method [10].Such strategy allows setting all computational windows before starting FD calculations.Thus we save computer time needed for the window identification and memory pointer shifts outside the window at each FD time step.Also the resulting windowed wavefield can be stored in RAM at every time step using the proposed strategy.

Hybrid Approach to Forward
Modeling.Let us list main steps of our hybrid modeling approach.For simplicity we consider the 2D wave equation (but it is straightforward to generalize it to the 3D case and elastic wave equation).
The windowed speedup technique proposed by Vidale [5] and illustrated by Kvasnička and Zahradnřandék [6] includes two steps: (i) Solving the eikonal equation to get traveltime at each grid node.(ii) Calculating the solution of the wave equation in the shifted window around the first-arrival front that is defined using the traveltime field from the previous step.
The mathematical background of the windowed calculation technique is obvious: the region of influence of the waveequation solution is determined by the characteristics.
During the second wave-equation-solution step at each time level all grid points can be updated or held unchanged.The FD calculations take place only for updated points that form the shifted zone around the first-arrival front.The solution in unchanged points is taken from the previous time level.The conventional way (Kvasnička and Zahradnřandék [6], Vidale [5]) to determine the updated zone is to check the condition: where  is the number of the time levels, Δ is the FD time spacing,   is the first-arrival traveltime, and  left ,  right are the constants.These constants can be defined, using the dominant wavelength.
The update condition (1) should be checked for every grid point at every time level during the wave-equation FD solution that may cost significant computer time.We would like to show how this checking can be avoided and propose more efficient new strategy for the windowed calculations that is based on traveltime increasing array reordering.

Numerical Grid Upwind Ordering (from Traveltime Solutions)
. As usual, the numerical solvers provide the viscosity solution of eikonal equation, which corresponds to the firstarrival traveltimes.At every node (, ) one has first-arrival traveltime  , .The resulting array can be sorted in the increasing order, so the result would be a vector (time increasing list):

𝑇 (𝑘) , 𝑇(𝑘 + 1) ≥ 𝑇 (𝑘) .
( In principle any eikonal FD solver or other computational techniques can be used to obtain traveltime field  , for the windowed calculations.We are going to illustrate the approach using the multistencils fast marching method, introduced by [17], which is the more recent and robust modification of classic FM. Here we give a very brief description of original fast marching method.The set of all grid points at each step of FM is divided into three sets: accepted points, narrowband, and far points.The main feature of FM algorithm is the special order of traveltimes computation in the grid nodes similar to Dijkstra's algorithm.If point (, ) is accepted, then the traveltime value in it  , becomes frozen and does not change any more.Every point that is not accepted but has an accepted point in its neighborhood should belong to the narrowband; the traveltimes in this set are recomputed.All other points are .At every iteration the point from narrowband with the minimum value of traveltime becomes accepted.For the efficiency narrowband is stored in the heap structure.
Note that if the fast marching solver is used, the rank of the given grid point in the time increasing list ( 2) is the number of acceptances coming out from the FM.
Let us consider the one-to-one mapping: where (, ) is the pair of Cartesian indexes and  is the rank of the  , in the traveltime increasing list (2);  = 1, . . .,   ;  = 1, . . .,   ;  = 1, . . .,   ×   .We use mapping (3) to reorder all the arrays, used for the wave-equation FD modeling.For example, let us consider velocity model matrix (, ); it can be stored in the form of a vector ().Note that mapping (3) should be remembered, as soon as there is a need for getting back to the Cartesian grid.For instance, one can use two integer vectors: () = , () = .Actually almost the same mapping (3) happens within the fast marching method for the narrowband, stored in the time-sorted heap structure.

Computing and Storing Widowed
Wavefield.The proposed time increasing reordering is natural for the windowed speedup calculations.Each time window is set up by its beginning number  begin and the end number  end that are defined using list (2).Thus the wave-equation FD calculations at every time level  are done in the cycle from  =  begin () to  =  end () and there is no need for checking the "update" condition (1).
In principle any time-domain wave-equation solver can be modified for windowed calculations, based on the introduced time increasing reordering.The space FD stencil for the node (, ) includes some neighborhood points, at least (± 1,  ± 1) (if one is using the second-order central differences).As soon as all the arrays are reordered, one should implement mapping (3) several times to get the ranks  for all the stencil nodes.The more efficient strategy is to introduce additional arrays for the ranks of neighborhood points.For instance, the rank of the "left" node (−1, ) for the node (, ) with the rank  can be remembered in the array  left ().
Note that one needs to store some additional indexationmapping arrays for the window and stencils neighborhoods.As an alternative, the first-arrival wavefront propagation should be done simultaneously with wave-equation solution just ahead of the computational window.The FM approach provides the possibility of such "on the fly" traveltime computations.However one should keep mapping (3) (at least backward).
The proposed windowed techniques, based on the time increasing array reordering, provide the opportunity of storing the "history" of the wavefield.At every time level  = 1, . . .,  the wavefield can be recorded just one by one at the current time level  from  =  begin () to  =  end () in the one large array, containing the entire "history" of the windowed wavefield.After the computations, one can get the access to the wavefield history by reading the wavefield in time level cycle, using the stored widowed limiting ranks  begin () and  end () and using mapping (3) for getting to the Cartesian grid.Note that one can produce the wavefield reading in a reversed time.

Fast Marching Method.
Consider the eikonal equation in isotropic medium: where (x) is the traveltime in the point x = (, ) and (x) is the seismic wave velocity.We will give here a brief description of original fast marching method in the two-dimensional case.
The set of all grid points at each step of FM is divided into three sets: accepted points, narrowband, and far points.The main feature of FM algorithm is the special order of traveltimes computation in the grid nodes similar to Dijkstra's algorithm.If point (, ) is accepted, then the traveltime value in it  , becomes frozen and does not change any more.Every point that is not accepted but has an accepted point in its neighborhood should belong to the narrowband; the traveltimes in this set are recomputed.All other points are far.At every iteration the point from narrowband with the minimum value of traveltime becomes accepted.For the efficiency narrowband is stored in the heap structure.

The Virieux Staggered
Grid FD Scheme.In the present paper 2D acoustic wave equation is considered: where (x) is pressure, (x) is the density, and (x) is the wave speed.
We use a staggered grid [18] to solve the wave equation (7).In order to implement this scheme the acoustic equation should be rewritten as two first-order equations: where (V  , V  ) is the particle velocity vector.The considered FD scheme uses the staggered grid, so that V  is attached to the points (, ) (  =  0 + ℎ,   =  0 + ℎ) of a Cartesian grid, V  is attached to the "half-integer" nodes ( + 1/2,  + 1/2), and pressure  is attached to the points ( + 1/2, ).Also the velocity and pressure are attached to different time layers.The staggered grid FD stencil for acoustic wave equation is shown in Figure 2. The FD equations are

Application to Wave-Equation Imaging and Tomography
We show the benefits of the proposed hybrid calculation-andstoring windowed technique for speeding up the reverse-time migration [13,14] and the wave-equation tomography [15].For both methods we consider shot gathers independently and for each of them we compute two wavefields: the socalled source field   (x, ) is a result of forward computation and the so-called adjoint field   (x, ) is a result of the reversetime computation; that is, time-reversed data is plugged in as source signature in the receiver positions.Interaction of these two wavefields is producing seismic image or the velocity model update.
In the reverse-time migration for each shot gather (corresponding to x  ) we get an image as a time integration of a product of these two wavefields: In the wave-equation tomography for each shot gather (corresponding to x  ) we also construct sensitivity kernels as an interaction of two wavefields: Here the adjoint wavefield p (x, ) is not just timereversed data but time-reversed pseudo-residual [15]: where  obs (x  , ) are observed data at receivers x  ,  (x  ) = ∫  obs (x  ,  + Δ)   (x  , ) , (13) and Δ are traveltime misfits estimated by cross-correlating observed data  obs (x  , ) and synthetic data   (x  , ) directly from waveforms (maximum of the cross-correlation): Δ (x  ) = arg min  ∫  obs (x  ,  + )   (x  , ) d.( 14) While summing up the results for all shot gathers (for different x  ), we get final migrated image: Similarly, we sum up all shot gather sensitivity kernels (for different x  ) in order to get the velocity model update: Iterative methods can be used to invert for the velocity model (x) in several steps: where γ(x)  is the sensitivity kernel at th iteration and   is updating step length.
The key step of constructing migrated image (10) or velocity model update (11) is to compute the source field   (x, ) and store it for further cross-correlation with the adjoint field   (x, ).For example, precomputed   (x, ) is required for estimating the pseudo-residual used for initiating the adjoint field (cf.(12)).
Thus there is a need to store the entire history of forward source field   (x, ) in the computer memory which is usually unfeasible.Alternatively, one can store the history of the forward field at the boundaries and then recalculate it backward in time simultaneously with the adjoint computation.This can reduce memory requirements for storing wavefields but require one extra waveform computation to be made.Computation of the source wavefield only in a narrowband around the first arrivals seems to be optimal for mentioned applications as it requires considerably less memory.Further we plan to apply the proposed approach only to the source field   (x, ) while computing the adjoint receiver field   (x, ) in a standard mode (note that it does not need to be stored in memory).Ricker wavelet with 20 Hz central frequency was used as a source function.Computed pressure data for source 1 are shown in Figure 4: "full" wavefield computed by standard FD scheme in panel (a) and "windowed" computation by the proposed hybrid kinematic-dynamic method in panel (b).Note that the first-arrival waveforms look the same.One can see differences for late arrivals where refracted waves can be seen on "full" data.For source 2 computed "full" and "windowed" data look pretty much the same because there are no refracted waves in this case (for this reason we do not show corresponding gathers).

Examples
For this example our hybrid computational approach is 20 times faster than the standard approach (full wavefield computation).Also it takes 17 times less memory to store the source wavefield within narrowband following first arrival (compared to storing full wavefield at each time step).These speedup and memory economy will increase with increasing velocity model size (in terms of dominant wavelength).

Wave-Equation
Tomography.Let us consider synthetic example motivated by near-surface applications.For this example we use gradient model with high velocity square anomaly as shown in Figure 5(a).For this true velocity we computed synthetic data using 60 Hz source wavelet (sources and receivers are placed every 5 m at the top of the model).For these data we performed velocity model update using wave-equation tomography (initial model is a model with linear gradient).Source wavefield was computed using our hybrid "windowed" approach.The result of velocity inversion ( 16) is shown in Figure 5(b).One can see that the velocity anomaly is recovered.
Next example is motivated by applications in seismology.We consider true gradient velocity model with two anomalies as shown in Figure 6(a).Source positions are shown by black stars and receivers were located regularly at the surface.Synthetic data were computed for 10 Hz Ricker wavelet.We performed velocity model update using waveequation tomography (initial model is a model with linear gradient).Source wavefield was computed using our hybrid "windowed" approach.The result of the velocity inversion is shown in Figure 6(b).One can see that the velocity anomalies can be recovered.We use our "windowed" hybrid approach to modeling the source field   (x, ); "windowed" modeling is 33 times faster than "full" waveform modeling (1.5 and 50 seconds  correspondingly on Intel Core i7-4770 3.4 GHz CPU).It also needs 20 times less memory to store the first-arrival wavefield, that is, 0.8 Gb versus 17 Gb.Reverse-time migrated images for 50 shot gathers (15) are shown in Figure 8: using "full" waveform computation (a) and using our "windowed" computation (b).Both migration results look identical.Note that the receiver wavefield   (x, ) was computed in a "full" waveform mode in both cases.

Discussion
Examples show that our "windowed" approach allows storing "windowed" source wavefield in memory for any 2D model while fitting "full" source wavefield in memory is not feasible in most cases.Overall speedup of using our "windowed" approach to wave-equation tomography and reverse-time migration is close to 3. For standard approach the source field does not fit memory and thus one requires three standard "full" wavefield computations [16] for implementing tomography update or migration.Our "windowed" approach dramatically reduces time required for the source field computation and it can be stored in memory for further use.Thus computational cost of our "windowed" approach is close to one "full" wavefield computation of the receiver field.
We want to mention potential problems in applying our "windowed" modeling approach for contrast models [6].The eikonal solver provides first arrivals which may correspond to weak head waves while energetic waves useful for imaging may arrive later.Computational window following the firstarrival front may miss these waves.However, the class of velocity models in which our approach works is large enough to make it useful for practical applications.

Conclusions
We have presented a new hybrid kinematic-dynamic finitedifference (FD) method for computing and storing the firstarrival seismic waveforms.In this "windowed" computational approach at every time step we perform FD computations only in a narrow region following the first-arrival wavefront.This wavefront is precomputed using fast marching eikonal solver.We have proposed a new effective strategy for computations and wavefield storing which is based on time increasing reordering.
We showed several examples of successful use of our "windowed" modeling approach in wave-equation tomography and migration.Proposed approach showed 20-30 times speedup in computing and 17-20 times less memory for storing the source wavefield when compared to standard "full" wavefield computing and storing.It allows overall 3 times speedup for implementing wave-equation tomography or reverse-time migration.
In this paper we consider only 2D isotropic models for illustrating the concept.Extension to 3D models is straightforward.Generalization of the approach to anisotropic models requires efficient anisotropic eikonal solver algorithm availability.It is a topic of our future research.

4. 1 .
Forward Modeling.Let us present an example of windowed wavefield computations.The smooth Marmousi velocity model is shown in Figure 3. Receivers are shown by white triangles and source locations are shown by stars.

4. 3 .
Reverse-Time Migration.Here we consider the reversetime migration example.Synthetic data was computed for Marmousi model shown in Figure 7(b); 50 shot gathers for 30 Hz Ricker sources were placed every 60 m.Migration velocity model is shown in Figure 7(a).

Figure 7 :
Figure 7: Smooth migration model (a) and true velocity model (b).

Figure 8 :
Figure 8: Reverse-time migrated images computing source field by standard full wavefield (a) and our "windowed" approach (b).