潘娟霞，鄒賢才,2

1. 武漢大學(xué)測(cè)繪學(xué)院，湖北 武漢 430079； 2. 武漢大學(xué)地球空間環(huán)境與大地測(cè)量教育部重點(diǎn)實(shí)驗(yàn)室，湖北 武漢 430079

重力測(cè)量衛(wèi)星GOCE(gravity field and steady-state ocean circulation explorer)是由歐空局(European Space Agency，ESA)研制和發(fā)射的低軌重力探測(cè)衛(wèi)星，其科學(xué)目標(biāo)是在100 km的空間分辨率內(nèi)(即恢復(fù)200階次以上的全球重力場(chǎng)模型)以2 cm的精度測(cè)定全球大地水準(zhǔn)面及10-5m/s2的精度測(cè)定重力異常。相比CHAMP(challenging minisatellite payload)和GRACE(gravity recovery and climate experiment)衛(wèi)星，GOCE關(guān)鍵載荷是一臺(tái)高精度的引力梯度儀，由6個(gè)加速度計(jì)對(duì)稱安置在3個(gè)正交軸上，3對(duì)加速度計(jì)基線長(zhǎng)約50 cm，衛(wèi)星質(zhì)心與梯度儀質(zhì)心重合，從而利用加速度計(jì)差分觀測(cè)值直接測(cè)定引力位的二階梯度張量，并以加速度計(jì)共模觀測(cè)值獲得作用于衛(wèi)星的非保守力[1]。GOCE還有精確的高低衛(wèi)星跟蹤數(shù)據(jù)及無(wú)阻尼控制，聯(lián)合引力梯度數(shù)據(jù)及高低衛(wèi)星跟蹤數(shù)據(jù)解算重力場(chǎng)是其獨(dú)特優(yōu)勢(shì)。

GOCE引力梯度儀受到測(cè)量誤差、外界觀測(cè)條件等因素影響，導(dǎo)致觀測(cè)值出現(xiàn)系統(tǒng)偏差、比例誤差及有色噪聲等，因此，反演重力場(chǎng)前針對(duì)GOCE觀測(cè)數(shù)據(jù)的校準(zhǔn)至關(guān)重要。安裝結(jié)構(gòu)不同導(dǎo)致GRACE衛(wèi)星的加速度計(jì)校準(zhǔn)方式不再適用于GOCE。國(guó)內(nèi)外諸多學(xué)者針對(duì)GOCE數(shù)據(jù)的使用及校準(zhǔn)做了大量研究，通常根據(jù)是否引入?yún)⒖贾亓?chǎng)模型等外部輔助數(shù)據(jù)將GOCE的校準(zhǔn)分為內(nèi)部校準(zhǔn)及外部校準(zhǔn)，校準(zhǔn)模型通常包括比例因子、偏差因子。文獻(xiàn)[2—3]提出用恒星敏感器測(cè)量數(shù)據(jù)及先驗(yàn)重力場(chǎng)估計(jì)比例因子，該模型還校準(zhǔn)了梯度儀坐標(biāo)系(gradiometer reference frame，GRF)與恒星敏感器坐標(biāo)系(star sensor reference frame，SSRF)間未配準(zhǔn)的問(wèn)題。文獻(xiàn)[4]提出了將恒星敏感器聯(lián)合梯度儀測(cè)量數(shù)據(jù)重建的衛(wèi)星角速度用于校準(zhǔn)模型中。文獻(xiàn)[5—6]提出類似方法用恒星敏感器數(shù)據(jù)及加速度計(jì)測(cè)量值進(jìn)行內(nèi)部校準(zhǔn)，并針對(duì)觀測(cè)值時(shí)間相關(guān)性提出解決方案，該方法即為ESA官方采用方法。文獻(xiàn)[7—8]根據(jù)精密軌道確定比例因子與偏差因子，利用恒星敏感器數(shù)據(jù)對(duì)校準(zhǔn)后的引力梯度數(shù)據(jù)進(jìn)行驗(yàn)證。文獻(xiàn)[9—10]發(fā)現(xiàn)在校準(zhǔn)模型中加入二次因子可以減弱地磁極周圍強(qiáng)風(fēng)的影響，于2018年引入外部重力場(chǎng)模型對(duì)GOCE數(shù)據(jù)重新處理。國(guó)內(nèi)諸多研究包括對(duì)GOCE衛(wèi)星數(shù)據(jù)的預(yù)處理研究、梯度數(shù)據(jù)確定地球重力場(chǎng)的理論方法研究[11-13]及基于加速度計(jì)數(shù)據(jù)的校準(zhǔn)方法研究，文獻(xiàn)[14—16]提出聯(lián)合幾何法精密軌道以動(dòng)力法完成單加速度計(jì)校準(zhǔn)及衛(wèi)星非保守力確定，同時(shí)解算重力場(chǎng)模型和校準(zhǔn)參數(shù)，從而降低參考重力場(chǎng)模型誤差對(duì)校準(zhǔn)結(jié)果的影響，并討論了衛(wèi)星無(wú)阻尼控制的補(bǔ)償效果。文獻(xiàn)[17]利用外部重力場(chǎng)及恒星敏感器數(shù)據(jù)初步驗(yàn)證對(duì)一定頻段內(nèi)校準(zhǔn)的有效性。文獻(xiàn)[18]討論利用不同外部重力場(chǎng)模型校準(zhǔn)及相同重力場(chǎng)模型不同階次對(duì)校準(zhǔn)結(jié)果的影響。

