重力測量衛星KBR系統相位中心在軌標定算法

辛寧 邱樂德 張立華 丁延衛

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

1 引言

低低跟蹤（SST-LL）重力測量衛星的核心載荷為K 頻段測距（KBR）系統，該系統包括兩副喇叭天線（分別安裝在兩顆衛星上），用于測量KBR 系統相位中心之間的距離及距離變化率，測量誤差可達微米級［1］。由于重力場反演過程中需要的數據是兩顆衛星質心之間的距離及距離變化率，因此還要確定兩顆衛星在軌運行過程中衛星的質心與KBR 系統相位中心在視線方向上的距離。在地面研制階段，KBR 系統相位中心與衛星質心的位置都是標定好的，但衛星在發射階段中產生的振動、空間溫度變化及電磁噪聲等，都會使KBR 系統相位中心與質心間的相對位置發生改變，影響KBR 系統的在軌測量誤差。為了反演高精度的地球重力場，進行KBR 系統相位中心在軌標定，進而為KBR 系統幾何修正提供依據，具有重要的研究意義。

文獻［2］提出了KBR 系統相位中心在軌標定的基本原理，將KBR 系統測量的距離及距離變化率與衛星旋轉角的關系作為切入點，設計了衛星姿態機動規律，討論KBR系統相位中心在軌標定的可行性。文獻［3］從工程角度出發，提出將磁力矩器及姿態控制噴氣發動機為執行部件設計衛星的姿態機動，以KBR系統、GPS接收機、星敏感器和磁強計為測量敏感元件的估計算法。該算法將KBR系統相位中心在軌標定問題轉換為對衛星姿態四元數的求解，并采用Batch估計算法進行嘗試性的仿真研究。不過，這種算法是建立在衛星動力學特性精確已知的基礎上，而動力學模型在實際工程中不可能保證精確，存在著不確定性的誤差，在這種情況下采用Batch估計算法很難保證標定精度。文獻［4］提出估計狀態的預測濾波法。在預測濾波法中，模型誤差無需已知，而是作為解的一部分被估計，有效地解決了存在顯著模型誤差情況下的非線性問題；但是，該算法對初始狀態信息的依賴程度很高，當初始狀態信息存在較大誤差時，其估計精度明顯下降，甚至發散。

綜合上述算法的特點，本文提出了一種利用預測卡爾曼濾波算法估計衛星的姿態四元數，以擴展卡爾曼濾波算法實現KBR 系統相位中心在軌標定的算法。當存在較大的衛星姿態動力學模型誤差時，該算法仍能精確地實現KBR 系統相位中心的標定。同時，該算法采用三維偏差四元數代替四元數作為狀態變量，避免了估計過程中四元數單位約束性被破壞。

2 KBR 系統相位中心在軌標定算法

KBR 系統相位中心在軌標定算法的基本過程（見圖1）為：首先，利用姿態控制器發送指令給磁力矩器及姿態控制噴氣發動機，使衛星按照設計的姿態機動規律進行姿態機動；然后，將星敏感器數據代入預測卡爾曼濾波算法中，估計出衛星的姿態四元數；最后，利用擴展卡爾曼濾波算法實現KBR 系統天線相位中心的在軌標定。

圖1 數值仿真算法Fig.1 Numerical simulation algorithm

2.1 衛星姿態機動規律

設重力測量衛星為A 星和B 星，則A 星俯仰／偏航機動規律為［3］

B星俯仰／偏航機動規律為

式中：φA和φB分別為A 星和B 星的俯仰角；φA和φB分別為A 星和B星的偏航角；C和f0分別為機動的振幅和頻率；t為機動時間；θ0為初始偏置角。

衛星進行姿態機動，當θ0為正值時，定義為衛星正像機動；當θ0為負值時，定義為衛星鏡像機動。A 星的姿態機動過程為：首先進行衛星俯仰機動，衛星的其他方向變化都須保持原標稱姿態；俯仰機動結束后，衛星調整為三軸穩定狀態，再進行偏航機動；在A 星機動過程中，B 星的姿態始終保持原標稱姿態。B 星的姿態機動過程同A 星一致。在相位中心在軌標定算法中，只取機動過程中的數據進行標定。

2.2 KBR系統模型

KBR系統輸出的雙星距離Robs的表達式為［5］

式中：RA和RB分別為A 星和B 星的相位中心在J2000慣性坐標系下的位置矢量，見式（4）和式（5）；Rbr為KBR 系統 的零偏 誤差；Rnr為KBR 系統的測量噪聲誤差。

式中：rA和rB分別為A 星和B 星的質心在J2000慣性坐標系下的位置向量；dpcA和dpcB分別為A 星和B星的相位中心在衛星本體坐標系下的位置矢量；qA和qB分別為A 星和B 星的姿態四元數；M為衛星本體坐標系轉換到地心慣性坐標系的轉換矩陣，其表達式為

式中：q0，q1，q2，q3為衛星姿態四元數。

顯然，在構建上述系統模型時，要利用到衛星姿態四元數q的信息，其計算公式為

