利用EEMD及超導重力數(shù)據(jù)探測超低頻地球自由振蕩模態(tài)分裂信號：探測實例

曾詩瑜　申文斌,2

1　武漢大學測繪學院，武漢市珞喻路129號，430079　2　武漢大學測繪遙感信息工程國家重點實驗室，武漢市珞喻路129號，430079

?

利用EEMD及超導重力數(shù)據(jù)探測超低頻地球自由振蕩模態(tài)分裂信號：探測實例

曾詩瑜1申文斌1,2

1武漢大學測繪學院，武漢市珞喻路129號，4300792武漢大學測繪遙感信息工程國家重點實驗室，武漢市珞喻路129號，430079

摘要：利用2004年蘇門答臘9.0級地震后GGP臺網(wǎng)中15臺超導重力儀的18組觀測數(shù)據(jù)，采用EEMD分解技術及疊積技術，檢測了頻率小于1 mHz的低頻地球自由振蕩的簡正模(3S1、0S4、0S5)及其頻譜分裂現(xiàn)象。通過對原始分鐘間隔的超導重力數(shù)據(jù)進行預處理得到重力殘差，然后運用EEMD方法將其分解到處于不同頻段上的IMF分量，提高了探測低頻自由振蕩信號的信噪比，使一些低頻自由振蕩信號的譜線分裂更清晰地顯現(xiàn)出來。結(jié)果表明，經(jīng)過EEMD分解過的重力殘差能更有效地檢測低頻自由振蕩信號，得到信噪比更高的低階球型振蕩分裂譜線。該研究進一步證實了EEMD數(shù)據(jù)處理方法的有效性以及超導重力儀在低階地球自由振蕩檢測中的優(yōu)越性。

關鍵詞：地球背景噪聲；低頻地球自由振蕩；頻率；振幅；品質(zhì)因子

地球上發(fā)生大地震時，不僅可以產(chǎn)生涉及地球局部運動的體波和面波，還能激發(fā)全球性的自由振蕩。通過研究地球自由振蕩，可以揭示地球的內(nèi)部結(jié)構(gòu)、動力學問題以及震源特性[1]。地球自由振蕩的頻率與地球的形狀、密度分布、剪切模量和體積模量等密切相關[2]。地球自轉(zhuǎn)和橢率的綜合作用使地球自由振蕩的本征頻率不是簡正的，從而產(chǎn)生相應的譜線(或譜峰)分裂耦合現(xiàn)象，同時地幔的橫向不均勻性和內(nèi)核的各向異性會導致自由振蕩譜的異常分裂[3]。此外，影響地球自由振蕩譜形態(tài)的因素還包括地球物理學結(jié)構(gòu)橫向和徑向的不均勻性、地球內(nèi)部結(jié)構(gòu)的各向異性、非完全彈性以及震源機制性質(zhì)等[2,4]。地球自由振蕩頻譜的振幅與震源機制有關[4]，我們將觀測到的自由振蕩的周期和振幅與PREM理論模型的振蕩周期和振幅相比較，就可以對地球內(nèi)部結(jié)構(gòu)進行不同于體波和面波的另一種約束[2,4]。由于地球長周期自由振蕩的振幅主要依賴于地震斷層的破裂方式及地震矩的大小，因此利用觀測記錄得到的地球自由振蕩波形資料，可以反演地震的震源破裂參數(shù)，從而約束地震的震源機制、地震大小及持續(xù)時間[2]，還可以對地球Q值進行約束[3]。例如，0S0(膨脹收縮模)的振幅大小與地震斷層的傾角、滑動方向角、震源深度及地震斷層的破裂時間有關[3]。

超導重力儀(SG)對地球自由振蕩本征周期的檢測能力可與彈簧重力儀和寬頻地震儀媲美，尤其是利用單臺SG觀測記錄就可以對低頻自由振蕩簡正模式的頻譜分裂現(xiàn)象進行高分辨率分析[5-6]，從而得到相應的本征頻率和自由振蕩參數(shù)。大量研究表明，SG還在構(gòu)制頻率小于1 mHz的長周期地震圖和研究地球深部構(gòu)造方面發(fā)揮著重要作用。

本文主要是通過對來自GGP臺網(wǎng)(http:∥www.eas.slu.edu/GGP/ggphome.html)2004年蘇門答臘地震的SG觀測數(shù)據(jù)進行預處理，從其殘差中提取出低頻率(0～1 mHz)自由振蕩模態(tài)，進而對地球深部結(jié)構(gòu)作出相應的解釋。

1　提取低頻自由振蕩信號的方法

1.1EEMD分解

本文運用EEMD方法[7]對每個臺站經(jīng)過氣壓改正和潮汐改正后的最終重力殘差時間序列進行分解(這里的EEMD方法起到了一個濾波器的作用)，得到一組有限的IMF分量(見圖1，數(shù)據(jù)長度為300 h，分解共得到9個IMF和1個殘余趨勢項,縱軸單位nm/s2)。由于不同的IMF處于不同的頻率域，因此不同范圍的重力殘差序列的信號能由合適的IMF分量重新組合。另外，EEMD能被用來調(diào)頻[7]，可用于檢波調(diào)頻時間序列。因此，EEMD可以同時用作濾波器和檢波器。

