姚培東，李星秀，吳盤龍，蘇 鼎

（1.南京理工大學(xué)自動(dòng)化學(xué)院， 南京 210094；2.南京理工大學(xué)理學(xué)院， 南京 210094；3.中國人民解放軍 66483 部隊(duì)， 北京 100093）

0 引言

隨著全球太空資源開發(fā)的提高和未來太空作戰(zhàn)趨勢(shì)的加劇，越來越多的國家擁有了進(jìn)入太空和利用太空的能力，太空變得更加具有競(jìng)爭(zhēng)性和對(duì)抗性。其中，空間監(jiān)視系統(tǒng)可以監(jiān)視空間目標(biāo)的運(yùn)行，確定潛在的威脅，還可以預(yù)測(cè)空間目標(biāo)的軌道，其水平的高低直接制約著空間對(duì)抗能力的發(fā)揮［1］。同時(shí)，空間監(jiān)視系統(tǒng)相比地面監(jiān)視系統(tǒng)具有全天候、視場(chǎng)廣、不受大氣環(huán)境影響等特點(diǎn)，發(fā)展天基監(jiān)視系統(tǒng)就顯得尤為重要。

天基觀測(cè)平臺(tái)上獲取目標(biāo)信息的方式包括有源主動(dòng)和無源被動(dòng)兩種［2］。其中，無源被動(dòng)的工作方式隱蔽性好、作用距離遠(yuǎn)，具有重要的實(shí)用價(jià)值，但是此種方法只能獲得目標(biāo)的角度測(cè)量信息，定位精度不高。為了提高定位精度，一些學(xué)者［3-4］通過引入改進(jìn)的迭代擴(kuò)展Kalman濾波（Iterated Extended Kalman Filter，IEKF）和迭代擴(kuò)展Kalman粒子濾波（Iterated Extended Kalman Particle Filter，IEKPF）來增加迭代次數(shù)。雖然能提高定位精度，但是這種方法計(jì)算量太大，實(shí)時(shí)性很難滿足，當(dāng)測(cè)量誤差較大時(shí)，效果反而更差。文獻(xiàn)［5］通過將初定軌（Initial Orbit Determination，IOD）和序貫估計(jì)相結(jié)合，利用改進(jìn)的Laplace法，通過在一定時(shí)間內(nèi)利用更多組測(cè)量信息提高初定軌的精度用于序貫估計(jì)。但是，這種方法需要在短時(shí)間內(nèi)獲取更多組測(cè)量值，對(duì)系統(tǒng)硬件要求很高，且單個(gè)平臺(tái)對(duì)空間目標(biāo)光學(xué)測(cè)角定軌可觀度和幾何約束較差，當(dāng)目標(biāo)與觀測(cè)星的距離較遠(yuǎn)時(shí)，短弧段的測(cè)量序列隨時(shí)間變化近似線性，Laplace方程接近病態(tài)，造成定軌偏差較大。因此，該方法不適合跟蹤遠(yuǎn)距離的目標(biāo)。為了解決此問題，文獻(xiàn)［6］通過對(duì)目標(biāo)長(zhǎng)時(shí)間的觀測(cè)數(shù)據(jù)采集，選取不同的可觀測(cè)弧段進(jìn)行數(shù)據(jù)融合，間接改善了觀測(cè)幾何條件，使定軌精度提高。但是，這種方法需要長(zhǎng)時(shí)間采集數(shù)據(jù)，實(shí)時(shí)性也比較差，且量測(cè)模型沒有考慮坐標(biāo)轉(zhuǎn)換，雖然公式簡(jiǎn)單，但是實(shí)用性不強(qiáng)。

