Pengliang Yang

S U M M A R Y In this paper, we study 3-D multiparameter full waveform inversion (FWI) in viscoelastic media based on the generalized Maxwell/Zener body including arbitrary number of attenu-ation mechanisms. We present a frequency-domain energy analysis to establish the stability condition of a full anisotropic viscoelastic system, according to zero-valued boundary condition and the elastic–viscoelastic correspondence principle: the real-valued stiffness matrix becomes a complex-valued one in Fourier domain when seismic attenuation is taken into account. We develop a least-squares optimization approach to linearly relate the quality factor with the anelastic coefficients by estimating a set of constants which are independent of the spatial coordinates, which supplies an explicit incorporation of the parameter Q in the general viscoelastic wave equation. By introducing the Lagrangian multipliers into the matrix expression of the wave equation with implicit time integration, we build a systematic formulation of multiparameter FWI for full anisotropic viscoelastic wave equation, while the equivalent form of the state and adjoint equation with explicit time integration is available to be resolved efficiently. In particular, this formulation lays the foundation for the inversion of the parameter Q in the time domain with full anisotropic viscoelastic properties. In the 3-D isotropic viscoelastic settings, the anelastic coefficients and the quality factors using bulk and shear moduli parametrization can be related to the counterparts using P and S velocity. Gradients with respect to any other parameter of interest can be found by chain rule. Pioneering numerical validations as well as the real applications of this most generic framework will be carried out to disclose the potential of viscoelastic FWI when adequate high-performance computing resources and the field data are available.

Three-dimensional implementations of reverse time migration (RTM) and full-waveform inversion (FWI) require efficient schemes to access the incident field to apply the imaging condition of RTM or build the gradient of FWI. Wavefield reconstruction by reverse propagation using final snapshot and saved boundaries appears quite efficient but unstable in attenuating media, whereas the checkpointing strategy is a stable alternative at the expense of increased computational cost through repeated forward modeling. We have developed a checkpointing-assisted reverse-forward simulation (CARFS) method in the context of viscoacoustic wave propagation with a generalized Maxwell body. At each backward reconstruction step, the CARFS algorithm makes a smart decision between forward modeling using checkpoints and reverse propagation based on the minimum time-stepping cost and an energy measure. Numerical experiments demonstrated that the CARFS method allows accurate wavefield reconstruction using less timesteppings than optimal checkpointing, even if seismic attenuation is very strong. For RTM and FWI applications involving a huge number of independent sources and/ or applications on architectures with limited memory, CARFS will provide an efficient tool with adequate accuracy in practical implementation.

Many practical seismic applications such as reverse time migration (RTM) and full-waveform inversion (FWI) are usually computation and memory intensive. To perform crosscorrelation in RTM or build the gradient for FWI, it is mandatory to access the forward and adjoint wavefields simultaneously. To do this, there are three methods: One is to read the stored forward wavefield from the disk, the second is using the final snapshot and the stored boundaries via reverse propagation, and the third is remodeling using checkpointing from stored state to another state. Among these techniques, wavefield reconstruction by reverse propagation appears to be a quite straightforward approach ; however, it suffers a stringent memory bottleneck for 3D large-scale imaging applications. The Courant-Friedrichs-Lewy (CFL) condition is a fundamental criterion to determine temporal sampling to achieve stable wavefield extrapolation. The injection of the boundary sequence in time is essentially determined by Nyquist sampling principle, rather than the time interval given by CFL, which is much smaller than the Nyquist requirement. Based on this recognition, we have developed three boundary interpolation techniques, such as the discrete Fourier transform (DFT) interpolation, Kaiser windowed sinc interpolation, and Lagrange polynomial interpolation, for wave-field reconstruction to move from CFL to the Nyquist limit. Wavefield reconstruction via DFT interpolation can be implemented by folding and unfolding steps in the forward simulation and backward reconstruction on the fly. Compared with the DFT interpolation, the wavefield reconstruction methods using Kaiser windowed sinc interpolation and Lagrange polynomial interpolation have better efficiency while remaining a competitive accuracy. These methods allow us to dramatically decimate the boundary without significant loss of information, and they nicely reconstruct the boundary elements in between the samples , making the in-core memory saving of the boundaries feasible in 3D large-scale imaging applications.

The graphics processing unit (GPU) has become a popular device for seismic imaging and inversion due to its superior speed-up performance. We implemented GPU-based full-waveform inversion using the wavefield reconstruction strategy. Because computation on the GPU was much faster than CPU-GPU data communication, in our implementation, the boundaries of the forward modeling were saved on the device to avert the issue of data transfer between the host and device. We adopted the Clayton-Enquist absorbing boundary to maintain the efficiency of the GPU computation. A hybrid nonlinear conjugate gradient algorithm combined with the parallel reduction scheme was used to do computation in GPU blocks. The numerical results confirmed the validity of our implementation.

In this paper, we propose an approximate shrinkage operator in the iterative shrinkage–thresholding (IST) algo-
rithm to interpolate the nonuniformly sampled seismic traces. The key to designing the approximate shrinkage
operator is the use of Taylor series. The redundant Fourier transform is employed to enhance the reconstruction
performance, inspired by the theory of spectral compressive sensing. Numerical experiments using 3D real data
and 5D synthetic data demonstrate the superiority of the proposed method.