圖1　Canberra臺重力殘差在2004年蘇門答臘地震5 h后的EEMD分解結(jié)果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臺站經(jīng)過EEMD分解后的IMF1～IMF6進行功率譜密度估計得到的功率譜結(jié)果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方法提出的頻率域自回歸估計方法對自由振蕩的相關參數(shù)進行估計。為了驗證信號經(jīng)過EEMD分解后被分解到相應的頻段上，我們對重力殘差的一組IMF進行功率譜密度估計，并將它轉(zhuǎn)化到頻率域。由于IMF7～IMF9所處的頻率幾乎為0，而地球自由振蕩模態(tài)所處的頻帶主要在0.2～1 mHz，所以只需分析IMF1～IMF6。如圖2(IMF7～IMF10沒有畫出，箭頭表示相關的自由振蕩信號模態(tài),縱軸單位nm2/s4)所示，0S0出現(xiàn)在IMF2～IMF4中，0S2僅出現(xiàn)在IMF5中，0S3出現(xiàn)在IMF4中，0S5和3S1出現(xiàn)在IMF2～IMF3中。在IMF5中可以觀測到0S2的5重分裂譜峰，因此IMF5可作為觀測0S2譜線分裂的新序列。此外，觀測3S1的譜線分裂時，IMF2和IMF3可以作為一個新序列，與Shen等[9]的結(jié)論一致。結(jié)果表明，信號的IMF嚴格按照分解順序從先到后頻率依次降低，而且沒有出現(xiàn)模態(tài)混疊的現(xiàn)象，可以達到信號分解的目的。

2　低頻地球自由振蕩模態(tài)及其頻譜分裂現(xiàn)象的檢測

選取GGP觀測臺網(wǎng)2004年蘇門答臘地震中的15個超導重力儀臺站，主要通過分鐘間隔的SG數(shù)據(jù)來求解0～1 mHz頻段內(nèi)的地球低頻自由振蕩模態(tài)，從而為評估其低頻簡正模和頻譜分裂現(xiàn)象提供高精度的觀測數(shù)據(jù)。1982年Dahlen[10]提出理論上最佳的數(shù)據(jù)記錄長度約為1.1Q個周期，以增大信號的SNR。本文對于不同的地球自由振蕩模態(tài)，選取數(shù)據(jù)的起始時間和長度需視情況而定。

圖3　由疊積方法得到18個SG記錄的模態(tài)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存在的頻段，多線態(tài)1S3對3S1的觀測會產(chǎn)生交叉耦合效應。前人的作法通常是采用震后一定時間(例如2 d)的數(shù)據(jù)，以盡量削弱1S3的強度[11]。本文對每個SG記錄取地震開始50 h之后，數(shù)據(jù)總長度為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]的結(jié)果接近；采用EEMD方法得到的各譜峰的估計值與PREM模型的理論預測值的相對誤差在0.04%以內(nèi)，并且與Ding等[13]的結(jié)果接近，m=0的譜峰估計值與Roult等[12]和Shen等[14]的結(jié)果接近(表1，單位mHz)。

圖4　在0.93～0.96 mHz頻段內(nèi)4個SG數(shù)據(jù)獲得的模態(tài)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)作疊積，數(shù)據(jù)長度為555 h，相應的歸一化功率譜密度如圖5

表1　3S1分裂譜線的模型預測值與頻率觀測值

注：上標a為通過疊積方法得到相應的功率譜密度(PSD)，用Lorentz諧振函數(shù)與觀測到的譜峰進行匹配計算求得；上標b為通過EEMD分解以及Fourier變換來觀測模態(tài)的譜峰，譜峰的頻率通過AR方法來估計，最后對臺站的結(jié)果進行加權(quán)平均得到最終的頻率估計值。表2、表3同。

(垂直藍色虛線對應于PREM模型中相關譜峰的值，曲線平滑采用了3次樣條插值)所示。由于耦合效應，很難從頻譜圖中剝離出各個譜峰。對4個SG記錄(h1、h2、cb、st)進行EEMD分解以探測0S4的譜線分裂，記錄開始于地震5 h后，數(shù)據(jù)長度約為706.68 h。如圖6 (a)～(d) 所示，在使用EEMD分解之后，0S4譜峰的相對信噪比有一定程度的提高(圖6 (a)、(d) )，0S4的9個譜峰可以區(qū)分開來，這進一步證實了EEMD的可行性。從表2可以看出，疊積結(jié)果的頻率估計值與PREM理論模型預測值的相對誤差除了m=-4為-0.116%外，其他譜峰的相對誤差均在0.02%左右；使用EEMD方法得到的各譜峰的頻率估計值與PREM理論模型的預測值的相對誤差除了m=-4譜峰為-0.103%外，其他譜峰的相對誤差均在0.05%左右，其估計值與Ding等[13]的結(jié)果較為接近。盡管估計的頻率序列不同于相關PREM模型的理論預測值以及Roult等[12]研究的觀測值(在他們的研究中，只有m=±4的譜峰是2004年Sumatra地震后在超過50臺寬頻地震儀的基礎上觀測得到的)，但是估計的頻率相互之間有很強的關聯(lián)性。可以看出，運用EEMD分解技術可以清楚地觀測到所有的分裂譜線。

