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

聶仁奇，章傳銀，秘金鐘

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

0 引 言

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

1 超導重力數(shù)據

1.1 數(shù)據獲取

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

1.2 重力潮汐改正

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

1.3 數(shù)據預處理

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

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

2 基本原理與方法

2.1 基于FFT的功率譜分析

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

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

式中：s為原始時間序列；Sdt為去除趨勢項后的時間序列。

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

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

2.2 最小二乘功率譜分析

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

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

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

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

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

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

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

3 實例分析與小波分析校核

3.1 基于FFT的功率譜分析

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

圖2 一月份重力值序列

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

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

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

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

最小二乘譜分析可以直接分析不等間隔不等權的數(shù)據序列。但對于較大的噪音會對譜分析結果產生扭曲，因此，在數(shù)據處理時必須將這部分較大的噪音刪除?；谧钚《俗V分析技術，無需內插缺失的數(shù)據，無需平滑尖峰和漂移，直接由加權最小二乘譜分析得到譜值。同樣也是以2008年1月份的數(shù)據為例用最小二乘功率譜分析進行了周期探測。

圖4 最小二乘譜分析結果

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

3.3 小波多辨分析校核

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

圖5 小波多分辨分析

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

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

表1 主要周期探測結果

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

4 結 論

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

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

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

［1]尹 暉.SPIROS D P.最小二乘譜及其在超導重力觀測數(shù)據分析中的應用［J].武漢大學學報·信息科學版，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]邢 喆.超導重力儀觀測數(shù)據的初步分析與信號檢測［D].武漢：武漢大學，2004.

［4]陳 暉，周繼慧，郭春亮，等.基于VB的FFT算法的設計與實現(xiàn)［J].華東交通大學學報，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.