利用EEMD及超导重力数据探测超低频地球自由振荡模态分裂信号：探测实例

曾诗瑜　申文斌,2

1　武汉大学测绘学院，武汉市珞喻路129号，430079　2　武汉大学测绘遥感信息工程国家重点实验室，武汉市珞喻路129号，430079



利用EEMD及超导重力数据探测超低频地球自由振荡模态分裂信号：探测实例

曾诗瑜1申文斌1,2

1武汉大学测绘学院，武汉市珞喻路129号，4300792武汉大学测绘遥感信息工程国家重点实验室，武汉市珞喻路129号，430079

摘要：利用2004年苏门答腊9.0级地震后GGP台网中15台超导重力仪的18组观测数据，采用EEMD分解技术及叠积技术，检测了频率小于1 mHz的低频地球自由振荡的简正模(3S1、0S4、0S5)及其频谱分裂现象。通过对原始分钟间隔的超导重力数据进行预处理得到重力残差，然后运用EEMD方法将其分解到处于不同频段上的IMF分量，提高了探测低频自由振荡信号的信噪比，使一些低频自由振荡信号的谱线分裂更清晰地显现出来。结果表明，经过EEMD分解过的重力残差能更有效地检测低频自由振荡信号，得到信噪比更高的低阶球型振荡分裂谱线。该研究进一步证实了EEMD数据处理方法的有效性以及超导重力仪在低阶地球自由振荡检测中的优越性。

关键词：地球背景噪声；低频地球自由振荡；频率；振幅；品质因子

地球上发生大地震时，不仅可以产生涉及地球局部运动的体波和面波，还能激发全球性的自由振荡。通过研究地球自由振荡，可以揭示地球的内部结构、动力学问题以及震源特性[1]。地球自由振荡的频率与地球的形状、密度分布、剪切模量和体积模量等密切相关[2]。地球自转和椭率的综合作用使地球自由振荡的本征频率不是简正的，从而产生相应的谱线(或谱峰)分裂耦合现象，同时地幔的横向不均匀性和内核的各向异性会导致自由振荡谱的异常分裂[3]。此外，影响地球自由振荡谱形态的因素还包括地球物理学结构横向和径向的不均匀性、地球内部结构的各向异性、非完全弹性以及震源机制性质等[2,4]。地球自由振荡频谱的振幅与震源机制有关[4]，我们将观测到的自由振荡的周期和振幅与PREM理论模型的振荡周期和振幅相比较，就可以对地球内部结构进行不同于体波和面波的另一种约束[2,4]。由于地球长周期自由振荡的振幅主要依赖于地震断层的破裂方式及地震矩的大小，因此利用观测记录得到的地球自由振荡波形资料，可以反演地震的震源破裂参数，从而约束地震的震源机制、地震大小及持续时间[2]，还可以对地球Q值进行约束[3]。例如，0S0(膨胀收缩模)的振幅大小与地震断层的倾角、滑动方向角、震源深度及地震断层的破裂时间有关[3]。

超导重力仪(SG)对地球自由振荡本征周期的检测能力可与弹簧重力仪和宽频地震仪媲美，尤其是利用单台SG观测记录就可以对低频自由振荡简正模式的频谱分裂现象进行高分辨率分析[5-6]，从而得到相应的本征频率和自由振荡参数。大量研究表明，SG还在构制频率小于1 mHz的长周期地震图和研究地球深部构造方面发挥着重要作用。

本文主要是通过对来自GGP台网(http:∥www.eas.slu.edu/GGP/ggphome.html)2004年苏门答腊地震的SG观测数据进行预处理，从其残差中提取出低频率(0～1 mHz)自由振荡模态，进而对地球深部结构作出相应的解释。

1　提取低频自由振荡信号的方法

1.1EEMD分解

本文运用EEMD方法[7]对每个台站经过气压改正和潮汐改正后的最终重力残差时间序列进行分解(这里的EEMD方法起到了一个滤波器的作用)，得到一组有限的IMF分量(见图1，数据长度为300 h，分解共得到9个IMF和1个残余趋势项,纵轴单位nm/s2)。由于不同的IMF处于不同的频率域，因此不同范围的重力残差序列的信号能由合适的IMF分量重新组合。另外，EEMD能被用来调频[7]，可用于检波调频时间序列。因此，EEMD可以同时用作滤波器和检波器。

图1　Canberra台重力残差在2004年苏门答腊地震5 h后的EEMD分解结果Fig.1　Results obtained after EEMD was applied to above Canberra station’s residual gravity sequences starting from 2004 Sumatra earthquake 5 h later

