日地L2點航天器編隊的分布式自主相對導航

葉子鵬，周慶瑞，王輝

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

日地L2點位于日地動力學平衡點，是進行天文觀測與研究的理想位置[1-2]。在日地L2點部署探測器編隊可實現高精度的光干涉測量，是深空探測的重要手段，如美國的TPF計劃[3-4]。L2點的深空探測編隊對相對控制和導航精度提出了很高要求。

對于地球軌道航天器編隊的相對導航，采用全球導航衛星系統(Global Navigation Satellite System，GNSS)與慣性導航系統(Inertial Navigation System，INS)結合的導航方法通常是最便捷有效的[5-7]，但并不適用于深空探測任務，一種工程可行的方法是通過測量航天器間的距離與方向信息，估計出航天器的相對狀態信息。Kim[8]與Alonso[9]等研究了基于視線矢量量測的航天器編隊相對導航，但是對視線矢量的量測需要光學傳感器，在有限的星載資源下，所能量測的視線矢量數量有限，量測過程也較為復雜。Kuang等[10-11]提出了一種通過測距進行航天器編隊相對導航的方法，該方法在工程上更可行，但未深入研究影響濾波精度的要素。Mcloughlin等[12-13]在濾波算法上進行了改進，首先利用自主編隊飛行(Autonomous Formation Flyer, AFF)傳感器測量航天器間的相對距離與方位，再結合改進的加權最小二乘法，對航天器間的相對狀態進行估計，同時對有限資源下AFF傳感器的調度問題進行了研究。Wang等[14-16]對分布式定位算法進行了研究，其算法可以實現高精度定位。他們基于對LEO軌道航天器的定位結果，實現了航天器高精度的協同控制。Ferguson和How[17]等針對編隊飛行相對導航中計算量較大的問題進行了研究。采用了Schmidt-Kalman濾波方法，通過降低所需估計的狀態維度，從而減少了計算量。但是隨著編隊成員數量增加，不同航天器之間的協方差矩陣處理會變得復雜。同時，該算法需要提前設置編隊中航天器的濾波順序，并且所有航天器需要按順序進行信息交互，因此本質上該算法是集中式的。

本文研究了日地L2點的探測器編隊僅通過星間測距進行相對導航的方法。由于實際工程中航天器的慣性姿態易于獲取，比如采用小型化的星敏感器即可獲得很高精度的慣性姿態，因此假設各航天器三軸慣性姿態已知，在此情況下研究相對導航問題，更具有實際價值。

為了避免集中濾波計算量過大，同時考慮在深空環境中保證較高濾波精度，本文提出了一種分布式的相對導航方法。并對影響導航精度的要素進行了分析，對某些可解耦項進行了精確地擬合，研究結果對實際工程具有一定參考價值。

1 分布式自主導航方法

基于AFF的方法可以實現一對地球軌道主從星之間可靠的相對導航。即使量測量較少，量測精度較低，濾波結果也能夠收斂。這是由于濾波結果受到非線性的地球軌道相對運動方程的約束，其運動軌跡不僅要符合量測值，同時也要符合相對運動方程的解。這樣的強約束導致濾波結果容易收斂。日地L2點主從星之間的相對運動方程是線性的，對濾波結果約束較弱，僅通過某一從星與主星之間測距進行相對導航的方法容易導致濾波誤差發散，尤其當從星與主星距離較遠時，可能會產生多組同時符合量測值與相對運動方程解的結果。

傳統的去中心化方法結合所有從星的相對運動方程，每個從星都對全局狀態信息進行估計，能夠有效解決濾波誤差容易發散的問題。但是在傳統的去中心化方法中，當編隊規模較大時，全局狀態維度過高，每一次估計的計算量較大，因此，本文提出了一種僅采用局部測量信息的分布式導航方法，不僅具有較低的狀態估計維度，同時可通過數據融合保持較高的估計精度。

不失一般性，以一個4星編隊為例介紹提出的分布式自主導航方法。在編隊中兩兩之間存在測量關系，其示意圖如圖1(a)所示。

本文主要研究航天器編隊的相對導航問題，假設主星(Chief)的狀態信息已知。在此假設下，如果采用集中式方法，每顆星都要進行18維全局狀態估計，且隨著編隊成員的增加，待估計狀態數線性增加。針對此問題，提出了一種采用局部測量信息的導航方法，如圖1(b)～圖1(d)所示，每顆從星選擇主星和另一顆從星組成三角量測構型，量測值為3顆星相互的測距值。由于主星的狀態信息可視為已知量，故每次只需估計2顆從星的狀態信息。在此模型中，每顆從星的狀態信息都可能被估計多次，可采用信息融合方法獲取最終狀態信息。

