Least Squares Migration: 2 of 2 Articles

The previous article in this two part series ("Least Squares Migration: 1 of 2 Articles") discussed important similarities and distinctions between RTM (Reverse Time Migration), LS-RTM (Least Squares RTM), and FWI (Full Waveform Inversion). It was noted that traditionally, LSM does not update the velocity model.

The discussion here describes an elegant process to recover high-resolution subsurface models of reflectivity, velocity and relative density from marine shot gathers with minimal pre-processing applied. The description here avoids mathematics and highlights the key breakthroughs that have enabled conventional FWI and data domain (LS-RTM to evolve into a methodology that can simultaneously recover several models of subsurface properties without significant leakage between parameters, make no assumptions about rock physics relationships, and bypass traditionally cascaded processing and imaging flows.

Introduction

Full Waveform Inversion (FWI) is a topic published and posed in many forms, many of them highly mathematical and therefore alien to a large proportion of practicing geoscientists. Below I attempt to frame FWI into more practical terms and explain why recent developments in seismic ‘imaging by inversion’ are exciting for geological interpreters. The most recent development takes marine seismic shot gathers with minimal processing and simultaneously produces both a velocity model and reflectivity image in vertical units of depth.

As discussed below, these results can also be used to recover a high-resolution resolution density model in one elegant step. Collectively, these three outputs present a highly efficient alternative to traditional cascaded processing and imaging workflows to produce migrated seismic stack images. The enabling technology is an evolution of FWI and data domain Least-Squares Reverse Time Migration (LS-RTM) algorithms.

Full Waveform Inversion: Building an Earth Model Directly from Observed Seismic Data

FWI typically begins with a very simplistic velocity model and compares synthetic shot gathers to observed data. Differences or ‘perturbations’ in modeled vs. observed seismic data are analyzed on a trace-by-trace basis and attributed to localized errors in the model. Several iterations of updates to the model as ‘pixels’ of differing scale are pursued until the modeled data is acceptably close to the observed data and the model is assumed to represent the true subsurface property distribution.

FWI is commonly implemented as an adaption of seismic imaging—and shares the modeling and migration engine of RTM. Several attempts have been made in recent years to adapt FWI to build interpretable models of the Earth, effectively bypassing traditional cascaded processing and imaging workflows with a highly abbreviated (but computationally rather expensive) methodology. Unfortunately, FWI is famously prone to instability and the generation of strong artifacts.

Platforms to Make FWI Stable and Robust

In a nutshell, FWI is more likely to work if an accurate ‘background’ velocity model represented by low wavenumbers can be accurately built before the high wavenumber details are recovered.

Note that ‘wavenumber’ is the reciprocal of wavelength and represents the ‘spatial frequency’ of the seismic wavefield. Just as broadband seismic imaging seeks to recover high-bandwidth temporal frequencies for each seismic event, FWI should obtain all spatial frequencies in the velocity spectrum except those that cannot be resolved due to ‘incompleteness’ in the data—limited source-receiver offsets and band-limited source frequencies.

The background model (low wavenumber values) can be robustly obtained during acquisition by decreasing the source frequency or by having access to long offset transmitted waves, i.e., diving waves.

Ideally, FWI can use both transmitted and scattered waves in a ‘full wavefield’ solution as reflections can provide much deeper model updates than diving waves allow for each offset. 

There are, however, several challenges to the pursuit of a true ‘full wavefield’ FWI:

• Diving waves are often only partly recorded for conventional streamer lengths of 8km. The pursuit of very long offset diving wave information has correspondingly driven interest in ocean bottom seismic (OBS) acquisition for regions affected by salt-related imaging challenges such as the Gulf of Mexico (GoM).
• Diving wave FWI requires offsets 3-6 times the maximum depth of interest for which model updates are desired. Model updates from diving wave information are therefore typically only shallow.
• Reflections can provide much deeper model updates than diving waves for conventional streamer lengths, but it is difficult to model reflections from a smooth initial model unless a density model corresponding to key impedance layer contrasts is used.
• FWI is highly sensitive to arrival time misalignments between modeled and observed data more than half a cycle (known as ‘cycle skipping’) and is prone to converging to local minima (incorrect solutions).
• High wavenumber model artifacts are easily generated when reflectivity due to density changes are incorrectly mapped as velocity updates.

Elegant solutions are now available to the challenges above.

Overall, the earth model can be represented as a smoothly varying (low-wavenumber) velocity macro-model onto which are superimposed sharp contrasts in acoustic properties (high-wavenumbers) associated with geological boundaries.

(upper row) Sensitivity kernels for full wavefield FWI. The left version shows the ‘banana’ and ‘rabbit ears’ kernels that represent the tomographic terms for diving wave and reflection information, respectively, and these recover the low wavenumber ‘background model’. The right version is the traditional full wavefield FWI version that also contains the migration isochron that provides high-wavenumber updates from reflections. (lower row) FWI model derived from low-wavenumber updates with transmission and reflection data (left); and the FWI model contaminated with high-wavenumber updates associated with the migration isochron (right). The ‘high resolution’ content of the model on the right represents high-wavenumber reflectivity erroneously mapping as velocity model updates. Courtesy of Øystein Korsmo (PGS).

The top row of the figure above shows where model updates can be expected for a given offset and reflector depth. Known by practitioners as ‘sensitivity kernels’, these plots also illustrate how the high-wavenumber migration component of FWI can easily contaminate the low-wavenumber tomography components. The modelled reflector depth is 4km, the source-receiver offset is 5km, and the background velocity is 4 km/s.

There are three features in the upper row:

1. The ‘banana’ velocity kernel corresponding to low wavenumber diving wave contributions (both panels);
2. The ‘rabbit ears’ velocity kernel corresponding to low wavenumber reflection contributions (both panels); and 
3. The ‘migration isochron’ reflectivity kernel corresponding to high wavenumber reflection contributions (upper-right panel only).

Both the banana and rabbit ears kernels represent tomographic terms within FWI, and their width represents the first Fresnel zone for the relevant frequency and velocity.

As the migration term is usually an order of magnitude larger in amplitude than the tomography terms, the high wavenumber components of the model are updated faster than the low wavenumbers.

If the starting velocity model is not accurate, or other incorrect assumptions mentioned below are made, the high-wavenumber content will be constructed in the wrong spatial locations, and it is very difficult to remove this content through later iterations.

Another challenge to the use of reflection information in FWI is the difficulty of generating reflections from a smooth initial velocity model without a knowledge of the density model. Formulae such as Gardner’s relation are commonly used to predict a density model, but such assumptions easily introduce high-wavenumber artifacts such as shown in the lower-right panel above.

This problem is also inherent in recent industry pursuits of so-called ‘FWI Imaging’ wherein spatial derivatives of a post-inversion FWI velocity model is used with a density model assumption to derive a pseudo-reflectivity model.

The approach will leak reflectivity into the velocity model when using a conventional imaging condition, i.e., the reflectivity kernel shown in the upper-right panel above contributes velocity model updates.

Even when a velocity model is employed, it only acts as a passive parameter during FWI, and all reflectivity updates are simply mapped as velocity updates. The lower-right panel above shows how such an approach can introduce an incorrect velocity perturbation when the reflectivity is driven by the density contrast and not the velocity.

PGS developed an elegant solution several years ago that allows the migration isochrone contribution to be separated from the tomographic terms in the manner of the upper-left panel above.

Furthermore, a reformulation of the variable density wave-equation that is often used for modeling wave propagation in FWI also helps avoid spurious high-wavenumber updates from reflectivity.

Shot gather modelled with the first order Born approximation (left); Modeled shot gather using the vector reflectivity implementation of FWI (right). Note the additional direct arrival and diving wave events on the right.

The ‘vector reflectivity’ formulation allows FWI to accurately model high-wavenumber reflections with a smooth initial velocity model, as well as modelling refractions and diving waves, which is beyond the capability of the currently available industry methods.

The right panel in the figure above shows modeled shot gathers with accurate amplitudes for all depths. No approximations are made within the modeling engine, in contrast to traditional ‘Born modeling’.

The left panel shows the Born modeling common to data domain LS-RTM. The nature of the approximations inherent within Born modeling result in amplitudes restricted to near-to-mid offsets/angles (left side).

In contrast, the full wavefield vector reflectivity modeling engine generates head waves, diving waves and refractions in addition to both low and high angle reflections (right side). Anisotropy is easily incorporated if necessary.

Collectively, the two aforementioned FWI improvements, combined with other features such as ‘dynamic time warping’ to avoid cycle skipping problems, allowed PGS to recast FWI into a robust solution that recovers accurate low-wavenumber velocity model updates from the full seismic wavefield for large depths with conventional streamer lengths.

Extending FWI to Simultaneous Multi-Parameter Inversion

The PGS simultaneous inversion solution uses diving waves and reflections to recover accurate low-wavenumber velocity model updates with the tomography terms of FWI (upper row). The method simultaneously uses the migration term within a non-linear least-squares framework to recover the high-wavenumber reflectivity model (lower row). The methodology avoids leakage between the model parameters being inverted. The vector reflectivity formulation of the acoustic wave equation also enables relative density to be accurately and robustly recovered from the two model outputs (refer to the following figure).

Depth slice from a velocity model (left) and the computed relative density model derived from both the velocity and reflectivity models (right). Courtesy of Yang Yang, PGS.

PGS extended the FWI methodology in 2021 into a multi-parameter inversion for both velocity and reflectivity models based upon a robust scale separation. An efficient non-linear implementation of LS-RTM that can be simultaneously run at different wavenumber scales within an FWI framework:

• An accurate low-wavenumber model is derived without leakage from spurious high-wavenumber artifacts (the upper row of the previous figure). High-waveumber updates are independently recovered in a simultaneous manner (the lower row of the previous figure).
• The input to PGS Ultima is the same as for traditional FWI: common shot gathers with minimal pre-processing and a starting velocity model. An initial reflectivity is computed during the first iteration as an RTM image.
• The high-wavenumber updates to the LS-RTM reflectivity model benefit the low-wavenumber velocity model updates without introducing high-wavenumber velocity model artifacts, and the final velocity model that contributes to the final reflectivity model should be better than that derived from traditional standalone FWI.
• Note that traditional LS-RTM does not update the velocity model, is only able to utilize near-to-mid reflection angles within the inversion updates and is therefore susceptible to errors in velocity when updating reflectivity.
• The velocity and reflectivity model deliverables are produced in ‘one step’ rather than a traditional multi-month cascaded processing and imaging workflow.

Recovering High Resolution Density Models

The final reflectivity model should implicitly contain a valid representation of the subsurface spatial density, so the reflectivity and velocity models can be used to directly recover high-resolution relative density models.

The results i the final figure above are an example from the Orphan Basin, offshore Newfoundland and Labrador, Canada. Multisensor GeoStreamer acquisition was used with 8.1km streamer length.

A traditional pre-stack simultaneous AVA inversion workflow to recover an estimate of density is as follows:

• Careful signal processing, multiple and noise removal, velocity model building, pre-stack migration and gather conditioning in a cascaded workflow. Attention must be paid to gather flatness, AVA (amplitude vs. angle) preservation and phase stability in a frequency-dependent manner. As pre-stack time migration (PSTM) is the common method, geostatistical time-to-depth correction on the gathers is important.
• An accurate low-frequency model (LFM) based on valid depth-dependent rock physics relationships for the study area is critical for the next step to be quantitatively accurate.
• Pre-stack simultaneous inversion of the condition gathers and LFM is used to estimate elastic impedances (Ip and Is).
• Density can also be estimated via a 3-term AVA equation such as Shuey or Aki & Richards. Note that both equations are a linear approximation of the Zoeppritz equations and assume weak contrasts in both velocity and density—neither of which apply to the methodology described here.

In contrast, the simultaneous inversion workflow described here requires no special gather conditioning or rock physics knowledge, is entirely data-driven, and incorporates multiples, diving waves and refractions into the final estimation. A remarkable result!

Summary

By reformulating the theory behind FWI, accurate low-wavenumber model updates can be derived for large depths with conventional streamer lengths and without contamination from spurious high-wavenumber artifacts.

The ‘vector reflectivity’ formulation also enables an efficient implementation of non-linear LS-RTM that can be simultaneously run at different wavenumber scales within an FWI framework.

The two outputs of the simultaneous inversion solution described are an accurate velocity model and a reflectivity model with enhanced spatial resolution and amplitude balance.

The velocity and reflectivity model deliverables are produced in ‘one step’ rather than a traditional multi-month cascaded processing and imaging workflow. The spatially varying reflectivity and velocity models can be used to directly recover high-resolution density models, thereby representing an exciting opportunity for subsurface characterization.

The new workflow requires no special gather conditioning or rock physics knowledge, is entirely data-driven, and incorporates multiples, diving waves and refractions into the final estimation.

Acknowledgements

Special thanks to Nizar Chemingui, Jaime Ramos-Martinez, Yang Yang, Sriram Arasanipalai, Neil Paddy, Øystein Korsmo, Yermek Balabekov and Cyrille Reiser (all of PGS). I also thank PGS, TGS, and OilCo NL.for access to the Blomidon MC3D data.

Interested in Hearing More?

• Canada | Appraising Hydrocarbon Potential in Orphan Basin: 30-minute Webinar.

Disclaimer

The content discussed here represents the opinion of Andrew Long only and may not be indicative of the opinions of Petroleum Geophysical AS or its affiliates ("PGS") or any other entity. Furthermore, material presented here is subject to copyright by Andrew Long, PGS, or other owners (with permission), and no content shall be used anywhere else without explicit permission. The content of this website is for general information purposes only and should not be used for making any business, technical or other decisions.