地球物理学报  2013, Vol. 56 Issue (8): 2790-2798   PDF    
Hilbert谱分析在探地雷达薄层识别中的应用
张先武1,2 , 高云泽1 , 方广有1     
1. 中国科学院电子学研究所, 电磁辐射与探测技术中国科学院重点实验室, 北京 100190;
2. 中国科学院大学, 北京 100049
摘要: 针对实际勘探中的需要, 本文对薄层的电磁波反射特性进行了深入研究.薄层在反射电磁波时, 会对入射电磁波进行滤波作用.薄层的类型和厚度会影响其滤波特性, 利用此特点, 本文将Hilbert谱分析方法引入到薄层识别中.该方法能将探地雷达对递变型薄层的垂直分辨率提升至λ/8.本文利用该方法对实测探地雷达资料进行了层位识别研究, 该方法能很好地提高探地雷达的层位识别能力.
关键词: 探地雷达      薄层识别      Hilbert谱分析     
Application of Hilbert spectrum analysis in Ground Penetrating Radar thin layer recognition
ZHANG Xian-Wu1,2, GAO Yun-Ze1, FANG Guang-You1     
1. Key Laboratory of Electromagnetic Radiation and Sensing Technology of CAS, Institute of Electronics, Chinese Academy of Sciences, Beijing 100190, China;
2. University of Chinese Academy of Sciences, Beijing 100049, China
Abstract: This paper makes a deep study of thin layer's electromagnetic reflection properties in accordance with the need of actual exploration. The thin layer function is a filter for incident electromagnetic wave when the thin layer reflects electromagnetic wave. For thin layer's type and thickness influence its filtering properties, the paper introduces Hilbert Spectrum Analysis method to recognize the thin layer by exploiting this character. The method can help improve Ground Penetrating Radar's vertical resolution on the graded thin layer up to λ/8.This paper employs such method to make a layer recognition study on actual Ground Penetrating Radar data, by which layer recognition ability of Ground Penetrating Radar can be well improved..
Key words: Ground Penetrating Radar      Thin layer recognition      Hilbert Spectrum Analysis     
1 引言

探地雷达是一种利用高频脉冲电磁波来探测介质内部物性分布规律的地球物理方法[1-4].由于探地雷达辐射出来的电磁波频率高,波长短,所以探地雷达在探测常见的分层介质(如探测冰川内部沉积结构[5-7]、公路路面层和基底结构[8-10]、近地表岩石分层[11-12])时能获得很高的垂直分辨率.借用地震反射理论,依据瑞雷标准,探地雷达所能识别最薄层状介质的厚度为λ/4(λ为介质中电磁波中心频率对应的波长)[13].实际勘探中,分层介质内部通常含有层厚小于λ/4的薄层,为了更精细地划分分层介质内部结构,需要对层厚小于λ/4的薄层进行识别.因此,研究层厚小于λ/4的薄层对电磁波的反射特性具有重要意义.

电磁波在分层介质中传播时,如果分界面两边介质的波阻抗不连续,电磁波在分界面上会产生反射.本文首先建立含有薄层的分层介质模型,推导出薄层反射系数.再通过分析薄层反射系数来分析薄层的电磁波反射特性.依据薄层反射系数分析结果,得出:薄层对电磁波的反射,可以等效为薄层对电磁波的滤波作用.薄层的滤波特性与薄层的类型和厚度密切相关.

由于薄层的滤波效应,电磁波经过薄层反射后,频谱会随着薄层厚度的变化而变化.根据此特点,本文将Hilbert谱分析方法[14-16]引入到薄层识别中.模型计算结果表明,Hilbert谱分析方法可以识别出层厚大于λ/8的递变型薄层.最后,本文利用该方法对实测探地雷达数据进行了层位识别处理,处理结果表明Hilbert谱分析方法能很好地提高探地雷达的层位识别能力.

2 薄层的电磁波反射特性 2.1 薄层反射系数分析

为了研究薄层的电磁波反射特性,建立如图 1所示的三层层状介质模型.三层介质的电导率都为0.第1层和第3层介质分别为半无限大空间,介电常数和磁导率分别为ε1μ1ε3μ3.第2层为薄层,介电常数和磁导率分别为ε2μ2.薄层顶界面深度为0,底界面深度为h,即薄层厚度为h.