目前，國(guó)內(nèi)針對(duì)GOCE數(shù)據(jù)的處理及使用以ESA內(nèi)部校準(zhǔn)后的數(shù)據(jù)為主，本文以ESA官方發(fā)布的方法為基礎(chǔ)，實(shí)現(xiàn)GOCE數(shù)據(jù)的內(nèi)部校準(zhǔn)，主要目的是為完善GOCE梯度儀觀測(cè)數(shù)據(jù)的校準(zhǔn)處理，并為討論內(nèi)部校準(zhǔn)、外部校準(zhǔn)方法及動(dòng)力法校準(zhǔn)的聯(lián)合、比較做初步準(zhǔn)備。ESA發(fā)布的內(nèi)部校準(zhǔn)方法將誤差因素模型化為3個(gè)加速度計(jì)對(duì)的逆校準(zhǔn)矩陣，并將該矩陣作用于加速度計(jì)實(shí)現(xiàn)校準(zhǔn)，本文在此基礎(chǔ)上提出了顧及SSRF與GRF間的旋轉(zhuǎn)矩陣的校準(zhǔn)參數(shù)的改進(jìn)模型，討論這一變化對(duì)衛(wèi)星引力梯度精度帶來(lái)的改變。涉及的數(shù)據(jù)包括衛(wèi)星搭載的梯度儀及恒星敏感器測(cè)量數(shù)據(jù)，利用加速度計(jì)共模觀測(cè)值及差分觀測(cè)值與恒星敏感器的姿態(tài)數(shù)據(jù)，根據(jù)衛(wèi)星引力梯度測(cè)量原理建立校準(zhǔn)模型，以獲取校準(zhǔn)后梯度儀坐標(biāo)系下的引力梯度。除此以外，恢復(fù)重力場(chǎng)還涉及梯度儀坐標(biāo)系與慣性系的高精度轉(zhuǎn)換，因此本文提出了關(guān)于衛(wèi)星姿態(tài)重建的討論。

1 理論與方法

1.1 內(nèi)部校準(zhǔn)

