深空探測航天器巡航段自主導航方法研究

葉子鵬，周慶瑞，王輝

中國空間技術研究院 錢學森空間技術實驗室，北京 100094，中國

21世紀以來，全球各國開始制定火星探測計劃，中國擬在2020年實現火星探測。航天器的精確導航將是該計劃最為關鍵的步驟之一。通常情況下，航天器的軌道參數可以通過地面測控站確定或者通過接收GNSS信號自行確定。然而對于深空探測任務，地面測控的方式不僅存在時延問題，對于一些無法在全球布置測控站的國家而言，還存在信號遮擋問題。GNSS則存在覆蓋問題、信號強度問題，以及對深空探測航天器的量測幾何差。因此，自主導航對航天器提高生存能力、有效完成多項任務具有重大意義。

由于圖像處理、敏感器等技術原因，直到1999年，NASA在“新千年計劃”中的第一顆探測器“深空一號”才完成了自主導航技術的在軌驗證。其方法是利用背景恒星和小行星進行自主導航，導航的位置精度為250 km(3σ)，速度精度0.5 m/s(3σ)。隨著技術發展，自主導航技術在深空探測中已具有一定的工程可行性，得到了越來越多研究者的關注。由于目前工程上自主導航的精度低于地面測控方式，美國的火星探測任務在巡航段主要采用了地面站測控的方法[1-2]，但同時也在大力發展自主導航技術，在多項火星探測任務中開展了天文光學導航的技術驗證。隨著自主導航技術逐漸成熟，未來的深空探測中自主導航必然會成為主要的導航方法。JPL已經開發了能夠實現行星際自主導航(圖像處理、軌道確定、機動計算)的軟件“AutoNav”[3]。

對于自主導航方法，雖然有學者提出了X射線脈沖星導航[4-6]的方法，但是由于可用脈沖星數量少，且脈沖信號抵達太陽系質心的時間數據庫[7]尚未完全建立，脈沖星導航精度無法保證，依然還處于研究階段。因此深空自主導航目前主要依賴于天文光學測量。Yim等[8]提出了利用航天器相對于太陽的徑向速度與航天器同太陽、地球視線矢量相結合的自主導航方法。該方法的仿真精度優于100 km，但是其計算過程中并未考慮動力學模型的噪聲，另外對測量精度要求非常高，速度量測精度為cm/s 量級，角度測量精度為亞角秒級，而且航天器與地球的位置有較強的約束。因此該仿真結果不適用于某些深空探測航天器的自主導航。Karimi等[9-10]提出了對量測量的修正方法，對行星視線矢量進行了光時延校正和光像差校正，但是未對濾波算法進行深入研究。

本文提出了一種用于火星巡航段軌道的自主導航方法。該方法首先測量航天器相對于太陽的徑向速度，利用擴展卡爾曼濾波估計出航天器狀態信息；同時測量多顆小行星的視線矢量，結合最小二乘法估計出航天器位置信息。最后利用一種改進的信息融合算法，對兩者的估計值進行有效融合。相比于其他深空探測自主導航方法，該方法對動力學模型精度的依賴度低，位置確定精度高，對于未來的深空探測具有一定借鑒意義。

1 系統狀態模型

本文采用的坐標系是J2000日心黃道坐標系，對于導航方案的總體設想如圖1所示。

圖1 導航方案流程Fig.1 Flow diagram of navigation

令航天器的狀態矢量為：

X(k)=[x,y,z,vx,vy,vz]T

(1)

式中：x,y,z分別為航天器三軸位置；vx,vy,vz分別為航天器三軸速度。以太陽為中心引力體的航天器動力學方程表示如下：

式中：r為航天器相對于太陽的矢徑，r為向量r的模；v為航天器相對于太陽的速度矢量；μs為太陽引力常數；μi為第i顆行星的引力常數；rip是航天器相對于第i顆行星的矢徑；rsi是第i顆行星相對于太陽的矢徑；a中包含太陽光壓等未建模項。

