- Email: [email protected]
ADCIGs extraction and reflection tomography modeling for elastic wave
Journal of Applied Geophysics 143 (2017) 203–211
Contents lists available at ScienceDirect
Journal of Applied Geophysics journal homepage: www.elsevier.com/locate/jappgeo
ADCIGs extraction and reflection tomography modeling for elastic wave Ruihe Liu a, Ning Qin b,⁎, Xingyao Yin a a b
China University of Petroleum (Huadong), Shandong, Qingdao 266580, China Shengli Geophysical Research Institute of Sinopec, Shandong, Dongying 257022, China
a r t i c l e
i n f o
Article history: Received 19 October 2016 Received in revised form 8 February 2017 Accepted 19 April 2017 Available online 23 April 2017 Keywords: Elastic wave Angle domain common imaging gathers Reflection tomography Velocity modeling Gaussian beam
a b s t r a c t Based on the theory of elastic wave, we derive the migration angle formula of P- and S-wave in the Gaussian beam migration and implement the ADCIGs (angle domain common imaging gathers) extraction for Gaussian beam of PP- and PS-wave firstly. Then, we derive the reflection tomography equations for elastic wave velocity analysis. They are used for the tomography velocity modeling base on the Gaussian beam ADCIGs. Making reflection tomography and depth migration into the same velocity modeling process, this method can not only use ADCIGs of elastic wave to update the velocities, but also can apply the migration results of tomography in each iteration as quality control, until the end of whole process in velocity modeling and finally we can get the prestack depth migration results after tomography. Synthetic and field examples validate the correctness and practicability of this method. © 2017 Elsevier B.V. All rights reserved.
1. Introduction In recent years, the exploration and exploitation of some mature oilfields become more and more difficult. Fracture and lithology reservoirs have been taken as the exploration targets. Unlike the use of the methods of P-wave data for the exploration of conventional structural reservoirs, processing and interpretation techniques of vector data are utilized for the exploration of fracture and lithology reservoirs. Usually, the vertical component is considered as P-wave and the rotated horizontal component is considered as converted wave in the process of traditional multi-wave and multi-component seismic data processing. The elastic and vector characteristics of seismic wave are not taken into account. Moreover, the precision of wave field separation influences the processing result greatly. If different waves were not separated completely, the residual energy of the other waves will damage the imaging results by causing plenty of noise. Actually, the underground media is inhomogeneous, inelastic, anisotropic and composed of multiphase components (Wu, 2006). The seismic wave propagating in the subsurface media is a kind of elastic vector wave that includes coupling P-wave and S-wave. The data processing methods based on acoustic media assumption are not available in the case of multi-wave and multi-component exploration. The related theory for elastic seismic data processing and its application has not been well developed (Pao and Varatharajulu, 1976; Mora, 1987; Alkhalifah and Tsvankin, 1996; Xue and McMechan, 2000; Duzhinin, 2003; Luth et al., 2005; Yan and Sava, 2010; Qin et al., 2012, 2013; Du and Li, 2012; Liu et al., 2013). Therefore, it is necessary to do more research on elastic media velocity modeling and imaging. ⁎ Corresponding author. E-mail address: [email protected] (N. Qin).
http://dx.doi.org/10.1016/j.jappgeo.2017.04.006 0926-9851/© 2017 Elsevier B.V. All rights reserved.
An efficient and accurate migration method that can output velocity analysis gathers is crucial for velocity modeling by reflection tomography. Currently, Kirchhoff migration algorithm has high computational efficiency, but it cannot work well in areas of complex structure due to its failure to deal with multipath problem (Kuo and Dai, 1984; Keho and Wu, 1987; Gherasim et al., 2005; Du et al., 2011). Although reverse time migration has high imaging precision, it takes large calculation costs and demands an accurate velocity model (Sun and McMechan, 1986; Li et al., 2012). Gaussian-beam is a particular description of propagating wave field based on ray theory. It is a representation of subsurface local wave field that is not only independent but also superimposed mutually. Consequently, compared to Kirchhoff migration, Gaussianbeam migration can handle the multipath problem and has higher imaging precision (Hill, 1990, 2001; Yue, 2011). Moreover, the initial wave-front of Gaussian-beam is plane. This can be used to perform plane wave decomposition by windowed local slant stack and the corresponding Gaussian-beam can be used for extrapolation and imaging. The calculation efficiency of Gaussian-beam migration is comparable with that of Kirchhoff migration and is much higher than that of reverse time migration. On the other hand, considering that the most calculation time of reflection tomography is spent on migration and angle domain gathers extraction, Gaussian-beam migration is an efficient and flexible multi-wave angle domain gathers extraction algorithm. So Gaussian-beam pre-stack depth migration is an appropriate tool for imaging and velocity modeling. In this paper, based on elastic wave theory, the P-wave and S-wave Gaussian-beam migration angle formula is deduced to extract PP-wave and PS-wave Gaussian-beam angle domain gathers. Then, a reflection tomography inversion equation is derived to perform elastic wave velocity modeling. Tomography modeling using Gaussian-beam
204
R. Liu et al. / Journal of Applied Geophysics 143 (2017) 203–211
angle domain gathers is achieved. This method combines reflection tomography and depth migration into a unified modeling procedure. It means that velocity model updating and imaging quality control can be implemented simultaneously until the entire velocity modeling is accomplished and the final pre-stack depth migration profile is obtained. Tests on synthetic and field data processing verify the validity and practicability of this method. 2. Extraction of elastic wave angle domain gathers 2.1. Migration angle of P-wave and S-wave The initial value of Gaussian-beam chosen by Hill can ensure that the travel time and amplitude vary slowly in a certain wave field scope. It means that the Gaussian-beam real value travel time information can be used directly to calculate the propagation angle. In a 2-D ray central coordinate system (Fig. 1), tA the Gaussian-beam real value travel time of point A can be described as: 1 t A ¼ t B þ n2 Re½M B 2
ð1Þ
where, MB is the second-order partial derivative of point B's travel time field with respect to ray central coordinate; n is point A's normal vector coordinate in the coordinate system taking point B as its origin; l is the ray unit tangent at point B; lx and lz are the x and z components respectively of l in an orthogonal coordinate system. Then, n ¼ ðx−xB Þlz −ðz−zB Þlx
ð2Þ
Take the derivation with respect to abscissa x on both sides of Eq. (1). ∂t A ∂t B ∂n ¼ þ n Re½M B ∂x ∂x ∂x
ð3Þ
Substitute Eq. (2) into Eq. (3) and assume pxA ¼ then
∂t A ∂t B ,p ¼ , and ∂x xB ∂x
h i 2 pxA ¼ pxB þ ðx−xB Þlz −ðz−zB Þlx lz Re½M B
ð4Þ
h i 2 pzA ¼ pzB − ðx−xB Þlx lz −ðz−zB Þlx Re½MB
ð5Þ
The grid point's propagation angle η (the angle between ray and positive direction of z axis) is: 8
px b 0; pz b 0 px N 0; pz b 0 others
domain gathers can be obtained by accumulating the imaging value in a common imaging gathers of a certain angle scope. For PS-wave imaging, the incident angle of P-wave and S-wave can be derived as follows. As shown in Fig. 2, assuming that the clockwise direction of reference axis is corresponding to positive angle, the angle of P-wave and S-wave Gaussian-beam relative to positive direction of Z axis is ϕ2 and ϕ1 respectively. n is the normal of reflection interface. And the incident angle of P-wave β and the reflection angle of S-wave θ are to be solved. Based on the geometrical relationship shown in Fig. 2, the following equation can be obtained: ϕ1 −ϕ2 ¼ θ−β According to Snell law, sinð−βÞ sin θ ¼ vp vs
ð8Þ
Substitute Eq. (7) to Eq. (8), sinð−βÞ sinðϕ1 −ϕ2 Þ cos β þ cosðϕ1 −ϕ2 Þ sin β ¼ vp vs
ð9Þ
Simplify the former equation, vp sinðϕ1 −ϕ2 Þ cos β þ vp cosðϕ1 −ϕ2 Þsin β þ vs sin β ¼ 0
ð10Þ
The incident angle of P-wave β can be obtained by solving Eq. (10), β ¼ arctan
− sinðϕ1 −ϕ2 Þ cosðϕ1 −ϕ2 Þ þ vs =vp
ð11Þ
In the same way, the reflection angle of S-wave θ is θ ¼ arctan
sinðϕ1 −ϕ2 Þ cosðϕ1 −ϕ2 Þ þ vp =vs
ð12Þ
After getting the incident angle of P-wave and the reflection angle of S-wave, the Gaussian-beam angle domain gathers extraction of elastic vector wave can be performed directly. 2.2. Test on angle domain gathers extraction Gaussian-beam angle domain gathers extraction of elastic wave is tested using a synthetic fault model shown in Fig. 3. The grid size of
ð6Þ
After getting the Gaussian-beam's propagation angle η of source and receiver, the migration angle of PP-wave can be obtained. Then, according to the migration angle, the PP-wave Gaussian-beam angle
Fig. 1. 2-D ray central coordinate system. (Blue line represents the ray and red line indicates the ray-central coordinate system).
ð7Þ
Fig. 2. The reflection and incident angles of conversion wave.
R. Liu et al. / Journal of Applied Geophysics 143 (2017) 203–211
(a) Velocity model of P-wave
205
(b) Velocity model of S-wave
Fig. 3. Velocity model of P- and S- wave.
this model is 601 ∗ 251. Both the vertical and horizontal grid intervals are 5 m. The velocity ratio of P-wave to S-wave is 1.732. The density is constant. The sample number is 3501 per trace in a shot record with 0.5 ms sample interval. The extracted angle domain gathers, angle range from −30° to 30°, at location 0.75 km and 1.50 km using 90%, 100% and 110% of the true velocity model respectively are shown in Figs. 4 and 5. It can be seen that events of both PP-wave and PSwave angle domain gathers are in a shape of a smile face while the used velocity is less than the true value. Additionally, the interfaces' imaging depth of P-wave is not consistent with that of S-wave. The events of both PP-wave and PS-wave angle domain gathers are flattened while the used velocity is equal to the true velocity and their imaging depth of the interface is consistent to each other. When the velocity is greater than the true value, events of both PP-wave and PS-wave angle domain gathers are shown in a shape of a weep face. Meanwhile, their imaging depths do not match. Some conclusions can be obtained from the test. The residual curvature of angle domain gathers is not very sensitive to velocity error. It means that the velocity error does not cause residual curvature definitely. The resolution limitation is a kind of inherent defect for the velocity analysis based on angle domain residual curvature. In addition, PP-wave angle domain gathers are only related to P-wave velocity. PS-wave angle domain gathers depend on both P-wave velocity and S-wave velocity. Therefore, PP-wave angle domain gathers are generally better than those of PS-wave.
(a) PP-wave angle domain gathers
3. Velocity model building based on elastic wave reflection tomography 3.1. Basic theory Based on the theory of acoustic wave tomography, the inversion equation of elastic wave tomography is derived. Taking sp1 and sp2 as the P-wave slowness before and after iteration respectively, Lp1 and Lp2 as the P-wave sensitivity matrix before and after iteration respectively, tpp1 and tpp2 as PP-wave ray travel-time before and after iteration respectively, the following equation can be obtained. Lp1 sp1 −Lp2 sp2 ¼ t pp1 −t pp2
ð13Þ
We assume that the variation of interface position causes weak depth iteration, which means approximately Lp1 = Lp2 = Lp, then Lp sp1 −sp2 ¼ t pp1 −t pp2
ð14Þ
The inversion equation of elastic PP-wave tomography can be written as follows: Lp Δsp ¼ Δt pp
ð15Þ
where, Lp is the sensitivity matrix of P-wave and its elements are the ray path length of P-wave velocity field within the grids. Δtpp is the residual
(b) PS-wave angle domain gathers
Fig. 4. Angle domain gathers using 90%, 100% and 110% of the true velocity at 0.75 km.
206
R. Liu et al. / Journal of Applied Geophysics 143 (2017) 203–211
(a) PP-wave angle domain gathers
(b) PS-wave angle domain gathers
Fig. 5. Angle domain gathers using 90%, 100% and 110% of the true velocity at 1.5 km.
vector of travel time. Δsp is the updated quantity of P-wave slowness which is to be inversed. Similarly, the inversion equation of elastic PS-wave tomography can be obtained. Ls Δss ¼ 2Δt ps −Δt pp
ð16Þ
where, Ls is the sensitivity matrix of S-wave and its elements are the ray path length of S-wave velocity field within the grids. Δtps is the residual vector of travel time. Δss is the updated quantity of S-wave slowness which is to be inversed. Therefore, elastic wave reflection tomography can be used to build a velocity model. The key point is the determination of sensitivity matrix and the calculation of travel time residual. The sensitivity matrix can be obtained by ray tracing in P-wave and S-wave velocity fields and its element represents the path length of ray i at grid point j. Travel time residual can be converted from angle domain residual curvature. The conversion equations of PP-wave and PS-wave are as follows: Δt pp ¼ 2sp Δzpp cos α cos β Δt ps
where, Δzpp and Δzps are the angle domain residual curvatures of PP-wave and PS-wave respectively. sp is the P-wave slowness at the imaging point. α is the dip angle of reflection interface. β is the incidence angle of P-wave or S-wave. θ is the reflection angle of S-wave. 3.2. The procedures First of all, it should be noted that converted the velocity of PS-wave depends on both P-wave velocity and S-wave velocity. The so-called compound velocity is not obtained by separation of P-wave and S-wave. The stack velocity Vps|stk of converted wave is achieved directly from velocity spectra or velocity scanning. With the assumption of small stratum dip angle, Vps|stk is the RMS (root-mean-square) velocity. The relationship between the compound velocity and the RMS velocity of P-wave and S-wave can be described as: V2ps|stk = Vp|stkVs|stk. The procedures for elastic wave reflection tomography velocity model building are as follows.
ð17Þ
tan β ¼ 1þ sp Δzps cos α cos β tan θ
(1) Input the original pre-stack multi-component seismic data. Stacking velocity analysis of P-wave and converted wave is performed to obtain their initial stacking velocity fields. The initial S-wave stacking velocity field can be calculated by the formula
ð18Þ
Final P -wave and S -wave velocity
Pre-stack mul-component seismic
Yes
No Initial interval velocity fields of P-wave and S-wave
Update interval
Meet the requirement
velocity fields of P -wave and S-wave
The principle of flattened angle domain gathers and imaging depth consistency of
Pre-stack depth migration
Ray tracing modeling
Angle domain gathers
Sensitivity matrix
multi-wave
Establish a regularization inversion equation set
Fig. 6. The workflow of elastic wave reflection tomography velocity model building.
R. Liu et al. / Journal of Applied Geophysics 143 (2017) 203–211
(a) Z component
207
(b) X component
Fig. 7. Z and X component of synthetic data.
(a) Initial velocity of P-wave
(b) Initial velocity of S-wave
Fig. 8. Initial velocity model of P- and S-wave.
V2ps | stk = Vp| stkVs| stk. And then the interval velocity field can be calculated with Dix Formula. (2) Choose the appropriate angle scope in the initial interval velocity fields of P-wave and S-wave according to the requirement for the prestack seismic processing. And then, ray tracing is implemented with a proper fixed step length which is determined by the precision requirement of velocity analysis in order to get the sensitivity matrix of P-wave and S-wave. (3) Perform Gaussian beam pre-stack depth migration to the multicomponent seismic data using the initial interval velocity fields of P-wave and S-wave to get the imaging results of PP-wave and PS-wave and their corresponding angle domain gathers.
(a) PP-wave
Afterwards, the travel time residual vector is calculated according to the conversion relationship between the travel time residual and the depth residual (residual curvature of angle domain gathers). (4) Taking the prior and constraint information into account, establish the corresponding regularization matrix and the joint inversion equation set of travel time tomography in the elastic vector wave imaging domain. Solve the equation set with linear optimization algorithm to obtain the slowness update quantity of P-wave and S-wave, and to update the velocity models. (5) The criteria of flattening the events in angle domain gathers, the depth consistency of multi-wave image, and the precision
(b) PS-wave
Fig. 9. Angle domain gathers of PP- and PS-wave for initial velocity model (left) and tomography velocity model (right).
208
R. Liu et al. / Journal of Applied Geophysics 143 (2017) 203–211
(a) Tomography velocity of P-wave
(b) Tomography velocity of S-wave
Fig. 10. Tomography velocity model of P- and S-wave.
(a) PP-wave migration
(b) PS-wave migration
Fig. 11. PP- and PS-wave migration profile with tomography velocity model.
requirement of tomography are used to determine whether more iterations are needed. If more iterations are needed, go back to step (2) to update tomography under the condition that the interval velocity and depth of overlying layers are kept unchanged. This loop will not be terminated until the results meet the criteria. Finally, output the final joint imaging results of elastic vector waves and flattened angle domain gathers.
4. Synthetic and field examples 4.1. Synthetic example The synthetic fault model shown in Fig. 3 is used to test the elastic wave joint tomography method. The grid size of this model is 601 ∗ 251. Both the vertical and horizontal grid intervals are 5 m. There are 126 shots of multi-wave seismic data totally, shown in Fig. 7.
(a) Comparison of P-wave velocity (b) Comparison of S-wave velocity Fig. 12. Comparison of true, initial and tomography P- and S-wave velocity.
R. Liu et al. / Journal of Applied Geophysics 143 (2017) 203–211
(a) Z component
209
(b) X component
Fig. 13. Z and X component of the seismic data.
Each shot is composed of 301 traces. The sample number per trace is 3500 in a shot record with 0.5 ms sample interval. The angles of angle domain gathers range from −30° to 30° with the interval of 1°. The velocity ratio of P-wave to S-wave is 1.732. The density is constant. Following the workflow described in Fig. 6, the initial velocity fields of P-wave and S-wave are obtained, shown in Fig. 8. The initial Gaussianbeam migration angle domain gathers of PP- and PS-wave are shown in the left hand side of Fig. 9. Because of the incorrect velocity, the events of angle domain gathers are in the shape of a smile face. The flattened angle domain gathers of reflection tomography is shown in the right hand side of Fig. 9. The imaging depth of PP-wave is consistent with
(a) Initial velocity of P-wave
that of and PS-wave. The built velocity model is shown in Fig. 10, the corresponding migration profiles are shown in Fig. 11. We compare the initial velocity fields of P-wave and S-wave, tomography velocity model, and true velocity model, and find that the tomography velocity model matches the synthetic velocity model very well except for the lowest layer caused by limited angle domain gathers. Moreover, the inverted P-wave velocity model is more reliable than that of S-wave as shown in Fig. 12. We can see the tomography model does a poor job of reconstructing the lowest velocity layer because the data do not contain information about the velocity in this layer which lies beneath the deepest reflector.
(b) Initial velocity of S-wave
Fig. 14. Initial velocity model of P- and S- wave.
(a) PP-wave
(b) PS-wave
Fig. 15. Tomography angle domain gathers of PP- and PS- wave.
210
R. Liu et al. / Journal of Applied Geophysics 143 (2017) 203–211
(a) Tomography velocity of P-wave
(b) Tomography velocity of S-wave
Fig. 16. Tomography velocity model of P- and S- wave.
(a) PP-wave migration
(b) PS-wave migration
Fig. 17. PP- and PS- wave migration profile with tomography velocity model.
4.2. Field example The proposed method is also applied to the multi-component seismic data in the oilfield. The recording time lengths of Z and X component for each trace are 6 s and 7 s, respectively. The sample interval is 4 ms. The trace interval is 12 m. The max offset is 2700 m. Fig. 13 shows that x component data have poor quality than z component data due to more missing traces and lower signal to noise ratio. The angles of gathers range from 0° to 34° with 1° interval. The density is constant. The initial velocity models of P-wave and S-wave are shown in Fig. 14. The initial Gaussian-beam migration angle domain gathers of PP-wave and PS-wave are shown in the left hand side of Fig. 15. Because of the incorrect velocity, the events of angle domain gathers are upward bending and contaminated seriously by shallow noise. The flattened angle domain gathers of reflection tomography are shown in the right hand side of Fig. 16. The imaging depth of PP-wave is almost consistent with that of PS-wave. The updated velocity models are shown in Fig. 16. And the corresponding migration profiles are shown in Fig. 17. It is obvious that the imaging quality of PP-wave is better than that of PS-wave. The first reason is that the quality of z component raw data is better than that of x component data due to higher signal to noise ratio and less missing traces. The second one is that PP-wave imaging only relies on P-wave velocity, and PS-wave imaging depends on both P-wave velocity and S-wave velocity. For field data processing, pre-stack pretreatment has not been done. Due to the resolution of angle domain gathers, low signal to noise ratio of multi-component data can have a great impact on the result. 5. Conclusions In the paper, the migration angle formulas of P- and S-wave in the Gaussian beam migration are derived based on the theory of elastic
wave, and the ADCIGs (angle domain common imaging gathers) extractions for Gaussian beam of PP- and PS-wave are provided. The reflection tomography equations for elastic wave velocity analysis and the conversion equations between travel time residual and ADCISs residual curvature are established. These equations are used to implement the velocity tomography base on the Gaussian beam ADCIGs. The method of velocity model building takes advantage both of the rich information of elastic wave ADCIGs and of the calculating stability and efficiency of Gaussian beam migration. It provides an accurate velocity model and high quality imaging result. Synthetic and field examples validate that the proposed method has higher accuracy for multi-wave joint inversion and better elastic wave pre-stack migration quality. Additionally, the inversion result of P-wave velocity has higher reliability than that of S-wave. Since the pre-stack data with low signal to noise ratio can influence the tomography results critically, preprocessing of pre-stack multi-wave data is necessary to guarantee the algorithm precision. References Alkhalifah, T., Tsvankin, 1996. Velocity analysis and imaging in transversely isotropic media: methodology and a case study (J). Lead. Edge 15 (2), 371–378. Du, Q.Z., Li, F., 2012. Multicomponent joint migration velocity analysis in the angle domain for PP-waves and PS-waves (J). Geophysics 77 (1), U1–U13. Du, Q.Z., Li, F., Hou, B., et al., 2011. Angle-domain migration velocity analysis based on elastic-wave Kirchhoff prestack depth migration. Chin. J. Geophys. 54 (5), 1327–1339 (in Chinese). Duzhinin, A., 2003. Decoupled elastic prestack depth migration (J). J. Appl. Geophys. 54 (1), 369–389. Gherasim, M., Hoelting, C., Marfurt, K., 2005. 3-D VSP elastic Kirchhoff pre-stack depth migration-Vinton Dome, Louisiana (C). SEG Expanded Abstracts, pp. 2649–2652. Hill, N.R., 1990. Gaussian beam migration. Geophysics 55 (11), 1416–1428. Hill, N.R., 2001. Prestack Gaussian-beam depth migration. Geophysics 66 (4), 1240–1250. Keho, K.H., Wu, R.S., 1987. Elastic Kirchhoff migration for vertical seismic profiles (C). SEG Expanded Abstracts, pp. 774–776. Kuo, J.T., Dai, T., 1984. Kirchhoff elastic wave migration for the case of noncoincident source and receiver (J). Geophysics 49 (5), 1223–1238.
R. Liu et al. / Journal of Applied Geophysics 143 (2017) 203–211 Li, W.J., Qu, S.L., Wei, X.C., et al., 2012. Discussion about reverse-time depth migration technology of prestack elastic wavefield for offset VSP data. Chin. J. Geophys. 55 (1), 238–251 (in Chinese). Liu, X.X., Yin, X.Y., Liu, B.H., 2013. Prestack reverse-time migration of elastic waves in porous media. Prog. Geophys. 28 (6), 3165–3173 (in Chinese). Luth, S., Buske, S., Giese, R., 2005. Fresnel volume migration of multi-component data (J). Geophysics 70 (6), S121–S129. Mora, P., 1987. Nonlinear 2-D elastic inversion of multi-offset seismic data (J). Geophysics 52 (7), 1200–1228. Pao, Y.H., Varatharajulu, V., 1976. Huygens' principle, radiation conditions and integral formulas for the scattering of elastic waves (J). J. Acoust. Soc. Am. 59 (4), 1361–1371. Qin, N., Li, Z.C., Zhang, K., 2012. A velocity tomography method in the angle domain for elastic vector wavefield (C). SEG Expanded Abstracts, pp. 1941–1945.
211
Qin, N., Li, Z.C., Zhang, K., et al., 2013. Joint tomography inversion of image domain for elastic vector wavefield. Chin. J. Geophys. 56 (9), 3109–3117 (in Chinese). Sun, R., McMechan, G.A., 1986. Pre-stack reverse-time migration for elastic waves with application to synthetic offset vertical seismic profiles (J). Proc. Inst. Electr. Electron. Eng. 74 (3), 457–465. Wu, G.C., 2006. Seismic Wave in Anisotropic Media Propagation and Imaging (M). China University of Petroleum Press, Dongying. Xue, A., McMechan, G.A., 2000. Prestack elastic Kirchhoff migration for multicomponent seismic data in variable velocity media (C). SEG Expanded Abstracts, pp. 449–452. Yan, J., Sava, P., 2010. Analysis of converted-wave extended images for migration velocity analysis (C). SEG Expanded Abstracts, pp. 1666–1669. Yue, Y.B., 2011. Study on Gaussian Beam Migration Methods in Complex Medium (D). China University of Petroleum (Huadong), Qingdao.
Contents lists available at ScienceDirect
Journal of Applied Geophysics journal homepage: www.elsevier.com/locate/jappgeo
ADCIGs extraction and reflection tomography modeling for elastic wave Ruihe Liu a, Ning Qin b,⁎, Xingyao Yin a a b
China University of Petroleum (Huadong), Shandong, Qingdao 266580, China Shengli Geophysical Research Institute of Sinopec, Shandong, Dongying 257022, China
a r t i c l e
i n f o
Article history: Received 19 October 2016 Received in revised form 8 February 2017 Accepted 19 April 2017 Available online 23 April 2017 Keywords: Elastic wave Angle domain common imaging gathers Reflection tomography Velocity modeling Gaussian beam
a b s t r a c t Based on the theory of elastic wave, we derive the migration angle formula of P- and S-wave in the Gaussian beam migration and implement the ADCIGs (angle domain common imaging gathers) extraction for Gaussian beam of PP- and PS-wave firstly. Then, we derive the reflection tomography equations for elastic wave velocity analysis. They are used for the tomography velocity modeling base on the Gaussian beam ADCIGs. Making reflection tomography and depth migration into the same velocity modeling process, this method can not only use ADCIGs of elastic wave to update the velocities, but also can apply the migration results of tomography in each iteration as quality control, until the end of whole process in velocity modeling and finally we can get the prestack depth migration results after tomography. Synthetic and field examples validate the correctness and practicability of this method. © 2017 Elsevier B.V. All rights reserved.
1. Introduction In recent years, the exploration and exploitation of some mature oilfields become more and more difficult. Fracture and lithology reservoirs have been taken as the exploration targets. Unlike the use of the methods of P-wave data for the exploration of conventional structural reservoirs, processing and interpretation techniques of vector data are utilized for the exploration of fracture and lithology reservoirs. Usually, the vertical component is considered as P-wave and the rotated horizontal component is considered as converted wave in the process of traditional multi-wave and multi-component seismic data processing. The elastic and vector characteristics of seismic wave are not taken into account. Moreover, the precision of wave field separation influences the processing result greatly. If different waves were not separated completely, the residual energy of the other waves will damage the imaging results by causing plenty of noise. Actually, the underground media is inhomogeneous, inelastic, anisotropic and composed of multiphase components (Wu, 2006). The seismic wave propagating in the subsurface media is a kind of elastic vector wave that includes coupling P-wave and S-wave. The data processing methods based on acoustic media assumption are not available in the case of multi-wave and multi-component exploration. The related theory for elastic seismic data processing and its application has not been well developed (Pao and Varatharajulu, 1976; Mora, 1987; Alkhalifah and Tsvankin, 1996; Xue and McMechan, 2000; Duzhinin, 2003; Luth et al., 2005; Yan and Sava, 2010; Qin et al., 2012, 2013; Du and Li, 2012; Liu et al., 2013). Therefore, it is necessary to do more research on elastic media velocity modeling and imaging. ⁎ Corresponding author. E-mail address: [email protected] (N. Qin).
http://dx.doi.org/10.1016/j.jappgeo.2017.04.006 0926-9851/© 2017 Elsevier B.V. All rights reserved.
An efficient and accurate migration method that can output velocity analysis gathers is crucial for velocity modeling by reflection tomography. Currently, Kirchhoff migration algorithm has high computational efficiency, but it cannot work well in areas of complex structure due to its failure to deal with multipath problem (Kuo and Dai, 1984; Keho and Wu, 1987; Gherasim et al., 2005; Du et al., 2011). Although reverse time migration has high imaging precision, it takes large calculation costs and demands an accurate velocity model (Sun and McMechan, 1986; Li et al., 2012). Gaussian-beam is a particular description of propagating wave field based on ray theory. It is a representation of subsurface local wave field that is not only independent but also superimposed mutually. Consequently, compared to Kirchhoff migration, Gaussianbeam migration can handle the multipath problem and has higher imaging precision (Hill, 1990, 2001; Yue, 2011). Moreover, the initial wave-front of Gaussian-beam is plane. This can be used to perform plane wave decomposition by windowed local slant stack and the corresponding Gaussian-beam can be used for extrapolation and imaging. The calculation efficiency of Gaussian-beam migration is comparable with that of Kirchhoff migration and is much higher than that of reverse time migration. On the other hand, considering that the most calculation time of reflection tomography is spent on migration and angle domain gathers extraction, Gaussian-beam migration is an efficient and flexible multi-wave angle domain gathers extraction algorithm. So Gaussian-beam pre-stack depth migration is an appropriate tool for imaging and velocity modeling. In this paper, based on elastic wave theory, the P-wave and S-wave Gaussian-beam migration angle formula is deduced to extract PP-wave and PS-wave Gaussian-beam angle domain gathers. Then, a reflection tomography inversion equation is derived to perform elastic wave velocity modeling. Tomography modeling using Gaussian-beam
204
R. Liu et al. / Journal of Applied Geophysics 143 (2017) 203–211
angle domain gathers is achieved. This method combines reflection tomography and depth migration into a unified modeling procedure. It means that velocity model updating and imaging quality control can be implemented simultaneously until the entire velocity modeling is accomplished and the final pre-stack depth migration profile is obtained. Tests on synthetic and field data processing verify the validity and practicability of this method. 2. Extraction of elastic wave angle domain gathers 2.1. Migration angle of P-wave and S-wave The initial value of Gaussian-beam chosen by Hill can ensure that the travel time and amplitude vary slowly in a certain wave field scope. It means that the Gaussian-beam real value travel time information can be used directly to calculate the propagation angle. In a 2-D ray central coordinate system (Fig. 1), tA the Gaussian-beam real value travel time of point A can be described as: 1 t A ¼ t B þ n2 Re½M B 2
ð1Þ
where, MB is the second-order partial derivative of point B's travel time field with respect to ray central coordinate; n is point A's normal vector coordinate in the coordinate system taking point B as its origin; l is the ray unit tangent at point B; lx and lz are the x and z components respectively of l in an orthogonal coordinate system. Then, n ¼ ðx−xB Þlz −ðz−zB Þlx
ð2Þ
Take the derivation with respect to abscissa x on both sides of Eq. (1). ∂t A ∂t B ∂n ¼ þ n Re½M B ∂x ∂x ∂x
ð3Þ
Substitute Eq. (2) into Eq. (3) and assume pxA ¼ then
∂t A ∂t B ,p ¼ , and ∂x xB ∂x
h i 2 pxA ¼ pxB þ ðx−xB Þlz −ðz−zB Þlx lz Re½M B
ð4Þ
h i 2 pzA ¼ pzB − ðx−xB Þlx lz −ðz−zB Þlx Re½MB
ð5Þ
The grid point's propagation angle η (the angle between ray and positive direction of z axis) is: 8
px b 0; pz b 0 px N 0; pz b 0 others
domain gathers can be obtained by accumulating the imaging value in a common imaging gathers of a certain angle scope. For PS-wave imaging, the incident angle of P-wave and S-wave can be derived as follows. As shown in Fig. 2, assuming that the clockwise direction of reference axis is corresponding to positive angle, the angle of P-wave and S-wave Gaussian-beam relative to positive direction of Z axis is ϕ2 and ϕ1 respectively. n is the normal of reflection interface. And the incident angle of P-wave β and the reflection angle of S-wave θ are to be solved. Based on the geometrical relationship shown in Fig. 2, the following equation can be obtained: ϕ1 −ϕ2 ¼ θ−β According to Snell law, sinð−βÞ sin θ ¼ vp vs
ð8Þ
Substitute Eq. (7) to Eq. (8), sinð−βÞ sinðϕ1 −ϕ2 Þ cos β þ cosðϕ1 −ϕ2 Þ sin β ¼ vp vs
ð9Þ
Simplify the former equation, vp sinðϕ1 −ϕ2 Þ cos β þ vp cosðϕ1 −ϕ2 Þsin β þ vs sin β ¼ 0
ð10Þ
The incident angle of P-wave β can be obtained by solving Eq. (10), β ¼ arctan
− sinðϕ1 −ϕ2 Þ cosðϕ1 −ϕ2 Þ þ vs =vp
ð11Þ
In the same way, the reflection angle of S-wave θ is θ ¼ arctan
sinðϕ1 −ϕ2 Þ cosðϕ1 −ϕ2 Þ þ vp =vs
ð12Þ
After getting the incident angle of P-wave and the reflection angle of S-wave, the Gaussian-beam angle domain gathers extraction of elastic vector wave can be performed directly. 2.2. Test on angle domain gathers extraction Gaussian-beam angle domain gathers extraction of elastic wave is tested using a synthetic fault model shown in Fig. 3. The grid size of
ð6Þ
After getting the Gaussian-beam's propagation angle η of source and receiver, the migration angle of PP-wave can be obtained. Then, according to the migration angle, the PP-wave Gaussian-beam angle
Fig. 1. 2-D ray central coordinate system. (Blue line represents the ray and red line indicates the ray-central coordinate system).
ð7Þ
Fig. 2. The reflection and incident angles of conversion wave.
R. Liu et al. / Journal of Applied Geophysics 143 (2017) 203–211
(a) Velocity model of P-wave
205
(b) Velocity model of S-wave
Fig. 3. Velocity model of P- and S- wave.
this model is 601 ∗ 251. Both the vertical and horizontal grid intervals are 5 m. The velocity ratio of P-wave to S-wave is 1.732. The density is constant. The sample number is 3501 per trace in a shot record with 0.5 ms sample interval. The extracted angle domain gathers, angle range from −30° to 30°, at location 0.75 km and 1.50 km using 90%, 100% and 110% of the true velocity model respectively are shown in Figs. 4 and 5. It can be seen that events of both PP-wave and PSwave angle domain gathers are in a shape of a smile face while the used velocity is less than the true value. Additionally, the interfaces' imaging depth of P-wave is not consistent with that of S-wave. The events of both PP-wave and PS-wave angle domain gathers are flattened while the used velocity is equal to the true velocity and their imaging depth of the interface is consistent to each other. When the velocity is greater than the true value, events of both PP-wave and PS-wave angle domain gathers are shown in a shape of a weep face. Meanwhile, their imaging depths do not match. Some conclusions can be obtained from the test. The residual curvature of angle domain gathers is not very sensitive to velocity error. It means that the velocity error does not cause residual curvature definitely. The resolution limitation is a kind of inherent defect for the velocity analysis based on angle domain residual curvature. In addition, PP-wave angle domain gathers are only related to P-wave velocity. PS-wave angle domain gathers depend on both P-wave velocity and S-wave velocity. Therefore, PP-wave angle domain gathers are generally better than those of PS-wave.
(a) PP-wave angle domain gathers
3. Velocity model building based on elastic wave reflection tomography 3.1. Basic theory Based on the theory of acoustic wave tomography, the inversion equation of elastic wave tomography is derived. Taking sp1 and sp2 as the P-wave slowness before and after iteration respectively, Lp1 and Lp2 as the P-wave sensitivity matrix before and after iteration respectively, tpp1 and tpp2 as PP-wave ray travel-time before and after iteration respectively, the following equation can be obtained. Lp1 sp1 −Lp2 sp2 ¼ t pp1 −t pp2
ð13Þ
We assume that the variation of interface position causes weak depth iteration, which means approximately Lp1 = Lp2 = Lp, then Lp sp1 −sp2 ¼ t pp1 −t pp2
ð14Þ
The inversion equation of elastic PP-wave tomography can be written as follows: Lp Δsp ¼ Δt pp
ð15Þ
where, Lp is the sensitivity matrix of P-wave and its elements are the ray path length of P-wave velocity field within the grids. Δtpp is the residual
(b) PS-wave angle domain gathers
Fig. 4. Angle domain gathers using 90%, 100% and 110% of the true velocity at 0.75 km.
206
R. Liu et al. / Journal of Applied Geophysics 143 (2017) 203–211
(a) PP-wave angle domain gathers
(b) PS-wave angle domain gathers
Fig. 5. Angle domain gathers using 90%, 100% and 110% of the true velocity at 1.5 km.
vector of travel time. Δsp is the updated quantity of P-wave slowness which is to be inversed. Similarly, the inversion equation of elastic PS-wave tomography can be obtained. Ls Δss ¼ 2Δt ps −Δt pp
ð16Þ
where, Ls is the sensitivity matrix of S-wave and its elements are the ray path length of S-wave velocity field within the grids. Δtps is the residual vector of travel time. Δss is the updated quantity of S-wave slowness which is to be inversed. Therefore, elastic wave reflection tomography can be used to build a velocity model. The key point is the determination of sensitivity matrix and the calculation of travel time residual. The sensitivity matrix can be obtained by ray tracing in P-wave and S-wave velocity fields and its element represents the path length of ray i at grid point j. Travel time residual can be converted from angle domain residual curvature. The conversion equations of PP-wave and PS-wave are as follows: Δt pp ¼ 2sp Δzpp cos α cos β Δt ps
where, Δzpp and Δzps are the angle domain residual curvatures of PP-wave and PS-wave respectively. sp is the P-wave slowness at the imaging point. α is the dip angle of reflection interface. β is the incidence angle of P-wave or S-wave. θ is the reflection angle of S-wave. 3.2. The procedures First of all, it should be noted that converted the velocity of PS-wave depends on both P-wave velocity and S-wave velocity. The so-called compound velocity is not obtained by separation of P-wave and S-wave. The stack velocity Vps|stk of converted wave is achieved directly from velocity spectra or velocity scanning. With the assumption of small stratum dip angle, Vps|stk is the RMS (root-mean-square) velocity. The relationship between the compound velocity and the RMS velocity of P-wave and S-wave can be described as: V2ps|stk = Vp|stkVs|stk. The procedures for elastic wave reflection tomography velocity model building are as follows.
ð17Þ
tan β ¼ 1þ sp Δzps cos α cos β tan θ
(1) Input the original pre-stack multi-component seismic data. Stacking velocity analysis of P-wave and converted wave is performed to obtain their initial stacking velocity fields. The initial S-wave stacking velocity field can be calculated by the formula
ð18Þ
Final P -wave and S -wave velocity
Pre-stack mul-component seismic
Yes
No Initial interval velocity fields of P-wave and S-wave
Update interval
Meet the requirement
velocity fields of P -wave and S-wave
The principle of flattened angle domain gathers and imaging depth consistency of
Pre-stack depth migration
Ray tracing modeling
Angle domain gathers
Sensitivity matrix
multi-wave
Establish a regularization inversion equation set
Fig. 6. The workflow of elastic wave reflection tomography velocity model building.
R. Liu et al. / Journal of Applied Geophysics 143 (2017) 203–211
(a) Z component
207
(b) X component
Fig. 7. Z and X component of synthetic data.
(a) Initial velocity of P-wave
(b) Initial velocity of S-wave
Fig. 8. Initial velocity model of P- and S-wave.
V2ps | stk = Vp| stkVs| stk. And then the interval velocity field can be calculated with Dix Formula. (2) Choose the appropriate angle scope in the initial interval velocity fields of P-wave and S-wave according to the requirement for the prestack seismic processing. And then, ray tracing is implemented with a proper fixed step length which is determined by the precision requirement of velocity analysis in order to get the sensitivity matrix of P-wave and S-wave. (3) Perform Gaussian beam pre-stack depth migration to the multicomponent seismic data using the initial interval velocity fields of P-wave and S-wave to get the imaging results of PP-wave and PS-wave and their corresponding angle domain gathers.
(a) PP-wave
Afterwards, the travel time residual vector is calculated according to the conversion relationship between the travel time residual and the depth residual (residual curvature of angle domain gathers). (4) Taking the prior and constraint information into account, establish the corresponding regularization matrix and the joint inversion equation set of travel time tomography in the elastic vector wave imaging domain. Solve the equation set with linear optimization algorithm to obtain the slowness update quantity of P-wave and S-wave, and to update the velocity models. (5) The criteria of flattening the events in angle domain gathers, the depth consistency of multi-wave image, and the precision
(b) PS-wave
Fig. 9. Angle domain gathers of PP- and PS-wave for initial velocity model (left) and tomography velocity model (right).
208
R. Liu et al. / Journal of Applied Geophysics 143 (2017) 203–211
(a) Tomography velocity of P-wave
(b) Tomography velocity of S-wave
Fig. 10. Tomography velocity model of P- and S-wave.
(a) PP-wave migration
(b) PS-wave migration
Fig. 11. PP- and PS-wave migration profile with tomography velocity model.
requirement of tomography are used to determine whether more iterations are needed. If more iterations are needed, go back to step (2) to update tomography under the condition that the interval velocity and depth of overlying layers are kept unchanged. This loop will not be terminated until the results meet the criteria. Finally, output the final joint imaging results of elastic vector waves and flattened angle domain gathers.
4. Synthetic and field examples 4.1. Synthetic example The synthetic fault model shown in Fig. 3 is used to test the elastic wave joint tomography method. The grid size of this model is 601 ∗ 251. Both the vertical and horizontal grid intervals are 5 m. There are 126 shots of multi-wave seismic data totally, shown in Fig. 7.
(a) Comparison of P-wave velocity (b) Comparison of S-wave velocity Fig. 12. Comparison of true, initial and tomography P- and S-wave velocity.
R. Liu et al. / Journal of Applied Geophysics 143 (2017) 203–211
(a) Z component
209
(b) X component
Fig. 13. Z and X component of the seismic data.
Each shot is composed of 301 traces. The sample number per trace is 3500 in a shot record with 0.5 ms sample interval. The angles of angle domain gathers range from −30° to 30° with the interval of 1°. The velocity ratio of P-wave to S-wave is 1.732. The density is constant. Following the workflow described in Fig. 6, the initial velocity fields of P-wave and S-wave are obtained, shown in Fig. 8. The initial Gaussianbeam migration angle domain gathers of PP- and PS-wave are shown in the left hand side of Fig. 9. Because of the incorrect velocity, the events of angle domain gathers are in the shape of a smile face. The flattened angle domain gathers of reflection tomography is shown in the right hand side of Fig. 9. The imaging depth of PP-wave is consistent with
(a) Initial velocity of P-wave
that of and PS-wave. The built velocity model is shown in Fig. 10, the corresponding migration profiles are shown in Fig. 11. We compare the initial velocity fields of P-wave and S-wave, tomography velocity model, and true velocity model, and find that the tomography velocity model matches the synthetic velocity model very well except for the lowest layer caused by limited angle domain gathers. Moreover, the inverted P-wave velocity model is more reliable than that of S-wave as shown in Fig. 12. We can see the tomography model does a poor job of reconstructing the lowest velocity layer because the data do not contain information about the velocity in this layer which lies beneath the deepest reflector.
(b) Initial velocity of S-wave
Fig. 14. Initial velocity model of P- and S- wave.
(a) PP-wave
(b) PS-wave
Fig. 15. Tomography angle domain gathers of PP- and PS- wave.
210
R. Liu et al. / Journal of Applied Geophysics 143 (2017) 203–211
(a) Tomography velocity of P-wave
(b) Tomography velocity of S-wave
Fig. 16. Tomography velocity model of P- and S- wave.
(a) PP-wave migration
(b) PS-wave migration
Fig. 17. PP- and PS- wave migration profile with tomography velocity model.
4.2. Field example The proposed method is also applied to the multi-component seismic data in the oilfield. The recording time lengths of Z and X component for each trace are 6 s and 7 s, respectively. The sample interval is 4 ms. The trace interval is 12 m. The max offset is 2700 m. Fig. 13 shows that x component data have poor quality than z component data due to more missing traces and lower signal to noise ratio. The angles of gathers range from 0° to 34° with 1° interval. The density is constant. The initial velocity models of P-wave and S-wave are shown in Fig. 14. The initial Gaussian-beam migration angle domain gathers of PP-wave and PS-wave are shown in the left hand side of Fig. 15. Because of the incorrect velocity, the events of angle domain gathers are upward bending and contaminated seriously by shallow noise. The flattened angle domain gathers of reflection tomography are shown in the right hand side of Fig. 16. The imaging depth of PP-wave is almost consistent with that of PS-wave. The updated velocity models are shown in Fig. 16. And the corresponding migration profiles are shown in Fig. 17. It is obvious that the imaging quality of PP-wave is better than that of PS-wave. The first reason is that the quality of z component raw data is better than that of x component data due to higher signal to noise ratio and less missing traces. The second one is that PP-wave imaging only relies on P-wave velocity, and PS-wave imaging depends on both P-wave velocity and S-wave velocity. For field data processing, pre-stack pretreatment has not been done. Due to the resolution of angle domain gathers, low signal to noise ratio of multi-component data can have a great impact on the result. 5. Conclusions In the paper, the migration angle formulas of P- and S-wave in the Gaussian beam migration are derived based on the theory of elastic
wave, and the ADCIGs (angle domain common imaging gathers) extractions for Gaussian beam of PP- and PS-wave are provided. The reflection tomography equations for elastic wave velocity analysis and the conversion equations between travel time residual and ADCISs residual curvature are established. These equations are used to implement the velocity tomography base on the Gaussian beam ADCIGs. The method of velocity model building takes advantage both of the rich information of elastic wave ADCIGs and of the calculating stability and efficiency of Gaussian beam migration. It provides an accurate velocity model and high quality imaging result. Synthetic and field examples validate that the proposed method has higher accuracy for multi-wave joint inversion and better elastic wave pre-stack migration quality. Additionally, the inversion result of P-wave velocity has higher reliability than that of S-wave. Since the pre-stack data with low signal to noise ratio can influence the tomography results critically, preprocessing of pre-stack multi-wave data is necessary to guarantee the algorithm precision. References Alkhalifah, T., Tsvankin, 1996. Velocity analysis and imaging in transversely isotropic media: methodology and a case study (J). Lead. Edge 15 (2), 371–378. Du, Q.Z., Li, F., 2012. Multicomponent joint migration velocity analysis in the angle domain for PP-waves and PS-waves (J). Geophysics 77 (1), U1–U13. Du, Q.Z., Li, F., Hou, B., et al., 2011. Angle-domain migration velocity analysis based on elastic-wave Kirchhoff prestack depth migration. Chin. J. Geophys. 54 (5), 1327–1339 (in Chinese). Duzhinin, A., 2003. Decoupled elastic prestack depth migration (J). J. Appl. Geophys. 54 (1), 369–389. Gherasim, M., Hoelting, C., Marfurt, K., 2005. 3-D VSP elastic Kirchhoff pre-stack depth migration-Vinton Dome, Louisiana (C). SEG Expanded Abstracts, pp. 2649–2652. Hill, N.R., 1990. Gaussian beam migration. Geophysics 55 (11), 1416–1428. Hill, N.R., 2001. Prestack Gaussian-beam depth migration. Geophysics 66 (4), 1240–1250. Keho, K.H., Wu, R.S., 1987. Elastic Kirchhoff migration for vertical seismic profiles (C). SEG Expanded Abstracts, pp. 774–776. Kuo, J.T., Dai, T., 1984. Kirchhoff elastic wave migration for the case of noncoincident source and receiver (J). Geophysics 49 (5), 1223–1238.
R. Liu et al. / Journal of Applied Geophysics 143 (2017) 203–211 Li, W.J., Qu, S.L., Wei, X.C., et al., 2012. Discussion about reverse-time depth migration technology of prestack elastic wavefield for offset VSP data. Chin. J. Geophys. 55 (1), 238–251 (in Chinese). Liu, X.X., Yin, X.Y., Liu, B.H., 2013. Prestack reverse-time migration of elastic waves in porous media. Prog. Geophys. 28 (6), 3165–3173 (in Chinese). Luth, S., Buske, S., Giese, R., 2005. Fresnel volume migration of multi-component data (J). Geophysics 70 (6), S121–S129. Mora, P., 1987. Nonlinear 2-D elastic inversion of multi-offset seismic data (J). Geophysics 52 (7), 1200–1228. Pao, Y.H., Varatharajulu, V., 1976. Huygens' principle, radiation conditions and integral formulas for the scattering of elastic waves (J). J. Acoust. Soc. Am. 59 (4), 1361–1371. Qin, N., Li, Z.C., Zhang, K., 2012. A velocity tomography method in the angle domain for elastic vector wavefield (C). SEG Expanded Abstracts, pp. 1941–1945.
211
Qin, N., Li, Z.C., Zhang, K., et al., 2013. Joint tomography inversion of image domain for elastic vector wavefield. Chin. J. Geophys. 56 (9), 3109–3117 (in Chinese). Sun, R., McMechan, G.A., 1986. Pre-stack reverse-time migration for elastic waves with application to synthetic offset vertical seismic profiles (J). Proc. Inst. Electr. Electron. Eng. 74 (3), 457–465. Wu, G.C., 2006. Seismic Wave in Anisotropic Media Propagation and Imaging (M). China University of Petroleum Press, Dongying. Xue, A., McMechan, G.A., 2000. Prestack elastic Kirchhoff migration for multicomponent seismic data in variable velocity media (C). SEG Expanded Abstracts, pp. 449–452. Yan, J., Sava, P., 2010. Analysis of converted-wave extended images for migration velocity analysis (C). SEG Expanded Abstracts, pp. 1666–1669. Yue, Y.B., 2011. Study on Gaussian Beam Migration Methods in Complex Medium (D). China University of Petroleum (Huadong), Qingdao.