The Offset-Domain Prestack Depth Migration with Optimal Separable Approximation

INTRODUCTION

Complicated geological structures are characterized by large dipping angles, faults, irregularities in a small range, large burial depths, and intense lateral variations of velocity, and the routine processing of seismic data cannot achieve clear information of subsurface media. In recent years, the prestack depth migration has been a critical technique to image complex subsurface structures such as thrust faults and deep ancient-buried hills. In the principle, this method can be classified into Kirchhoff integration prestack depth migration and wave-equation-based prestack depth migration. Although the Kirchhoff prestack depth migration requires less computation, it is suitable for various geometry and is easy to construct migration velocity fields; some problems are not treated well, such as multi-traveltime and amplitude, scattered points, and blind areas. Considerable efforts were put to improve this technique, but the results are not satisfying.

Based on the wave equation, the prestack depth migration of the wavefield extrapolation has the advantages of complete theory, high precision of imaging, amplitude preservation, and being applicable to complex structures. This method has been focused on imaging in complicated media. There are two kinds of approaches to achieve the wavefield extrapolation operator. One is the finite difference (Ristow and Ruhl, 1994), in which a large number of calculation is required to consider the high-angle imaging and 3D azimuth anisotropy (Zhang et al., 2000; Ristow and Ruhl, 1997). The other is the explicit method by Fourier transform in the frequency-wavenumber domain. The method is exact for homogeneous media (Gazdag, 1978). However, in the presence of lateral velocity variations, making approximation to the operator is required. In recent years, the approximation method has been greatly developed, such as the splitting step phase shift (SSF), generalized screen (GSP), Fourier finite difference (FFD), optimum split-step Fourier (OSP), and optimal separable approximation (OSA). Stoffa et al. (1990) presented the SSF, which decomposes the media slowness (velocity) into background slowness and perturbation component, and generalizes the phase shift method to media of laterally varying velocity based on the small perturbation theory. This method is suitable for media of weakly lateral velocity variations. The GSP algorithm (Le Rousseau and De Hoop, 2001) generalized the split-step Fourier methods. By two Taylor expansions, the accuracy of the operator approximation was greatly improved over the SSF method. However, the method suffers a problem of singular points in the complex plane. The OSP (Liu and Zhang, 2006) achieves optimum approximation to the Taylor expansion coefficients, and the operator accuracy and numerical stability are enhanced. The OSA (Zhang et al., 2005; Chen and Liu, 2004; Song, 2001) uses the product summation of two single-variable complex functions to represent the symbol of the one-way wave operator, so that the operator is separated in the wavenumber domain and the space domain. The phase shift operation is implemented in the wavenumber domain and the time shift caused by the lateral velocity variation is corrected in the space domain. The background velocity is not used to calculate the phase shift. Hence, the wave fields of each velocity value in the media of laterally strong varying velocity can be located and focused very well. The migration results of 2D and 3D seismic data show that this method is capable of imaging small faults and lenses in complex media.

A large amount of computation is required by the wave equation prestack depth migration, especially for migration velocity analysis, which leads to an unacceptable amount. The method to reduce prestack numerical calculation to raise the speed of this method remains to be explored. One approach is the area shot technique that the single shot data are synthesized to records of receivers, which reduces the data amount for migration. The method is implemented by the single square root equation. Another alternative approach is to use the double square root equation on the middle point-offset domain data to achieve migration. Popovici (1996) and Jin et al. (2002) employed the OSP and GSP to implement migration with the double square root equation, respectively. On 3D data, Biondi and Palacharla (1996) completed double square root equation migration in the common azimuth domain. These results show that the offset-domain method of the double square root equation can reduce the computation amount of migration with respect to the shot domain method.

In this article, we begin with the double square root equation under the midpoint-offset coordinates, and apply the optimal separable approximation to approximate the symbol of the operator. The forward and inverse Fourier transforms are adapted to construct the continuation operator of the wave field and to implement imaging. In this way, we calculate the impulse response of the non-zero offset to examine the precision of the operator constructed by the separable approximation method. Finally, we achieve migration of the model data.

CONTINUATION OF THE OFFSET-DOMAIN WAVE FIELD

Double Square Root Operator and Its Symbol

The double square root equation in the frequency-wavenumber domain inhomogeneous media can be expressed as (Claerbout, 1985)

