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

黃靜琪，何雨帆，孫山鵬

(1. 西安衛(wèi)星測(cè)控中心宇航動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，西安 710043；2. 火箭軍工程大學(xué)，西安 710025)

0 引 言

航天器實(shí)時(shí)定軌在空間態(tài)勢(shì)感知、航天器實(shí)時(shí)跟蹤監(jiān)視、碰撞預(yù)警、航天器機(jī)動(dòng)效果快速評(píng)估等方面有著重要作用。實(shí)時(shí)定軌的測(cè)量數(shù)據(jù)一般由多個(gè)相同或者不同類型的傳感器組成的傳感器網(wǎng)絡(luò)提供，因此如何融合多傳感器的測(cè)量數(shù)據(jù)快速準(zhǔn)確地進(jìn)行軌道確定十分關(guān)鍵。

目前，我國(guó)航天測(cè)控網(wǎng)主要采用集中式濾波算法[1-3]進(jìn)行實(shí)時(shí)定軌，各測(cè)量數(shù)據(jù)均傳輸至軌道計(jì)算中心進(jìn)行融合計(jì)算。隨著航天技術(shù)應(yīng)用的不斷發(fā)展，對(duì)實(shí)時(shí)定軌系統(tǒng)可靠性的要求也越來(lái)越高。集中式算法過(guò)分依靠中心節(jié)點(diǎn)使的魯棒性差。同時(shí)，各傳感器測(cè)量數(shù)據(jù)均要傳輸至中心節(jié)點(diǎn)，導(dǎo)致通信量大、中心節(jié)點(diǎn)計(jì)算負(fù)擔(dān)高、擴(kuò)展性較差，分布式算法可有效解決以上問(wèn)題[4-6]。

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

上述分布式算法針對(duì)不同的問(wèn)題提出，常常有特定的假設(shè)，對(duì)節(jié)點(diǎn)及其鄰居節(jié)點(diǎn)的聯(lián)合可觀性、通信拓?fù)涞膹?qiáng)連通性、動(dòng)力學(xué)模型和通信拓?fù)涞臅r(shí)不變性等都有嚴(yán)格的要求，有的還要求局部節(jié)點(diǎn)需要知道與通信拓?fù)浣Y(jié)構(gòu)相關(guān)的全局信息，不利于實(shí)際應(yīng)用。文獻(xiàn)[11]在上述研究的基礎(chǔ)上，提出了一種分布式混合信息濾波算法(Distributed hybrid information filtering,DHIF)，該算法以擴(kuò)展信息濾波為基礎(chǔ)，與協(xié)方差交叉算法(CI)結(jié)合，是一種普適性較強(qiáng)的分布式濾波算法，通信代價(jià)低，并且可以使每個(gè)節(jié)點(diǎn)的估計(jì)值逼近全局最優(yōu)。

本文以DHIF算法為基礎(chǔ)，結(jié)合我國(guó)航天器測(cè)控系統(tǒng)現(xiàn)狀，設(shè)計(jì)了一種分布式實(shí)時(shí)軌道確定算法，并對(duì)測(cè)量網(wǎng)絡(luò)中不同傳感器類型和不同網(wǎng)絡(luò)拓?fù)浣Y(jié)構(gòu)進(jìn)行了仿真計(jì)算，結(jié)果表明該算法使得各測(cè)站定軌精度均優(yōu)于單站濾波精度，并逼近集中式濾波定軌精度。最后，本文分析了該分布式實(shí)時(shí)定軌算法收斂速度與通信拓?fù)涞年P(guān)系。

1 信息濾波算法

1.1 擴(kuò)展信息濾波

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

設(shè)有非線性系統(tǒng)：

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

(1)

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

(2)

式中:x(k)為k時(shí)刻的狀態(tài)向量;z(k)為k時(shí)刻傳感器的測(cè)量值。w(k)和v(k)是零均值不相關(guān)白噪聲，且協(xié)方差矩陣分別為Q(k)和R(k)。

1) 狀態(tài)預(yù)測(cè)

(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) 狀態(tài)更新

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)右側(cè)，整理后可得：