图 1 三层层状介质模型 Fig. 1 Three layers layered medium model

假定电磁波以TE平面波入射到薄层顶界面上,入射波的表达式为:Ey=E0e-ikzz+ikxxkzkx分别为入射波沿z方向和x方向的波矢量).薄层的反射系数R[17]

(1)

(1)式中,k1zk2zk3z分别表示的是TE平面波在第1层、薄层和第3层中沿z方向的波矢量.

,则(1)式可以简化为:

(2)

R12R23的表达式可知,R12表示的是薄层顶界面在顶层中引起的波的反射系数,R23表示的是薄层底界面在薄层中引起的波的反射系数.

令薄层顶界面在顶层中引起的波的透射系数为T12,在薄层中引起的波的反射系数和透射系数分别为R21T21.则R12T12R21T21满足以下关系:

(3)

(4)

(5)

利用(3)、(4)、(5)式,可将(2)式改写为:

(6)

因为|R21R23| < 1,所以

(7)

由(6)、(7)式,得到R的级数表达式为:

(8)

R的级数表达式可以看出,薄层反射电磁波是薄层顶界面的单次反射和薄层内部多次反射综合作用的结果(图 2).

图 2 薄层反射电磁波示意图 Fig. 2 Diagram of thin layer reflecting electromagnetic wave

实际勘探中,|T12R232R21T21|一般远小于|T12R23T21|,R可近似表示为:RR12 +T12R23T21ei2k2zh.假定入射的平面电磁波为单色平面波(频率为f),因为薄层为理想电介质,所以k2z=2πf/vv为电磁波在薄层中的传播速度,.R的近似表达式可改写为:RR12+T12R23T21ei2π0,其中τ0=2h/v,为电磁波在薄层中的双程走时.

R的近似表达式进行化简,化简结果为:

(9)

(9)式中:

(10)

(11)

(12)

由反射系数R,得到薄层反射的电磁波Ery为:

(13)

由(13)式可以看出,薄层反射的电磁波的振幅和相位与Pθ)、Qθ)密切相关.Pθ)、Qθ)由薄层顶、底界面的反射系数、薄层厚度、介电常数、磁导率、入射波频率这些参数决定.当薄层顶、底界面反射系数、薄层厚度、介电常数、磁导率确定时,反射波的振幅会随着入射波频率的变化而变化.当入射波不是单色波时,入射波经过薄层反射后,反射波的频谱会发生改变.因此,薄层在反射电磁波时,可以将薄层等效为一个滤波器.薄层反射电磁波,可以看作薄层对电磁波进行了滤波.

2.2 薄层滤波特性分析

由前面分析可知,在反射电磁波时,薄层可以等效为一个滤波器,Pθ)、Qθ)可以视为等效滤波器的振幅特性和相位特性.本文主要研究的是等效滤波器的振幅特性Pθ).

借用地震勘探中对薄层类型的划分标准[18],可以将薄层划分为递变型薄层(R12R23>0)和韵律型薄层(R12R23 < 0)两种类型.

2.2.1 递变型薄层滤波特性分析

因为R12R23 >0,(1-R122)≥0,所以当cosθ=1,即0=mm为整数)时,Pθ)取极大值Pmax

当cosθ=-1,即0=m+1/2(m为整数)时,Pθ)取极小值Pmin

PmaxPmin对应的频点会随着τ0的变化而变化.当薄层介电常数和磁导率确定时,τ0与薄层厚度成线性关系,此时,等效滤波器的振幅特性由薄层厚度决定.图 3为不同厚度递变型薄层滤波特性.递变型薄层对应的模型参数为:ε1=ε0ε2=2ε0ε3=4ε0μ1=μ2=μ3=μ0.ε0μ0分别为真空中的介电常数和磁导率.薄层厚度h分别设置为1.4、0.7、0.35、0.175m.

图 3 不同厚度递变型薄层的滤波特性 (a)h=1.4m;(b)h=0.7m;(c)h=0.35m;(d)h=0.175m. Fig. 3 Filter property of graded thin layer with different thickness

图 3中可以看出,等效滤波器的零频对应的为最大振幅值点Pmax.Pmax和最小振幅值点Pmin周期性出现.随着薄层厚度的减小,Pmin对应的频点逐渐变大,PmaxPmin之间对应的频宽fB逐渐变宽.当fB与入射波截止频率fC接近时,反射波的高频成分会被压制,薄层滤波特性与低通滤波器相似.当fB接近于fC/2时,反射波的中间频率成分会被压制,薄层滤波特性类似于带阻滤波器.

