Stretching correction for amplitude-preserving vector wavefield reverse-time migration

Jia-jia Yang , Bing-shou He, Hua-ning Xu , Jun Pan , Jun Liu , Hong Liu

a Qingdao Institute of Marine Geology, China Geological Survey, Ministry of Natural Resources, Qingdao 266071, China

b Laboratory for Marine Mineral Resources, Qingdao National Laboratory for Marine Science and Technology, Qingdao 266071, China

c Key Laboratory of Submarine Geosciences and Prospecting Techniques, Ministry of Education, Ocean University of China, Qingdao 266100, China

Keywords:

Reverse-time migration

Multi-wave and multi-component

Migration stretch correction

Amplitude-preserving imaging

Seismic migration imaging

A B S T R A C T

The migration of multi-wave seismic data is aimed at obtaining the P- and S-wave imaging results of the amplitude preserving. But the P- and S-wave stretching effect produced by the reverse time migration of the elastic wave equation will not only reduce the vertical resolution of the migration results and the amplitude preserving of the large reflection angle. In this paper, the reverse time migration technique of amplitude preserving vector wave-field separating is used. Based on the analysis of the stretch mechanism and the influencing factors of stretch magnitude, the paper gave the stretch correcting factors. Then, realize the stretch correction method at the time that after the reverse extrapolation and before the imaging by solving the problem which is how to calculate the P-wave and Ps-wave propagation directions of imaging points at different times. The stretch correction method can improve the vertical resolution and amplitude fidelity of the imaging results and provide high fidelity input data for seismic data interpretation and inversion.

1. Introduction

Prestack depth migration is an important means of accurate seismic imaging (Chang and McMechan, 1986,1987, 1990; Dong and McMechan, 1993; Sun R and McMechan, 2006; He BH et al., 2010), the significance of prestack depth migration imaging is: (1) accurate imaging of underground structures requires high-precision prestack migration technology; (2) the characteristics of seismic wave mechanics contained in multiwave migration profiles can be used to identify lithology, fluid and microstructural characteristics of rocks (cracks, etc.); (3) angle domain common imaging gathers (ADCIGs) from prestack depth migration can provide input data for prestack inversion. At the same time, seismic data migration should achieve two targets:(1) accurate imaging of P- and S-waves and precise characterization of stratum trend and geological body shape;(2) acquisition of fidelity P- and S-wave prestack inversion gathers and migration stack profiles. At present, the first target of multi-wave seismic data migration has been basically achieved (Liu SW et al., 2013; Clapp RG, 2008; Guo P et al.,2011; Chen K and Wu GZ, 2012; Ding L and Liu SW, 2012,Chen K and He BS, 2014), the second target has been slow to develop. The main reason is a misunderstanding in cognition.It is generally believed that the reverse-time migration of the elastic wave is based on the two-way wave equation, so the amplitude is fidelity after migration. This is a misunderstanding, because in the existing elastic wave reverse-time migration algorithms, no matter what imaging conditions are used, the dynamic characteristics of seismic waves after migration can’t accurately and fully reflect lithology and fluid-bearing properties.

The reverse-time migration based on the elastic wave equation can theoretically achieve amplitude-preserving imaging of P- and S-waves, but the P- and S-wave fields are included in the horizontal and vertical components of the elastic wave equation. It is necessary to decouple P- and S-wavefields during elastic wavefield extrapolation to achieve independent imaging of the P- and S-waves (Sun R et al.,2004). The traditional method is to obtain divergence and curl of horizontal and vertical components to separate P- and S-waves (Sun R, 1999), but this method will lead to distortion of amplitude and phase (Sun R et al., 2001, 2004, 2011),which violates the goal of amplitude preservation reverse-time migration. At the same time, migration stretching correction is also an important part of the research on the fidelity of reverse-time migration. This is because the stretching effect of reverse-time migration will lead to the inconsistency on the spectrum information of seismic data with different reflection angles at the same reflection point, which directly reduces the amplitude fidelity of large reflection angle data.

