辛寧邱樂(lè)德張立華劉乃金

（1中國(guó)空間技術(shù)研究院通信衛(wèi)星事業(yè)部，北京 100094）（2航天東方紅衛(wèi)星有限公司，北京 100094）

一種重力測(cè)量衛(wèi)星靜電加速度計(jì)在軌標(biāo)定算法

辛寧1邱樂(lè)德1張立華2劉乃金1

（1中國(guó)空間技術(shù)研究院通信衛(wèi)星事業(yè)部，北京 100094）（2航天東方紅衛(wèi)星有限公司，北京 100094）

靜電加速度計(jì)是低低跟蹤（SST－LL）重力測(cè)量衛(wèi)星的關(guān)鍵載荷之一，其性能直接影響地球重力場(chǎng)空間變化率的測(cè)定結(jié)果。為了確保靜電加速度計(jì)長(zhǎng)期在軌工作，結(jié)合擴(kuò)展卡爾曼濾波估計(jì)算法，提出了一種應(yīng)用動(dòng)力學(xué)方法確定靜電加速度計(jì)校準(zhǔn)參數(shù)的算法。首先建立靜電加速度計(jì)及K頻段測(cè)距（KBR）系統(tǒng)的量測(cè)模型；然后將高精度地球重力場(chǎng)模型和靜電加速度計(jì)觀測(cè)數(shù)據(jù)代入擴(kuò)展卡爾曼濾波算法的狀態(tài)方程中，將KBR系統(tǒng)觀測(cè)數(shù)據(jù)代入觀測(cè)方程中，建立靜電加速度計(jì)在軌標(biāo)定模型。數(shù)學(xué)仿真結(jié)果表明：靜電加速度計(jì)的標(biāo)度因子和零偏估計(jì)誤差均在0.2%以內(nèi)，實(shí)現(xiàn)了衛(wèi)星靜電加速度計(jì)較為精確的標(biāo)定。

重力測(cè)量衛(wèi)星；靜電加速度計(jì)；K頻段測(cè)距系統(tǒng)；在軌標(biāo)定

1 引言

靜電加速度計(jì)是低低跟蹤（SST－LL）重力測(cè)量衛(wèi)星的關(guān)鍵載荷，用于精確測(cè)量作用于衛(wèi)星上的非保守力加速度，其測(cè)量精度達(dá)到3×10－10m／s2。衛(wèi)星在軌運(yùn)行過(guò)程中，由于受空間環(huán)境輻射、各元件性能退化的影響，加速度計(jì)的零偏誤差、標(biāo)度因子具有不確定性［1］，直接影響重力場(chǎng)反演的精度。因此，為了反演高精度和高分辨率的地球重力場(chǎng)模型，提高重力測(cè)量衛(wèi)星的精密定軌精度，進(jìn)行加速度計(jì)的在軌標(biāo)定研究具有重要的意義。