图 5为不同厚度递变型薄层(薄层厚度分别取:1.4、0.7、0.35、0.175 m)反射电磁波的一维FDTD模拟结果.模拟时,入射电磁波选择的是一阶高斯脉冲源(图 4),源的中心频率f0为150 MHz.源的空间位置与薄层顶界面的垂直距离为2.5m,波场记录点与源点置于同一空间位置.空间网格ds取0.01m,时间网格dt取0.01ns.边界采用的是PML边界.

图 4 一阶高斯脉冲源的(a)时域波形和(b)频谱 Fig. 4 Time domain waveform (a) and frequency spectrum (b) of Gaussian pulse's first derivative
图 5 不同厚度递变型薄层反射电磁波的一维FDTD模拟结果 (a)h=1.4m;(b)h=0.7m;(c)h=0.35m;(d)h=0.175m. Fig. 5 One dimensional FDTD simulation results of graded thin layer reflecting electromagnetic wave with different thickness

图 5c图 5d对应的频谱图为图 6a图 6b.当薄层厚度取0.35、0.175 m时,fB分别接近于fC/2、fC.对比分析图 4b图 6,在图 6a中,反射波的中间频率成分被压制.而在图 6b中,反射波的高频成分被压制.FDTD模拟结果与递变型薄层滤波特性理论分析结果一致.

图 6 不同厚度递变型薄层反射电磁波频谱图 (a)h=0.35m;(b)h=0.175m. Fig. 6 Reflected electromagnetic wave frequency spectrum of graded thin layer with different thickness
2.2.2 韵律型薄层滤波特性分析

因为R12R23 < 0,(1-R122)≥0,所以

当cosθ=-1,即0=m+1/2(m为整数)时,Pθ)取极大值Pmax

当cosθ=1,即0=mm为整数)时,Pθ)取极小值Pmin

图 7为不同厚度韵律型薄层滤波特性,韵律型薄层对应的模型参数为:ε1=2ε0ε2=ε0ε3=4ε0μ1=μ2=μ3=μ0.薄层厚度h分别设置为2、1、0.5、0.25m.当薄层介电常数和磁导率确定时,等效滤波器的振幅特性也是由薄层厚度决定.等效滤波器的零频对应的为最小振幅值点Pmin.Pmin和最大振幅值点Pmax也是周期性出现.随着薄层厚度的减小,Pmax对应的频点逐渐变大,PminPmax之间对应的频宽fB逐渐变宽.当fB与入射波截止频率fC接近时,反射波的低频成分被压制,薄层滤波特性与高通滤波器相似.当fB接近于fC/2时,反射波的高频和低频成分被压制,薄层滤波特性类似于带通滤波器.图 8a图 8b为FDTD模拟的韵律型薄层(薄层厚度分别取0.5、0.25m)反射电磁波频谱图.图 4为模拟时选用的脉冲源.

图 7 不同厚度韵律型薄层的滤波特性 (a)h=2m;(b)h=1m;(c)h=0.5m;(d)h=0.25m. Fig. 7 Filter property of rhythmicity thin layer with different thickness
图 8 不同厚度韵律型薄层反射电磁波频谱图 (a)h=0.5m;(b)h=0.25m. Fig. 8 Reflected electromagnetic wave frequency spectrum of rhythmicity thin layer with different thickness

当薄层厚度取0.5 m、0.25 m时,fB分别接近于fC/2、fC.对比分析图 4b图 8,在图 8a中,反射波的低频成分和高频成分都被压制.而在图 8b中,反射波的低频成分被压制.FDTD模拟结果验证了韵律型薄层滤波特性理论分析的正确性.

3 Hilbert谱分析识别薄层

由第2节分析结果,薄层反射波频谱会随着薄层厚度的变化而变化.而薄层厚度与薄层顶、底界面反射波的双程走时密切相关.根据此特点,可以通过分析薄层反射波频谱随时间的变化来分析薄层顶、底界面反射波的双程走时,进而实现对薄层厚度和深度的识别.