As for the decoupling of the elastic wavefield, since P-wave is a non-rotational field and S-wave is a non-divergent field in isotropic media, the separation of P- and S-waves can be achieved by divergence and curl calculation of elastic wavefield (Aki K and Richards PG, 1980) in the process of wavefield extrapolation (Sun R, 1999). The divergence and curl of the elastic wavefield is actually the spatial derivatives of the elastic wavefield, which will lead to the amplitude distortion and phase shift of 90° (Sun R et al., 2001, 2004,2011). Sun R et al. (2001) made Hilbert transform to correct the phase shift of the separated wavefield in view of the phase change after wavefield separation. Sun R et al. (2011)proposed a method to correct the amplitude ratio of P- and S-waves by using P- and S-wave velocities on the basis of Sun R (1999) wavefield separation algorithm. Then, based on the view of Sun R et al. (2011), Li ZY (2013) proposed a method for accurately estimating near-surface velocity values, in order to completely separate P- and S-waves from a corrected amplitude ratio of P-and S-waves, at the same time, he corrected the phase in the frequency domain. Du QZ et al.(2015) proposed separating P- and S-wavefield snapshots in the wavenumber domain by divergence and curl operators.When cross-correlation imaging conditions are used, phase correction is realized by multiplying the imaging amplitude by -1. Wang WL et al. (2015) proposed a reverse-time migration algorithm for vector elastic wave. In the process of elastic wave propagation, the P-wave pressure, the velocity equations of P- and S-waves are added to obtain the decomposed vector P- and S-waves. The imaging conditions are applied to the modes of vector P- and S-waves. The polarity is based on the relationship between the direction of P- and S-wave vibrations and the reflection layer, and the polarity of the imaging results (especially if the reflection is greater than the critical angle) is uncertain. Yang JJ et al.(2016a) based on Wang WL et al. (2015) decoupling of P-and S-wave vectors, proposed a method of scaling the vector P- and S-waves and then applying imaging conditions to avoid the problem of uncertainty caused by calculating the polarity of imaging results.

Migration stretching correction is to suppress the low frequency noise interference in migration results with the increase of offset, which is common in prestack depth migration. Tygel M et al. (1994) studied this effect and obtained the conclusion that migration stretching is affected by reflection angle of seismic wave, dip angle of stratum and velocity of medium, and gave the mathematical expression of stretching factor. Levin SA (1998) takes inclined medium as its main target, and focuses on the influence of the reflection angle of seismic waves on migration stretching. Because migration stretching and normal move out correction (NMO)stretching have similar characteristics in performance, that is,with the increase of offset, the frequency moves to the low direction. Many scholars use this characteristic to eliminate the frequency imbalance on common offset gathers or common imaging gathers. Lazaratos S et al. (2004) proposed a method of frequency compensation for stretched wavelet in frequency domain based on unstretched wavelet to eliminate the effect of NMO stretching. Based on the assumption of horizontal layered media, Roy B et al. (2005) designed a deconvolution filter to correct the stretch in the common offset gathers by using the stationary wavelet with a certain incidence angle. Perez extended the former algorithm to inclined media. Fixed spectral morphology operator was applied to data with different incident angles to realize the wavelet stretching correction in ADCIGs (Perez G and Marfurt K, 2006). Canning A et al. (2008) assumes that the zero offset seismic wavelet is known, and proposes a method to realize the migration stretching correction in the angle domain gathers. All the above methods are proposed based on the characteristics of stretching noise, however the pre-stack depth migration stretching is produced when the seismic wave in the time domain is compared to that in the depth domain.For pre-stack reverse-time depth migration, the noise is generated when the imaging condition is applied. Therefore,Zhu XF and McMechan GA (2013) proposed a method of migration stretching correction after wavefield extrapolation and before imaging. The advantage of this method is to choose a better time for stretching correction, that is, between wavefield extrapolation and imaging, this is the time that migration stretching is producing. But his calculation method of reflection angles and formation dip angles is complex,which requires huge calculation and storage cost. On the basis of Zhu XF and McMechan GA (2013), Yang JJ et al. (2016a)extended the imaging conditions to the cross-correlation category, simplified the calculation method of the stretching factor parameters, and realized the stretching correction of the reverse-time migration of the acoustic equation by increasing a small amount of calculation. Yang SW et al. (2016) applied this method to the elastic wave equation and realized the migration stretching correction of P- and S-waves preliminarily.