(10)

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

對(duì)于多傳感器的集中式信息融合濾波，每個(gè)傳感器節(jié)點(diǎn)有各自的測(cè)量方程，k時(shí)刻第i個(gè)節(jié)點(diǎn)的量測(cè)方程可以表示為：

(11)

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

1)信息狀態(tài)預(yù)測(cè)

狀態(tài)、協(xié)方差矩陣和信息矩陣的預(yù)測(cè)同式(3)、(4)、(5)。信息狀態(tài)為：

(12)

2)信息狀態(tài)更新

根據(jù)文獻(xiàn)[2]的集中式融合方法，系統(tǒng)信息狀態(tài)及信息矩陣的更新為：

(13)

(14)

其中,i=1，…，N，N為網(wǎng)絡(luò)中傳感器的數(shù)量

(15)

(16)

2 基于傳感器網(wǎng)絡(luò)的分布式濾波

在航天器實(shí)時(shí)軌道確定這類問(wèn)題中，往往有多個(gè)測(cè)站，各測(cè)站可能使用不同類型的傳感器對(duì)目標(biāo)進(jìn)行跟蹤測(cè)量，可以歸為基于傳感器網(wǎng)絡(luò)的融合濾波定軌問(wèn)題。目前航天器實(shí)時(shí)軌道確定系統(tǒng)中使用的集中式融合濾波方法依賴于計(jì)算中心，一方面需要各傳感器向中心傳遞大量數(shù)據(jù)，增加了通信負(fù)載和計(jì)算量，另一方面中心節(jié)點(diǎn)故障將直接導(dǎo)致整個(gè)系統(tǒng)癱瘓。使用分布式濾波方法定軌，各測(cè)站既是測(cè)量節(jié)點(diǎn)也是計(jì)算節(jié)點(diǎn)，均可以對(duì)目標(biāo)狀態(tài)進(jìn)行估計(jì)。

現(xiàn)有的大多數(shù)分布式算法中，針對(duì)不同背景問(wèn)題，有較為嚴(yán)格的應(yīng)用條件。例如，單個(gè)節(jié)點(diǎn)需要知道某些全局信息(整個(gè)網(wǎng)絡(luò)中的節(jié)點(diǎn)數(shù)量或通信拓?fù)涞淖畲笕攵鹊?、通信拓?fù)洳荒馨l(fā)生改變、每個(gè)節(jié)點(diǎn)和鄰居節(jié)點(diǎn)是聯(lián)合可觀的等。DHIF算法則克服了上述缺點(diǎn)：1)每次局部狀態(tài)更新前，鄰居節(jié)點(diǎn)之間只需進(jìn)行一次信息交換，通信量低；2)節(jié)點(diǎn)間交換的信息只有局部狀態(tài)估計(jì)及量測(cè)信息，單個(gè)節(jié)點(diǎn)無(wú)需知道整個(gè)網(wǎng)絡(luò)的全局信息，方便隨時(shí)增減傳感器節(jié)點(diǎn)，實(shí)現(xiàn)了傳感器在網(wǎng)絡(luò)中的即插即用；3)DHIF算法應(yīng)用條件寬松，不需要每個(gè)節(jié)點(diǎn)及其鄰居都是聯(lián)合可觀的，只要滿足傳感器網(wǎng)絡(luò)是聯(lián)合可觀的，即可保證各節(jié)點(diǎn)的局部估計(jì)誤差有統(tǒng)一下界，并且誤差的期望值漸進(jìn)為零[10]；4)算法支持時(shí)變的通信拓?fù)洹?/p>

DHIF算法如下：

1) 局部狀態(tài)預(yù)測(cè)

(17)

(18)

2) 局部狀態(tài)更新

在k時(shí)刻，節(jié)點(diǎn)i將自身以及接收到的鄰居節(jié)點(diǎn)的信息融合進(jìn)行局部狀態(tài)更新。為了保證各節(jié)點(diǎn)只通過(guò)接收鄰居節(jié)點(diǎn)信息就可以獲取到網(wǎng)絡(luò)中其他節(jié)點(diǎn)的信息，網(wǎng)絡(luò)的通信拓?fù)湟鬄檫B通圖。