GOCE衛(wèi)星上搭載了幾類重要載荷，本文重點(diǎn)關(guān)注用于確定中短波重力場(chǎng)的靜電引力梯度儀(electrostatic gravity gradiometer，EGG)，以及提供衛(wèi)星姿態(tài)的3個(gè)恒星敏感器。GOCE是第一顆采用無(wú)阻尼控制技術(shù)的衛(wèi)星，沿軌方向大氣阻力得到持續(xù)補(bǔ)償。衛(wèi)星實(shí)際運(yùn)行過(guò)程中，存在諸多誤差因素，為此GOCE設(shè)計(jì)了每?jī)蓚€(gè)月執(zhí)行一次特定的持續(xù)1 d的振動(dòng)，用于衛(wèi)星的校準(zhǔn)，該振動(dòng)期間的數(shù)據(jù)稱為shaking數(shù)據(jù)，其他時(shí)段稱為nominal數(shù)據(jù)。圖1為梯度儀坐標(biāo)系中3對(duì)加速度計(jì)安置結(jié)構(gòu)，實(shí)線代表超敏感軸，虛線代表非敏感軸[6]。

圖1 梯度儀坐標(biāo)系中3對(duì)加速度計(jì)安置結(jié)構(gòu)Fig.1 Arrangement of the three pairs of accelerometers in the GRF

1.1.1 共模加速度與差分加速度

圖2是梯度儀中3對(duì)加速度計(jì)的特殊結(jié)構(gòu)，根據(jù)衛(wèi)星引力梯度測(cè)量原理，第i個(gè)加速度計(jì)測(cè)得的加速度值為

(1)

梯度儀兩類觀測(cè)模式為差分模式加速度(differential mode acceleration，DM)與共模加速度(common mode acceleration，CM)[6]

(2)

(3)

ac,ij≈d

(4)

(5)

式(5)表達(dá)為矩陣形式

(6)

式中，Ad=(ad,14,ad,25,ad,36)，而矩陣L[6]

(7)

(8)

AdL-1+(AdL-1)T=-V+Ω2

(9)

1.1.2 校準(zhǔn)模型

衛(wèi)星實(shí)際運(yùn)行過(guò)程中，通?？紤]以下誤差因素：①加速度計(jì)質(zhì)心偏離標(biāo)稱位置；②加速度計(jì)各坐標(biāo)軸未嚴(yán)格與梯度儀相應(yīng)坐標(biāo)軸對(duì)齊；③加速度計(jì)3軸不完全正交；④由于數(shù)據(jù)輸出增益的不確定性，產(chǎn)生的加速度計(jì)比例因子[10]。將誤差參數(shù)化為比例因子、偏差因子及二次項(xiàng)，其中二次項(xiàng)以物理振動(dòng)的方式予以消除[19]，因此內(nèi)部校準(zhǔn)矩陣應(yīng)包括比例因子與偏差因子。圖3表示除非線性二次項(xiàng)以外的加速度計(jì)誤差因素，每個(gè)加速度計(jì)包括6個(gè)角度校準(zhǔn)參數(shù)、3個(gè)比例因子，6個(gè)加速度計(jì)共需要確定54個(gè)校準(zhǔn)參數(shù)，定義校準(zhǔn)矩陣[6]

(10)

圖2 單個(gè)加速度計(jì)的誤差來(lái)源Fig.2 Error sources of a single accelerometer

(11)

(12)

式(12)給出了測(cè)量值、真實(shí)值及逆校準(zhǔn)矩陣間的關(guān)系。

式(4)表示由3對(duì)加速度計(jì)觀測(cè)值給出的共模加速度CM即是衛(wèi)星在3個(gè)方向所受到的非保守力，有以下6個(gè)獨(dú)立條件

E{ac,ij-ac,kl}=0ij=14,25,36,ij≠kl

(13)

(14)

類似地，聯(lián)合恒星敏感器角速度ΩS，式(9)可變換為

(15)

圖3 0.05～0.1 Hz內(nèi)離心加速度、差分加速度、引力梯度PSD1/2比較Fig.3 Comparison of PSD1/2 of centrifugal acceleration, differential acceleration and gravity gradients in the 0.05～0.1 Hz

(16)

根據(jù)式(13)、式(14)與式(16)，以ESA官方發(fā)布的衛(wèi)星nominal時(shí)段的梯度儀測(cè)量數(shù)據(jù)EGG_NOM_L1b及恒星敏感器數(shù)據(jù)STR_VC2_1b與STR_VC3_1b作為輸入數(shù)據(jù)以實(shí)現(xiàn)內(nèi)部校準(zhǔn)。

