-
PDF
- Split View
-
Views
-
Cite
Cite
Siamak Abolhassani, Leo Hoogerbrugge, Dirk Jacob Verschuur, One-way reflection waveform inversion with depth-dependent gradient pre-conditioning, Geophysical Journal International, Volume 240, Issue 1, January 2025, Pages 652–672, https://doi.org/10.1093/gji/ggae397
- Share Icon Share
SUMMARY
Reflection waveform inversion (RWI) is a technique that uses pure reflection data to estimate subsurface background velocity, relying on evolving seismic images. Conventional RWI operates in a cyclic workflow, with two key components in each cycle—migration and reflection tomography. Conventional RWI may result in suboptimal background velocity estimation, partly due to limited or unresolved resolution within each component in each cycle. While gradient pre-conditioning with the reciprocal of Hessian information helps resolve this issue in both components of RWI, it becomes impractical for a large number of model parameters. One-way reflection waveform inversion (ORWI) is a reflection waveform inversion technique in which the forward modelling scheme operates in one direction (downward and then upward) via virtual parallel depth levels within the medium. Leveraging the ORWI framework, we decompose and reduce the linear Hessian operator (also known as the approximate Hessian or Gauss–Newton Hessian) into multiple smaller suboperators. In particular, the diagonal blocks of the monofrequency approximate Hessian operators, each corresponding to a single depth level within the medium, are extracted and inverted to pre-condition the corresponding monofrequency gradients in both the migration and reflection tomography components of ORWI. This depth-dependent gradient pre-conditioning transforms standard ORWI into a high-resolution, yet computationally feasible version aimed at addressing suboptimal velocity estimation, referred to as high-resolution ORWI. The effectiveness of the proposed approach is demonstrated through successful applications to synthetic data examples.
1 INTRODUCTION
Tarantola (1984) introduced the conventional full waveform inversion (FWI) in acoustic approximation, a well-known non-linear data-fitting technique to estimate high-resolution subsurface velocity models. Conventional FWI has become widely adopted for evaluating shallow subsurface targets. Indeed, incorporating refracted and diving waves, it offers a detailed subsurface model and has proven to be highly effective for this purpose (e.g. Vigh et al. 2011; Zou et al. 2014). However, when attempting to map the targets that lie beyond the reach of refracted and diving waves (“deep targets”), its performance may falter (e.g. Plessix et al. 2013; Brittan & Jones 2019). This limitation arises because, while the sampling of high-reflective model wavenumbers is nonlinearly dependent on the sampling of low-propagative model wavenumbers in FWI (Audebert & Ortiz-Rubio 2018), conventional FWI fails to actively sample low-propagative model wavenumbers within the deep targets because of weak transmission-after-reflection wave paths within the FWI gradient, commonly known as “rabbit-ear” wave paths. Scaled by reflection coefficients, the rabbit-ear wave paths exhibit weaker amplitudes compared to other wave paths within the FWI gradient, that is the transmission wave path and migration isochrone (Fig. 1). This failure to sample low-propagative model wavenumbers significantly limits the capture of high-reflective model wavenumbers by conventional FWI, highlighting the significance of model scale separation for reflection tomography (Mora 1989).

