ГЕОФИЗИКА, ГЕОФИЗИЧЕСКИЕ МЕТОДЫ ПОИСКОВ ПОЛЕЗНЫХ ИСКОПАЕМЫХ A Comparison of the Classical OWWE Migration Methods: Analysis of Feasibility in Subsalt Plays Exploration

Subsalt seismic imaging is particularly difficult due to the complex overburdens, caused by the movement of salt that usually results in steeply dipping structures, which along with strong lateral velocity contrast at the sediments-salt interface distort the structural position and the stratigraphic resolution of the subsalt reservoirs. Despite it was proven that two-way wave equation migration provides the best results illuminating these reservoirs, it has huge computational cost, especially in three dimensions. One-way wave equation (OWWE) migration techniques are good alternative in this case as providing the acceptable quality of seismic sections with an appropriate computational cost. To know the advantages and limitations of the OWWE techniques in subsalt imaging, three classical OWWE algorithms were evaluated for depth migration of prestack data (PSDM): phase-shift-plus-interpolation (PSPI), split-step Fourier (SSF), and Fourier finite-difference (FFD). These algorithms were tested with three different fully elastic synthetic models which simulates the structural complexity showed in subsalt plays. It was concluded that FFD gives very accurate results when the lateral velocity variation was strong with acceptable computational cost. The PSPI provided the best quality results but required about twice the computer time needed for FFD, and SSF was the fastest but clearly the least accurate.


Introduction
Salt flow frequently leads to strongly deformed formations with steeply dipping structures that present significant challenges for the seismic migration methods. Prestack depth migration methods are the most commonly used in structurally complex areas. These are usually put into two separate categories -ray-based methods and wave field extrapolation-based methods (Han, 1998), both derived from the acoustic or elastic wave equation. The ray-based prestack Kirchhoff migration is still the most popular as it is costeffective and can image most of the complex structures. However, this technique has some intrinsic limitations preventing it from handling extremely complex structures characterized by steep dips caused by salt movement (Mulder and Plessix, 2004). Its main limiting factor is the dependence on accuracy of the ray tracing used to get the travel times needed for the migration. The multi-path ray tracing in salt structures may yield improper travel times, causing poor subsurface imaging.
On the other hand, wave field extrapolationbased methods can handle this problem allowing more accurate imaging of the complex structures. In the wave field extrapolation techniques, there are two options, derived from the two-way wave equation and from the one-way wave equation.
Two-way techniques are most accurate but demand a greater computation time. Therefore, the one-way techniques may be a good choice. Claerbout (1971) formulated one-way wave equation (OWWE) for seismic migration involving up and downgoing waves. It does not provide true amplitudes for reflectors since the one-way propagator does not properly model the wave amplitudes and computes only wave fields with propagation angles less than 90°. However, with some modification to the OWWE like proposed by Vivas et al. (2010) and Zhang et al. (2006), this problem can be solved.
Some OWWE migration methods use the Fourier transform to solve the wave equation in the frequency-space and frequency-wavenumber domain called usually as dual-domain migration methods. They incorporate the Laplacian operators in the wave equation to reduce the numerical dispersion effects that usually affect the migration methods based on finite differences. Gazdag (1978) introduced the phase shift migration method, which works in the frequencywavenumber domain, and requires a homogeneous velocity model to correctly image the reflectors. In order to solve the problem with strong lateral velocity variation, several methods were introduced. Gazdag and Sguazzero (1984) proposed the phase-shift-plus-interpolation (PSPI), Stoffa et al. (1990) introduced the split-step Fourier (SSF), and Ristow and Rühl (1994) developed a split-step Fourier and finite-difference hybrid method called Fourier finite-difference (FFD).
In this paper, three OWWE prestack commonshot migration methods were evaluated on 2D synthetic data, which are closely related to each other (Rühl and Ristow, 1995), to compare their imaging accuracy in salt structures. Finally, the results of migrating of the 2D seismic records obtained from the full elastic synthetic models are shown.

