1 INTRODUCTION

Seismic migration plays a key role in oil and gas exploration to obtain the subsurface structure at present,especially in complex exploration area (Gray et al. 2001; Yilmaz,2001; Etgen et al., 2009; Leveille et al., 2011) . However,with the increasing complexity and dem and for precise imaging of the exploration target,seismic migration result is required not only to satisfy the imaging of complicated structure,but also to show lithological changes. It is hoped to give similar result as full waveform inversion in lithological research,so as to satisfy the lithology interpretation of seismic data,especially the lithology interpretation of composite complex reservoir. Thus it is necessary for seismic migration to have a more solid theoretical foundation.

Originated from early 1970s,the modern seismic migration method initially tried to back propagate the reflected seismic wave observed on the surface to its starting location through mathematical method, and gain the image which can indicate the location of the subsurface reflectivity (Claerbout, 1971,1976; Ma et al., 1989; Robein,2010; Stolt,1986; Bleistein et al., 2000; Li et al., 2011) . After 40 years development,seismic migration has developed from post-stack to pre-stack migration,from two-dimension to three-dimension,from time domain migration to depth domain migration,from isotropic medium to anisotropic medium,from scalar acoustic wave equation to vector elastic wave equation (Du,2012; Chang,1994) ,from conventional migration without iteration to true amplitude migration and iterative least squares migration (Zhang,2013; Yao, 2012,2012; Dong,2012; Gray,1997; Sun et al., 2002; Cui et al., 2004; Liu et al., 2006; Zhang,2006; Li,2008; Zhou et al., 2014) , and it also tried to realize velocity and impedance inversion (Zhang, 2013b) . The calculation of seismic migration can be divided into two broad categories: 1) is based on Kirchhoff integral,2) is based on solving differential equation. The Kirchhoff integral migration includes general ray Kirchhoff integral method (French,1974; Schneider,1978; Sun,2000) and Gaussian beam Kirchhoff method (Červený,1982; Costa,1989; Hill,1990) . Differential equation migration includes one-way wave equation method based on wavefield decomposition (Gazdag,1978; Berkhout,1979; Huang,1996; Ristow,1994; Chen,2006; Chen,2007) and full waveform equation method (Hemon,1978; McMechan,1983; Baysal,1983; Whitmore,1983) .

The propagation of reflected seismic wave has three processes including incident wave propagation,reflection, and reflected wave propagation. It can be formally denoted as WRW concept model proposed by Berkhout (1982) . Currently,differential equation migration methods are developed according to this reflected wave propagation. However,as the lack of specific mathematical physics equation that describes the propagation of reflected wave,those migration methods are not established on the basis of the reflected wave propagation equation. The Claerbout's time consistency imaging principle,which is the basis of the imaging formula of the differential equation,is a conceptual description about the position of subsurface reflection point by the relationship of the incident wave's and reflected wave's travel time. The imaging equation proposed by Claerbout states the relationship between reflected wavefield and incident wavefield (that means the reflected wave on the reflection is equal to the multiplication of reflectivity and incident wave) . But this quantitative relationship can't satisfy reflected wave propagation equation,which is the main error of the Claerbout's imaging method. When applying Claerbout's imaging method,the migration result may have phase error,leading to inaccurate imaging position. As we know,the reflection coefficient is defined as the reflected wave amplitude divided by the incident wave amplitude,not the reflected wavefield divided by the incident wavefield. So we should derive the imaging formula by the relationship between reflected wavefield and incident wavefield,which is established according to the mathematical and physical formula of seismic wave, and correct the Claerbout migration error based on this. Based on the Huygens principle and the Kirchhoff integral formula of diffraction wave propagation,Kircchoff integral migration only considers seismic wave propagation, and does not consider the interaction of incident wave and uneven velocity model during the seismic wave propagation. However,scattering wave and reflected wave are generated by this interaction,so we think that Kirchhoff integral migration can't match with the generation and propagation of scattering wave and reflected wave. By using the scattering theory of seismic wave propagation,quasi differential operator theory,Fourier integral operator theory and applying the Kirchhoff integral formula to reflected wavefield ,Bleistein (2000) proposed the modified Kirchhoff integral migration method-true-amplitude Kirchhoff integral migration method. But he still used Claerbout's imaging formula to express the relation between reflected wave and incident wave. Mathematically,seismic wave migration belongs to the solution of the linear inversion problem,so we should build the migration method according to the forward propagation equation of seismic wave. While in the current seismic migration,instead of taking the difference between scattering wave and reflected wave into consideration,it's commonly to share the same migration method with scattering wave and reflected wave,even directly use the scattering wave propagation equation to derive reflected wave migration equation. That is very inappropriate. In the process of seismic wave propagation,the incident wave interacts with the heterogeneities to produce scattering wave. If the size of the heterogeneity is relatively large compared to the wavelength of the seismic waves,the interaction of incident wave and heterogeneities generates reflected wave. So scattering seismic migration should be established based on scattering wave propagation equation, and reflected seismic migration should be established based on reflected wave propagation equation.