文獻(xiàn)［2］中利用旋轉(zhuǎn)衛(wèi)星法開(kāi)展了加速度計(jì)的在軌標(biāo)定研究，對(duì)于評(píng)估加速度計(jì)在軌性能具有一定的指導(dǎo)意義，但是其中需要衛(wèi)星的姿態(tài)機(jī)動(dòng)、質(zhì)心調(diào)節(jié)等一系列復(fù)雜的操作，而且標(biāo)度因子檢驗(yàn)和坐標(biāo)軸偏差測(cè)量無(wú)法滿足預(yù)期精度。文獻(xiàn)［3］中利用加速度計(jì)在軌測(cè)量非保守力數(shù)據(jù)，并與標(biāo)準(zhǔn)的非保守力模型進(jìn)行比對(duì)給出標(biāo)定結(jié)果，但是非保守力模型（如大氣模型）往往很難準(zhǔn)確描述衛(wèi)星所受的非保守力，因此標(biāo)定結(jié)果精度較差。文獻(xiàn)［4］中將GPS精密定軌數(shù)據(jù)與根據(jù)重力場(chǎng)模型和加速度計(jì)實(shí)時(shí)測(cè)量數(shù)據(jù)計(jì)算出的軌道對(duì)比，利用Batch估計(jì)方法給出加速度計(jì)的校準(zhǔn)參數(shù)。該方法建立了GPS與加速度計(jì)之間的數(shù)據(jù)聯(lián)系，對(duì)于重力測(cè)量衛(wèi)星的總體設(shè)計(jì)具有一定的指導(dǎo)意義。但由于GPS精密定軌數(shù)據(jù)中并未包含衛(wèi)星的加速度信息，因此Batch估計(jì)方法中的系統(tǒng)狀態(tài)轉(zhuǎn)移矩陣非常復(fù)雜，同時(shí)采用積分的方法對(duì)軌道進(jìn)行估計(jì)會(huì)使觀測(cè)噪聲累積，且一般采用一天的數(shù)據(jù)估算一次校準(zhǔn)參數(shù)，使Batch估計(jì)方法的估計(jì)值可能無(wú)法收斂。文獻(xiàn)［5］中利用重力場(chǎng)反演方法（能量法）確定加速度計(jì)校準(zhǔn)參數(shù)，基本原理是利用現(xiàn)有的重力場(chǎng)模型計(jì)算衛(wèi)星高度處的擾動(dòng)位，與利用加速度計(jì)觀測(cè)值計(jì)算的擾動(dòng)位進(jìn)行比較，通過(guò)最小二乘法確定校準(zhǔn)參數(shù)，但該方法依賴于先驗(yàn)重力場(chǎng)模型與數(shù)據(jù)處理方法，根據(jù)不同的重力場(chǎng)模型計(jì)算出的加速度計(jì)標(biāo)度因子和零偏誤差并不相同。總體說(shuō)來(lái)，以上靜電加速度計(jì)的標(biāo)定方法均僅基于高低跟蹤（高軌GPS衛(wèi)星跟蹤低軌衛(wèi)星）的重力場(chǎng)測(cè)量模式，因此無(wú)法有效消除由非保守力模型、GPS定軌誤差、重力場(chǎng)模型等引起的系統(tǒng)誤差，加速度計(jì)標(biāo)定精度有待提高。

本文從重力測(cè)量衛(wèi)星工程設(shè)計(jì)出發(fā)，在文獻(xiàn)［4］的基礎(chǔ)上，利用低低跟蹤（低軌衛(wèi)星跟蹤低軌衛(wèi)星）的重力場(chǎng)測(cè)量模式，將重力測(cè)量衛(wèi)星另一個(gè)關(guān)鍵載荷K頻段測(cè)距（KBR）系統(tǒng)觀測(cè)到的星間加速度［6］代替雙星差分GPS數(shù)據(jù)作為系統(tǒng)觀測(cè)量，并利用高精度的地球重力場(chǎng)模型，基于動(dòng)力學(xué)方法建立靜電加速度計(jì)的量測(cè)模型，并利用擴(kuò)展卡爾曼濾波算法實(shí)現(xiàn)靜電加速度計(jì)參數(shù)的在軌標(biāo)定。由于星間加速度可由KBR系統(tǒng)直接輸出，因此本文算法無(wú)須進(jìn)行積分，可大大簡(jiǎn)化估計(jì)算法中的系統(tǒng)狀態(tài)轉(zhuǎn)移矩陣，建模簡(jiǎn)便，同時(shí)采用低低跟蹤的重力場(chǎng)測(cè)量模式進(jìn)行數(shù)據(jù)差分可有效消除GPS定軌誤差、先驗(yàn)重力場(chǎng)模型引起的系統(tǒng)誤差。此算法的可行性通過(guò)數(shù)學(xué)仿真進(jìn)行了驗(yàn)證，可為重力測(cè)量衛(wèi)星工程應(yīng)用提供參考。

2 靜電加速度計(jì)在軌標(biāo)定算法

靜電加速度計(jì)參數(shù)標(biāo)定算法的基本原理為：根據(jù)靜電加速度計(jì)和KBR系統(tǒng)的量測(cè)模型構(gòu)建靜電加速度計(jì)在軌標(biāo)定的觀測(cè)方程，再結(jié)合高精度重力場(chǎng)模型，利用擴(kuò)展卡爾曼濾波算法分離出具有線性特征的靜電加速度計(jì)誤差分量，由此來(lái)確定靜電加速度計(jì)的零偏誤差和標(biāo)度因子，從而實(shí)現(xiàn)靜電加速度計(jì)的在軌標(biāo)定。