針對(duì)上述問題，結(jié)合現(xiàn)有分布式衛(wèi)星快速發(fā)展的背景，利用編隊(duì)衛(wèi)星對(duì)目標(biāo)進(jìn)行多視線測(cè)量，改善觀測(cè)的幾何條件，融合觀測(cè)數(shù)據(jù)，在保證實(shí)時(shí)性的同時(shí)，提高了對(duì)非合作目標(biāo)的跟蹤精度和穩(wěn)定性。此外，受空間幾何因素和光學(xué)因素的影響［7］，利用光電傳感器對(duì)目標(biāo)進(jìn)行觀測(cè)時(shí)會(huì)出現(xiàn)某段時(shí)間內(nèi)測(cè)量值丟失和測(cè)量值誤差過大的現(xiàn)象，這些測(cè)量故障會(huì)造成跟蹤精度不高甚至跟蹤結(jié)果發(fā)散?；诖嗽O(shè)計(jì)自適應(yīng)濾波，該濾波基于新息構(gòu)建自適應(yīng)因子，設(shè)計(jì)增益調(diào)節(jié)矩陣，提高了濾波的魯棒性。最后，通過仿真對(duì)提出方法的有效性進(jìn)行了驗(yàn)證和分析。

1 空間目標(biāo)狀態(tài)模型

如圖1所示，在理想情況下，衛(wèi)星和地球組成一個(gè)二體系統(tǒng)，衛(wèi)星沿固定的軌道圍繞地球運(yùn)動(dòng)。但實(shí)際情況下，衛(wèi)星在圍繞地球運(yùn)動(dòng)的過程中會(huì)受到很多干擾，其中最主要的是地球非球形引力攝動(dòng)、大氣阻力攝動(dòng)、日月引力攝動(dòng)、太陽光壓攝動(dòng)等。不同高度軌道的衛(wèi)星所受的攝動(dòng)影響不同，本文選取的非合作目標(biāo)運(yùn)行在高度較低的軌道，因此所受的主要攝動(dòng)為地球非球形引力攝動(dòng)。

圖1 衛(wèi)星編隊(duì)非合作目標(biāo)跟蹤系統(tǒng)Fig.1 Satellite formation for NCST tracking system

在僅考慮地球J2攝動(dòng)時(shí)，衛(wèi)星的軌道方程如下

式（1）中，f為地球非球形J2攝動(dòng)下所產(chǎn)生的加速度。f在地心慣性坐標(biāo)系（Earth Centered Inertial，ECI）下可表示為［6］

則式（2）可改寫為

式（1）～式（5）中，r＝［xyz］T為衛(wèi)星在ECI下的位置矢量；為衛(wèi)星在ECI下的速度矢量；為衛(wèi)星距地心的距離；μ為地球引力常數(shù)，其值為3.986×1014m3/s2；J2為二階帶諧項(xiàng)系數(shù)，其值為1.08263×10-3；Re為地球赤道半徑；k＝［0 0 1］T。

將r和v定義為狀態(tài)矢量，即X＝［r；v］，則式（1）可寫為

式（6）離散化后可進(jìn)一步寫為

利用Euler法對(duì)式（7）進(jìn)行求解

式（8）中，Δt＝tk+1-tk，wk為過程噪聲矩陣，假定Gauss白噪聲均值為0，協(xié)方差為Qk。

2 測(cè)量模型

2.1 觀測(cè)方程中的坐標(biāo)轉(zhuǎn)換

在典型的僅測(cè)角無源定位中，方位角和高低角的計(jì)算定義為目標(biāo)與觀測(cè)器位置矢量之差。以往的文獻(xiàn)多是直接在ECI下建立觀測(cè)方程，公式簡(jiǎn)單，但是觀測(cè)數(shù)據(jù)是在主星和從星各自的軌道坐標(biāo)系下測(cè)量的，而雙星觀測(cè)器和目標(biāo)的位置矢量都建立在ECI下，而且需要建立編隊(duì)雙星之間的坐標(biāo)關(guān)系。因此，需要進(jìn)行坐標(biāo)轉(zhuǎn)換才更接近真實(shí)系統(tǒng)，Euler角的坐標(biāo)轉(zhuǎn)換如圖2所示。

圖2 Euler角坐標(biāo)轉(zhuǎn)換示意圖Fig.2 Schematic diagram of Euler angles coordinate transformation

式（11）中，Gi為坐標(biāo)轉(zhuǎn)換矩陣。以觀測(cè)主星為例，G1＝R3（u）R1（i）R3（Ω），Ω、i、u分別為觀測(cè)主星的升交點(diǎn)赤經(jīng)、軌道傾角和緯度幅角［8］，R3（Ω）、R3（u）為繞z軸旋轉(zhuǎn)角度Ω和u的變換矩陣，R1（i）為繞x軸旋轉(zhuǎn)角度i的變換矩陣，其矩陣形式為

