申文斌，欒威

1 武漢大學(xué)測(cè)繪學(xué)院，武漢　430079　2　武漢大學(xué)測(cè)繪遙感信息工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，武漢　430079

?

利用超導(dǎo)重力數(shù)據(jù)探測(cè)Slichter模三重分裂信號(hào)

申文斌1,2，欒威1

1 武漢大學(xué)測(cè)繪學(xué)院，武漢4300792武漢大學(xué)測(cè)繪遙感信息工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，武漢430079

摘要本文基于2004蘇門答臘地震后全球地球動(dòng)力學(xué)計(jì)劃觀測(cè)網(wǎng)中9臺(tái)超導(dǎo)重力儀13200 h的連續(xù)重力觀測(cè)數(shù)據(jù)，聯(lián)合利用整體經(jīng)驗(yàn)?zāi)B(tài)分解(EEMD)、最優(yōu)序列估計(jì)(OSE)、積譜密度分析(PSA)和自回歸估計(jì)(AR)方法，探測(cè)地球固態(tài)內(nèi)核平動(dòng)振蕩模態(tài)(Slichter模)的三重分裂信號(hào).利用EEMD提取常規(guī)預(yù)處理后的重力殘差中包含目標(biāo)頻段的本征模態(tài)函數(shù)(IMFs)作為最終的重力殘差之后，將每個(gè)臺(tái)站的重力殘差平均分成無(wú)重疊的3個(gè)子數(shù)據(jù)塊，再應(yīng)用OSE和PSA方法獲取Slichter模三重分裂譜線的積譜密度，并采用AR方法估計(jì)各積譜中弱共振信號(hào)的中心頻率及其誤差.結(jié)果表明，結(jié)合OSE和PSA方法探測(cè)到三個(gè)清晰的譜峰，周期分別為5.8307±7.1×10-4h、5.2161±8.1×10-4h和4.7536±5.2×10-4h，分別對(duì)應(yīng)Slichter模三重分裂譜線m=-1、m=0和m=+1，且三個(gè)譜峰的周期與Crossley(1992)、Rochester和Peng(1993)、Peng(1997)和Rogister(2003)基于PREM地球模型給出的Slichter模理論周期非常接近，也與Ding和Shen(2013)建議的“可能的觀測(cè)結(jié)果”較為接近.參考現(xiàn)有的理論和觀測(cè)結(jié)果，本文認(rèn)為這三個(gè)信號(hào)可能是Slichter模三重分裂譜線.

關(guān)鍵詞Slichter模三重分裂； 超導(dǎo)重力儀； EEMD； 最優(yōu)序列估計(jì)； 積譜密度

1引言

地球固態(tài)內(nèi)核在其平衡位置附近的平動(dòng)振蕩是地球自由振蕩基本核模，也是基本簡(jiǎn)正模之一，或稱之為Slichter模(Slichter,1961)，包括赤道面上的正向運(yùn)動(dòng)、地球自轉(zhuǎn)軸方向上的軸向運(yùn)動(dòng)和赤道面上的逆向運(yùn)動(dòng)(Slichter,1961; Smith,1976).對(duì)于非旋轉(zhuǎn)、球?qū)ΨQ、理想彈性和各向同性地球模型(SNREI)，簡(jiǎn)并的Slichter模是周期最長(zhǎng)的一階球型振蕩模態(tài)(1S1)，它的本征位移為一階球型位移.與其他的地球自由振蕩簡(jiǎn)正模一樣，地球的自轉(zhuǎn)、橢率和3-D結(jié)構(gòu)將會(huì)引起Slichter模的譜線分裂，形成一組方位角級(jí)數(shù)m分別為-1,0和1的三條譜線，分別對(duì)應(yīng)地球固態(tài)內(nèi)核平動(dòng)振蕩的三重分裂譜峰(王迪晉,2012; Ding and Shen,2013; Shen and Ding,2013).Slichter模的本征周期與內(nèi)核邊界(inner core boundary,ICB)處的密度跳躍ΔρICB有密切關(guān)聯(lián)(Crossley et al.,1992; Rieutord,2002)，因此Slichter模三重分裂譜線的本征周期觀測(cè)有利于約束地球內(nèi)外核邊界處的密度跳躍和黏滯度.

研究表明，重力測(cè)量可有效應(yīng)用于探測(cè)地球內(nèi)部構(gòu)造特征和動(dòng)力學(xué)現(xiàn)象(Melchior and Ducarme,1986).然而地球固體內(nèi)核平動(dòng)振蕩引起的地表重力變化信號(hào)非常微弱，易淹沒(méi)于Slichter模所在的亞潮汐頻段的臺(tái)站背景噪聲中，探測(cè)難度很大.研究表明(Smith,1976; Crossley et al.,1991)，大地震激發(fā)的Slichter模導(dǎo)致地表重力信號(hào)的變化幅度可達(dá)nGal量級(jí)，與超導(dǎo)重力儀(superconducting gravimeter,SG)的觀測(cè)精度(10-11m·s-2)相當(dāng).SG是目前觀測(cè)精度最高、穩(wěn)定性和連續(xù)性最好的地球重力場(chǎng)信號(hào)觀測(cè)儀器，對(duì)重力場(chǎng)變化有很高的靈敏度，在低于0.8 mHz頻段優(yōu)于當(dāng)前最好的寬頻地震儀；它在低頻段具有極低的噪聲水平，特別適用于長(zhǎng)周期信號(hào)的觀測(cè)，在研究地球自由振蕩模態(tài)方面優(yōu)勢(shì)明顯(Hinderer et al.,1995; Rosat et al.,2004,2006).相關(guān)學(xué)者在利用SG數(shù)據(jù)檢測(cè)內(nèi)核平動(dòng)振蕩方面做了許多有益的嘗試，以Smylie為代表的少數(shù)科學(xué)家基于亞地震波近似的理論計(jì)算結(jié)果，聲稱在SG的長(zhǎng)期觀測(cè)資料的亞潮汐頻段(0.172～0.285cph)發(fā)現(xiàn)了“Slichter?！比豕舱裥盘?hào)(Smylie,1992; Smylie et al.,1993; Courtier et al.,2000)，并在此基礎(chǔ)上研究了內(nèi)核邊界附近的密度和黏滯度，但他們的結(jié)果受到了很多科學(xué)家的質(zhì)疑(Hinderer et al.,1995; Jensen et al.,1995; Rosat et al.,2003).至今，是否已觀測(cè)到甚至能否觀測(cè)到Slichter模三重分裂，爭(zhēng)議很大，這也成為當(dāng)今國(guó)際地學(xué)界的一個(gè)熱點(diǎn)研究課題.