图2　对Canberra台站经过EEMD分解后的IMF1～IMF6进行功率谱密度估计得到的功率谱结果Fig.2　Power spectra of the IMF1-IMF6 of the SG residual record obtained at the Canberra station after using EEMD

1.2谱估计

本文采用Chao等[8]基于Prony方法提出的频率域自回归估计方法对自由振荡的相关参数进行估计。为了验证信号经过EEMD分解后被分解到相应的频段上，我们对重力残差的一组IMF进行功率谱密度估计，并将它转化到频率域。由于IMF7～IMF9所处的频率几乎为0，而地球自由振荡模态所处的频带主要在0.2～1 mHz，所以只需分析IMF1～IMF6。如图2(IMF7～IMF10没有画出，箭头表示相关的自由振荡信号模态,纵轴单位nm2/s4)所示，0S0出现在IMF2～IMF4中，0S2仅出现在IMF5中，0S3出现在IMF4中，0S5和3S1出现在IMF2～IMF3中。在IMF5中可以观测到0S2的5重分裂谱峰，因此IMF5可作为观测0S2谱线分裂的新序列。此外，观测3S1的谱线分裂时，IMF2和IMF3可以作为一个新序列，与Shen等[9]的结论一致。结果表明，信号的IMF严格按照分解顺序从先到后频率依次降低，而且没有出现模态混叠的现象，可以达到信号分解的目的。

2　低频地球自由振荡模态及其频谱分裂现象的检测

选取GGP观测台网2004年苏门答腊地震中的15个超导重力仪台站，主要通过分钟间隔的SG数据来求解0～1 mHz频段内的地球低频自由振荡模态，从而为评估其低频简正模和频谱分裂现象提供高精度的观测数据。1982年Dahlen[10]提出理论上最佳的数据记录长度约为1.1Q个周期，以增大信号的SNR。本文对于不同的地球自由振荡模态，选取数据的起始时间和长度需视情况而定。

图3　由叠积方法得到18个SG记录的模态3S1归一化功率谱密度(PSD)Fig.3　The normalized product spectra density ( PSD ) of the mode 3S1 from 18 SG records with stacking method

2.1检测3S1的谱线分裂

理论上在3S1存在的频段，多线态1S3对3S1的观测会产生交叉耦合效应。前人的作法通常是采用震后一定时间(例如2 d)的数据，以尽量削弱1S3的强度[11]。本文对每个SG记录取地震开始50 h之后，数据总长度为650 h。首先对18个SG记录(h1、h2、cb、es、ka、ma、mb、mc、me、m1、ny、st、s1、s2、w1、w2、vi、tc)作叠积，得到3S1的归一化功率谱密度(见图3，垂直蓝色虚线对应于PREM模型中相关谱峰的值，曲线平滑采用了3次样条插值。图4同)。随后选取4个SG记录(h1、h2、m1、 st)进行EEMD分解来探测3S3的谱线分裂，如图4(a)～(d) 所示。EEMD可以使谱峰m=0被探测到，且谱峰m=-1的相对振幅会有一定提高，说明EEMD可以将信号分解到不同频段上，且提高相对信噪比。采用叠积方法得到的各谱峰的频率估计值与PREM模型的理论预测值的相对误差在0.03%左右，m=±1的谱峰估计值与Roult等[12]的结果接近；采用EEMD方法得到的各谱峰的估计值与PREM模型的理论预测值的相对误差在0.04%以内，并且与Ding等[13]的结果接近，m=0的谱峰估计值与Roult等[12]和Shen等[14]的结果接近(表1，单位mHz)。

图4　在0.93～0.96 mHz频段内4个SG数据获得的模态3S1的归一化振幅谱Fig.4　The normalized amplitude of the mode 3S1obtained from 4 SG records during 0.93-0.96 mHz

2.2检测0S4的谱线分裂

观测0S4的谱线分裂，同时要考虑到信噪比和频率精度，本文选取各记录开始于地震5 h后。首先对18个SG记录(h1、h2、cb、es、ka、ma、mb、mc、me、m1、ny、st、s1、s2、w1、w2、vi、tc)作叠积，数据长度为555 h，相应的归一化功率谱密度如图5

表1　3S1分裂谱线的模型预测值与频率观测值

注：上标a为通过叠积方法得到相应的功率谱密度(PSD)，用Lorentz谐振函数与观测到的谱峰进行匹配计算求得；上标b为通过EEMD分解以及Fourier变换来观测模态的谱峰，谱峰的频率通过AR方法来估计，最后对台站的结果进行加权平均得到最终的频率估计值。表2、表3同。

