SMIwiz-2.0: Extended functionalities for wavefield decomposition, linearized and nonlinear inversion

We extend the functionalities of SMIwiz open source software to include up-down wavefield separation, reflection waveform inversion, as well as linearized waveform inversion in data and image domain. The fundamental functionalities for 2D/3D wave modelling and imaging (reverse time migration and nonlinear full waveform inversion) are backward compatible with improvements in seismic imaging processing. Reproducible examples are supplied to verify these developments.

New version program summary

Program Title: SMIwiz

CPC Library link to program files: https://doi.org/10.17632/tygszns27k.2

Developer's repository link: https://github.com/yangpl/SMIwiz

Licensing provisions: GNU General Public License v3.0

Programming language: C, Shell, Fortran

Software dependencies: MPI [1], FFTW [2]

Journal reference of previous version: Comput. Phys. Commun. 295 (2024) 109011. https://doi.org/10.1016/j.cpc.2023.109011

Does the new version supersede the previous version?: Yes.

Nature of problem: Seismic modelling and imaging (linearized and nonlinear waveform inversion).

Solution method: Conjugate gradient (CGNR) method for linearized inversion, quasi-Newton LBFGS and line search for nonlinear optimization.

Summary of revisions: The following new features (specified by a mode parameter) are added to extend the functionalities of SMIwiz:

• 1.
  Up-going and down-going wavefield separation (mode=10)
  We have implemented decomposition of up-going and down-going wavefield following [3,4]. The method relies on the use of Hilbert transform for constructing analytic wavefield. To avoid huge storage requirement and expensive computation, the Hilbert transform was applied to the wavelet or source time function. By simultaneously propagating the wavefield excited by a complex-valued analytic source time function (its imaginary part is the Hilbert transform of the real part), the analytic wavefield (with real and imaginary parts stored separately) is then obtained by time stepping algorithms. The up-going wavefield and down-going wavefields are then computed by a filtering process in the frequency-wavenumber domain (based on the sign of the frequency and the wavenumber) [3,4].
• 2.
  Improved reverse time migration (RTM, mode=2)
  We have improved RTM by switching from cross-correlation imaging condition to an impedance kernel. This change was motivated by the fact that classic cross-correlation imaging condition suffers from the low-frequency noises. These low-frequency noises are created along the wave path, which are unwanted image artefacts [5]. Removing these low-frequency artefacts normally requires a Laplace filtering, which is mathematically equivalent to using impedance kernel [6]. Switching to the imaging condition with impedance kernel directly gives an artefact-free image without Laplace filtering.
• 3.
  Reflection waveform inversion (RWI, mode=1 & 4, rwi=1)
  Classical full waveform inversion mainly updates the model using direct arrivals and far-offset refractions, while near-offset reflections are not fully explored. RWI is therefore proposed since the concept of migration-based traveltime tomography [7]. The construction of a sensitivity kernel for background model in [8] is in the shape of a pair of rabbit ear clearly demonstrate the capability of RWI in updating the deep part of the velocity model, which goes very much beyond the depth that first arrivals and diving waves can reach. Our implementation inherits the first-order velocity-pressure formulation of acoustic wave equation embedded in SMIwiz [9], and employs the velocity-impedance parametrization [10]. Compared to [10], the distinctive feature of our implementation is the flexibility to switch between simultaneous multiparameter inversion and alternating optimization over velocity and impedance parameters using variable projection method [11].
• 4.
  Multiparameter linearized waveform inversion (mode=3)
  Linearized waveform inversion [12], often referred to as least-squares reverse time migration (LSRTM), intends to retrieve the perturbation of the model parameters, assuming the background of these parameters is known. Using density ρ and bulk modulus κ (linked with P wavespeed $V_p$ via $\kappa =\rho V_p^2$) as the basic parameter family, we implemented LSRTM in the data domain, allowing the linearized inversion under different families of parametrization with family 1 $(V_p,\rho )$ or family 2 $(V_p,I_p)$ (where the impedance is $I_p=\rho V_p$). Switching parametrization from $(\rho ,\kappa )$ to $(V_p,\rho )$ and $(V_p,I_p)$ is trivial thanks to the chain rule. To manipulate these relations conveniently, a log scaling has been applied to convert multiplication into addition to facilitate efficient computation. The CGNR method [13, Algorithm 9.7] is used to solve the problem. It requires repetitive Born modelling and migration at each iteration to form a Hessian in Krylov space, which is rather computationally expensive.