On the basis of the above knowledge and the current problems existing in the seismic migration,this paper use the scattered wave propagation theory as a starting point,firstly set up the scattering seismic migration method according to the linear equation of the primary scattered wave,then under the condition of high frequency approximation,approximate the spatial variation of the velocity perturbation to gain the subsurface reflectivity distribution,thus degrade the scattered wave equation based on perturbation into the reflected wave equation based on reflectivity, and establish the reflected seismic migration method by this linear primary reflected wave equation,finally use simulated seismic data to prove the validity of the scattered and reflected wave migration method proposed in this paper.

2 THEORY

Scalar seismic wave equation is generally used to realize migration,that is using the following wave equation to describe the propagation of seismic wave,

where p (x_s,x,t) is the pressure wavefield,v (x) is the velocity model; f (t) is the source wavelet; x = (x,y,z) is the scattering point,x_s is the source location; ▽² denotes Laplacian operator with the form ▽² =$\frac{{{\partial }^{2}}}{\partial {{x}^{2}}}+\frac{{{\partial }^{2}}}{\partial {{y}^{2}}}+\frac{{{\partial }^{2}}}{\partial {{z}^{2}}}$. Given velocity model v (x) ,Eq. (1) describes the nonlinear propagation of the seismic data,that means the wavefield p (x_s,x,t) has nonlinear relationship with v (x) . The Green function g (x_s,x,t) of Eq. (1) can be described as Transfer Eq. (1) and (2) into frequency domain,respectively, and obtain In order to gain the corresponding linear forward modeling of seismic migration,perturbation method is used to establish the seismic wave propagation equation.

2.1 Forward Propagation Equation of Primary Scattered Wave

In order to depict the effect arising from velocity inhomogeneity on the wave propagation,perturbation method can be used to denote velocity model. That means:

where c (x) is the background velocity model,α (x) is the heterogeneous distributed velocity model,also can be named as perturbation model. Substitute this expression into Eq. (3) to obtain the scattered wave propagation equation, where P_s (x, (ω) ) is the scattered wavefield produced by perturbation; P (x, (ω) ) is the full wavefield,P (x, (ω) ) = P_i (x, (ω) ) + P_s (x, (ω) ) ,here P_s (x, (ω) ) is the background incident wavefield,it can be expressed as Transfer Eq. (7) to the time domain, The full wavefield on the right side of Eq. (6) can be expressed as the summation of incident wavefield and scattered wavefield, The right side of Eq. (8) is the secondary source term that can excite scattered wavefield,which is composed of the interaction between incident wavefield and perturbation and the interaction between scattered wavefield and perturbation,so Eq. (8) also describes the nonlinear propagation of scattered wave. That means scattered wavefield P_s (x, (ω) ) has nonlinear relationship with α (x) . When P_s (x, (ω) ) $\ll $ P_i (x, (ω) ) ,under the Born approximation,we can neglect the interaction between scattered wavefield and perturbation,then Eq. (8) can be rewritten as Eq. (9) also can be transferred into time domain, Eq. (9) expresses a linear seismic propagation,incident wavefield interact with perturbation to produce primary scattered wavefield,so Eq. (9) can be named as linear primary scattered wave propagation equation. In Eq. (9) ,P_i (x, (ω) ) has the linear relationship with both P_i (x, (ω) ) and α (x) .

Given background velocity model c (x) ,velocity perturbation α (x) and source wavelet f (t) ,we can simulate scattered wave in the time domain with Eqs. (7-1) and (9-1) ,