The wave paths constructing the FWI gradient: the transmission wave path, labelled as number 1, the pair of transmission-after-reflection wave paths (rabbit-ear wave paths), labelled as number 2, and the migration isochrone, labelled as number 3. As can be observed, the pair of transmission-after-reflection wave paths appears in a lighter shade, illustrating its much weaker amplitude compared to the other existing wave paths. This occurs as a consequence of scaling by the reflection coefficient of the reflector within the medium. Hence, it can be concluded that the FWI gradient is dominated by the transmission wave path and migration isochrone, while the pair of transmission-after-reflection wave paths remains relatively inactive.
Over the years, several variants of migration-based velocity analysis (MVA), each built upon the principle of model scale separation based on Born modelling, have been developed. Among them, the ones incorporating wave-equation forward modelling have stood out due to the improved handling of wave propagation in complex media. Wave-equation MVA (WEMVA) algorithms (Sava & Biondi 2004) typically optimize an image-domain error function measured on the common image gathers (CIGs) either via a classical semblance or a differential semblance (DS) function, measuring the flatness (coherency property) in CIGs. While with the classical semblance function, the flatness measure of CIGs is given by the stack power of image amplitudes over surface offsets (Chauris & Noble 2001), the flatness measure of CIGs with the DS function is given by squaring the L2 norm of the derivative of image amplitudes with respect to surface offsets (Shen et al. 2003). WEMVA may also be formulated in the subsurface offset domain, where, at the correct velocity, the image is expected to be focused at zero offset. In this domain, measuring the focusing of subsurface-offset CIGs mirrors the use of the classical semblance and DS function with surface CIGs (Shen 2005; Shen & Symes 2008; Lameloise et al. 2015). Alternatively, a scattering angle can serve as the extension parameter, with flatness measured on angle-domain CIGs (Biondi & Symes 2004).
Expanding on the foundations of extended-domain WEMVA, Symes (2008) formulated a theory that laid the groundwork for integrating extended-domain WEMVA and FWI, leveraging their synergy. This theory was subsequently validated by Almomin & Biondi (2012). In recent initiatives to enhance MVA, Chauris & Cocher (2017) introduced the notion of inversion velocity analysis (IVA) by substituting standard migration with iterative migration to map cleaner CIGs, reducing migration artifacts and multiple crosstalk. Moreover, they adopted a pseudo-inverse of the Born modelling operator, rather than an adjoint, to reduce the migration iterations. In the same vein, Assis et al. (2022) introduced tomographic Hessian-based IVA.
Employing migrated images for velocity model building also triggered another trend of development known as reflection waveform inversion (RWI). RWI represents a wave-equation reflection tomography tool with a data-domain error function (Xu et al. 2012; Zhou et al. 2012; Berkhout 2012; Wang et al. 2013). Conventional RWI in acoustic approximation alternately solves a two-parameter minimization problem for the background velocity and “reflectivity” of the subsurface. Although RWI and the migration-based traveltime waveform inversion (MBTT; Clément et al. 2001) share some similarities, it is important to note that the migration step in MBTT operates on reflectivities in time rather than depth—images in time are largely invariant to velocity, unlike images in depth, but this assumes a 1-D velocity model.
Unlike FWI, RWI is specifically tailored to sample deep targets tomographically harnessing exclusively reflection data through the pair of transmission-after-reflection wave paths. RWI involves estimating the background velocity of the subsurface in a cyclic process as displayed in Fig. 2, relying on an evolving stacked image that fits near-offset observed and modelled “phase information” within each cycle. Therefore, it is expected that RWI naturally suppresses any ambiguities initiated by the phase mismatch in the near-offset reflection data (Provenzano et al. 2023), particularly within the deep areas or the areas with challenging data acquisition conditions where FWI suffers severely. To improve the RWI tomogram, specifically within reach of refractions and diving waves, some studies have also investigated the impact of integrating early-arrival FWI and RWI in order to steer the search direction toward recovering a full (horizontal and vertical) low local-model-wavenumber spectrum for the background velocity model (Zhou et al. 2015; Wu & Alkhalifah 2015). For a full review of the other aspects of RWI in recent years, the reader can refer to Yao et al. (2020).

