GOCE引力梯度數據濾波方法的設計與結果分析*

楊洪國 黃 強

1)青島農業大學，青島 266109

2)西南交通大學地球科學與環境工程學院，成都 610097

GOCE引力梯度數據濾波方法的設計與結果分析*

楊洪國1)黃 強2)

1)青島農業大學，青島 266109

2)西南交通大學地球科學與環境工程學院，成都 610097

為了濾掉GOCE引力梯度數據的低頻誤差，基于引力梯度數據的頻譜特性設計了3種濾波方法，并選取其中濾波效果最優的一種方法處理了2009-11-02～2010-01-10共70 d的引力梯度數據。利用濾波后的70 d引力梯度數據恢復的1個200階次的重力場模型結果顯示，引力梯度數據恢復的模型必須進行正則化處理，施加了正則化的模型與GO_CONS_GCF_2_TIM_R3在160階前具有一致的精度。

GOCE;引力梯度;低頻誤差;快速傅里葉變換;帶通濾波器

GOCE衛星的目的是恢復高精度靜態地球重力場模型。ESA預期可恢復大地水準面精度為1 cm、重力異常精度為10－5ms－2、空間分辨率為100 km(半波長)的地球重力場模型［1－3］。由于GOCE衛星自身的技術特性，導致6個梯度分量中的4個(Vxx，Vyy，Vzz，Vxz)精度較高，另外兩個(Vxy，Vyz)精度較差［4］。此外，引力梯度分量含有巨大的低頻誤差，低頻誤差的濾波處理直接影響到恢復重力場模型的計算。目前，ESA已發布了大量的引力梯度數據，包括兩種數據產品:EGG_TRF_2_和EGG_NOM_2_。其中，EGG_TRF_2_為經過高通濾波處理、在局部指北坐標系LNOF(local north oriented frame)下的引力梯度數據;EGG_NOM_2_為未經濾波處理，在梯度儀坐標系GRF(gradiometer reference frame)下的引力梯度數據。由于EGG_TRF_2_在處理時引入了大量的外部模型信息，并在坐標系之間進行了轉換，受到姿態誤差的影響，不能用于恢復重力場模型。EGG_NOM_2_可用于恢復重力場模型的計算，但必須進行濾波處理。本文主要對引力梯度數據的噪聲特性和濾波處理方法進行研究。

1 濾波方法的設計

為了避免引力梯度數據濾波后失真，將濾波對象定為觀測值減去模擬值的剩余部分OMC(observation minus computation)。模擬引力梯度數據所用的模型為EGM2008的前300階次，模擬引力梯度數據的Laplace算子的絕對值達到了10－21數量級。以2009-11-02的引力梯度數據為例，4個精度較高的梯度分量的OMC時間序列見圖1。

根據引力梯度數據的頻譜特性，本文設計了3種引力梯度數據濾波方法。

方法1:利用FFT技術把時域下的引力梯度數據轉換到頻域，將頻帶外的各頻率點賦0值，觀測頻帶為5 mHz～0.1 Hz，并利用FFT技術將頻域下的序列值轉換到時域。

方法2:直接采用FIR(finite impulse response)濾波器對引力梯度序列進行處理。帶通濾波器的階數設定為1 000階，通帶設定為5 mHz～0.1 Hz，窗函數為漢寧窗。為了避免濾波結果出現相位漂移現象，將引力梯度序列分別按正向和反向通過FIR濾波器，即

式中，h(n)為帶通濾波器，x0(n)為原始梯度數據序列，x1(n)為正向濾波后的梯度序列，x2(n)為反向的梯度數據序列，x3(n)為反向濾波后的梯度數據序列，x4(n)為經過正反向濾波后的梯度數據序列。

方法3:在傅里葉多項式的基礎上，以每個頻率點處的相位值作為未知參數，建立引力梯度觀測值關于各頻率點相位的觀測方程，利用最小二乘法解算出各頻率點的相位值，然后利用各頻點的相位值計算出各歷元的低頻誤差，在觀測值中扣除各歷元的低頻誤差。利用帶通濾波器或低通濾波器處理高頻誤差，即

圖1 OMC的時間序列Fig.1 The time series of OMC data