where x_g is the receiver's location,d_s (x_s,x_g,t) is the primary scattered data recorded on the receivers. From Eq. (10) ,the simulation of primary scattered wave includes the propagation of incident wavefield,the second source due to the interaction of the second order temporal derivative of the incident wavefield with the perturbation and thus generated scattered wave, and the propagation of the scattered wave. According to the linearity of the wave equation,Eq. (10) also can be expressed as

2.2 The Migration Theory of Scattered Seismic Data

Under the Born approximation,wave Eq. (9) , (9-1) not only can be used to simulate the primary scattered wave,but also is the foundation of scattered seismic inversion and migration. With the Green function of the background velocity model,incident wavefield and primary scattered wavefield in the frequency domain can respectively be expressed as

let then Eq. (15) can be expressed in the matrix form as where u is the matrix of the scattered seismic data D_s (x_s,x_g, (ω) ) ,W is the propagated matrix with G (x_g,x, (ω) ) as integral kenel; m is the unknown matrix. The least square solution of the matrix Eq. (16) is where W^-g is the least square generalized inverse matrix of W,W^* is the adjoint maxtrix of W.

Put themobtained from Eq. (17) and the incident wavefield P_i (x_s,x, (ω) ) obtained from Eq. (12) into Eq. (14) to gain the evaluation of velocity perturbation α (x) ,

where Re means taking the real part.

We replace the least square generalized inverse matrix W^-g with its adjoint matrix W^* to ensure the stability of numerical computation and reduce computational complexity in Eqs. (17) and (18) ,

Virtually,the calculation of Eq. (19) is equal to the adjoint propagation of the wavefield u,then Eq. (18) can be rewritten as where P_ic (x_s,x, (ω) ) is the complex conjugate of P_i (x_s,x, (ω) ) .

Eqs. (12) , (19) and (20) together constitute the approximate solution of velocity perturbation α (x) ,which in essence is the migration image of α (x) . Thus in time domain,the migration calculation of α (x) is,

1) Forward propagation of the incident wavefield

2) Adjoint or reverse time reconstruction of the scattered wavefield

3) The image of the velocity perturbation

According to the simulation of primary scattered wave in Eq. (11) ,the calculation of Eq. (21) can be modified to the following equation to calculate the second time derivative of incident wavefield,

And the corresponding imaging of velocity perturbation is

The above Eqs. (21) or (24) , (22) , (23) or (25) constitute the primary scattered wave migration method proposed in this paper.

2.3 Forward Propagation Equation of Primary Reflected Wave

High frequency approximation method has been widely applied in seismic research. When the seismic wavelength is relatively small compared with velocity perturbation,high frequency approximation can be used to approximate the wave phenomenon in seismic propagation and the spatial variation of the velocity perturbation. Under the high frequency approximation,geometrical optical theory can approximately replace wave theory and then the scattered wave generated by perturbation degrades into reflection and refraction. What's more,the spatial variation of the perturbation also can be viewed as relatively slow-varying or even constant compared with seismic wavelength.

Base on such recognition and under the high frequency approximation,defining the wavenumber of background velocity model along the incident propagation direction n,k =$\frac{\omega }{c\left (x \right) }$ and defining the reflectivity as the derivation of the perturbation α (x) along the incident propagation direction in the wavenumber domain,

Then under the high frequency approximation,scattered wave propagation Eq. (6) degrades into reflected wave equation, The corresponding primary scattered wave propagation Eq. (9) degrades into primary reflected wave equation, This equation describes the linear propagation of the reflected wave P_r (x, (ω) ) ,which has a linear relationship with both incident wavefield P_i (x, (ω) ) and reflectivity r (x) . Its time domain form is Eqs. (28) and (28-1) respectively are the primary reflected wave propagation equation proposed by this paper in the frequency domain and the time domain. Reflectivity r (x) also can be defined as the spatial derivation of the perturbation α (x) along the incident direction n, Reflectivity r (x) defined by Eq. (26) or (29) is a distribution function along the velocity interface. It has a reciprocal length dimension.

Given background velocity model c (x) ,reflectivity distribution function r (x) and source function f (t) ,primary reflected wave simulation can be done with Eqs. (7) and (28) as follows

where x_g is the receiver location,d_r (x_s,x_g,t) is the reflected wave data recorded on the receivers. From Eq. (30) ,primary reflected wave simulation includes three processes: the propagation of incident wave,the interaction of the first order time derivative of the incident wave and the reflectivity and then generation of reflected wave,the propagation of the reflected wave. Similar to Eq. (11) ,Eq. (30) can be rewritten as