However, in order to reach the industrial application level,the above achievements need to be further improved,manifested in: (1) correct the calculation method of the reflection angle. In the above algorithms, the reflection angle of the seismic wave is calculated by wavefield propagation direction and stratum dip angle, which is feasible, but there is an error in the algorithm which only uses the travel time of direct waves to calculate the direction of seismic wave propagation, which affects the accuracy of the calculation of the reflection angle; (2) Zhu XF and McMechan GA (2013)and Yang JJ et al. (2016a) can only be applied to reverse-time migration of the acoustic wave equation, meanwhile, although Yang SW et al. (2016) extended the algorithm to the elastic wave equation, he did not give a method for calculating the S-wave reflection angle. If the S-wavefield is corrected by the P-wave reflection angle, the problem of over-correction will be introduced; (3) traditional reverse-time migration of elastic wave equation uses divergence and the curl operator to separate P- and S-waves and then image them, which results in distortion of amplitude and phase of P- and S-waves.Therefore, the same problem will arise in the stretch correction of P- and S-waves based on traditional reverse time migration of elastic wave equation.

Aiming at the realization of amplitude-preserving elastic wave reverse-time migration, this paper uses the reverse-time migration algorithm of amplitude-preserving wavefield separation (Yang JJ et al., 2016b), then eliminate the stretching effect of prestack depth migration and further improve the amplitude-preserving imaging effect of elastic wave reverse-time migration. Firstly, the stretching mechanism and influencing factors of reverse-time migration are analyzed, and a method of calculating the stretching factor of each grid point in space by using the elastic wave equation of wavefield decoupling is proposed. The stretching correction of imaging wavelets in time domain after wavefield extrapolation is carried out. Then, the amplitude-preserving P-and S-wave imaging results after stretching correction are obtained by applying imaging conditions. After correction, the vertical resolution of seismic imaging results with a large reflection angle is effectively improved, and the amplitude and phase characteristics of imaging results remain unchanged, which provides high fidelity input data for seismic data interpretation and inversion and precise exploration of special targets.

2. Generation mechanism of reverse-time migration stretching

Because seismic records are acquired at equal time intervals and seismic waves have different propagation velocities in different media, multiplying equal time intervals with different propagation velocities will lead to distortion in depth domain imaging of seismic waves, which is a common problem in pre-stack depth migration. Tygel et al. has made a thorough study on this problem and deduced an analytical formula for the local stretching factor of prestack depth migration in theory (Tygel M et al., 1994), as follows:

Which,m0is stretching factor, which is defined asΔtis shift in time, Δzis corresponding shift in depth,vis migration velocity, θ is reflection angle, β is the local dip angle of the layer. The meaning of the stretching factor is the amount of variation with a reflection angle, dip angle and migration velocity when time interval is compared to the corresponding space interval. Therefore, pre-stack depth migration stretching is not the propagation effect of seismic waves. For reverse-time depth migration, the stretching is caused at the time between wavefield extrapolation and imaging, that is, it can be corrected before imaging condition is applied.

Tygel’s pre-stack depth migration stretching formula(Tygel M et al., 1994) is extended to elastic wave migration according to the principle of this formula. The analytical formulas of local stretching factors of P- and S-waves with elastic wave pre-stack depth migration are as follows:

Which,mP0andmS0are stretching factors of P- and S-waves, respectively, which are defined asandrespectively; Δtis shift in time; ΔzPand ΔzSare corresponding shifts of P- and S-waves in depth, respectively;vPandvSare migration velocity of P-and S-waves,respectively; θPand θSare the reflection angle of P- and S-waves, respectively; β is the local dip angle of the layer.

Formula (2a) shows that the stretching factor of the P-wave increases with the increase of the reflection angle and the local dip angle of the P-wave, and increases with the decrease of migration velocity of the P-wave. The same conclusion can be obtained by analyzing the S-wave stretching factor of formula (2b). For a particular reflector, the migration velocities above and below can be regarded as constants respectively. The influence of migration velocities on different offsets of a particular reflector is the same. It will not cause the frequency changes of different offsets of the reflector, but will only cause certain changes in imaging width between different reflectors. In this paper, the effect of different offset stretching on the same reflector is studied and corrected, so the effect of migration velocity on P- and S-wave migration stretching is not considered for the time being.