式中，s(t)為時刻 t處的引力梯度值，A0、Aic和 Ais為傅里葉多項式的系數，fi為第i點的頻率，fi的最大值為低頻帶的上限0.005 Hz。利用式(2)通過最小二乘法即可得到低頻誤差。

2 數據計算與分析

利用3種濾波方法對2009-11-02～2010-01-10的引力梯度數據進行濾波處理，方法1濾波后的功率譜密度見圖2。

圖2 方法1濾波后的結果Fig.2 The filtered results with method 1

從圖2可見，引力梯度在頻帶(5 mHz～0.1 Hz)內的信號被提取出來，頻帶外的梯度信號被全部去掉。由于頻帶外的信號為0，在FFT逆變換時出現了混頻現象，并導致頻帶內梯度信號的譜密度降低。這種濾波方法降低了引力梯度在頻帶內的精度，最終將影響重力場模型的恢復。

方法2濾波后的功率譜密度見圖3。從圖3可見，FIR帶通濾波器能夠完整地提取出頻帶內的梯度信號，并且保留頻帶外的引力梯度信號，濾波效果較好。

方法3濾波后的功率譜密度見圖4。圖4中，OMCxx0、OMCyy0、OMCzz0 和 OMCxz0 為減去低頻誤差 后 的 結 果;OMCxx1、OMCyy1、OMCzz1 和OMCxz1為最終的濾波結果。可見，濾波后的結果能較好地提取出觀測頻帶內的信號。

在引力梯度張量中，Laplace算子是評價引力梯度質量的重要指標。為了更好地評估3種方法，計算了3種方法濾波后梯度的Laplace算子的絕對值(表1)。從表1可見，方法1的平均值和標準差都為最小，但是個別值的精度極差，大于100×10－12s－2，遠低于梯度測量儀在頻帶內的精度。由此可見，方法1出現的混頻現象使濾波后的引力梯度數據失真。方法2和方法3濾波后的結果精度接近，但總體上看，方法3的精度最高。從實測的引力梯度數據的濾波結果可知，實測的引力梯度數據精度略低于標稱精度［5］。

表1 3種濾波結果的統計(單位:10－12s－2)Tab.1 Statistics of the results from three filtered methods(unit:10 －12s－2)

圖3 方法2濾波后的結果Fig.3 The filtered results with method 2

采用第3種方法對2009-11-02～2010-01-10共70 d的引力梯度數據進行濾波，并以濾波后的引力梯度數據為基礎利用時域法恢復200階次的地球重力場模型(參與計算的數據包括精度較高的4個分量，精度較差的兩個分量不參與反演計算)。由于GOCE衛星的太陽同步軌道特性，導致觀測數據無法完全覆蓋地球，在地球兩極分別有半徑約為6.5°的空白區域［6－7］。兩極空白將會導致法方程嚴重病態，解算的模型的帶諧位系數附近會形成一個嚴重的低精度帶(圖5(a))。為了抑制法方程的病態，對法方程進行Kaula正則化，最優正則化參數α=10－12.3，選取的原則是以反演模型的大地水準面誤差最小為標準。施加了Kaula正則化后恢復的模型能有效體現出引力梯度數據在中短波信息方面的優勢(圖5(b))。

從圖5(a)可見，無正則化的模型在帶諧位系數附近誤差很大，而施加了Kaula正則化之后帶諧位系數精度得到了明顯改善(圖5(b))，但個別低階次位系數精度仍然偏低。從梯度數據反演的模型可見，梯度數據的優勢在重力場模型的高階次部分。

圖4 方法3的濾波結果Fig.4 The filtered results with method 3

為了比較恢復的重力場模型精度，計算了3個模型的階方差，包括利用70 d梯度數據無正則化反演的模型Model1，利用70 d梯度數據施加Kaula正則化反演的模型Model2，ICGEM公布的利用近18個月的GOCE觀測數據反演的模型GO_CONS_GCF_2_TIM_R3，用Model3表示。階方差比較見圖6。

圖5 無正則化和Kaula正則化恢復的模型Fig.5 Gravity field model without regularization and with Kaula regularization

圖6 3個模型的階方差比較Fig.6 Comparison of degreevariance of three models