建立以上模型時(shí)，做了以下假定：①ESA給定的引力梯度儀基線長(zhǎng)Lx、Ly、Lz是準(zhǔn)確且為固定常數(shù)；②ESA給定的SSRF與GRF間的相對(duì)關(guān)系是準(zhǔn)確且為常數(shù)矩陣；本文探索考慮這兩個(gè)假定是否準(zhǔn)確以及其對(duì)引力梯度精度的影響，在基礎(chǔ)模型上加上新的參數(shù)重新建立模型。用微小旋轉(zhuǎn)角rx、ry、rz描述恒星敏感器與梯度儀坐標(biāo)系間旋轉(zhuǎn)矩陣隨時(shí)間的變化對(duì)角加速度的綜合影響，用GRF′表示標(biāo)稱梯度儀坐標(biāo)系，GRF表示存在偏差的坐標(biāo)系，兩坐標(biāo)系中角加速度的關(guān)系表達(dá)為

(17)

(18)

(19)

寫成分量形式為

(20)

由于式(17)在0.05～0.1 Hz頻段內(nèi)作了近似處理，且分母上省略了基線長(zhǎng)Lx、Ly、Lz，導(dǎo)致基線長(zhǎng)的估計(jì)少了完整約束，因此本文暫不討論基線長(zhǎng)隨時(shí)間的變化，仍然認(rèn)為其數(shù)值為常數(shù)。本文在ESA方法的基礎(chǔ)上，利用式(13)、式(20)及式(16)完成內(nèi)部校準(zhǔn)。

1.1.3 協(xié)方差矩陣處理

本文采用2 d的數(shù)據(jù)估計(jì)一組校準(zhǔn)參數(shù)，根據(jù)GOCE采樣率及有色噪聲的特點(diǎn)，難以實(shí)現(xiàn)矩陣的存儲(chǔ)與求逆運(yùn)算，參考文獻(xiàn)[6,10]，利用對(duì)稱的滑動(dòng)平均去相關(guān)濾波，實(shí)現(xiàn)去相關(guān)處理的同時(shí)實(shí)現(xiàn)相應(yīng)頻段濾波，避免大型矩陣的求逆，式(21)為去相關(guān)濾波表達(dá)式

(21)

1.2 恒星敏感器的聯(lián)合

圖4 梯度儀坐標(biāo)系與3個(gè)星敏感器坐標(biāo)系的相對(duì)關(guān)系[6]Fig.4 Relative orientation of the GRF and the SSRFs

1.3 角速度與姿態(tài)重建

角速度與姿態(tài)四元數(shù)的精度直接影響到引力梯度與地球重力場(chǎng)反演精度，對(duì)角速度及姿態(tài)的重建是必要的。恒星敏感器觀測(cè)的四元數(shù)變換到角速度需經(jīng)過(guò)數(shù)值微分的過(guò)程，導(dǎo)致高頻段噪聲被放大，加速度計(jì)觀測(cè)值經(jīng)積分導(dǎo)出的角速度及四元數(shù)則放大了低頻噪聲。文獻(xiàn)[19]采用卡爾曼濾波實(shí)現(xiàn)角速度與姿態(tài)重建。文獻(xiàn)[20]利用恒星敏感器角速度與梯度儀角速度建立頻域內(nèi)的權(quán)重模型，以維納濾波重建角速度與姿態(tài)。文獻(xiàn)[10]提出采用最小二乘擬合以重建姿態(tài)?？柭鼮V波與最小二乘擬合均在時(shí)域內(nèi)展開(kāi)而未使用頻域特點(diǎn)，且卡爾曼濾波的瞬態(tài)效應(yīng)嚴(yán)重造成數(shù)據(jù)利用率不高，本文參考[20]利用維納濾波在頻域內(nèi)聯(lián)合兩類數(shù)據(jù)的優(yōu)點(diǎn)實(shí)現(xiàn)角速度及姿態(tài)重建。