2.1 靜電加速度計(jì)量測(cè)模型

靜電加速度計(jì)在動(dòng)態(tài)設(shè)計(jì)良好并進(jìn)入穩(wěn)態(tài)后，輸出值可表示為［7］

式中：靜電加速度計(jì)的標(biāo)度因子s＝diag（s1，s2，s3），s1，s2，s3為靜電加速度計(jì)標(biāo)度因子在3軸方向的分量；為在衛(wèi)星本體坐標(biāo)系下的非引力加速度，衛(wèi)星本體坐標(biāo)系的原點(diǎn)為衛(wèi)星質(zhì)心，x軸垂直于衛(wèi)星軌道平面，y軸為衛(wèi)星飛行方向，z軸為衛(wèi)星對(duì)地指向；b為零偏誤差；An為觀測(cè)高斯白噪聲。

由式（1）可得

式中：q＝［q0q1q2q3］T，為星敏感器觀測(cè)的姿態(tài)四元數(shù)。

根據(jù)地球重力場(chǎng)模型可確定衛(wèi)星所受的引力加速度f(wàn)g［5］，則衛(wèi)星加速度可表示為

2.2 KBR系統(tǒng)量測(cè)模型

KBR系統(tǒng)觀測(cè)值包括星間距離、星間速度及星間加速度。設(shè)衛(wèi)星A和衛(wèi)星B在地心慣性坐標(biāo)系中的質(zhì)心位置分別為rA和rB，則星間距離rAB的大小為式中：和為衛(wèi)星A和衛(wèi)星B的加速度；和為GPS數(shù)據(jù)精密定軌后得到的衛(wèi)星A和衛(wèi)星B的速度。

以上觀測(cè)值均可通過(guò)文獻(xiàn)［6］中構(gòu)建的KBR仿真系統(tǒng)中的衛(wèi)星軌道姿態(tài)模塊獲取。

2.3 靜電加速度計(jì)參數(shù)的擴(kuò)展卡爾曼濾波估計(jì)

靜電加速度計(jì)參數(shù)標(biāo)定模型的系統(tǒng)狀態(tài)變量為X＝［τAmAτBmB］T，則其系統(tǒng)狀態(tài)方程為

式中：τA，τB和mA，mB是計(jì)算得到的衛(wèi)星A和衛(wèi)星B的τ和m。

將式（2）和式（4）代入式（9），得系統(tǒng)觀測(cè)方程為

式中：H（t）為狀態(tài)轉(zhuǎn)移矩陣，如式（12）所示，t為積分時(shí)間；和分別為衛(wèi)星A和衛(wèi)星B根據(jù)地球重力場(chǎng)模型計(jì)算的引力加速度；Zn為測(cè)量誤差。

式中：MA和MB分別為衛(wèi)星A和衛(wèi)星B本體坐標(biāo)系到地心慣性坐標(biāo)系的轉(zhuǎn)換矩陣；NA和NB分別為衛(wèi)星A和衛(wèi)星B的靜電加速度計(jì)輸出。

給定初始狀態(tài)估計(jì)值＾X（t0）及初始估計(jì)均方差誤差P（t0），利用擴(kuò)展卡爾曼濾波算法即可估計(jì)出兩顆衛(wèi)星靜電加速度計(jì)的標(biāo)度因子和零偏誤差，其迭代過(guò)程如式（13）～（17）所示。

3 數(shù)學(xué)仿真

靜電加速度計(jì)在軌標(biāo)定算法的仿真流程見(jiàn)圖1。

圖1 靜電加速度計(jì)在軌標(biāo)定算法仿真過(guò)程Fig.1 Simulation procedure of electrostatic accelerometer on－orbit calibration algorithm