• 5.
  Migration deconvolution using point spread function (mode=7 & 8)
  Assuming the Hessian H is a priori known or accessible. The above data-domain LSRTM for finding reflectivity images can be reformulated into a new linear optimization problem, often called migration deconvolution,(1)$\underset{\delta m}{\min }{\Vert H\delta m-m_{rtm}\Vert }^2,$ where $m_{rtm}$ is the RTM image by migrating the seismic reflections into image domain. For a model of size N, the Hessian H is normally a very large and dense matrix consisting of $N^2$ elements. Fully storing the Hessian is therefore impractical. We have adopted the concept of point spread function (PSF) [14,15] to sparsely sample the Hessian using evenly distributed point scatters (spacings are specified by control parameters n1win, n2win and n3win with mode=7 and mdopt=1 being activated). The CGNR algorithm is then utilized to optimize the problem in equation (1) (mode=8). Each iteration requires only the matrix-vector product, which is constructed on-the-fly using trilinear interpolation [16], assuming the Hessian is approximated by a sparse banded matrix. The interpolated Hessian therefore loses symmetric positive definite property, requiring the evaluation of the product of both Hessian H and its transpose $H^T$ with an input vector (cf. the subroutine matmul_Hv()). Our implementation has symmetrized the interpolated Hessian, i.e. using $(H+H^T)/2$.
• 6.
  Migration deconvolution using non-stationary Wiener filter (mode=7 & 9)
  The problem above may be solved in a non-iterative fashion, assuming the Hessian is locally invariant in space. This corresponds to a circulant Toeplitz structure of Hessian matrix which may be diagonalized in the Fourier domain. Our implementation is thus performed patch by patch with domain overlapping. In each patch, we use Wiener filtering (frequency-domain division) to obtain a local solution for model perturbations. Finally, they are stacked together following the principle of partition of unity [17]. In the first step, the effect of Hessian in mode=7 is captured using RTM image as the input (corresponding to mdopt=2). At the second stage, migration deconvolution using non-stationary deblurring filtering is carried out using mode=9.

References

• [1]
  https://www.mpich.org/.
• [2]
  http://fftw.org/.
• [3]
  F. Liu, G. Zhang, S.A. Morton, J.P. Leveille, An effective imaging condition for reverse-time migration using wavefield decomposition, Geophysics 76 (1) (2011) S29–S39.
• [4]
  T.W. Fei, Y. Luo, J. Yang, H. Liu, F. Qin, Removing false images in reverse time migration: the concept of de-primary, Geophysics 80 (6) (2015) S237–S244.
• [5]
  Y. Zhang, J. Sun, Practical issues in reverse time migration: true amplitude gathers, noise removal and harmonic source encoding, First Break 27 (2009) 53–59.
• [6]
  H. Douma, D. Yingst, I. Vasconcelos, J. Tromp, On the connection between artifact filtering in reverse-time migration and adjoint tomography, Geophysics 75 (6) (2010) S219–S223.
• [7]
  G. Chavent, F. Clément, S. Gòmez, Automatic determination of velocities via migration-based traveltime waveform inversion: a synthetic data example, SEG Tech. Program Expand. Abstr. 1994 (1994), https://doi.org/10.1190/1.1822731.
• [8]
  S. Xu, D. Wang, F. Chen, G. Lambaré, Y. Zhang, Inversion on reflected seismic wave, SEG Tech. Program Expand. Abstr. 2012 (2012) 1–7, https://doi.org/10.1190/segam2012-1473.1.
• [9]
  P. Yang, SMIwiz: an integrated toolbox for multidimensional seismic modelling and imaging, Comput. Phys. Commun. 295 (2024) 109011, https://doi.org/10.1016/j.cpc.2023.109011.
• [10]
  W. Zhou, R. Brossier, S. Operto, J. Virieux, Full waveform inversion of diving & reflected waves for velocity model building with impedance inversion based on scale separation, Geophys. J. Int. 202 (3) (2015) 1535–1554.
• [11]
  G. Golub, V. Pereyra, Separable nonlinear least squares: the variable projection method and its applications, Inverse Probl. 19 (2) (2003) R1.
• [12]
  A. Tarantola, Linearized inversion of seismic reflection data, Geophys. Prospect. 32 (1984) 998–1015.
• [13]
  Y. Saad, Iterative Methods for Sparse Linear Systems, SIAM, Philadelphia, 2003.
• [14]
  I. Lecomte, Resolution and illumination analyses in PSDM: a ray-based approach, Lead. Edge 27 (5) (2008) 650–663, https://doi.org/10.1190/1.2919584, https://pubs.geoscienceworld.org/tle/article-pdf/27/5/650/3086199/650.pdf.
• [15]
  R.P. Fletcher, D. Nichols, R. Bloor, R.T. Coates, Least-squares migration—data domain versus image domain using point spread functions, Lead. Edge 35 (2) (2016) 157–162.
• [16]
  L.N. Osorio, B. Pereira-Dias, A. Bulcão, L. Landau, Migration deconvolution using domain decomposition, Geophysics 86 (3) (2021) S247–S256.
• [17]
  P. Yang, Z. Ji, A comparative study of data- and image-domain LSRTM under velocity-impedance parametrization, Comput. Geosci. (2025), submitted for publication.