一種分布式信息融合的航天器實時軌道確定算法

黃靜琪，何雨帆，孫山鵬

(1. 西安衛星測控中心宇航動力學國家重點實驗室，西安 710043；2. 火箭軍工程大學，西安 710025)

0 引 言

航天器實時定軌在空間態勢感知、航天器實時跟蹤監視、碰撞預警、航天器機動效果快速評估等方面有著重要作用。實時定軌的測量數據一般由多個相同或者不同類型的傳感器組成的傳感器網絡提供，因此如何融合多傳感器的測量數據快速準確地進行軌道確定十分關鍵。

目前，我國航天測控網主要采用集中式濾波算法[1-3]進行實時定軌，各測量數據均傳輸至軌道計算中心進行融合計算。隨著航天技術應用的不斷發展，對實時定軌系統可靠性的要求也越來越高。集中式算法過分依靠中心節點使的魯棒性差。同時，各傳感器測量數據均要傳輸至中心節點，導致通信量大、中心節點計算負擔高、擴展性較差，分布式算法可有效解決以上問題[4-6]。

分布式狀態估計方法主要分為兩類：基于一致性算法和基于融合算法。對于一致性算法，Olfati-Saber在文獻[7]中將平均一致性算子加入節點的局部卡爾曼濾波中，提出了卡爾曼一致性濾波(Kalman consistency filter,KCF)算法，該算法要求每個傳感器節點和鄰居節點對目標有聯合可觀性。隨后，Olfati-Saber在文獻[4]中首次提出了基于狀態一致性的最優分布式卡爾曼濾波，但該方法要求每個節點已知鄰居的鄰居節點信息，且計算量大，無法應用到實際中。為了降低該方法的通信和計算復雜度，文獻[4]中給出了一種次優的分布式卡爾曼濾波器。以文獻[4]為基礎，針對不同問題發展出許多分布式卡爾曼濾波算法。上述一致性算法通常要求各節點在局部更新前與鄰居節點進行多次通信。除了一致性算法外，文獻[8]提出了基于融合方法的卡爾曼濾波(Distributed Kalman filtering,DKF)，文獻[9]將協方差交叉算法[10](Consensus-plus-innovations,CI)與DKF算法結合，提出了一種基于CI算法的DKF算法(CI-DKF)，但該方法普適性不強，只能用來解決文中提到的兩種特定場景下的分布式濾波問題。

上述分布式算法針對不同的問題提出，常常有特定的假設，對節點及其鄰居節點的聯合可觀性、通信拓撲的強連通性、動力學模型和通信拓撲的時不變性等都有嚴格的要求，有的還要求局部節點需要知道與通信拓撲結構相關的全局信息，不利于實際應用。文獻[11]在上述研究的基礎上，提出了一種分布式混合信息濾波算法(Distributed hybrid information filtering,DHIF)，該算法以擴展信息濾波為基礎，與協方差交叉算法(CI)結合，是一種普適性較強的分布式濾波算法，通信代價低，并且可以使每個節點的估計值逼近全局最優。

本文以DHIF算法為基礎，結合我國航天器測控系統現狀，設計了一種分布式實時軌道確定算法，并對測量網絡中不同傳感器類型和不同網絡拓撲結構進行了仿真計算，結果表明該算法使得各測站定軌精度均優于單站濾波精度，并逼近集中式濾波定軌精度。最后，本文分析了該分布式實時定軌算法收斂速度與通信拓撲的關系。

1 信息濾波算法

1.1 擴展信息濾波

在實施擴展卡爾曼濾波(Extended Kalman filter,EKF)的過程中，若迭代更新的不是協方差矩陣P，而是協方差矩陣P的逆矩陣，則該濾波過程稱作擴展信息濾波(Extended information filter,EIF)。

設有非線性系統：

x(k)=f(x(k-1))+w(k-1)

(1)

z(k)=h(x(k))+v(k)

(2)

式中:x(k)為k時刻的狀態向量;z(k)為k時刻傳感器的測量值。w(k)和v(k)是零均值不相關白噪聲，且協方差矩陣分別為Q(k)和R(k)。

1) 狀態預測

(3)

Pk|k-1=F(k-1)Pk-1|k-1(F(k-1))T+Q(k-1)

(4)

Ik|k-1=(Pk|k-1)-1

(5)

2) 狀態更新

Ik|k=Ik|k-1+(H(k))T(R(k))-1H(k)

(6)

K(k)=(Ik|k)-1(H(k))T(R(k))-1

(7)

(8)

(9)

式(9)乘Ik|k后，將(6)代入式(9)右側，整理后可得：

(10)

1.2 多傳感器的集中式信息濾波

對于多傳感器的集中式信息融合濾波，每個傳感器節點有各自的測量方程，k時刻第i個節點的量測方程可以表示為：

(11)

集中式信息濾波算法如下：

1)信息狀態預測

狀態、協方差矩陣和信息矩陣的預測同式(3)、(4)、(5)。信息狀態為：

(12)

2)信息狀態更新

根據文獻[2]的集中式融合方法，系統信息狀態及信息矩陣的更新為：