在衛(wèi)星軌道和姿態(tài)模擬器中［8］，地球引力場(chǎng)模型采用120階EMG96模型；大氣阻力模型采用大氣阻力溫度模型（Drag Temperature Model，DTM），大氣阻力系數(shù)為2.2；在太陽(yáng)光壓攝動(dòng)模型中，太陽(yáng)光壓反照系數(shù)為0.5；地球磁場(chǎng)模型采用13階IGRF2005模型。根據(jù)衛(wèi)星軌道和姿態(tài)模擬器輸出的衛(wèi)星位置、速度和姿態(tài)四元數(shù)，模擬生成KBR系統(tǒng)觀測(cè)值、靜電加速度計(jì)觀測(cè)值、GPS觀測(cè)值及星敏感器觀測(cè)值，將這些觀測(cè)值代入到擴(kuò)展卡爾曼濾波算法中，得到靜電加速度計(jì)的在軌標(biāo)定結(jié)果。

衛(wèi)星的飛行姿態(tài)為三軸對(duì)地定向模式，衛(wèi)星A和衛(wèi)星B初始軌道參數(shù)見(jiàn)表1。

表1 兩顆衛(wèi)星的初始軌道參數(shù)Table 1 Initial orbit parameters of two satellites

KBR系統(tǒng)仿真值如圖2所示，由于軌道周期為5600s，并且軌道偏心率不為零，因此星間距離序列表現(xiàn)為周期為5600s、振幅為2km的余弦波形。星間加速度值如圖3所示，其表現(xiàn)為周期為5600s、振幅為0.002 5m／s2的余弦波形。根據(jù)非保守力模型，衛(wèi)星A和衛(wèi)星B的靜電加速度計(jì)的仿真值如圖4和圖5所示，由于在飛行過(guò)程中為衛(wèi)星B領(lǐng)飛，衛(wèi)星A跟飛，因此衛(wèi)星B在y軸的非保守力加速度要略大于衛(wèi)星A，其余兩個(gè)方向則基本保持一致。

圖2 KBR系統(tǒng)測(cè)距觀測(cè)值Fig.2 KBR system range observation

圖3 衛(wèi)星A與衛(wèi)星B星間加速度觀測(cè)值Fig.3 Range acceleration observation between satellite A and B

圖4 衛(wèi)星A加速度計(jì)輸出值Fig.4 Accelerometer observation of satellite A

圖5 衛(wèi)星B加速度計(jì)輸出值Fig.5 Accelerometer observation of satellite B

靜電加速度計(jì)標(biāo)定算法仿真中，地球引力場(chǎng)模型采用120階TJGRACE02S模型［9］和120階EGM96模型［10］。衛(wèi)星A靜電加速度計(jì)的三軸標(biāo)度因子為1.0，零偏誤差為5×10－5m／s2；衛(wèi)星B靜電加速度計(jì)的三軸標(biāo)度因子為1.0，零偏誤差為5×10－5m／s2；星敏感器的姿態(tài)確定精度為0.05mrad；KBR系統(tǒng)的偏置誤差為1μm，系統(tǒng)噪聲誤差為0.5μm。GPS接收機(jī)的定軌精度為1cm，測(cè)速精度為0.1mm／s；衛(wèi)星質(zhì)量為500kg；轉(zhuǎn)動(dòng)慣量為diag（80，420，470）kg·m2；仿真中以1d的軌道弧長(zhǎng)為單位，對(duì)10d的軌道弧長(zhǎng)進(jìn)行10次標(biāo)定，仿真結(jié)果中給出了x軸方向的標(biāo)定結(jié)果（y軸方向和z軸方向與x軸方向具有相同的結(jié)論），如圖6～9所示。標(biāo)度因子和零偏誤差平均統(tǒng)計(jì)結(jié)果如表2所示。

由圖6～9可以看出：在不同時(shí)間，靜電加速度計(jì)的標(biāo)度因子和零偏誤差的標(biāo)定值是波動(dòng)的，這說(shuō)明標(biāo)度因子和零偏誤差并不是固定不變的。因此，對(duì)靜電加速度計(jì)數(shù)據(jù)精度要求較高的用戶，建議每天或半天甚至數(shù)小時(shí)估計(jì)一次標(biāo)度因子和零偏誤差，并根據(jù)其時(shí)間變化規(guī)律進(jìn)行曲線擬合。同時(shí)，可以看出：參考地球重力場(chǎng)模型的精度對(duì)校準(zhǔn)參數(shù)仍會(huì)有一定的影響，但總體來(lái)說(shuō)，計(jì)算結(jié)果相差不是很大，因而可認(rèn)為此算法校準(zhǔn)加速度數(shù)據(jù)對(duì)先驗(yàn)地球重力場(chǎng)模型的選擇并不是很敏感。