維納濾波的原理是根據(jù)兩類觀測(cè)值精度實(shí)現(xiàn)頻域內(nèi)的加權(quán)平均，頻域內(nèi)某頻率的功率譜密度值(power spectral density,PSD)代表該處的精度

(22)

式中，HSTR(f)與HGRAD(f)分別代表敏感器與梯度儀的權(quán)重；PSTR(f)、PGRAD(f)分別表示頻域內(nèi)兩類角速度的平方根功率譜密度，模型參考文獻(xiàn)[20]。圖5為兩類角速度噪聲PSD1/2，其中，G代表梯度儀，S代表恒星敏感器。

圖5 兩類角速度噪聲PSD1/2Fig.5 PSD1/2 of two types of noise

角速度最終通過(guò)頻域內(nèi)乘積與傅里葉正反變換計(jì)算得到

(23)

(24)

式中，t代表觀測(cè)時(shí)刻。

2 數(shù)據(jù)分析

圖6 部分參數(shù)序列Fig.6 Part parameters series of

圖7 校準(zhǔn)后引力梯度張量跡的PSD1/2Fig.7 PSD1/2 of calibrated gravity gradient trace

表1 加速度計(jì)對(duì)比例因子Tab.1 Scale factors of accelerometer pairs

圖8 恒星敏感器聯(lián)合前后引力梯度張量跡的PSD1/2Fig.8 PSD1/2 of calibrated gravity gradient trace before and after combination of star sensors

圖9 恒星敏感器聯(lián)合前后對(duì)應(yīng)引力梯度分量與模型參考值差值的PSD1/2Fig.9 PSD1/2 of differences of gravity gradients to EIGEN-5C model before and after combination of star sensors

在第1類ESA模型基礎(chǔ)上，本文增加了SSRF與GRF之間旋轉(zhuǎn)矩陣隨時(shí)間變化的參數(shù)ΔR，這一變化將對(duì)角速度及姿態(tài)重建過(guò)程產(chǎn)生影響，引力梯度分量也將產(chǎn)生相應(yīng)改變。這里給出新增參數(shù)ΔR的數(shù)據(jù)序列及線性擬合序列，圖10中存在幾個(gè)明顯異常值，分析相應(yīng)時(shí)間段內(nèi)3個(gè)恒星敏感器工作狀態(tài)，統(tǒng)計(jì)發(fā)現(xiàn)在異常值存在的時(shí)段內(nèi)，僅有單個(gè)恒星敏感器數(shù)據(jù)可用的時(shí)刻數(shù)遠(yuǎn)大于其他時(shí)段相應(yīng)值，這也體現(xiàn)了聯(lián)合多個(gè)恒星敏感器的優(yōu)勢(shì)和必要性。rx、ry、rz數(shù)值絕對(duì)值在100″左右且該月數(shù)值呈現(xiàn)約3～30″的線性趨勢(shì)，因而建立校準(zhǔn)模型時(shí)該校準(zhǔn)參數(shù)不應(yīng)被忽略。

圖10 參數(shù)ΔR序列Fig.10 Time series of ΔR

將ΔR作用到恒星敏感器角速度及四元數(shù)，再重建角速度與姿態(tài)以定量分析新增校準(zhǔn)參數(shù)對(duì)引力梯度精度的影響。圖11給出了加入?yún)?shù)ΔR前后引力梯度分量與參考模型引力梯度差值的平方根功率譜分析，表2為低于0.005 Hz頻率的兩類校準(zhǔn)模型對(duì)應(yīng)引力梯度分量與模型參考值差值的平方根功率譜統(tǒng)計(jì)結(jié)果，其中“方法1”表示與ESA模型一致，“方法2”代表新增ΔR的結(jié)果，4個(gè)引力梯度分量Vxx、Vxz、Vyy、Vzz對(duì)應(yīng)的差值平方根功率譜均顯示增加新的參數(shù)后引力梯度分量的精度有一定提升，其中分量Vxz精度提升最大，這是由于“方法2”中新增的參數(shù)改變了衛(wèi)星的姿態(tài)，相對(duì)其他分量而言，Vxz分量對(duì)姿態(tài)的變化更敏感，從而其精度提升最大。由于新增參數(shù)主要影響到姿態(tài)四元數(shù)，而梯度張量跡不會(huì)隨四元數(shù)改變[21-27]，這里未提供相應(yīng)對(duì)比圖。