在實際場景中，當存在更多的從星時，可采用如下策略進行量測構型生成：每個從星與自己最近的從星，以及主星生成三角量測構型。如果2顆從星互相選擇，則其中一顆選擇其他從星生成量測構型。如從星i與從星j互為最近鄰星，當i選擇了j后，從星j則不能選擇i，只能選擇離自己次近的從星k。以上策略可以保證所有航天器都至少有一次狀態估計。當完成濾波后，每個從星收集關于自身的狀態估計信息，如果信息組數大于等于2組，則對多組信息進行融合，從而完成對自身狀態的估計；如果某航天器收集的關于自身的信息組數小于2組，則該航天器在自身未選擇的航天器中再選一顆最近的，生成量測構型，對自身再進行一次狀態估計，從而得到2組估計值，進行信息融合。

本文的算法可總結如下：

圖1 編隊量測構型示意圖Fig.1 Diagram of formation measurement configuration

步驟1通過通信協商形成三角測量劃分，原則是各個航天器不選擇重復的三角形。

步驟2根據測量構型，分別進行測量和狀態估計。

步驟3對多組估計數據融合，得到最終結果。

步驟4與鄰居的拓撲關系是否改變，如改變轉到步驟 1，如未改變轉到步驟 2。

以圖1所示的構型為例，提出的航天器編隊分布式自主導航的全部過程如圖2所示。每顆從星都得到了多組狀態估計，通過信息融合，可以得到每顆從星較高精度的狀態估計。

圖2 編隊狀態估計流程圖Fig.2 Flow of formation state estimation

2 相對導航的數學模型

地球軌道航天器編隊通常采用CW方程[8,18]作為相對運動方程。但是日地L2點附近的Halo軌道屬于三體軌道[19]，CW方程不適用。由于L2點的重力梯度較小，可假設航天器編隊間的相對運動不受約束，只受牛頓第二定律的作用。

2.1 相對運動方程

令從星i相對于主星的狀態信息為

i=1,2,3

(1)

令間隔時間為ΔT，受牛頓第二定律作用，短時相對勻速運動可描述為

Xi(k+1)=ΦXi(k)+ΓWi(k),i=1,2,3

(2)

式中：Wi(k)為噪聲；Φ與Γ為狀態轉移矩陣，可以表示為

(3)

(4)

其中：I3×3為三階單位陣。令本體坐標系下從星i的轉動速率為ωi，主星的轉動速率為ωc，則從星i相對于主星有

(5)

Mi=Rz(ψi)Ry(θi)Rx(φi)

(6)

根據式(5)和式(6)可以得到從星i相對于主星的歐拉角速率為

(7)

式中：f(φi,θi,ψi)是與從星i相對于主星本體坐標系的歐拉角相關的矩陣，可由式(5)求得。

2.2 量測模型

量測模型的示意圖如圖3所示，每2顆航天器之間通過收發天線進行測距，結合天線的本體系安裝位置，可以建立量測方程。圖中：Rcd為慣性系下從星Deputy 1相對于主星的位置；Rct為主星本體系下發射天線的安裝位置；Rdr為從星Deputy 1本體系下某一接收天線的安裝位置；Rct與Rdr均為已知量；Rtr為從星接收天線與主星發射天線在主星本體坐標系下的距離矢量。

量測量為發射天線與接收天線之間的距離，可以得到量測量為

(8)

式中：Mci為主星相對于慣性系的姿態矩陣；Mcd為主星相對于從星的姿態矩陣。其他量測量計算方式同理，不同航天器間相對的姿態矩陣可通過轉換計算得到。

量測量可用式(9)描述：

Z=h(Xform,φform,θform,ψform)+V

(9)

式中：Xform為構成量測幾何中的2顆從星狀態信息；φform、θform、ψform為構成量測幾何中的3顆航天器姿態角信息；V為量測噪聲；h(·)為量測方程。

圖3 量測模型示意圖Fig.3 Diagram of measurement model

3 分布式相對導航的狀態估計算法

根據慣性姿態已知的假設，將姿態角信息作為相對導航的輸入量。

令由Deputy 1，Deputy 2，Chief組成的編隊狀態為

(10)

式中：X1(k)與X2(k)如式(1)所示。由式(2)得到編隊的運動學方程為

(11)

式中：Φ、Γ含義如式(3)和式(4)所示；W12為Deputy 1，Deputy 2運動中的過程噪聲。

