Dong-Joo Min

Waveform inversion using a logarithmic wavefield

Although waveform inversion has been studied extensively since its beginning 20 years ago, applications to seismic field data have been limited, and most of those applications have been for global-seismology- or engineering-seismology-scale problems, not for exploration-scale data. As an alternative to classical waveform inversion, we propose the use of a new, objective function constructed by taking the logarithm of wavefields, allowing consideration of three types of objective function, namely, amplitude only, phase only, or both. In our waveform inversion, we estimate the source signature as well as the velocity structure by including functions of amplitudes and phases of the source signature in the objective function. We compute the steepest-descent directions by using a matrix formalism derived from a frequency-domain, finite-element/finite-difference modeling technique. Our numerical algorithms are similar to those of reverse-time migration and waveform inversion based on the adjoint state of the wave equation. In order to demonstrate the practical applicability of our algorithm, we use a synthetic data set from the Marmousi model and seismic data collected from the Korean continental shelf. For noisefree synthetic data, the velocity structure produced by our inversion algorithm is closer to the true velocity structure than that obtained with conventional waveform inversion. When random noise is added, the inverted velocity model is also close to the true Marmousi model, but when frequencies below 5 Hz are removed from the data, the velocity structure is not as good as those for the noise-free and noisy data. For field data, we compare the time-domain synthetic seismograms generated for the velocity model inverted by our algorithm with real seismograms and find that the results show that our inversion algorithm reveals short-period features of the subsurface. Although we use wrapped phases in our examples, we still obtain reasonable results. We expect that if we were to use correctly unwrapped phases in the inversion algorithm, we would obtain better results.

Improved frequency‐domain elastic wave modeling using weighted‐averaging difference operators

We develop a new finite-difference scheme that reduces the number of grid points per wavelength required in frequency-domain elastic modeling. Our approach computes weighted averages of the spatial second-order derivative and the mass acceleration terms using a 25-point computational stencil. By determining the weighting coefficients to minimize numerical dispersion and numerical anisotropy, we reduce the number of grid points to 3.3 per shear wavelength, with a resulting error in velocities smaller than 1%. Our choice of grid points reduces the computer memory needed to store the complex impedance matrix to 4% of that for a conventional second-order scheme and to 54% of that for a combined second-order scheme. The 25-point weighted averaging scheme of this paper makes it possible to accurately simulate realistic models. Numerical examples show that this technique can achieve the same accurate solutions with fewer grid points than those from previous frequency-domain second-order schemes. Our technique can be extended directly to 3-D elastic modeling; the computational efficiency will be even greater than that realized for 2-D models.

Frequency-Domain Elastic Full Waveform Inversion Using the New Pseudo-Hessian Matrix: Experience of Elastic Marmousi-2 Synthetic Data

Abstract A proper scaling method allows us to find better solutions in waveform inversion, and it can also provide better images in true-amplitude migration methods based on a least-squares method. For scaling the gradient of a misfit function, we define a new pseudo-Hessian matrix by combining the conventional pseudo-Hessian matrix with amplitude fields. Because the conventional pseudo-Hessian matrix is assumed to neglect the zero-lag autocorrelation terms of impulse responses in the approximate Hessian matrix of the Gauss–Newton method, it has certain limitations in scaling the gradient of a misfit function relative to the approximate Hessian matrix. To overcome these limitations, we introduce amplitude fields to the conventional pseudo-Hessian matrix, and the new pseudo-Hessian matrix is applied to the frequency-domain elastic full waveform inversion. This waveform inversion algorithm follows the conventional procedures of waveform inversion using the backpropagation algorithm. A conjugate-gradient method is employed to derive an optimized search direction, and a backpropagation algorithm is used to calculate the gradient of the misfit function. The source wavelet is also estimated simultaneously with elastic parameters. The new pseudo-Hessian matrix can be calculated without the extra computational costs required by the conventional pseudo-Hessian matrix, because the amplitude fields can be readily extracted from forward modeling. Synthetic experiments show that the new pseudo-Hessian matrix provides better results than the conventional pseudo-Hessian matrix, and thus, we believe that the new pseudo-Hessian matrix is an alternative to the approximate Hessian matrix of the Gauss–Newton method in waveform inversion.

