一種重力衛星質心在軌標定算法

辛寧 邱樂德 張立華 丁延衛 王大雷

（1航天東方紅衛星有限公司，北京100094）（2中國空間技術研究院，北京100094）

1 引言

重力場探測衛星的核心有效載荷為測量非保守力的靜電加速度計，用于確定星間距離的K波段微波星間雙向測距儀（K Band Ranging System，KBR）。當衛星質心發生漂移時，靜電加速度計內部懸浮質量塊質心與衛星質心無法重合，嚴重影響了加速度計的測量精度［1］，同時衛星質心的偏差對KBR相位中心的幾何標定也有一定的影響［2］。然而，由于衛星在發射前質心校準存在誤差，發射中經受振動環境，衛星入軌后姿態調整飛行中燃料損失等原因均會引起衛星質心發生漂移。為了使衛星在軌飛行測量中，達到懸浮質量塊質心和衛星質心的最大重合，提高靜電加速度計和KBR的測量精度，進而反演出高精度的地球重力場信息，必須進行重力衛星質心的在軌標定研究［3］。

文獻［4］提出了一種利用衛星自旋晃動，并通過靜電加速度計實現衛星質心在軌標定的方案。但是該在軌自旋方案控制難度較大，工程上難以實現。文獻［5］從工程實現的角度出發，提出了一種利用周期性磁力矩使衛星產生姿態機動，利用經典的Kalman濾波算法對衛星質心進行檢測的方案。此方案將衛星質心在軌標定問題轉換為對衛星角加速度的求解，并進行了嘗試性的仿真研究。但該算法需假設衛星加速度模型近似線性，且過程噪聲為高斯白噪聲，如果實際的角速度響應中存在著非線性、非平穩過程，則會降低質心估計精度。文獻［6］提出了一種利用最小模型誤差準則（MME）估計衛星角加速度，利用EKF估計衛星質心的算法。由于該算法的測量元件中增加了陀螺儀，因此不需要再進行衛星角速度的估計，節省了一定的計算量。MME算法雖然無需利用假設的系統模型即可獲得角速度估計，但是MME估計器是一種離線估計算法，無法對質心進行實時估計。為了解決實時性問題，在經典的卡爾曼濾波和MME估計器的基礎上，文獻［7］提出了一種應用于無陀螺角速度確定預測濾波算法，這種算法不再局限于系統模型誤差只是高斯白噪聲的情況，能夠對空間飛行姿態進行實時估計，有效解決了存在顯著動態模型誤差情況下的非線性濾波問題。

在上述文獻的基礎上，提出了一種使用靜電加速度計和陀螺數據估計衛星質心的算法，其將預測濾波與擴展卡爾曼濾波相結合，以衛星角加速度和質心偏移量為狀態估計量，將衛星角加速度作為系統模型誤差可以實時進行估計。該算法中的濾波器不需要衛星的控制輸入以及動力學特性信息，針對衛星質心及角加速度能夠實時估計，具有快速的收斂性能和良好的穩定精度。

2 衛星質心在軌標定算法

衛星質心的在軌標定基本原理如圖1所示，根據地球磁場模型及衛星的軌道位置，利用磁力矩器給衛星施加一定頻率的力矩，使衛星產生相應的轉動運動。采用加速度計、陀螺儀作為檢測傳感器，對陀螺儀觀測矢量進行預測濾波處理，估計出角加速度；對靜電加速度觀測矢量進行擴展卡爾曼濾波處理，分離出具有線性特征的非重力加速度分量，由此來確定衛星質心的位置。

圖1 重力衛星質心在軌標定原理Fig.1 Principle of on-orbit calibration of center of mass for gravity satellite

2.1 系統模型

靜電加速度計在衛星本體坐標系下輸出的觀測矢量表達式可表示為［5］

式中Aobs為靜電加速度計在衛星本體坐標系下的輸出；RAB為安裝誤差矩陣；s1為加速度計的標度因數；Ab為零偏誤差；An為加速度計的噪聲誤差；aout為靜電加速度計的實際輸出：

式中ω為衛星角速度觀測數據；為衛星角加速度；l為衛星質心偏差；gg為由于質心偏差引起的重力梯度加速度；ang為非保守力加速度。

陀螺儀的觀測數據可表示為：

式中k1為標度因數；ωs為加速度計輸出；ωb為零偏誤差；ωn為陀螺的噪聲誤差。

由文獻［5］可知，由重力梯度和非保守力因素引起的加速度分量在較短的標定時間t內基本保持線性關系，因此質心在軌標定問題的測量方程可表示為

式中α和β為aout、gg以及Ab線性化近似后的斜率估計量和截距估計量；為等效矩陣：

顯然，在構建上述測量方程時，需要利用衛星的角加速度信息。剛體衛星角加速度為

