Local events-based fast RTM surface-offset gathers via dip-guided interpolation

Reverse Time Migration (RTM) Surface Offset Gathers (SOGs) are demonstrated to deliver more superior residual dip information than ray-based approaches. It appears more powerful in complex geological settings, such as salt areas. Still, the computational cost of constructing RTM SOGs is a big challenge in applying it to 3D field data. To tackle this challenge, we propose a novel method using dips of local events as a guide for RTM gather interpolation. The residual-dip information of the SOGs is created by connecting local events from depth-domain to time-domain via ray tracing. The proposed method is validated by a synthetic experiment and a field example. It mitigates the computational cost by an order of magnitude while producing comparable results as fully computed RTM SOGs.


Introduction
Classic ray-based migration velocity analysis (MVA) mainly involves three major steps: (1) constructing offset or angle domain gathers from Kirchhoff migration, (2) picking residual moveouts from these migrated gathers, and (3) inverting the updated velocity model using ray tracing. In complicated geological regions, ray-based imaging methods such as Kirchhoff migration usually fail to provide high-quality gathers. In contrast, reverse time migration, a method based on the two-way wave equation, can handle complex velocity overhang, has no dip limitations, thus can provide highquality gathers for accurate estimation of residual moveout. Ehinger et al. (1996) initially performed the common-offset wave-equation migration by generating Green's functions from the surface and then used them in a frequency-dependent one-way migration. Etgen (2012) extended this method from 2D to 3D and optimized the computational efficiency by recycling one-way Green's functions via extrapolating them from the surface. Giboli et al. (2012) upgraded Green's function from one-way to two-way and created SOGs via stationary-phase extraction of an encoded attribute. Yang et al. (2015) proposed a simple but expensive method. This method divided a common-shot gather into offset sets (Duquet and Lailly 2006) and back-propagated the wavefield for each offset set separately. Zhao et al. (2020a) extended it to the elastic cases and perform a comprehensive elastic migration velocity analysis. Each offset gathers can survive the destructive interference from neighboring offsets; therefore RTM SOGs can provide reliable moveouts since each source-receiver section is migrated independently (Zhao et al. 2018).
All previous works were suffered from significant computational cost, which prevents this technology from being widely used in practice, even for 2D applications. For Edited by Jie Hao and Chun-Yan Tang 1 3 example, for one shot, RTM SOGs with 20 offsets are about 10 times more expensive than regular RTM with crosscorrelation imaging conditions. As a result, mitigating the SOGs cost is particularly essential and is the main objective of this article. One potential way to mitigate this problem is to use sparse offsets. A sparse offset system can significantly reduce the computational cost but can produce migration operators. To compensate for the sparse offset, an effective interpolation is critical to make this proposed workflow feasible in practice. Various approaches were established for data interpolations in the geophysical community (Yilmaz 2001), in which a variety of data volumes to interpolate data at the target positions to regularize and increase the spatial density of the geophysical datasets. For example, the frequency-space (FX) domain (Spitz 1991) interpolation, based on complex spatial prediction filters (Ronen 1987), does not require any attempt to determine the geological dips. In contrast, many interpolation methods were developed along with geologic structures, such as the energy-scanning local dip method (Marfurt 2006), local structure tensors (Fehmers and Höcker 2003), and plane-wave destruction (PWD) filters (Fomel 2002b). In this study, we select PWD as the interpolation method since most gathers present a non-complex curvature structure.
In summary, the image gathers at additional offsets are generated by structural-dip-guided interpolation. The proposed workflow consists of the following five steps: a SOGs computations. Each shot-domain seismic records is divided into different groups with sparse offsets to mitigate the cost. RTM is performed on each offset group independently. b Ray tracing in the depth domain. Rays are shooting from migrated events up to the corresponding shot/receiver locations. These parameters, such as geological dip, reflection angle, locations, and ray travel time are determined by this ray-tracing process. c Dip calculations in the time domain. The process of depth-to-time mapping provides time-domain parameters such as recorded time, source and receiver locations, ray-traced emerging angles at source, and receiver. We may, therefore, estimate the observed dip of local events at each sparse offset by PWD. d Residual dip calculation. We then calculate the misfits of surface angles between raytraced and observed results and yield the residual dips at these sparsely offsets. e Gather interpolation. Finally, we perform the PWD interpolation using the obtained residual dips as a guide and rebuild the gathers at denser spatial samplings.