令編隊中的航天器姿態矩陣為M1i，M2i，Mci，分別代表Deputy 1，Deputy 2，Chief的本體坐標系相對于慣性系的姿態。由M1i，M2i，Mci可以得到

(12)

式中：M21為Deputy 2本體系相對與Deputy 1本體系的姿態矩陣；M1c為Deputy 1本體系相對與Chief本體系的姿態矩陣；M2c為Deputy 2本體系相對與Chief本體系的姿態矩陣。

因此，結合式(9)，編隊的量測方程可以表示為

Z12=h(X12,M1i,M2i,Mci)+V12

(13)

式中：V12為Deputy 1，Deputy 2，Chief所組成編隊的量測噪聲。

采用無跡卡爾曼濾波[20](UKF)進行狀態估計，UKF處理非線性模型時精度更高，計算量更小。

(14)

(15)

設計信息融合方式為

(16)

(17)

式中：W1為Deputy 1的過程噪聲。Deputy 1的整個濾波過程可以用圖4流程圖表示。Deputy 2和Deputy 3的計算過程與Deputy 1相同。

圖4 Deputy 1濾波算法流程示意圖Fig.4 Flow of Deputy 1 filtering algorithm

4 算法的能觀性與精度分析

本節結合仿真實驗，對導航算法的能觀性以及對影響導航精度的因素進行了分析。

4.1 仿真條件

在仿真實驗與分析中，將過程噪聲設為1 mm/s2(1σ)，將三軸姿態角估計精度設置為10″，(1σ)，分別將量測噪聲設置為5 mm(1σ)和1 mm(1σ)。 仿真中濾波周期為1 s，共仿真10 000個步長。

同時，本文對安裝天線的位置進行了分析，令天線安裝示意圖如圖5所示。

令衛星模型為立方體，邊長為R，則接收天線在本體系下的安裝位置分別為

圖5 天線安裝示意圖Fig.5 Diagram of antenna installation

(18)

發射天線在本體系下的安裝位置為

(19)

在仿真分析時，通過改變R的值，從而改變各測量點之間的基線距離。同時，通過減少接收天線數量，研究量測值數量對相對導航精度的影響。

在靜態條件下，當航天器的慣性姿態已知時，在式(8)中，唯一的未知項為三維列向量Rcd，因此求解該三維向量至少需要3個量測值。而在本文所使用的量測構型中，靜態條件下共有2個航天器共六維未知項，因此至少需要6個量測值。即在量測構型的編隊中，每2顆航天器之間只要有2個量測值，就可以實現靜態條件下的相對位置確定。因此，在以下仿真分析中選擇航天器兩兩間產生的量測值個數大于2。

在仿真中分別令Chief，Deputy 1，Deputy 2，Deputy 3的旋轉角速度為

(20)

令Deputy 1，Deputy 2，Deputy 3相對于主星的位置、速度初始值為

(21)

仿真中令初始位置誤差為5 m，初始速度誤差為0.01 m/s。

4.2 仿真過程

為驗證在量測構型中每對航天器之間僅有2個量測量的編隊相對導航情況，設置量測精度為1 mm(1σ)，設置R為12 m，進行濾波得到在第10 000個仿真步長時3顆從星相對位置的估計誤差均方差如表1所示。可以看出該條件下雖然可以實現相對導航，但是相對導航的精度較低。

表1 從星相對位置誤差均方差Table 1 Standard deviation of deputies

令量測精度為1 mm(1σ)，R為12 m，每個航天器安裝4根接收天線，每2顆衛星之間產生8個量測量，濾波結果和局部放大圖如圖6所示。由圖6可以得到，位置估計誤差大約為0.05 m， 速度估計誤差約為5 mm/s。隨著時間增加，由于航天器之間的距離增加，濾波精度下降。

圖6 量測噪聲1 mm時的從星估計誤差Fig.6 Deputies estimation error with 1 mm measurement noise

在實際編隊中，航天器之間的距離通常保持在100 m到1 000 m之間。仿真中，在初始相對速度與噪聲的影響下，航天器間的距離變化如圖7所示。

圖7 從星與主星間的距離變化Fig.7 Distance changing between deputies and chief

第10 000個仿真步長時航天器間的距離如表2所示。可以得到仿真模型與實際編隊情形相似，因此圖6中的濾波精度具有參考價值。

令量測精度為5 mm(1σ)，每個航天器安裝4根接收天線，濾波結果和局部放大圖如圖8所示由圖8可以得到，位置估計誤差大約為0.1 m，速度估計誤差大約為10 mm/s，相較于圖6，濾波精度有所下降。