根據動力學方程，可以將狀態模型用下式表示：

式中：f[X(t)]為狀態方程；W(t)為過程噪聲模型。

2 系統量測模型

2.1 徑向速度量測

當航天器相對于某恒星在徑向上靜止時，接收到的恒星光譜恒定不變。而當航天器存在徑向速度時，接收到的恒星光譜將會產生漂移，其漂移量與徑向速度成正比。因此可通過測量光譜漂移量來獲得航天器相對于某一恒星的徑向速度。

分光計基于這一原理，能夠達到cm/s量級[8]的徑向速度測量精度。因此在進行自主導航研究時，多數學者都選擇將徑向速度作為確定航天器軌道的量測量。但是只有徑向速度作為量測量的模型是不可觀的，為了提高自主導航精度，通常還需要加入其他信息量。軌道模型示意如圖2所示。

圖2 航天器軌道示意Fig.2 Diagram of spacecraft orbit

徑向速度的大小vr可以表示為：

量測模型可表示為：

Z(t)=h[X(t)]+V(t)

式中：h[X(t)]為模型量測方程；V(t)為量測噪聲。

2.2 小行星視線矢量量測

在巡航段時，小行星帶的小行星數量多，距離較地球等大行星離航天器更近，對航天器位置確定的能觀度較好。因此選擇量測小行星帶的小行星視線矢量來確定航天器位置。

星光方向測量多用由鏡頭光學系統和CCD光敏元件組成的星敏感器來完成，它可以建立如圖3所示的測量坐標系。

圖3 視線矢量測量模型Fig.3 Diagram of sight vector measurement model

由CCD陣列的測量數據經數據處理可以得到圖像中心位置(Px,Py)，從而可以得到視線矢量在測量坐標系上的坐標Vp為：

(2)

式中：f為光學系統的焦距。

由已知星敏感器在衛星本體上的安裝矩陣M，結合式(2)可得到視線矢量在衛星本體坐標系中的坐標為：

VB=MTVp

式中：M為星敏感器安裝矩陣。

通過獲取視線矢量在衛星本體坐標系中的方向以及衛星在J2000日心黃道坐標系下的姿態角，就可以求出視線矢量在日心黃道坐標系中的方向。

2.3 系統的可觀測度分析

當僅使用徑向速度作為量測量時，系統是不可觀的。但是結合動力學模型進行卡爾曼濾波時，徑向速度的量測可以提高狀態估計的精度。

通過測量多顆小行星視線矢量，航天器的位置信息是可觀的，利用多個小行星視線矢量可確定航天器的位置。

本文在仿真中對系統的可觀測度進行了計算分析。對于線性模型，系統的可觀測度計算過程如下：

式中：Φi,k為第k步到第i步的狀態轉移矩陣；Hi為第i個量測值的量測矩陣；Ri為該量測值的噪聲方差陣；n為量測值的維度；W為系統的能觀度矩陣，當其非奇異時，系統可觀；η為系統的可觀測度，系統不可觀時為0。對于系統可觀測度的分析有助于在小行星星歷表中選擇合適的小行星進行導航。

由于本文中的量測模型為非線性模型，因此采用的是基于Fisher信息陣[11-13]的系統可觀測度分析，令觀測向量的概率密度為Pz，同時令量測噪聲服從高斯分布，對于一組觀測序列：

Z=[Z1,Z2,…,Zk]T

其中Z1,Z2,…,Zk分別為各量測量。其概率密度Pz(Z|x)的對數表示形式為：

lnLx(Z)=ln[Pz(Z|x)]=

根據Fisher信息陣計算方法得信息陣：

式中：X為對象狀態信息。根據Cramer-Rao不等式，估計均方誤差陣P有：

P≥J-1

因此，可以用信息陣的逆來估計可觀測度，其表達形式為：

式中：J為信息陣。

對選擇不同數量導航小行星的模型可觀度進行計算，得到結果如圖4所示。