Local reflection events from RTM SOGs
Similar to Kirchhoff migration, offset-domain partial images are straightforward to produce by RTM. The migration integral is performed on separate subsets of the data and the partial images that are kept separately. In 2D, it becomes the following integral over time: where I(x, z, h) is the surface-offset image gathers. x and z are the image-point location. src and rec denote the source wavefield and receiver wavefield, respectively. T is the record duration. rec (x + h, z, t) gives the receiver wavefields using data at each offset h . Equation (1) is performed independently on each offset h and the size of I(x, z, h) depends on the offset selection. SOGs are well-known essential outputs of pre-stack depth migration used for velocity estimation. The local reflected events are selected on the migrated gather volume. They are characterized by their dips in the commonoffset panels and residual moveout curve. Migration velocity analysis aims to invert the velocity model by flattening the residual moveouts of picked local events. Definitions of parameters associated with events x, z, h to ray trajectories and residual moveouts in the context of migration velocity analysis are summarized in Fig. 1. As observed in Fig. 1, the following requirement of the travel time needs to be satisfied  (src,x,z) T calc (rec,x,z)

Fig. 1
Geometry of a local event with geological dip and its residual moveout from common-image gather via ray tracing. The depth perturbation δz is exaggerated to depth z, offset h for illustration purpose. The far-offset source (src) ray, with the travel time T calc (src, x, z) (the left red solid line), forms the incident angle θ to the normal direction (the purple dot line), and so as the receiver (rec) ray, with the travel time T calc (rec, x, z) . In contrast, the far-offset ray of the observed seismic data (the blue lines), with the travel time T obs (src, rec) , shoot from a shallower reflector and reach the src and rec at the same h. The near-offset ray path (yellow solid line) forms similarly and share the same reflector with the observed ray path. The corresponding CIGs range from the near-offset (yellow lines) to the far-offset ray (red lines) presenting a residual dip ϕ (the solid green lines) is shown on the right. By performing tomographic updates for slowness χ, the local event is shifted vertically by δz to the location of the observed event by flattening the residual dips with a local event T obs (src, rec) migrated to a 2D position x, z, h: where T calc (src, x, z, ) and T calc (rec, x, z, ) are the calculated one-way travel time (dependent of slowness ) from the source and the receiver locations to the image location x, z , respectively. Similarly, T obs (src, rec) is the observed two-way travel time (independent of slowness ) for the given event from the seismic records in the time domain. The depth perturbation z is applied to depth z , offset h for illustration purpose. Figure. 2a shows a synthetic example of collected SOGs from Fig. 1 of different CDPs (common-depth points). It shows the relationship of different dips of one local event in three dimensions (CDP, offsets, and depth), specifically, is apparent geological migrated dip measured in the common offset axis and is the residual dip in the SOG. and are defined as: For the correct velocity, the SOG should be flat and should be zero. The dip information (blue lines in Fig. 1) is given by In tomography, Stork (1992) established that the relationships mapping traveltime perturbations to depth deviations are T calc (src, x, z) + T calc (rec, x, z) = 2 z cos cos . The source and receiver locations are switched to the bout 1985). Equation (1) can be rewritten as where (p obs src − p calc src ) − (p obs rec − p calc rec ) is the surface angle misfit caused by the velocity errors, and 2 cos cos is a stretching factor at the place of the local event (x, z, h). Chauris et al. (2002) derived this connection between local events in the pre-stack time domain and the events in the depth domain. A full RTM SOG using Equation (5) is fairly expensive. Figure. 2b shows a sparse and muted version of Fig. 2a with selected offsets. The gathers are generated with only 1/10 cost comparing to a full operation but an effective interpolation is required to compensate for the sparsity. A dip-guided interpolation is selected for the relatively simple structure as CIGs. Equation (5) is the key equation providing the residual dip to perform gather interpolation.
These measurements p calc src , p calc rec , cos , cos are determined by ray tracing in depth-domain (red lines in Fig. 1). Specifically, rays are traced from each image point to the source and receiver locations for the given dip and reflection angle. Each ray possesses five parameters: apparent geological dip cos , reflection angle cos , horizontal and vertical positions x, z , and ray travel time T calc (src, rec) . Source-side dip p calc src and receiver-side dip p calc rec can be determined during ray tracing and used for the depth delay estimation. In contrast, these time-domain measurements p obs src , p obs rec require designed techniques to be directly measured in the common shot and receiver domains, respectively (red lines in Fig. 2a). In the next section, we will discuss the dip estimation and how to use the dip to perform the dip-guided interpolation.  Fig. 2 Geometry relationship of the local migrated event from Fig. 1 mapping from time to depth time. a Definitions of the residual dip and apparent geological dip associated with a local event in CIGs in depth domain. The Y axis represents the CDP location, whereas the XZ plane is equivalent to Fig. 1, with is apparent geological migrated dip measured in the common offset (Y) axis, and is the residual dip in the CIGs. b Sparse SOG with 1/8 of the offset groups as in (a) The depth perturbation z is exaggerated to depth z , offset h for illustration purpose. The far-offset source (src ) ray, with the travel time T calc (src, x, z) (the left red solid line), forms the incident angle to the normal direction (the purple dot line), and so as the receiver ( rec) ray, with the travel time T calc (rec, x, z) . In contrast, the far-offset ray of the observed seismic data (the blue lines), with the travel time T obs (src, rec) , shoot from a shallower reflector and reach the src and rec at the same h . The near-offset ray path (yellow solid line) forms similarly and share the same reflector with the observed ray path. The corresponding CIGs range from the near-offset (yellow lines) to the far-offset ray (red lines) presenting a residual dip (the solid green lines) is shown on the right. By performing tomographic updates for slowness , the local event is shifted vertically by z to the location of the observed event by flattening the residual dips.