分別在1 mm(1σ)和5 mm(1σ)量測精度下，設置不同的天線距離與不同數量的接收天線，3顆從星相對位置估計的標準差分別如表3、表4和圖9、圖10所示，其中DN表示每2顆航天器之間的量測值個數，R表示天線之間的距離，表中數組(a,b,c)分別表示Deputy 1，Deputy 2，Deputy 3的位置估計均方差。

圖8 量測噪聲5 mm時的從星估計誤差Fig.8 Deputies estimation error with 5 mm measurement noise

表2 第10 000個仿真步長時航天器間距Table 2 Distance between spacecraft in step 10 000

4.3 仿真分析

從仿真結果可以看出，隨著編隊中航天器間的量測數據增加，航天器上測量點間基線距離增加，濾波精度將會增加，同時量測精度直接影響狀態估計精度。

在文中仿真條件下，當編隊規模在1 000 m以內，量測精度達到1 mm(1σ)時，相位位置估計的均方誤差通常可以達到厘米量級。當量測精度為5 mm(1σ)時，相對位置估計的均方誤差通常在0.1 m到0.3 m之間。

同時，對大量仿真數據進行分析得到，接收天線間的距離R與量測值量測精度為濾波誤差的可解耦項。其他條件保持不變，僅改變接收天線間的距離R，可以得到如表5所示的結果。

表5中n1為相應R值下的濾波誤差與R值為1 m下的濾波誤差的比值。通過擬合，可以得到

(22)

表5中的原始數據和擬合數據的曲線圖如圖11所示。

其他條件保持不變，僅改變量測精度，可以得到如表6所示的結果。

表6中σ為量測噪聲標準差，n2為相應σ值下的濾波誤差與σ值為1 mm下的濾波誤差的比值。

表3 量測噪聲1 mm時位置估計均方差Table 3 Std. deviation with 1 mm measurement noise

表4 量測噪聲5 mm時位置估計均方差Table 4 Std. deviation with 5 mm measurement noise

通過擬合，可以得到

n2=σ0.76

(23)

表6中的原始數據和擬合數據的曲線圖如圖12所示。

結合式(22)和式(23)，可以得到濾波誤差同天線距離R及量測噪聲的標準差σ之間的關系為

(24)

圖9 量測噪聲1 mm時估計均方差曲線圖Fig.9 Diagram of std. deviation with 1 mm measurement noise

圖10 量測噪聲5 mm時估計均方差曲線圖Fig.10 Diagram of std. deviation with 5 mm measurement noise

表5 不同天線距離下的濾波誤差比例Table 5 Proportion under different antenna distances

圖11 濾波誤差與天線距離的關系擬合圖Fig.11 Relationship between filtering error and antenna distance

表6 不同量測噪聲下的濾波誤差比例Table 6 Proportion with different measurement noise

圖12 濾波誤差與量測噪聲的關系擬合圖Fig.12 Relationship between filtering error and measurement noise

式中：ε為估計誤差；F(others)為與其他項(如衛星運動狀態，天線數量，天線安裝構型)有關的數值。通過式(24)，可以快速求出在其他項不變時，不同量測精度及不同天線距離下的估計誤差。

5 結 論

研究了日地L2點深空探測航天器編隊的相對導航問題，提出了一種基于局部測量信息的分布式自主導航方法，通過生成測量構型、濾波、融合3個步驟來確定從星相對于主星的狀態信息。只采用局部測量信息降低了估計狀態的維度，提高了計算效率，同時通過信息的融合提高狀態估計的精度。還對影響相對導航精度的因素進行了研究分析，得到以下結論：

1) 顯然量測值的數量越多，濾波誤差越小，相對導航精度越高。但誤差的下限值受多方面因素影響，當每對航天器間的測量量大于4個的時候，對于精度的提升不明顯，因此，使用本文的方法進行相對導航時，建議安裝4個接收天線。

2) 在相同安裝構型下，接收天線間的距離越遠，即測量點間的基線越大，會有更高的相對導航精度，基線距離為濾波誤差的可解耦項，文中給出了擬合經驗公式。

3) 隨著編隊中航天器的距離增加，相對導航精度將會下降。當航天器之間的距離在500 m以內時，通常可以保持0.01 m(1σ)左右的位置估計誤差；當航天器之間的距離增加到1 000 m左右時，位置估計誤差逐漸逐漸上升至0.03 m(1σ)。

4) 測距精度將直接影響最終的相對狀態估計精度，測距精度為濾波誤差的可解耦項，文中給出擬合經驗公式。