Theoretical framework
The phase shift migration method is the basis for the three methods evaluated here. Gazdag (1978) proposed it as a more efficient alternative to the finite-difference methods in general. The 2D acoustic wave equation in the ( , ) domain, given a homogeneous velocity medium, is where = ( , , ) is the pressure wave field, is the midpoint variable, is the depth, is the velocity of the medium, and is two-way travel time. Taking the temporal Fourier transform and spatial Fourier transform in the -direction, to get the wave equation in the frequency-wavenumber domain − = . (2) Here is the wavenumber in the lateral direction and is the frequency. Now, taking the spatial Fourier transform in the -direction, one obtains the dispersion relationship Using equation (3) in equation (2), we obtain the one-way wave equations = ± .
The analytic solution of eq. (5) is the wave field extrapolation equation: where ± signs stand for wave field extrapolation direction and involve the phase shift in the frequency-wavenumber domain. When the velocity changes with depth, phase shift migration provides accurate images, but it works only in mediums with weak lateral velocity variation.

Phase-shift-plus-interpolation migration
Phase-Shift-Plus-Interpolation (PSPI) method is an updated version of the phase-shift method, developed by Gazdag and Sguazzero (1984) to handle the strong lateral velocity variations. PSPI is based in the idea of introducing several reference velocities to consider the lateral velocity variation by interpolation among wave fields obtained when they are propagated in the ( , ) domain for each of those reference velocities. These wave fields are brought to the ( , ) domain by the inverse Fourier transform, where the reference wave fields are interpolated to obtain the true wave field based on the relationship of the local and reference velocity.
To avoid decreasing accuracy for small dips, when / ≪ 1, Gazdag and Sguazzero (1984) split the phase shift in equation (6) and Here is the reference velocity ≠ ( , ) and is an approximation to ( , ), * is the Fourier transform of * from ( , ) to ( , ). Clearly, the reference velocities chosen are a crucial factor in this migration method because the computation cost is directly proportional to the number of reference velocities. To handle the computation cost and preserve accuracy was realized by the adaptive choice of Bagaini et al. (1995) for the reference velocities.

Split-Step Fourier migration
Split-Step Fourier (SSF) method was introduced into seismic imaging by Stoffa et al. (1990) for migration of the stacked seismic data, and was implemented in 2-D prestack migration by Roberts et al. (1997).
The SSF method, which is implemented in both the ( , ) and ( , ) domain, is based on separation the velocity field in two parts: a background velocity (which can be some average of the medium velocities) and a perturbation term: Here is the background velocity, which is used to the wave field extrapolation in ( , ) domain, * = ( , , ) exp ± ( ) − .
After taking the inverse Fourier transform, * is transformed into the ( , ) domain, where is applied a phase shift to take into account the lateral velocity variation The SSF method involves a wave field extrapolation with the background velocity in the ( , ) and a phase shift in the ( , ) domain using perturbation terms, just like the PSPI method, but its workflow is just opposite to that of PSPI (for PSPI, the phase shift is done before to the wave field extrapolation). For strong lateral velocity variations, this perturbation term is not enough, and more than one reference velocity is required for SSF, which increases the computational cost to that of PSPI method (Han, 1998).

Fourier finite-differences migration
Though PSPI and SSF methods can deal with lateral velocity variations, they are less accurate when these lateral variations are very strong, which are usual in salt basins, where salt bodies usually have much larger velocities than surrounding sediments. Ristow and Rühl (1994) presented an extension of the Split-Step Fourier method to handle strong lateral velocity variations, which uses a more precise approximation for the dispersion equation by adding an adaptive finite-difference (FD) operator. The expression for the wave field propagator for the Fourier finite-differences (FFD) method is: where symbolizes the minimum velocity at layer [ , + ], and are coefficients with = 2 and = [ + + 1]. The combination of І and ІІ is implemented when the SSF method is used. ІІІ is the finite difference operator.