Traveltime and amplitude calculations using the damped wave solution

Because of its computational efficiency, prestack Kirchhoff depth migration remains the method of choice for all but the most complicated geological depth structures. Further improvement in computational speed and amplitude estimation will allow us to use such technology more routinely and generate better images. To this end, we developed a new, accurate, and economical algorithm to calculate first-arrival traveltimes and amplitudes for an arbitrarily complex earth model. Our method is based on numerical solutions of the wave equation obtained by using well-established finite-difference or finite-element modeling algorithms in the Laplace domain, where a damping term is naturally incorporated in the wave equation. We show that solving the strongly damped wave equation is equivalent to solving the eikonal and transport equations simultaneously at a fixed reference frequency, which properly accounts for caustics and other problems encountered in ray theory. Using our algorithm, we can easily calculate first-arrival traveltimes for given models. We present numerical examples for 2-D acoustic models having irregular topography and complex geological structure using a finite-element modeling code.

Weighted-Averaging Finite-Element Method for 2D Elastic Wave Equations in the Frequency Domain

We present a weighted-averaging frequency-domain finite-element method for an accurate and efficient 2D elastic wave modeling technique. Our method introduces three kinds of supplementary element sets in addition to a basic element set that is used in the standard finite-element method. By constructing global stiffness and mass matrices for four kinds of element sets and then averaging them with weighting coefficients, we obtain a new global stiffness and mass matrix. With optimal weighting coefficients determined by a Marquardt–Levenberg method to minimize grid dispersion and grid anisotropy, we can reduce the number of nodal points per shear wavelength from 33.3 (using the standard finite-element method) and 20 (using the eclectic method) to 5, with the errors of group velocities no larger than 1%. By reducing the number of grid points per wavelength, we achieve a 97% and 75% reduction of computer memory required to store the complex impedance matrix for a band-type matrix solver and a nested dissection method, respectively, compared with those of the eclectic method. Our method gives approximate solutions compatible with exact solutions for an infinite homogeneous, a semi-infinite homogeneous (Lamb9s problem), and a horizontal two-layer model with fewer grid points than the standard and the eclectic method. A major advantage of the weighted-averaging finite-element method for the elastic wave equation is that it provides solutions very close to correct solutions for Lamb9s problem economically, unlike most of the displacement approaches. In addition, our scheme makes the complex impedance matrix symmetric, which satisfies reciprocity. Seismic forward modeling techniques that satisfy reciprocity are of critical importance in seismic imaging and inversion because we can economically calculate a Jacobian matrix using the reciprocity. Successful simulation of a large-size model shows that our method can be used for the simulation of wave propagation in the geological model needed in the reverse-time migration or seismic inversion.

Refraction tomography using a waveform-inversion back-propagation technique

One of the applications of refraction-traveltime tomography is to provide an initial model for waveform inversion and Kirchhoff prestack migration. For such applications, we need a refraction-traveltime tomography method that is robust for complicated and high-velocity-contrast models. Of the many refraction-traveltime tomography methods available, we believe wave-based algorithms to be best suited for dealing with complicated models. We developed a new wave-based, refraction-tomography algorithm using a damped wave equation and a waveform-inversion back-propagation technique. The imaginary part of a complex angular frequency, which is generally introduced in frequency-domain wave modeling, acts as a damping factor. By choosing an optimal damping factor from the numerical-dispersion relation, we can suppress the wavetrains following the first arrival. The objective function of our algorithm consists of residuals between the respective phases of first arrivals in field data and in forward-modeled data. The...

Free-Surface Boundary Condition in Finite-Difference Elastic Wave Modeling

