超導(dǎo)重力儀觀測數(shù)據(jù)分析①

聶仁奇，章傳銀，秘金鐘

（1.中國測繪科學(xué)研究院，北京100830；2.山東科技大學(xué)，山東 青島266510）

0 引 言

地球時(shí)變重力場包含地球系統(tǒng)物質(zhì)分布及其運(yùn)移的豐富信息，直接反映了地球各圈層最基本的物質(zhì)及其變化特性，是各種環(huán)境變化導(dǎo)致地球動(dòng)力學(xué)特征最基本、最直接和最重要的物理量之一。時(shí)變重力場觀測信號(hào)是深入研究、了解全球和局部動(dòng)力學(xué)特征、固體地球內(nèi)部的構(gòu)造和運(yùn)動(dòng)特征、各圈層的相互作用和耦合機(jī)制的重要信息源［1]。目前超導(dǎo)重力儀具有很高的靈敏度和穩(wěn)定性、很低的噪音水平和漂移以及很寬的動(dòng)態(tài)頻率響應(yīng)范圍，是目前觀測精度最高的重力儀。因此，超導(dǎo)重力儀作為研究時(shí)變重力場的關(guān)鍵儀器。探測時(shí)間序列中隱含的周期信號(hào)是超導(dǎo)重力儀觀測數(shù)據(jù)分析的主要目標(biāo)之一，通過快速傅里葉變換（FFT）和最小二乘譜分析方法分析超導(dǎo)重力數(shù)據(jù)中隱含信號(hào)的檢測，兩種方法結(jié)果進(jìn)行了分析比較。最終用小波分析的方法進(jìn)行了檢校。

1 超導(dǎo)重力數(shù)據(jù)

1.1 數(shù)據(jù)獲取

獲取的數(shù)據(jù)是一臺(tái)GWR臺(tái)站式超導(dǎo)重力儀的數(shù)據(jù)，以1min的采樣頻率進(jìn)行數(shù)據(jù)采集，并以1d作為一個(gè)記錄文件，觀測數(shù)據(jù)為電壓（V），這些數(shù)據(jù)乘以一個(gè)系數(shù)（不同的測站其系數(shù)不同）即可轉(zhuǎn)化為重力值（單位μGal）。這個(gè)超導(dǎo)重力儀的系數(shù)為－76.504μGal／V.

1.2 重力潮汐改正

受日月引潮位作用影響重力產(chǎn)生周期性的潮汐變化即為重力潮汐。在重力觀測中包含了非常顯著的重力潮汐變化，其幅度可以超過200，且表現(xiàn)為非常強(qiáng)烈的緯度依賴特征，如果不精確改正，就會(huì)嚴(yán)重影響對(duì)地球的動(dòng)力學(xué)和環(huán)境變化的研究。由重力潮汐預(yù)測程序G－wave計(jì)算得到重力潮汐理論結(jié)果，可以用來剔除超導(dǎo)重力儀原始觀測數(shù)據(jù)中重力潮汐信號(hào)的影響。G－wave程序的預(yù)測精度（在日月引力作用下）為40μGal.在計(jì)算前，應(yīng)輸入測站的經(jīng)度和緯度，根據(jù)理論計(jì)算結(jié)果，可以從原始數(shù)據(jù)中剔除重力潮汐的影響［2]。

1.3 數(shù)據(jù)預(yù)處理

由于超導(dǎo)重力儀受到環(huán)境的影響（地震，地下水變化）和儀器因素（氦注入、故障、維修、校準(zhǔn)、系統(tǒng)升級(jí)等），超導(dǎo)重力儀的數(shù)據(jù)是不連續(xù)不穩(wěn)定的。在原始數(shù)據(jù)中存在著數(shù)據(jù)的重疊、斷鏈、尖峰、突跳、漂移及噪聲干擾。故在分析前必須刪除重疊數(shù)據(jù)，使數(shù)據(jù)序列的時(shí)間單位呈單調(diào)遞增。如圖1所示，三月份數(shù)據(jù)中出現(xiàn)了四次比較大的噪音，對(duì)結(jié)果會(huì)產(chǎn)生很大的扭曲，因此，必須將這部分?jǐn)?shù)據(jù)刪除［3]。

圖1 三月份重力數(shù)據(jù)

2 基本原理與方法

2.1 基于FFT的功率譜分析

FFT是進(jìn)行超導(dǎo)重力觀測數(shù)據(jù)中隱含信號(hào)檢測的常用經(jīng)典方法。FFT算法提供了一種快速頻譜分析方法，可直接用來處理離散信號(hào)數(shù)據(jù)。