EEMD是Wu和Huang(2009)基于噪聲輔助數(shù)據(jù)分析(NADA)理論，對(duì)EMD方法(Huang et al.,1998)做出的進(jìn)一步改進(jìn)方法，其目的是減弱或消除模態(tài)混疊效應(yīng)(mode mixing effect)及端點(diǎn)效應(yīng)(end effect)，使分解得到的IMF具有更集中的頻率信息，有助于提取數(shù)據(jù)中的真實(shí)信號(hào).應(yīng)用EEMD處理SG重力殘差序列，可有效過(guò)濾重力殘差中非目標(biāo)頻段的信號(hào)，增強(qiáng)目標(biāo)模態(tài)的信噪比，有利于目標(biāo)微弱信號(hào)的提取,其有效性和優(yōu)越性參見(jiàn)Shen和Ding (2014).基于最小二乘原理的OSE方法是對(duì)MSE方法的改進(jìn)，提供了以更高的信噪比剝離單線態(tài)信號(hào)的手段，可有效應(yīng)用于探測(cè)Slichter模三重分裂.本文基于全球地球動(dòng)力學(xué)計(jì)劃(GGP)觀測(cè)網(wǎng)提供的超導(dǎo)重力儀觀測(cè)資料，利用EEMD方法提取經(jīng)常規(guī)預(yù)處理(包括潮汐改正、氣壓改正和極移改正)后的重力殘差中包含目標(biāo)頻段(0.15～0.30cph)的IMFs作為最終的重力殘差；結(jié)合OSE和PSA方法探測(cè)相應(yīng)頻段內(nèi)與地球固態(tài)內(nèi)核平動(dòng)振蕩有關(guān)的重力信號(hào)；最后應(yīng)用基于Prony方法的頻率域自回歸估計(jì)(AR)方法(Prony,1795; Chao and Gilbert,1980)，估算可能的弱共振信號(hào)的中心頻率，并結(jié)合最小二乘誤差傳播定律估算其精度.

2重力殘差序列

2.1數(shù)據(jù)準(zhǔn)備

為探測(cè)Slichter模三重分裂信號(hào)，本文選取了全球地球動(dòng)力學(xué)計(jì)劃(GGP)觀測(cè)網(wǎng)中的9個(gè)超導(dǎo)重力儀臺(tái)站在2004蘇門答臘大地震發(fā)生5小時(shí)后連續(xù)13200 h、采樣率為1h的觀測(cè)數(shù)據(jù)，這些臺(tái)站分別是德國(guó)的Bad Homburg、Moxa、Wettzell，加拿大的Canberra，奧地利的Vienna，比利時(shí)的Membach，意大利的Medicina，南非的Southland，日本的Matsushiro，相關(guān)信息見(jiàn)表1.從臺(tái)站地理位置來(lái)看，除美洲和南極洲以外，其他大洲均有臺(tái)站，排除了臺(tái)站局部區(qū)域分布的影響，在一定程度上保證了探測(cè)結(jié)果的可靠性.

2.2基于EEMD的數(shù)據(jù)預(yù)處理

為減少數(shù)據(jù)預(yù)處理帶來(lái)的人為影響，所有9個(gè)觀測(cè)序列采用相同的預(yù)處理和調(diào)和分析方法.其處理過(guò)程主要包括以下3個(gè)步驟.

(1) 利用德國(guó)Wenzel 研制的ETERNA標(biāo)準(zhǔn)潮汐分析軟件(Wenzel,1996)，精密確定重力潮汐參數(shù)(包括振幅因子和相位滯后)，并從各原始觀測(cè)序列中扣除相應(yīng)的重力潮汐信號(hào)(包括固體潮和海潮負(fù)荷效應(yīng))；采用國(guó)際上通用的大氣重力導(dǎo)納值法，得到氣壓觀測(cè)數(shù)據(jù)與其重力影響之間的對(duì)應(yīng)關(guān)系，計(jì)算出氣壓重力影響并實(shí)施改正；根據(jù)IERS(International Earth Rotation Service)提供的地球自轉(zhuǎn)參數(shù)計(jì)算極移效應(yīng)并加以扣除，得到重力殘差序列.