角動量H＝Jpω，Jp為衛星的轉動慣量矩陣；T為作用于衛星上的總力矩。

衛星運行期間，其轉動慣量矩陣Jp存在不確定性，而且，作用在衛星上的總外部力矩受外部干擾力矩和其他攝動因素的影響，往往無法精確已知［6-7］。如果直接利用衛星姿態動力學方程式（7）計算衛星的四元數信息，其估計結果未必可靠，這是KBR 系統在軌標定的一個難點所在。本文引入預測卡爾曼濾波算法來解決這一問題，首先利用預測濾波算法估計外部力矩誤差，然后利用卡爾曼濾波算法對衛星姿態進行估計，獲得高精度的坐標系轉換矩陣M，最后利用擴展卡爾曼濾波算法獲得KBR 系統相位中心的估計值。

2.3 應用預測卡爾曼濾波算法的姿態四元數估計

設狀態變量為X1＝［qH］T，在式（7）中增加力矩模型誤差d（tk），得到姿態估計的狀態方程為

式中：ωk，qk，Tk，Hk分別為衛星在tk時刻的角速度、姿態四元數、力矩和角動量，k為時間序號；I為單位陣；f（X1）為系統狀態模型函數；G為模型誤差分布矩陣。

姿態估計的觀測方程為

式中：和為tk＋1時刻的2個測量矢量；B1和B2為2個不平行參考矢量；h為系統觀測模型函數；vk＋1為觀測噪聲。

將觀測方程式（9）泰勒級數展開得［4］

式中：Δt為采樣時間間隔。

式中：為f（X1）的一階李導數；f（X1）的二階李導數

式中：W為3×3階加權矩陣；M為測量噪聲協方差矩陣。

將狀態方程式（8）離散化，保留二階項，可得到離散的狀態估計傳播方程如下。

式中：A［X1（tk）］＝

將tk＋1時刻星敏感器輸出測量值，tk時刻估計值代入式（13），求得模型誤差d（tk）；然后把d（tk）代入式（14）中，計算出tk＋1時刻的一步狀態估計值X1（tk＋1｜tk）。

為防止四元數的正交性奇異問題，對四元數誤差進行降階處理。設四元數誤差為矢部誤差，δq3為標部誤差。給定姿態初始狀態估計值及初始估計均方差誤差P1（t0），利用擴展卡爾曼濾波算法［8］就可以得到衛星姿態四元數的估計值。其中，擴展卡爾曼濾波算法中的姿態狀態轉移矩陣為

量測轉移矩陣為

則矢部誤差估計為

式中：Kq（tk＋1）為姿態反饋增益矩陣。

標部誤差估計為

姿態四元數的估計值為

2.4 應用擴展卡爾曼濾波算法的KBR 系統相位中心在軌標定

獲得衛星姿態四元數的估計值后，便可精確計算兩顆衛星的坐標轉換矩陣M，在此基礎上，設狀態變量X2＝［dpcAdpcB］T，建立式（20）中的狀態方程和觀測方程，利用擴展卡爾曼濾波算法遞推公式［8］即可估計出KBR 系統相位中心的偏差值。

式中：傳遞函數矩陣H2（t）＝［eTM1-eTM2］，e＝

3 數值仿真

設A 星和B星均運行在高度為450km、偏心率為0.001、傾角為89.9°的軌道上，兩星相距239km，它們的初始軌道和姿態參數如表1所示。

表1 雙星初始參數Table 1 Initial parameters of two satellites

仿真中采用的衛星攝動模型如表2所示，仿真中衛星的測量誤差參數如表3所示［9］。

表2 雙星動力學模型Table 2 Dynamic model of two satellites

表3 雙星測量誤差Table 3 Measurement errors of two satellites

衛星姿態機動過程中的機動規律設置為φ／φ＝±2＋sin（2πt／250）；為了減小多徑誤差對KBR 系統測距誤差的影響，衛星每種姿態機動時間取機動周期的整數倍，這里取1000s。為了分析得到數據采集的軌道相位與相位中心標定誤差的關系，分別選在赤道附近、北極附近進行衛星KBR 系統相位中心在軌標定。

定義重力測量衛星的本體坐標系為OXYZ，原點O位于衛星質心，X軸指向KBR 系統的相位中心，Z軸指向地心，Y軸與Z軸、X軸構成右手直角坐標系，則在衛星本體坐標系下赤道上空力矩模型估計誤差曲線如圖2所示。

圖2 力矩模型誤差估計曲線Fig.2 Torques model estimation errors

對圖2中的估值序列進行統計平均，得到力學模型誤差的穩態估值為［0.96 18.8 0.97］×10-6N·m。標定過程中，KBR系統輸出的距離觀測值見圖3，雙星距離大概變化560m。標定過程中，衛星質心與KBR系統相位中心之間的距離變化見圖4，其距離變化周期與姿態角機動周期相同，為250s，變化幅值為1.6mm。KBR系統的多徑傳播誤差見圖5，機動過程中多徑傳播誤差最大為1.45mm，變化周期為250s。北極附近的分析結果與赤道上空的一致。

