Mathematical Problems in Engineering

Volume 2015 (2015), Article ID 543540, 8 pages

http://dx.doi.org/10.1155/2015/543540

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

^{1}Institute of Petroleum Geology and Geophysics SB RAS, Pr. Ac. Koptyuga 3, Novosibirsk 630090, Russia^{2}Novosibirsk State University, Street Pirogova 2, Novosibirsk 630090, Russia

Received 15 April 2015; Revised 1 July 2015; Accepted 6 July 2015

Academic Editor: Yuming Qin

Copyright © 2015 Alexandr S. Serdyukov and Anton A. Duchkov. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

Estimation of the structure response to seismic motion is an important part of structural analysis related to mitigation 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 to make 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.

#### 1. 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 finite-difference 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 so-called 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.

#### 2. 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.

##### 2.1. 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 wave-equation 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 , 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*.

##### 2.2. Numerical Grid Upwind Ordering (from Traveltime Solutions)

As usual, the numerical solvers provide the* viscosity* solution of eikonal equation, which corresponds to the first-arrival 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):

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); ; ; .

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.

##### 2.3. 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 and the end number that are defined using list (2). Thus the wave-equation FD calculations at every time level are done in the cycle from to 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 (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 for the node with the rank can be remembered in the array .

Note that one needs to store some additional indexation-mapping 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 the wavefield can be recorded just one by one at the current time level from to 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 and and using mapping (3) for getting to the Cartesian grid. Note that one can produce the wavefield reading in a reversed time.

##### 2.4. Fast Marching Method

Consider the eikonal equation in isotropic medium: where is the traveltime in the point and is the seismic wave velocity. We will give here a brief description of original* fast marching method* in the two- dimensional case. Let us consider the rectangular Cartesian grid: , , , and . The original FM method [10] is based on the first-order “upwind” FD approximation of eikonal equation (4): where The essence of “upwind” approximation (5) is shown in Figure 1. The grid nodes for FD stencil are choosing principle that is based on traveltime increasing “causality”; that is, the direction of “wind” is the direction of wavefront expansion.