(2) 選擇0.15～0.30cph頻段作為目標(biāo)頻段，將EEMD作為一種并矢濾波器，提取含有目標(biāo)頻段的IMFs作為最終的重力殘差.為描述該過(guò)程，我們以意大利 Medicina臺(tái)站重力殘差序列為例，對(duì)重力殘差序列實(shí)施EEMD分解，可得到12個(gè)IMF和一個(gè)殘余趨勢(shì)項(xiàng)(圖1).為詳細(xì)了解各IMF的頻率信息，我們對(duì)每個(gè)IMF作傅利葉變換得到各自的振幅譜.圖2給出了重力殘差經(jīng)過(guò)EEMD分解后得到的前6個(gè)IMF的振幅譜.可以看到，按照分解的順序IMF的能量分別集中在由高到低分布的頻段上；在目標(biāo)頻段內(nèi)，除IMF1和IMF2以外，其他IMF的振幅幾乎為零，因此我們選擇IMF1+IMF2作為最終的重力殘差.圖3給出了EEMD處理前后的重力殘差在目標(biāo)頻段的振幅譜，圖中黑色曲線表示原始重力殘差的振幅譜，紅色虛線表示經(jīng)過(guò)EEMD處理后的重力殘差振幅譜.可以看出，原始重力殘差經(jīng)EEMD處理后，目標(biāo)頻段的噪聲水平并沒(méi)有比較明顯的降低；盡管Shen and Ding(2014)已驗(yàn)證了EEMD在探測(cè)微弱信號(hào)時(shí)能相對(duì)增強(qiáng)它們的信噪比(SNR)，但對(duì)于某些未知頻率的弱信號(hào) (如Slichter模，我們只能結(jié)合前人的理論研究結(jié)果確定一個(gè)寬泛的目標(biāo)頻段)，EEMD的效果并不明顯，但這并不妨礙我們利用EEMD過(guò)濾某些殘余的潮汐信號(hào).確切地說(shuō)，EEMD 有助于提高SNR.

表1　9個(gè)超導(dǎo)重力儀臺(tái)站的相關(guān)信息

圖1　Medicina臺(tái)站長(zhǎng)度為13200 h的重力殘差序列及其EEMD分解結(jié)果，

圖2　Medicina臺(tái)站重力殘差經(jīng)過(guò)EEMD分解后得到的前六個(gè)IMF的振幅譜(縱坐標(biāo)為振幅，單位為nm·s-2)

圖3　Medicina臺(tái)站重力殘差的振幅譜黑色實(shí)線為沒(méi)有實(shí)施EEMD處理的重力殘差振幅譜曲線，紅色虛線為EEMD處理后的重力殘差振幅譜曲線.

(3) 將EEMD處理后的各臺(tái)站最終重力殘差序列分為3個(gè)無(wú)重疊的子數(shù)據(jù)塊S1i(t)、S2i(t)、S3i(t)(i=1、2、…、9)，即每個(gè)臺(tái)站重力殘差按照時(shí)間段被分成3等份，每份數(shù)據(jù)長(zhǎng)度為4400 h，則可得到27個(gè)新的殘差序列，供后續(xù)使用.

3基于OSE方法探測(cè)Slichter模三重分裂

3.1OSE的基本原理

Ding和Shen(2013)基于最小二乘原理，對(duì)多臺(tái)站試驗(yàn)(MSE)技術(shù)作了進(jìn)一步改進(jìn)，提出了一種時(shí)間域的最優(yōu)序列估計(jì)(OSE)方法，主要解決了MSE方法中隨機(jī)噪聲影響問(wèn)題.MSE方法和OSE方法，二者均是基于目標(biāo)模態(tài)的特性而提出的方法，旨在將多重態(tài)自由振蕩信號(hào)的單線態(tài)逐個(gè)剝離，達(dá)到識(shí)別信號(hào)的目的(丁浩,申文斌,2013).

任一臺(tái)站的觀測(cè)重力殘差gj(t)包含三個(gè)平動(dòng)振蕩模信號(hào)，分別為赤道面正向、自轉(zhuǎn)軸方向和赤道面逆向振蕩(即Slichter模的三個(gè)單線態(tài)信號(hào))，可表示為(Courtier et al.,2000; Ding and Shen,2013)gj(t)=apei(ωpt-φj)sinθj+aaeiωatcosθj

+arei(ωrt+φj)sinθj+nj(t)，

(1)

其中ap、aa和ar分別表示Slichter模三重分裂譜峰的振幅，ωp、ωa和ωr為相應(yīng)的角頻率，nj(t)表示線性無(wú)關(guān)的臺(tái)站背景噪聲，θj、φj分別表示第j個(gè)臺(tái)站所在的余緯和經(jīng)度.

當(dāng)臺(tái)站數(shù)N>3時(shí)，可將式(1)表達(dá)為矩陣形式為

(2)

(3)

由方程(2)可求解未知時(shí)間序列矩陣S為

(4)

由方程(4)可以看出，B-1V是未知的.MSE方法忽略了B-1V項(xiàng)，將觀測(cè)重力殘差序列視為實(shí)數(shù)值，方程(1)的實(shí)數(shù)部分可表示觀測(cè)值，則方程(4)有唯一解(GilbertandBackus,1965;WoodhouseandGirnius,1982).盡管線性無(wú)關(guān)的臺(tái)站背景噪聲nj(t)可以通過(guò)疊積方法減弱，然而方程(4)中的B-1V項(xiàng)無(wú)法估計(jì)，因此應(yīng)用MSE方法探測(cè)Slichter模的觀測(cè)結(jié)果依然無(wú)法避免主要誤差項(xiàng)的影響.

考慮到MSE方法的限制，Ding和Shen(2013)采用最小二乘方法來(lái)估計(jì)未知時(shí)間序列apeiωpt、aaeiωat和areiωrt，方程(2)的最小二乘解為

(5)

采用OSE方法可以剝離出Slichter模的三個(gè)單線態(tài)信號(hào)，且三個(gè)獨(dú)立的時(shí)間序列可以表示為

(6)

針對(duì)三個(gè)獨(dú)立的簡(jiǎn)正模，利用基于Prony方法的頻率域的自回歸估計(jì)(AR)技術(shù)估計(jì)Slichter模三重分裂譜線的復(fù)頻率σ和復(fù)振幅A，記復(fù)頻率σ=ω+iα，其中α為衰減因子，則根據(jù)最小二乘誤差傳播定律,可估算ω、α的精度.