Conventional RWI cycle in which background velocity estimation and image reconstruction alternate.
Although RWI is built upon the notion of model scale separation, making this separation has been a struggle. To address the issue, several approaches have been put forth. Among others, we reference Alkhalifah (2014) for introducing the scattering-angle filtering solution, Wang et al. (2013) for adopting the wavefield decomposition solution, and Xu et al. (2012) for leveraging Born modelling as the solution. Additionally, Berkhout (2012, 2014) introduced the joint migration inversion (JMI) technique that makes use of a one-way forward modelling scheme called full-wavefield modelling (FWMod), parametrized based on reflectivity and acoustic background velocity, to fulfill the model scale separation by wavefield decomposition into upgoing and downgoing wavefields. FWMod models full (both primaries and multiples) downgoing and upgoing reflection wavefields within the JMI framework, in which zero-lag cross-correlations of same-direction source and receiver wavefields (either both upwards or both downwards) reconstruct the reflectivity model update (small-scale medium variations), and correlations of opposite-traveling source and receiver wavefields (one upgoing, one downgoing) estimate the background velocity model update (large-scale medium variations).
Adopting the concept of JMI for reflection waveform inversion and constraining FWMod to primary wavefield modelling (PWMod), assuming internal multiples are weak or resolved via pre-processing, this paper represents and investigates upon a one-way reflection waveform inversion technique, hereafter referred to as standard ORWI. The approporiate parameterization in PWMod leads to a natural scale separation—separation of migration isochrones and transmission-after-reflection wavepaths—in standard ORWI. This separation allows for the independent calculation of migration and tomography gradients within each ORWI cycle, thereby freeing tomograms from high-reflective model wavenumbers. Algorithm 1 shows the cyclic workflow of standard ORWI, which involves an alternating sequence between the image and background velocity reconstructions. Standard ORWI cyclically minimizes the error function
where
Suboptimal background velocity estimation in RWI can be attributed, among others, to the use of seismic images with limited or unresolved resolution as well as suboptimal preservation of true amplitudes in each cycle (e.g. Hou & Symes 2015; Chauris & Cocher 2017; Gomes & Yang 2017). With the same line of reasoning, one can also link suboptimal background velocity estimation in RWI to the use of tomograms with limited or unresolved resolution. Pratt et al. (1998) showed that the Newton optimization method in waveform inversion estimates accurate velocity models with higher resolution than the gradient decent method since it uses the reciprocal of second-order derivatives of the error function (inverse Hessian) to pre-condition the gradient. However, the calculation demand of the Hessian matrix-operator for large-scale problems in seismology renders the Newton method impossibly expensive in computing terms. Several attempts, either in the data or image domain, have been made thus far to approximate an effective yet cost-effective pre-conditioner: migration deconvolution filters (Hu et al. 2001; Yu et al. 2006), matching filters (Guitton 2004; Aoki & Schuster 2009; Guitton 2017; Guo & Wang 2020; Yang et al. 2021), point spread functions (Lecomte 2008; Fletcher et al. 2016; Yang et al. 2022), target-oriented solutions (Valenciano et al. 2006; Tang 2009), asymptotic pseudo inverses to the Born modelling operator (ten Kroode 2014; Hou & Symes 2015; Chauris & Cocher 2017), linear Hessian (also known as the Gauss–Newton Hessian or approximate Hessian) (Pratt et al. 1998), pseudo-Hessian (Choi et al. 2008; Jun et al. 2015), diagonal approximate Hessian (Chavent et al. 1999; Shin et al. 2001; Plessix & Mulder 2004), and approximating the application of the Hessian inverse operator on the search direction (Brossier et al. 2009; Asnaashari et al. 2013; Métivier et al. 2013; Assis & Schleicher 2021). Within the context of iterative least-squares one-way wave-equation migration (LS-WEM), while Lu et al. (2018) obtained an approximate Hessian by implicitly formulating the modelling operator and its adjoint, Abolhassani & Verschuur (2024) introduced a computationally affordable, depth-dependent approximate Hessian, referred to henceforth as pre-conditioned LS-WEM (PLS-WEM).
This paper aims to increase the resolution power of standard ORWI tomograms, expanding upon the findings of Abolhassani & Verschuur (2024) on depth-dependent pre-conditioning for image reconstruction. While employing depth-dependent pre-conditioners for high-resolution image reconstruction with standard ORWI, we concurrently develop the required mathematical groundwork to integrate the depth-dependent pre-conditioning concept for background velocity estimation. This two-fold effort facilitates more accurate velocity estimation with ORWI. We refer to it as high-resolution ORWI (HR-ORWI).
In the following sections, we will first review the theory of HR-ORWI. Then, we will present a brief overview of the depth-dependent gradient pre-conditioning concept for migration (PLS-WEM). Following that, we will mathematically demonstrate its relevance for tomography within the context of ORWI. Next, we will provide three numerical examples, where the first two compare standard ORWI and HR-ORWI, and the third validates that our cost-friendly high-resolution updates in both migration and tomography loops of HR-ORWI result in a reliable background velocity estimation for the Marmousi2 model. Lastly, we end with discussions and conclusions.
2 HR-ORWI THEORY
The theoretical aspects of HR-ORWI are presented here. We rely on Berkhout (1981, 1982) for the forward problem theory.
2.1 Forward problem
The two-way acoustic wave equation for a homogeneous medium in the space-frequency domain can be written as
in which u represents the monochromatic pressure wavefield, and c is the velocity.
Replacing the differentiations in eq. (2) with spatial convolutions (Berkhout 1982) gives
where
After splitting the total pressure wavefield into up- and downgoing wavefields (
where z is pointing downward,
Inserting the derived partial derivative wavefields into the Taylor expansion for
where
When dealing with extrapolation distances with lateral velocity changes, a matrix formulation is adopted to represent eq. (7) as
where

(a) Schematic representation of downward wavefield extrapolation between two virtual depth levels, considering lateral velocity variation. This graphic represents one single row of
By discretizing the subsurface into N virtual depth levels separated by
with
where
2.2 Optimization problem
HR-ORWI minimizes a two-parameter error function by alternating updates of the two parameters: angle-independent upward reflectivity and background velocity (see Algorithm 1). The error function through a linearized forward problem approximation is defined as a quadratic function of frequency-dependent model perturbation vectors
with
and
where
Setting the derivative of eq. (13) with respect to
which leads to
with
and
where
Following Jang et al. (2009) and Oh & Min (2013), the total model perturbation vector, after accounting for the negative sign indicating the descent direction, can be expressed as the linear combination of all the monofrequency model perturbation vectors
where
In the
where
resulting from minimizing the residual data with respect to the optimization step length.
To reduce the computational load of solving the multiparameter inverse problem, we first discard the off-diagonal blocks of the monofrequency approximate Hessian, which represent the coupling between parameter classes. This is achieved by setting the off-diagonal elements of
where
In the following section, we present how to calculate
2.3 Depth-dependent gradient pre-conditioning
2.3.1 Reflectivity
Since PWMod allows access to the acoustic wavefield at different virtual depth levels,
and
where
Within the LS-WEM context, Abolhassani & Verschuur (2024) showed that the sensitivity matrix with respect to the upward reflectivity model at a given virtual depth level (
in which
Integrating eqs (24), (25), and (26) into eq. (17) and then into eq. (20) yield the total angle-independent upward reflectivity perturbation vector
2.3.2 Background velocity
By similarity, we will here derive the sensitivity matrix for the background velocity model at a given virtual depth level (
We rewrite the monochromatic primary reflection data recorded at the Earth’s surface (eq. 10) into the superposition of two upgoing wavefields, one travelling from
The data sensitivity with respect to a single background velocity model parameter located at
After reintroducing
and plugging them into eq. (28), we get
in which given the assumption that
and
Fig. 4 provides a schematic representation of the physical interpretation of eqs (31) and (32) for a medium consisting of six virtual depth levels, including
Now, after eq. (24), we can write,
and
where
Incorporating eqs (34), (35), and (28) into eq. (17) and subsequently into eq. (20) give the total background velocity perturbation vector
According to eqs (25) and (35), the depth-marching nature of HR-ORWI framework facilitates constructing monofrequency, depth-dependent approximate Hessian operators. Each monofrequency depth-dependent approximate Hessian operator is calculated for model parameters located on a single virtual depth level (with
2.4 HR-ORWI versus standard ORWI
Standard ORWI leverages the same forward problem theory as HR-ORWI and similarly updates reflectivity and background velocity by minimizing the errors in primary reflection waveforms for each class of parameters alternately, meaning that both standard ORWI and HR-ORWI follow the same inversion workflow, i.e. Algorithm 1.
Unlike HR-ORWI, standard ORWI minimizes eq. (13) as a quadratic function of
where
We now use a two-layer model to evaluate the fidelity of the reflectivity and background velocity gradients in HR-ORWI versus Standard ORWI. The upper layer has a homogeneous background velocity of 3500 m s−1, while the bottom layer has a homogeneous background velocity of 3000 m s−1. The interface is at a depth of 800 m. 41 shots are triggered every 75 m on the surface. A 10 Hz Ricker wavelet is used as the source function with a fixed-spread acquisition geometry (the maximum available offset is 3000 m). We fix the migration offset—data offsets used for seismic migration—at 400 m (near-offset imaging) and use an initial background velocity model of 3400 m s−1 (i.e.

Stacked reflectivity perturbation and one-source–receiver-pair background velocity perturbation in a two-layer model after one cycle (1× LS migration and 1× LS tomography). (a, b) Standard ORWI: (a) Stacked reflectivity perturbation and (b) One-source–receiver-pair background velocity perturbation. (c, d) HR-ORWI: (c) Stacked reflectivity perturbation and (d) One-source–receiver-pair background velocity perturbation.
Comparing Figs 5(a) and (b) (standard ORWI) with Figs 5(c) and (d) (HR-ORWI) reveals that HR-ORWI provides significantly better resolution in both the image and tomogram than standard ORWI. In Figs 5(a) and (b) (standard ORWI), the signature of the source function is visible in both the image and tomogram. The pair of transmission-after-reflection wave paths in Fig. 5(b) (standard ORWI) are unfocused and band-limited, as manifested by the scattering of reflected energies over a larger area compared to Fig. 5(d) (HR-ORWI), which makes information extraction from the data more challenging. In contrast, Figs 5(c) and (d) (HR-ORWI) show superior resolution in both the image and tomogram. This improvement is attributed to the source signature deconvolution effect embedded in the proposed depth-dependent approximate Hessian operators. With better-focused transmission-after-reflection wave paths, HR-ORWI is better positioned to handle the complex interference of the wave paths in challenging geological settings, supporting more reliable search directions than standard ORWI.
With acquisition configurations identical to those used in the experiment associated with Fig. 5, we use the same two-layer model to compare the stacked background velocity perturbations estimated through HR-ORWI and standard ORWI. In the experiment, we use an initial background velocity model of 3000 m s−1 (i.e.

Background velocity perturbation corresponding to the two-layer model, estimated after one cycle (1× LS migration and 1× LS tomography) for (a) standard ORWI and (b) HR-ORWI.
3 NUMERICAL EXAMPLES
First, we validate the effectiveness of HR-ORWI compared to standard ORWI using two synthetic examples. Second, we evaluate the performance of HR-ORWI with synthetic data on the marine-environment Marmousi model (Martin et al. 2006). In all the examples, we adopt the offset-selection strategy presented in Appendix A.
3.1 Fault model
In the fault model example (Fig. 7a), the observed data is generated via PWMod (only primary reflections), employing a 10 Hz Ricker wavelet. The fault model is discretized with 201 gridpoints in the horizontal direction (20 m interval) and 201 gridpoints in the vertical direction (5 m interval). 41 shots with 100 m spacing are used, and each shot is recorded by 201 receivers with a fixed-spread acquisition geometry at the surface of the Earth. The maximum available offset in the data is limited to 3000 m. The data is recorded for a duration of 1.6 s. The initial background velocity model is a 1-D linearly increasing gradient model (Fig. 7b), and the initial reflectivity model is zero. We fix the migration offset at 1000 m, use the mid-to-far-offset residual data (

Fault model example. (a) True background velocity model. (b) Initial background velocity model.
The 0–30 Hz frequency band data is inverted at once, that is the inversion process does not involve a multiscaling strategy (Bunks et al. 1995). Each cycle of ORWI includes 1× LS migration and 1× LS tomography. While Fig. 8 provides a comparison of the estimated background velocity perturbations (

Background velocity perturbations associated with the fault model example estimated after one cycle (1× LS migration and 1× LS tomography) of ORWI. (a) Background velocity perturbation estimated via standard ORWI. (b) Background velocity perturbation estimated via HR-ORWI.

Estimated tomograms associated with the fault model example. (a) and (c) Estimated tomograms after 25 and 50 cycles (each cycle involves 1× LS migration and 1× LS tomography) of standard ORWI. (b) and (d) Estimated tomograms after 25 and 50 cycles (each cycle involves 1× LS migration and 1× LS tomography) of HR-ORWI.
Comparing Figs 8(a) and (b) reveals a clear optimization direction with HR-ORWI shortly after the first cycle, where the fault line and balanced-amplitude background velocity perturbations are visible from shallow to deep. Emphasizing the reliability of HR-ORWI, this comparison also highlights the potential confusion that might be introduced with standard ORWI.
By comparing Figs 9(a) and (c) to Figs 9(b) and (d), respectively, it becomes evident that the estimated tomogram via HR-ORWI after 25 cycles closely aligns with, or even excel, the estimated tomogram via standard ORWI after 50 cycles. Further comparison also confirms the existence of persistent and non-physical artifacts in tomograms estimated via standard ORWI from the first to the fiftieth cycle. The accuracy of the tomogram estimated via HR-ORWI after 50 cycles (Fig. 9d) is self-evident.
Fig. 10 (convergence history) again validates that HR-ORWI exhibits faster convergence and is able to achieve a better data fit in fewer cycles compared to standard ORWI.

Convergence history (data error curve) associated with the fault model example.
3.2 Reservoir model
In the reservoir model example, a buried low velocity reservoir lies beneath a lens-shaped geological formation with high velocity (Fig. 11a), discretized into a grid of 251 points horizontally and 121 points vertically, both spaced at 10 m intervals. The observed data set is acquired with a fixed-spread acquisition geometry, including 26 shots with 100 m spacing and 251 receivers, both distributed at the surface of the Earth. The data set contains 1.4 s records, generated with a 15 Hz Ricker wavelet as the source function using PWMod, including only primary reflections. The maximum available offset in the data set is limited to 2000 m. The initial background velocity model is a 1-D linearly increasing gradient model, illustrated in Fig. 11(b), and the initial reflectivity model is zero. We fix the migration offset at 1000 m (MEMO), use the mid-to-far-offset residual data (

Reservoir model example. (a) True background velocity model. (b) Initial background velocity model.

Estimated tomograms associated with the reservoir model example. (a) and (c) Estimated tomograms after 25 and 50 cycles (each cycle involves 1× LS migration and 1× LS tomography) of standard ORWI. (b) and (d) Estimated tomograms after 25 and 50 cycles (each cycle involves 1× LS migration and 1× LS tomography) of HR-ORWI.
Comparing Fig. 12(a) with Fig. 12(b) reveals that standard ORWI, after 25 cycles, still struggles to distinguish between the lens body and the underlying flat layer, whereas HR-ORWI reliably resolves this. Furthermore, while standard ORWI completely misestimates the background velocity of the reservoir body, HR-ORWI shows an accurate retrieval.
Comparing Fig. 12(c) with Fig. 12(b) shows that the tomogram estimated with HR-ORWI after 25 cycles not only aligns with the tomogram estimated with standard ORWI after 50 cycles but also surpasses in accuracy, especially at the reservoir level; this is observed while both share nearly identical data errors according to Fig 13. The accuracy of the tomogram estimated via HR-ORWI after 50 cycles (Fig. 12d) is self-evident.

Convergence history (data error curve) associated with the reservoir model example.
Fig. 13 (convergence history) again verifies that HR-ORWI comes with faster convergence and can reach a better data fit in fewer cycles compared to standard ORWI.
3.3 Marmousi2 model
The Marmousi2 model example (Fig. 14) contains 334 × 103 gridpoints with intervals of 22 m. The observed data is generated via PWMod (only primary reflections) with a 10 Hz Ricker wavelet as the source function. With a shooting interval of nearly 200 m, 41 shots are triggered at the surface of the Earth, and each shot is recorded by 334 receivers planted on the surface of the Earth (fixed-spread acquisition). The data is recorded for 4.096 s, and the maximum available offset is limited to 4 km. To build the initial background velocity model, the true background velocity model is smoothed out by applying a square Gaussian kernel with a standard deviation of 1100 m. Subsequently, a single vertical profile is chosen to construct the initial 1D model with a water layer on top, as shown in Fig. 14(b). The initial reflectivity model is set to zero. We fix the migration offset at 2 km (MEMO), use the mid-to-far-offset residual data (

Marmousi2 model example. (a) True background velocity model. (b) Initial background velocity model.
To invert the data within the 0–18 Hz frequency band, a multiscaling strategy is employed, involving three frequency bands: 0–5, 0–11, and 0–18 Hz. For both 0–50 and 0–11 Hz, the inversion process involves 20 cycles, each with 1 least-squares LS migration iteration and 1 LS tomography iteration. The image does not reset to zero after each cycle. For the 0–18 Hz band, the inversion process involves a total of 60 cycles, each with 1× least-squares LS migration iteration and 1× LS tomography iteration, and no image reset after each cycle. Note that the reflectivity model is reset back to zero after each frequency scale. The estimated background velocity perturbations (cycles 1 and 100) and tomogram (cycle 100) through HR-ORWI are displayed in Figs 15 and 16, respectively.

Background velocity perturbations associated with the Marmousi2 model example. (a) Background velocity perturbation estimated after the first cycle of HR-ORWI. (b) Background velocity perturbation estimated after the 100th cycle of HR-ORWI.

(a) Smoothed true background velocity model of Marmousi2 generated via a square Gaussian kernel with a standard deviation of 75 m and a water layer on top. (b) Estimated tomogram associated with the Marmousi2 model example after 100 cycles of HR-ORWI.
The estimated background velocity perturbation after the first cycle (Fig. 15a) illustrates that the pre-conditioned tomographic gradient reliably points in the right direction even from the outset. The main reflections have been correctly interpreted. The tomographic update focuses primarily on reflected events within the depth range of 0–1.5 km. Fig. 15(b) shows the estimated background velocity perturbation after the 100th cycle. It is evident that the tomographic gradient of HR-ORWI successfully samples the low-to-intermediate model wavenumbers of the true background velocity properly in 100 cycles from shallow to deep, though with a depth-dependent decreasing resolution (Williamson 1991).
The estimated tomogram after 100 cycles is shown in Fig. 16(a). In addition, a smoothed version of the true background velocity model with the aid of a square Gaussian kernel with a standard deviation of 75 m and a water layer on top is generated and shown in Fig. 16(a). A comparison of the two models presented in Fig. 16 indicates that the estimated tomogram has a resolution approximately similar to that of the smoothed true model.
To check the quality of the estimated tomogram, different vertical profiles located at the lateral locations of 3000, 3750, and 4600 m are extracted and plotted out in Fig. 17. It can be seen that the estimated background velocity in all extracted lateral locations is smooth but significantly close to the true velocity.

Vertical background velocity profiles associated with the Marmousi2 model example. A vertical profile of the estimated tomogram at the lateral location of (a) 3000 m, (b) 3750 m, and (c) 4600 m.
To assess the extent of uplift in the background velocity from cycle 1 to cycle 100, various images using different background velocity models, including true, initial, smoothed and estimated, are presented in Fig. 18. Upon comparing the images mapped with the initial and estimated background velocity models (Figs 18b and d), a significant kinematic uplift can be observed in Fig. 18(d), attributed to the enhanced focus of the reflectors. After comparing the images mapped with the true and estimated background velocity models (Figs 18a and d), a significant kinematic uplift is again confirmed in Fig. 18(d) based on the reasonably accurate positioning of the reflectors. Comparing the images mapped with the smoothed and estimated background velocity models (Figs 18c and d) also suggests that a “maximum” kinematic uplift within the achievable resolution of RWI has occurred in Fig. 18(d).

Images using different background velocity models associated with the Marmousi2 model example. The mapped images correspond to the (a) true, (b) initial, (c) smoothed and (d) estimated background velocity models.
4 DISCUSSION
The one-way (up and down) wave-equation extrapolation technique has persisted in the seismic imaging industry despite its weaknesses (for details on weaknesses, see Mulder & Plessix 2004), as it offers computational efficiency when contrasted with direct solutions to the Helmholtz equation. In this research, we demonstrated the potential use of one-way wave-equation forward modelling, that is PWMod, in efficient gradient pre-conditioning to enhance suboptimal images and tomograms resulting from RWI. In terms of constructing, inverting or storing, depth-dependent approximate Hessians computed based on PWMod are more feasible than doing so for pre-conditioners computed at once for all model parameters within the medium. Despite being feasible, it remains time-intensive compared to standard ORWI (i.e. HR-ORWI is six to nine times more expensive but requires half the iterations of standard ORWI, effectively making it three to five times more expensive overall), given the current implementation and the capabilities of the current computing resources. However, there are ways to reduce the time intensity. Among others, source subsampling in the computation of
Selecting one technology over another often involves evaluating the cost-versus-quality trade-off. While cost-effectiveness must remain a key priority, compromising accuracy due to budget constraints in exploration seismology can raise significant concerns. In the reservoir model example presented in this paper, we clearly observed how choosing a budget-conscious technology (i.e. standard ORWI) can lead to incorrect subsurface illumination, resulting in misinterpretations of subsurface properties and potentially far-reaching consequences. Using cheaper technology can indeed be a double-edged sword; while it saves money upfront, the outcome could be irreversible. This portrays the complexity of decision-making dynamics in exploration seismology.
Although we did not use any model-based regularization in this study, integrating approaches like structure-oriented regularization guided by the reflectivity model may further enhance the fidelity of tomographic updates in HR-ORWI (Masaya & Verschuur 2018; Yao et al. 2019; Provenzano et al. 2023).
Real data applications require careful consideration of several additional factors. The negative potential impact of source wavelet estimation errors and the absence of low-frequency data on the stability of HR-ORWI should be carefully explored. Furthermore, addressing amplitude versus offset (AVO) effects in the ORWI process is crucial. Qu & Verschuur (2021) recommend minimizing AVO effects before applying ORWI. In this study, through inverse-crime numerical examples, we assumed reduced AVO effects in the data, facilitating a thorough examination of ORWI convergence. Challenges may also arise in accommodating anisotropy effects and integrating the Q-effect into ORWI, as highlighted in other studies (Alshuhail & Verschuur 2019; Safari & Verschuur 2023).
5 CONCLUSIONS
Limited or unresolved resolution in both migration and reflection tomography components of RWI can lead to suboptimal velocity estimation. To address this within the context of ORWI, we introduced a pre-conditioning strategy that enhanced both the migration and tomography loops, resulting in high-resolution updates for both reflectivity and background velocity. We named this approach HR-ORWI, where diagonal blocks of the monofrequency Gauss–Newton Hessian operators, each corresponding to a specific virtual depth level, are extracted and inverted to pre-condition their corresponding monofrequency gradients. This approach performs a partial deconvolution of the gradient while keeping memory and computational requirements manageable. Through three synthetic examples, we demonstrated that HR-ORWI generates superior background velocity estimations compared to standard ORWI.
ACKNOWLEDGMENTS
The authors express their gratitude to the sponsors of the Delphi consortium for their stimulating discussions, valuable feedback and financial support.
DATA AVAILABILITY
The code of this research is confidential. All data generated in this research can be shared upon a reasonable request directed to the corresponding author.
REFERENCES
APPENDIX A: OFFSET-SELECTION STRATEGY
Despite its conceptual appeal and aside from the resolution challenges, ORWI could face other limitations due to other factors: (1) Full-wave inconsistency in the short-offset residual data for tomography due to inconsistent reflectivity and velocity models (Valensi & Baina 2021; Liang et al. 2022; Provenzano et al. 2023). (2) The standard practice in the RWI community for mitigating the adverse imprint of inconsistent reflectivity and velocity models often involves using short-offset or near-zero offset data for migration (Provenzano et al. 2023). However, this can degrade the illumination and signal-to-noise ratio of the reflectivity model and result in an increased number of least-squares migration iterations. (3) Adverse contribution of cycle-skipped long-offset reflection data in tomography (e.g. Abolhassani & Verschuur 2022; Provenzano et al. 2023).
A1 Tomography offset: short offsets
To minimize the impact of inconsistent reflectivity and velocity models on the tomographic gradient of ORWI, we propose muting their relevant erroneous tomographic wave paths in the residual data gathers, which are typically established at short offsets.
A2 Migration offset
Building upon the muting of short-offset residual waveforms for tomography, we extend the migration offset to improve both the illumination and the signal-to-noise ratio of the reflectors. We follow Algorithm 2, aiming for an extended offset for migration (rather than short offset) where a significant reduction in the misfit function is impracticable. The misfit indeed remains unchanged as the amplitude fit fails to better itself while mapping more out-of-phase reflection data into the model space. This denotes the MEMO, lying averagely between the short offset and the maximum uncycle-skipped offset in the data. Algorithm 2 is applied as a quality control measure prior to ORWI. This algorithm calculates the misfit function value for various offset ranges up to the maximum offset within the acquisition. Following this, on the misfit-offset graph, we pick an offset right before the curve starts to become nearly horizontal on a logarithmic scale.
A3 Tomography offset: long offsets
A given pair of modelled and observed seismic traces in the time domain are called cycle-skipped traces if their time distance is larger than half a cycle time shift (Virieux & Operto 2009). Inspired by van Leeuwen & Mulder (2010), we defined Algorithm 3 as a time-domain two-step data-selection algorithm excluding the contribution of the cycle-skipped long-offset residual data in the tomography loop of ORWI to obtain a tomographic update free of the damaging effect of cycle-skipping in each iteration. In Algorithm 3, we identify a given pair of time-domain modelled and observed seismic traces as cycle-skipped if their maximum correlation lag exceeds a reference lag determined by the dominant period of the observed trace. Based on this criterion, our time-domain data-selection algorithm relies on the cross-correlation of these traces within a sliding window (referred to as local cross-correlation), safeguarded by a preliminary global cross-correlation of the given traces. This global cross-correlation functions as a protection, protecting subsequent local cross-correlations from irregularities in the modelled waveform. Without this safeguard, the outputs of the local cross-correlations would not meet the required accuracy.