- Email: [email protected]
Resolution equivalence of dispersion-imaging methods for noise-free high-frequency surface-wave data
Journal of Applied Geophysics 122 (2015) 167–171
Contents lists available at ScienceDirect
Journal of Applied Geophysics journal homepage: www.elsevier.com/locate/jappgeo
Resolution equivalence of dispersion-imaging methods for noise-free high-frequency surface-wave data Chao Shen a,b,⁎, Ao Wang a, Limin Wang a,b,c, Zongbo Xu a,b, Feng Cheng a,b a b c
Institute of Geophysics and Geomatics, China University of Geosciences (Wuhan), 388 Lumo Rd., Wuhan, Hubei 430074, China Subsurface Imaging and Sensing Laboratory, Institute of Geophysics and Geomatics, China University of Geosciences, 388 Lumo Rd., Wuhan, Hubei 430074, China Hubei Subsurface Multi-scale Imaging Key Laboratory, China University of Geosciences (Wuhan), 388 Lumo Rd., Wuhan, Hubei 430074, China
a r t i c l e
i n f o
Article history: Received 19 July 2015 Received in revised form 2 September 2015 Accepted 22 September 2015 Available online 25 September 2015 Keywords: Surface wave Dispersion image Resolution Equivalence
a b s t r a c t A key step in high-frequency surface-wave methods is to acquire the dispersion properties of surface waves. To pick dispersion curves, surface-wave data is required to be transformed from the time–space (t–x) domain to the frequency–phase velocity (f–v) domain. In constructing accurate multi-mode dispersion curves, it is necessary to generate a reliable and high-resolution image of dispersion energy in the f–v domain. At present, there are five methods available for imaging dispersion energy: the τ-p transformation, the f-k transformation, the phase shift, the frequency decomposition and slant stacking, and the high-resolution linear Radon transformation. Among them, resolution of the high-resolution linear Radon transformation is the highest. We proposed a processing step to enhance resolution of the other four methods. Resolution of the four enhanced methods and high-resolution linear Radon transformation is compared by imaging dispersion energy of the same synthetic data. As a result, the four enhanced methods generated dispersion images with the same resolution, which revealed that the resolution of the five dispersion-imaging methods is essentially equivalent for noise-free data. © 2015 Elsevier B.V. All rights reserved.
1. Introduction Surface waves are dispersive for layered earth models. Particular modes of surface waves will possess unique phase velocities at different frequencies. A key step in high-frequency surface-wave methods is to acquire the dispersion properties of surface waves. Generating a reliable and high-resolution image of dispersion energy in frequency–phase velocity (f–v) domain is of essential necessity (Xia, 2014). From a highresolution dispersion image, accurate dispersion curves can be picked by following peaks of dispersion energy at different frequencies, leading to an effective inversion. In high-frequency surface-wave methods, there are five methods available for imaging dispersion energy at present: the τ-p transformation (McMechan and Yedlin, 1981), the f-k transformation (e.g., Yilmaz, 1987), the phase shift (Park et al., 1998), the frequency decomposition and slant stacking (Xia et al., 2007), and the high-resolution linear Radon transformation (Luo et al., 2008). Among them, high-resolution linear Radon transformation usually generates higher-resolution dispersion images than the other four methods do (Luo et al., 2008). Since the method has been utilized in many applications (Herrmann et al., 2000; Foster and Mosher, 1992; Thorson and Claerbout, 1985), it has received more and more attentions. Based on our current works, ⁎ Corresponding author at: Institute of Geophysics and Geomatics, China University of Geosciences (Wuhan), 388 Lumo Rd., Wuhan, Hubei 430074, China. E-mail address: [email protected] (C. Shen).
http://dx.doi.org/10.1016/j.jappgeo.2015.09.019 0926-9851/© 2015 Elsevier B.V. All rights reserved.
we find that there are ways to enhance resolution of the four other methods, so that the five methods can possess the same resolution power in the f–v domain. One simple practice is the numeric processing that based on power operation. In this paper, we first introduce principles of the five dispersionimaging methods and are solution-enhancing processing. Then dispersion energy of the same synthetic surface-wave shot gather is imaged using all the five methods. Besides, the processing step is executed to enhance resolution of dispersion images that are generated by τ-p transformation, f-k transformation, phase shift, and frequency decomposition and slant stacking. Finally, we conduct a comparison to summarize resolution of the five dispersion-imaging methods.
2. Methods To acquire dispersion energy in the f–v domain, a shot gather in the time–space (t–x) domain is necessary to go through two steps: one is to transform from the time domain to the frequency domain; the other one is to transform from the space domain to the phase velocity domain by velocity scanning and slant stacking. These two steps are independent and they can switch their transformation sequences with each other during calculation. Let s(x, t) be surface-wave wavefield, where x is offset (distance between a source and a receiver) and t is travel time, all receivers are evenly arranged along a survey line with a spacing Δx.
168
C. Shen et al. / Journal of Applied Geophysics 122 (2015) 167–171
2.1. τ-p transformation (McMechan and Yedlin, 1981) A slant stack operator is first applied to linearly transform wavefield s(x, t) into τ-p domain data: U ðp; τÞ ¼
N X
sðxi ; t ¼ τ þ pxi Þ;
ð1Þ
i¼1
where p is the slowness and τ is the intercept on the time axis. If the scanning slowness p is close to the real slowness of surface wave, U(p, τ) will reach a maximum. A Fourier transformation is then applied to the data in the τ domain to transform them into the frequency domain. As p is simply the reciprocal of phase velocity (p = 1/v), the f–v domain result can be easily converted and achieved.
Z
∞ −∞
∞ −∞
sðx; t Þe−2πiðftþkxÞ dtdx:
ð2Þ
This procedure could be replaced with a one-dimensional Fourier transformation in the time domain and a one-dimensional inverse Fourier transformation in the space domain. Then the f–v domain result can be transformed according to the relation: v = f/k. 2.3. Phase shift (Park et al., 1998) Phase information P( x,ω) of a wavefield can be calculated by a Fourier transformation in the time domain. It is expressed as ω
P ð x; ωÞ ¼ e−i v xþφ0 ;
ð3Þ
where ω is the circular frequency, v is the phase velocity and φ0 is the initial phase. The phase is linearly correlated to offset. Integrate the phase using a scanning φ: Z U ðφ; ωÞ ¼
Z eiφx P ðx; ωÞdx ¼
eiðφ−
ω ν
Þxþφ0
dx ¼
N X
eiðφ−
ω ν
Þxi þφ0
:
When the scanning φ gets close to ω/v, phase angle of the complex ω numbereiðφ− v Þxþφ0 is the same in the complex plane, making the integral reaching a maximum value. At last, the f–v domain result can be achieved according to the following relations: f ¼ 2πω; and ω 2πf 2π ¼ ¼ f: φ φ φ
ð5Þ ð6Þ
Phase shift method utilizes the normalized spectrum and can perform well in generating dispersion image (Moro et al., 2003). Nevertheless, due to the uniqueness of calculated phase at each frequency, dispersion energy of higher modes may not be completely imaged. 2.4. Frequency decomposition and slant stacking (Xia et al., 2007) Shot gather s(x, t) is first stretched into pseudo-vibroseis data d(x,t) using a frequency decomposition technique (Coruh, 1985): dðx; t Þ ¼ yðt Þ⊕sðx; t Þ;
n X
mðpi ; f Þei2πfpi x ;
ð8Þ
where d(x,f) is a shot gather in the frequency domain and m(p, f) is the Radon matrix. So Eq. (8) can be written as a matrix form: d ¼ Lm;
ð9Þ
where L = ei2πfpx is the linear Radon transformation (LRT) operator. Eq. (9) can be rearranged as madj ¼ LT d;
ð10Þ
where madj represents a low-resolution solution using the transpose of L. An appropriate weighting and damping factor will lead to a sparse inversion solution and improve the resolution. Thus, the objective function becomes 2 J ¼ ‖W d d−LW −1 m W m m ‖ þ λ‖W m m‖ ;
ð11Þ
and the problem turns into solving the following equation (Ji, 2006) T T −1 −T T T W −T m L W d W d LW m þ λI W m m ¼ W m L W d W d d;
ð12Þ
ð4Þ
i¼1
v¼
Standard linear Radon transformation can be expressed as the following equation
i¼1
Wavefield s(x, t) is first transformed to the f-k domain data S(k, f) by a two-dimensional Fourier transformation: Z
2.5. High-resolution linear Radon transformation (Luo et al., 2008)
dðx; f Þ ¼
2.2. F-k transformation (e.g., Yilmaz, 1987)
Sðk; f Þ ¼
So d(x, t) is also a frequency-swept data d(x, f) and each t corresponds with a frequency. A slant stacking (Yilmaz, 1987, p. 430) is then performed to obtain f–p domain data U(p, f). The f–v domain result can be eventually achieved by inverting p. The calculated result of this method is better than that of τ-p transformation because frequencies are decomposed before slant stacking. Furthermore, compared with phase shift, this algorithm could extract stronger energy for higher-mode surface waves.
ð7Þ
where ⊕ stands for the convolution operator and y(t) is a sweep function. Generally, y(t) is usually chosen as a linear function y(t) = a + bt.
where I is the identity matrix, Wd is a data weighting matrix and Wm is a model weighting matrix, λ is a damping factor and controls the compromise between the error misfit and model length. Eq. (12) can be solved by the weighted preconditioned conjugate gradient (CG) algorithm (e.g., Sacchi and Ulrych, 1995). The resolution of high-resolution LRT depends on the iteration in the CG algorithm. At the beginning of iteration, this method equals with a standard LRT (only iterate once). With the number of iteration growing, the resolution will be greatly improved (normally iterate 3 times). 2.6. A resolution-enhancing processing We define the resolution of a dispersion image as the width between the two half-values of dispersion energy at a given frequency in the f–v domain (Xia et al., 2006). τ-p transformation, f-k transformation, phase shift, frequency decomposition and slant stacking are actually kinds of a standard discretized LRT and they generate dispersion images with a low resolution (Eq. (10)). To compare them with high-resolution LRT, we deem them regular dispersionimaging methods in this paper for convenience. They suffer from incomplete information (Trad et al., 2003) and could hardly generate high-resolution dispersion images. The mathematical power operation is able to amplify the difference between large and small values, outstanding the maxima. Computing a
C. Shen et al. / Journal of Applied Geophysics 122 (2015) 167–171
169
power of the originally generated dispersion image will enhance its resolution. The explicit formula of our proposed power operation is Aenhanced ¼ Aαorignal ;
ð13Þ
where Aoriginal is a matrix of the original dispersion image generated by the four regular methods, Aenhanced is a matrix of the enhanced dispersion image with a higher resolution, and α is the exponent of the power. Resolution of the enhanced dispersion image can be controlled through adjusting α in this step. The stronger the power operation is applied to a dispersion image, the higher-resolution of dispersion energy will be imaged. High-resolution LRT images multi-mode dispersion energy simultaneously. Since the power operation will only reserve the most dominant peaks (the maximum value), an automatic normalization of dispersion energy on each frequency is preferred if higher-mode dispersion curves (normally not exhibiting the strongest energy, and being suppressed by the power operation) are needed. In the phase-velocity direction (yaxis) of a dispersion image at a certain frequency, the fundamental mode, the first higher-mode, the second higher-mode and so on, will successively appear. Accordingly, the strength of the imaged dispersion energy will ascend and descend repeatedly (like peaks and troughs). Thus, a series of independent phase-velocity zones could be easily divided for each mode by the energy minima (troughs). The automatic normalization is similar to the processing of automatic gain control (AGC). It searches phase-velocity ranges of every mode of dispersion energy and normalize them separately: C ð f ; vÞnorm ¼
Aoriginal ð f ; vÞ ; k ¼ 1; 2; 3…; max Aoriginal ð f ; vÞ
ð14Þ
vlow ðkÞbvbvhigh ðkÞ
where Cnorm is the normalized coefficient matrix of Aoriginal, k represents the fundamental (k = 1), the first higher (k = 2) and the second higher mode (k = 3) and so forth, vlow(k) and vhigh(k) are the lower and higher limits of phase-velocity range of mode k. In Cnorm, energy peaks of different modes are all raised to 1 (although their original values are different), so they will always remain and equally highlighted during the power operation. By multiplying the sparse coefficient matrix Cαnorm, the resolution of dispersion image can also be enhanced and sharper dispersion energy is achieved: Aenhanced ¼ Aoriginal C αnorm :
ð15Þ
Since the dispersion energy in a low-frequency range possesses lower resolution than in a high-frequency range (Forbriger, 2003), a variant exponent that is correlated with frequency should be chosen so as to generate a dispersion image with the same resolution in all frequency ranges. The exponent α can be written as b αð f Þ ¼ a þ ; f
ð16Þ
Fig. 1. Synthetic Love wave data due to a five-layer model simulated by a finite-difference method (Luo et al., 2010). We selected a spatial grid size of 0.5 m and a time step of 0.1 ms in the modeling. The source was described by a Ricker wavelet with 20 Hz dominant frequency (with a 50 ms delay).
decomposition and slant stacking, respectively (Fig. 2a, d, g, j). Dispersion images that generated by standard (iterate 1 time) and highresolution (iterate 2 and 3 times) LRT are appended (Fig. 2m–o) as a comparison. They present strong and clear dispersion energy of the fundamental and higher-mode Love waves. The energy peaks of Love waves in these images are continuous and match the theoretical multi-mode dispersion curves (Schwab and Knopoff, 1972) perfectly. Multi-mode dispersion energy is sufficiently separated and can be easily distinguished. Thus, multi-mode dispersion curves can be accurately picked. Resolution of dispersion energy that imaged by the four lowresolution methods (Fig. 2a, d, g, j) and the standard LRT (Fig. 2m) is very close. But specifically, resolution of high-resolution LRT (Fig. 2n–o) is evidently higher than them (Table 2). After performing the resolution-enhancing step, eight enhanced results (Fig. 2b, c, e, f, h, i, k, l) are generated. With appropriate strength of power operation, resolution of the four enhanced methods turns out to be virtually the same (Fig. 2b, e, h, k compared with Fig. 2n; Fig. 2c, f, i, l compared with Fig. 2o) with that of high-resolution LRT (it is meaningless to compare the detailed value of the enhanced resolution because it could be endlessly enhanced if we use a larger exponent). Consequently, combined with the enhancing processing, the four methods and high-resolution LRT generated dispersion images with equivalent resolution. 4. Discussion
where a and b are two parameters that control the strength of power operation. The first parameter a is the minimum exponent of the power that is applied in a high frequency range. The second parameter b decides the level of the extra power in a low frequency range. The higher the requested resolution is, the larger value of a should be chosen. If higher-resolution is required at a low frequency range, b is ought to be raised. 3. Results High frequency Love waves (Fig. 1) due to a five-layer model (Table 1) are simulated using a finite-difference method (Luo et al., 2010). Dispersion images of this noise-free data are generated by τ-p transformation, f-k transformation, phase shift, and frequency
In fact, the five kinds of algorithms (the four regular methods and standard LRT) own the same capability to obtain the dispersion kernel
Table 1 Parameters of a five-layer model. Layer
Shear wave velocity (m/s)
Density (g/cm3)
Thickness (m)
1 2 3 4 Half-space
200 300 400 500 600
1.75 1.80 1.85 1.90 1.95
2.0 3.0 4.0 5.0 ∞
170
C. Shen et al. / Journal of Applied Geophysics 122 (2015) 167–171
Fig. 2. Dispersion images generated by (a–c) τ-p transformation, (d–f) f-k transformation, (g–i) phase shift, (j–l) frequency decomposition and slant stacking, (m) standard LRT and (n–o) high-resolution LRT. Fig. 2b, e, h, k (choosing a = 4 and b = 25 in Eq. (16), the exponent of power gradually changes from 4.25 at 100 Hz to 6.5 at 10 Hz) and Fig. 2c, f, i, l (choosing a = 8 and b = 200 in Eq. (16), the exponent of power gradually changes from 10 at 100 Hz to 28 at 10 Hz) are all enhanced result of Fig. 2a, d, g, j, respectively. The iteration number of LRT differs in Fig. 2m–o (corresponding iteration numbers are 1, 2, and 3, respectively). The black dots represent theoretical dispersion curves calculated by a Knopoff's method (Schwab and Knopoff, 1972). In each column, resolution of dispersion energy that is imaged by the five kinds of methods is virtually the same.
(exact phase-velocity dispersion curves) for noise-free surface-wave data. These original algorithms, however, merely have elementary resolving ability to image the dispersion kernel, resulting in their low resolution. The sparse inversion technique makes the resolving ability of high-resolution LRT stronger than standard LRTs, so it generates higher-resolution dispersion images. The resolution-enhancing
processing is able to enhance resolution as well because the power operation is akin to iteration and weighting in high-resolution LRT. With the exponent of power increasing, its resolving ability gets stronger and the dispersion kernel is better highlighted. Therefore, the four enhanced methods are fully capable of generating competing highresolution dispersion images as high-resolution LRT.
C. Shen et al. / Journal of Applied Geophysics 122 (2015) 167–171 Table 2 Quantitative statistics of resolution (in m/s) of dispersion energy in Fig. 2. The resolution results of τ-p transformation, f-k transformation, phase shift, frequency decomposition and slant stacking refer to Fig. 2a, d, g, j, respectively. The resolution result of high-resolution linear Radon transformation refers to Fig. 2o. Dispersion-imaging method
Fundamental mode at 20 Hz
First higher mode at 50 Hz
Second higher mode at 80 Hz
τ-p transformation F-k transformation Phase shift Frequency decomposition and slant stacking High-resolution linear Radon transformation
74.5 76.8 72.1 76.9
37.2 37.1 34.3 36.6
25.4 26.2 26.7 27.3
11.5
10.1
9.3
For real-world data, which inevitably contains noise, and the response of the five dispersion- imaging methods could be different. A low signal to noise ratio will give rise to the improper functioning of the high-resolution LRT. In this situation, it is obliged to sacrifice resolution to ensure the accuracy due to its inversion methodology. So highresolution LRT is not always the best choice in the real world. The enhanced methods based on numeric processing are relatively more stable and could still generate high-resolution dispersion images. The detailed comparison of resolution of the five methods for real-world data is on the way. 5. Conclusions The algorithms of τ-p transformation, f-k transformation, phase shift, frequency decomposition and slant stacking, and high-resolution LRT are compared in terms of resolution in the f–v domain. Though the original resolution of the four regular methods is lower than the high-resolution LRT, their enhanced approaches that combined with an additional processing step, are able to generate high-resolution dispersion images as well. For noise-free (synthetic) data, resolution of the five methods is equivalent in substance. Acknowledgments The first author thanks Jianghai Xia of China University of Geosciences (Wuhan), Hubei, and Chin-Ping Lin of National Chiao Tung University,
171
Taiwan for initiating discussion on this topic. This work is supported by the National Natural Science Foundation of China (NSFC), under Grant No. 41274142. References Coruh, C., 1985. Stretched automatic amplitude adjustment of seismic data. Geophysics 50 (2), 52–256. Forbriger, T., 2003. Inversion of shallow-seismic wavefields: I. Wavefield transformation. Geophys. J. Int. 153, 719–734. Foster, D.J., Mosher, C.C., 1992. Suppression of multiple reflections using the Radon transform. Geophysics 57, 386–395. Herrmann, P., Mojesky, T., Hugonnet, P., 2000. Dealiased high-resolution Radon transforms. 70th AnnInternat. Meeting, SEG, Expanded Abstracts, pp. 1953–1956. Ji, J., 2006. CGG method for robust inversion and its application to velocity-stack inversion. Geophysics 71 (4), R59–R67. Luo, Y., Xia, J., Miller, R.D., et al., 2008. Rayleigh-wave dispersive energy imaging by highresolution linear Radon transform. Pure Appl. Geophys. 165 (5), 903–922. Luo, Y., Xia, J., Xu, Y., et al., 2010. Finite-difference modeling and dispersion analysis of high-frequency Love waves for near-surface applications. Pure Appl. Geophys. 167 (12), 1525–1536. McMechan, G.A., Yedlin, M.J., 1981. Analysis of dispersive waves by wave field transformation. Geophysics 46, 869–874. Moro, G.D., Pipan, M., Forte, E., et al., 2003. Determination of Rayleigh wave dispersion curves for near surface applications in unconsolidated sediments. Technical Program with Biographies, SEG, 73rd Annual Meeting, Dallas, TX, pp. 1247–1250. Park, C.B., Miller, R.D., Xia, J., 1998. Imaging dispersion curves of surface waves on multichannel record. Technical Program with Biographies, SEG, 68th Annual Meeting, New Orleans, Louisiana, pp. 1377–1380. Sacchi, M., Ulrych, T., 1995. High resolution velocity gathers and offset space reconstruction. Geophysics 60, 1169–1177. Schwab, F.A., Knopoff, L., 1972. Fast surface wave and free mode computations: in methods in computational physics, edited by B.A. Bolt: Academic Press, New York, 87–180 Thorson, J.R., Claerbout, J.F., 1985. Velocity stack and slant stochastic inversion. Geophysics 50, 2727–2741. Trad, D., Ulrych, T., Sacchi, M., 2003. Latest views of the sparse Radon transform. Geophysics 68, 386–399. Xia, J., Xu, Y., Chen, C., et al., 2006. Simple equations guide high-frequency surface-wave investigation techniques. Soil Dyn. Earthq. Eng. 26 (5), 395–403. Xia, J., Xu, Y., Miller, R.D., 2007. Generating image of dispersive energy by frequency decomposition and slant stacking. Pure Appl. Geophys. 164 (5), 941–956. Xia, J., 2014. Estimation of near-surface shear-wave velocities and quality factors, using multichannel analysis of surface-wave methods. J. Appl. Geophys. 103, 140–151. Yilmaz, Ö., 1987. Seismic Data Processing. Society of Exploration Geophysicists, Tulsa, OK, p. 526.
Contents lists available at ScienceDirect
Journal of Applied Geophysics journal homepage: www.elsevier.com/locate/jappgeo
Resolution equivalence of dispersion-imaging methods for noise-free high-frequency surface-wave data Chao Shen a,b,⁎, Ao Wang a, Limin Wang a,b,c, Zongbo Xu a,b, Feng Cheng a,b a b c
Institute of Geophysics and Geomatics, China University of Geosciences (Wuhan), 388 Lumo Rd., Wuhan, Hubei 430074, China Subsurface Imaging and Sensing Laboratory, Institute of Geophysics and Geomatics, China University of Geosciences, 388 Lumo Rd., Wuhan, Hubei 430074, China Hubei Subsurface Multi-scale Imaging Key Laboratory, China University of Geosciences (Wuhan), 388 Lumo Rd., Wuhan, Hubei 430074, China
a r t i c l e
i n f o
Article history: Received 19 July 2015 Received in revised form 2 September 2015 Accepted 22 September 2015 Available online 25 September 2015 Keywords: Surface wave Dispersion image Resolution Equivalence
a b s t r a c t A key step in high-frequency surface-wave methods is to acquire the dispersion properties of surface waves. To pick dispersion curves, surface-wave data is required to be transformed from the time–space (t–x) domain to the frequency–phase velocity (f–v) domain. In constructing accurate multi-mode dispersion curves, it is necessary to generate a reliable and high-resolution image of dispersion energy in the f–v domain. At present, there are five methods available for imaging dispersion energy: the τ-p transformation, the f-k transformation, the phase shift, the frequency decomposition and slant stacking, and the high-resolution linear Radon transformation. Among them, resolution of the high-resolution linear Radon transformation is the highest. We proposed a processing step to enhance resolution of the other four methods. Resolution of the four enhanced methods and high-resolution linear Radon transformation is compared by imaging dispersion energy of the same synthetic data. As a result, the four enhanced methods generated dispersion images with the same resolution, which revealed that the resolution of the five dispersion-imaging methods is essentially equivalent for noise-free data. © 2015 Elsevier B.V. All rights reserved.
1. Introduction Surface waves are dispersive for layered earth models. Particular modes of surface waves will possess unique phase velocities at different frequencies. A key step in high-frequency surface-wave methods is to acquire the dispersion properties of surface waves. Generating a reliable and high-resolution image of dispersion energy in frequency–phase velocity (f–v) domain is of essential necessity (Xia, 2014). From a highresolution dispersion image, accurate dispersion curves can be picked by following peaks of dispersion energy at different frequencies, leading to an effective inversion. In high-frequency surface-wave methods, there are five methods available for imaging dispersion energy at present: the τ-p transformation (McMechan and Yedlin, 1981), the f-k transformation (e.g., Yilmaz, 1987), the phase shift (Park et al., 1998), the frequency decomposition and slant stacking (Xia et al., 2007), and the high-resolution linear Radon transformation (Luo et al., 2008). Among them, high-resolution linear Radon transformation usually generates higher-resolution dispersion images than the other four methods do (Luo et al., 2008). Since the method has been utilized in many applications (Herrmann et al., 2000; Foster and Mosher, 1992; Thorson and Claerbout, 1985), it has received more and more attentions. Based on our current works, ⁎ Corresponding author at: Institute of Geophysics and Geomatics, China University of Geosciences (Wuhan), 388 Lumo Rd., Wuhan, Hubei 430074, China. E-mail address: [email protected] (C. Shen).
http://dx.doi.org/10.1016/j.jappgeo.2015.09.019 0926-9851/© 2015 Elsevier B.V. All rights reserved.
we find that there are ways to enhance resolution of the four other methods, so that the five methods can possess the same resolution power in the f–v domain. One simple practice is the numeric processing that based on power operation. In this paper, we first introduce principles of the five dispersionimaging methods and are solution-enhancing processing. Then dispersion energy of the same synthetic surface-wave shot gather is imaged using all the five methods. Besides, the processing step is executed to enhance resolution of dispersion images that are generated by τ-p transformation, f-k transformation, phase shift, and frequency decomposition and slant stacking. Finally, we conduct a comparison to summarize resolution of the five dispersion-imaging methods.
2. Methods To acquire dispersion energy in the f–v domain, a shot gather in the time–space (t–x) domain is necessary to go through two steps: one is to transform from the time domain to the frequency domain; the other one is to transform from the space domain to the phase velocity domain by velocity scanning and slant stacking. These two steps are independent and they can switch their transformation sequences with each other during calculation. Let s(x, t) be surface-wave wavefield, where x is offset (distance between a source and a receiver) and t is travel time, all receivers are evenly arranged along a survey line with a spacing Δx.
168
C. Shen et al. / Journal of Applied Geophysics 122 (2015) 167–171
2.1. τ-p transformation (McMechan and Yedlin, 1981) A slant stack operator is first applied to linearly transform wavefield s(x, t) into τ-p domain data: U ðp; τÞ ¼
N X
sðxi ; t ¼ τ þ pxi Þ;
ð1Þ
i¼1
where p is the slowness and τ is the intercept on the time axis. If the scanning slowness p is close to the real slowness of surface wave, U(p, τ) will reach a maximum. A Fourier transformation is then applied to the data in the τ domain to transform them into the frequency domain. As p is simply the reciprocal of phase velocity (p = 1/v), the f–v domain result can be easily converted and achieved.
Z
∞ −∞
∞ −∞
sðx; t Þe−2πiðftþkxÞ dtdx:
ð2Þ
This procedure could be replaced with a one-dimensional Fourier transformation in the time domain and a one-dimensional inverse Fourier transformation in the space domain. Then the f–v domain result can be transformed according to the relation: v = f/k. 2.3. Phase shift (Park et al., 1998) Phase information P( x,ω) of a wavefield can be calculated by a Fourier transformation in the time domain. It is expressed as ω
P ð x; ωÞ ¼ e−i v xþφ0 ;
ð3Þ
where ω is the circular frequency, v is the phase velocity and φ0 is the initial phase. The phase is linearly correlated to offset. Integrate the phase using a scanning φ: Z U ðφ; ωÞ ¼
Z eiφx P ðx; ωÞdx ¼
eiðφ−
ω ν
Þxþφ0
dx ¼
N X
eiðφ−
ω ν
Þxi þφ0
:
When the scanning φ gets close to ω/v, phase angle of the complex ω numbereiðφ− v Þxþφ0 is the same in the complex plane, making the integral reaching a maximum value. At last, the f–v domain result can be achieved according to the following relations: f ¼ 2πω; and ω 2πf 2π ¼ ¼ f: φ φ φ
ð5Þ ð6Þ
Phase shift method utilizes the normalized spectrum and can perform well in generating dispersion image (Moro et al., 2003). Nevertheless, due to the uniqueness of calculated phase at each frequency, dispersion energy of higher modes may not be completely imaged. 2.4. Frequency decomposition and slant stacking (Xia et al., 2007) Shot gather s(x, t) is first stretched into pseudo-vibroseis data d(x,t) using a frequency decomposition technique (Coruh, 1985): dðx; t Þ ¼ yðt Þ⊕sðx; t Þ;
n X
mðpi ; f Þei2πfpi x ;
ð8Þ
where d(x,f) is a shot gather in the frequency domain and m(p, f) is the Radon matrix. So Eq. (8) can be written as a matrix form: d ¼ Lm;
ð9Þ
where L = ei2πfpx is the linear Radon transformation (LRT) operator. Eq. (9) can be rearranged as madj ¼ LT d;
ð10Þ
where madj represents a low-resolution solution using the transpose of L. An appropriate weighting and damping factor will lead to a sparse inversion solution and improve the resolution. Thus, the objective function becomes 2 J ¼ ‖W d d−LW −1 m W m m ‖ þ λ‖W m m‖ ;
ð11Þ
and the problem turns into solving the following equation (Ji, 2006) T T −1 −T T T W −T m L W d W d LW m þ λI W m m ¼ W m L W d W d d;
ð12Þ
ð4Þ
i¼1
v¼
Standard linear Radon transformation can be expressed as the following equation
i¼1
Wavefield s(x, t) is first transformed to the f-k domain data S(k, f) by a two-dimensional Fourier transformation: Z
2.5. High-resolution linear Radon transformation (Luo et al., 2008)
dðx; f Þ ¼
2.2. F-k transformation (e.g., Yilmaz, 1987)
Sðk; f Þ ¼
So d(x, t) is also a frequency-swept data d(x, f) and each t corresponds with a frequency. A slant stacking (Yilmaz, 1987, p. 430) is then performed to obtain f–p domain data U(p, f). The f–v domain result can be eventually achieved by inverting p. The calculated result of this method is better than that of τ-p transformation because frequencies are decomposed before slant stacking. Furthermore, compared with phase shift, this algorithm could extract stronger energy for higher-mode surface waves.
ð7Þ
where ⊕ stands for the convolution operator and y(t) is a sweep function. Generally, y(t) is usually chosen as a linear function y(t) = a + bt.
where I is the identity matrix, Wd is a data weighting matrix and Wm is a model weighting matrix, λ is a damping factor and controls the compromise between the error misfit and model length. Eq. (12) can be solved by the weighted preconditioned conjugate gradient (CG) algorithm (e.g., Sacchi and Ulrych, 1995). The resolution of high-resolution LRT depends on the iteration in the CG algorithm. At the beginning of iteration, this method equals with a standard LRT (only iterate once). With the number of iteration growing, the resolution will be greatly improved (normally iterate 3 times). 2.6. A resolution-enhancing processing We define the resolution of a dispersion image as the width between the two half-values of dispersion energy at a given frequency in the f–v domain (Xia et al., 2006). τ-p transformation, f-k transformation, phase shift, frequency decomposition and slant stacking are actually kinds of a standard discretized LRT and they generate dispersion images with a low resolution (Eq. (10)). To compare them with high-resolution LRT, we deem them regular dispersionimaging methods in this paper for convenience. They suffer from incomplete information (Trad et al., 2003) and could hardly generate high-resolution dispersion images. The mathematical power operation is able to amplify the difference between large and small values, outstanding the maxima. Computing a
C. Shen et al. / Journal of Applied Geophysics 122 (2015) 167–171
169
power of the originally generated dispersion image will enhance its resolution. The explicit formula of our proposed power operation is Aenhanced ¼ Aαorignal ;
ð13Þ
where Aoriginal is a matrix of the original dispersion image generated by the four regular methods, Aenhanced is a matrix of the enhanced dispersion image with a higher resolution, and α is the exponent of the power. Resolution of the enhanced dispersion image can be controlled through adjusting α in this step. The stronger the power operation is applied to a dispersion image, the higher-resolution of dispersion energy will be imaged. High-resolution LRT images multi-mode dispersion energy simultaneously. Since the power operation will only reserve the most dominant peaks (the maximum value), an automatic normalization of dispersion energy on each frequency is preferred if higher-mode dispersion curves (normally not exhibiting the strongest energy, and being suppressed by the power operation) are needed. In the phase-velocity direction (yaxis) of a dispersion image at a certain frequency, the fundamental mode, the first higher-mode, the second higher-mode and so on, will successively appear. Accordingly, the strength of the imaged dispersion energy will ascend and descend repeatedly (like peaks and troughs). Thus, a series of independent phase-velocity zones could be easily divided for each mode by the energy minima (troughs). The automatic normalization is similar to the processing of automatic gain control (AGC). It searches phase-velocity ranges of every mode of dispersion energy and normalize them separately: C ð f ; vÞnorm ¼
Aoriginal ð f ; vÞ ; k ¼ 1; 2; 3…; max Aoriginal ð f ; vÞ
ð14Þ
vlow ðkÞbvbvhigh ðkÞ
where Cnorm is the normalized coefficient matrix of Aoriginal, k represents the fundamental (k = 1), the first higher (k = 2) and the second higher mode (k = 3) and so forth, vlow(k) and vhigh(k) are the lower and higher limits of phase-velocity range of mode k. In Cnorm, energy peaks of different modes are all raised to 1 (although their original values are different), so they will always remain and equally highlighted during the power operation. By multiplying the sparse coefficient matrix Cαnorm, the resolution of dispersion image can also be enhanced and sharper dispersion energy is achieved: Aenhanced ¼ Aoriginal C αnorm :
ð15Þ
Since the dispersion energy in a low-frequency range possesses lower resolution than in a high-frequency range (Forbriger, 2003), a variant exponent that is correlated with frequency should be chosen so as to generate a dispersion image with the same resolution in all frequency ranges. The exponent α can be written as b αð f Þ ¼ a þ ; f
ð16Þ
Fig. 1. Synthetic Love wave data due to a five-layer model simulated by a finite-difference method (Luo et al., 2010). We selected a spatial grid size of 0.5 m and a time step of 0.1 ms in the modeling. The source was described by a Ricker wavelet with 20 Hz dominant frequency (with a 50 ms delay).
decomposition and slant stacking, respectively (Fig. 2a, d, g, j). Dispersion images that generated by standard (iterate 1 time) and highresolution (iterate 2 and 3 times) LRT are appended (Fig. 2m–o) as a comparison. They present strong and clear dispersion energy of the fundamental and higher-mode Love waves. The energy peaks of Love waves in these images are continuous and match the theoretical multi-mode dispersion curves (Schwab and Knopoff, 1972) perfectly. Multi-mode dispersion energy is sufficiently separated and can be easily distinguished. Thus, multi-mode dispersion curves can be accurately picked. Resolution of dispersion energy that imaged by the four lowresolution methods (Fig. 2a, d, g, j) and the standard LRT (Fig. 2m) is very close. But specifically, resolution of high-resolution LRT (Fig. 2n–o) is evidently higher than them (Table 2). After performing the resolution-enhancing step, eight enhanced results (Fig. 2b, c, e, f, h, i, k, l) are generated. With appropriate strength of power operation, resolution of the four enhanced methods turns out to be virtually the same (Fig. 2b, e, h, k compared with Fig. 2n; Fig. 2c, f, i, l compared with Fig. 2o) with that of high-resolution LRT (it is meaningless to compare the detailed value of the enhanced resolution because it could be endlessly enhanced if we use a larger exponent). Consequently, combined with the enhancing processing, the four methods and high-resolution LRT generated dispersion images with equivalent resolution. 4. Discussion
where a and b are two parameters that control the strength of power operation. The first parameter a is the minimum exponent of the power that is applied in a high frequency range. The second parameter b decides the level of the extra power in a low frequency range. The higher the requested resolution is, the larger value of a should be chosen. If higher-resolution is required at a low frequency range, b is ought to be raised. 3. Results High frequency Love waves (Fig. 1) due to a five-layer model (Table 1) are simulated using a finite-difference method (Luo et al., 2010). Dispersion images of this noise-free data are generated by τ-p transformation, f-k transformation, phase shift, and frequency
In fact, the five kinds of algorithms (the four regular methods and standard LRT) own the same capability to obtain the dispersion kernel
Table 1 Parameters of a five-layer model. Layer
Shear wave velocity (m/s)
Density (g/cm3)
Thickness (m)
1 2 3 4 Half-space
200 300 400 500 600
1.75 1.80 1.85 1.90 1.95
2.0 3.0 4.0 5.0 ∞
170
C. Shen et al. / Journal of Applied Geophysics 122 (2015) 167–171
Fig. 2. Dispersion images generated by (a–c) τ-p transformation, (d–f) f-k transformation, (g–i) phase shift, (j–l) frequency decomposition and slant stacking, (m) standard LRT and (n–o) high-resolution LRT. Fig. 2b, e, h, k (choosing a = 4 and b = 25 in Eq. (16), the exponent of power gradually changes from 4.25 at 100 Hz to 6.5 at 10 Hz) and Fig. 2c, f, i, l (choosing a = 8 and b = 200 in Eq. (16), the exponent of power gradually changes from 10 at 100 Hz to 28 at 10 Hz) are all enhanced result of Fig. 2a, d, g, j, respectively. The iteration number of LRT differs in Fig. 2m–o (corresponding iteration numbers are 1, 2, and 3, respectively). The black dots represent theoretical dispersion curves calculated by a Knopoff's method (Schwab and Knopoff, 1972). In each column, resolution of dispersion energy that is imaged by the five kinds of methods is virtually the same.
(exact phase-velocity dispersion curves) for noise-free surface-wave data. These original algorithms, however, merely have elementary resolving ability to image the dispersion kernel, resulting in their low resolution. The sparse inversion technique makes the resolving ability of high-resolution LRT stronger than standard LRTs, so it generates higher-resolution dispersion images. The resolution-enhancing
processing is able to enhance resolution as well because the power operation is akin to iteration and weighting in high-resolution LRT. With the exponent of power increasing, its resolving ability gets stronger and the dispersion kernel is better highlighted. Therefore, the four enhanced methods are fully capable of generating competing highresolution dispersion images as high-resolution LRT.
C. Shen et al. / Journal of Applied Geophysics 122 (2015) 167–171 Table 2 Quantitative statistics of resolution (in m/s) of dispersion energy in Fig. 2. The resolution results of τ-p transformation, f-k transformation, phase shift, frequency decomposition and slant stacking refer to Fig. 2a, d, g, j, respectively. The resolution result of high-resolution linear Radon transformation refers to Fig. 2o. Dispersion-imaging method
Fundamental mode at 20 Hz
First higher mode at 50 Hz
Second higher mode at 80 Hz
τ-p transformation F-k transformation Phase shift Frequency decomposition and slant stacking High-resolution linear Radon transformation
74.5 76.8 72.1 76.9
37.2 37.1 34.3 36.6
25.4 26.2 26.7 27.3
11.5
10.1
9.3
For real-world data, which inevitably contains noise, and the response of the five dispersion- imaging methods could be different. A low signal to noise ratio will give rise to the improper functioning of the high-resolution LRT. In this situation, it is obliged to sacrifice resolution to ensure the accuracy due to its inversion methodology. So highresolution LRT is not always the best choice in the real world. The enhanced methods based on numeric processing are relatively more stable and could still generate high-resolution dispersion images. The detailed comparison of resolution of the five methods for real-world data is on the way. 5. Conclusions The algorithms of τ-p transformation, f-k transformation, phase shift, frequency decomposition and slant stacking, and high-resolution LRT are compared in terms of resolution in the f–v domain. Though the original resolution of the four regular methods is lower than the high-resolution LRT, their enhanced approaches that combined with an additional processing step, are able to generate high-resolution dispersion images as well. For noise-free (synthetic) data, resolution of the five methods is equivalent in substance. Acknowledgments The first author thanks Jianghai Xia of China University of Geosciences (Wuhan), Hubei, and Chin-Ping Lin of National Chiao Tung University,
171
Taiwan for initiating discussion on this topic. This work is supported by the National Natural Science Foundation of China (NSFC), under Grant No. 41274142. References Coruh, C., 1985. Stretched automatic amplitude adjustment of seismic data. Geophysics 50 (2), 52–256. Forbriger, T., 2003. Inversion of shallow-seismic wavefields: I. Wavefield transformation. Geophys. J. Int. 153, 719–734. Foster, D.J., Mosher, C.C., 1992. Suppression of multiple reflections using the Radon transform. Geophysics 57, 386–395. Herrmann, P., Mojesky, T., Hugonnet, P., 2000. Dealiased high-resolution Radon transforms. 70th AnnInternat. Meeting, SEG, Expanded Abstracts, pp. 1953–1956. Ji, J., 2006. CGG method for robust inversion and its application to velocity-stack inversion. Geophysics 71 (4), R59–R67. Luo, Y., Xia, J., Miller, R.D., et al., 2008. Rayleigh-wave dispersive energy imaging by highresolution linear Radon transform. Pure Appl. Geophys. 165 (5), 903–922. Luo, Y., Xia, J., Xu, Y., et al., 2010. Finite-difference modeling and dispersion analysis of high-frequency Love waves for near-surface applications. Pure Appl. Geophys. 167 (12), 1525–1536. McMechan, G.A., Yedlin, M.J., 1981. Analysis of dispersive waves by wave field transformation. Geophysics 46, 869–874. Moro, G.D., Pipan, M., Forte, E., et al., 2003. Determination of Rayleigh wave dispersion curves for near surface applications in unconsolidated sediments. Technical Program with Biographies, SEG, 73rd Annual Meeting, Dallas, TX, pp. 1247–1250. Park, C.B., Miller, R.D., Xia, J., 1998. Imaging dispersion curves of surface waves on multichannel record. Technical Program with Biographies, SEG, 68th Annual Meeting, New Orleans, Louisiana, pp. 1377–1380. Sacchi, M., Ulrych, T., 1995. High resolution velocity gathers and offset space reconstruction. Geophysics 60, 1169–1177. Schwab, F.A., Knopoff, L., 1972. Fast surface wave and free mode computations: in methods in computational physics, edited by B.A. Bolt: Academic Press, New York, 87–180 Thorson, J.R., Claerbout, J.F., 1985. Velocity stack and slant stochastic inversion. Geophysics 50, 2727–2741. Trad, D., Ulrych, T., Sacchi, M., 2003. Latest views of the sparse Radon transform. Geophysics 68, 386–399. Xia, J., Xu, Y., Chen, C., et al., 2006. Simple equations guide high-frequency surface-wave investigation techniques. Soil Dyn. Earthq. Eng. 26 (5), 395–403. Xia, J., Xu, Y., Miller, R.D., 2007. Generating image of dispersive energy by frequency decomposition and slant stacking. Pure Appl. Geophys. 164 (5), 941–956. Xia, J., 2014. Estimation of near-surface shear-wave velocities and quality factors, using multichannel analysis of surface-wave methods. J. Appl. Geophys. 103, 140–151. Yilmaz, Ö., 1987. Seismic Data Processing. Society of Exploration Geophysicists, Tulsa, OK, p. 526.