3.2Slichter模的檢測(cè)結(jié)果

為作比較，我們?cè)O(shè)計(jì)了以下2個(gè)方案來(lái)處理相同的3個(gè)數(shù)據(jù)塊，方案中所有功率譜密度均采用“l(fā)og”形式，得到的積譜密度單位相同，均為nm2·s-4/cph(dB).

圖4為采用方案1處理后的積譜密度.類似于Smylie(1992)、Courtier等(2000)和Ding和Shen(2013)， 我們采用五次多項(xiàng)式擬合3個(gè)分裂譜線的積譜密度來(lái)評(píng)估背景噪聲水平.從圖中可以看出，一些已知信號(hào)(包括大氣潮汐信號(hào)S4以及一些殘余潮汐信號(hào)MN4、M4、SN4、N4、MS4、K3)對(duì)目標(biāo)頻段的噪聲水平有影響.通過(guò)剔除這些已知信號(hào)，背景噪聲水平明顯降低，一些微弱的共振信號(hào)有所增強(qiáng).

圖4　采用方案1處理后的原始積譜密度(藍(lán)色曲線)與剔除某些已知信號(hào)后的積譜密度(紅色曲線)積譜密度曲線中央的藍(lán)色虛線為原始積譜密度的背景噪聲水平估計(jì)，兩側(cè)的藍(lán)色虛線為相應(yīng)的95%置信水平(CI)曲線；積譜密度曲線中央的紅色虛線為剔除某些已知信號(hào)后的積譜密度的背景噪聲水平估計(jì)，兩側(cè)的紅色虛線為相應(yīng)的95% CI.垂直虛線表示大氣潮信號(hào)S4—S7，箭頭表示某些已知的殘余潮汐信號(hào).

在檢測(cè)1S1三重分裂信號(hào)時(shí)，可能會(huì)出現(xiàn)若干組譜峰(分別對(duì)應(yīng)m=-1、m=0、m=+1的譜線分裂)，它們均有可能是Slichter 模1S1的三重分裂譜線.為判斷某一譜峰是否為可能的1S1單線態(tài)，我們參照Ding和Shen(2013)基于OSE的特征設(shè)置了兩個(gè)準(zhǔn)則： (1)該譜峰振幅Ps(ωi) (ωi為該譜峰對(duì)應(yīng)的頻率)要高于95%的置信水平(CI)曲線； (2)其他兩個(gè)積譜密度中相應(yīng)頻率的譜峰振幅Pq(ωi)(q=p，a，r; q≠s)要小于80%的CI曲線.嚴(yán)格滿足以上兩個(gè)準(zhǔn)則，則認(rèn)為該譜峰可能是某一1S1單線態(tài).

圖5給出了在0.16～0.28cph頻段內(nèi)采用方案1計(jì)算剔除某些已知信號(hào)后的積譜密度.從圖中可以清楚地看到，有8個(gè)譜峰同時(shí)滿足兩個(gè)準(zhǔn)則，圖中用不同顏色的豎線和箭頭表示出了各譜峰在該頻段內(nèi)的位置和相應(yīng)頻率，其中對(duì)應(yīng)于m=-1單線態(tài)信號(hào)的可能譜峰有1個(gè)，對(duì)應(yīng)于m=0單線態(tài)信號(hào)的可能譜峰有5個(gè)，對(duì)應(yīng)于m=+1單線態(tài)信號(hào)的可能譜峰有2個(gè)；有5個(gè)譜峰雖然滿足第一準(zhǔn)則，但不滿足第二準(zhǔn)則；因此，盡管它們的信噪比較高，我們認(rèn)為這幾個(gè)譜峰可能是臺(tái)站的背景噪聲信號(hào)，也可能是某些與其他地球物理和地球動(dòng)力學(xué)因素有關(guān)的信號(hào)(孫和平等,2004; 徐建橋等,2005).此外，在三個(gè)積譜中未發(fā)現(xiàn)任何譜峰的周期對(duì)應(yīng)于Smylie等 (1992)聲稱檢測(cè)到的1S1三重分裂譜線的周期.

我們重點(diǎn)關(guān)注可能頻段(0.17～0.26cph)內(nèi)的積譜密度，并采用Savitaky-Golay平滑濾波(Savitzky and Golay,1964)方法來(lái)處理相應(yīng)頻段的積譜密度，以凸顯信噪比高的信號(hào)，而壓制信噪比低的信號(hào).圖6給出了經(jīng)過(guò)平滑處理后的積譜密度及其相 應(yīng)的CI曲線.對(duì)比圖5和圖6，平滑處理前后積譜密度曲線整體上變化不大；圖5中與噪聲信號(hào)非常接近的一個(gè)可能譜峰信號(hào)經(jīng)過(guò)平滑后，基本與噪聲信號(hào)的頻率一致，說(shuō)明該譜峰實(shí)際上是由于該噪聲信號(hào)譜峰寬度較大，剝離在不同單線態(tài)頻譜中尖峰值可能發(fā)生偏移.圖6表明，共有7個(gè)譜峰同時(shí)滿足上述兩個(gè)準(zhǔn)則，也即這7個(gè)譜峰可能為1S1三重分裂信號(hào).用P、A、R分別表示1S1三重分裂譜線，采用AR方法估計(jì)7個(gè)譜峰的中心頻率ω0，表2給出了經(jīng)過(guò)平滑處理后檢測(cè)到的7個(gè)可能譜峰相應(yīng)的中心頻率ω0和其對(duì)應(yīng)的周期T0.雖然檢測(cè)到了可能為1S1單線態(tài)的譜峰，但沒(méi)有找到一組譜峰符合1S1三重分裂的分裂寬度，說(shuō)明這7個(gè)譜峰中可能只有2個(gè)或更少的譜峰為可能的1S1單線態(tài)，其他幾個(gè)信噪比較高的背景噪聲信號(hào)對(duì)1S1單線態(tài)干擾較大.

