Gaussian beam migration for sparse common-shot data - CiteSeerX
Gaussian beam migration for sparse common-shot and common-receiver data Robert L. Nowack*, Purdue University, USA, Mrinal K. Sen and Paul L. Stoffa, The University of Texas at Austin, USA. Summary We investigate the Gaussian beam migration of commonshot and common-receiver data. The imaging of commonshot data is useful for seismic data where the receiver coordinates are well sampled, but the source coordinates are less well sampled. By reciprocity, this approach can also be applied to common-receiver data, such as from OBS experiments where the source locations are dense but the receiver locations are sparse. Since Gaussian beam migration uses smoothed, localized windowing of the data, it can provide more stable results for the inversion of sparsely sampled data. In order to test the Gaussian beam imaging, pre-stack data sets have been computed using the finite difference method for slices obtained from the SEG/EAGE salt model which has significant lateral velocity variations. Prior to tracing of the beams, smoothing of the velocity was performed. The phase times were then corrected using first-order perturbation theory. A receiver spacing of 20 m was used with an offset range of 5000 m, and the shot spacing was allowed to vary. Comparing the Gaussian beam migration images with the true model, most features were well imaged. In particular, features beneath the salt were imaged including steeply dipping events and a lower flat reflector. The Gaussian beam images also compare well with results from other migration methods. Introduction The inversion of common-shot data is important for seismic data where the receiver coordinates are well sampled, but the source coordinates are less well sampled. Since Gaussian beam migration uses smoothed localized windowing of the data, it can provide stable results for the inversion of sparsely sampled data in one coordinate. By reciprocity, this can also be applied to common-receiver data, such as from OBS or onshore/offshore experiments where the source coordinates are dense but the receiver locations are sparse. Theory The Gaussian beam method is an asymptotic method for the computation of wavefields in smoothly varying inhomogeneous media (Popov, 1982; Červeny, et al. 1982; Nowack and Aki, 1984; Babich and Popov, 1989; Nowack, 2003) that describes high-frequency seismic wavefields by the summation of paraxial Gaussian beams. The
advantages of Gaussian beam synthesis over other asymptotic methods are: • Individual Gaussian beams have no singularities either at caustics in the spatial domain or at pseudo-caustics in the wavenumber domain. • The Gaussian beam method naturally introduces smoothing effects and is therefore not as sensitive to model parameterizations as the ray method. • The Gaussian beam method does not require two-point ray tracing and naturally allows for the incorporation of multiple arrivals in either forward or inverse modeling. In the 2-D case, the Green’s function can be written in terms of a Gaussian beam summation as
g ( x, x ' , ω ) =
−i 4π
⎛ p0 ⎞ ⎜ ⎟ ⎜v ⎟ ⎝ 0 ⎠
1/ 2
∫ dγ u
gb
γ
( x, x', ω )
where
⎡ v( s) ⎤ u γgb ( x, x', ω ) = ⎢ ⎥ ⎣ q( s) ⎦
1/ 2
iω ⎧ ⎫ M ( s ) n2 ⎬ exp⎨i ω τ ( s ) + 2 ⎩ ⎭
and (s,n) corresponds to the point x in ray centered coordinates, τ(s) is the travel time along the central ray, v(s) is the velocity along the center ray, p0 = p(s=0), and M ( s ) = p( s ) / q ( s ) . The take-off angle γ at the source can be written as dγ
dp1′ . p 3′
=
The migration of the individual shot gathers can be written as, dω dp1g L g L s g I s ( x) = ∑ ∞ ∫ A(ω ) U ( x, x , p , ω ) D( x , x , p , ω ) 2 π ∫ p3g m = −∞ where A ( ω ) is an amplitude factor, U ( x, x , p , ω ) is a L
g
propagation integral and D ( x , x , p , ω ) is the local L
s
g
slant-stack of the data that has been Gaussian windowed about a center distance and dip angle. This is similar to Hill (2001) where we have used a common-shot instead of a common-offset geometry. To summarize the computational aspects, first the complex travel-times and amplitudes at subsurface points x from the Gaussian beams originating from (xL, pg) for each data window centered on xL are computed. This can be done on
Gaussian beam migration for sparse data
∆L ~
⎛ω 2 ⎜⎜ r ⎝ ωH
1/ 2
⎞ ⎟⎟ W0 ⎠
where ωr is the lowest frequency in the data and ωH is the highest (Hale, 1992). W0 is the width of the Gaussian beam at the reference frequency ω r . p1g must also be adequately sampled with
∆p1g
≅
1 1 , W0 (ω r ω H )1 / 2
(Hale, 1992 and Hill, 1990). Using the values of W0 ,
∆ L , and ∆ p1g above, the image for each shot point can be numerically evaluated. The final image is then the sum of the individual common-shot images. Alternative forms of asymptotic synthesis and inversion of seismic data include wavepath migration (Sun and Schuster, 2001), Coherent states (Thomson, 2001; Albertin, 2001) and the Maslov method (Chapman and Drummond, 1982; Xu et al., 1998). Also, other early applications of Gaussian beam migration include Costa et al. (1989) and Lazaratos and Harris (1990) and Alkhalifah (1995). Examples In order to test the Gaussian beam algorithm, several preliminary synthetic applications are performed. To compare with the Gaussian beam inversion results, a ray/Born inversion based on the formulation of Bleistein et al. (1987) is also used (see also, Bleistein et al., 2001). Figure 1: The top plot shows a slice from the SEG/EAGE salt model. The bottom plot shows the smoothed velocity model used in the beam tracing for the Gaussian beam migration.
a course grid and interpolated to a finer grid. Also, only subsurface points within the support of the beams need to be retained and this provides added computational savings, particularly in 3-D. As an approximation, the propagation components from the source are computed using just the first arrivals. Beam tracing is used to compute the paraxial Gaussian beams for the receiver propagation components. The triplicated arrivals are then accounted for in the propagation from the receivers, but only the first arrivals for the propagation from the source. To retain all the multipathing, the complete integral for U would need to be evaluated by Gaussian beams, but this would require more computations and is not considered here. The spacing of the beams in x and p depends on the frequency of the data. In 2-D, one choice is
A realistic synthetic example is taken from a prestack data set computed using the finite difference method for a 2-D slice of the SEG/EAGE salt model (Aminzadeh, et al., 1996; O’Brien and Gray, 1996). The model has a 20 m spacing in x and z and has 716 points in x and 250 in z which includes a 20 point padding on the sides. The receivers have a spacing of 20 m along the surface. Figure 1 shows the slice through the model with a salt body in the center as well as different reflectors. The model also has significant lateral velocity variations. For the generation of the synthetic data, no multiples from the free surface were included. The top plot of Figure 1 shows the original 2-D structure taken from the SEG/EAGE salt model. In order to trace rays and beams through the model, a smoothed velocity model is needed, and an example is shown at the bottom of Figure 1. If no prior information is known, a highly smoothed velocity model is justified and also gives a stable propagation of rays and beams through the model. For the simple case here, first-order perturbation theory is applied to obtain perturbed phase
Gaussian beam migration for sparse data times between the smoothed velocity model and the true model and these are used to correct the migration image. In this example, a sparse sampling of shots is used with a spacing alternating between 240 m and 260 m. The receiver spacing is 20 m and the receiver offsets range from plus and minus 2500 m except near the ends of the profile where the total range of receiver coverage is taken to be 5000 m. Figure 2 shows the results of inversions with a sparse shot array. The final inversion images are generally well migrated even with the sparse shot spacing. Comparing the Gaussian beam image in the top plot with the true model in Figure 1, several features can be seen beneath the salt body including some steeply dipping events and the lower flat reflector, but are not as well imaged. As shown by O’Brien and Gray (1996), subsalt imaging can be improved by the use of larger offset. Similar to the Gaussian beam inversion, the common-shot ray/Born images are stacked to obtain a final image. The results of the inversion of the sparse shot gathers with the ray/Born formulation are shown in the bottom plot of Figure 2. Compared to the Gaussian beam inversion and the true model in Figure 1, the steeply dipping faults above the salt body are not as well imaged. Also, the structures beneath the salt are much less clear. This is due to the restricted offset range used, as well as the use of only first arrivals for the propagation from both sources and receivers in the Born inversion. In addition, some artifacts for both inversions may result from inter-bed multiples in the synthetic data. Thus for the data set used, the inversion from the Gaussian beam method gives a much clearer image than that from the common-shot ray/Born inversion. An even better image of the sub-salt structure could be obtained from the Gaussian beam migration if a beam decomposition were used for both the source and receiver propagation terms, but this would require more extensive computations. Nonetheless, the Gaussian beam results are still better than those obtained using the ray/Born inversion using only first arrivals. This also results from the smoothing effects of the beams for the sparse source coverage. Still, improved subsalt images should be obtained from both methods if greater offsets were used. Conclusions We have investigated common-shot Gaussian beam inversion which by reciprocity can also be applied to common-receiver inversions, such as from OBS or onshore/offshore experiments. The Gaussian beam method
Figure 2: These plots show the inversion results for receivers every 20 m over a 5000 m aperture centered about each shot. The source array is sparse with a spacing alternating between 240 and 260 m. The top plot shows the results from the Gaussian beam inversion and the bottom plot shows the ray/Born inversion.
is tested using the SEG/EAGE salt model and compared with the results from a ray/Born inversion. The results for the Gaussian beam inversion compare well with the true model and give better images for structures directly above the salt body as well as below the salt body compared with the Ray/Born inversion. References Albertin, U., Yingst, D. and Jaramillo, H., 2001, Comparing common-offset Maslov, Gaussian Beam, and Coherent state migration, 71th Ann. Internat. Mtg. Soc. Expl. Geophys., Expanded Abstracts, 913-916. Alkhalifah, T., 1995, Gaussian beam depth migration for anisotropic media, Geophysics, 60, 1474-1484.
Gaussian beam migration for sparse data
Aminzadeh, F., Brac, J. and Kunz, T., 1996, Threedimensional SEG/EAGE salt model – an update, The Leading Edge, 15, no. 2, 131-134.
Popov, M.M., 1982, A new method of computation of wave fields using Gaussian beams, Wave Motion, 4, 85-97. Sun, H. and Schuster, G., 2001, 2-D wavepath migration, Geophysics, 66, 1528-1537.
Babich, V.M and Popov, M.M., 1989, Gaussian summation method (review), translated in Radiophysics and Quantum Electronics 32, 1063-1081.
Thomson, C.J., 2001, Seismic coherent states and ray geometrical spreading, Geophys. J. Int., 144, 320-342.
Bleistein, N., Cohen, J.K., Hagin, F.G., 1987, Two and one half dimensional inversion with an arbitrary reference, Geophysics, 52, 26-36.
Xu, S. and Lambere, G., 1998, Maslov+Born migration /inversion in complex media, 68th Ann. Internat. Mtg. Soc. Expl. Geophys., Expanded Abstracts, 1704-1707.
Bleistein, N., Cohen, J.K. and Stockwell, Jr., J.W., 2001, Mathematics of Multichannel Seismic Imaging, Migration and Inversion, Springer, New York. Červeny, V., Popov, M.M. and Pšencík, I., 1982, Computation of wavefields in inhomogeneous media – Gaussian beam approach, Geophys. J.R. Astr. Soc., 70, 109-128. Chapman, C.H. and Drummond, R., 1982, Body-wave seismograms in inhomogeneous media using Maslov asymptotic theory, Bull. Seism. Soc. Am., 72, S277-S317. Costa, C.A., Raz, S. and Kosloff, D., 1989, Gaussian Beam migration, 59th Ann. Internat. Mtg. Soc. Expl. Geophys., Expanded Abstracts, 1169-1171. Hale, D., 1992, Migration by the Kirchhoff, slant stack and Gaussian beam methods, CWD-121, Center for Wave Phenomena, Colorado School of Mines, Golden, CO. Hill, N. R., 1990, Gaussian beam migration, Geophysics, 55, 1416-1428. Hill, N. R., 2001, Prestack Gaussian-beam depth migration, Geophysics, 66, 1240-1250. Lazaratos, S.K. and. Harris, J.M., 1990, Radon Transform /Gaussian Beam migration, 60st Ann. Internat. Mtg. Soc. Expl. Geophys., Expanded Abstracts, 1314-1317. Nowack, R.L., 2003, Calculation of synthetic seismograms with Gaussian Beams, Pure and Applied Geophys., 160, 487507. Nowack, R.L. and Aki, K., 1984, The two-dimensional Gaussian beam synthetic method: testing and application, J. Geophys., 89, 7797-7819. O’Brien, M. J. and Gray, A.H., 1996, Can we image beneath salt? The Leading Edge, 15, no. 1, 17-22.
advantages of Gaussian beam synthesis over other asymptotic methods are: • Individual Gaussian beams have no singularities either at caustics in the spatial domain or at pseudo-caustics in the wavenumber domain. • The Gaussian beam method naturally introduces smoothing effects and is therefore not as sensitive to model parameterizations as the ray method. • The Gaussian beam method does not require two-point ray tracing and naturally allows for the incorporation of multiple arrivals in either forward or inverse modeling. In the 2-D case, the Green’s function can be written in terms of a Gaussian beam summation as
g ( x, x ' , ω ) =
−i 4π
⎛ p0 ⎞ ⎜ ⎟ ⎜v ⎟ ⎝ 0 ⎠
1/ 2
∫ dγ u
gb
γ
( x, x', ω )
where
⎡ v( s) ⎤ u γgb ( x, x', ω ) = ⎢ ⎥ ⎣ q( s) ⎦
1/ 2
iω ⎧ ⎫ M ( s ) n2 ⎬ exp⎨i ω τ ( s ) + 2 ⎩ ⎭
and (s,n) corresponds to the point x in ray centered coordinates, τ(s) is the travel time along the central ray, v(s) is the velocity along the center ray, p0 = p(s=0), and M ( s ) = p( s ) / q ( s ) . The take-off angle γ at the source can be written as dγ
dp1′ . p 3′
=
The migration of the individual shot gathers can be written as, dω dp1g L g L s g I s ( x) = ∑ ∞ ∫ A(ω ) U ( x, x , p , ω ) D( x , x , p , ω ) 2 π ∫ p3g m = −∞ where A ( ω ) is an amplitude factor, U ( x, x , p , ω ) is a L
g
propagation integral and D ( x , x , p , ω ) is the local L
s
g
slant-stack of the data that has been Gaussian windowed about a center distance and dip angle. This is similar to Hill (2001) where we have used a common-shot instead of a common-offset geometry. To summarize the computational aspects, first the complex travel-times and amplitudes at subsurface points x from the Gaussian beams originating from (xL, pg) for each data window centered on xL are computed. This can be done on
Gaussian beam migration for sparse data
∆L ~
⎛ω 2 ⎜⎜ r ⎝ ωH
1/ 2
⎞ ⎟⎟ W0 ⎠
where ωr is the lowest frequency in the data and ωH is the highest (Hale, 1992). W0 is the width of the Gaussian beam at the reference frequency ω r . p1g must also be adequately sampled with
∆p1g
≅
1 1 , W0 (ω r ω H )1 / 2
(Hale, 1992 and Hill, 1990). Using the values of W0 ,
∆ L , and ∆ p1g above, the image for each shot point can be numerically evaluated. The final image is then the sum of the individual common-shot images. Alternative forms of asymptotic synthesis and inversion of seismic data include wavepath migration (Sun and Schuster, 2001), Coherent states (Thomson, 2001; Albertin, 2001) and the Maslov method (Chapman and Drummond, 1982; Xu et al., 1998). Also, other early applications of Gaussian beam migration include Costa et al. (1989) and Lazaratos and Harris (1990) and Alkhalifah (1995). Examples In order to test the Gaussian beam algorithm, several preliminary synthetic applications are performed. To compare with the Gaussian beam inversion results, a ray/Born inversion based on the formulation of Bleistein et al. (1987) is also used (see also, Bleistein et al., 2001). Figure 1: The top plot shows a slice from the SEG/EAGE salt model. The bottom plot shows the smoothed velocity model used in the beam tracing for the Gaussian beam migration.
a course grid and interpolated to a finer grid. Also, only subsurface points within the support of the beams need to be retained and this provides added computational savings, particularly in 3-D. As an approximation, the propagation components from the source are computed using just the first arrivals. Beam tracing is used to compute the paraxial Gaussian beams for the receiver propagation components. The triplicated arrivals are then accounted for in the propagation from the receivers, but only the first arrivals for the propagation from the source. To retain all the multipathing, the complete integral for U would need to be evaluated by Gaussian beams, but this would require more computations and is not considered here. The spacing of the beams in x and p depends on the frequency of the data. In 2-D, one choice is
A realistic synthetic example is taken from a prestack data set computed using the finite difference method for a 2-D slice of the SEG/EAGE salt model (Aminzadeh, et al., 1996; O’Brien and Gray, 1996). The model has a 20 m spacing in x and z and has 716 points in x and 250 in z which includes a 20 point padding on the sides. The receivers have a spacing of 20 m along the surface. Figure 1 shows the slice through the model with a salt body in the center as well as different reflectors. The model also has significant lateral velocity variations. For the generation of the synthetic data, no multiples from the free surface were included. The top plot of Figure 1 shows the original 2-D structure taken from the SEG/EAGE salt model. In order to trace rays and beams through the model, a smoothed velocity model is needed, and an example is shown at the bottom of Figure 1. If no prior information is known, a highly smoothed velocity model is justified and also gives a stable propagation of rays and beams through the model. For the simple case here, first-order perturbation theory is applied to obtain perturbed phase
Gaussian beam migration for sparse data times between the smoothed velocity model and the true model and these are used to correct the migration image. In this example, a sparse sampling of shots is used with a spacing alternating between 240 m and 260 m. The receiver spacing is 20 m and the receiver offsets range from plus and minus 2500 m except near the ends of the profile where the total range of receiver coverage is taken to be 5000 m. Figure 2 shows the results of inversions with a sparse shot array. The final inversion images are generally well migrated even with the sparse shot spacing. Comparing the Gaussian beam image in the top plot with the true model in Figure 1, several features can be seen beneath the salt body including some steeply dipping events and the lower flat reflector, but are not as well imaged. As shown by O’Brien and Gray (1996), subsalt imaging can be improved by the use of larger offset. Similar to the Gaussian beam inversion, the common-shot ray/Born images are stacked to obtain a final image. The results of the inversion of the sparse shot gathers with the ray/Born formulation are shown in the bottom plot of Figure 2. Compared to the Gaussian beam inversion and the true model in Figure 1, the steeply dipping faults above the salt body are not as well imaged. Also, the structures beneath the salt are much less clear. This is due to the restricted offset range used, as well as the use of only first arrivals for the propagation from both sources and receivers in the Born inversion. In addition, some artifacts for both inversions may result from inter-bed multiples in the synthetic data. Thus for the data set used, the inversion from the Gaussian beam method gives a much clearer image than that from the common-shot ray/Born inversion. An even better image of the sub-salt structure could be obtained from the Gaussian beam migration if a beam decomposition were used for both the source and receiver propagation terms, but this would require more extensive computations. Nonetheless, the Gaussian beam results are still better than those obtained using the ray/Born inversion using only first arrivals. This also results from the smoothing effects of the beams for the sparse source coverage. Still, improved subsalt images should be obtained from both methods if greater offsets were used. Conclusions We have investigated common-shot Gaussian beam inversion which by reciprocity can also be applied to common-receiver inversions, such as from OBS or onshore/offshore experiments. The Gaussian beam method
Figure 2: These plots show the inversion results for receivers every 20 m over a 5000 m aperture centered about each shot. The source array is sparse with a spacing alternating between 240 and 260 m. The top plot shows the results from the Gaussian beam inversion and the bottom plot shows the ray/Born inversion.
is tested using the SEG/EAGE salt model and compared with the results from a ray/Born inversion. The results for the Gaussian beam inversion compare well with the true model and give better images for structures directly above the salt body as well as below the salt body compared with the Ray/Born inversion. References Albertin, U., Yingst, D. and Jaramillo, H., 2001, Comparing common-offset Maslov, Gaussian Beam, and Coherent state migration, 71th Ann. Internat. Mtg. Soc. Expl. Geophys., Expanded Abstracts, 913-916. Alkhalifah, T., 1995, Gaussian beam depth migration for anisotropic media, Geophysics, 60, 1474-1484.
Gaussian beam migration for sparse data
Aminzadeh, F., Brac, J. and Kunz, T., 1996, Threedimensional SEG/EAGE salt model – an update, The Leading Edge, 15, no. 2, 131-134.
Popov, M.M., 1982, A new method of computation of wave fields using Gaussian beams, Wave Motion, 4, 85-97. Sun, H. and Schuster, G., 2001, 2-D wavepath migration, Geophysics, 66, 1528-1537.
Babich, V.M and Popov, M.M., 1989, Gaussian summation method (review), translated in Radiophysics and Quantum Electronics 32, 1063-1081.
Thomson, C.J., 2001, Seismic coherent states and ray geometrical spreading, Geophys. J. Int., 144, 320-342.
Bleistein, N., Cohen, J.K., Hagin, F.G., 1987, Two and one half dimensional inversion with an arbitrary reference, Geophysics, 52, 26-36.
Xu, S. and Lambere, G., 1998, Maslov+Born migration /inversion in complex media, 68th Ann. Internat. Mtg. Soc. Expl. Geophys., Expanded Abstracts, 1704-1707.
Bleistein, N., Cohen, J.K. and Stockwell, Jr., J.W., 2001, Mathematics of Multichannel Seismic Imaging, Migration and Inversion, Springer, New York. Červeny, V., Popov, M.M. and Pšencík, I., 1982, Computation of wavefields in inhomogeneous media – Gaussian beam approach, Geophys. J.R. Astr. Soc., 70, 109-128. Chapman, C.H. and Drummond, R., 1982, Body-wave seismograms in inhomogeneous media using Maslov asymptotic theory, Bull. Seism. Soc. Am., 72, S277-S317. Costa, C.A., Raz, S. and Kosloff, D., 1989, Gaussian Beam migration, 59th Ann. Internat. Mtg. Soc. Expl. Geophys., Expanded Abstracts, 1169-1171. Hale, D., 1992, Migration by the Kirchhoff, slant stack and Gaussian beam methods, CWD-121, Center for Wave Phenomena, Colorado School of Mines, Golden, CO. Hill, N. R., 1990, Gaussian beam migration, Geophysics, 55, 1416-1428. Hill, N. R., 2001, Prestack Gaussian-beam depth migration, Geophysics, 66, 1240-1250. Lazaratos, S.K. and. Harris, J.M., 1990, Radon Transform /Gaussian Beam migration, 60st Ann. Internat. Mtg. Soc. Expl. Geophys., Expanded Abstracts, 1314-1317. Nowack, R.L., 2003, Calculation of synthetic seismograms with Gaussian Beams, Pure and Applied Geophys., 160, 487507. Nowack, R.L. and Aki, K., 1984, The two-dimensional Gaussian beam synthetic method: testing and application, J. Geophys., 89, 7797-7819. O’Brien, M. J. and Gray, A.H., 1996, Can we image beneath salt? The Leading Edge, 15, no. 1, 17-22.