融合的信息包括兩部分，第一部分為自身及鄰居節(jié)點(diǎn)信息矩陣和信息狀態(tài)的局部先驗(yàn)估計(jì)，第二部分為自身及鄰居節(jié)點(diǎn)的y(k)和s(k)值。

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

(19)

(20)

(21)

(22)

3 分布式濾波實(shí)時(shí)軌道確定

3.1 動(dòng)力學(xué)模型

(23)

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

(24)

將其離散化，其他攝動(dòng)看作系統(tǒng)噪聲w，T為離散時(shí)間間隔，可得

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

(25)

式(25)可寫為狀態(tài)方程(1)的形式，線性化可得狀態(tài)轉(zhuǎn)移矩陣為

(26)

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

(27)

3.2 量測(cè)模型

(28)

(29)

在測(cè)站坐標(biāo)系中有：

(30)

(31)

(32)

(33)

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

(34)

3.3 分布式實(shí)時(shí)軌道確定算法

在分布式實(shí)時(shí)軌道確定中，所有測(cè)站構(gòu)成一個(gè)傳感器網(wǎng)絡(luò)，每個(gè)觀測(cè)站看作一個(gè)節(jié)點(diǎn)，算法采用DHIF，各測(cè)站分別對(duì)航天器進(jìn)行實(shí)時(shí)定軌。計(jì)算流程如下：

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

1)步驟1:測(cè)站狀態(tài)初始化

2)步驟2:狀態(tài)預(yù)測(cè)

3)步驟3:獲取量測(cè)值

4)步驟4:測(cè)站i計(jì)算狀態(tài)信息和量測(cè)信息

(35)

(36)

(37)

(38)

5)步驟5:測(cè)站信息交換

7)步驟7:量測(cè)更新

(39)

(40)

4 仿真結(jié)果與分析

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

表1 目標(biāo)航天器初始軌道Table 1 Target spacecraft initial orbit

表2 測(cè)站坐標(biāo)Table 2 Coordinate of the observation stations

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

計(jì)算采用200次蒙特卡洛仿真，分別使用圖2兩種通信拓?fù)浣Y(jié)構(gòu)計(jì)算。

圖2 兩種不同的通信拓?fù)浣Y(jié)構(gòu)Fig.2 Two different communication topologies

4.1 測(cè)站測(cè)量方式相同

設(shè)4個(gè)測(cè)站均采用統(tǒng)一S頻段系統(tǒng)(USB)設(shè)備測(cè)量，觀測(cè)數(shù)據(jù)間隔1 s。

設(shè)初始狀態(tài)為

初始狀態(tài)協(xié)方差矩陣為

系統(tǒng)噪聲方差矩陣為

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

量測(cè)噪聲方差矩陣為

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

試驗(yàn)結(jié)果見圖3～圖4。其中，圖3 (a)、(b)分別為在圖2兩種通信拓?fù)湎拢鳒y(cè)站80 s內(nèi)分布式濾波定軌位置和速度的均方根誤差(Root mean squared error,RMSE)[1]曲線。

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

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

從式(39)、(40)可以看出，對(duì)某一測(cè)站，其相鄰測(cè)站的狀態(tài)信息和量測(cè)信息直接影響該測(cè)站的濾波結(jié)果，不同的通信拓?fù)鋾?huì)導(dǎo)致測(cè)站的鄰居節(jié)點(diǎn)發(fā)生變化，繼而影響該測(cè)站計(jì)算收斂的速度。對(duì)比圖3(a)、(b)，在各測(cè)站觀測(cè)弧段相同的情況下，拓?fù)鋋的收斂速度要快于拓?fù)鋌。

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