(垂直蓝色虚线对应于PREM模型中相关谱峰的值，曲线平滑采用了3次样条插值)所示。由于耦合效应，很难从频谱图中剥离出各个谱峰。对4个SG记录(h1、h2、cb、st)进行EEMD分解以探测0S4的谱线分裂，记录开始于地震5 h后，数据长度约为706.68 h。如图6 (a)～(d) 所示，在使用EEMD分解之后，0S4谱峰的相对信噪比有一定程度的提高(图6 (a)、(d) )，0S4的9个谱峰可以区分开来，这进一步证实了EEMD的可行性。从表2可以看出，叠积结果的频率估计值与PREM理论模型预测值的相对误差除了m=-4为-0.116%外，其他谱峰的相对误差均在0.02%左右；使用EEMD方法得到的各谱峰的频率估计值与PREM理论模型的预测值的相对误差除了m=-4谱峰为-0.103%外，其他谱峰的相对误差均在0.05%左右，其估计值与Ding等[13]的结果较为接近。尽管估计的频率序列不同于相关PREM模型的理论预测值以及Roult等[12]研究的观测值(在他们的研究中，只有m=±4的谱峰是2004年Sumatra地震后在超过50台宽频地震仪的基础上观测得到的)，但是估计的频率相互之间有很强的关联性。可以看出，运用EEMD分解技术可以清楚地观测到所有的分裂谱线。

图5　由叠积方法得到18个SG记录的模态0S4的归一化功率谱密度(PSD)Fig.5　The normalized product spectra density ( PSD ) of the mode 0S4 from 18 SG records with stacking method

图6　在0.635～0.660 mHz频段内由4个SG数据获得的模态0S4的归一化振幅谱Fig.6　The normalized amplitude of the mode 0S4obtained from 4 SG records during 0.635～0.660 mHz

模态m=-4/+4m=-3/+3m=-2/+2m=-1/+1m=0PREM0.642098/0.6516050.643401/0.6505310.644671/0.6494240.645908/0.6482850.647113Roult等(2010)0.64147±2.562×10-4/0.65106±1.908×10-4Ding等(2013)(GJI)0.641642±4.5×10-5/0.651738±4.0×10-50.643336±2.0×10-5/0.650723±1.7×10-50.644585±1.7×10-5/0.649551±1.3×10-50.645858±3.8×10-5/0.648325±9.4×10-60.647057±3.1×10-5本文a0.641352±1.1×10-60.643284±1.3×10-6/0.649692±9.6×10-60.644401±1.2×10-6/0.649406±1.3×10-60.646195±5.1×10-60.647323±1.3×10-6本文b0.641438±9.8×10-7/0.651594±1.1×10-60.643370±1.2×10-6/0.650923±4.9×10-70.644464±8.7×10-7/0.649357±9.7×10-70.646374±6.9×10-7/0.648323±9.6×10-70.647124±7.0×10-7

2.3检测0S5的谱线分裂

观测0S5的谱线分裂，首先对18个SG记录(h1、h2、cb、es、ka、ma、mb、mc、me、m1、ny、st、s1、s2、w1、w2、vi、tc)作叠积，选取各记录开始于地震5 h后，数据长度为698 h，得到0S5的归一化功率谱密度(见图7，垂直蓝色虚线对应于PREM模型中相关谱峰的值，曲线平滑采用了3次样条插值)。由于耦合重叠的作用，很难从频谱图中辨别出所有谱峰。随后选取4个SG记录(h2、m1、st、w1)，通过EEMD分解来探测0S5的谱线分裂，记录开始于地震5 h后，数据长度约为706.685 h。如图8所示，通过对比使用EEMD分解前后的结果发现，通过对信号进行EEMD分解以及目标频段的IMF叠加，可以清楚地观测到大约10个谱峰，且各谱峰的相对信噪比都有所提高。使用叠积方法得到0S5的谱峰频率估计结果与PREM模型理论预测值的相对误差在0.04%左右。由于耦合效应，谱峰m=+3与m=+5的估计值没有估计出来，使用EEMD分解后各谱峰的频率估计值与PREM理论模型预测值的相对误差也在0.04%左右，只有谱峰m=+3没有估计出来(表3)。

图7　由叠积方法得到18个SG记录的模态0S5的归一化功率谱密度(PSD)Fig.7　The normalized product spectra density ( PSD ) of the mode 0S5 from 18 SG records with stacking method

图8　在0.832～0.848 mHz频段内由4个SG数据获得的模态0S5的归一化振幅谱Fig.8　The normalized amplitude of the mode 0S5obtained from 4 SG records during 0.832～0.848 mHz

