基于新息趨勢的納衛星地磁定軌*

陳貴芳 郁 豐 王 潤

1.“空間光電探測與感知”工信部重點實驗室，南京 210016 2.南京航空航天大學航天學院，南京 210016

0 引言

作為一種自主的、無長期積累誤差且具有較強抗干擾能力的導航技術[1]，地磁定軌技術適合對儀器復用率高的納衛星使用。

當前，地磁定軌多使用批處理算法和卡爾曼濾波算法進行數據融合。前者經多年發展，已有基于無跡變換的批處理算法(unscented batch filter, UBF)提出并應用于地磁定軌中[2]。但批處理算法的數據需提前采集，不適合衛星的實時定位。后者包含擴展卡爾曼濾波器(extended Kalman filter，EKF)和無跡卡爾曼濾波器(unscented Kalman filter，UKF)。EKF引入泰勒近似實現了對非線性系統的狀態估計。但地磁模型復雜，求導困難。文獻[3-4]建立了多項式擬合的地磁模型，文獻[5]則在計算觀測矩陣時采用低階地磁模型，均減小了系統的計算復雜度，但對定位精度提升不多。由于沒有忽略高階項，UKF的精度通常高于EKF。通過引入地平儀[6]，太陽敏感器[7-8]等傳感器，基于UKF的地磁定軌系統的定位精度可達到千米級。此外，上述傳感器也是常用的姿態敏感器，文獻[9-10]利用衛星的位置信息估計在慣性系下的地磁矢量信息和地心矢量(或太陽方向矢量)信息，再結合衛星的姿態信息將其轉換到衛星本體坐標系下，和傳感器的輸出值進行比較，實現了衛星的位置和姿態一體化估計。文獻[11-12]研究了地磁場幅值和衛星軌道根數在頻域上的關系，據此實現了衛星半長軸的高精度估計。文獻[13]在此基礎上實現了碟狀飛行器的姿態估計。

定義磁強計誤差和外源磁場擾動為地磁偏差，地磁定軌系統受其影響精度不佳。隨科技進步，后者的影響逐漸減小[14]。文獻[15]建立了磁強計誤差和地理緯度間的函數，在衛星工作時實時補償地磁信息，減小了定位誤差。但磁強計誤差還會隨著時間而變化，因此該方法不適合衛星長期使用。文獻[16]推導并建立了外源磁場的模型并將地磁場模型的系數作為狀態向量參與濾波，提升了地磁定軌的精度。但該方法計算量大，不適合納衛星使用。文獻[17]將地磁主磁場模型中構成磁耦極子模型的3個系數作為狀態向量參與濾波，仿真證明了該方法的可行性。外源磁場在磁暴時變化劇烈，文獻[18]對地磁偏差中的時變部分建模并改進了EKF算法，減小了磁暴下地磁定軌系統的位置誤差，但對無磁暴影響時的系統性能提升不大。

本文以地磁場幅值作為觀測量，推導了自適應卡爾曼濾波算法(adaptive Kalman filter，AKF)，并以此實現了納衛星的自主導航。該算法將地磁偏差近似建模為隨機游走模型并作為狀態變量參與濾波。然后，依據新息的變化趨勢，實時調整濾波器的參數。以Swarm-A的觀測數據的實驗表明，該算法有效提升了地磁偏差的估計精度，從而提升了系統的性能。和傳統濾波算法對比，AKF的收斂速度更快，其定位精度優于EKF，和UKF相當，最大定位誤差約6 km，但計算量卻遠小于UKF，適合納衛星等使用。

1 地磁定軌模型建立

1.1 狀態方程

以衛星的軌道動力學方程作為系統的狀態方程。為簡化地磁模型的計算，該方程描述于地心地固坐標系下。由于地磁定軌精度不高，本文只考慮了J2攝動項。

(1)

式中：μ為地球引力常數，Re為地球的參考半徑，ω為地球的自轉角速度，J2為攝動項系數，f=[fx,fy,fz]為攝動力。

對該方程離散化，得到離散型干擾式：

δXk,k-1=φk,k-1δXk-1,k-1+Γk,k-1Wk-1

(2)

1.2 觀測方程

選取國際地磁參考場(international geomagnetic reference field，IGRF)作為地球主磁場模型，其磁位勢函數如下：

(3)

將磁位勢對北向、東向和徑向求導，則地磁場在北東地坐標系下的分量可表示為：

(4)

選擇地磁場幅值作為觀測量，觀測方程為：

(5)

式中：Bb為地磁偏差，將其建模為隨機游走模型，并擴充到狀態變量中，η為觀測噪聲。

(6)

式中：ab為地磁偏差的驅動噪聲。

對觀測方程離散化，有：

δB=H·δX+η

(7)

(8)

式中：

地磁模型的高階導數對軌道修正能力弱[5]，為節約星上計算量，該偏微分陣僅考慮地磁模型的低階項：