在赤道上空，當初始狀態估計出現很大的誤差，X1（t0）＝［0.5 -0.5 0.5 0.5］T時，姿態估計結果見圖6和圖7。

圖3 KBR 系統距離觀測值Fig.3 KBR system range observation

圖4 衛星質心與KBR 系統相位中心之間的距離變化Fig.4 Distance changing between satellite centers of mass and KBR system phase centers

圖5 KBR 系統相位中心標定過程中多徑傳播誤差值Fig.5 Mutipath error during KBR system phase center calibration

圖6 俯仰機動姿態估計誤差Fig.6 Attitude estimation errors at pitch maneuver

圖7 偏航機動姿態估計誤差Fig.7 Attitude estimation errors at yaw maneuver

單獨使用卡爾曼濾波算法的估計值并不能收斂于真值，估計結果嚴重發散。預測卡爾曼濾波算法可以實現姿態四元數的準確估計，三軸估計誤差均在20μrad以內。這是由于預測卡爾曼濾波算法采用更新量測矩陣信息，保證了模型誤差估計值的準確性，因而其魯棒性更強。需要指出的是，預測卡爾曼濾波算法的上述性能是以犧牲運算量為代價的，其運算量為4430次乘法和2986次加法，約是預測濾波算法的3倍，但是要小于通常采用的卡爾曼濾波算法。

設相位中心估計值為dePCi，它與理論值dPCi之間的夾角θier作為標定誤差。根據文獻［10］，如果要將由KBR 系統相位中心偏移引起的測距誤差控制到1μm，則標定誤差θier必須在0.3mrad以內。分別采用文獻［3］中的Batch估計算法及本文提出的算法，在相同的輸入條件下得到相位中心的標定結果，如表4所示。當初始狀態估計出現很大的誤差時，使用Batch估計算法結果無法收斂，而采用本文提出的算法，可以達到很好的估計結果。從標定結果可以看出，無論在赤道地區還是在北極地區，單獨進行正像俯仰／偏航機動或者鏡像俯仰／偏航機動，相位中心的標定結果均不能滿足0.3mrad的精度要求，必須分別進行正像／鏡像機動，然后將數據進行聯合處理，才可達到精度要求。其主要原因是：單獨進行正像機動或鏡像機動，受多徑傳播誤差影響較大；而將正像機動與鏡像機動數據相結合，可有效抵消多徑傳播誤差影響，使相位中心估計精度大幅度提高。總體而言，在該算法中，標定軌道相位的選取，對最終標定誤差的影響并不是特別顯著，赤道地區略優。從激勵時長看，單次機動1000s的總機動時間是可以接受的，激勵時間再長，噴氣發動機燃料消耗較大，不利于衛星的長壽命運行。

表4 KBR 系統相位中心在軌標定誤差Table 4 KBR system phase center in-orbit calibration errors mrad

4 結束語

KBR 系統相位中心在軌標定的一個關鍵問題，是如何高精度地獲取衛星在軌機動時的姿態轉換矩陣。針對這一問題，本文提出了一種應用預測卡爾曼濾波算法和擴展卡爾曼濾波算法的KBR 系統相位中心在軌標定算法，使KBR 系統相位中心標定誤差優于0.3mrad。仿真結果表明：該算法魯棒性強，對初始姿態誤差不敏感，能有效利用預測卡爾曼濾波算法的優點。目前，濾波器的設計還有些復雜，在實際工程應用還要結合KBR 系統的具體工作特性進行深入研究，進一步優化濾波器結構，為重力測量衛星的工程化設計服務。

（References）

［1］佘世剛，王鍇，周毅，等.高精度星間微波測距技術［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（in Chinese）

［2］Romans L.Optimal determination quaternion from multiple star camera sensors［R］.Pasadena，California：NASA JPL，2003

［3］Wang Furun.Antenna phase center determination of inter-communicating satellites［C］／／Proceedings of the AAS／AIAA Astrodynamics Conference. Washington D.C.：AIAA，2002：212-218

［4］John L.Predictive filtering for attitude estimation without rate sensors［J］.Journal of Guidance，Control and Dynamic，1997，15（3）：522-527

［5］張紅軍，趙艷彬，孫克新.星間微波測距系統相位中心在軌標定研究［J］.上海航天，2010，27（4）：1-5 Zhang Hongjun，Zhao Yanbin，Sun Kexin.Study on phase centers calibration of K-band ranging system during in-flight phase［J］.Aerospace Shanghai，2010，27（4）：1-5（in Chinese）

［6］Hall C D，Tsiotras P，Shen H.Tracking rigid body motion using thrusters and momentum wheels［J］.The Journal of the Astronautical Sciences，2002，50（3）：311-323

［7］Bang H，Choi H D.Attitude control of a bias momentum satellite using momentum of inertia ［J］.IEEE Transactions on Aerospace and Electronic System，2002，38（1）：243-250

［8］王本利，廖鶴，韓毅.基于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（in Chinese）

［9］Kim J R.Simulaton study of a low-low satellite-to-satellite tracking mission［D］.Austin：University of Texas at Austin，2000

［10］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