圖6 衛(wèi)星A靜電加速度計(jì)x軸方向標(biāo)度因子標(biāo)定結(jié)果Fig.6 Calibration results of scale factor along x axis of satellite A electrostatic accelerometer

圖7 衛(wèi)星A靜電加速度計(jì)x軸方向零偏誤差標(biāo)定結(jié)果Fig.7 Calibration results of bias along xaxis of satellite A electrostatic accelerometer

圖8 衛(wèi)星B靜電加速度計(jì)x軸方向標(biāo)度因子標(biāo)定結(jié)果Fig.8 Calibration results of scale factor along x axis of satellite B electrostatic accelerometer

圖9 衛(wèi)星B靜電加速度計(jì)x軸方向零偏誤差標(biāo)定結(jié)果Fig.9 Calibration results of bias along xaxis of satellite B electrostatic accelerometer

表2 雙星靜電加速度計(jì)標(biāo)度因子和零偏誤差均方差統(tǒng)計(jì)結(jié)果Table 2 Scale factor and bias calibration error of satellite A and B

由表2可以看出，利用EGM96和TJGRACE02S模型計(jì)算出的標(biāo)度因子和零偏誤差，估計(jì)誤差均在0.2%以內(nèi)，從統(tǒng)計(jì)角度而言，并無(wú)顯著差別，而從標(biāo)度因子和偏差參數(shù)的變化規(guī)律來(lái)看，也十分相似（個(gè)別地方有差別），這說(shuō)明本文給出的算法穩(wěn)定性好，也說(shuō)明EGM96模型與TJGRACE02S模型有較好的一致性。

4 結(jié)論

（1）本文提出的在軌標(biāo)定算法，充分考慮到靜電加速度計(jì)和KBR系統(tǒng)觀測(cè)數(shù)據(jù)之間的數(shù)據(jù)聯(lián)系，建立了靜電加速度計(jì)參數(shù)標(biāo)定的狀態(tài)方程和觀測(cè)方程，并采用擴(kuò)展卡爾曼濾波對(duì)雙星的加速度計(jì)進(jìn)行實(shí)時(shí)在軌補(bǔ)償。

（2）本文算法重點(diǎn)解決了動(dòng)力學(xué)方法中由于要進(jìn)行軌道積分而導(dǎo)致?tīng)顟B(tài)轉(zhuǎn)移矩陣十分復(fù)雜的問(wèn)題，有效簡(jiǎn)化了在軌標(biāo)定算法結(jié)構(gòu)，且計(jì)算量有所減少。

（3）利用不同的參考地球重力場(chǎng)模型分別進(jìn)行數(shù)學(xué)仿真驗(yàn)證，數(shù)值分析結(jié)果說(shuō)明該方法是有效的。此外，該方法對(duì)先驗(yàn)地球重力場(chǎng)模型的選擇不太敏感，但與所選參考地球重力場(chǎng)模型的精度仍有一定關(guān)系，因而在靜電加速度計(jì)標(biāo)定時(shí)應(yīng)盡量選擇精度較高的模型。

（4）本文算法也適用于單獨(dú)對(duì)尺度因數(shù)或零偏誤差進(jìn)行在軌標(biāo)定的情況。為了進(jìn)一步提高估計(jì)精度，可以充分利用該算法對(duì)尺度因數(shù)和零偏誤差的標(biāo)定結(jié)果，以及更多的可以估計(jì)到的誤差結(jié)果對(duì)系統(tǒng)進(jìn)行優(yōu)化，使系統(tǒng)達(dá)到最優(yōu)。

（References）

［1］周澤兵，白彥崢，祝竺，等.衛(wèi)星重力測(cè)量中加速度計(jì)在軌參數(shù)校準(zhǔn)方法研究［J］.中國(guó)空間科學(xué)技術(shù)，2009，29（6）：74－80 Zhou Zebing，Bai Yanzheng，Zhu Zhu，et al.In－orbit calibration methods of accelerometer parameters on satellite－borne gravimetry［J］.Chinese Space Science and Technology，2009，29（6）：74－80（in Chinese）