FFT在處理時(shí)間序列數(shù)據(jù)時(shí)要求數(shù)據(jù)必須等間隔、等權(quán)、連續(xù)和平穩(wěn)。但在超導(dǎo)重力數(shù)據(jù)采集工程中難免出現(xiàn)斷鏈、尖峰、噪音、漂移等特點(diǎn)。造成了數(shù)據(jù)的不等間隔、不等權(quán)、不連續(xù)以及不平穩(wěn)的情況。所以，在用FFT譜分析處理超導(dǎo)重力數(shù)據(jù)時(shí)，必須要對(duì)超導(dǎo)重力事先處理的觀測數(shù)據(jù)進(jìn)行數(shù)據(jù)預(yù)處理。數(shù)據(jù)預(yù)處理的目的在于保證所需信號(hào)不受損失，消除干擾因素，從而得到高精度的觀測資料。譜分析和時(shí)間序列分析都要求原始數(shù)據(jù)為零均值，沒有線性趨勢(shì)項(xiàng)。為此，需要對(duì)時(shí)間序列進(jìn)行線性擬合，即去除時(shí)間序列的趨勢(shì)項(xiàng)。

式中：s為原始時(shí)間序列；Sdt為去除趨勢(shì)項(xiàng)后的時(shí)間序列。

利用FFT進(jìn)行譜分析時(shí)要求均勻采樣，但在超導(dǎo)重力時(shí)間序列中難免出現(xiàn)斷鏈的情況。為了獲得均勻采樣的時(shí)間序列，需要對(duì)時(shí)間序列進(jìn)行插值。插值采用拉格朗日插值進(jìn)行。由此獲得了經(jīng)預(yù)處理的時(shí)間序列［4]。

對(duì)于數(shù)據(jù)預(yù)處理獲得的時(shí)間序列，達(dá)到了等間隔、等權(quán)、連續(xù)和平穩(wěn)的要求。用離散傅里葉變換計(jì)算功率譜

2.2 最小二乘功率譜分析

用FFT進(jìn)行頻譜分析，研究時(shí)間序列中可能存在的周期信號(hào)，并比較各種頻率波動(dòng)的振幅或功率，可以確定主要周期。但FFT的不足主要表現(xiàn)在：它要求數(shù)據(jù)必須等間隔、等權(quán)、連續(xù)和平穩(wěn)，對(duì)于不滿足這些條件的數(shù)據(jù)處理，F(xiàn)FT采用數(shù)據(jù)內(nèi)插（如線性內(nèi)插）使間斷的數(shù)據(jù)連續(xù)，使不等間距數(shù)據(jù)等間距以及提取趨勢(shì)項(xiàng)、修正數(shù)據(jù)漂移等。針對(duì)FFT的不足和局限性，Vanicek于1971年提出的最小二乘譜分析（LSSA），將最小二乘逼近原理用于譜分析中，解決了測量數(shù)據(jù)采集中經(jīng)常出現(xiàn)的諸如數(shù)據(jù)序列不等間隔、不等權(quán)、不平穩(wěn)、有斷鏈、尖峰、漂移等需要進(jìn)行數(shù)據(jù)預(yù)處理的問題。Pagiatakis給出了數(shù)據(jù)序列不等權(quán)情況下的LSSA公式，提出了用數(shù)理統(tǒng)計(jì)理論檢驗(yàn)最小二乘譜的顯著性方法［5－6]，并總結(jié)出該方法具有如下優(yōu)越性：①系統(tǒng)誤差可以得到嚴(yán)格抑制，而不產(chǎn)生譜峰的漂移；②不等間距的序列不需要預(yù)處理；③可以處理具有協(xié)方差矩陣的時(shí)間序列；④能進(jìn)行譜峰顯著性的統(tǒng)計(jì)檢驗(yàn)［7]。

在譜分析中，通常將尋找的隱含周期信號(hào)表示為一組以正弦和余弦為基函數(shù)的線性組合。給定一組觀測值向量f＝f（t）＝｛fi｝，i＝1，2，…，n，若設(shè)定基函數(shù)形式為一組頻率ωi（i＝1，2，…，m）的三角函數(shù)，則有

其中，Ω 為未知參數(shù)向量，設(shè) Ω＝［Ω1i，Ω2i]T，Φ＝［cos（ωit），sin（ωit）]，其中Ω 由公式

其中，Φ為已知基函數(shù)向量。

對(duì)于不同的頻率ωi（i＝1，2，…，m ），可以得到不同的頻譜值

這就是最小二乘譜分析的基本思想。顯然，f的最小二乘譜是所有頻率ωi（i＝1，2，…，m）下的譜值集合，ωi譜值越大，表明f在該頻率的功率越大。

具有協(xié)方差函數(shù)Cf時(shí)的最小二乘譜計(jì)算公式為

3 實(shí)例分析與小波分析校核

3.1 基于FFT的功率譜分析

以2008年1月份的數(shù)據(jù)為例，進(jìn)行功率譜的計(jì)算。如圖2示出的是2008年1月經(jīng)過去粗差、去趨勢(shì)項(xiàng)、插值形成的重力值序列。

圖2 一月份重力值序列

基于FFT技術(shù)，對(duì)預(yù)處理后的時(shí)間序列進(jìn)行譜分析得到譜值。圖3詳細(xì)示出了2008年1月重力數(shù)據(jù)功率譜圖。

圖3 一月份重力數(shù)據(jù)功率譜