圖5　采用方案1計(jì)算且剔除某些已知信號(hào)后的功率譜密度黑色虛曲線為背景噪聲水平估計(jì)，藍(lán)色和紅色虛曲線分別對(duì)應(yīng)95%和80% CI曲線.黑色豎直虛線表示只滿足第一個(gè)準(zhǔn)則的相應(yīng)譜峰；不同顏色的豎直虛線表示同時(shí)滿足兩個(gè)準(zhǔn)則的譜峰，相應(yīng)顏色的箭頭表示這些可能信號(hào)在頻譜密度中出現(xiàn)的位置.其中，紅色豎直虛線對(duì)應(yīng)1S1單線態(tài)信號(hào)m=-1，藍(lán)色豎直虛線對(duì)應(yīng)1S1單線態(tài)信號(hào)m=0，綠色豎直虛線對(duì)應(yīng)1S1單線態(tài)信號(hào)m=+1.

圖6　經(jīng)過(guò)平滑后的積譜密度(只關(guān)注0.17～0.26cph的可能頻段)檢測(cè)到的譜峰與圖5相應(yīng)頻率的譜峰基本一致.黑色豎直虛線表示只滿足第一個(gè)準(zhǔn)側(cè)的相應(yīng)譜峰；不同顏色的豎直虛線表示同時(shí)滿足兩個(gè)準(zhǔn)側(cè)的譜峰，相應(yīng)顏色的箭頭表示這些可能信號(hào)在積譜密度中出現(xiàn)的位置.

譜峰ω0(cph)T0(h)P0.171505±2.09×10-55.83073±7.10×10-4AA10.175431±7.85×10-55.70026±2.55×10-3A20.177047±6.99×10-55.64822±2.23×10-3A30.191714±2.97×10-55.21611±8.08×10-4A40.218056±2.57×10-54.58597±5.39×10-4RR10.247609±2.55×10-54.03863±4.15×10-4R20.258569±2.92×10-53.86743±4.37×10-4

為完善1S1三重分裂譜線的觀測(cè)結(jié)果，我們基于Smylie等(1993)給出的PSA方法(即方案2)來(lái)檢測(cè)可能的弱共振信號(hào)，并以五次多項(xiàng)式擬合平滑后的積譜密度作為相應(yīng)頻段上背景噪聲，以譜峰峰值大于95%的CI作為可檢測(cè)信號(hào)的最低能量限制.按照此準(zhǔn)則，我們?cè)?.17～0.26cph頻段內(nèi)一共檢測(cè)到了13個(gè)譜峰(Spectral Peak，SP)，并用SPi(i=1,2,… ,13)按照頻率由高到低表示相應(yīng)的譜峰.圖7給出了平滑前后剔除某些已知信號(hào)后的積譜密度，表3給出了采用方案2檢測(cè)到的13個(gè)譜峰的中心頻率ω0及其對(duì)應(yīng)的周期T0.從圖7可以看出，經(jīng)過(guò)平滑處理后，一些振幅較高而信噪比不大的譜峰得到明顯壓制；而一些信噪比不高的信號(hào)，經(jīng)過(guò)平滑處理后依然可能達(dá)到我們?cè)O(shè)定的大于95%的CI可檢測(cè)準(zhǔn)則.比較兩種方案的檢測(cè)情況，可以看出基于OSE方法檢測(cè)到的譜峰與基于PSA方法檢測(cè)到的譜峰中心頻率相差較大.

圖7　采用方案2計(jì)算的積譜密度黑色曲線為剔除了某些已知信號(hào)后的積譜密度，紅色曲線為剔除了某些已知信號(hào)和平滑后的積譜密度；積譜密度曲線中央的黑色虛線為背景噪聲水平估計(jì)，上下兩側(cè)的藍(lán)色虛線為相應(yīng)的95% CI曲線.

譜峰ω0(cph)T0(h)SP10.170826±1.93×10-55.85391±6.60×10-4SP20.174681±2.02×10-55.72473±6.63×10-4SP30.175258±2.28×10-55.70588±7.43×10-4SP40.178269±1.94×10-55.60950±6.10×10-4SP50.180781±2.40×10-55.53157±7.34×10-4SP60.184193±3.82×10-55.42909±1.12×10-3SP70.195345±1.95×10-55.11914±5.12×10-4SP80.197837±2.19×10-55.05466±5.61×10-4SP90.198754±1.93×10-55.03135±4.89×10-4SP100.210366±2.28×10-54.75363±5.15×10-4SP110.241721±2.03×10-54.13700±3.47×10-4SP120.249664±2.13×10-54.00538±3.41×10-4SP130.256722±1.96×10-53.89527±2.97×10-4