［2］祝竺，張曉敏，周澤兵.利用旋轉(zhuǎn)衛(wèi)星法開(kāi)展加速度計(jì)在軌檢驗(yàn)研究［J］.宇航學(xué)報(bào)，2010，31（5）：1362－1367 Zhu Zhu，Zhang Xiaomin，Zhou Zebing.In－orbit verification for accelerometers using rotating spacecraft method［J］.Journal of Astronautics，2010，31（5）：1362－1367（in Chinese）

［3］Visser Pnam，Ijssel J.Verification of CHAMP accelerometer observations［J］.Advance in Space Research，2003，31（6）：1905－1910

［4］Tom V，Eelco D，Pieter V.CHAMP and GRACE accelerometer calibration by GPS－based orbit determination［J］.Advance in Space Research，2009，24（8）：45－50

［5］徐新禹，李建成，王正濤，等.利用參考重力場(chǎng)模型基于能量法確定GRACE加速度計(jì)校準(zhǔn)參數(shù)［J］.武漢大學(xué)學(xué)報(bào)（信息科學(xué)版），2008，33（1）：72－75 Xu Xinyu，Li Jiancheng，Wang Zhengtao，et al.Calibration of GRACE accelerometer using reference gravity field model based on energy balance approach［J］.Geomatics and Information Science of Wuhan University，2008，33（1）：72－75（in Chinese）

［6］辛寧，邱樂(lè)德，張立華，等.USO－KBR測(cè)距系統(tǒng)建模與仿真［J］.航天器工程，2012，21（5）：91－95 Xin Ning，Qiu Lede，Zhang Lihua，et al.Modeling and simulation for USO－KBR ranging system［J］.Spacecraft Engineering，2012，21（5）：91－95（in Chinese）

［7］Kim J R.Simulation study of a low－low satellite－to－satellite tracking mission［D］.Austin：University of Texas at Austin，2000

［8］Furun Wang.Study on center of mass calibration and K－band ranging system calibration of the GRACE mission［D］.Austin：University of Texas at Austin，2003

［9］Flury J，Bettadpur S，Tapley B D.Precise accelerometry onboard the GRACE gravity field satellite mission［J］.Advances in Space Research，2008，42（8）：1414－1423

［10］Tapley B D，Ries J C，Bettadpur S，et al.Neutral density measurements from the gravity recovery and climate experiment accelerometers［J］.Journal of Spacecraft and Rockets，2007，44（6）：1220－1225

（編輯：夏光）

On－orbit Calibration Algorithm of Electrostatic Accelerometer for Gravity Measurement Satellite

XIN Ning1QIU Lede1ZHANG Lihua2LIU Naijin1
（1Institute of Telecommunication Satellite，China Academy of Space Technology，Beijing 100094，China）（2DFH Satellite Co.，Ltd.，Beijing 100094，China）

The electrostatic accelerometer is one of the most important instruments of the SST－LL（low－low satellite－satellite tracing）gravity measurement satellite which directly affects the measurement of the earth gravity field.In order to ensure long－term on－orbit working of electrostatic accelerometer and in combination with extended Kalman filter，a calibration algorithm based on dynamics is proposed for the accelerometer.Firstly，the observation model of the electrostatic accelerometer and KBR（K band ranging）system are obtained.Then，the accurate gravitational model and the measurements of electrostatic accelerometer are fitted into KBR－based extended Kalman algorithm.The associated simulation is performed to verify the availability and feasibility of the proposed calibration algorithm.The calibration accuracy of scale factors and biases can be better than 0.2%，which achieves the accurate calibration of satellite electrostatic accelerometer.

gravity measurement satellite；electrostatic accelerometer；K band ranging system；onorbit calibration

U666.1

：ADOI：10.3969／j.issn.1673－8748.2016.01.004

2015－03－26；

：2015－06－01

國(guó)家自然科學(xué)基金（91438205）資助項(xiàng)目

辛寧，男，博士，研究方向?yàn)樾l(wèi)星系統(tǒng)總體設(shè)計(jì)。Email：xinning7＠sina.com。