We design a new frequency-domain, finite-difference approach, based on a displacement formulation, which correctly describes the stress-free conditions at a free surface. In the conventional, displacement-based finite-difference method, we assign both displacements and material properties such as density and Lame constants to nodal points (a node-based grid set), whereas in our new finite-difference method, displacements are still defined at nodal points but material properties within cells (a cell-based grid set). In our new finite-difference technique using the cell-based grid set, stress-free conditions at the free surface are described by the changes in the material properties without any additional stress-free boundary condition. By conducting accuracy tests, we confirmed that the new second-order finite differences formulated in the cell-based grid set generate numerical solutions compatible with analytic solutions within the range of errors determined by dispersion analysis; the new, cell-based, weighted-averaging finite-difference method also yields better solutions than the old, node-based, weighted-averaging finite-difference method. The cell-based finite-difference method does not violate the reciprocity.

Traveltime calculations from frequency‐domain downward‐continuation algorithms

We present a new, fast 3D traveltime calculation algorithm that employs existing frequency‐domain wave‐equation downward‐continuation software. By modifying such software to solve for a few complex (rather than real) frequencies, we are able to calculate not only the first arrival and the approximately most energetic traveltimes at each depth point but also their corresponding amplitudes. We compute traveltimes by either taking the logarithm of displacements obtained by the one‐way wave equation at a frequency or calculating derivatives of displacements numerically. Amplitudes are estimated from absolute value of the displacement at a frequency.By using the one‐way downgoing wave equation, we also circumvent generating traveltimes corresponding to near‐surface upcoming head waves not often needed in migration. We compare the traveltimes computed by our algorithm with those obtained by picking the most energetic arrivals from finite‐difference solutions of the one‐way wave equation, and show that our trave...

Frequency-domain reverse-time migration with source estimation

Although artificially generated seismic sources such as dynamite, vibroseis, and air guns are used in seismic exploration, it is not easy to exactly recover the source wavelet in field recording or in data processing. For this reason, seismic data processing often assumes that an explosive-source wavelet can be described by a well-known function (e.g., a Ricker wavelet), a near-offset trace, or a deconvolved wavelet. In frequency-domain waveform inversion, it has been proven that a source wavelet can be estimated by an optimization method, and incorporating the source wavelet estimation into an inversion algorithm yields better inversion results. We have developed source wavelet estimation into 2D two-way frequency-domain reverse-time migration. The source wavelet is first estimated independently of reverse-time migration by an optimization method such as the full Newton method. It is then used in reverse-time migration. This source-wavelet-incorporated reverse-time migration algorithm is applied to three...

3D Laplace-domain full waveform inversion using a single GPU card

Abstract The Laplace-domain full waveform inversion is an efficient long-wavelength velocity estimation method for seismic datasets lacking low-frequency components. However, to invert a 3D velocity model, a large cluster of CPU cores have commonly been required to overcome the extremely long computing time caused by a large impedance matrix and a number of source positions. In this study, a workstation with a single GPU card (NVIDIA GTX 580) is successfully used for the 3D Laplace-domain full waveform inversion rather than a large cluster of CPU cores. To exploit a GPU for our inversion algorithm, the routine for the iterative matrix solver is ported to the CUDA programming language for forward and backward modeling parts with minimized modification of the remaining parts, which were originally written in Fortran 90. Using a uniformly structured grid set, nonzero values in the sparse impedance matrix can be arranged according to certain rules, which efficiently parallelize the preconditioned conjugate gradient method for a number of threads contained in the GPU card. We perform a numerical experiment to verify the accuracy of a floating point operation performed by a GPU to calculate the Laplace-domain wavefield. We also measure the efficiencies of the original CPU and modified GPU programs using a cluster of CPU cores and a workstation with a GPU card, respectively. Through the analysis, the parallelized inversion code for a GPU achieves the speedup of 14.7 – 24.6 x compared to a CPU-based serial code depending on the degrees of freedom of the impedance matrix. Finally, the practicality of the proposed algorithm is examined by inverting a 3D long-wavelength velocity model using wide azimuth real datasets in 3.7 days.