綜合考慮前人的理論計(jì)算結(jié)果，Slichter模的周期約為3～6 hours，分裂寬度約為15～40 minutes(Dahlen and Sailor,1979; Crossley,1992; Hinderer and Crossley,1993; Rochester and Peng,1993; Peng,1995,1997; Denis et al.,1997; Rogister,2003; Smylie et al.,2009; Rosat,2011)；而對(duì)于PREM模型，Slichter模的中心周期約為5.2～5.3 h(Peng,1995,1997; Denis et al.,1997; Rogister,2003; Smylie et al.,2009; Rosat,2011).鑒于OSE方法檢測(cè)1S1三重分裂信號(hào)更為嚴(yán)格，我們以O(shè)SE方法檢測(cè)的結(jié)果為基礎(chǔ)，結(jié)合兩種方案檢測(cè)結(jié)果，并參照理論模型結(jié)果，找到了唯一一組可能的三重分裂信號(hào)(見(jiàn)表4)，它們是基于OSE方法檢測(cè)到的2個(gè)譜峰P(m=-1：5.8307±7.1×10-4h)、A3(m=0：5.2161±8.1×10-4h)和基于PSA方法檢測(cè)到的1個(gè)譜峰SP10(m=+1：4.7536±5.2×10-4h)，這三個(gè)譜峰所對(duì)應(yīng)的周期與Crossley(1992)、Rochester和Peng(1993)、Peng(1997)和Rogister(2003)基于PREM模型給出的理論周期之間的相對(duì)誤差約1.48%，與Ding和Shen(2013) 建議的“可能的觀測(cè)結(jié)果”之間的相對(duì)誤差約為0.59%；其中，m=+1信號(hào)的強(qiáng)度明顯高于m=0和m=-1信號(hào).

4結(jié)論及討論

(1) 基于2004年蘇門答臘地震后GGP觀測(cè)網(wǎng)中9臺(tái)超導(dǎo)重力儀13200h的連續(xù)重力觀測(cè)數(shù)據(jù)，將EEMD作為濾波器處理常規(guī)預(yù)處理后的重力殘差， 提取其中包含目標(biāo)頻段的IMFs，作為最終各臺(tái)站的重力殘差，有利于提高信噪比.結(jié)合時(shí)間域內(nèi)的OSE方法和頻率域內(nèi)PSA方法，迭積了9個(gè)超導(dǎo)重力臺(tái)站重力殘差的功率譜密度，并對(duì)積譜采用剔除某些已知信號(hào)和平滑處理；參照Ding和Shen(2013)設(shè)置的兩個(gè)檢測(cè)準(zhǔn)則，并依據(jù)Slichter模的分裂特征選取了兩個(gè)候選譜峰，對(duì)應(yīng)于Slichter模的兩個(gè)單線態(tài)信號(hào)m=-1和m=0.只采用頻率域內(nèi)PSA方法處理相同的重力殘差序列，我們找到13個(gè)微弱譜峰，對(duì)照基于OSE方法得到的兩個(gè)單線態(tài)并兼顧理論分裂寬度，從中選取了一個(gè)候選譜峰，對(duì)應(yīng)于Slichter模的單線態(tài)信號(hào)m=+1.結(jié)合兩種處理過(guò)程的檢測(cè)結(jié)果，我們探測(cè)到了一組結(jié)果(5.8307±7.1×10-4h，5.2161±8.1×10-4h，4.7536±5.2×10-4h)，不僅與Crossley(1992)、Rochester和Peng(1993)、Peng(1997)和Rogister(2003)基于PREM模型給出的理論周期很接近，而且也與 Ding和Shen(2013)的探測(cè)結(jié)果非常接近.需要指出，本研究與Ding和Shen(2013)的研究完全獨(dú)立，盡管檢測(cè)Slichter模的方法相同，但采用的數(shù)據(jù)源不同，且數(shù)據(jù)預(yù)處理的方法也不同.這一獨(dú)立探測(cè)結(jié)果，從某種意義上說(shuō)，證實(shí)了Ding和Shen(2013)探測(cè)到的“可能的Slichter 模三重分裂值”.

表4　Slichter模的理論預(yù)測(cè)周期及本文觀測(cè)周期(單位：h)

(2) 最后需要指出，Rosat和Rogister (2012)認(rèn)為，只有震級(jí)超過(guò)MW9.7的大地震才能激發(fā)出足以觀測(cè)得到的Slichter模三重分裂信號(hào)，但他們的結(jié)論有待商榷，因?yàn)橐粋€(gè)顯然的事實(shí)是，通過(guò)多臺(tái)站疊積技術(shù)可有效放大全球性信號(hào)，但他們忽視了這一點(diǎn).不過(guò)，毫無(wú)疑問(wèn)，由于我們并不能象探測(cè)如2S1模態(tài)那樣看到清晰、明顯的三個(gè)譜峰，因此，我們尚無(wú)法判定本文所檢測(cè)到的三個(gè)譜峰是否一定是地球內(nèi)核的平動(dòng)振蕩信號(hào).此外，還有一些譜峰信號(hào)其原因不明，這些信號(hào)是否由地球內(nèi)核平動(dòng)振蕩所致，還是與其他一些地球物理過(guò)程(或事件)如海浪、潮汐效應(yīng)等有關(guān)，還有待進(jìn)一步研究.

致謝GGP臺(tái)網(wǎng)提供了高精度SG數(shù)據(jù)，丁浩博士在數(shù)據(jù)處理方面提供了寶貴意見(jiàn)，審稿人對(duì)初稿提出了有價(jià)值的評(píng)論和修改建議,致使本文更加完善，在此向他們表示衷心感謝.

References

Buland R,Berger J,Gilbert F.1979.Observations from the IDA network of attenuation and splitting during a recent earthquake.Nature,277(5695): 358-362,doi: 10.1038/277358a0.

Chao F B,Gilbert F.1980.Autoregressive estimation of complex eigenfrequencies in low frequency seismic spectra.Geophys.J.Int.,63(3): 641-657.

Courtier N,Ducarme B,Goodkind J,et al.2000.Global superconducting gravimeter observations and the search for the translational modes of the inner core.Phys.Earth Planet.Int.,117(1-4): 3-20,doi: 10.1016/S0031-9201 (99)00083-7.

Crossley D J.1992.Eigensolutions and seismic excitation of the Slichter mode triplet for a fully rotating Earth model.EOS,73,60.Crossley D J,Hinderer J,Legros H.1991.On the excitation,detection and damping of core modes.Phys.Earth Planet.Int.,68(1-2): 97-116.