圖11 兩類校準(zhǔn)模型對(duì)應(yīng)引力梯度分量與模型參考值差值的PSD1/2Fig.11 PSD1/2 of differences of gravity gradients to EIGEN-5C model for the two calibration models

表2 兩類引力梯度分量與模型參考值差值的PSD1/2統(tǒng)計(jì)Tab.2 PSD1/2 statistics of differences of gravity gradients to EIGEN-5C model for the two calibration models

3 結(jié)論與展望

本文從GOCE衛(wèi)星引力梯度測(cè)量原理著手，利用L1b數(shù)據(jù)中nominal時(shí)段的引力梯度儀觀測(cè)數(shù)據(jù)及恒星敏感器姿態(tài)數(shù)據(jù)對(duì)加速度計(jì)測(cè)量值做內(nèi)部校準(zhǔn)，以最小二乘聯(lián)合多個(gè)恒星敏感器觀測(cè)值以避免SSRF及GRF間轉(zhuǎn)換導(dǎo)致的誤差傳播；為得到慣性坐標(biāo)系下更為精確的引力梯度張量以恢復(fù)重力場(chǎng)，采用維納濾波重建了衛(wèi)星角速度及姿態(tài)四元數(shù)。在ESA校準(zhǔn)方法的基礎(chǔ)上，提出顧及SSRF與GRF間旋轉(zhuǎn)矩陣校準(zhǔn)參數(shù)的內(nèi)部校準(zhǔn)模型，結(jié)果表明該組校準(zhǔn)參數(shù)的數(shù)值絕對(duì)值在100″左右，且在該月呈現(xiàn)約3～30″的線性趨勢(shì)，同時(shí)解算逆校準(zhǔn)矩陣及SSRF與GRF間旋轉(zhuǎn)矩陣校準(zhǔn)參數(shù)的內(nèi)部校準(zhǔn)模型可改進(jìn)低于0.005 Hz頻段內(nèi)引力梯度各分量的精度，對(duì)于梯度分量Vxz改進(jìn)最大，證實(shí)了這一校準(zhǔn)參數(shù)存在的必要性，對(duì)GOCE自身數(shù)據(jù)處理及后續(xù)重力衛(wèi)星數(shù)據(jù)處理帶來(lái)一定的參考價(jià)值。

最后在本文的基礎(chǔ)上提出關(guān)于GOCE數(shù)據(jù)處理可能的改進(jìn)方向：①由于無(wú)外部數(shù)據(jù)的約束，內(nèi)部校準(zhǔn)方法不能避免GOCE自身系統(tǒng)偏差造成的影響，因此考慮比較外部校準(zhǔn)方法及動(dòng)力法以聯(lián)合各方法的優(yōu)勢(shì)，可推動(dòng)梯度數(shù)據(jù)處理問(wèn)題的討論；②關(guān)于基線長(zhǎng)Lx、Ly、Lz及ΔR的確定及二次項(xiàng)K2的補(bǔ)償，考慮引入先驗(yàn)重力場(chǎng)并聯(lián)合恒星敏感器數(shù)據(jù)解算；③在對(duì)協(xié)方差矩陣進(jìn)行對(duì)稱滑動(dòng)平均去相關(guān)濾波過(guò)程中，本文統(tǒng)一選擇濾波器階數(shù)等于軌道周期長(zhǎng)度，但實(shí)際解算時(shí)該長(zhǎng)度估計(jì)的功率譜并不是最優(yōu)，針對(duì)去相關(guān)過(guò)程中功率譜分辨率及精度的折中選擇值得深入討論。