2.4 Reflected Seismic Data Migration Theory

The wave Eqs. (28) and (28-1) obtained under the high frequency approximation condition,can not only be used to simulate primary reflected wave,but also is the foundation of seismic inversion and migration. Similar to the scattered wave migration method,we calculate the imaging of the reflectivity in the time domain.

1) The incident wavefield

2) The adjoint reconstruction of the reflected wavefield

3) The imaging of the reflectivity

According to Eq. (31) ,the above first step can be changed to calculate the first derivative of the incident wave versus time,

The corresponding reflectivity imaging is

Eq. (32) or (35) , (33) , (34) or (36) form the primary reflected seismic migration theory proposed in this paper. For comparison,the following is the conventional reserve time migration equation,which is based on the Claerbout's imaging theory.

1) The incident wavefield calculation

2) The adjoint reconstruction of the reflected wavefield

3) The reflectivity imaging

By comparison,the main difference between the reflected data migration proposed in this paper and the conventional migration is the imaging equation,which is the difference between Eqs. (34) and (39) . Compared with conventional migration image,reflected migration image has 90° phase shift due to the first order time derivative terms in the source wavefield. What's more,derivative operation can improve the image's resolution. Thus for reflected wave migration,we suggest to apply the reflected wave migration method proposed by this paper,because Claerbout's migration method not only has 90° phase error,leading to inaccurate imaging position,but also results in lower-resolution image.

Comparing the scattered imaging Eqs. (21) , (22) and (23) with the conventional migration Eqs. (37) , (38) and (39) ,it can be found that the main difference also is the imaging equation. Because imaging Eq. (23) contains the second order time derivation terms,its imaging result has an 180° phase shift or reversed polarity compared with conventional migration. What's more,as derivation can improve resolution,the scattered seismic migration has higher resolution than conventional migration,so we suggest applying the scattered wave migration method proposed in this paper to scattered seismic migration. But if applying this scatter migration method to reflected data,the reflectivity result may have -90° phase error. Conversely,if applying the reflection migration method to scattered data,the scattered result may have 90° phase error.

3 NUMERICAL TESTS

Above theory and formulas are derived in the form of three-dimension,for convenience,we use twodimensional model to verify the accuracy and efficiency of the proposed scattering seismic data and reflected seismic data. For comparison,conventional migration results using Eqs. (37) , (38) and (39) are also given in the experiment,so finally including conventional migration result,scatter migration result and reflection migration result in the model test. As all the seismic simulations and migrations are in two-dimensional form,the test result is reliable.

3.1 Simple Model Test

Figure 1 is the velocity model. This model contains nine different size high-speed and low-speed scatterers and one horizontal reflector,which can be used to verify the effect of scattered wave migration and reflected wave migration respectively. The model is 5 km in length and 2 km in depth,grid spacing is 10 m in both directions. There are 5 high speed scatterers whose size from left to right respectively is three,four,five,six and eight grid points (number 1,3,5,7,9 black scatterer in Fig. 1) . All the high-speed scatterers' velocity is 2500 m·s^-1. Four low-speed white scatterers with velocity 1600 m·s^-1 are arranged left-to-right,whose size respectively is three,four,five and six grid points. The first layer is 2000 m·s^-1, and the second layer is 3000 m·s^-1,the velocity interface is at 150 grids point in depth. 250 shots are fired with a middle shooting and both sides receiving spread geometry. The shot interval is 20 m and each shot is recorded by 401 receivers with a 10 m receiver interval and 0 m minimum offset. A Ricker wavelet with a 20 Hz peak frequency is used as the source wavelet, and the record length is 2.5 s,the time sampling rate is 1 ms.

Applying the conventional migration,scattered wave migration and reflected wave migration method to this model (Fig. 1) ,then he migration result is shown in Figs. 2,3, and 4,respectively.