为了分析薄层反射波频谱随时间的变化,需要借助时频分析方法.常用的时频分析方法有短时Fourier变换[19],小波变换[20-21],S变换[22-24],Hilbert谱分析等.同其它时频分析方法相比,Hilbert谱分析能在时频域获得较高的时频分辨率.本文选用Hilbert谱分析方法对薄层进行了识别研究,该方法能识别出层厚大于λ/8的递变型薄层.

3.1 Hilbert谱分析

对时变信号xt)进行Hilbert变换,得到解析信号qt),

式中,ht)为xt)的Hilbert变换.

xt)的瞬时振幅at)、瞬时相位θt)、瞬时频率ft)分别为:

在时间-频率平面内,利用at)和ft)来表示xt),得到xt)的Hilbert谱Htf),

3.2 Hilbert谱分析识别递变型薄层

对递变型薄层反射电磁波的一维FDTD模拟结果(图 5)进行Hilbert谱分析.图 9为Hilbert谱分析结果.

图 5中,入射波的中心频率f0为150 MHz,f0在薄层中对应的波长λ为1.4 m,薄层厚度分别为λλ/2、λ/4、λ/8.随着薄层厚度的减小,薄层顶、底界面的反射波相互叠加.这种叠加会影响薄层顶、底界面双程走时的拾取,进而影响薄层厚度的识别.当薄层厚度减小为λ/8时,在图 5中拾取不出薄层底界面的双程走时.此时,已经无法识别出薄层厚度.

对比分析图 9图 5图 9中的直达波、薄层反射波能很好地与图 5中对应.通过拾取时频平面上的能量极大值点,可以准确地拾取薄层顶、底界面的双程走时.随着薄层厚度的减小,薄层顶、底界面反射波的叠加会使薄层反射波在时频平面上的能量轨迹发生改变.当薄层厚度接近λ/4时,薄层反射波在时频平面上的频率值会异常增高.当薄层厚度减小为λ/8时,在时频平面上,依然可以准确地拾取出薄层顶、底界面的双程走时.同时在图 9中能清晰地看见FDTD模拟时出现的噪声和多次波的能量轨迹,但是它们在时频平面上的能量很弱.从时频平面上的能量角度看,相对于直达波和薄层反射波,噪声和多次波可以被忽略.

图 9 不同厚度递变型薄层反射电磁波的Hilbert谱分析结果 (a)h=1.4m;(b)h=0.7m;(c)h=0.35m;(d)h=0.175m. Fig. 9 Reflected electromagnetic wave Hilbert Spectrum Analysis results of graded thin layer with different thickness
4 Hilbert谱分析在实测探地雷达资料层位识别中的应用

图 10为某一湖区实测探地雷达资料.测量时,天线中心频率选用150 MHz,时间窗口为500ns,采样点数为512,总的记录道数为141道.

图 10 实测探地雷达资料 Fig. 10 Ground Penetrating Radar data

选取第80道数据(图 11),对其进行Hilbert谱分析.在时频平面内,利用Hilbert谱分析结果,对能量轨迹进行曲线拟合,得到第80道数据时频曲线(图 12).

图 11 第80道实测探地雷达数据 Fig. 11 The 80th trace Ground Penetrating Radar data
图 12 第80道实测探地雷达数据时频曲线 Fig. 12 Time-frequency curve of the 80th trace Ground Penetrating Radar data

图 11中可以看出,雷达子波时域波形存在两个较强的旁瓣.在利用湖底淤泥沉积层回波信号(图 11中A区域)划分湖底淤泥沉积层层位时,旁瓣会给层位的划分带来很大困难.而在图 12中A区域,我们可以清晰地将湖底淤泥沉积层划分为四个独立反射层.同时,依据湖底淤泥沉积层底部回波信号频率异常增高这个特点,可以预测湖底淤泥沉积层底层厚度接近1/(2εr1/2)(εr为湖底淤泥的相对介电常数).

图 12中B区域,可以清晰地看见湖水回波信号的能量轨迹.但是在时频平面内,湖水回波信号的能量比湖底淤泥沉积层回波信号的能量小很多.利用Hilbert谱分析方法对实测雷达信号进行层位划分时,能量轨迹能很好地识别出弱信号,但是对噪声也很敏感.所以,对实测雷达信号进行Hilbert谱分析后,在时频平面内需要综合利用信号的能量轨迹和能量大小来划分层位.