從圖3中可以看出四個(gè)比較清晰的功率譜峰值，即具有四個(gè)很明顯的周期項(xiàng)。分別為以67.635h、23.951 5h、11.999 9h為周期和以8.001 2為周期的四個(gè)周期項(xiàng)。

3.2 基于最小二乘功率譜分析

最小二乘譜分析可以直接分析不等間隔不等權(quán)的數(shù)據(jù)序列。但對(duì)于較大的噪音會(huì)對(duì)譜分析結(jié)果產(chǎn)生扭曲，因此，在數(shù)據(jù)處理時(shí)必須將這部分較大的噪音刪除。基于最小二乘譜分析技術(shù)，無需內(nèi)插缺失的數(shù)據(jù)，無需平滑尖峰和漂移，直接由加權(quán)最小二乘譜分析得到譜值。同樣也是以2008年1月份的數(shù)據(jù)為例用最小二乘功率譜分析進(jìn)行了周期探測。

圖4 最小二乘譜分析結(jié)果

從圖4中可以看出四個(gè)比較清晰的功率譜峰值，即具有四個(gè)很明顯的周期項(xiàng)，這一點(diǎn)和FFT譜分析得到的結(jié)果是一致的。從左至右依次為以68.133h、23.991 7h、以12.162 3h為周期和以8.183 1為周期的四個(gè)周期項(xiàng)。

3.3 小波多辨分析校核

為了再次說明兩種方法的可行性和準(zhǔn)確性，利用目前成熟的小波多分辨分析方法對(duì)兩種計(jì)算結(jié)果也進(jìn)行了校核。

圖5 小波多分辨分析

圖5是超導(dǎo)重力數(shù)據(jù)時(shí)間序列的小波多分辨分析得到的結(jié)果。圖5中第一行、第二行、第三行、第四行表示的是周期分別為8.030 2h、12.000 0 h、23.936 7h、68h的信號(hào)。和由FFT譜分析方法和最小二乘譜分析方法計(jì)算得到的結(jié)果十分吻合。

分別利用FFT譜分析方法和最小二乘譜方法對(duì)2008年1月的超導(dǎo)重力數(shù)據(jù)進(jìn)行了分析。主要的周期探測結(jié)果如表1所示。

表1 主要周期探測結(jié)果

從表1可以看出，利用兩種方法計(jì)算得到的周期結(jié)果與理論結(jié)果基本一致，說明了這兩種周期探測方法在超導(dǎo)重力數(shù)據(jù)探測周期中的正確性。

4 結(jié) 論

基于超導(dǎo)重力數(shù)據(jù)研究時(shí)變重力場、重力潮汐、自由核章動(dòng)參數(shù)測定、海潮模型有效性檢驗(yàn)、大氣負(fù)荷信號(hào)改正、地球自由振蕩檢測等方面發(fā)揮了重要作用。超導(dǎo)重力儀觀測序列具有豐富的信息，蘊(yùn)藏了地球內(nèi)部變化的所有動(dòng)態(tài)變量的痕跡。分別利用FFT譜分析方法和最小二乘譜分析方法對(duì)2008年1月份的數(shù)據(jù)進(jìn)行了周期探測。兩種方法的計(jì)算結(jié)果十分吻合，且與理論值也十分接近。總體來看：

1）快速傅里葉變換（FFT）譜分析，計(jì)算簡單，處理海量數(shù)據(jù)時(shí)有較大優(yōu)勢(shì)，可快速探測出觀測序列中的隱含周期信號(hào)，其計(jì)算結(jié)果與理論值基本吻合。

2）最小二乘譜分析，對(duì)存在斷鏈、不等間隔和不等權(quán)數(shù)據(jù)的頻譜分析有較大優(yōu)勢(shì)。

［1]尹 暉.SPIROS D P.最小二乘譜及其在超導(dǎo)重力觀測數(shù)據(jù)分析中的應(yīng)用［J].武漢大學(xué)學(xué)報(bào)·信息科學(xué)版，2005，30（7）：613－616.

［2]MERRIAM J B.An ephemeris for gravity tide predic－tions at the Nanogal Level［J].Geophys.J.Int.，1992（108）：415－422.

［3]邢 喆.超導(dǎo)重力儀觀測數(shù)據(jù)的初步分析與信號(hào)檢測［D].武漢：武漢大學(xué)，2004.

［4]陳 暉，周繼慧，郭春亮，等.基于VB的FFT算法的設(shè)計(jì)與實(shí)現(xiàn)［J].華東交通大學(xué)學(xué)報(bào)，2003，20（1）：94－96.

［5]PAGIATAKIS S D.Stochastic significance of peaks in the least－squares spectrum［J].Journal of Geodesy，1999（73）：67－78.

［6]CHEN Y Q.Deformation observation data process［M].Beijing：Press of Surveying and Mapping，1988.

［7]WELL D E，VANICEK P，PAGIATAKIS S D.Least squares spectral analysis revisited.Tech.Rep.84，Department of Surveying Engineering［R].University of New Brunswick，F(xiàn)redericton，1985.