黃靜琪，何雨帆，孫山鵬

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

0 引 言

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

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

分布式狀態(tài)估計方法主要分為兩類：基于一致性算法和基于融合算法。對于一致性算法，Olfati-Saber在文獻(xiàn)[7]中將平均一致性算子加入節(jié)點(diǎn)的局部卡爾曼濾波中，提出了卡爾曼一致性濾波(Kalman consistency filter,KCF)算法，該算法要求每個傳感器節(jié)點(diǎn)和鄰居節(jié)點(diǎn)對目標(biāo)有聯(lián)合可觀性。隨后，Olfati-Saber在文獻(xiàn)[4]中首次提出了基于狀態(tài)一致性的最優(yōu)分布式卡爾曼濾波，但該方法要求每個節(jié)點(diǎn)已知鄰居的鄰居節(jié)點(diǎn)信息，且計算量大，無法應(yīng)用到實(shí)際中。為了降低該方法的通信和計算復(fù)雜度，文獻(xiàn)[4]中給出了一種次優(yōu)的分布式卡爾曼濾波器。以文獻(xiàn)[4]為基礎(chǔ)，針對不同問題發(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)，只能用來解決文中提到的兩種特定場景下的分布式濾波問題。

上述分布式算法針對不同的問題提出，常常有特定的假設(shè)，對節(jié)點(diǎn)及其鄰居節(jié)點(diǎn)的聯(lián)合可觀性、通信拓?fù)涞膹?qiáng)連通性、動力學(xué)模型和通信拓?fù)涞臅r不變性等都有嚴(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é)點(diǎn)的估計值逼近全局最優(yōu)。

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

1 信息濾波算法

1.1 擴(kuò)展信息濾波

在實(shí)施擴(kuò)展卡爾曼濾波(Extended Kalman filter,EKF)的過程中，若迭代更新的不是協(xié)方差矩陣P，而是協(xié)方差矩陣P的逆矩陣，則該濾波過程稱作擴(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時刻的狀態(tài)向量;z(k)為k時刻傳感器的測量值。w(k)和v(k)是零均值不相關(guān)白噪聲，且協(xié)方差矩陣分別為Q(k)和R(k)。

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

(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 多傳感器的集中式信息濾波

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

(11)

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

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

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

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

DHIF算法如下：

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

(17)

(18)

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

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

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

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

(19)

(20)

(21)

(22)

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

3.1 動力學(xué)模型

(23)

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

(24)

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

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 量測模型

(28)

(29)

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

(30)

(31)

(32)

(33)

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

(34)

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

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

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

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

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

3)步驟3:獲取量測值

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

(35)

(36)

(37)

(38)

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

7)步驟7:量測更新

(39)

(40)

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

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

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

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

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

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

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

4.1 測站測量方式相同

設(shè)4個測站均采用統(tǒng)一S頻段系統(tǒng)(USB)設(shè)備測量，觀測數(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)

量測噪聲方差矩陣為

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

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

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

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

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

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

4.2 測站測量方式不同

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

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

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

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

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

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

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

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

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

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

5 結(jié) 論

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