对141道数据分别进行Hilbert谱分析.在时频平面内,利用每道数据的Hilbert谱分析结果,对能量轨迹进行曲线拟合,得到每道数据的时频曲线.最后将141道数据的时频曲线显示在同一平面内(图 13).

图 13 实测探地雷达资料Hilbert谱分析结果 Fig. 13 Hilbert Spectrum Analysis result of Ground Penetrating Radar data

对比分析图 13图 10.在图 13中,湖底淤泥沉积层(A区域)层位更加清晰.多道数据Hilbert谱分析结果也表明Hilbert谱分析方法能很好地提高探地雷达的层位识别能力.

5 结论

本文对薄层的电磁波反射特性进行了深入的研究.理论分析和FDTD模拟结果都表明,不同厚度、不同类型的薄层在反射电磁波时,会表现出不同滤波特性.

本文利用薄层厚度与薄层滤波特性相关这个特点,将Hilbert谱分析引入到薄层识别中.模型计算和实测探地雷达数据处理结果证明该方法能很好地提高探地雷达的层位识别能力.

参考文献
[1] Slob E, Sato M, Olhoeft G. Surface and borehole ground-penetrating-radar developments. Geophysics , 2010, 75(5): 75A103-75A120. DOI:10.1190/1.3480619
[2] 李华, 鲁光银, 何现启, 等. 探地雷达的发展历程及其前景探讨. 地球物理学进展 , 2010, 25(4): 1492–1502. Li H, Lu G Y, He X Q, et al. The progress of the GPR and discussion on its future development. Progress in Geophysics (in Chinese) (in Chinese) , 2010, 25(4): 1492-1502.
[3] 李嘉, 郭成超, 王复明, 等. 探地雷达应用概述. 地球物理学进展 , 2007, 22(2): 629–637. Li J, Guo C C, Wang F M, et al. The summary of the surface ground penetrating radar applied in subsurface investigation. Progress in Geophysics (in Chinese) (in Chinese) , 2007, 22(2): 629-637.
[4] 底青云, 许琨, 王妙月. 衰减雷达波有限元偏移. 地球物理学报 , 2000, 43(2): 257–263. Di Q Y, Xu K, Wang M Y. The attenuated radar wave migration with finite element method. Chinese J. Geophys. (in Chinese) , 2000, 43(2): 257-263.
[5] 王璞玉, 李忠勤, 吴利华, 等. 探地雷达在冰川厚度及冰下地形探测中的应用. 吉林大学学报(地球科学版) , 2011, 41(Suppl.1): 393–400. Wang P Y, Li Z Q, Wu L H, et al. Application of ground penetrating radar to the survey of glacier thickness and bedrock topography. Journal of Jilin University (Earth Science Edition) (in Chinese) (in Chinese) , 2011, 41(Suppl.1): 393-400.
[6] 武震, 张世强, 刘时银, 等. 祁连山老虎沟12号冰川冰内结构特征分析. 地球科学进展 , 2011, 26(6): 631–641. Wu Z, Zhang S Q, Liu S Y, et al. Structural characteristics of the No. 12 glacier in Laohugou Valley, Qilian Mountain based on the ground penetrating radar combined with FDTD simulation. Advances in Earth Science (in Chinese) (in Chinese) , 2011, 26(6): 631-641.
[7] Plewes L A, Hubbard B. A review of the use of radio-echo sounding in glaciology. Progress in Physical Geography , 2001, 25(2): 203-236. DOI:10.1191/030913301668581943
[8] 杨天春, 吕绍林, 王齐仁. 探地雷达检测道路厚度结构的应用现状及进展. 物探与化探 , 2003, 27(1): 79–82. Yang T C, Lü S L, Wang Q R. The application and development of GPR in detection of pavement thickness and highway structures. Geophysical and Geochemical Exploration (in Chinese) (in Chinese) , 2003, 27(1): 79-82.
[9] 卢成明, 秦臻, 朱海龙, 等. 探地雷达检测公路结构层隐含裂缝实用方法研究. 地球物理学报 , 2007, 50(5): 1558–1568. Lu C M, Qin Z, Zhu H L, et al. Practical methods for detection of concealed cracks in highway pavement using ground penetration radar data. Chinese J. Geophys. (in Chinese) , 2007, 50(5): 1558-1568.
[10] 冯德山, 戴前伟. 高速公路路面厚度探地雷达检测. 地球物理学进展 , 2008, 23(1): 289–294. Feng D S, Dai Q W. Application of ground penetrating radar in the survey of the pavement thickness in highway. Progress in Geophysics (in Chinese) (in Chinese) , 2008, 23(1): 289-294.
[11] 刘四新, 曾昭发. 频散介质中地质雷达波传播的数值模拟. 地球物理学报 , 2007, 50(1): 320–326. Liu S X, Zeng Z F. Numerical simulation for Ground Penetrating Radar wave propagation in the dispersive medium. Chinese J. Geophys. (in Chinese) , 2007, 50(1): 320-326.
[12] 周立军, 董荣伟, 徐波, 等. 探地雷达在城市地质调查中的应用. 工程地球物理学报 , 2009, 6(5): 632–635. Zhou L J, Dong R W, Xu B, et al. The application of ground penetrating radar in city geological investigation. Chinese Journal of Engineering Geophysics (in Chinese) (in Chinese) , 2009, 6(5): 632-635.
[13] Widess M B. How thin is a thin bed. Geophysics , 1973, 38(6): 1176-1180. DOI:10.1190/1.1440403
[14] Huang N E, Shen Z, Long S R, et al. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis. Proc. R. Soc. London , 1998, 454(1971): 903-995. DOI:10.1098/rspa.1998.0193
[15] Huang N E, Shen Z, Long S R. A new view of nonlinear water waves-Hilbert Spectrum. Ann. Rev. Fluid Mech. , 1999, 31: 417-457. DOI:10.1146/annurev.fluid.31.1.417
[16] 汤井田, 化希瑞, 曹哲民, 等. Hilbert-Huang变换与大地电磁噪声压制. 地球物理学报 , 2008, 51(2): 603–610. Tang J T, Hua X R, Cao Z M, et al. Hilbert-Huang transformation and noise suppression of magnetotelluric sounding data. Chinese J. Geophys. (in Chinese) , 2008, 51(2): 603-610.
[17] 孔金欧.电磁波理论.吴季译.北京:电子工业出版社, 2003: 210-220. Kong J O. Electromagnetic Wave Theory (in Chinese). Wu J Trans. Beijing: Electronic Industry Press, 2003: 210-220.
[18] 何樵登, 熊维纲. 应用地球物理教程--地震勘探. 北京: 地质出版社, 1991 : 27 -30. He Q D, Xiong W G. The Applied Geophysics Course-Seismic Exploration (in Chinese). Beijing: Geological Publishing House, 1991 : 27 -30.
[19] Gabor D. Theory of communication, Part Ⅰ, Ⅱ, Ⅲ. Journal of the IEEE , 1946, 93(26): 429-457.
[20] Gressman A, Morlet J. Decomposition of hardy functions into square integrable wavelets of constant shape. SIAM J. Math. Anal. , 1984, 15(4): 723-736. DOI:10.1137/0515056
[21] 高静怀, 汪文秉, 朱光明. 小波变换与信号瞬时特征分析. 地球物理学报 , 1997, 40(6): 821–832. Gao J H, Wang W B, Zhu G M. Wavelet transform and instantaneous attributes analysis. Chinese J. Geophys. (in Chinese) , 1997, 40(6): 821-832.
[22] Stockwell R G, Mansinha L, Lowe R P. Localization of the complex spectrum: the S transform. IEEE Transactions on Signal Processing , 1996, 44(4): 998-1001. DOI:10.1109/78.492555
[23] 高静怀, 陈文超, 李幼铭, 等. 广义S变换与薄互层地震响应分析. 地球物理学报 , 2003, 46(4): 526–532. Gao J H, Chen W C, Li Y M, et al. Generalized S transform and seismic response analysis of thin interbeds. Chinese J. Geophys. (in Chinese) , 2003, 46(4): 526-532.
[24] 刘喜武, 刘洪, 李幼铭, 等. 基于广义S变换研究地震地层特征. 地球物理学进展 , 2006, 21(2): 440–451. Liu X W, Liu H, Li Y M, et al. Study on characteristics of seismic stratigraphy by generalized S-transform. Progress in Geophysics (in Chinese) (in Chinese) , 2006, 21(2): 440-451.