從圖6可見，無正則化的模型Model1不可用;施加了Kaula正則化的模型Model2在160階以前與Model3具有相同的精度，在160階以后，Model2的精度低于Model3。在160階以后精度略低的原因可能是恢復 Model2僅用了70 d數據，而恢復Model3用了近18個月的數據。從恢復的模型可知，本文所用的引力梯度數據能有效地濾掉低頻誤差，濾波后的引力梯度數據能滿足恢復重力場模型的需要。

3 結論

1)從3種濾波方法的濾波結果可知，GOCE實測引力梯度數據的精度并未達到發射前設計的標稱精度3.2×10－12s－2，略低于 ESA 的標稱精度。由于低頻誤差的量級較大，未經濾波的梯度數據不能用于恢復重力場模型。

2)利用濾波方法3處理的70 d引力梯度數據恢復了1個200階次的重力場模型，施加正則化后能體現出GOCE引力梯度數據在中短波方面的優勢。施加了正則化的模型在160階以前與ICGEM公布的GOCE時域法模型精度一致;160階以后本文恢復的模型精度略低，原因可能是本文參與計算的數據量較少。

3)從恢復的模型可見，本文所采用的引力梯度數據濾波方法可用于恢復地球重力場模型。GOCE衛星的遠期目標是恢復高精度的靜態重力場模型，本文的工作為進一步聯合GOCE軌道和引力梯度數據恢復高精度重力場模型作好了準備。

1 Pail R，et al.First GOCE gravity field models derived by three different approaches［J］.Journal of Geodesy，2011，85(11):819－843.

2 萬曉云，于錦海，曾艷艷.GOCE引力梯度的頻譜分析及濾波［J］.地球物理學報，2012，55(9):2 909 －2 916.(Wan Xiaoyun，Yu Jinhai，Zeng Yanyan.Frequency analysis and filtering processing of gravity gradients data from GOCE［J］.Chinese J Geophys，2012，55(9):2 909 －2 916)

3 Floberghagen R，et al.Mission design，operation and exploitation of the gravity field and steady-state ocean circulation explorer mission［J］.Journal of Geodesy，2011，85(11):749 －758.

4 徐新禹，等.GOCE衛星重力梯度測量誤差分析及其模擬研究［J］.大地測量與地球動力學，2010(2):71 －75.(Xu Xinyu，et al.Research on analysis and simulation of gravity gradiometry error of GOCE satellite［J］.Journal of Geodesy and Geodynamic，2010(2):71 －75)

5 鐘波.基于GOCE衛星重力測量技術確定地球重力場的研究［D］.武漢:武漢大學，2010.(Zhong Bo.Determination of earth’s gravitational field based on the GOCE gravity measurement technigue［D］.Wuhan:Wuhan Univercity，2010)

6 Metzler B，Pail R.GOCE data processing:the spherical cap regularization approach［J］.Studia Geophysica et Geodaetica，2005(4):441 －462.

7 Ditmar P，Kusche J，Klees R.Computation of spherical harmonic coefficients from gravity gradiometry data to be acquired by the GOCE satellite:regularization issues［J］.Journal of Geodesy，2003，77(7 －8):465－477.

DESIGN OF FILTERING METHOD OF GOCE GRAVITY GRADIENT AND ANALYSIS OF RESULT

Yang Hongguo1)and Huang Qiang2)
1)Qingdao Agriculture University，Qingdao 266109
2)Faculty of Geoscience and Environment Engineering，Southwest Jiaotong University，Chengdu610097

For the reason that GOCE gravity gradients is impacted badly by low frequency error，three filtered methods were designed based on the spectral characteristics of GOCE gravity gradients.The best method was choosed through comparison of the results with the three methods.The method was applied to filter 70 days gravity gradients between 2009-11-02—2010-01-10.An Earth’s gravity field model with 200 degree and order was reconstructed based on 70 days filtered gravity gradients.The results show that reconstruction of model have to be processed with regularization.The precision of the regularized model is consistant with that of GO_CONS_GCF_2_TIM_R3 with 160 degrees.The results prove that the designed method can filter gravity gradients successfully.

GOCE;gravity gradients;low frequency error;fast Fourier transform;band pass filter

P223.0

A

1671-5942(2014)03-0146-05

2013-08-23

中央高校基本科研業務費專項(20120184120006)。

楊洪國，男，1976年生，工程師，主要從事大地測量學與測量工程研究。E－mail:874143839@.qq.com。