2.2 目標(biāo)運(yùn)動(dòng)觀測(cè)方程

設(shè)觀測(cè)主星對(duì)目標(biāo)測(cè)得的高低角和方位角分別為αk1、βk1，觀測(cè)從星對(duì)目標(biāo)測(cè)得的高低角和方位角分別為αk2、βk2，其定義如下

式（13）中，ηα1、ηα2為高低角的測(cè)量噪聲，ηβ1、ηβ2為方位角的測(cè)量噪聲。

將獲得的方位角和高低角投影到觀測(cè)主星和觀測(cè)從星的軌道坐標(biāo)系下，則測(cè)量值可表示為

從星的測(cè)量值通過星間鏈路傳輸給主星，zk2為目標(biāo)在從星軌道系下的量測(cè)值，應(yīng)轉(zhuǎn)換到主星軌道系下。將zk2從從星軌道系轉(zhuǎn)換到ECI，其坐標(biāo)轉(zhuǎn)換矩陣為；再從ECI轉(zhuǎn)換到主星軌道坐標(biāo)系，其坐標(biāo)轉(zhuǎn)換矩陣為G1。因此，為從星軌道坐標(biāo)系到主星軌道坐標(biāo)系的轉(zhuǎn)換矩陣。最終，雙星編隊(duì)對(duì)非合作目標(biāo)的量測(cè)值為

觀測(cè)雙星對(duì)目標(biāo)的觀測(cè)量均為單位矢量，故量測(cè)方程的右端也應(yīng)為單位矢量。在慣性坐標(biāo)系下，觀測(cè)雙星到NCST的單位觀測(cè)矢量為

式（15）建立在主星軌道坐標(biāo)系下，式（16）建立在ECI下。為了方便濾波時(shí)協(xié)方差矩陣等的獲取，需進(jìn)行坐標(biāo)轉(zhuǎn)換，最終可寫為

式（18）中，v為測(cè)量噪聲矢量。假設(shè)為Gauss白噪聲序列，其均值為0，協(xié)方差為Rk。

3 跟蹤濾波算法

3.1 平方根容積Kalman濾波

天基僅測(cè)角跟蹤是非線性狀態(tài)估計(jì)，平方根容積Kalman濾波（Square Root Cubature Kalman Filter，SRCKF）使用基于容積規(guī)則的數(shù)值積分方法直接計(jì)算非線性函數(shù)的均值和方差，只需計(jì)算函數(shù)值，避免求導(dǎo)計(jì)算。同時(shí)，SRCKF利用協(xié)方差的平方根直接參與遞推運(yùn)算，確保了協(xié)方差矩陣的對(duì)稱性和正定性，能保證濾波的精度和穩(wěn)定性。

SRCKF濾波的步驟如下：

（1）初始化

（2）時(shí)間更新

計(jì)算容積點(diǎn)，有

式（21）中，N＝2m，m為狀態(tài)的維數(shù)；S0＝chol（P0），即。記m維單位向量為e＝［1 0 … 0］T，使用符號(hào)［1］表示對(duì)e的元素進(jìn)行全排列和改變?cè)胤?hào)所產(chǎn)生的點(diǎn)集，稱為完整全對(duì)稱點(diǎn)集。

（3）量測(cè)更新

3.2 自適應(yīng)SRCKF

在復(fù)雜的太空環(huán)境中，由于地球遮擋等幾何因素、目標(biāo)輻射強(qiáng)度以及目標(biāo)背景較亮等影響光電傳感器工作的光學(xué)因素，難免會(huì)出現(xiàn)量測(cè)信息不準(zhǔn)確或者信息傳輸?shù)臅簳r(shí)中斷，導(dǎo)致不能提供正確的量測(cè)信息。此時(shí)，傳統(tǒng)的Kalman濾波精度會(huì)大大降低，甚至導(dǎo)致濾波發(fā)散。因此，需要設(shè)計(jì)自適應(yīng)濾波來提高其魯棒性。自適應(yīng)濾波一般是基于新息和基于殘差設(shè)計(jì)的［9-10］，本文采用的是基于新息的方法來設(shè)計(jì)自適應(yīng)平方根容積Kalman濾波（Adaptive Square Root Cubature Kalman Filter，ASRCKF）。