圖4中系統可觀測度大于0，說明系統的狀態是可觀測的。同時通過觀察可觀測度變化的規律可以得到，選擇的導航星數越多，系統可觀性越好。

在選擇相同數量的導航小行星時，可觀測度的上下波動可以反映航天器與導航小行星之間的距離變化：可觀測度越高，航天器與導航小行星之間的距離越近；可觀測度越低，航天器與導航小行星之間的距離越遠。

選擇不同星數的最小二乘濾波結果的誤差平滑后如圖5所示。

根據圖4與圖5，可以得到在利用較少數量導航小行星的情形下，當可觀度下降時，濾波誤差呈現上升趨勢；當可觀度上升時，濾波誤差呈現下降趨勢。

同時可以得到，隨著導航小行星數量增加，系統可觀測度上升，濾波誤差下降。但是當導航小行星數提升至12顆以后，濾波精度提升不大。綜合考慮航天器的量測代價，故選擇導航小行星數量為12顆。

圖4 不同數量導航小行星的系統可觀度Fig.4 Observability with different number of navigation asteroids

3 濾波算法

3.1 擴展卡爾曼濾波

擴展卡爾曼濾波方法計算簡便，而且可以獲得較好的濾波結果，廣泛應用于航天器自主導航中。因此本文采用擴展卡爾曼濾波對航天器狀態進行估計，離散的狀態模型可以簡單化地表示如下：

X(k)=f[X(k-1)]+W(k-1)

式中：X(k)如式(1)所示；W(k)為過程噪聲；V(k)為量測噪聲。令：

E[W(k)W(k+i)T]=Q(k)δ(i)

E[V(k)V(k+i)T]=R(k)δ(i)

式中：Q，R分別為過程噪聲與量測噪聲的均方差陣。當且僅當i=0時，δ(i)為單位陣，否則為0矩陣。

則擴展卡爾曼濾波過程如下所示：

3.2 最小二乘法

最小二乘法的優點是方法簡單，使用過程中不需要任何先驗的狀態信息，也不需要動力學模型的參與，因此可以避免因動力學模型不準確而造成濾波誤差的問題。同時，最小二乘法得到的位置精度會隨導航小行星數量的增加而提高，而小行星帶恰好能夠提供足夠數量的導航小行星,并且當航天器越靠近小行星，自主導航的精度會越高。最小二乘法小行星視線矢量的分布如圖6所示。

對于小行星的矢量距離Sn有：

Sn=lnρn

式中：ln為第n顆小行星的視線矢量；ρn為航天器同第n顆小行星的標量距離。可以得到如下等式：

圖6 小行星量測示意Fig.6 Diagram of asteroid measurement

r=r1+l1ρ1=r2+l2ρ2=…=rn+lnρn

(3)

式中：l1到ln為量測的視線矢量，其形式為3維的列向量；r1到rn可通過查詢小行星的星歷表獲取，均為3維的列向量。小行星的星歷誤差通常為100 km(3σ)，在仿真中設置小行星視線矢量量測誤差為角秒級時，可以滿足仿真過程中所需的星歷誤差條件。根據式(3)可以建立如下方程：

一是現行法律規定不明導致訴訟主體不清。民事訴訟法將公益訴訟主體限定在“法律規定的機關和有關組織”，但未對其具體范圍作出明確規定。在水資源保護公益訴訟領域，“法律規定的機關”是否包括水行政主管部門，抑或還包括其他機關，尚須進一步明確；至于“有關組織”，其范圍更加廣泛，是僅指所有在民政部門登記的組織（合法社會組織），還是指須經相關法律特別授權的社會組織，目前也不明確。此外，“法律規定的機關”和“其他組織”之間的訴權順序，目前也缺乏規定。

(4)

ri可視作已知量，故式(4)可表示為：

li=hi(r)

基于遞推最小二乘算法，狀態r的迭代可用算法中的X等價表示，計算過程如下：

式中：Wk+1為量測量的權重系數；Hk+1為

由于量測方程存在線性化誤差，仿真試驗發現，通過遞推最小二乘法所得結果的精度略差于通過下式批處理法所求結果的精度：