Formula (2) is an analytical formula for P- and S-wave migration stretching, which can be used to correct P- and S-waves migration stretching. The following is the explanation by taking cross-correlation imaging conditions as an example,in conjunction with Fig.1. Fig. 1a-c is the imaging result of a P-wave migration with stretch. From this section, it can be seen that the imaging width increases with the increase of offset. The source wavefield wavelet (blue line in Fig. 1b) and the receiver wavefield wavelet (red line in Fig. 1b) at the depth of 0.5 km are extracted. It can be seen that the source wavefield and the receiver wavefield coincide well in the middle of reflection layer. Extraction of the source wavefield wavelet (blue line in Fig. 1c) and the receiver wavefield wavelet (red line in Fig. 1c) at the depth of 0.54 km shows that the source wavefield does not coincide with the receiver wavefield well. Fig. 1d-f is the imaging result after the stretching correction of the P-wave migration. From this section, it can be seen that the imaging width does not change with the increase of offset. Source wavefield wavelet (blue line in Fig. 1e and Fig. 1f) and receiver wavefield wavelet(red line in Fig. 1e and Fig. 1f) at depths of 0.5 km and 0.54 km are also extracted. They are all corrected by migration stretching. Compared with Fig.1c and Fig. 1f, it can be seen that after migration stretching correction, the P-wave wavelet of the source wavefield and the receiver wavefield are reduced, and the overlapping part of the two fields is reduced,which leads to the reduction of the imaging width (compared with Fig.1a and Fig. 1d). Similarly, the same conclusion can be obtained by comparing the results of S-wave imaging in Fig. 1g-l.

Fig. 1. P-wave (a-f) and S-wave (g-l) RTM images before and after migration stretching correction. a, b, c-P-wave RTM image before migration stretching correction and the extracted imaging wavelet; d, e, f-P-wave RTM image after migration stretching correction and the extracted imaging wavelet; g, h, i-S-wave RTM image before migration stretching correction and the extracted imaging wavelet; j, k, l-S-wave RTM image after migration stretching correction and the extracted imaging wavelet.

Fig.1 illustrates the effectiveness of migration stretching correction using Formula (2). How to realize migration stretching correction is described below.

3. Realization of migration stretching correction

Based on the elastic wave RTM algorithm with amplitudepreserved (Yang JJ et al., 2016b), this paper illustrates the realization of stretching correction for P- and S-waves migration by taking cross-correlation imaging conditions as an example. The specific steps are as follows:

(i) The excitation time of the downward P-wave is calculated and recorded by the forward extrapolation of the elastic wavefield. In this paper, the maximum amplitude time of the downward wave is taken as the excitation time.

(ii) When reconstructing the source wavefield and the receiver wavefield, extracting the propagate waveletWAVSp(x,z,t)of the source P-wave field in the time domain at the time ofT(x,z) (the time window is centered on excitation timeT(x,z), including (n+1) time samples, n is greater than or equal to the length of the source wavelet), using a decoupled vector at excitation timeT(x,z) calculate and record the propagation angle α (x,z) of the source P-wavefield by measuring the P-wave (see Appendix I); extracting the propagated waveletWAVRp(x,z,t) andWAVRs(x,z,t) of the receiver wavefield in the time domain during the receiver wave field extrapolation, and calculating the angle between P-and S-wave by using the Poynting vectors of the decoupled vector P- and S-wave at the time of excitationT(x,z).

(iii) P-wave migration profileIP(x,z) obtained by crosscorrelation imaging conditions.

(iv) Repeat steps (i)-(iii) for each shot, then stack all shots to obtain the P-wave migration stack profileISP(x,z). The signal-to-noise ratio of the P-wave migration stack profile is generally higher than that of the S-wave migration stack profile, which can be used to calculate the local dip angle.

(v) Calculate the instantaneous phase gradient of the complex signal ofISP(x,z) to obtain the local dip angleβ(x,z)( Appendix I).

(vi) The reflection angle of the P-wave (θP(x,z)) is the difference of the propagation angle (α (x,z)) of every shot from step (ii) and the local dip angle (β (x,z)) from step (v).The reflection angle of the S-wave (θS(x,z)) is the difference of the propagation angle (α (x,z)) of every shot from step (ii)and the intersection angle of the receiver P- and S-waves(ϕ (x,z)).