Dip estimation and dip-guided interpolation
Local reflected events, either in time-domain or depthdomain gathers, contain not only temporal and spatial information but also dip information. Seismic dips can be used for velocity model building, gather interpolation, and conditioning. These measurements of seismic dips can be completed by PWD filters and are valuable tools for numerous geophysical applications (Fomel 2002a(Fomel , 2010Fomel and Guitton 2006). PWD predicts the next trace by phaseshifting from the previous trace along with the dominant event dip while maintaining the amplitude. The estimated dominant dip is then determined by minimizing the prediction error (Fomel 2010;Xue et al. 2019). Shaping regularization limits the assessed dips to change smoothly in neighboring traces for the interpolated records. We can formularize it as a regularized least-squares problem as obs ≈ , where is the dip field between traces, associated with the PWD operator, the delta is the prediction error, and obs is the original seismic trace in time-domain. This approximated equality can be achieved by solving the least-squares problem ≈ (( obs ) T obs ) −1 ( obs ) T .
Dip calculation of time-domain gathers is challenging due to its low S/N ratio. Preprocessing steps should be applied before the estimation for accurate dip extraction. The preprocessing steps include low-pass filtering, removal of coherent noise, and estimation in the common-offset domain with dip constraints within a realistic range. Figure 3a shows local migrated events in the corresponding time domain with record time T obs (src, rec) . Events are identified by the soft thresholding of the semblance within sliding windows. A simple example of the resultant picking of a common-shot gather is illustrated in Fig. 3a. Figure 3b illustrates the calculated dip from the input data and provides p obs src , p obs rec for Equation (5). Other measurements p calc src , p calc rec , cos , cos can be achieved via ray tracing in the depth domain. We can then obtain tan using Equation (5). Parameter calculations are the most crucial step in the entire workflow.
The next task is to use tan to perform dip-guided interpolation to fill in the missing offset. We again need to solve the regularized least-squares problem by switching matrix position ≈ , where is still the dip field containing tan but is sparse SOGs. Instead of solving can be computed similarly but in common receiver gathers. b Dip estimation and event identification using PWD. Red and blue represent the slope of local events dipping in different directions for , we solve for this interpolation purpose. This linear problem can be solved efficiently by ≈ ( T ) −1 . To demonstrate this dip estimation and interpolation process, a classic synthetic image (Claerbout and Fomel 2006) is selected to validate our proposed dip-guided interpolation procedure. This classic image contains dipping beds, an unconformity, and a fault. Figure 4a shows the subsampled version (1/4 sampling rate), and we want to recover the original image from this image. Figure 4b shows the dip estimated from the input image, and it has been smoothed similarly to generate a dip that has the same dimension as the target image (Zhao et al. 2020c). We can easily obtain the interpolation result following the derived equations above. Figure 4c and d show the interpolation images and their difference from the true reference image, respectively. The interpolated results act as a smoother version of the inputs, and residuals only around the unconformity and fault can be observed. A comprehensive analysis can be found in Zhao et al. (2020b).
As we discussed earlier, PWD is suitable for the gathers with simple curvatures since this technique assumes smooth amplitude variances and a single plane-wave direction. Other advanced interpolation methods, such as cubic spline interpolation (Zhang and Zheng 2014), can also be used. Cubic splines offer more stable and smooth results as well as continuous second-order derivatives for given input points. Cubic splines thus produce continuous results for smoothly varying shapes. The shape of migration gathers along offsets are generally varying smoothly and generate continuous spatial derivatives. Taking advantage of highorder nonlinearity, the cubic-spline interpolation may serve as a great option to interpolate gathers with complex curvatures (e.g., mixed up-down curve). Figure 5 demonstrates this dip-guided interpolation process, where Fig. 5a is a selected sparse SOG from Fig. 2b. The ray-tracing procedure described in the previous section was carried out for each local event individually. The resulting values overlapping with its associated event in sparse  Fig. 5b. The sparse needs to be populated with the desired offset grid. We applied linear interpolation to fill in in the offset gaps, followed by horizontal linear smoothing. Figure 5c shows the interpolated results using in the sparse offset grid as the initial dips with the PWD. A total of 81 offsets of SOG have been recovered. The reference SOG is shown in Fig. 5d. The interpolated gather shows a comparable quality and preserves the slope comparing with full SOGs. The proposed method successfully recovers the missing traces and removes the aliased events, but with only one-eighth of the computation cost.

Examples
In this section, we aim to demonstrate the proposed method using synthetic and field data examples. The first example is a transition-zone model with a shallow water environment. The model contains sediments with strong velocity anomalies such as salt bodies and islands. 250 shot gathers were preprocessed with refraction-wave muting and FK noise removal. We separated the common-shot data with an offset ranging from -8 to 8 km into 11 groups with a spacing of 1.6 km. Similar to the example in Fig. 5, the migrated CIGs are shown in Fig. 6 generated the input SOGs by subsampling the original offsets shown in Fig. 6a, which has severe aliasing issues in the offset dimension; Fig. 6b shows the dip estimated with the PWD filter. We used it to perform dip-oriented interpolation, resulting in the result shown in Fig. 6c. To perform high-quality picking of residual moveout for tomography, the target output gathers have 81 offsets with 0.2 km spacing. The difference between the full offsets in Fig. 6d is very small. The interpolation has reconstructed the original SOGs with a much smaller cost. In the second example, a field data experiment with a real transition-zone environment was performed, and the results are shown in Fig. 7. We used the same parameters as the previous synthetic example to examine the robustness of our proposed method. Comparing the interpolated SOGs ( Fig. 7c) with full SOGs (Fig. 7d), the recovery quality is very comparable to the synthetic example, which further demonstrates the robustness of our workflow. The computation time of these two datasets is listed in Table 1. Compared to the computation time for standard GOSs, our proposed method indicates a speedup factor of about 5 times over standard full-volume SOGs. The efficiency is mainly attributed to three factors: sparse offsets we selected, the cheap cost of 3D ray-tracing, and fast dip-guided interpolation.

Conclusions
We proposed a dip-based method to gather interpolations to mitigate the heavy cost of generating RTM SOGs. The main idea of our proposed method is mapping the local events between time and depth domains via ray tracing. The local events are defined by these depth-domain parameters: geological dip, reflection angle, locations, and ray travel time; and time-domain parameters: recorded time, source and receiver locations, ray-traced emerging angles at source, and receiver. Compared with traditional techniques, our proposed workflow accomplishes comparable gather quality, but with considerably lower cost as demonstrated by these synthetic and field data examples. The proposed workflow offers a practically realistic approach for velocity model building in the presence of complex geological settings. The computation cost ratio highly depends on the percentage of sparse offsets to the full, and the optimum number of the sparse offsets are related to curvatures caused by velocity inaccuracy. In the simplest case of a constant velocity error in a homogeneous medium, we may only need two offsets to reconstruct residual moveouts because the curvature is controlled by the tangent function of the incident ray angle.
In the presence of complex velocity error in a heterogeneous medium, however, the gather may experience up-down curvatures rather than a single curve-down (up) moveout. These complicated scenarios require dense offsets among each peak-trough to fully reconstruct the correctness of kinematic information from these gathers. We are currently in a preliminary interpolation stage, in which the offset selection is determined in a heuristic way instead of through systematic optimization. A more robust and self-adaptive approach is one of our ongoing research directions. Dip-guided PWD interpolation process of synthetic transition zone example. Y-axis represents the CDP location in addition to SOG shown in Fig. 5. a RTM SOGs with sparse offsets. Comparing this with Fig. 2, only 10 offsets are selected for RTM SOG computations. Other traces are filled with zeros. b the calculated dip of the sparse offsets using equation (5); c Interpolated results. Seven additional traces were inserted between each of the neighboring input traces. d Original SOG with dense 81 offset traces. e a waveform comparison between our proposed method (green lines) and the original SOG (blue lines) at a selected section. The detailed comparison (e) demonstrates that our method provides a highly comparable quality to the original SOGs and the kinematic moveout matches closely with the proposed method