自適應(yīng)濾波通過采集過去一定時(shí)間的新息，利用采集到的新息實(shí)時(shí)調(diào)整濾波器的誤差矩陣［11-12］。具體調(diào)整策略是通過比較實(shí)時(shí)估計(jì)的矩陣值和設(shè)定的矩陣值，當(dāng)相差大于一定的閾值時(shí)就進(jìn)行調(diào)整。

真實(shí)測(cè)量值和一步預(yù)測(cè)值計(jì)算的殘差公式如下

式（37）中，Zk為k時(shí)刻的真實(shí)測(cè)量值，為k時(shí)刻的量測(cè)估計(jì)值。

若系統(tǒng)量測(cè)的真實(shí)誤差統(tǒng)計(jì)特性與濾波遞推誤差模型一致時(shí)，以下等式成立

式（38）中，n為滑動(dòng)窗口寬度，表示歷元新息值的采集個(gè)數(shù)。

當(dāng)量測(cè)異常時(shí)，量測(cè)噪聲與濾波遞推模型不符，需要引入一個(gè)由比例因子組成的自適應(yīng)矩陣Sk到算法中來調(diào)整測(cè)量協(xié)方差矩陣和讓理論值接近真實(shí)值，關(guān)系式如下

其中，式（39）左邊為實(shí)際的新息協(xié)方差，右邊為經(jīng)調(diào)整后的理論新息協(xié)方差，則Sk可推導(dǎo)如下

由以上公式可知，當(dāng)測(cè)量噪聲異常時(shí)，Sk是一個(gè)單位矩陣，對(duì)濾波無影響。當(dāng)量測(cè)出現(xiàn)異常時(shí)，由比例因子組成的自適應(yīng)矩陣對(duì)應(yīng)的項(xiàng)將變得相對(duì)更大，會(huì)消除測(cè)量值異常的影響。

矩陣Sk在計(jì)算過程中，由于誤差等因素的影響，其對(duì)角線元素可能小于1或者為負(fù)數(shù)。然而，自適應(yīng)矩陣Sk不能有元素小于1，故需要對(duì)自適應(yīng)矩陣做進(jìn)一步處理

（Sk）jj為Sk的第j個(gè)主對(duì)角線元素，當(dāng)量測(cè)噪聲突變時(shí)，會(huì)對(duì)濾波增益產(chǎn)生影響，式（34）變?yōu)?/p>

4 仿真驗(yàn)證

觀測(cè)雙星編隊(duì)平臺(tái)處于高軌道，非合作目標(biāo)衛(wèi)星處于低軌道，編隊(duì)衛(wèi)星和目標(biāo)衛(wèi)星的初始軌道參數(shù)如表1所示。其中，a為軌道半常軸，e為偏心率，i為軌道傾角，Ω為升交點(diǎn)赤經(jīng)，ω為近地點(diǎn)幅角，f為真近點(diǎn)角。

表1 編隊(duì)衛(wèi)星和空間非合作目標(biāo)初始軌道參數(shù)Table 1 Initial orbit parameters of satellite formation and NCST

利用STK11衛(wèi)星仿真軟件高精度軌道預(yù)報(bào)模型（High-precision Orbit Propagator，HPOP）生成NCST的星歷作為標(biāo)稱數(shù)據(jù)，生成觀測(cè)雙星的星歷作為雙星的狀態(tài)量，利用Access功能產(chǎn)生仿真的數(shù)據(jù)弧度，可得到一天之內(nèi)非合作目標(biāo)對(duì)觀測(cè)編隊(duì)衛(wèi)星的可見時(shí)間表。本文選定仿真時(shí)間為2019年3月11日07：56：50～09：03：23約4000s，采樣間隔為1s。

仿真的初始狀態(tài)矢量X0＝［-6535.038km-2142.106km 2329.429km-0.730km/s-4.327km/s -5.905km/s］，初始狀態(tài)誤差P0＝diag（25km，25km，25km，10m/s，10m/s，10m/s），方位角和高低角測(cè)量誤差均為0.0001rad，其他的仿真參數(shù)在上文已給出，仿真結(jié)果如圖3所示。

圖3 空間非合作目標(biāo)狀態(tài)估計(jì)誤差Fig.3 State estimation error of NCST