(13)

(14)

其中,i=1，…，N，N為網絡中傳感器的數量

(15)

(16)

2 基于傳感器網絡的分布式濾波

在航天器實時軌道確定這類問題中，往往有多個測站，各測站可能使用不同類型的傳感器對目標進行跟蹤測量，可以歸為基于傳感器網絡的融合濾波定軌問題。目前航天器實時軌道確定系統中使用的集中式融合濾波方法依賴于計算中心，一方面需要各傳感器向中心傳遞大量數據，增加了通信負載和計算量，另一方面中心節點故障將直接導致整個系統癱瘓。使用分布式濾波方法定軌，各測站既是測量節點也是計算節點，均可以對目標狀態進行估計。

現有的大多數分布式算法中，針對不同背景問題，有較為嚴格的應用條件。例如，單個節點需要知道某些全局信息(整個網絡中的節點數量或通信拓撲的最大入度等)、通信拓撲不能發生改變、每個節點和鄰居節點是聯合可觀的等。DHIF算法則克服了上述缺點：1)每次局部狀態更新前，鄰居節點之間只需進行一次信息交換，通信量低；2)節點間交換的信息只有局部狀態估計及量測信息，單個節點無需知道整個網絡的全局信息，方便隨時增減傳感器節點，實現了傳感器在網絡中的即插即用；3)DHIF算法應用條件寬松，不需要每個節點及其鄰居都是聯合可觀的，只要滿足傳感器網絡是聯合可觀的，即可保證各節點的局部估計誤差有統一下界，并且誤差的期望值漸進為零[10]；4)算法支持時變的通信拓撲。

DHIF算法如下：

1) 局部狀態預測

(17)

(18)

2) 局部狀態更新

在k時刻，節點i將自身以及接收到的鄰居節點的信息融合進行局部狀態更新。為了保證各節點只通過接收鄰居節點信息就可以獲取到網絡中其他節點的信息，網絡的通信拓撲要求為連通圖。

融合的信息包括兩部分，第一部分為自身及鄰居節點信息矩陣和信息狀態的局部先驗估計，第二部分為自身及鄰居節點的y(k)和s(k)值。

首先，第一部分的融合采用了CI[10]算法，如下：

(19)

(20)

(21)

(22)

3 分布式濾波實時軌道確定

3.1 動力學模型

(23)

式中：μ為地球引力常數，Re為地球半徑。式(23)可以表示為如下形式

(24)

將其離散化，其他攝動看作系統噪聲w，T為離散時間間隔，可得

x(k)=x(k-1)+A(x(k-1))T+w(k-1)

(25)

式(25)可寫為狀態方程(1)的形式，線性化可得狀態轉移矩陣為

(26)

其中，I為3×3的單位矩陣，且

(27)

3.2 量測模型

(28)

(29)

在測站坐標系中有：

(30)

(31)

(32)

(33)

式(30)～(33)可以表示為量測方程(2)的形式，線性化后可得量測矩陣為：

(34)

3.3 分布式實時軌道確定算法

在分布式實時軌道確定中，所有測站構成一個傳感器網絡，每個觀測站看作一個節點，算法采用DHIF，各測站分別對航天器進行實時定軌。計算流程如下：

圖1 分布式實時軌道確定算法流程Fig.1 Distributed real-time orbit determination algorithm

1)步驟1:測站狀態初始化

2)步驟2:狀態預測

3)步驟3:獲取量測值

4)步驟4:測站i計算狀態信息和量測信息

(35)

(36)

(37)

(38)

5)步驟5:測站信息交換

7)步驟7:量測更新

(39)

(40)

4 仿真結果與分析

本節對分布式實時定軌算法進行仿真計算。假設有4個測站構成的傳感器網絡對一航天器進行分布式實時軌道確定。航天器初始軌道見表1，定軌弧段為2007-7-1 12∶32∶19至2007-7-1 12∶43∶59，使用STK高精度軌道外推模型(HPOP)將初軌外推700 s數據作為精密星歷，觀測站坐標見表2，各地面站跟蹤航天器弧段見表3，可以看出航天器不同時對所有測站可見。

表1 目標航天器初始軌道Table 1 Target spacecraft initial orbit

表2 測站坐標Table 2 Coordinate of the observation stations

表3 各測站跟蹤弧段Table 3 The observation arc of each station s

計算采用200次蒙特卡洛仿真，分別使用圖2兩種通信拓撲結構計算。

圖2 兩種不同的通信拓撲結構Fig.2 Two different communication topologies

4.1 測站測量方式相同

設4個測站均采用統一S頻段系統(USB)設備測量，觀測數據間隔1 s。

設初始狀態為

初始狀態協方差矩陣為

系統噪聲方差矩陣為

Q=diag(0.12,0.12,0.12,0.012,0.012,0.012)

量測噪聲方差矩陣為

Ri=diag(102,0.052,0.022,0.022)

試驗結果見圖3～圖4。其中，圖3 (a)、(b)分別為在圖2兩種通信拓撲下，各測站80 s內分布式濾波定軌位置和速度的均方根誤差(Root mean squared error,RMSE)[1]曲線。