Crossley D,Rochester M G,Peng Z R.1992.Slichter modes and Love numbers.Geophys.Res.Lett.,19(16):1679-1682.

Dahlen F A,Sailor R V.1979.Rotational and elliptical splitting of the free oscillation of the Earth.Geophys.J.Int.,58(3): 609-623,doi: 10.1111/j.1365-246X.1979.tb04797.x.

Denis C,Rogister Y,Amalvict M,et al.1997.Hydrostatic flattening,core structure,and translational mode of the inner core.Phys.Earth Planet.Int.,99(3-4): 195-206.

Ding H,Shen W B.2013.Search for the Slichter modes based on a new method: Optimal sequence estimation.J.Geophys.Res.,118(9): 5018-5029,doi: 10.1002/jgrb.50344.

Ding H,Shen W B.2013.Comparative study of the MSE technique and SHS method used for detectionnS1mode triplet.Chinese J.Geophys.(in Chinese),56(10): 3313-3323,doi: 10.6038/cjg20131008.

Gilbert F.1971.The diagonal sum rule and averaged eigenfrequencies.Geophys.J.Int.,23(1): 119-123,doi: 10.1111/j.1365-246X.1971.tb01806.x.Gilbert F,Backus G.1965.The rotational splitting of the free oscillations of the Earth 2.Rev.Geophys.,3(1): 1-9.

Hinderer J,Crossley D.1993.Core dynamics and surface gravity changes.∥ Dynamics of Earth′s Deep Interior and Earth Rotation.AGU,Geophysical Monograph,72: 1-16.Hinderer J,Crossley D,Jensen O.1995.A search for the Slichter triplet in superconducting gravimeter data.Phys.Earth Planet.Int.,90(3-4): 183-195.

Huang N E,Shen Z,Long S R,et al.1998.The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis.Proc.Roy.Soc.Lond.A.,454(1971): 903-995.

Jensen O,Hinderer J,Crossley D J.1995.Noise limitations in the core-mode band of superconducting gravimeter data.Phys.Earth Planet.Int.,90(3-4): 169-181.

Masters G,Laske G,Gilbert F.2000.Matrix autoregressive analysis of free-oscillation coupling and splitting.Geophys.J.Int.,143(2): 478-489.

Melchior P,Ducarme B.1986.Detection of inertial gravity oscillations in the Earth′s core with a superconducting gravimeter at Brussels.Phys.Earth Planet.Int.,42(3): 129-134.

Peng Z R.1995.The Slichter modes in a realistic Earth model [PhD thesis].Memorial University of Newfoundland.

Peng Z R.1997.Effects of a mushy transition zone at the inner core boundary on the Slichter modes.Geophys.J.Int.,131(3): 607-617,doi: 10.1111/j.1365-246X.1997.tb06602.x.

Prony R.1795.Essai experimental et analytique sur les lois de la dilatabilité des fluids élastiques et sur celles de la force expansive de la vapeur de l′eau et de la vapeur de l′alkool,à différentes températures.J.Ec.Polytech.(Paris),1(2): 24-76.

Rieutord M.2002.Slichter modes of the Earth revisited.Phys.Earth Planet.Int.,131(3-4): 269-278.

Rochester M G,Peng Z R.1993.The slichter modes of the rotating Earth: A test of the subseismic approximation.Geophys.J.Int.,113(3): 575-585.

Rogister Y.2003.Splitting of seismic-free oscillations and of the Slichter triplet using the normal mode theory of a rotating,ellipsoidal Earth.Phys.Earth Planet.Int.,140(1-3): 169-182.

Rosat S.2011.A review of the Slichter modes: an observational challenge.∥ Phillips J M.The Earth′s Core: Structure,Properties and Dynamics.Nova Science Publishers,Inc.,63-77.Rosat S,Fukushima T,Sato T,et al.2008.Application of a non-linear damped harmonic analysis method to the normal modes of the Earth.J.Geodyn.,45(1): 63-71.

Rosat S,Hinderer J,Crossley D,et al.2003.The search for the Slichter mode: comparison of noise levels of superconducting gravimeters and investigation of a stacking method.Phys.Earth Planet.Int.,140(1-3): 183-202,doi: 10.1016/j.pepi.2003.07.010.

Rosat S,Hinderer J,Crossley D,et al.2004.Performance of superconducting gravimeters from long-period seismology to tides.J.Geodyn.,38(3-5): 461-476.

Rosat S,Rogister Y.2012.Excitation of the Slichter mode by collision with a meteoroid or pressure variations at the surface and core boundaries.Phys.Earth Planet.Int.,190-191: 25-33,doi: 10.1016/j.pepi.2011.10.007.

Rosat S,Rogister Y,Crossley D,et al.2006.A search for the Slichter triplet with superconducting gravimeters: Impact of the density jump at the inner core boundary.J.Geodyn.,41(1-3): 296-306,doi: 10.1016/j.jog.2005.08.033.

Savitzky A,Golay M J E.1964.Smoothing and differentiation of data by simplified least squares procedures.Anal.Chem.,36(8): 1627-1639.

Shen W B,Ding H.2013.Detection of the inner core translational triplet using superconducting gravimetric observations.Journal of Earth Science,24(5): 725-735,doi: 10.1007/s12583-013-0369-3.

Shen W B,Ding H.2014.Observation of spheroidal normal mode multiplets below 1 mHz using ensemble empirical mode decomposition.Geophys.J.Int.,196(3): 1631-1642.Slichter L B.1961.The fundamental free mode of the Earth′s inner core.Proc.Natl.Acad.Sci.USA,47(2): 186-190.

Smith M L.1976.Translational inner core oscillations of a rotating,slightly elliptical Earth.J.Geophys.Res.,81(17): 3055-3065.