Comparing Figs. 2,3, and 4 with Fig. 1,it is clear that as for the imaging of scatterers,the imaging quality in Fig. 3 is the best,it not only has high resolution,but also correct phase position. Although the result in Fig. 2 got the right location,it has low resolution and 180° phase error,or reversed polarity,which means that the five high-speed scatterers 1,3,5,7, and 9 in Fig. 1 are imaged as low-speed scatterers and the four low-speed scatterers are imaged as high-speed scatterers. The imaging result in Fig. 4 has 90° phase error,leading to incorrect imaging location. However,for the imaging of reflection,the result in Fig. 4 is the best,not only has high resolution,but also right phase and position,the reflection imaging result of conventional migration in Fig. 2 has a 90° phase error, and the result in Fig. 3 also has a -90° phase error,both resulting in wrong positions. In order to further compare the resolution,phase and imaging position of those three methods,we plot the 232nd trace of above three migration results,which crosses the fifth scatterer and the reflector. In Fig. 5,the left picture is conventional migration result,the middle is scattered wave migration result and the right is reflected wave migration result. From the waveforms in Fig. 5,it is clear that the result from scattering migration method can get the right position and correct polarity of the crest indicating a high-speed scatterer. Meanwhile,the reflection migration method also does well in reflector imaging. It can image at the correct position with high resolution and the image has right polarity of the crest indicating a high-speed reflector. But the conventional migration image has phase error for both scatter and reflector. Its image result has 180° phase error for scatterer,which changes the high-speed scatterer into low-speed scatterer, and has 90° phase error for reflector,which results in incorrect image position and low resolution. The numerical result of this model consists with the theoretical result of this paper.

3.2 Sigbee2A Model Test

An international st and ard Sigbee2A model data is used to further verify the proposed scattering migration and reflection migration method. This model contains high-speed salt,subtle reflector and many high-speed scatterers,which is suitable for our test. Fig. 6 is the background velocity model with 3201 grid points horizontally and 1201 grid points in depth. The grid interval of both directions is 7.62 m. Marine towing observation system is used to generate synthetic seismic data. The seismic data contains 500 shots with 45.72 m shot-interval and each shot is recorded by different number of receivers with 22.86 m interval,at least 57,up to 348. The towing length is 7932.42 m. The time interval is 8 ms and time sample points is 1500.

Applying the conventional migration,the proposed scattering migration and reflection migration method to this Sigbee2A model and obtain the migration result shown in Figs. 7,8,9.

From Figs. 7,8, and 9,it is clear that for high-speed scatterers image,scattering migration gets better result than conventional migration and reflection migration. While for the imaging of reflector,reflection migration method gains the best result. Enlarging the left part of the model and the above three images to carefully compare the imaging result,as shown in Figs. 10,11,12, and 13. Fig. 10 is the magnification of the true model,Figs. 10,11,12, and 13 are respectively the magnification of conventional reverse time migration,scattered migration and reflection migration results. It is shown that the proposed scattered migration and reflection migration method separately get an ideal result for high-speed scatterers' imaging and reflector's imaging,not only have higher resolution,but also have correct imaging position and phase information. That means scattering migration benefits the imaging of scatters and reflection migration benefits the imaging of reflector.

4 CONCLUSIONS

The study of seismic migration method should be based upon the mathematical physical equation. The migration of different wave should adopt different kind of wave propagation equation. Scattering migration method established on the basis of scattering wave propagation equation should be used for scattering data generated by scatterers. As for the reflector imaging,reflection migration method established on the basis of reflection wave propagation equation,should be applied to the reflection data. Based on the scattering wave propagation equation,we use high frequency approximation to derive the reflected wave propagation equation. This is a development of the reflected wave propagation theory. There is no problem about the forward extrapolation and backward extrapolation in reverse time migration,the main problem is imaging formula,which is derived from Claerbout's time consistency imaging principle,not from scattering equation or reflection equation. The conventional migration has a 90° phase error for reflector imaging,leading to incorrect imaging position, and has an 180° phase error for scatterer imaging,resulting in reverse polarity of the image result. This paper derives different imaging formula for scatterer and reflector with wave propagation equation,. The migration result of scattering imaging is a relatively dimensionless velocity perturbation,while the migration result of reflection imaging is the reflectivity with the reciprocal of the length as its dimension. The reflection migration method proposed in this paper has better image result for reflector than conventional migration method,such as correct phase and imaging position,high resolution. For scatterer image,the scattering method proposed by this paper also has correct phase information and higher resolution than conventional migration method.

5 ACKNOWLEDGMENTS

This work was supported by the National Natural Science Foundation of China (41074133,41374001) .