圖3 各測站80 s內分布式濾波定軌的位置、速度RMSE曲線Fig.3 80 s position and velocity RMSE of DHIF for each station

結果表明，在初軌誤差較大的情況下(位置誤差13.69 km，速度誤差52.1 m/s)，該分布式濾波定軌算法可以使各測站的計算結果快速收斂，且定軌結果趨于一致。根據圖3，在前65 s只有測站S1、S2有量測數據(見表3)的情況下，各測站50 s定軌位置誤差達到200 m，速度誤差下降到10 m/s。在65 s時，測站S3開始跟蹤目標，此時傳感器網絡中三個測站具有量測數據，這些量測信息通過相鄰測站之間的通信在網絡中傳播。從RMSE曲線斜率的變化趨勢可以看出，各測站定軌收斂速度加快。80 s時，定軌位置誤差達到10 m以內，速度誤差在0.1 m/s以內。

從式(39)、(40)可以看出，對某一測站，其相鄰測站的狀態信息和量測信息直接影響該測站的濾波結果，不同的通信拓撲會導致測站的鄰居節點發生變化，繼而影響該測站計算收斂的速度。對比圖3(a)、(b)，在各測站觀測弧段相同的情況下，拓撲a的收斂速度要快于拓撲b。

圖4為各測站使用不同濾波方法在整個定軌弧段的RMSE曲線。圖4中，DHIF-a，DHIF-b為拓撲a，b下的分布式濾波，IF為測站單站信息濾波，CIF為4個測站集中式信息濾波。表4和表5列出了集中式濾波和分布式濾波RMSE的均值。

結果表明，使用該分布式算法，各測站的定軌精度均優于單站信息濾波定軌精度，且在不同的通信拓撲下均逼近于集中式信息濾波定軌精度。該算法中，測站每次狀態更新前與相鄰測站進行信息傳遞，單個測站融合了其他測站的信息，因此定軌結果優于單站信息濾波。同時，雖然不像集中式算法計算中心可以直接獲取所有測站的測量值，但通信拓撲的連通性使得各個測站間接獲得了其他測量信息，得到了與集中式定軌精度相當的結果。

4.2 測站測量方式不同

設4個測站采用不同測量設備，見表6，不同測量設備的觀測量見表7。

R1=R4=diag(202,0.022,0.022)，

R2=diag(0.0022,0.0022)，

R3=diag(102,0.052,0.022,0.022)。

圖5為兩種拓撲下，各測站80 s內分布式濾波定軌位置和速度的RMSE曲線。

表5 各測站集中式與分布式濾波定軌速度平均RMSETable 5 Velocity average RMSE of DHIF and CIF for each station

表6 各測站測量設備Table 6 Measure equipment of each observation station

80 s時，定軌位置誤差在10 m以內，速度誤差在0.1 m/s以內。圖6為各測站不同拓撲結構下定軌比較。

表7 不同設備的觀測量Table 7 The observation type of different equipment

圖5 各測站80 s內分布式濾波定軌的位置、速度RMSE曲線Fig.5 80 s position and velocity RMSE of DHIF for each station

圖6(a)、(b)中，測站S1和S2在拓撲b下收斂速度快于拓撲a。這是因為拓撲b的結構使得測站S1和S2的1和2類信息能更快的傳輸至整個網絡中。S1和S2在起始階段，收到更多的是1和2類信息。而在拓撲a中，起始階段S1和S2很快收到了來自S3和S4的第3類信息，特別對于S2，開始幾秒只接收到了第3類信息，因此收斂速度最慢。

圖6(c)、(d)中，結果則相反，測站S3，S4在拓撲a下收斂速度快于拓撲b。因為雖然通信拓撲a的信息交換速度比拓撲b慢，但拓撲a中的S3和S4在前幾秒收到的更多是第1和2類信息，尤其S4，開始只收到了第1和2類信息，因此很快收斂。而拓撲b中的S3和S4開始時收到的1和2類信息少，第3類信息多，因此收斂速度相對較慢。

圖5也驗證了上述結論，圖5 (a)中，根據拓撲a中信息的傳遞方向，收斂速度最快的為S4，最慢的為S2。圖5(b)中，因為拓撲b使得各類信息在整個網絡中可以快速傳遞，四個測站的收斂速度基本接近。

圖6 測站在不同拓撲下80 s內分布式實時定軌位置、速度RMSE曲線對比Fig.6 Comparison of 80 s position and velocity RMSE for each station in different topologies

5 結 論

針對我國航天器測控系統，本文結合DHIF設計了一種分布式實時軌道確定算法。該算法使得單站分布式定軌精度優于單站濾波定軌精度并逼近集中式濾波方法，可以擴展至天地基多種傳感器聯合觀測的分布式實時定軌中。該算法要求各測站的量測數據在同一時刻，在實際工程中，由于不同傳感器的量測模型和數據修正模型不同，存在不同測站修正后的量測數據不在同一時刻的問題，因此如何進行多傳感器異步數據的分布式融合，是需進一步研究的課題。