Smylie D E.1992.The inner core translational triplet and the density near Earth′s center.Science,255(5052): 1678-1682,doi: 10.1126/science.255.5052.1678.

Smylie D E,Brazhkin V V,Palmer A.2009.Direct observations of the viscosity of Earth′s outer core and extrapolation of measurements of the viscosity of liquid iron.Phys.-Usp.,52(1): 79-92.Smylie D E,Hinderer J,Richter B,et al.1993.The product spectra of gravity and barometric pressure in Europe.Phys.Earth Planet.Int.,80(3-4): 135-157.

Smylie D E,Jiang X H,Brennan B J,et al.1992.Numerical calculation of modes of oscillation of the Earth′s core.Geophys.J.Int.,108(2): 465-490.

Sun H P,Xu J Q,Ducarme B.2004.Detection of the translational oscillation of the Earth′s solid inner core based on the international superconducting gravimeter observations.Chinese Science Bulletin,49(11): 1165-1176.

Wang D J.2012.Study of low-frequency earth free oscillations based on global superconducting gravimeter data [Ph.D.thesis] (in Chinese).Wuhan: Wuhan University,68-112.Wenzel H G.1996.The nanogal software: earth tide data processing package ETERNA3.30.Bull.Inf.Marees Terrestres,124: 9425-9439.

Woodhouse J H,Girnius T P.1982.Surface waves and free oscillations in a regionalized Earth model.Geophys.J.Int.,68(3): 653-673.

Wu Z H,Huang N E.2009.Ensemble empirical mode decomposition: a noise-assisted data analysis method.Adv.Adapt.Data Anal.,1(1): 1-41.

Xu J Q,Sun H P,Fu R S.2005.Detection of long-period core modes by using the data from global superconducting gravimeters.Chinese J.Geophys.(in Chinese),48(1): 69-77.

附中文參考文獻(xiàn)

丁浩,申文斌.2013.探測(cè)一階模態(tài)三重分裂的MSE技術(shù)和SHS方法的比較研究.地球物理學(xué)報(bào),56(10): 3313-3323,doi: 10.6038/cjg20131008.

孫和平,徐建橋,Ducarme B.2004.基于國(guó)際超導(dǎo)重力儀觀測(cè)資料檢測(cè)地球固態(tài)內(nèi)核的平動(dòng)振蕩.科學(xué)通報(bào),49(8): 803-813.

王迪晉.2012.基于全球超導(dǎo)重力數(shù)據(jù)的低頻地球自由振蕩研究[博士論文].武漢: 武漢大學(xué),68-112.

徐建橋,孫和平,傅容珊.2005.利用全球超導(dǎo)重力儀數(shù)據(jù)檢測(cè)長(zhǎng)周期核模.地球物理學(xué)報(bào),48(1): 69-77.

(本文編輯張正峰)

基金項(xiàng)目國(guó)家重點(diǎn)基礎(chǔ)研究發(fā)展計(jì)劃(973)(2013CB733305)和國(guó)家自然科學(xué)基金項(xiàng)目(41574007，41210006,41174011,41429401)聯(lián)合資助.

作者簡(jiǎn)介申文斌,博士,武漢大學(xué)教授,博士生導(dǎo)師,主要從事物理大地測(cè)量、相對(duì)論大地測(cè)量和地球物理研究.E-mail:wbshen@sgg.whu.edu.cn

doi:10.6038/cjg20160307 中圖分類號(hào)P312

收稿日期2014-07-14，2015-09-08收修定稿

Detection of the Slichter mode triplet using superconducting gravimetric observations

SHEN Wen-Bin1,2,LUAN Wei1

1SchoolofGeodesyandGeomatics,WuhanUniversity,Wuhan430079,China2StateKeyLaboratoryofInformationEngineeringinSurveying,MappingandRemoteSensing,WuhanUniversity,Wuhan430079,China

AbstractBased on 13200-hour records by nine superconducting gravimeters (SGs) under a worldwide network of the Global Geodynamics Project after the 2004 Sumatra earthquake,a series methods including the ensemble empirical mode decomposition (EEMD),the optimal sequence estimation (OSE),the product spectrum analysis (PSA) and the autoregressive (AR) are applied jointly to detect the Earth′s inner core translational mode triplet (Slichter mode triplet).We first obtain the intrinsic mode functions (IMFs) of the target frequency band extracted from the residual sequence after the routine pretreatment as the final gravity residual,and each of the nine residual gravity series is divided into three segments of the same length without overlap,and the product spectra of Slichter triplet are obtained by using OSE and PSA.Then,the AR method is used to estimate the central frequencies of the week resonance signals and the corresponding errors.The results show that three clear spectral peaks are detected using OSE and PSA,which correspond respectively to the singlets m=-1,m=0 and m=+1.The periods corresponding to the three spectral peaks are close to the prediction of the PREM model given by Crossley (1992),Rochester and Peng (1993),Peng (1997) and Rogister (2003),as well as the “possible observation results” suggested by Ding and Shen (2013).Referring to the existing theory and observed results of the Slichter triplet,we conclude that the three observed signals in this study might correspond to the Slichter triplet.

KeywordsSlichter mode triplet; Superconducting gravimeter; EEMD; Optimal sequence estimation; Product spectrum analyze

申文斌，欒威.2016.利用超導(dǎo)重力數(shù)據(jù)探測(cè)Slichter模三重分裂信號(hào).地球物理學(xué)報(bào)，59(3):840-851,doi:10.6038/cjg20160307.

Shen W B,Luan W.2016.Detection of the Slichter mode triplet using superconducting gravimetric observations.Chinese J.Geophys.(in Chinese),59(3):840-851,doi:10.6038/cjg20160307.