式中Jp為轉動慣量；M為作用于衛星上的合力矩。由于Jp和M在質心標定過程中存在不確定性，因此利用公式（4）無法實時準確的得到衛星角加速度。為了彌補上述公式的缺點，本文首先將陀螺數據代入預測濾波器中，實時獲取到角加速度信息，然后利用擴展卡爾曼濾波對質心位置進行標定。

將式（4）右邊的狀態信息全部定義為模型誤差d（t），由此獲得的系統狀態方程為

式中G為d（t）的系數矩陣。

觀測方程為

式中為陀螺儀輸出；為靜電加速度計輸出。觀測值

2.2 角加速度預測濾波估計

預測濾波器的物理概念：設已經獲得了tk時刻的狀態估值，采樣間隔時間為Δt，在接到tk＋Δt時刻的測量信息后，利用模型最小誤差逼近可預測 ［tk，tk＋1］間隔內的模型誤差d（t），然后將d（t）的預測值代入系統狀態方程和觀測輸出方程中，將狀態估值傳播到tk＋1時刻，從而得到較高的估計精度。

將觀測方程（7）泰勒級數展開得到［7］

將式（9）代入式（10）后求導，得到最小模型誤差：

其中

式中ci為系統觀測模型函數c的第i行；Λ（Δt）為對角矩陣，其對角元素為λii＝Δtpi／pi！，pi是對公式（7）微分時中首次出現d（t）中元素的最低階數；是ci關于f的k階李導數；si為S的第i行元素；j為d（t）的行索引；gj是矩陣G的第j列元素；R為量測噪聲；W為半正定模型誤差加權矩陣。

利用上述迭代公式便可實時精確地獲取衛星角加速度˙ω（t）的估計值。

2.3 質心的擴展卡爾曼濾波估計

在衛星角速度及角加速度估計值已知的條件下，根據系統狀態方程式（6）及觀測方程式（8）對X1＝［lαβ］進行估計。

待估狀態變量一步預測均值估計為

式中Φ（tk）＝Ⅰ9×9為狀態轉移矩陣。

一步估計均方誤差陣為P（tk＋1／k）＝Φ（tk）P（tk）ΦT（tk）

反饋增益矩陣為K（tk＋1）＝P（tk＋1／k）H（tk＋1）T（H（tk＋1）P（tk＋1／k）H（tk＋1）＋Q（tk＋1））-1，其 中，為量測轉移矩陣；Q（tk＋1）為量測噪聲矩陣。

估計均方差陣為P（tk＋1）＝（Ⅰ-K（tk＋1）H（tk＋1））P（tk＋1／k）

狀態估計值：

上述質心標定過程的實質是先利用預測濾波預測下一時刻的模型誤差校正量d（t），然后再對狀態估值進行傳播估計，從而有效地提高了角加速度的估計精度。然后將角速度及高精度的角加速度估計值代入中更新量測轉移矩陣H，根據公式（12）實現了衛星質心的高精度估計。

3 數值仿真

衛星初始時刻：2015年1月1日正午（UTC）。

衛星動力學參數［8］：10階EMG96重力場模型；DTM大氣模型，阻力系數Cd＝2.2；太陽光壓反照系數η＝0.5；地球磁場模型，13階IGRF2005模型。

衛星仿真參數［9］：質量m＝500kg；初始質心位置三軸偏差磁力矩器力偶mm＝ ［0 1 -sgn（ByBz）］T·30sin（2π×0.1t）（A·m2），By、Bz為地球磁場在衛星本體坐標系下的分量；靜電加速度計標度因數s1＝ ［1＋0.01% 1＋0.01% 1＋0.01%］T；零偏誤差Ab＝ ［2×10-62×10-6-2×10-6］Tm／s2，測量誤差An＝ ［1×10-101×10-101×10-10］Tm／s2；陀螺儀 標 度 因 數k1＝ 1＋0.1% ，三軸零偏誤差ωb＝［0.3 0.3 0.3］／［（°）／h］，測量誤差ωn＝［0.1 0.1 0.1］／［（°）／h］。

加速度計與陀螺儀的輸出周期為0.1s，預測估計濾波算法的周期為0.1s，為了檢驗所設計的磁力矩器力偶信號的持續時間與質心標定精度的關系，預測估計濾波算法的運行時間分別設定為1min、3min、5min、10min。為了分析得到數據采集的軌道相位與質心標定精度的關系，分別選在赤道附近及南極附近進行衛星質心標定。

