- Email: [email protected]
Vector-based excitation amplitude imaging condition for elastic RTM
Vector-based excitation amplitude imaging condition for elastic RTM Jinju Zhou, Deli Wang PII: DOI: Reference:
S0926-9851(17)30377-4 doi:10.1016/j.jappgeo.2017.10.003 APPGEO 3349
To appear in:
Journal of Applied Geophysics
Received date: Revised date: Accepted date:
16 April 2017 10 July 2017 10 October 2017
Please cite this article as: Zhou, Jinju, Wang, Deli, Vector-based excitation amplitude imaging condition for elastic RTM, Journal of Applied Geophysics (2017), doi:10.1016/j.jappgeo.2017.10.003
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT
Vector-based excitation amplitude imaging condition for elastic RTM
SC
130026, China
RI
College of Geoexploration Science and Technology, Jilin University, Changchun
TE
D
MA
NU
E-mail: [email protected]; [email protected]
AC CE P
1
PT
Jinju Zhou1, Deli Wang1
ACCEPTED MANUSCRIPT
ABSTRACT In recent years, many studies have focused on elastic reverse time migration (RTM). In
PT
response to the problems associated with elastic RTM, we propose a new procedure for 2D
RI
elastic multicomponent RTM. In this new method, decomposed P- and S-wave components are
SC
obtained from the decoupled propagation of the source and receiver wavefields, which allows the expedient calculation of the Poynting vectors and the incident and reflection angles of the P- and
NU
S-waves. In addition, we deduce the vector-based excitation amplitude imaging condition. This
MA
process automatically accounts for the particle vibration directions when determining the angledependent signed reflection coefficients, and does not require the sign to be determined apart
D
from the value of the reflection coefficients. This concept was further extended to the source-
TE
normalized crosscorrelation imaging condition. The reflection coefficient of the layered model
AC CE P
test was in agreement with the Zoeppritz theory, the PP and PS wave images of the Marmousi II model were clear, and the PS wave images had higher resolution and richer details. In addition, since the calculated reflection coefficients are angle-dependent, they can be easily used for the extraction of angle-domain common-image gathers. Moreover, the imaging condition avoids the polarization reversal in PS wave images and does not require all of the source wavefield data. Consequently, the computation and storage requirements are significantly reduced, which will facilitate the use of the elastic RTM in practice. Keywords: elastic reverse time migration; vector-based; angle-dependent; imaging condition.
ACCEPTED MANUSCRIPT
1. Introduction
PT
The theory of reverse time migration (RTM) was first proposed over forty years ago, and
RI
has since been examined by many scholars (Claerbout, 1971; Baysal, 1983; Chang and
SC
McMechan, 1987; Wapenaar and Haime, 1990). The concept of RTM involves imaging the interface using the source and receiver wavefields. Compared to acoustic RTM, elastic RTM is
NU
more representative of the actual situation and can handle multicomponent seismic data. In
MA
addition, converted-wave RTM results have higher resolution due to their shorter wavelengths compared to those of P-waves (Wang and McMechan, 2015).
D
Elastic RTM involves three main problems: elastic simulation, wavefield decomposition,
TE
and imaging conditions. The staggered grid (Virieux, 1986; Levander, 1988) and optimization of
AC CE P
the difference coefficients based on the method of least-squares (Liu et al., 2014; Cai et al., 2015) can improve the accuracy of elastic simulations. In terms of the boundary condition, Hastings et al. (1996) extended the perfectly matched layer (PML) from electromagnetic to elastic, which enables good boundary absorption with only a few grid points. In addition, the application of parallel computing technology can greatly accelerate the computation of the RTM using tools such as OpenMP, MPI, etc. (Komatitsch, 2010). Since the elastic RTM method must accommodate multicomponent data, and its wavefield contains both P- and S-waves, we need to separate the P- and S-waves before imaging to avoid crosstalk. However, traditional decomposition methods based on the divergence and curl will lose the vector information of the elastic wavefield and cause a polarization reversal in the converted wave imaging results (Sun and McMechan, 2001; Wang and McMechan, 2015; Li et
ACCEPTED MANUSCRIPT al., 2016). Many studies have been conducted to identify a solution for the polarization reversal problem (Yan and Sava, 2008; Yan and Xie, 2012; Du et al., 2015; Li et al., 2016). In contrast to
PT
traditional methods, Wang et al. (2015) tested two methods for isotropic P and S vector decomposition, both of which preserved the amplitude and phase information, and the decoupled
SC
RI
propagation resulted in faster computations and lower memory requirements. Nguyen and McMechan (2015) discussed five ways to avoid the storage of source
NU
wavefield slices, three of which were wavefield reconstruction methods that have since been
MA
adopted by many scholars (Dussaud et al., 2008; Feng et al., 2012; Sun and Fu, 2013; Liu et al., 2015). The remaining two imaging conditions, namely, the excitation amplitude (EA) imaging
D
condition and the sparse crosscorrelation imaging condition, did not require wavefield
TE
reconstruction, which greatly reduced the calculation or storage requirements. Nguyen and
AC CE P
McMechan (2013) used the EA imaging condition for acoustic RTM. Wang and McMechan (2015) extended this imaging condition to elastic RTM, and termed it the vector-based (VB) imaging condition. The procedure they established avoids polarization reversal in the PS wave imaging results, but requires the sign of the reflection coefficients to be calculated separately. In this study, we extend the EA imaging condition to elastic RTM based on the ratio of the particle velocity vector inner product, and propose a new procedure for 2D elastic multicomponent RTM. The proposed imaging condition is different from the VB imaging condition in that there is no need to separately calculate the sign of the reflection coefficients. We first present the method for elastic modelling and wavefield decomposition. Next, we derive the vector-based excitation amplitude (VEA) imaging condition and introduce a detailed process for applying the new imaging condition. Then, we show a simple layered model test and a
ACCEPTED MANUSCRIPT Marmousi II model test to demonstrate the efficiency of the new imaging condition and procedure.
RI
PT
2. Methodology
SC
2.1 Wavefield modelling and decomposition
First-order stress-velocity equations are often used for elastic simulation because they are
NU
suitable for staggered-grid finite-difference methods (Virieux, 1986; Levander, 1988). In this
MA
study, two-dimensional (2D) equations are adopted for wavefield modelling. However, to solve the equations, the spatial difference order is increased from second-order (Virieux, 1986) and
D
fourth-order (Levander, 1988) to 18th-order (Dong et al., 2000) or higher. Increasing the order
TE
reduces the spatial dispersion and allows the use of a larger time step in the simulation. However, the calculations become more complicated as the difference order increases. Considering
AC CE P
accuracy and efficiency requirements, we choose a second-order time difference and a 12th-order spatial difference. To improve the calculation accuracy, many scholars optimize the difference coefficients. For example, Xin et al. (2014) compared four different coefficient solving methods, including Taylor's expansion, least squares optimization (Yang et al., 2014), simulated annealing, and linearization, and the latter three methods provided better results than those obtained from the Taylor's expansion method. In this paper, we applied the least squares optimization method of calculating the coefficients. In addition, we also adapted the PML boundary condition and OpenMP parallel algorithm to wavefield modelling. In terms of the wavefield decomposition method, Wang et al. (2015) demonstrated that the wavefield decomposition methods used by Zhang et al. (2007) and Xiao and Leaney (2010)
ACCEPTED MANUSCRIPT are mathematically equivalent, and both can be used to obtain the decomposed P- and S-waves and preserve the vector information. The latter used a P-wave stress p to reconstruct the
PT
decoupled elastic wave equation, which is more efficient method. We therefore adopted this
RI
method in this study.
SC
2.2 The VEA imaging condition
NU
The proposed imaging condition needs to calculate the imaging time T at each grid point in the forward propagation, which is the time when the amplitude of the P-wave particle velocity
MA
vector at each grid point is at its maximum. Thus,
（1）
TE
D
p T ( x, z ) argmax[ v src x, z , t ],
AC CE P
p where v src x, z is the amplitude of the P-wave particle velocity vector, and the argmax
operator provides the imaging time T when the amplitude is at its maximum over all moments. In general, it is possible to accurately determine the imaging time T; however, when the media is very complex or the velocity model has been smoothed inappropriately, the calculated imaging time may be uncorrected. Gu et al. (2015) studied this problem, and determined that this could be addressed by restricting the imaging time calculation based on the ray travel time. After the imaging time has been calculated, the P-wave particle velocity vector and Pwave stress at this instant are saved. Next, the Poynting vector and incident angle are calculated. For an isotropic medium, since the P-wave stress is a scalar field, the P-wave Poynting vectors of the source and receiver wavefields can be calculated using the following equation,
ACCEPTED MANUSCRIPT s p x, z p x, z v p ( x, z),
（2）
PT
where s p is the P-wave Poynting vector; p is the P-wave stress; and v p is the decomposed Pwave particle velocity vector. However, the calculation of the S-wave Poynting vector is more
SC
RI
complicated (Wang and McMechan, 2015), as shown by the following:
NU
sxs x, z xx x, z p x, z vxs x, z xz x, z vzs x, z ,
（4）
MA
szs x, z zz x, z p x, z vzs x, z xz x, z vxs x, z ,
（3）
D
where vxs and vzs are the horizontal and vertical components of the decomposed S-wave particle
TE
velocity vector; xx , zz and xz are the stress tensors; and sxs and szs are the horizontal and
AC CE P
vertical components of the S-wave Poynting vector. Then, the incident angle at each grid point can be calculated based on the geometric relationship (Wang and McMechan, 2015). When the interface is tilted (Figure 1), we propose the new incident angle calculation method shown in Eq. 5, in which the left incident angle is a positive value and the right incident angle is a negative value.
szp, src x, z szp,rec x, z, T x, z, T 0.5 actan p actan p , s x , z s x , z , T x , src x , rec
(5)
where is the incident angle; szp,src and sxp,src are the vertical and horizontal components of the P-wave Poynting vector of the source wavefield at imaging time T; and the subscript ‘rec’ denotes the corresponding components of the receiver wavefield.
SC
RI
PT
ACCEPTED MANUSCRIPT
Figure 1. Geometric relationship between the incident angle and Poynting vectors of the source
NU
p p ( s src ) and receiver ( s rec ) wavefields.
MA
Similar to the EA imaging condition, the reflection coefficient can be obtained by dividing the receiver wavefield by the source wavefield at the imaging time. For elastic media,
TE
AC CE P
the vector modes,
D
both the receiver and source wavefields are vectors, and the reflection coefficient is the ratio of
rpp ( x, z , )
rps ( x, z , )
p v rec x, z , T p v src x, z
s v rec x, z , T p v src x, z
（6）
,
,
（7）
where rpp and rps are the calculated reflection coefficients of the PP and PS waves, respectively. p The saved P-wave particle velocity vector is represented by v src , and the decomposed P- and S-
p s wave particle velocity vectors at imaging time T are represented by v rec and v rec , respectively.
Since Equations 6 and 7 are module ratios, only the magnitude of the particle vibration is considered, and the direction of the particle vibration is not taken into account. Consequently, the
ACCEPTED MANUSCRIPT calculated reflection coefficients are always positive and do not fit the actual situation. In this paper, we propose a new VEA imaging condition,
PT
p v src x, z vrecp x, z, T , p v src x, z v srcp x, z
（8）
RI
rpp ( x, z, )
SC
p s v src x, z vrec x, z, T . p p v src x, z v src x, z
（9）
NU
rps ( x, z, )
MA
Equations 8 and 9 replace the vector module ratio with the vector inner product ratio, in which the particle vibration direction is taken into account. However, there is a problem in
TE
D
Equations. 8 and 9. Taking the PP imaging condition as an example,
AC CE P
p p p x, z, T cos 2 v src x, z v rec x, z, T v srcp x, z v rec 2 p v src x, z v srcp x, z v P x, z
.or.
src
p p v src x, z v rec x, z, T cos 2
v
p src
x, z
2
. （10）
Compared with Eq. 6, Eq. 10 is multiplied by cos 2 , which results in a deviation from the calculated result as the incident angle increases. However, this also has a benefit. As shown in Figure 2, the incident angle is 40°. Taking the particle vibration direction as an example, which is the same as the direction of wave propagation, when the reflection coefficient is positive, the corresponding particle vibration direction of the reflected P-wave is similar to that represented by the solid arrow; however, when the reflection coefficient is negative, the corresponding particle vibration direction of the reflected P-wave is similar to that indicated by the dotted arrow. It can be seen that the sign of the reflection coefficient calculated using Eq. 10
ACCEPTED MANUSCRIPT is opposite to that of the true reflection coefficient. These two problems can be solved simultaneously. Thus,
PT
p v p x, z v rec x, z, T , 1 srcp p cos 2 v src x, z v src x, z
（11）
RI
rpp ( x, z, )
SC
where 45 . If theta equals 45 degrees, the denominator cos 2 will be equal to zero, in
NU
which case, the equation is meaningless. In practice, the incident angle used in the seismic migration is often less than 40 degrees. In theory, when the incident angle is essentially 45
MA
degrees, we can approximate the denominator by replacing it with 0.001 to avoid any instability. With this modification, the angle limit no longer poses a significant problem. Equation 11 can be
AC CE P
TE
D
used to directly obtain the signed accurate reflection coefficient.
Figure 2. Relationship between the P-wave particle vibration direction and PP wave reflection coefficient.
ACCEPTED MANUSCRIPT Figure 3. Geometric relationship between the vectors and angles. The reflection angle of the converted S-wave is represented by ' , and the angle between the incident P-wave particle
PT
velocity vector and the converted S-wave particle velocity vector is represented by .
RI
Similarly, we can modify Eq. 9 to obtain:
(12)
NU
SC
p s v src x, z v rec x, z, T 1 . rps ( x, z, ) p p v src x, z v src x, z cos
MA
According to the geometric relationship shown in Figure 3, and referring to Eq. 5, szp, src x, z szs,rec x, z , T actan p actan s , 2 sx , src x, z sx ,rec x, z , T
(13)
TE
D
AC CE P
where szs,rec and sxs,rec are the vertical and horizontal components of the decomposed S-wave Poynting vector at imaging time T. Since the denominator is cos , it cannot be zero, that is
0 . In theory, when the incident angle is vertical, a converted S-wave is not produced, and when the incident angle is very small, the converted shear wave energy is extremely weak. Equations 11 and 12 are complete representations of the PP and PS wave VEA imaging condition, respectively. When multi-shot RTM is applied using the VEA imaging condition, each grid point in each shot migration result corresponds to a reflection coefficient and incident angle. Then, the reflection coefficient of each grid point is sorted according to the angle of incidence and angle-domain common-image gathers (ADCIGs) are I x, , z .
ACCEPTED MANUSCRIPT For the source-normalized crosscorrelation (SNC) imaging condition, the application to elastic RTM based on vector decomposition also involves vector issues. Referring to Equations
PT
11 and 12, the vector-based source-normalized crosscorrelation (VSNC) imaging conditions of
RI
the PP and PS waves can be represented as,
nt
1 cos 2
p src
p ( x, z, t ) v rec ( x, z, t )
0
nt
v
p src
nt
s ( x, z, t ) v rec ( x, z, t )
0
nt
MA
I PS x, z
p src
( x, z , t )
,
(14)
2
NU
0
v
SC
I PP x, z
v
v
p src
( x, z, t )
2
1 , cos
(15)
TE
D
0
where I pp and I ps are the imaging results of the PP and PS waves, respectively. However, one
AC CE P
difficulty in the VSNC imaging condition is that it is not easy to determine the corresponding incident angle because the imaging value of each grid point is the sum of the imaging results at all times. In this paper, in order to compare the imaging results of the VEA and VSNC imaging conditions, the incident angle in the VSNC imaging condition also accepts the incident angle of the former. To summarize, a flow chart of the new imaging method is shown in Figure 4. The source and receiver wavefields are decomposed using the decoupled propagation. The P-wave particle velocity vector, P-wave stress, and imaging time of the source wavefield are saved. In the process of reverse time propagation, PP images can be obtained at imaging time with the PP imaging condition, and PS images can be similarly obtained. At the same time, the incident angle
ACCEPTED MANUSCRIPT and angle are calculated using the corresponding Poynting vectors. After the image of each grid point is sorted according to the incident angle, the ADCIGs are obtained. To reduce the
PT
artifacts, the PP and PS images are stacked with the restriction of the incident angle. Note that
RI
because model smoothing has been discussed in many articles about RTM, we do not provide the details here. Before RTM processing, the velocity models tested in this article were properly
AC CE P
TE
D
MA
NU
SC
smoothed.
Figure 4. Flow chart of elastic RTM with the VEA imaging condition.
3. Model test 3.1 Layered model test
ACCEPTED MANUSCRIPT To test the procedure, we first designed a simple three-layer model (Figure 5). The grid spacing in two directions was 5 m, and the PML boundary condition was added with 10 grid
PT
points. The source wavelet was a Ricker wavelet with a main frequency of 30 Hz. Before applying elastic RTM, we first tested the efficiency of the wavefield decomposition method.
RI
Based on the results shown in Figure 6, the wavefield decomposition method is valid and the
TE
D
MA
NU
SC
Poynting vectors calculated using the P- and S-wave are acceptable.
AC CE P
Figure 5. Parameters of the layered model
AC CE P
TE
D
MA
NU
SC
RI
PT
ACCEPTED MANUSCRIPT
Figure 6. Pictures of the x- (a) and z-components (b) of the particle velocity; x- (c) and zcomponents (d) of the decomposed P-wave; x- (e) and z-components (f) of the decomposed Swave, respectively. Panels (g and h) show the calculated Poynting vectors of the P- and S-waves, respectively. The single-shot migration results are shown in Figure 7. The black boxes in the figures show the effective migration results, which indicates that the migration aperture of the PS wave
ACCEPTED MANUSCRIPT is much larger than that of the PP wave. Note that the PS images are poorly illuminated near the centre of the panels due to the weak energy of the converted wave when the incident angle was
PT
small. These results demonstrate that the VEA and VSNC imaging conditions yield imaging results similar to the EA and SNC conditions, and that the directly calculated reflection
RI
coefficients are representative. In addition, it can be seen that the PS wave migration result with
AC CE P
TE
D
MA
NU
SC
the VEA imaging condition has the highest spatial resolution.
ACCEPTED MANUSCRIPT Figure 7. Single-shot migration results calculated using the EA, VEA, SNC, and VSNC imaging conditions are represented in four columns from top to bottom, respectively. The left are the PP
PT
wave images and the right are the PS wave images.
RI
Figure 8 provides a comparison of the first layer reflection coefficients obtained using the
SC
VEA and VSNC imaging conditions and the Zoeppritz equation. The PP wave reflection coefficients obtained from the two imaging conditions agree with that calculated from Zoeppritz
NU
theory when the offset is between 0–0.7 km (corresponding to an incidence angle of 0–35°). The
MA
PS wave reflection coefficients are in good agreement with the Zoeppritz theory when the offset
AC CE P
TE
D
is between 0.1 km and 1.5 km (corresponding to an incidence angle of 5–56°).
Figure 8. Comparison of the first layer reflection coefficients obtained from the two imaging conditions and the Zoeppritz equation. Panel (a) shows the PP wave reflection coefficients, and panel (b) shows the PS wave reflection coefficients. 3.2 Marmousi II test
ACCEPTED MANUSCRIPT We also performed a complex model test based on the Marmousi II model (Figure 9). The grid spacing in both directions was 5 m, the time sampling interval was 0.5 ms, and the total
PT
recording time was 4 s. A total of 101 shots were distributed over the surface with 50 m spacing. A total of 1,001 receivers were located on the surface. The water layer was replaced by a rock
RI
layer with a non-zero S-wave velocity. The wavefield modelling method and boundary were the
AC CE P
TE
D
MA
NU
SC
same as previous, and the source frequency was reduced to 20 Hz.
Figure 9. A portion of the Marmousi II model. Panels a, b, and c represent the P-wave velocity, S-wave velocity, and density distributions, respectively.
ACCEPTED MANUSCRIPT The ADCIGs extracted from the PP and PS wave migration results based on the VEA imaging condition are shown in Figures 10a and 10c (only the angle gathers corresponding to the
PT
black dashed line in Figure 9a are shown). From the figures, many reflection layers can be distinguished. Compared with PP ADCIGs, the PS ADCIGs have a higher resolution and can
RI
therefore provide a more detailed depiction. As mentioned earlier, each grid point corresponds to
SC
only one reflection coefficient and incident angle pair in each shot migration result.
NU
Consequently, the number of effective reflection coefficients of ADCIG for each grid point is less than the number of shots, which results in some discrete points in the ADCIGs rather than
MA
continuous events. However, as the number of shots increases, the events in the ADCIGs will become clearer. If a velocity analysis is performed using ADCIGs, continuous events may be
D
necessary, which would require a sufficient number of shots or interpolating ADCIGs. In
TE
addition, when the velocity model has –10% errors, the events shown in Figures 10b and 10d are
AC CE P
no longer horizontal, but curve upward.
AC CE P
TE
D
MA
NU
SC
RI
PT
ACCEPTED MANUSCRIPT
Figure 10. The columns in (a) and (c) are PP and PS ADCIGs with the true velocity model, respectively. The columns in (b) and (d) are similar to those in (a) and (c), except that the velocity model has –10% errors.
Figure 11 show the stacked results of 101 shots calculated using the VEA and VSNC imaging conditions with the true velocity model. Overall, these imaging results are very clear and many details of the subsurface structure are recognizable. By comparing the detailed areas (indicated by the black arrows), the results obtained with the VEA imaging condition have a higher imaging resolution than those obtained with the VSNC imaging condition. The PS wave migration results are imaged with higher resolution and richer imaging detail than the PP wave migration results.
AC CE P
TE
D
MA
NU
SC
RI
PT
ACCEPTED MANUSCRIPT
ACCEPTED MANUSCRIPT Figure 11. Panels (a) and (b) are the stacked results of PP wave imaging with the VEA and VSNC imaging conditions, respectively. Panels (c) and (d) are the stacked results of PS wave
PT
imaging with the VEA and VSNC imaging conditions, respectively.
RI
4. Discussion
SC
The VEA imaging condition is an extension of the EA imaging condition. First, only the
NU
P-wave particle velocity vector and corresponding P-wave stress are saved at the imaging time of each grid point in the forward propagation, instead of the complete wavefield at each moment.
MA
Compared with the traditional crosscorrelation imaging conditions and wavefield reconstruction method, this approach greatly reduces the storage and computing requirements and allows the
D
elastic RTM to be applied to actual data. In addition, only the maximum amplitude is used for
TE
imaging, which is the largest contribution to the image, and results in a higher resolution. Second,
AC CE P
as was the case in the research of Wang and McMechan (2015), the wavefield decomposition method proposed in this paper is based on vector decomposition. The PS wave images do not involve polarization reversal, so multi-shot imaging results can be directly stacked without the need for polarization correction.
Although there are some common points between the imaging conditions, there are also several differences. First, the VEA imaging condition directly includes information regarding the particle vibration directions in the process of imaging, which means that the signed angledependent reflection coefficients can be directly obtained without calculating the sign. Second, a new equation to calculate the incident angle is proposed that can distinguish between the left or right incidence. Third, the VSNC imaging condition is deduced by extending the concept to the SNC imaging condition.
ACCEPTED MANUSCRIPT In addition, there are some things that need to be explained. First, the ADCIGs can be extracted directly from the multi-shot imaging results and then used in the velocity analyses. At
PT
the same time, the ADCIGs can also constrain the superposition process and allow the results to be superimposed with the effective imaging angle to reduce the imaging artifacts. Second, the
RI
process of wave propagation and imaging uses the OpenMP parallel algorithm. In actual data
SC
calculations or three-dimensional (3D) RTM calculations, the MPI parallel technique can also be
NU
used, which will greatly improve the computational efficiency.
MA
Although the methodology and tests shown above are about 2D isotropic media, it is possible to extend them to 2D anisotropic media, 3D isotropic media, and attenuation. If the 2D
D
model contains anisotropic media, taking VTI media as an example, the imaging conditions can
TE
also be applied after separating the P- and S-waves well (Xiao and Leaney, 2010; Wang et al.,
AC CE P
2015). For 3D, if the medium is isotropic, the particle vibration directions of the incident P-wave, the reflected P- and S-waves are in the same plane, except that the particle velocity vector contains three components. On this basis, it is easy to extend this method to 3D isotropic media. However, for anisotropic media, there is a change in the azimuth angle, and the reflected S-wave may be split into SV and SH waves. Thus, this method cannot be simply extended to 3D anisotropic media, which requires more studies. The maximum amplitude is used to calculate the correct imaging time. The imaging time is the time when the particle amplitude of each grid point at its largest over all times, similar to the one-way time from the source to that grid point. When the attenuative anomalies exist, although the amplitude will be attenuated, the amplitude of the wavefront is still the largest over all times. Therefore, we infer that the attenuative anomalies do not affect the calculation of the imaging time, and it will not affect the application of this method.
ACCEPTED MANUSCRIPT Finally, the multipath problem cannot be solved (Wang and McMechen, 2015). In the original velocity model, there is a low-velocity body (marked with the hollow arrow in Figure
PT
9a). The P-wave velocity of the surrounding rocks is about 1,700 m/s and the P-wave velocity of the low-velocity body is about 1,100 m/s. There will be strong interface multiples in this low-
RI
velocity body, and a strong diffraction wave at both ends. The imaging results will appear in a
SC
strong illusion under the low-velocity body (Du et al., 2015). Therefore, we set the P-wave
NU
velocity of this low-velocity body to 1,500 m/s, so that the difference between the surrounding rocks was not very large. Consequently, there is no obvious illusion in the imaging results.
MA
Another noteworthy problem is the jaggedness in the PS imaging results (Figures 11c and 11d) when the interface is tilted, which also occurs in the original velocity model due to the regular
AC CE P
5. Conclusion
TE
D
grid. This shows that the PS imaging results have a high resolution.
In this paper, we propose a new process for 2D elastic RTM. The decomposed P- and Swaves are obtained by the decoupled propagation of the source and receiver wavefields. In the process of forward modelling, the imaging time of each grid point is calculated, and the corresponding P-wave particle velocity vector and P-wave stress are saved. During the process of reverse time propagation, each grid point is imaged using the VEA imaging condition at the imaging time. The P- and S-wave Poynting vectors are calculated to obtain the incident and reflection angles, which are then used to extract the ADCIGs and amplitude corrections. The migration results of the simple model and the Marmousi II model show that the VEA imaging condition can be used to obtain signed angle-dependent reflection coefficients directly with high computational efficiency and strong adaptability. The PS wave images have high imaging
ACCEPTED MANUSCRIPT resolution and rich imaging details. In addition, by extending the concept to the SNC imaging condition, the VSNC imaging condition is deduced and the synthetic results prove its validity.
PT
Acknowledgments
RI
The research work in this paper was supported by the Natural Science Foundation of
SC
China (NSFC) [Grant No. 41374108] and Major Projects of the National Science and
AC CE P
TE
D
MA
NU
Technology of China [Grant No. 2016ZX05026-002-003].
ACCEPTED MANUSCRIPT References Baysal, E., Kosloff, D., Sherwood, J., 1983. Reverse time migration. Geophysics 48 (11),
RI
PT
1514-1524, doi: 10.1190/1.1441434.
Claerbout, J.F., 1971. Towards a unified theory of reflector mapping. Geophysics 36 (3),
SC
467-481, doi: 10.1190/1.1440185.
NU
Chang, W., McMechan, G.A., 1987. Elastic reverse-time migration. Geophysics 52 (10),
MA
1365-1375, doi: 10.1190/1.1442249.
Cai, X., Liu, Y., Ren, Z., 2015. Elastic Reverse-Time Migration Based on the Optimal
AC CE P
TE
10.1190/segam2015-5906925.1.
D
Staggered-Grid Finite-Difference Scheme. SEG, Expanded Abstracts, 4148-4152, doi:
Dong, L.G., Ma, Z.T., Cao, J.Z., Wang, H.Z., Geng, J.H., Lei, B., 2000. A staggered-grid high-order difference method of one-order elastic wave equation. Chinese Journal of Geophysics 43 (3), 411-419.
Dussaud, E., Symes, W.W., Williamson, P., Lemaistre, L., Singer, P., Denel, B., 2008. Computational strategies for reverse-time migration. SEG, Expanded Abstracts, 2267-2271, doi: 10.1190/1.3059336. Du, Q., Zhang, M., Gong, X., Chen, X., 2015. Polarity-consistent excitation amplitude imaging condition for elastic reverse time migration. Journal of Geophysics and Engineering 12 (1), 33-44, doi: 10.1088/1742-2132/12/1/33.
ACCEPTED MANUSCRIPT Gu, B.L., Liu, Y.S., Ma, X.N., Li, Z.Y., Liang, G.H., 2015. A modified excitation amplitude imaging condition for prestack reverse time migration. Exploration Geophysics 46 (4),
PT
359-370, doi: 10.1071/EG14039.
RI
Hastings, F.D., Schneider, J. B., Broschat, S. L., 1996. Application of the perfectly
SC
matched layer (PML) absorbing boundary conditions to elastic wave propagation. Journal of the
NU
Acoustical Society of America 100 (5), 3061-3069, doi: 10.1121/1.417118. Komatitsch, D., Goddeke, D., Erlebacher, G., Michea, D., 2010. Modeling the
MA
propagation of elastic waves using spectral elements on a cluster of 192 GPUs. Computer
D
Science - Research and Development 25, 75-82, doi: 10.1007/s00450-010-0109-1.
TE
Levander, A.R., 1988. Fourth-Order Finite-Difference P-SV Seismograms. Geophysics,
AC CE P
53 (11), 1425-1436, doi: 10.1190/1.1442422. Liu, Y., 2013. Globally optimal finite-difference schemes based on least squares. Geophysics 78 (4), T113-T132, doi: 10.1190/geo2012-0480.1. Liu, S., Li, X., Wang, W., Zhu, T., 2015. Source wavefield reconstruction using a linear combination of the boundary wavefield in reverse time migration. Geophysics 80 (6), S203-S212. doi: 10.1190/geo2015-0109.1. Li, Z., Ma, X., Fu, C., Liang, G., 2016. Wavefield separation and polarity reversal correction in elastic reverse time migration. Journal of Applied Geophysics 127, 56-67, doi: 10.1016/j.jappgeo.2016.02.012.
ACCEPTED MANUSCRIPT Nguyen, B.D., McMechan, G. A., 2013. Excitation amplitude imaging condition for prestack reverse-time migration. Geophysics 78 (1), S37-S46, doi: 10.1190/geo2012-0079.1.
PT
Nguyen, B.D., McMechan, G. A., 2015. Five ways to avoid storing source wavefield
RI
snapshots in 2d elastic prestack reverse time migration, Geophysics 80 (1), S1-S18, doi:
SC
10.1190/geo2014-0014.1.
NU
Sun, R., McMechan, G. A., 2001. Scalar reverse-time depth migration of prestack elastic
MA
seismic data. Geophysics 66 (5), 1519-1527, doi: 10.1190/1.1487098. Sun, W., Fu, L.Y., 2013. Two effective approaches to reduce data storage in reverse time
D
migration. Computers and Geosciences 56, 69-75, doi: 10.1016/j.cageo.2013.03.013.
TE
Tang, C., Wang, D., 2012. Reverse Time Migration with Source Wavefield
AC CE P
Reconstruction in Target Imaging Region. EAGE, Extended Abstracts, P272, doi: 10.3997/22144609.20148620.
Virieux, J., 1986. P-SV wave propagation in heterogeneous media: velocity-stress finitedifference method, Geophysics 51 (4), 889–901. doi: 10.1190/1.1442147. Wapenaar, C.P., Haime, G.C., 1990. Elastic extrapolation of primary seismic p- and swaves, Geophysical Prospecting 38 (1), 23-60. doi: 10.1111/j.1365-2478.1990.tb01833.x. Wang, W., McMechan, G.A., Zhang, Q., 2015. Comparison of two algorithms for isotropic elastic P and S vector decomposition. Geophysics 80 (4), T147-T160, doi: 10.1190/geo2014-0563.1.
ACCEPTED MANUSCRIPT Wang, W., McMechan, G.A., 2015. Vector-based elastic reverse time migration. Geophysics 80 (6), S245-S258, doi: 10.1190/geo2014-0620.1.
PT
Xiao, X., Leaney, W.S., 2010. Local vertical seismic profiling (VSP) elastic reverse-time
SC
Geophysics 75 (2), S35-S49, doi: 10.1190/1.3309460.
RI
migration and migration resolution: salt-flank imaging with transmitted P-to-s waves.
NU
Xin, W., Yan, Z.C., Liang, W.Q., Chen, Y.H., Yang, C.C., 2015. Methods to determine the finite difference coefficients for elastic wave equation modeling. Chinese Journal of
MA
Geophysics- Chinese Edition 58 (7), 2486-2495, doi: 10.6038/cjg20150724.
D
Yan, J., Sava, P., 2008. Isotropic angle-domain elastic reverse-time migration.
TE
Geophysics 73 (6), S229-S239, doi: 10.1190/1.2981241.
AC CE P
Yan, R., Xie, X., 2012. An angle-domain imaging condition for elastic reverse time migration and its application to angle gather extraction. Geophysics 77 (5), S105-S115, doi: 10.1190/geo2011-0455.1.
Yang, L., Yan, H., Liu, H., 2014. Least squares staggered-grid finite-difference for elastic wave modeling. Exploration Geophysics 45 (4), 255-260, doi: 10.1071/EG13087. Zhang, J., Tian, Z., Wang, C., 2007. P- and S-wave-separated elastic wave-equation numerical modeling using 2D staggered grid. SEG, Expanded Abstracts, 2104-2109, doi: 10.1190/1.2792904.
ACCEPTED MANUSCRIPT
AC CE P
TE
D
MA
NU
SC
RI
PT
Highlights: - A new method for 2D elastic multicomponent reverse time migration (RTM) is proposed. - The vector-based excitation amplitude imaging condition is deduced. - The signed angle-dependent reflection coefficients are directly obtained. - The concept is extended to the source-normalized crosscorrelation imaging condition.
S0926-9851(17)30377-4 doi:10.1016/j.jappgeo.2017.10.003 APPGEO 3349
To appear in:
Journal of Applied Geophysics
Received date: Revised date: Accepted date:
16 April 2017 10 July 2017 10 October 2017
Please cite this article as: Zhou, Jinju, Wang, Deli, Vector-based excitation amplitude imaging condition for elastic RTM, Journal of Applied Geophysics (2017), doi:10.1016/j.jappgeo.2017.10.003
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT
Vector-based excitation amplitude imaging condition for elastic RTM
SC
130026, China
RI
College of Geoexploration Science and Technology, Jilin University, Changchun
TE
D
MA
NU
E-mail: [email protected]; [email protected]
AC CE P
1
PT
Jinju Zhou1, Deli Wang1
ACCEPTED MANUSCRIPT
ABSTRACT In recent years, many studies have focused on elastic reverse time migration (RTM). In
PT
response to the problems associated with elastic RTM, we propose a new procedure for 2D
RI
elastic multicomponent RTM. In this new method, decomposed P- and S-wave components are
SC
obtained from the decoupled propagation of the source and receiver wavefields, which allows the expedient calculation of the Poynting vectors and the incident and reflection angles of the P- and
NU
S-waves. In addition, we deduce the vector-based excitation amplitude imaging condition. This
MA
process automatically accounts for the particle vibration directions when determining the angledependent signed reflection coefficients, and does not require the sign to be determined apart
D
from the value of the reflection coefficients. This concept was further extended to the source-
TE
normalized crosscorrelation imaging condition. The reflection coefficient of the layered model
AC CE P
test was in agreement with the Zoeppritz theory, the PP and PS wave images of the Marmousi II model were clear, and the PS wave images had higher resolution and richer details. In addition, since the calculated reflection coefficients are angle-dependent, they can be easily used for the extraction of angle-domain common-image gathers. Moreover, the imaging condition avoids the polarization reversal in PS wave images and does not require all of the source wavefield data. Consequently, the computation and storage requirements are significantly reduced, which will facilitate the use of the elastic RTM in practice. Keywords: elastic reverse time migration; vector-based; angle-dependent; imaging condition.
ACCEPTED MANUSCRIPT
1. Introduction
PT
The theory of reverse time migration (RTM) was first proposed over forty years ago, and
RI
has since been examined by many scholars (Claerbout, 1971; Baysal, 1983; Chang and
SC
McMechan, 1987; Wapenaar and Haime, 1990). The concept of RTM involves imaging the interface using the source and receiver wavefields. Compared to acoustic RTM, elastic RTM is
NU
more representative of the actual situation and can handle multicomponent seismic data. In
MA
addition, converted-wave RTM results have higher resolution due to their shorter wavelengths compared to those of P-waves (Wang and McMechan, 2015).
D
Elastic RTM involves three main problems: elastic simulation, wavefield decomposition,
TE
and imaging conditions. The staggered grid (Virieux, 1986; Levander, 1988) and optimization of
AC CE P
the difference coefficients based on the method of least-squares (Liu et al., 2014; Cai et al., 2015) can improve the accuracy of elastic simulations. In terms of the boundary condition, Hastings et al. (1996) extended the perfectly matched layer (PML) from electromagnetic to elastic, which enables good boundary absorption with only a few grid points. In addition, the application of parallel computing technology can greatly accelerate the computation of the RTM using tools such as OpenMP, MPI, etc. (Komatitsch, 2010). Since the elastic RTM method must accommodate multicomponent data, and its wavefield contains both P- and S-waves, we need to separate the P- and S-waves before imaging to avoid crosstalk. However, traditional decomposition methods based on the divergence and curl will lose the vector information of the elastic wavefield and cause a polarization reversal in the converted wave imaging results (Sun and McMechan, 2001; Wang and McMechan, 2015; Li et
ACCEPTED MANUSCRIPT al., 2016). Many studies have been conducted to identify a solution for the polarization reversal problem (Yan and Sava, 2008; Yan and Xie, 2012; Du et al., 2015; Li et al., 2016). In contrast to
PT
traditional methods, Wang et al. (2015) tested two methods for isotropic P and S vector decomposition, both of which preserved the amplitude and phase information, and the decoupled
SC
RI
propagation resulted in faster computations and lower memory requirements. Nguyen and McMechan (2015) discussed five ways to avoid the storage of source
NU
wavefield slices, three of which were wavefield reconstruction methods that have since been
MA
adopted by many scholars (Dussaud et al., 2008; Feng et al., 2012; Sun and Fu, 2013; Liu et al., 2015). The remaining two imaging conditions, namely, the excitation amplitude (EA) imaging
D
condition and the sparse crosscorrelation imaging condition, did not require wavefield
TE
reconstruction, which greatly reduced the calculation or storage requirements. Nguyen and
AC CE P
McMechan (2013) used the EA imaging condition for acoustic RTM. Wang and McMechan (2015) extended this imaging condition to elastic RTM, and termed it the vector-based (VB) imaging condition. The procedure they established avoids polarization reversal in the PS wave imaging results, but requires the sign of the reflection coefficients to be calculated separately. In this study, we extend the EA imaging condition to elastic RTM based on the ratio of the particle velocity vector inner product, and propose a new procedure for 2D elastic multicomponent RTM. The proposed imaging condition is different from the VB imaging condition in that there is no need to separately calculate the sign of the reflection coefficients. We first present the method for elastic modelling and wavefield decomposition. Next, we derive the vector-based excitation amplitude (VEA) imaging condition and introduce a detailed process for applying the new imaging condition. Then, we show a simple layered model test and a
ACCEPTED MANUSCRIPT Marmousi II model test to demonstrate the efficiency of the new imaging condition and procedure.
RI
PT
2. Methodology
SC
2.1 Wavefield modelling and decomposition
First-order stress-velocity equations are often used for elastic simulation because they are
NU
suitable for staggered-grid finite-difference methods (Virieux, 1986; Levander, 1988). In this
MA
study, two-dimensional (2D) equations are adopted for wavefield modelling. However, to solve the equations, the spatial difference order is increased from second-order (Virieux, 1986) and
D
fourth-order (Levander, 1988) to 18th-order (Dong et al., 2000) or higher. Increasing the order
TE
reduces the spatial dispersion and allows the use of a larger time step in the simulation. However, the calculations become more complicated as the difference order increases. Considering
AC CE P
accuracy and efficiency requirements, we choose a second-order time difference and a 12th-order spatial difference. To improve the calculation accuracy, many scholars optimize the difference coefficients. For example, Xin et al. (2014) compared four different coefficient solving methods, including Taylor's expansion, least squares optimization (Yang et al., 2014), simulated annealing, and linearization, and the latter three methods provided better results than those obtained from the Taylor's expansion method. In this paper, we applied the least squares optimization method of calculating the coefficients. In addition, we also adapted the PML boundary condition and OpenMP parallel algorithm to wavefield modelling. In terms of the wavefield decomposition method, Wang et al. (2015) demonstrated that the wavefield decomposition methods used by Zhang et al. (2007) and Xiao and Leaney (2010)
ACCEPTED MANUSCRIPT are mathematically equivalent, and both can be used to obtain the decomposed P- and S-waves and preserve the vector information. The latter used a P-wave stress p to reconstruct the
PT
decoupled elastic wave equation, which is more efficient method. We therefore adopted this
RI
method in this study.
SC
2.2 The VEA imaging condition
NU
The proposed imaging condition needs to calculate the imaging time T at each grid point in the forward propagation, which is the time when the amplitude of the P-wave particle velocity
MA
vector at each grid point is at its maximum. Thus,
（1）
TE
D
p T ( x, z ) argmax[ v src x, z , t ],
AC CE P
p where v src x, z is the amplitude of the P-wave particle velocity vector, and the argmax
operator provides the imaging time T when the amplitude is at its maximum over all moments. In general, it is possible to accurately determine the imaging time T; however, when the media is very complex or the velocity model has been smoothed inappropriately, the calculated imaging time may be uncorrected. Gu et al. (2015) studied this problem, and determined that this could be addressed by restricting the imaging time calculation based on the ray travel time. After the imaging time has been calculated, the P-wave particle velocity vector and Pwave stress at this instant are saved. Next, the Poynting vector and incident angle are calculated. For an isotropic medium, since the P-wave stress is a scalar field, the P-wave Poynting vectors of the source and receiver wavefields can be calculated using the following equation,
ACCEPTED MANUSCRIPT s p x, z p x, z v p ( x, z),
（2）
PT
where s p is the P-wave Poynting vector; p is the P-wave stress; and v p is the decomposed Pwave particle velocity vector. However, the calculation of the S-wave Poynting vector is more
SC
RI
complicated (Wang and McMechan, 2015), as shown by the following:
NU
sxs x, z xx x, z p x, z vxs x, z xz x, z vzs x, z ,
（4）
MA
szs x, z zz x, z p x, z vzs x, z xz x, z vxs x, z ,
（3）
D
where vxs and vzs are the horizontal and vertical components of the decomposed S-wave particle
TE
velocity vector; xx , zz and xz are the stress tensors; and sxs and szs are the horizontal and
AC CE P
vertical components of the S-wave Poynting vector. Then, the incident angle at each grid point can be calculated based on the geometric relationship (Wang and McMechan, 2015). When the interface is tilted (Figure 1), we propose the new incident angle calculation method shown in Eq. 5, in which the left incident angle is a positive value and the right incident angle is a negative value.
szp, src x, z szp,rec x, z, T x, z, T 0.5 actan p actan p , s x , z s x , z , T x , src x , rec
(5)
where is the incident angle; szp,src and sxp,src are the vertical and horizontal components of the P-wave Poynting vector of the source wavefield at imaging time T; and the subscript ‘rec’ denotes the corresponding components of the receiver wavefield.
SC
RI
PT
ACCEPTED MANUSCRIPT
Figure 1. Geometric relationship between the incident angle and Poynting vectors of the source
NU
p p ( s src ) and receiver ( s rec ) wavefields.
MA
Similar to the EA imaging condition, the reflection coefficient can be obtained by dividing the receiver wavefield by the source wavefield at the imaging time. For elastic media,
TE
AC CE P
the vector modes,
D
both the receiver and source wavefields are vectors, and the reflection coefficient is the ratio of
rpp ( x, z , )
rps ( x, z , )
p v rec x, z , T p v src x, z
s v rec x, z , T p v src x, z
（6）
,
,
（7）
where rpp and rps are the calculated reflection coefficients of the PP and PS waves, respectively. p The saved P-wave particle velocity vector is represented by v src , and the decomposed P- and S-
p s wave particle velocity vectors at imaging time T are represented by v rec and v rec , respectively.
Since Equations 6 and 7 are module ratios, only the magnitude of the particle vibration is considered, and the direction of the particle vibration is not taken into account. Consequently, the
ACCEPTED MANUSCRIPT calculated reflection coefficients are always positive and do not fit the actual situation. In this paper, we propose a new VEA imaging condition,
PT
p v src x, z vrecp x, z, T , p v src x, z v srcp x, z
（8）
RI
rpp ( x, z, )
SC
p s v src x, z vrec x, z, T . p p v src x, z v src x, z
（9）
NU
rps ( x, z, )
MA
Equations 8 and 9 replace the vector module ratio with the vector inner product ratio, in which the particle vibration direction is taken into account. However, there is a problem in
TE
D
Equations. 8 and 9. Taking the PP imaging condition as an example,
AC CE P
p p p x, z, T cos 2 v src x, z v rec x, z, T v srcp x, z v rec 2 p v src x, z v srcp x, z v P x, z
.or.
src
p p v src x, z v rec x, z, T cos 2
v
p src
x, z
2
. （10）
Compared with Eq. 6, Eq. 10 is multiplied by cos 2 , which results in a deviation from the calculated result as the incident angle increases. However, this also has a benefit. As shown in Figure 2, the incident angle is 40°. Taking the particle vibration direction as an example, which is the same as the direction of wave propagation, when the reflection coefficient is positive, the corresponding particle vibration direction of the reflected P-wave is similar to that represented by the solid arrow; however, when the reflection coefficient is negative, the corresponding particle vibration direction of the reflected P-wave is similar to that indicated by the dotted arrow. It can be seen that the sign of the reflection coefficient calculated using Eq. 10
ACCEPTED MANUSCRIPT is opposite to that of the true reflection coefficient. These two problems can be solved simultaneously. Thus,
PT
p v p x, z v rec x, z, T , 1 srcp p cos 2 v src x, z v src x, z
（11）
RI
rpp ( x, z, )
SC
where 45 . If theta equals 45 degrees, the denominator cos 2 will be equal to zero, in
NU
which case, the equation is meaningless. In practice, the incident angle used in the seismic migration is often less than 40 degrees. In theory, when the incident angle is essentially 45
MA
degrees, we can approximate the denominator by replacing it with 0.001 to avoid any instability. With this modification, the angle limit no longer poses a significant problem. Equation 11 can be
AC CE P
TE
D
used to directly obtain the signed accurate reflection coefficient.
Figure 2. Relationship between the P-wave particle vibration direction and PP wave reflection coefficient.
ACCEPTED MANUSCRIPT Figure 3. Geometric relationship between the vectors and angles. The reflection angle of the converted S-wave is represented by ' , and the angle between the incident P-wave particle
PT
velocity vector and the converted S-wave particle velocity vector is represented by .
RI
Similarly, we can modify Eq. 9 to obtain:
(12)
NU
SC
p s v src x, z v rec x, z, T 1 . rps ( x, z, ) p p v src x, z v src x, z cos
MA
According to the geometric relationship shown in Figure 3, and referring to Eq. 5, szp, src x, z szs,rec x, z , T actan p actan s , 2 sx , src x, z sx ,rec x, z , T
(13)
TE
D
AC CE P
where szs,rec and sxs,rec are the vertical and horizontal components of the decomposed S-wave Poynting vector at imaging time T. Since the denominator is cos , it cannot be zero, that is
0 . In theory, when the incident angle is vertical, a converted S-wave is not produced, and when the incident angle is very small, the converted shear wave energy is extremely weak. Equations 11 and 12 are complete representations of the PP and PS wave VEA imaging condition, respectively. When multi-shot RTM is applied using the VEA imaging condition, each grid point in each shot migration result corresponds to a reflection coefficient and incident angle. Then, the reflection coefficient of each grid point is sorted according to the angle of incidence and angle-domain common-image gathers (ADCIGs) are I x, , z .
ACCEPTED MANUSCRIPT For the source-normalized crosscorrelation (SNC) imaging condition, the application to elastic RTM based on vector decomposition also involves vector issues. Referring to Equations
PT
11 and 12, the vector-based source-normalized crosscorrelation (VSNC) imaging conditions of
RI
the PP and PS waves can be represented as,
nt
1 cos 2
p src
p ( x, z, t ) v rec ( x, z, t )
0
nt
v
p src
nt
s ( x, z, t ) v rec ( x, z, t )
0
nt
MA
I PS x, z
p src
( x, z , t )
,
(14)
2
NU
0
v
SC
I PP x, z
v
v
p src
( x, z, t )
2
1 , cos
(15)
TE
D
0
where I pp and I ps are the imaging results of the PP and PS waves, respectively. However, one
AC CE P
difficulty in the VSNC imaging condition is that it is not easy to determine the corresponding incident angle because the imaging value of each grid point is the sum of the imaging results at all times. In this paper, in order to compare the imaging results of the VEA and VSNC imaging conditions, the incident angle in the VSNC imaging condition also accepts the incident angle of the former. To summarize, a flow chart of the new imaging method is shown in Figure 4. The source and receiver wavefields are decomposed using the decoupled propagation. The P-wave particle velocity vector, P-wave stress, and imaging time of the source wavefield are saved. In the process of reverse time propagation, PP images can be obtained at imaging time with the PP imaging condition, and PS images can be similarly obtained. At the same time, the incident angle
ACCEPTED MANUSCRIPT and angle are calculated using the corresponding Poynting vectors. After the image of each grid point is sorted according to the incident angle, the ADCIGs are obtained. To reduce the
PT
artifacts, the PP and PS images are stacked with the restriction of the incident angle. Note that
RI
because model smoothing has been discussed in many articles about RTM, we do not provide the details here. Before RTM processing, the velocity models tested in this article were properly
AC CE P
TE
D
MA
NU
SC
smoothed.
Figure 4. Flow chart of elastic RTM with the VEA imaging condition.
3. Model test 3.1 Layered model test
ACCEPTED MANUSCRIPT To test the procedure, we first designed a simple three-layer model (Figure 5). The grid spacing in two directions was 5 m, and the PML boundary condition was added with 10 grid
PT
points. The source wavelet was a Ricker wavelet with a main frequency of 30 Hz. Before applying elastic RTM, we first tested the efficiency of the wavefield decomposition method.
RI
Based on the results shown in Figure 6, the wavefield decomposition method is valid and the
TE
D
MA
NU
SC
Poynting vectors calculated using the P- and S-wave are acceptable.
AC CE P
Figure 5. Parameters of the layered model
AC CE P
TE
D
MA
NU
SC
RI
PT
ACCEPTED MANUSCRIPT
Figure 6. Pictures of the x- (a) and z-components (b) of the particle velocity; x- (c) and zcomponents (d) of the decomposed P-wave; x- (e) and z-components (f) of the decomposed Swave, respectively. Panels (g and h) show the calculated Poynting vectors of the P- and S-waves, respectively. The single-shot migration results are shown in Figure 7. The black boxes in the figures show the effective migration results, which indicates that the migration aperture of the PS wave
ACCEPTED MANUSCRIPT is much larger than that of the PP wave. Note that the PS images are poorly illuminated near the centre of the panels due to the weak energy of the converted wave when the incident angle was
PT
small. These results demonstrate that the VEA and VSNC imaging conditions yield imaging results similar to the EA and SNC conditions, and that the directly calculated reflection
RI
coefficients are representative. In addition, it can be seen that the PS wave migration result with
AC CE P
TE
D
MA
NU
SC
the VEA imaging condition has the highest spatial resolution.
ACCEPTED MANUSCRIPT Figure 7. Single-shot migration results calculated using the EA, VEA, SNC, and VSNC imaging conditions are represented in four columns from top to bottom, respectively. The left are the PP
PT
wave images and the right are the PS wave images.
RI
Figure 8 provides a comparison of the first layer reflection coefficients obtained using the
SC
VEA and VSNC imaging conditions and the Zoeppritz equation. The PP wave reflection coefficients obtained from the two imaging conditions agree with that calculated from Zoeppritz
NU
theory when the offset is between 0–0.7 km (corresponding to an incidence angle of 0–35°). The
MA
PS wave reflection coefficients are in good agreement with the Zoeppritz theory when the offset
AC CE P
TE
D
is between 0.1 km and 1.5 km (corresponding to an incidence angle of 5–56°).
Figure 8. Comparison of the first layer reflection coefficients obtained from the two imaging conditions and the Zoeppritz equation. Panel (a) shows the PP wave reflection coefficients, and panel (b) shows the PS wave reflection coefficients. 3.2 Marmousi II test
ACCEPTED MANUSCRIPT We also performed a complex model test based on the Marmousi II model (Figure 9). The grid spacing in both directions was 5 m, the time sampling interval was 0.5 ms, and the total
PT
recording time was 4 s. A total of 101 shots were distributed over the surface with 50 m spacing. A total of 1,001 receivers were located on the surface. The water layer was replaced by a rock
RI
layer with a non-zero S-wave velocity. The wavefield modelling method and boundary were the
AC CE P
TE
D
MA
NU
SC
same as previous, and the source frequency was reduced to 20 Hz.
Figure 9. A portion of the Marmousi II model. Panels a, b, and c represent the P-wave velocity, S-wave velocity, and density distributions, respectively.
ACCEPTED MANUSCRIPT The ADCIGs extracted from the PP and PS wave migration results based on the VEA imaging condition are shown in Figures 10a and 10c (only the angle gathers corresponding to the
PT
black dashed line in Figure 9a are shown). From the figures, many reflection layers can be distinguished. Compared with PP ADCIGs, the PS ADCIGs have a higher resolution and can
RI
therefore provide a more detailed depiction. As mentioned earlier, each grid point corresponds to
SC
only one reflection coefficient and incident angle pair in each shot migration result.
NU
Consequently, the number of effective reflection coefficients of ADCIG for each grid point is less than the number of shots, which results in some discrete points in the ADCIGs rather than
MA
continuous events. However, as the number of shots increases, the events in the ADCIGs will become clearer. If a velocity analysis is performed using ADCIGs, continuous events may be
D
necessary, which would require a sufficient number of shots or interpolating ADCIGs. In
TE
addition, when the velocity model has –10% errors, the events shown in Figures 10b and 10d are
AC CE P
no longer horizontal, but curve upward.
AC CE P
TE
D
MA
NU
SC
RI
PT
ACCEPTED MANUSCRIPT
Figure 10. The columns in (a) and (c) are PP and PS ADCIGs with the true velocity model, respectively. The columns in (b) and (d) are similar to those in (a) and (c), except that the velocity model has –10% errors.
Figure 11 show the stacked results of 101 shots calculated using the VEA and VSNC imaging conditions with the true velocity model. Overall, these imaging results are very clear and many details of the subsurface structure are recognizable. By comparing the detailed areas (indicated by the black arrows), the results obtained with the VEA imaging condition have a higher imaging resolution than those obtained with the VSNC imaging condition. The PS wave migration results are imaged with higher resolution and richer imaging detail than the PP wave migration results.
AC CE P
TE
D
MA
NU
SC
RI
PT
ACCEPTED MANUSCRIPT
ACCEPTED MANUSCRIPT Figure 11. Panels (a) and (b) are the stacked results of PP wave imaging with the VEA and VSNC imaging conditions, respectively. Panels (c) and (d) are the stacked results of PS wave
PT
imaging with the VEA and VSNC imaging conditions, respectively.
RI
4. Discussion
SC
The VEA imaging condition is an extension of the EA imaging condition. First, only the
NU
P-wave particle velocity vector and corresponding P-wave stress are saved at the imaging time of each grid point in the forward propagation, instead of the complete wavefield at each moment.
MA
Compared with the traditional crosscorrelation imaging conditions and wavefield reconstruction method, this approach greatly reduces the storage and computing requirements and allows the
D
elastic RTM to be applied to actual data. In addition, only the maximum amplitude is used for
TE
imaging, which is the largest contribution to the image, and results in a higher resolution. Second,
AC CE P
as was the case in the research of Wang and McMechan (2015), the wavefield decomposition method proposed in this paper is based on vector decomposition. The PS wave images do not involve polarization reversal, so multi-shot imaging results can be directly stacked without the need for polarization correction.
Although there are some common points between the imaging conditions, there are also several differences. First, the VEA imaging condition directly includes information regarding the particle vibration directions in the process of imaging, which means that the signed angledependent reflection coefficients can be directly obtained without calculating the sign. Second, a new equation to calculate the incident angle is proposed that can distinguish between the left or right incidence. Third, the VSNC imaging condition is deduced by extending the concept to the SNC imaging condition.
ACCEPTED MANUSCRIPT In addition, there are some things that need to be explained. First, the ADCIGs can be extracted directly from the multi-shot imaging results and then used in the velocity analyses. At
PT
the same time, the ADCIGs can also constrain the superposition process and allow the results to be superimposed with the effective imaging angle to reduce the imaging artifacts. Second, the
RI
process of wave propagation and imaging uses the OpenMP parallel algorithm. In actual data
SC
calculations or three-dimensional (3D) RTM calculations, the MPI parallel technique can also be
NU
used, which will greatly improve the computational efficiency.
MA
Although the methodology and tests shown above are about 2D isotropic media, it is possible to extend them to 2D anisotropic media, 3D isotropic media, and attenuation. If the 2D
D
model contains anisotropic media, taking VTI media as an example, the imaging conditions can
TE
also be applied after separating the P- and S-waves well (Xiao and Leaney, 2010; Wang et al.,
AC CE P
2015). For 3D, if the medium is isotropic, the particle vibration directions of the incident P-wave, the reflected P- and S-waves are in the same plane, except that the particle velocity vector contains three components. On this basis, it is easy to extend this method to 3D isotropic media. However, for anisotropic media, there is a change in the azimuth angle, and the reflected S-wave may be split into SV and SH waves. Thus, this method cannot be simply extended to 3D anisotropic media, which requires more studies. The maximum amplitude is used to calculate the correct imaging time. The imaging time is the time when the particle amplitude of each grid point at its largest over all times, similar to the one-way time from the source to that grid point. When the attenuative anomalies exist, although the amplitude will be attenuated, the amplitude of the wavefront is still the largest over all times. Therefore, we infer that the attenuative anomalies do not affect the calculation of the imaging time, and it will not affect the application of this method.
ACCEPTED MANUSCRIPT Finally, the multipath problem cannot be solved (Wang and McMechen, 2015). In the original velocity model, there is a low-velocity body (marked with the hollow arrow in Figure
PT
9a). The P-wave velocity of the surrounding rocks is about 1,700 m/s and the P-wave velocity of the low-velocity body is about 1,100 m/s. There will be strong interface multiples in this low-
RI
velocity body, and a strong diffraction wave at both ends. The imaging results will appear in a
SC
strong illusion under the low-velocity body (Du et al., 2015). Therefore, we set the P-wave
NU
velocity of this low-velocity body to 1,500 m/s, so that the difference between the surrounding rocks was not very large. Consequently, there is no obvious illusion in the imaging results.
MA
Another noteworthy problem is the jaggedness in the PS imaging results (Figures 11c and 11d) when the interface is tilted, which also occurs in the original velocity model due to the regular
AC CE P
5. Conclusion
TE
D
grid. This shows that the PS imaging results have a high resolution.
In this paper, we propose a new process for 2D elastic RTM. The decomposed P- and Swaves are obtained by the decoupled propagation of the source and receiver wavefields. In the process of forward modelling, the imaging time of each grid point is calculated, and the corresponding P-wave particle velocity vector and P-wave stress are saved. During the process of reverse time propagation, each grid point is imaged using the VEA imaging condition at the imaging time. The P- and S-wave Poynting vectors are calculated to obtain the incident and reflection angles, which are then used to extract the ADCIGs and amplitude corrections. The migration results of the simple model and the Marmousi II model show that the VEA imaging condition can be used to obtain signed angle-dependent reflection coefficients directly with high computational efficiency and strong adaptability. The PS wave images have high imaging
ACCEPTED MANUSCRIPT resolution and rich imaging details. In addition, by extending the concept to the SNC imaging condition, the VSNC imaging condition is deduced and the synthetic results prove its validity.
PT
Acknowledgments
RI
The research work in this paper was supported by the Natural Science Foundation of
SC
China (NSFC) [Grant No. 41374108] and Major Projects of the National Science and
AC CE P
TE
D
MA
NU
Technology of China [Grant No. 2016ZX05026-002-003].
ACCEPTED MANUSCRIPT References Baysal, E., Kosloff, D., Sherwood, J., 1983. Reverse time migration. Geophysics 48 (11),
RI
PT
1514-1524, doi: 10.1190/1.1441434.
Claerbout, J.F., 1971. Towards a unified theory of reflector mapping. Geophysics 36 (3),
SC
467-481, doi: 10.1190/1.1440185.
NU
Chang, W., McMechan, G.A., 1987. Elastic reverse-time migration. Geophysics 52 (10),
MA
1365-1375, doi: 10.1190/1.1442249.
Cai, X., Liu, Y., Ren, Z., 2015. Elastic Reverse-Time Migration Based on the Optimal
AC CE P
TE
10.1190/segam2015-5906925.1.
D
Staggered-Grid Finite-Difference Scheme. SEG, Expanded Abstracts, 4148-4152, doi:
Dong, L.G., Ma, Z.T., Cao, J.Z., Wang, H.Z., Geng, J.H., Lei, B., 2000. A staggered-grid high-order difference method of one-order elastic wave equation. Chinese Journal of Geophysics 43 (3), 411-419.
Dussaud, E., Symes, W.W., Williamson, P., Lemaistre, L., Singer, P., Denel, B., 2008. Computational strategies for reverse-time migration. SEG, Expanded Abstracts, 2267-2271, doi: 10.1190/1.3059336. Du, Q., Zhang, M., Gong, X., Chen, X., 2015. Polarity-consistent excitation amplitude imaging condition for elastic reverse time migration. Journal of Geophysics and Engineering 12 (1), 33-44, doi: 10.1088/1742-2132/12/1/33.
ACCEPTED MANUSCRIPT Gu, B.L., Liu, Y.S., Ma, X.N., Li, Z.Y., Liang, G.H., 2015. A modified excitation amplitude imaging condition for prestack reverse time migration. Exploration Geophysics 46 (4),
PT
359-370, doi: 10.1071/EG14039.
RI
Hastings, F.D., Schneider, J. B., Broschat, S. L., 1996. Application of the perfectly
SC
matched layer (PML) absorbing boundary conditions to elastic wave propagation. Journal of the
NU
Acoustical Society of America 100 (5), 3061-3069, doi: 10.1121/1.417118. Komatitsch, D., Goddeke, D., Erlebacher, G., Michea, D., 2010. Modeling the
MA
propagation of elastic waves using spectral elements on a cluster of 192 GPUs. Computer
D
Science - Research and Development 25, 75-82, doi: 10.1007/s00450-010-0109-1.
TE
Levander, A.R., 1988. Fourth-Order Finite-Difference P-SV Seismograms. Geophysics,
AC CE P
53 (11), 1425-1436, doi: 10.1190/1.1442422. Liu, Y., 2013. Globally optimal finite-difference schemes based on least squares. Geophysics 78 (4), T113-T132, doi: 10.1190/geo2012-0480.1. Liu, S., Li, X., Wang, W., Zhu, T., 2015. Source wavefield reconstruction using a linear combination of the boundary wavefield in reverse time migration. Geophysics 80 (6), S203-S212. doi: 10.1190/geo2015-0109.1. Li, Z., Ma, X., Fu, C., Liang, G., 2016. Wavefield separation and polarity reversal correction in elastic reverse time migration. Journal of Applied Geophysics 127, 56-67, doi: 10.1016/j.jappgeo.2016.02.012.
ACCEPTED MANUSCRIPT Nguyen, B.D., McMechan, G. A., 2013. Excitation amplitude imaging condition for prestack reverse-time migration. Geophysics 78 (1), S37-S46, doi: 10.1190/geo2012-0079.1.
PT
Nguyen, B.D., McMechan, G. A., 2015. Five ways to avoid storing source wavefield
RI
snapshots in 2d elastic prestack reverse time migration, Geophysics 80 (1), S1-S18, doi:
SC
10.1190/geo2014-0014.1.
NU
Sun, R., McMechan, G. A., 2001. Scalar reverse-time depth migration of prestack elastic
MA
seismic data. Geophysics 66 (5), 1519-1527, doi: 10.1190/1.1487098. Sun, W., Fu, L.Y., 2013. Two effective approaches to reduce data storage in reverse time
D
migration. Computers and Geosciences 56, 69-75, doi: 10.1016/j.cageo.2013.03.013.
TE
Tang, C., Wang, D., 2012. Reverse Time Migration with Source Wavefield
AC CE P
Reconstruction in Target Imaging Region. EAGE, Extended Abstracts, P272, doi: 10.3997/22144609.20148620.
Virieux, J., 1986. P-SV wave propagation in heterogeneous media: velocity-stress finitedifference method, Geophysics 51 (4), 889–901. doi: 10.1190/1.1442147. Wapenaar, C.P., Haime, G.C., 1990. Elastic extrapolation of primary seismic p- and swaves, Geophysical Prospecting 38 (1), 23-60. doi: 10.1111/j.1365-2478.1990.tb01833.x. Wang, W., McMechan, G.A., Zhang, Q., 2015. Comparison of two algorithms for isotropic elastic P and S vector decomposition. Geophysics 80 (4), T147-T160, doi: 10.1190/geo2014-0563.1.
ACCEPTED MANUSCRIPT Wang, W., McMechan, G.A., 2015. Vector-based elastic reverse time migration. Geophysics 80 (6), S245-S258, doi: 10.1190/geo2014-0620.1.
PT
Xiao, X., Leaney, W.S., 2010. Local vertical seismic profiling (VSP) elastic reverse-time
SC
Geophysics 75 (2), S35-S49, doi: 10.1190/1.3309460.
RI
migration and migration resolution: salt-flank imaging with transmitted P-to-s waves.
NU
Xin, W., Yan, Z.C., Liang, W.Q., Chen, Y.H., Yang, C.C., 2015. Methods to determine the finite difference coefficients for elastic wave equation modeling. Chinese Journal of
MA
Geophysics- Chinese Edition 58 (7), 2486-2495, doi: 10.6038/cjg20150724.
D
Yan, J., Sava, P., 2008. Isotropic angle-domain elastic reverse-time migration.
TE
Geophysics 73 (6), S229-S239, doi: 10.1190/1.2981241.
AC CE P
Yan, R., Xie, X., 2012. An angle-domain imaging condition for elastic reverse time migration and its application to angle gather extraction. Geophysics 77 (5), S105-S115, doi: 10.1190/geo2011-0455.1.
Yang, L., Yan, H., Liu, H., 2014. Least squares staggered-grid finite-difference for elastic wave modeling. Exploration Geophysics 45 (4), 255-260, doi: 10.1071/EG13087. Zhang, J., Tian, Z., Wang, C., 2007. P- and S-wave-separated elastic wave-equation numerical modeling using 2D staggered grid. SEG, Expanded Abstracts, 2104-2109, doi: 10.1190/1.2792904.
ACCEPTED MANUSCRIPT
AC CE P
TE
D
MA
NU
SC
RI
PT
Highlights: - A new method for 2D elastic multicomponent reverse time migration (RTM) is proposed. - The vector-based excitation amplitude imaging condition is deduced. - The signed angle-dependent reflection coefficients are directly obtained. - The concept is extended to the source-normalized crosscorrelation imaging condition.