(5)

式中：I3×3為3階單位陣。對于不相容的方程：

Ax=b

式(5)中有3n個方程，需要解得n+3個解，當n≥2時，存在矩陣A列滿秩，則有最小二乘解：

x=(ATA)-1ATb

(6)

利用式(6)即可得到航天器位置信息的最小二乘解。

3.3 信息融合

(7)

式中：

βi(k)為兩種估計值的權重；P1(k)為擴展卡爾曼濾波在第k步估計后的均方誤差陣；P2(k)為最小二乘估計均方誤差陣，由于在本文所仿真的場景中最小二乘估計所得誤差變化不大，因此可以將其設置為一個不變的矩陣。本文仿真中將P2(k)的2-范數值設置為 106～107之間的固定值，該值大約為最小二乘估計所得位置誤差的平方值。若是P2(k)變化較大，則可以用最小二乘估計均方差陣的迭代算法進行計算。

4 仿真試驗

4.1 仿真條件

導航小行星的選擇方法來源于文獻[14]，小行星的星歷數據來源于JPL的HORIZONS System。

相關學者對深空探測轉移軌道[15-16]進行了深入研究，本文采用的火星巡航段軌道是雙切軌道，考慮在J2000日心黃道坐標系下航天器Z軸狀態分量較小，為簡化計算，只考慮X、Y軸的位置和速度狀態信息。由于初始的位置誤差不影響最終的濾波結果，多數相關文獻的仿真試驗中將其設置為幾十到幾百千米，因此本文將其設置為500 km，并且將速度誤差設置為0.01 km/s，徑向速度量測誤差設為10-5km /s，小行星的視線矢量量測誤差設為4″(3σ) 。當選擇濾波周期過短，不僅會占用過多星上資源，而且相鄰兩步中量測量變化過小，甚至被量測噪聲覆蓋；當濾波周期過長，仿真精度下降。文獻[9]中將濾波周期設置為50 min左右，為了便于比較，本文設置仿真中的濾波周期為3 600 s，仿真總時長大約為140 d。

4.2 仿真結果

計算可以得到太陽對火星巡航段航天器的加速度為10-6km/s2量級，航天器在火星巡航段受到擾動總和約為10-10km/s2量級，因此對動力學模型添加10-9km/s2量級的噪聲。

最小二乘法的位置估計結果如圖5所示。單獨以徑向速度值為量測量進行卡爾曼濾波的結果如圖7所示，信息融合的濾波結果如圖8所示。

由于單獨以徑向速度值作為量測量的系統是不可觀的，因此圖7中位置估計誤差無法收斂至0附近。而多次仿真證明，僅使用該方法時最終收斂到的結果是一個無法預測的數值。

由圖8可以看出，信息融合濾波結果中位置估計誤差在80 km以內，遠優于單獨使用最小二乘法和卡爾曼濾波的結果。速度估計精度優于1.2×10-4km/s，與單獨使用卡爾曼濾波的結果相近。

我們也分別做了提高小行星視線矢量與徑向速度量測精度的仿真試驗，都得到了更高精度的試驗結果。

圖7 卡爾曼濾波結果Fig.7 Results of Kalman filtering

圖8 信息融合濾波結果Fig.8 Results of information fusion filtering

5 結束語

本文對深空探測航天器巡航段的自主導航方法開展了研究，提出了一種結合相對太陽徑向速度測量與小行星視線測量的導航方法和數據融合算法。在理論上分析了系統的能觀性，定量計算了能觀度，給出了小行星的選取建議。以實際工程可行的量測精度和動力學模型精度為輸入條件，開展了仿真分析，結果表明位置估計誤差為80 km(3σ)，速度估計誤差為0.12 m/s(3σ)。該方法計算簡單，精度可滿足巡航段的工程需求，具有很強的實用性。本文只對導航小行星提出了部分選取建議，對于具體的導航小行星選擇未做深入研究，該方面研究將作為本課題的未來工作展開。