GOCE引力梯度數(shù)據(jù)濾波方法的設(shè)計與結(jié)果分析*

楊洪國 黃 強

1)青島農(nóng)業(yè)大學(xué)，青島 266109

2)西南交通大學(xué)地球科學(xué)與環(huán)境工程學(xué)院，成都 610097

GOCE引力梯度數(shù)據(jù)濾波方法的設(shè)計與結(jié)果分析*

楊洪國1)黃 強2)

1)青島農(nóng)業(yè)大學(xué)，青島 266109

2)西南交通大學(xué)地球科學(xué)與環(huán)境工程學(xué)院，成都 610097

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

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

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

1 濾波方法的設(shè)計

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

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

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

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

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

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

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

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

2 數(shù)據(jù)計算與分析

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

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

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

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

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

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

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

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

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

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

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

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

圖5 無正則化和Kaula正則化恢復(fù)的模型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階以后精度略低的原因可能是恢復(fù) Model2僅用了70 d數(shù)據(jù)，而恢復(fù)Model3用了近18個月的數(shù)據(jù)。從恢復(fù)的模型可知，本文所用的引力梯度數(shù)據(jù)能有效地濾掉低頻誤差，濾波后的引力梯度數(shù)據(jù)能滿足恢復(fù)重力場模型的需要。

3 結(jié)論

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

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

3)從恢復(fù)的模型可見，本文所采用的引力梯度數(shù)據(jù)濾波方法可用于恢復(fù)地球重力場模型。GOCE衛(wèi)星的遠期目標(biāo)是恢復(fù)高精度的靜態(tài)重力場模型，本文的工作為進一步聯(lián)合GOCE軌道和引力梯度數(shù)據(jù)恢復(fù)高精度重力場模型作好了準(zhǔn)備。

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］.地球物理學(xué)報，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衛(wèi)星重力梯度測量誤差分析及其模擬研究［J］.大地測量與地球動力學(xué)，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衛(wèi)星重力測量技術(shù)確定地球重力場的研究［D］.武漢:武漢大學(xué)，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

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

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