在赤道附近進行姿態機動，圖2為地球磁場強度仿真結果，可以看出三軸磁場強度在短時間內可認為是近似線性變化。圖3為施加在衛星上的周期性磁力矩仿真結果，三軸磁力矩的最大幅值分別為0.5mNm、0.75mNm和0.75mNm，形式為正弦，頻率為0.1Hz。線性加速度的仿真時序如圖4所示，其中X方向的線性加速度近似效果相對較差，主要原因為X方向上大氣阻力的線性度比Y方向和Z方向要差很多，因此線性度較差。圖5為陀螺儀輸出的角速度數據，由圖5中的角速度數據計算的角加速度估計值如圖6所示，預測濾波算法具有很好的估計精度和可靠性，三軸估計精度分別為2.939×10-6rad／s2，0.975×10-6rad／s2，0.858×10-6rad／s2；該算法僅利用了陀螺提供的角速度便能實時地估計出高精度的衛星角加速度，南極上空的現象與赤道上空一致，這里不再贅述。通過圖7中的靜電加速度計的輸出仿真值，利用擴展卡爾曼濾波得到的質心標定結果如圖8所示，由于X方向的線性加速度線性度較其他兩個方向相對較差，因此質心X方向的估計時間較長，Y方向和Z方向的估計時間較短。衛星質心在軌標定的軌道相位、精度和標定時間的關系如表1所示。

圖2 赤道上空三軸磁場強度分布Fig.2 Magnetic flux density over equator

圖3 赤道上空三軸磁力矩Fig.3 Magnetic torque simulation over equator

圖4 赤道上空線性加速度Fig.4 Linear acceleration over equator

圖5 赤道上空陀螺儀觀測值Fig.5 Gyro observation over equator

圖6 赤道上空角加速度估計值Fig.6 Estimation of angle acceleration over equator

圖7 赤道上空靜電加速度計觀測值Fig.7 Accelerometer observation over equator

圖8 赤道上空質心卡爾曼濾波估計值Fig.8 Kalman estimation of center of mass over equator

表1 衛星質心在軌標定誤差精度Tab.1 Center of mass calibration error of satellite mm

由表1可以看出，X軸的最佳標定精度為0.062 8mm，Y軸的最佳標定精度為0.032 4 mm，Z軸的最佳精度為0.041 4mm。本文提出的質心標定算法滿足標定精度需求。在相同條件下，在赤道附近進行姿態機動，質心標定結果要略優于在南極附近機動的結果。其主要原因是赤道附近激勵起的三軸角加速度較為均勻，這樣利于辨識精度的提高；而南極激勵的三軸角加速度幅值差異較大。通過不同的激勵時長看，3min或5min的總機動時間是可以接受的，再長的激勵時間，則質心測量精度反而降低。其主要原因是，辨識中假設非重力加速度項在辨識期間保持線性關系，而隨著辨識時間的增加，這種假設的誤差就越大，反而不利于質心辨識精度的提高。

4 結束語

綜上所述，目前關于重力衛星質心的標定方法均為離線估計算法，無法實現質心的高精度實時估計，針對這一問題，本文通過設計估計角加速度的預測濾波器與質心擴展卡爾曼濾波器實現了衛星質心的實時估計，仿真結果表明該算法具有較高的精度，魯棒性能強，可為質心在軌修正及KBR天線相位中心在軌標定提供設計依據。同時考慮該算法對陀螺儀的精度要求較高，下一步將研究在陀螺儀精度不高的情況下，如何確保算法的精度。

［1］周澤兵，白彥崢，祝竺，等.衛星重力測量中加速度計在軌參數校準方法研究 ［J］.中國空間科學技術，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.

［2］佘世剛，王鍇，周毅，等.高精度星間微波測距技術 ［J］.宇航學報，2006，27（3）：403-406.SHE SHIGANG，WANG KAI，ZHOU YI，et al.The technology of high accuracy inter-satellite microwave ranging ［J］.Journal of Astronautics，2006，27（3）：403-406.

［3］陳光鋒，唐富榮，薛大同.重力衛星在軌質心修正原理 ［J］.宇航學報，2005，26（5）：567-570.CHEN GUANGFENG，TANG FURONG，XUE DATONG.The trim principle of center of mass of gravity satellite during orbit flight［J］.Journal of Astronautics，2005，26（5）：567-570.

［4］TANGYIN S，WILLIAMS T.Mass property estimation using coasting maneuvers［J］.Journal of Guidance，Control and Dynamics，1997，20（4）：625-632.

［5］WANG F，TAPLEY B D.Precise determination of accelerometer proof mass ［C］.AAS／AIAA Space Flight Mechanics Meeting，Texas，2002：739-752.

［6］王本利，廖鶴，韓毅.基于MME／EKF算法的衛星質心在軌標定 ［J］.宇航學報，2010，31（9）：2150-2156.WANG BENLI，LIAO HE，HAN YI.On-orbit calibration of satellite center of mass based on MME／EKF algorithm ［J］.Journal of Astronautics，2010，31（9）：2150-2156.

［7］CRASSIDIS J L.Predictive filtering for attitude estimation without rate sensors ［J］.Journal of Guidance，Control and Dynamics，1997，20（3）：522-527.

［8］KIM J R.Simulation study of a low-low satellite-to-satellite tracking mission ［D］.Austin：University of Texas at Austin，2000.

［9］WANG FURUN.Study on center of mass calibration and K-band ranging system calibration of the GRACE mission［D］.Austin：University of Texas at Austin，2003.