where X_s and X_r are the positions of the shot and the receiver, respectively; ω is the circular frequency; ; v(x_s, z) and v(x_r, z) denote the velocities of the shot and receiver at depth z, respectively; k_s and k_r are wave numbers, and P(k_s, k_r, z, ω) is the wave field in the frequency-wavenumber domain. Equation (1) is the square root equation under the shot-receiver coordinates, which can be transformed into the equation in midpoint-offset coordinates. Let

where m is the midpoint, and h is the half offset. By substituting equation (2) into equation (1) and transforming, the double square root equation is obtained in the midpoint-offset coordinates

When the medium is homogeneous, the wave field continuation operator is a simple operation of phase shift. In this case, the operator consists of Fourier transform, product of phase factors, and inverse Fourier transform, which is implemented by fast Fourier transform. For the laterally heterogeneous medium, the double square root operator of equation (3) is a pseudo-differential operator. However, the continuation operator and the phase shift operator are almost similar. The difference is that the phase factor depends on spatial coordinates, and Fourier transform is performed for various velocity values in space, which leads to a large amount of computation. Based on the theory of propagator (Le Rousseau and De Hoop, 2001), we can provide the symbol of the wave field continuation operator for the double square root equation

where, E = exp[zΔzk_z(k_m, k_h, z)] is called the symbol of the wave field continuation for the double square root equation, and k_z(k_m, k_h, z) is the vertical wave number associated with space coordinates. When the velocity of the media varies laterally, the processing of the wave field continuation operator turns to the symbol approximation. The algorithms SSF, GSP, and FFD use the Taylor expansion to approximate the symbol that is decomposed into the phase shift operation containing the background velocity and the time shift operation to velocity variations, in which the precision of wave field imaging to various velocity values is related to the background velocity.

Optimal Separable Approximation of the Symbol

Based on the ideas of the approximations stated above, a wave field continuation operator of the double square root equation was constructed in the media of lateral velocity variation, and the optimal separable approximation was used to approximate the symbol of the operator. The method has high precision at each velocity value of the media, and is not influenced by the background velocity. Meanwhile, it is still implemented by fast Fourier transform. Equation (4) is written as

where, . Equation (5) is a function with four variables, which is very difficult to solve. Let u = (u_r, u_s), and k = (k_m, k_h). Analogous to the expression of the optimal separable approximation of the symbol for the one-way wave operator, the symbol of the wave field continuation operator can be expressed as (Chen et al., 2007; Chen and Liu, 2004)

For β^(l)(k) = β^(l)(k_m, k_h), the same optimal separable approximation is derived

Substituting equation (7) into equation (6), we obtain

where λ =λ₁λ₂.Since the velocities of the shot and the receiver are the same at the corresponding space positions, the velocity terms in space variables are also the same for wave field continuation. We rewrite equation (8) as

The references of the optimal separable approximation of the symbol for the one-way wave operator show that the symbol approximation is of high precision with small value for s and t.

Figure 1 shows the comparison of the numerical approximations to the operator symbol for the double square root equation using three algorithms (SSF, GSP, and OSA). Since the offset is zero, the offset item function can be omitted, and the symbol is simplified as E(u, k_m). The velocity range of the medium is 1 500-2 500 m/s, km is 0–π/25, and ω is 40π, and we make a uniform sampling to u and k_m (or v and k) by 40 and 100, respectively. As seen in Fig. 1, the phase errors of the three algorithms are all very small at the narrow propagation angle. With increasing propagation angle, the phase errors of SSF and GSP become gradually large, and GSP has a singular point near the propagation angle 90. By the optimal separable approximation of 5, the symbol of the operator can be approximated well at wide propagation angles, except for a small loss of amplitude at the large angle. Therefore, the symbol of the operator can be approximated with optimal separable approximation very well, while the computational precision and efficiency remain good.

Figure  1.  Comparison of the numerical approximations to the symbol using different approximation methods.

DownLoad: Full-Size Img PowerPoint

Wave Field Continuation of Offset Data

The offset data migration of the double square root equation is implemented by wave field continuation. The optimal separable approximation is used to express the symbol of the wave field continuation operator in a thin-layer medium, and by substituting it in equation (3), we obtain

where denote the forward and inverse Fourier transforms, respectively. By continuing the wave field of all frequencies, and then performing summation, the wave fields of the double square root equation are achieved at t = 0 and h = 0.

Impulse Response of the Operator

To test the migration effect with optimal separable approximation in a medium of laterally strong varying velocity, we calculate the impulse response in the case of laterally varying velocity. In this model, there are 200 points in the X direction with sampling interval 10 m and 200 points in the y direction with sampling interval 5 m. The Ricker wavelet of main frequency 15 Hz is placed at 0.4 s as a point impulse, and the time sampling interval is 4 ms.