圖5　由疊積方法得到18個SG記錄的模態(tài)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頻段內(nèi)由4個SG數(shù)據(jù)獲得的模態(tài)0S4的歸一化振幅譜Fig.6　The normalized amplitude of the mode 0S4obtained from 4 SG records during 0.635～0.660 mHz

模態(tài)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后，數(shù)據(jù)長度為698 h，得到0S5的歸一化功率譜密度(見圖7，垂直藍色虛線對應于PREM模型中相關譜峰的值，曲線平滑采用了3次樣條插值)。由于耦合重疊的作用，很難從頻譜圖中辨別出所有譜峰。隨后選取4個SG記錄(h2、m1、st、w1)，通過EEMD分解來探測0S5的譜線分裂，記錄開始于地震5 h后，數(shù)據(jù)長度約為706.685 h。如圖8所示，通過對比使用EEMD分解前后的結(jié)果發(fā)現(xiàn)，通過對信號進行EEMD分解以及目標頻段的IMF疊加，可以清楚地觀測到大約10個譜峰，且各譜峰的相對信噪比都有所提高。使用疊積方法得到0S5的譜峰頻率估計結(jié)果與PREM模型理論預測值的相對誤差在0.04%左右。由于耦合效應，譜峰m=+3與m=+5的估計值沒有估計出來，使用EEMD分解后各譜峰的頻率估計值與PREM理論模型預測值的相對誤差也在0.04%左右，只有譜峰m=+3沒有估計出來(表3)。

圖7　由疊積方法得到18個SG記錄的模態(tài)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頻段內(nèi)由4個SG數(shù)據(jù)獲得的模態(tài)0S5的歸一化振幅譜Fig.8　The normalized amplitude of the mode 0S5obtained from 4 SG records during 0.832～0.848 mHz

模態(tài)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　結(jié)　語

本文利用EEMD技術，采用較少的SG觀測記錄，以較高的信噪比剝離了3個低頻自由振蕩振型3S1、0S4、0S5幾乎所有的分裂譜線。隨后，對各個臺站的數(shù)據(jù)進行疊積。得到的疊積譜密度圖顯示，疊積技術能很好地壓制周圍環(huán)境的噪聲，提高自由振蕩信號的相對信噪比，能有效觀測到微弱的低頻信號，但其在譜線分裂檢測方面的效果比不上EEMD方法。研究表明，一個譜峰的信噪比與其相對應的估計頻率的誤差呈負相關關系，這與H?fner等[15]及Ding等[13]得到的結(jié)果相一致。此外，本文給出的對應于這3個模態(tài)全部分裂譜峰的頻率估計與PREM模型的理論預測值有不同程度的偏離，可為構(gòu)建3D地球模型提供有效約束。

盡管典型的多重態(tài)剝離法能提供譜峰的頻率和振幅信息，但它極大地依賴于來源于不同測站接收器的結(jié)構(gòu)和記錄，有時還需要結(jié)合不同的地震事件。為了提高頻率估計的精度，我們可以采用其他方法對自由振蕩信號進行分離或?qū)EMD技術與疊積技術相結(jié)合，并改進頻率估計方法，從而既能壓制周圍噪聲的干擾，又能觀測到更高精度的簡正模的譜線分裂。

地球的自由振蕩參數(shù)，特別是振蕩信號的頻率和品質(zhì)因子Q值等是直接反映地球內(nèi)部介質(zhì)粘彈性質(zhì)的重要物理量。由于目前還很難準確描述地球自由振蕩的高階譜線分裂現(xiàn)象，這在很大程度上限制了研究進展。因此，有必要綜合利用多種數(shù)據(jù)源(包括長周期應變儀、地震儀和重力儀臺網(wǎng)資料等)深入研究高階環(huán)型模態(tài)和球型模態(tài)的譜線分裂及耦合問題，有效約束地幔結(jié)構(gòu)的橫向不均勻性和地球內(nèi)核各向異性，深化對地球內(nèi)部動力學現(xiàn)象的認識。

致謝:感謝丁浩和欒威提出有價值的建議，促使本文得以順利完成。

參考文獻

[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

第一作者簡介：曾詩瑜，碩士生，主要從事地球自由振蕩、超導重力數(shù)據(jù)分析及應用研究，E-mail: syzeng@whu.edu.cn。 通訊作者：申文斌，教授，博士，博士生導師，主要從事重力理論及應用、地球自轉(zhuǎn)、地球自由振蕩、相對論大地測量研究，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)。