楊洪國 黃 強(qiáng)

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

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

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

楊洪國1)黃 強(qiáng)2)

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

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

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

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

GOCE衛(wèi)星的目的是恢復(fù)高精度靜態(tài)地球重力場(chǎng)模型。ESA預(yù)期可恢復(fù)大地水準(zhǔn)面精度為1 cm、重力異常精度為10－5ms－2、空間分辨率為100 km(半波長)的地球重力場(chǎng)模型［1－3］。由于GOCE衛(wèi)星自身的技術(shù)特性，導(dǎo)致6個(gè)梯度分量中的4個(gè)(Vxx，Vyy，Vzz，Vxz)精度較高，另外兩個(gè)(Vxy，Vyz)精度較差［4］。此外，引力梯度分量含有巨大的低頻誤差，低頻誤差的濾波處理直接影響到恢復(fù)重力場(chǎng)模型的計(jì)算。目前，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_在處理時(shí)引入了大量的外部模型信息，并在坐標(biāo)系之間進(jìn)行了轉(zhuǎn)換，受到姿態(tài)誤差的影響，不能用于恢復(fù)重力場(chǎng)模型。EGG_NOM_2_可用于恢復(fù)重力場(chǎng)模型的計(jì)算，但必須進(jìn)行濾波處理。本文主要對(duì)引力梯度數(shù)據(jù)的噪聲特性和濾波處理方法進(jìn)行研究。

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

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

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

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

方法2:直接采用FIR(finite impulse response)濾波器對(duì)引力梯度序列進(jìn)行處理。帶通濾波器的階數(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)為正向?yàn)V波后的梯度序列，x2(n)為反向的梯度數(shù)據(jù)序列，x3(n)為反向?yàn)V波后的梯度數(shù)據(jù)序列，x4(n)為經(jīng)過正反向?yàn)V波后的梯度數(shù)據(jù)序列。

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

圖1 OMC的時(shí)間序列Fig.1 The time series of OMC data

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

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

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

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

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

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

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

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

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

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

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

為了比較恢復(fù)的重力場(chǎng)模型精度，計(jì)算了3個(gè)模型的階方差，包括利用70 d梯度數(shù)據(jù)無正則化反演的模型Model1，利用70 d梯度數(shù)據(jù)施加Kaula正則化反演的模型Model2，ICGEM公布的利用近18個(gè)月的GOCE觀測(cè)數(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個(gè)模型的階方差比較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個(gè)月的數(shù)據(jù)。從恢復(fù)的模型可知，本文所用的引力梯度數(shù)據(jù)能有效地濾掉低頻誤差，濾波后的引力梯度數(shù)據(jù)能滿足恢復(fù)重力場(chǎng)模型的需要。

3 結(jié)論

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

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

3)從恢復(fù)的模型可見，本文所采用的引力梯度數(shù)據(jù)濾波方法可用于恢復(fù)地球重力場(chǎng)模型。GOCE衛(wèi)星的遠(yuǎn)期目標(biāo)是恢復(fù)高精度的靜態(tài)重力場(chǎng)模型，本文的工作為進(jìn)一步聯(lián)合GOCE軌道和引力梯度數(shù)據(jù)恢復(fù)高精度重力場(chǎng)模型作好了準(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é)報(bào)，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)星重力梯度測(cè)量誤差分析及其模擬研究［J］.大地測(cè)量與地球動(dòng)力學(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)星重力測(cè)量技術(shù)確定地球重力場(chǎng)的研究［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

中央高?；究蒲袠I(yè)務(wù)費(fèi)專項(xiàng)(20120184120006)。

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