For a laterally heterogeneous medium with a velocity variation range of 1 500 m/s to 2 500 m/s, we sample to u at interval 40 m. The optimal separable approximation of 5 expansions is used to approximate the symbol of the operator for different velocities, and the wave field continuation operator of the double square root equation is constructed. The calculation results of the impulse responses in the medium of velocity 2 500 m/s using the three methods are shown in Fig. 2. It indicates that the results of both the SSF and OSA methods approximate to the analytic solution very well at small propagation angles. For the medium of laterally varying velocity, these two methods can image the wave field well at small angles. However, the impulse response by the SSF method deviates from the theoretical values at large angles, because there is an error in the approximation of the symbol. The impulse response by the OSA method at large propagation angles coincides well with the theoretical solution, which indicates that the method is better than that of SSF in wave field imaging for large angle phases. Figure 2 also shows that the response impulse of the common midpoint-offset domain is an ellipse, and its ellipticity becomes larger with increasing offset.

Figure  2.  Comparison of the calculated impulse responses. Upper, middle, and lower are responses of theoretical solution, SSF, and OSA, respectively. Left: offset=200 m; Right: offset=400 m.

DownLoad: Full-Size Img PowerPoint

CALCULATION OF MODEL AND FIELD DATA

To further examine the imaging method presented in this article, the SEG/Marmousi model data are used to test the prestack offset-domain depth migration for the double square root equation using the OSA method. The velocity of the model is from IFP based on a geological profile in the Cuanza basin. The single-shot records are simulated by the 2D acoustic finite difference method. As seen in Fig. 3a, there are 737 samples in the space direction and 750 samples in the depth direction, respectively. The space sampling interval is 12.5 m for the velocity field, with a maximum depth of 3 000 m and a sampling interval of 4 m in the depth direction. For the data, the number of time samples is 750 with a sample interval of 4 ms, and a sampling length of 3 000 ms. Totally, there are 240 shots, and 96 traces for each shot. Figure 3b shows the simulated seismic records of the 100th and 101st shots. Since the model is a complex structure with several faults and tilting strata, the obvious diffractive waves are shown in records of single shots, and the events of reflective waves are not of hyperbolic shape and are discontinuous.

Figure  3.  Two single-shot records of the Marmousi model.

DownLoad: Full-Size Img PowerPoint

Using the optimal separable approximation of 5 expansions, the imaging result of the prestack shot and the offset depth migration is shown in Fig. 4. The section exhibits clear images of three distinct faults and roughly convergent breaking points. The deep gas-bearing low-velocity bodies are imaged exactly. For event with small dip angles, the method is well capable of imaging, and the effect of the shot-domain method is better than that of the offset-domain method, since for tilting strata, the true CDP is not completely consistent with that of the extracted offset data. Nevertheless, the offset method can choose any record of offset and objective data for migration and imaging, which can be used for velocity analysis of prestack depth migration and for reducing the computation amount.

Figure  4.  Migration results of the Marmousi model with the OSA method.

DownLoad: Full-Size Img PowerPoint

CONCLUSIONS

The optimal separable approximation (OSA) is a regional approximation method and can approach well to the symbol of the continuation operator for the double square root equation at each velocity. The wave field continuation operator constructed by the OSA method has relatively high precision, and the variables in the wavenumber domain are separated from those in the space (velocity) domain. Phase shift is performed in the wave number domain, and the time delay, which is caused by the lateral velocity variation in the media, is corrected in the space domain. In the midpoint-offset coordinates, the wave fields of shots and receivers are continued downward simultaneously, and seismic data do not need to be migrated shot by shot. Therefore, the method has fairly high efficiency. For 3D data, the common azimuth data sets are generated first and the operator can be constructed in a similar way, of which the computation amount can be reduced greatly (Biondi and Palacharla, 1996).

The impulse responses of the operator are achieved on laterally variable velocity, and are roughly in agreement with the theoretical solution. Especially, the impulse responses can approximate well the theoretical wave fronts at the large propagation angles, which indicates that the offset domain OSA prestack depth migration is suitable for imaging complex structure in media of strongly varying velocity. The 2D Marmousi model was tested, and the result shows that the method can be applied to image complicated structure. The used data are midpoint-offset gathers in the method, and input data can be selected according to imaging objects. Hence, the method can be used to iterate the migration and estimation of the prestack migration velocity for the wave equation.