結(jié)果表明，使用該分布式算法，各測(cè)站的定軌精度均優(yōu)于單站信息濾波定軌精度，且在不同的通信拓?fù)湎戮平诩惺叫畔V波定軌精度。該算法中，測(cè)站每次狀態(tài)更新前與相鄰測(cè)站進(jìn)行信息傳遞，單個(gè)測(cè)站融合了其他測(cè)站的信息，因此定軌結(jié)果優(yōu)于單站信息濾波。同時(shí)，雖然不像集中式算法計(jì)算中心可以直接獲取所有測(cè)站的測(cè)量值，但通信拓?fù)涞倪B通性使得各個(gè)測(cè)站間接獲得了其他測(cè)量信息，得到了與集中式定軌精度相當(dāng)?shù)慕Y(jié)果。

4.2 測(cè)站測(cè)量方式不同

設(shè)4個(gè)測(cè)站采用不同測(cè)量設(shè)備，見表6，不同測(cè)量設(shè)備的觀測(cè)量見表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為兩種拓?fù)湎拢鳒y(cè)站80 s內(nèi)分布式濾波定軌位置和速度的RMSE曲線。

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

表6 各測(cè)站測(cè)量設(shè)備Table 6 Measure equipment of each observation station

80 s時(shí)，定軌位置誤差在10 m以內(nèi)，速度誤差在0.1 m/s以內(nèi)。圖6為各測(cè)站不同拓?fù)浣Y(jié)構(gòu)下定軌比較。

表7 不同設(shè)備的觀測(cè)量Table 7 The observation type of different equipment

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

圖6(a)、(b)中，測(cè)站S1和S2在拓?fù)鋌下收斂速度快于拓?fù)鋋。這是因?yàn)橥負(fù)鋌的結(jié)構(gòu)使得測(cè)站S1和S2的1和2類信息能更快的傳輸至整個(gè)網(wǎng)絡(luò)中。S1和S2在起始階段，收到更多的是1和2類信息。而在拓?fù)鋋中，起始階段S1和S2很快收到了來(lái)自S3和S4的第3類信息，特別對(duì)于S2，開始幾秒只接收到了第3類信息，因此收斂速度最慢。

圖6(c)、(d)中，結(jié)果則相反，測(cè)站S3，S4在拓?fù)鋋下收斂速度快于拓?fù)鋌。因?yàn)殡m然通信拓?fù)鋋的信息交換速度比拓?fù)鋌慢，但拓?fù)鋋中的S3和S4在前幾秒收到的更多是第1和2類信息，尤其S4，開始只收到了第1和2類信息，因此很快收斂。而拓?fù)鋌中的S3和S4開始時(shí)收到的1和2類信息少，第3類信息多，因此收斂速度相對(duì)較慢。

圖5也驗(yàn)證了上述結(jié)論，圖5 (a)中，根據(jù)拓?fù)鋋中信息的傳遞方向，收斂速度最快的為S4，最慢的為S2。圖5(b)中，因?yàn)橥負(fù)鋌使得各類信息在整個(gè)網(wǎng)絡(luò)中可以快速傳遞，四個(gè)測(cè)站的收斂速度基本接近。

圖6 測(cè)站在不同拓?fù)湎?0 s內(nèi)分布式實(shí)時(shí)定軌位置、速度RMSE曲線對(duì)比Fig.6 Comparison of 80 s position and velocity RMSE for each station in different topologies

5 結(jié) 論

針對(duì)我國(guó)航天器測(cè)控系統(tǒng)，本文結(jié)合DHIF設(shè)計(jì)了一種分布式實(shí)時(shí)軌道確定算法。該算法使得單站分布式定軌精度優(yōu)于單站濾波定軌精度并逼近集中式濾波方法，可以擴(kuò)展至天地基多種傳感器聯(lián)合觀測(cè)的分布式實(shí)時(shí)定軌中。該算法要求各測(cè)站的量測(cè)數(shù)據(jù)在同一時(shí)刻，在實(shí)際工程中，由于不同傳感器的量測(cè)模型和數(shù)據(jù)修正模型不同，存在不同測(cè)站修正后的量測(cè)數(shù)據(jù)不在同一時(shí)刻的問(wèn)題，因此如何進(jìn)行多傳感器異步數(shù)據(jù)的分布式融合，是需進(jìn)一步研究的課題。