表2給出了兩種跟蹤方法的狀態(tài)估計(jì)均方根誤差（Root Mean Square Error，RMSE）結(jié)果對(duì)比，Monte Carlo仿真次數(shù)為100次。從結(jié)果中可以看出，在相同的仿真條件下，相比單星跟蹤，雙星跟蹤的收斂位置誤差減少了27.06%，收斂速度誤差減少了26.96%。

表2 兩種跟蹤方法均方根誤差對(duì)比Table 2 RMSE comparison of two tracking methods

4.1 量測(cè)故障

為了驗(yàn)證所設(shè)計(jì)的自適應(yīng)濾波算法的性能，跟蹤方式為雙星跟蹤，分別給出如下兩個(gè)算例。

算例1：在保持以上仿真條件不變的基礎(chǔ)上，模擬測(cè)量數(shù)據(jù)丟失的場(chǎng)景。在歷元2000s～2050s，方位角和高低角輸出都設(shè)置為0，ASRCKF和SRCKF的仿真結(jié)果如圖4、圖5所示。

圖4 算例1下的ASRCKF狀態(tài)估計(jì)誤差Fig.4 State estimation errors of ASRCKF in case1

圖5 算例1下的SRCKF狀態(tài)估計(jì)誤差Fig.5 State estimation errors of SRCKF in case1

算例2：模擬測(cè)量數(shù)據(jù)不準(zhǔn)確時(shí)的場(chǎng)景。在歷元2500s～2510s 高低角和方位角的角度測(cè)量基礎(chǔ)上加入0.0004rad的常值粗差，ASRCKF和SRCKF的仿真結(jié)果如圖6、圖7 所示。

圖6 算例2下的ASRCKF狀態(tài)估計(jì)誤差Fig.6 State estimation errors of ASRCKF in case2

圖7 算例2下SRCKF狀態(tài)估計(jì)誤差Fig.7 State estimation errors of SRCKF in case2

在此基礎(chǔ)上求解兩種濾波算法的均方根誤差，如表3 所示。

表3 兩種濾波算法均方根誤差對(duì)比Table 3 RMSE comparison of two filter algorithms

對(duì)比表2和表3中算例1的數(shù)據(jù)可以看到，當(dāng)出現(xiàn)量測(cè)缺失擾動(dòng)時(shí)，自適應(yīng)濾波結(jié)果幾乎無變化，依然能夠不受量測(cè)異常擾動(dòng)的影響。濾波結(jié)果趨于平穩(wěn)，而常規(guī)濾波則直接發(fā)散，濾波結(jié)果不正確。

對(duì)比表2和表3中算例2的數(shù)據(jù)可以看到，ASRCKF面對(duì)量測(cè)值不準(zhǔn)確擾動(dòng)時(shí)，濾波結(jié)果出現(xiàn)小幅度范圍的波動(dòng)，但是和正常量測(cè)相比，濾波結(jié)果變化不大且很快重新收斂；而SRCKF則出現(xiàn)了較大的波動(dòng)，濾波精度較差，雖然最終依然能夠收斂，但是需要較長(zhǎng)的時(shí)間。

5 結(jié)論

本文提出了一種雙星編隊(duì)對(duì)非合作目標(biāo)的僅測(cè)角跟蹤方法，建立了跟蹤的狀態(tài)模型和量測(cè)模型，利用ASRCKF濾波算法進(jìn)行狀態(tài)估計(jì)，并且完成了不同量測(cè)條件下的仿真估計(jì)。仿真結(jié)果表明，雙星編隊(duì)的聯(lián)合跟蹤相比單星具有更高的跟蹤精度，位置誤差和速度誤差分別減少了27.06%和26.96%。針對(duì)量測(cè)異常的情況，ASRCKF相比SRCKF能夠很好地抑制量測(cè)異常擾動(dòng)對(duì)濾波結(jié)果的影響。本文提出的方法在實(shí)際應(yīng)用中還存在一定的局限性，例如雙星編隊(duì)量測(cè)數(shù)據(jù)不同步對(duì)跟蹤精度的影響、測(cè)量角度象限的跳變對(duì)跟蹤精度的影響，這些問題都需要后續(xù)進(jìn)一步研究解決。