(vii) Compress the waveletWAVSp(x,z,t),WAVRp(x,z,t)andWAVRs(x,z,t) to obtain a new waveletWAVSSSp(x,z,t),WAVSSRp(x,z,t)andWAVSSRs(x,z,t). Compress principle is according to the stretching factor, compressing the wavelet length from (n+1) to (n+1)·cosθP·cosβ and (n+1)·cosθS·cosβ, and keeping the amplitude and phase characteristics unchanged.

(viii) Apply cross-correlation imaging conditions toWAVSSSP(x,z,t),WAVSSRP(x,z,t) andWAVSSRS(x,z,t),obtain the migration results after migration stretching correctionISSP(x,z) andISSS(x,z).

Because the cross-correlation imaging condition is to calculate the cross-correlation between the source wavefield and the receiver wavefield, the propagation wavelet of the source wavefield and the receiver wavefield must be compressed simultaneously. The excitation time imaging condition only uses the receiver wavefield for imaging, so only receiver wavelet compression processing is needed.

4. Numeral experiences

4.1. Inclined layered medium model

An inclined layered medium model (shown in Fig. 2) is used to verify the algorithm. The elastic parameters of each layer are shown in Table. 1. The size of the model is 5.0 km ×2.5 km, the spatial interval is 5.0 m, and the time interval is 0.5 ms, with 6001 samples per trace. The source of seismic record is Ricker wavelet, the main frequency is 30 Hz; the ground shot interval is 25 m; the ground receiver interval is 5 m.

Table 1. Elastic parameters for inclined layered medium.

Fig. 2. Inclined layered medium model.

Fig. 3 shows the P-wave and converted S-wave migration profiles before and after migration stretching correction with the shot located in (2.5 km, 0 km). It can be seen from Fig. 3a that before the migration stretching correction, the vertical width of the reflection layer on the P-wave profile increases significantly with the increase of offset, while the imaging width of the P-wave (Fig. 3c) after correction does not change with offset; compared with Fig. 3a and Fig. 3b, the S-wave migration stretching noise is significantly less than the P-wave, because at the same imaging point, the reflection angle of the converted S-wave is less than that of the P-wave, so the S-wave migration stretching strength is smaller.

Fig. 3. P- and S-wave RTM images before and after migration stretching correction with the source located in (2.5 km, 0 km). a-P-wave RTM image before migration stretching correction; b-S-wave RTM image before migration stretching correction; c-P-wave RTM image after migration stretching correction; d-S-wave RTM image after migration stretching correction.

The P- and S-wave stacked profiles before and after migration stretcing correction are processed with the same mute showed in Fig. 4. Compared with Fig. 4a and 4c, it can be seen that P-wave has obvious low-frequency noise before migration stretching correction (Fig. 4a). This is because migration stretching makes the imaging result move to the low-frequency direction. The migration stretching correction(Fig. 4c) corrects the frequency of the imaging result,improves the vertical resolution of the result and suppresses the low-frequency noise caused by migration stretching.Compared with the stacked profiles of converted S-wave migration in Fig. 4b and Fig. 4d, it can be seen that the converted S-wave is less affected by migration stretching noise, but the imaging width of each reflector after stretching correction decreases, which improves the vertical resolution of the profiles.

Fig. 4. Stacked P- and S-wave RTM images before and after migration stretching correction. a-Stacked P-wave RTM image before migration stretching correction; b-stacked S-wave RTM image before migration stretching correction; c-stacked P-wave RTM image after migration stretching correction; d-stacked S-wave RTM image after migration stretching correction.

Fig. 5 extracts ADCIGs of P-waves and converted S-waves before and after migration stretching correction, which gives the gathers at 1.25 km, 1.875 km, 2.5 km, 3.125 km and 3.75 km, respectively with the angle ranges from -60° to 60°.It can be seen that the events of ADCIGs before the stretching correction of P-wave migration (Fig. 5a) widening with the increase of the incident angle, and the obvious stretching effect occurs. After the stretching correction of the P-wave migration (Fig. 5b), the events of the large angle get convergence, and the events continuity of different incident angles are the same. The angle gathers before and after the stretching correction of S-wave migration in Fig. 5c and Fig.6d are the same. From Fig. 5a and Fig. 6c, it can be seen that the stretching effect of the P-wave is stronger than that of the S-wave at the same incident angle. This is because the same incident angle of the P-wave corresponds to different reflection angles of the P-wave and the S-wave. Meanwhile,the larger reflection angle, the more obvious the migration stretching is. Therefore, the stretching magnitude of the S-wave at the same imaging point is smaller than that of the P-wave.