Methodology
To evaluate the three OWWE migration methods we created three synthetic elastic models with different salt bodies, like salt diapirs, pillows, and nappes, which are overlaying oil and gas reservoirs. The models (including Vp, Vs, density) are 20 km long and 4 km depth. They contain 800 x 4000 grid points with a lateral and vertical spacing of 5 m. The sources (air guns) and receivers (hydrophones) were located along one horizontal line at a depth of 10m below sea level. The data set contains 199 shots; each one contains 500 traces with a recording length of 4 s and sampling interval of 2 ms. The shot and receiver intervals are 100 m and 20 m, respectively. A 24.5 Hz Gaussian derivative wavelet was used as source function. Three models have a water layer of 550 m depth. To reduce the "hard water bottom" reflection and the multiples generated on the migrated sections, one transitional layer with a thickness of 50 m over the original water bottom was added. The structural configurations proposed for the models are the compressional, extensional and un-faulted settings. Several reservoirs were inserted below the salt layers. Their detection was a key parameter when evaluating the migration methods.  (Greenberg and Castagna, 1992;Castagna et al., 1993) All the models are fully elastic, allowing the propagation of compressional, shear, converted and all manner of guided waves by using the module suea2df (Seismic Unix version of (an)elastic anisotropic 2D finite difference forward modeling (Cohen and Stockwell, 2017)), which uses finite differences to solve the wave equation. Create of these elastic models needs to assign compressional velocity, shear velocity, and density to each grid point. Before assigning these properties, the lithology was given to every layer, among which are: water, sandstone, shale, salt, limestone, and dolomite. All layers with hydrocarbons are sandstone. After assigning the lithology, the P-wave velocity was defined considering the compaction gradient k: where is a depth, is the velocity at datum (sea level) and is the P-wave velocity.
Given the P-wave velocity for each layer, the S-wave velocity and density were calculated applying the transforms from (Greenberg and Castagna, 1992) and (Castagna et al., 1993) ( Table 1). Water and salt P-wave velocities are 1.480 and 4.550 km/s, respectively.

Results of study
Elastic models were created as close to real ones as possible to properly evaluate the quality of the results obtained with each method. The program Gassmann's Fluid Substitution v1.0 (Al-Khateb, 2015), which uses the Gasmann's equation for fluid substitution, was implemented to assign the elastic properties of the hydrocarbon bearing layers. Results obtained for each model are shown below.

Model 1
The first model is a representation of a salt structure with a basal thrust fault defining its rising compressional nose was partially adapted from (Harrison and Patton, 1995). The p-wave velocity model is shown in Fig. 1A. The salt body overlays several channel sands reservoirs, located at a depth of 2300 m to 2900 m. The model of oil reservoir trapped under the salt formation is located at 2500 m depth. Here, the salt offers a weak detachment layer of the thrust fault. Over the salt layer, there are several normal faults generated by stretching on the hanging wall. Granitic dike with a top located at 3400 m depth is shown in the bottom of the model.
In Fig. 1, the migrated results for all three algorithms are shown. Specifically, PSPI (Fig. 1B) and FFD (Fig. 1D) methods provided good imaging of the salt shape and subsalt layers. The SSF section (Fig. 1C) does not image properly the salt structure, because the top of salt is not clearly identified, and the base salt reflection is located higher than its real position.
Geological section in Fig. 2 shows the channel sand reservoirs below the salt formation and a fault trap reservoir, which are enclosed in the yellow square in Fig. 1A. Fig. 6, 7 and 8 show the migration results for this cropped section, where the dashed red line indicates the correct position of the salt layer.
PSPI and FFD methods provide acceptable imaging of the channel sands, different from the SSF method, which has problem imaging the smaller channels, and positioning the salt layer 100 m above its real position like indicates the dashed blue line (  Fig. 4) this depth gap occurs in all migrated sections when using the SSF method. Fig.6. shows the fault trap reservoir enclosed in the yellow square in ed by conformity with the local, which can be identified by conformity with the local structures and the polarity change of the reflector presented by oil-water interface. PSPI (  Fig. 7) and FFD (  Fig. 8) methods allow us to see the polarity change. The SSF (  Fig. 9) method masks it due to low signal-tonoise ratio (SNR).