模态m=-5/+5m=-4/+4m=-3/+3m=-2/+2m=-1/+1m=0PREM0.836467/0.8435560.837398/0.8430690.838279/0.8425330.839111/0.8419470.839895/0.8413120.840628本文a0.836059±3.6×10-60.836782±1.2×10-6/0.842712±6.0×10-60.837720±3.6×10-60.838627±1.0×10-6/0.841671±7.0×10-60.839233±3.9×10-6/0.840815±4.7×10-60.840109±7.2×10-6本文b0.836147±4.7×10-7/0.843738±5.2×10-70.837189±5.4×10-7/0.843007±6.5×10-70.838167±5.6×10-70.838799±7.5×10-7/0.842163±6.0×10-70.839472±6.0×10-7/0.841295±5.3×10-70.840534±5.7×10-7

3　结　语

本文利用EEMD技术，采用较少的SG观测记录，以较高的信噪比剥离了3个低频自由振荡振型3S1、0S4、0S5几乎所有的分裂谱线。随后，对各个台站的数据进行叠积。得到的叠积谱密度图显示，叠积技术能很好地压制周围环境的噪声，提高自由振荡信号的相对信噪比，能有效观测到微弱的低频信号，但其在谱线分裂检测方面的效果比不上EEMD方法。研究表明，一个谱峰的信噪比与其相对应的估计频率的误差呈负相关关系，这与Häfner等[15]及Ding等[13]得到的结果相一致。此外，本文给出的对应于这3个模态全部分裂谱峰的频率估计与PREM模型的理论预测值有不同程度的偏离，可为构建3D地球模型提供有效约束。

尽管典型的多重态剥离法能提供谱峰的频率和振幅信息，但它极大地依赖于来源于不同测站接收器的结构和记录，有时还需要结合不同的地震事件。为了提高频率估计的精度，我们可以采用其他方法对自由振荡信号进行分离或将EEMD技术与叠积技术相结合，并改进频率估计方法，从而既能压制周围噪声的干扰，又能观测到更高精度的简正模的谱线分裂。

地球的自由振荡参数，特别是振荡信号的频率和品质因子Q值等是直接反映地球内部介质粘弹性质的重要物理量。由于目前还很难准确描述地球自由振荡的高阶谱线分裂现象，这在很大程度上限制了研究进展。因此，有必要综合利用多种数据源(包括长周期应变仪、地震仪和重力仪台网资料等)深入研究高阶环型模态和球型模态的谱线分裂及耦合问题，有效约束地幔结构的横向不均匀性和地球内核各向异性，深化对地球内部动力学现象的认识。

致谢:感谢丁浩和栾威提出有价值的建议，促使本文得以顺利完成。

参考文献

[1]Gilbert F, Dziewonski A M. An Application of Normal Mode Theory to the Retrieval of Structural Parameters and Source Mechanisms from Seismic Spectra[J]. Phil Trans R Soc London A, 1975, 278: 187-269

[2]Derr J S. Internal Structure of the Earth Inferred from Free Oscillations[J].J Geophys Res, 1969, 74 (5): 202-219

[3]薛秀秀，胡小刚，郝晓光,等. 利用重力观测约束2011日本Tohoku大地震的震源机制[J]．地球物理学报，2012, 55 (9): 3 006-3 015(Xue Xiuxiu, Hu Xiaogang, Hao Xiaoguang, et al. Constraining Focal Mechanism of the 2011 Tohoku Earthquake by Gravity Observations[J]. Chinese J Geophys, 2012, 55(9): 3 006-3 015)

[4]Dziewonski A M, Anderson D L. Preliminary Reference Earth Model (PREM)[J]. Phys Earth Planet Interi, 1981, 25: 297-356

[5]Crossley D, Hinderer J. Report of GGP Activities to Commission 3,Completing 10 Years for the Worldwide Network of Superconducting Gravimeters[C]. IAG Symposia, 2008

[6]Rosat S, Sato T, Imanishi Y, et al. High Resolution Analysis of the Gravest Seismic Normal Modes after the 2004Mw=9 Sumatra Earthquake Using Superconducting Gravimeter Data[J]. Geophys Res Lett, 2005,32(13):313-320

[7]Huang N E, Wu Z. A Review on Hilbert-Huang Transform Method and Its Applications to Geophysical Studies[J].Rev Geophys, 2008, 46(2):2008

[8]Chao B F, Gilbert F. Autoregressive Estimation of Complex Eigenfrequencies in Low Frequency Seismic Spectra[J].Geophys J R Astron Soc, 1980, 63:641-657

