DOI: 10.11817/j.issn.1672-7207.2015.08.027
基于多模式分离的S变换瑞雷波频散曲线提取
张大洲1, 2，练小聪1，杨威1，熊章强1, 2
(1. 中南大学 地球科学与信息物理学院，湖南 长沙，410083；
2. 中南大学 有色金属成矿预测教育部重点实验室，湖南 长沙，410083)
摘要：对多阶模式下的瑞雷波信号进行分离，利用S变换计算2道数据单模式的频散曲线并进行反演得到二维横波速度剖面，以提高瞬态瑞雷波法横向分辨率。对理论模型和实际瑞雷波数据进行处理和分析，研究结果表明：利用S变换计算2道数据的频散曲线时，采用F-K分离后的单模式瑞雷波数据可有效提高计算精度，在此基础上对多条频散曲线进行反演得到的二维横波速度剖面可提高瑞雷波法勘探的横向分辨率；采用该方法还可计算2道数据的高阶模式频散曲线。
关键词：横向分辨率；频散曲线；F-K分离；S变换
中图分类号：P631 文献标志码：A 文章编号：1672-7207(2015)08-2950-07
Extract dispersion curves of Rayleigh wave through S transform based on F-K separation
ZHANG Dazhou1, 2, LIAN Xiaocong1, YANG Wei1, XIONG Zhangqiang1, 2
(1. School of Geosciences and Info-Physics, Central South University, Changsha 410083, China;
2. Key Laboratory of Metallogenic Prediction of Nonferrous Metals, Ministry of Education,
Central South University, Changsha 410083, China)
Abstract: The accuracy of transient Rayleigh wave lateral resolution was improved according to the two-dimensional S-wave velocity profile which was gotten by inversing multiple dispersion curves. Those dispersion curves were calculated by S transform using two channel single mode data which separate the multi-mode Rayleigh wave signal. The data from the theoretical model and the actual Rayleigh wave were processed and analyzed. The results show that the precision of dispersion curves calculated by S transforms method can be improved if the single mode Rayleigh wave data separated by F-K method are used. Based on this, S-wave velocity profile can be obtained from inversing those dispersion curves and the lateral resolution precision of Rayleigh wave prospecting can be improved. The high mode Rayleigh wave dispersion curves of two channel data can be obtained by using this method.
Key words: lateral resolution precision; dispersion curve; F-K separation; S transform
利用瑞雷波的频散特性可直接对场地进行分层，并可计算出各层的横波速度[1-3]。但目前瑞雷波勘探中的数据处理仍以多道面波为主，而多道面波相对于2道面波来说，虽然其频散曲线的精度较高，但横向分辨率较低。为了提高瑞雷波法的横向分辨率，可采用2道瑞雷波分析方法，但由于提取的频散曲线精度较低，在瑞雷波勘探中的效果较差。如何提高2道瑞雷波频散曲线的计算精度，进而提高瑞雷波勘探的精度是瞬态瑞雷波法应用中一个非常重要的问题。计算2道瑞雷波数据的频散曲线目前常采用的是相位差法[4]，这种方法利用傅里叶变换得到相邻2道瑞雷波信号的频率谱和相位谱，再计算相同频率下的相位差，从而得到该频率下的相速度。但该方法受到傅里叶变换本身的限制，其相位变化范围为[-π，+π]，这就要求在数据采集时，2个检波器之间的距离不得大于接收波的最小波长，否则，对于相位差大于2π的信号，该方法提取的相位差是不真实的。为了提高2道瑞雷波数据频散曲线的计算精度，采用时频分析方法。这种方法是通过对2道瑞雷数据进行时频变换，将时间域的信号转化到时频域，再利用2道信号之间相同频率能量极值的时间差计算相速度，从而克服相位差法对于道间距的限制[5]。目前所采用的时频分析方法主要有短时傅里叶变换[6]、小波变换[7]和S变换[8]。其中短时傅里叶变换只能有单一的分辨率，对于同时具有高、低频的地震信号不太适用；小波变换具有多分辨率，但由于尺度因子与频率无直接联系，且小波选取的多可能性使得结果具多解性，因此，计算相速度时其精度也不高。S变换是一种线性、自适应的时频分析方法，它具有小波变换的多分辨率，同时又如傅里变换一样具有无损可逆性，并且S变换对于多分量信号具有可加性，是一种线性变换，因此，使用S变换计算频散曲线时可以得到更好的效果。周竹生等[9]利用含可变因子的广义S变换提取瑞雷面波频散曲线，但该研究没有涉及瑞雷波的多模态特征。在进行频散曲线计算时，由于原始地震数据中包含各种不同模态的瑞雷波[10]，而不同模态瑞雷波的能量极值点出现时间不同，这样就会影响相速度的计算精度。另外，虽然S变换计算频散曲线效果较好，但它无法直接计算高阶模式[11]的频散曲线。为此，本文作者通过F-K分离法[12]，先将不同模态的瑞雷波数据进行分离，避免各阶模态对极值点的影响，再使用S变换准确提取基阶和高阶瑞雷波频散曲线，这样就可以得到较准确的2道瑞雷波的频散曲线，从而提高瑞雷波的勘探精度。
1 方法原理
基于多模式分离的S变换瑞雷波频散曲线的提取原理主要分为瑞雷波F-K域的分离和S变换计算频散曲线2部分。瑞雷波F-K域的分离原理见文献[13]。下面就S变换计算频散曲线的原理进行介绍。S变 换是Stockwell等[14]提出的一种时频分析方法，其正变换表达式为
(1)
其离散形式的表达式为
(2)
式中：u为h(t)的S变换；f为频率；τ和t为时间；τ可控制时间轴上高斯窗的位置。其反变换为
(3)
需说明的是：在对瑞雷波数据进行时频分析时，所计算的幅值就是该时刻瑞雷波能量，而根据2道最大幅值之差计算的速度为瑞雷波传播的群速度，而不是相速度。但在瑞雷波勘探中必须获得瑞雷波传播的相速度，因此，需要将时频分析中计算得到的群速度转换为相速度，其转换关系可表示为[15]
vg=dω/dk，vR=ω/k (4)
其中：vg为群速度；vR为相速度；ω为角频率；k为波数。联立式(3)和(4)，得到vR关于vg和ω的关系式如下：
(5)
S变换提取瑞雷波频散曲线的步骤如下：第1步，对原始瑞雷波数据利用F-K法进行多模式分离；第2步，利用分离后的数据分别对相距为d的2道信号进行S变换，得到2个时频谱；第3步，固定频率f，分别在2个时频谱上寻找该频率能量极值点所对应的时间t1和t2，并计算出时差△t，△t=|t2-t1|；第4步，计算群速度，vg(f)=d/△t；第5步，利用式(5)计算相速度vR(f)；第6步，对每一个频率重复第3~5步就可得到1条频散曲线。
综上所述，将采集到的瑞雷波记录先进行F-K分离提取单一模态的记录，然后用S变换就可计算得到较高精度的频散曲线。
2 理论模型计算
2.1 水平2层介质模型
为了验证以上方法提取频散曲线的效果，特设计1个2层介质模型。模型长×宽为50 m×50 m。模型参数如下：第1层，厚度为5 m，纵波速度vp=500 m/s，模波波数vs=200 m/s，密度ρ=1.9 g/cm3；第2层，vp=1 000 m/s，vs=400 m/s，ρ=2.0 g/cm3。利用高精度交错网格有限差分法对该模型进行全波场正演模拟[16]，模拟采用的参数如下：道间距为1.0 m，采样间隔为0.2 ms，采样长度为0.6 s，震源是频率为20 Hz的雷克子波。所得到的地震记录如图1(a)所示。
对模拟得到的地震记录进行F-K变换，结果如图1(b)所示。从图1(b)可以看到：频率-波数域中除了基阶模式外还有其他高阶模式能量团，但各模式的瑞雷波在F-K域中是相互分开的，这为各模式瑞雷波的分离提供了基础。为了获得基阶模式瑞雷波，在F-K域中剔除高阶模式能量团，保留基阶模式的能量团，然后进行反变换，得到分离后的基阶瑞雷波记录，如图2所示。
图1 模式分离前的瑞雷波记录和频率波数谱
Fig. 1 Synthetic data and frequency wave-number spectrum before separated
图2 基阶模式瑞雷波记录
Fig. 2 Record of fundamental model
在完成瑞雷波各模态分离后就可采用S变换来计算频散曲线。若采用未分离的地震记录进行S变换，则由于不同模态间极大值出现的时间不同，导致其极大值出现偏差。任意取1道模拟的地震记录(现取第5道)，分别对基阶、一阶和混合模式进行S变换，并计算30 Hz对应的振幅，如图3所示。由图3可见：混合模式的振幅谱由于受基阶和一阶模式振幅的共同影响，极大值出现的时间为0.061 6 s，位于基阶模式峰值和一阶模式峰值之间。若采用混合模式的最大值所对应的时间计算瑞雷波相速度则会产生误差，因此，利用多模式分离后的数据进行S变化进而计算频散曲线是非常必要的。
选取图2中基阶模式瑞雷波记录第5和第6道进行S变换，得到其时频谱如图4所示。固定频率为30 Hz，在第5道和第6道振幅值最大值分别对应时间为0.059 3 s和0.065 2 s，如图5所示，计算时间差为△t=0.0059 s。
图3 第5道30 Hz下不同模式的瑞雷波振幅对比图
Fig. 3 Amplitude of different models for the fifth trace in 30 Hz
用道间距d除以△t即为30 Hz频率所对应的群速度，可得vg=169.5 m/s，转换则可得相速度vR=197.2 m/s。与30 Hz的理论相速度195.8 m/s相比相对误差仅为0.72%。同理，若为混合模式瑞雷波，则计算所得时间差为0.006 8 s，群速度为vg=147.05 m/s，则对应的相速度vR=205.5 m/s，相对误差为4.95%。
同理可计算出其他频率的相速度，最后得到1条频率-相速度曲线，并将它与理论值对比，如图6所示。
图4 30 Hz时的基阶瑞雷波S变换时频谱
Fig. 4 Time-frequency spectrum of fundamental-model By S transform in 30 Hz
图5 第5道和第6道30 Hz时的基阶瑞雷波振幅值对比图
Fig. 5 Comparison of fundamental model’s amplitude spum of the fifth trace and the sixth trace in 30 Hz
从图6可以看出：进行模式分离后得到的基阶频散曲线与理论值基本重合，相对误差较大的低频段最大误差也仅为2%。低频段误差较大的原因是为了获得完整的低频信号，使用宽度较大的高斯窗时降低了时间分辨率。
采用相同的方法对未分离的瑞雷波数据计算其频散曲线，并与理论值及分离后的基阶瑞雷波频散曲线进行对比，如图7所示。从图7可以看出：未进行模式分离的瑞雷波其相速度在频率为20~30 Hz段误差较大，最大达到5.4%，由此可知瑞雷波各模态的分离对于提高频散曲线精度是非常必要的。
对瑞雷波原始数据进行模态分离，一方面可提高计算频散曲线的精度，另一方面可通过分离后的高模式瑞雷波数据计算出高模式瑞雷波频散曲线。图8所示为计算得到的一阶模式频散曲线与理论值的对比结果，可见2条频散曲线基本吻合。同理可计算其他高阶模式的频散曲线。因此，利用此方法计算2道间不同模式的频散并进行多模式联合反演，可提高瑞雷波的反演精度。
图6 基阶瑞雷波群速度、相速度及理论值对比
Fig. 6 Comparison of velocity curves of fundamental model surface wave among group-velocity, phase-velocity and theory values
图7 分离前后瑞雷波相速度频散曲线与理论值对比
Fig. 7 Comparison of dispersion curves between phase-velocity before and after separation and theory values
图8 一阶瑞雷波相速度频散曲线与理论值对比
Fig. 8 Comparison of dispersion curves of fundamental model surface wave between calculated results and theory values
2.2 倾斜2层介质模型
为了进一步说明S变换法在实际应用中的优点及瑞雷波分离后计算频散曲线的效果，设计1个倾斜界面的2层介质模型，如图9所示。模型长×宽为20 m×20 m。模型参数如下：第1层，厚度5 m，纵波波速vp=800 m/s，横波波速vs=200 m/s，密度ρ= 1.9 g/cm3；第2层，vp=1 000 m/s，vs=400 m/s，ρ=2.0g/cm3。利用高精度交错网格有限差分法对该模型进行正演模拟，采样间隔为0.2 ms，采样时长为0.3 s，震源采用20 Hz的雷克子波。图10所示为模拟得到的地震合成记录。若对模拟得到的合成记录采用多道瑞雷波分析方法，则对整个剖面可得到如图11所示的1条频散曲线，最终也只能得到界面的平均深度，无法获得界面的倾斜特性。若对分离后的地震记录利用S变换法，则只需要其中的2道数据就可计算得到1条频散曲线。这样通过对多道地震记录可以得到多条频散曲线，从而提高了横向分辨率。对图10所示的合成记录先进行S变换，然后计算相速度，按半波长理论，可得到相速度-半波长剖面图，如图12所示。从图12可以看出界面的倾斜特征较明显，这也充分体现了瑞雷波的横向分辨率。
图9 倾斜界面模型
Fig. 9 Model of slant interface
图10 倾斜界面模型合成地震记录(垂直分量)
Fig. 10 Synthetic data (vertical-component) of slant interface model
图11 相移法计算得到的频散曲线
Fig. 11 Dispersion curves calculated by phase-shift method
图12 利用S变换计算得到的相速度剖面图
Fig. 12 Phase velocity profile calculated by S transform method
3 实际算例
为了分析S变换提取实际瑞雷波数据频散曲线的应用效果，给出1个实例进行说明。某场地地层自上而下分为粉质黏土、强风化泥质砂岩及中风化泥质砂岩。3个地质层位存在着密度和速度差异，这为瑞雷波勘探提供了良好的地球物理条件。瑞雷波野外数据采集参数设置如下：偏移距为5 m，采样间隔为0.5 ms，采样长度为0.6 s，道间距为1 m，共48道。图13所示为原始瑞雷波记录，图14所示为利用F-K变换方法分离出的基阶瑞雷波记录。对该记录通过S变换计算每2道间的频散曲线，然后用最小二乘法进行反演，最后得到如图15所示的横波速度断面图。分析图15可知：随着深度的增加横波速度逐渐增大，0~6.0 m处为粉质黏土层，横波速度为150~200 m/s；6.0~14.0 m处为强风化泥岩层，横波速度为200~400 m/s，强风化层厚度在横向上变化较大；14.0~25.0 m处为中风化泥岩层，横波速度为400~800 m/s。
通过以上实例可知：根据不同岩土层的横波速度差异，利用单个排列的数据就可得到二维横波速度剖面，结合钻孔资料就可以解释测线下方的二维地质结构。该方法相对于一维多道瑞雷波法(该方法对1个排列仅得到1条频散曲线)，大大提高了瑞雷波勘探的横向分辨率。
图13 实际瑞雷波记录
Fig. 13 Real seismic record
图14 进行F-K分离后的瑞雷波记录
Fig. 14 Record of Rayleigh wave after F-K separated
图15 横波速度剖面图
Fig. 15 Inversion contour map of shear wave velocity
4 结论
1) 利用2道地震记录计算频散曲线时，通过将F-K分离后的瑞雷波进行S变换从而计算得到频散曲线，可有效抑制时差计算时瑞雷波多模式之间的相互干扰，从而提高频散曲线的计算精度；同时，利用S变换可以得到高模式瑞雷波频散曲线。
2) 对单炮瑞雷波记录中每2道计算其频散曲线，进而反演得到横波速度分布，获得二维地质结构，从而提高瑞雷波勘探的横向分辨率。
参考文献：
[1] 曹圣华, 杨晓东, 朱正坤. 瑞雷面波法用于软土地基勘察的试验研究[J]. 河海大学学报(自然科学版), 2003, 31(4): 419-423.
CAO Shenghua, YANG Xiaodong, ZHU Zhengkun. Experimental study on reconnaissance of soft foundations by Reyleigh surface wave method[J]. Journal of Hehai University (Natural Sciences), 2003, 31(4): 419-423.
[2] 童立元, 陈征宙, 刘定恩. 瞬态瑞雷面波勘探技术在公路工程勘察中的应用[J]. 地质与勘探, 2000, 36(5): 73-75.
DONG Liyuan, CHEN Zhengzhou, LIU Dingen. The application of transient rayleigh surface wave prospecting technique in engineering exploration of highway[J]. Geology and Prospecting, 2000, 36(5): 73-75.
[3] McMechan G A, Yedlin M J. Analysis of dispersive waves by wave-field[J]. Geophysics, 1981, 46(6): 101-113.
[4] 张大洲, 顾汉明, 熊章强, 等. 基于多模态分离的面波谱分析方法[J]. 中国地质大学学报, 2009, 34(6): 1012-1018.
ZHANG Dazhou, GU Hanming, XIONG Zhangxiang, et al. Spectrum analysis of surface wave method based on multi-mode separation[J]. Journal of China University of Geosciences, 2009, 34(6): 1012-1018.
[5] 彭文, 王亮. 瑞雷面波频散特征的时频分析方法及应用[J]. 物化探计算技术, 2006, 28(3): 233-239.
PENG Wen, WANG Liang. Theory and application of time-frequency analysis to rayleigh waves prospecting[J]. Computing Techniques for Geophysical and Geochemical Exploration, 2006, 28(3): 233-239.
[6] Garbo D. Theory of communications[J]. Journal of Institute for Electrical Engineering, 1946, 93(3): 429-457.
[7] Gressmann A, Morlet J. Decomposition of hardy functions into square integrable wavelets of constant shape[J]. SIAM Journal on Mathematical Analysis, 1984, 15(4): 723-736.
[8] Pinnegar C R, Mansinha L. The S-transform with windows of arbitrary and varying shape[J]. Geophysics, 2003, 68(1): 381-385.
[9] 周竹生, 杨明凯, 陈友良. 基于含可变因子广义S变换的瑞雷面波频散曲线提取方法[J]. 物探化探计算技术, 2012, 34(5): 518-522.
ZHOU Zhusheng, YANG Pengkai, CHEN Youliang. The method of Rayleigh wave dispersion curve extraction based on generalized S-transform with variable-factors[J]. Computing Techniques for Geophysical and Geochemical Exploration, 2012, 34(5): 518-522.
[10] 肖博勋. 高模式瑞雷面波及其正反演研究[D]. 长沙: 中南大学信息物理工程学院, 2000: 50-65.
XIAO Boxun. High Models rayleigh wave and forward modelling and inversion of its dispersion curves[D]. Changsha: School of Info-physics and Geomatics Engineering. Central South University, 2000: 50-65.
[11] Gabriels P, Snieder R, NoletG.Insitu measurements of shear-wave velocity in sediments with higher-mode Rayleigh waves[J]. Geophysical Prospecting, 1987, 35（2）: 187-196.
[12] Capon J. High-resolution frequency-wavenumber spectrum analysis[J]. Proceedings of the IEEE, 1969, 57(8): 1408-1418.
[13] 张大洲, 熊章强, 秦臻. 基于傅立叶变换的瑞雷面波分离提取及实例分析[J]. 中南大学学报(自然科学版), 2010, 41(2): 643-648.
ZHANG Dazhou, XION Zhangqiang, QIN Zhen. Separation and extraction of Rayleigh wave based on Fourier transform and case analysis[J]. Journal of Central South University (Science and Technology), 2010, 41(2): 643-648.
[14] Stockwell R G, Mansinha L, Lowe R P. Localization of the complex spectrum:the S transform[J]. IEEE Transactions on Signal Processing, 1996, 44(4): 998-1001.
[15] YAO Huajian, van der Hilst R D, de Hoop M V. Surface-wave array tomography in SE Tibet from ambient seismic noise and two-station analysis: I. Phase velocity maps[J]. Geophysical Journal International, 2006, 166(2): 732-744.
[16] 熊章强, 张大洲, 秦臻. 瑞雷波数值模拟中的边界条件及模拟实例分析[J]. 中南大学学报(自然科学版), 2008, 39(4): 824-830.
XIONG Zhangqiang, ZHANG Dazhou, QIN Zhen. Boundary conditions and case analysis of numerical modeling of Rayleigh wave[J]. Journal of Central South University (Science and Technology), 2008, 39(4): 824-830.
(编辑 陈灿华)
收稿日期：2014-09-10；修回日期：2014-11-20
基金项目(Foundation item)：国家自然科学基金资助项目(41274123)；教育部博士点基金资助项目(20130162110066)(Project (41274123) supported by the National Natural Science Foundation of China；Project (20130162110066) supported by Doctoral Fund of Ministry of Education)
通信作者：张大洲，博士，讲师，从事地震波正反演及偏移成像研究；E-mail：dazhou2005@163.com