Model 2
The second model is a representation of a rift valley overlaid by a salt layer, which has been deformed to form four diapirs and two salt pillows, partially adapted from (John Perez, 2017). The P-wave velocity model is shown in Fig. 14A.
The salt bodies overlay several channel sand reservoirs located from 2000 m to 2350 m depth in the post-rift rocks. Several oil reservoirs are trapped against the post-rift unconformity in the syn-rift rocks at 2800 m depth. In addition, there are reservoirs in the syn-rift rocks, trapped against the horst and the channel sands in the fault blocks. The pre-rift sediments provides the source rocks and there are no oil or gas reservoirs. No reservoirs were included over the salt layer because it is not an objective of this study.
The migrated sections in Fig. 14 show obtained results for all three algorithms. PSPI method provides the best images of the subsalt structures with the best signal-to-noise ratio. FFD method also demonstrates acceptable results.
The PSPI and FFD methods properly located the base salt, but the salt flanks were poorly illuminated (see Fig. 11 and Fig. 12). The SSF method provides the least accurate results seriously affected by seismic noise, which does not allow identifying some post-rift and syn-rift sandchannel reservoirs. In addition, the SSF method does not locate correctly the base salt, which is almost flat in the model (Fig. 10) but appeared discontinuous and strongly bended in the migrated seismic section (Fig. 13).
The PSPI and FFD sections are affected by artifacts slightly masking some features below the vertical salt bodies. This masking is stronger on the SSF section.

Model 3
The third model is a representation of the Kuqa foreland fold-and-thrust belt in the Tarim basin, where the salt layer is overlaying a series of imbricated thrust faults, pop-up and duplex structures. The model was constructed based on Xiuxiang et al., 2000 (Fig. 19A). In this basin, there are many trap structures formed in the subsalt sequences because of the intense compression of the Kuqa foreland fold belt. The thrust faults served as good routes for hydrocarbon migration, and the dense halite layers acted as good regional seals. The oil reservoirs trapped on top of the folds are in a depth of 1800 m to 1900 m. Some channel sands reservoir were included in the suprasalt layers but they are not an objective of this study.
The migrated sections obtained using the three methods ae shown in Fig. 19. The main reservoir in the section located in the top of a pop-up are presented in Fig. 15, where the dashed green line represents the oil-water contact. PSPI and FFD methods provide acceptable results in general, but they are affected by multiples below the salt nappe, like the arrow 1 shows. The oil reservoir can be identified by conformity with the local structures in PSPI and FFD sections ( Fig. 16 and Fig. 17). SSF method has the lower noise-to-signal ratio and is seriously affected by multiples below the salt nappe (arrows 1 and 2), which mask the other reflections. These seismic multiples also affect the structures into and above the pop-up located in the middle of the model causing a confusion about the real geometry of the layers (arrow 2 in Fig. 18).

Computing time
Computation times for all three methods are shown in Fig. 20, where it is clear that image quality and computing time are inversely proportional to each other. The PSPI algorithm is about seven times more expensive than SSF and about twice than FFD. The FFD algorithm is around three times more expensive than SSF.

Discussion
All the evaluated migration methods showed the variation in the accuracy of reflectors positioning, as well as the effectiveness of collapsing the energy of the reflectors into a single point. The OWWE migration methods have an intrinsic limitation due to the use of wave field propagation algorithms, which implies the knowledge of very accurate velocity models to reach the reliable results. In this work, the velocity models used for migration were the same as used to acquire the synthetic seismograms. Therefore, the migration methods were correctly evaluated. The analysis of the results showed that SSF method has serious problems in the presence of strong lateral velocity variations, and especially when layers have a dip angle greater than 45 degrees (e.g., model 2). The PSPI and FFD methods provided the more quality images, accurately positioning the subsalt layers, and showing important feasibility to be used in the improvement of seismic imaging in basins with the presence of salt bodies or other highvelocity structures (e.g. igneous structures).

Conclusions
The performance of OWWE methods is seriously affected by the steeply dipping structures generated by the salt movement. It is very difficult, maybe impossible, for these methods to eliminate the unwanted artifacts generated by the vertical salt structures and the subsalt multiples generated by the salt-sediment impedance contrast. Although it was possible to identify almost all reservoirs, the interpreters should consider that in any setting with salt bodies present, the seismic artifacts and multiples could give rise to misleading images. However, the analysis of the obtained seismic sections showed the feasibility of using the OWWE methods in subsalt provinces, as long as the interpreters take into account these issues during the prospects identification.    It was established that PSPI provides the most accurate results. Despite its high computational cost, this method recovers most of the events and structures overlaid by salt layers. Certainly, using a more accurate velocity model increases the chances of imaging the reservoir targets. The FFD method also provides good results, with the best ratio quality/computational time. The SSF method was the fastest, but also less accurate than the other methods.