(9)

式中：k為低階截斷階次。

2 自適應濾波算法推導

圖1是用Swarm-A衛星在2015年4月19日至20日的實測地磁幅值數據和解算IGRF模型獲得的地磁幅值數據做差獲得的地磁偏差曲線。由圖1可知，地磁偏差隨時間和衛星位置的變化而變化，因此常規卡爾曼濾波器的性能會受影響。本文提出了一種自適應濾波算法，依據新息序列的變化趨勢，實時調節濾波器的協方差矩陣和狀態噪聲矩陣，可提升地磁偏差的估計精度，從而提升地磁導航的精度。

定義系統的新息如下：

(10)

式中：Zk為系統觀測量。

新息vk的協方差矩陣可由式(11)獲得[20]：

(11)

式中：Pk為系統的協方差矩陣，R為系統的噪聲矩陣。

圖1 地磁偏差

在系統協方差矩陣時間更新方程中引入變量S(S≥1)，有：

(12)

則，

(13)

其中，tr(O)表示求矩陣O的跡。

因S≥1，有：

(14)

新息vk的協方差矩陣可用實際測量的最新N個新息的協方差均值近似表示：

(15)

式中：M為估計窗大小，由經驗給定，i0=k-M。

結合式(14)和式(15)，有：

(16)

本文認為，若式(16)的右式小于1，則S=1，否則，式(16)取等號。

綜上，AKF可表示為：

(17)

3 仿真校驗

首先分析各濾波器的定位性能。圖2給出了常規濾波器和自適應濾波器的位置誤差曲線。地磁場的可觀測性弱，濾波器收斂較慢。在0～20000 s間，各濾波器均有不同程度的超調，UKF最大，達到近44 km，AKF最小，僅17 km，并迅速減小到極小值。在系統運行20000 s后，其最大位置誤差小于5 km時濾波器收斂。從收斂速度上看，AKF最快，需要約40000 s，而EKF和UKF分別需要50000 s和60000 s才能收斂。從定位精度上看，AKF和UKF的精度最高，而使用EKF的地磁定軌系統的位置誤差最大，約10 km。因此，AKF收斂速度最快，精度和UKF相當，最大誤差約為6km。

接下來比較各濾波器的計算量。從內存占用上看。UKF通過sigma點集來逼近非線性函數，每次迭代相較于EKF和自適應算法需分別保存2n+1個狀態量和觀測量的sigma點集(n為狀態向量維數)，內存需求更大。從運行速度上看。使用一臺系統為Windows 10，安裝有matlab 2019a的電腦分別以這3種濾波器對同一組數據進行仿真，其仿真時間如表1。可見，UKF的運行速度遠低于EKF和AKF。因此，EKF和AKF對計算資源的需求遠低于UKF，更適合納衛星使用。

圖2 地磁導航位置誤差

表1 各濾波器仿真時間

常規濾波器的位置誤差曲線有類周期特征。圖3給出了地磁定軌系統(EKF)在僅做時間更新時和正常工作時的位置誤差曲線。當系統僅做時間更新時，曲線無類周期性，說明該波動與狀態方程無關。引入量測修正后，系統的位置誤差減小，同時引入了類周期波動。取120000～140000 s的位置誤差，在該區間內出現了約3.5個震蕩，每次用時約95 min，和Swarm-A的周期相同[22]。由于Swarm-A為極軌衛星，在不同周期，其在地磁場的位置不同。地磁場不同區域的可觀測性不同，因此在不同的衛星周期內，其定位性能不一，從而導致了該類周期性波動的出現。

圖3 常規濾波(EKF)和僅時間更新的位置誤差曲線

最后分析各濾波器對地磁偏差的估計性能。由圖4，EKF和UKF的地磁偏差估計曲線近乎重合，說明UKF并不能提升地磁偏差的估計性能。圖5是濾波器的地磁偏差的估計誤差曲線，這里只比較EKF和AKF。當系統開始工作時，AKF的地磁偏差的估計誤差迅速達到一個極小值，在圖2(a)表現為位置誤差的超調量遠小于常規算法。在20000～40000 s之間，AKF的地磁偏差的誤差和EKF相近，在位置誤差中表現為二者的最大誤差接近。40000s后2個濾波器逐步收斂，此后AKF對地磁偏差的估計性能始終優于EKF，其定位誤差也更小。因此，本文提出的自適應濾波算法有效地提升了地磁偏差的估計性能，進而提升了地磁定軌的性能。

圖4 常規濾波器的地磁偏差

圖5 地磁偏差的估計誤差

4 結論

對于地磁定軌易受地磁偏差影響的問題，針對性地將其建模為隨機游走模型并使之參與濾波。提出了自適應濾波算法，有效克服了空間磁場擾動和磁強計誤差的不確定性。實驗表明，該算法計算量低，定位精度高，優于傳統濾波算法，適合納衛星使用。