[9]Shen W B, Ding H. Observation of Spheroidal Normal Mode Multiplets below 1 mHz Using Ensemble Empirical Mode Decomposition[J]. Geophysical Journal International, 2014, 196: 1 631-1 642

[10]Dahlen F A. The Effect of Data Windows on the Estimation of Free Oscillations Parameters[J].Geophys J R Astron Soc, 1982，69: 537-549

[11]Masters G, Gilbert F. Attenuation in the Earth at Low Frequencies[J]. Phil Trans R Soc A, 1983, 308 (1 504): 479-522

[12]Roult G, Roch J, Clevede E. Observation of Split Modes from the 26th December 2004 Sumatra-Andaman Mega-Event[J]. Phys Earth Planet Inter, 2010, 179:45-59

[13]Ding H, Shen W B. Comparative Study of the MSE Technique and SHS Method Used for Detection ofnS1Mode Triplet[J]. Chinese Journal Geophysics, 2013,56 (10): 3 313-3 323

[14]Shen W B, Wu B. A Case Study of Detecting the Triplet of3S1Using Superconducting Gravimeter Records with an Alternative Data Preprocessing Technique[J]. Ann Geophys, 2012, 55 (2): 293-300

[15]Häfner R, Widmer-Schnidrig R. Signature of 3-D Density Structure in Spectra of the Spheroidal Free Oscillation0S2[J]. Geophys J Int, 2013, 192 (1):285-294

Foundation support:National Key Basic Research Program of China,No. 2013CB733305; National Natural Science Foundation of China,No.41128003, 41174011, 41021061, 40974015.

About the first author:ZENG Shiyu, postgraduate, majors in the free oscillation of the earth, the analysis of superconducting gravity data and its application,E-mail: syzeng@whu.edu.cn.

收稿日期：2015-07-01

第一作者简介：曾诗瑜，硕士生，主要从事地球自由振荡、超导重力数据分析及应用研究，E-mail: syzeng@whu.edu.cn。 通讯作者：申文斌，教授，博士，博士生导师，主要从事重力理论及应用、地球自转、地球自由振荡、相对论大地测量研究，E-mail: wbshen@sgg.whu.edu.cn。

DOI：10.14075/j.jgg.2016.08.010

文章编号：1671-5942(2016)08-0698-05

中图分类号：P312

文献标识码：A

Corresponding author:SHEN Wenbin, professor, PhD, PhD supervisor, majors in gravity theory and application, earth rotation, earth’s free oscillation, relativistic geodesy, E-mail: wbshen@sgg.whu.edu.cn.

Observation of Ultralow-Frequency Earth Free Oscillation Mode Multi-Splitting Based on EEMD and Global Superconducting Gravimeter Data：A Case Study

ZENGShiyu1SHENWenbin1,2

1School of Geodesy and Geomatics,Wuhan University, 129 Luoyu Road,Wuhan 430079, China2State Key Laboratory of Information Engineering in Surveying, Mapping and Remote Sensing,Wuhan University,129 Luoyu Road, Wuhan 430079, China

Abstract:In this study, we select 15 sets of 18 of superconducting gravimeter (SG) records from GGP (global geodynamics project) stations after the 2004 Sumatra Mw 9.0 earthquake and use EEMD (ensemble empirical

mode decomposition) and singlet stacking technology to detect low-frequency earth’s free oscillations mode (3S1,0S4,0S5) and their spectral splitting with frequencies less than 1 mHz. After removing the tidal and local atmospheric pressure effects from the original minute-interval SG records, we obtain a residual gravity data set. Then, EEMD is applied to this SG time-series to obtain different IMF (intrinsic mode function) on different frequencies. This will significantly reduce the possibility of mode mixing and end effect, and it could improve some low-frequency seismic signals’ SNR (signal-to-noise ratio). Therefore, EEMD could enable some splitting spectral of low-frequency free oscillation signals to be observed more clearly. Through comparisons of the normalized amplitude of the results obtained without using EEMD and after using EEMD of the residual gravity records, the experimental results show that when EEMD is applied to residual gravity records, it can be more effective for observation of the Earth’s low-frequency signals and will obtain higher resolution of low-order spherical oscillations’ singlets. This study demonstrates that EEMD is effective in data-processing and that the superconducting gravimeter is superior in detecting the low-order earth’s free oscillations.

Key words:earth’s background noise; low-frequency earth’s free oscillations; frequency; amplitude; QF (quality factor)

项目来源：国家973计划(2013CB733305)；国家自然科学基金(41128003,41174011,41021061,40974015)。