Fig. 5. P- and S-wave ADCIGs before and after migration stretching correction with angle range -60°-60°. a-P-wave ADCIGs before migration stretching correction; b-P-wave ADCIGs after migration stretching correction; c-S-wave ADCIGs before migration stretching correction;d-S-wave ADCIGs after migration stretching correction.

Fig. 6. The stacked P-wave RTM images before and after migration stretching correction. a-Stacked P-wave RTM images before migration stretching correction; b-stacked P-wave RTM images after migration stretching correction.

4.2. Measured data of area

The test of the measured data in a certain area is carried out. Firstly, simple post-processing is implemented to the single shot migration data which is before and after stretching correction; then, extract the common imaging gathers to obtain the stacked profiles. The P-wave migration stacked profiles before and after the migration stretch correction are shown in Fig.6. It can be seen that the signal-to-noise ratio of the profile after migration stretching is better than that of the profile before migration stretching. Because the missing source is located in the horizontal position of 2540-3160 m and 6040-6650 m of the original data (the yellow line in Fig.6a-b). The original stacked profiles have problems in that the events are discontinuous and the events’ energy is unbalanced. Compared with the original profile, the formation information with the large reflection angle can be retained after migration stretching correction. Therefore, in theory,structural imaging can be obtained in the case of sparse shot points.

Blocks A, B and C of Fig.6a-b are partially magnified as shown in Fig. 7. Block A in Fig. 6a-b is shown in Fig. 7a(block A in Fig. 6a) and 7b (block A in Fig. 6b), and Block B in Fig. 6a and Fig. 6b is shown in Fig. 7c (block B in Fig. 6a)and 7d (block B in Fig. 6b). By comparison, we can see that the results of migration stretching correction have higher resolution. Moreover, the stacked profile events before migration stretching correction is discontinuous in the area where the shot point is missing, and the events continuity is effectively improved after migration stretching correction.Block C in Fig. 6a-b is shown in Fig.7e (block C in Fig. 6a)and 7f (block C in Fig. 6b). The contrast shows that migration stretch correction improves the resolution of imaging results.

Fig. 7. Local enlargement of block A, B and C of Fig. 6. a-Local enlargement of block A of Fig. 6a; b-local enlargement of block A of Fig.6b; c-local enlargement of block B of Fig. 6a; d-local enlargement of block B of Fig. 6b; e-local enlargement of block C of Fig. 6a; f-local enlargement of block C of Fig. 6b.

5. Conclusions

In this paper, based on the elastic wave equation with amplitude-preserving separation of P- and S-waves, the generation mechanism and influencing factors of migration stretching are analyzed, and the stretching correction algorithm of P- and S-wave reverse-time migration is studied.This method guarantees the fidelity of amplitude of P- and S-wave imaging results calculated by multi-component seismic wave reverse-time migration which is based on elastic theory.The processing of amplitude fidelity includes two aspects:Firstly, using the elastic wave equation of amplitudepreserving separation of P- and S-waves, the amplitude and phase of the separated P- and S-waves are kept unchanged in the process of wavefield extrapolation. Secondly, the time of migration stretching correction is between wavefield extrapolation and imaging. It guarantees the consistency of the frequency spectrum of seismic data with different reflection angles at the same reflection point after migration processing, and the amplitude fidelity of large reflection angle data. In additional, it also provide high fidelity input data for seismic data interpretation and inversion.

Acknowledgement

The article is financially supported by Qingdao National Laboratory for Marine Science and Technology(QNLM2016ORP0206), National Science and Technology Major Project (2016ZX05027-002), National Key R&D Plan(2017YFC0306706-04, 2017YFC0307400) Qingdao National Laboratory for Marine Science and Technology(QNLM201708), China Postdoctoral Science Foundation (No.2017M612219), and Natural Science Foundation of Shandong Province (No. ZR2016DB10). The anonymous reviewers provided valuable suggestions and comments on the manuscript, which is greatly appreciated.

Appendix I