雙星編隊對空間非合作目標的僅測角跟蹤方法

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

（1.南京理工大學自動化學院， 南京 210094；2.南京理工大學理學院， 南京 210094；3.中國人民解放軍 66483 部隊， 北京 100093）

0 引言

隨著全球太空資源開發的提高和未來太空作戰趨勢的加劇，越來越多的國家擁有了進入太空和利用太空的能力，太空變得更加具有競爭性和對抗性。其中，空間監視系統可以監視空間目標的運行，確定潛在的威脅，還可以預測空間目標的軌道，其水平的高低直接制約著空間對抗能力的發揮［1］。同時，空間監視系統相比地面監視系統具有全天候、視場廣、不受大氣環境影響等特點，發展天基監視系統就顯得尤為重要。

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

針對上述問題，結合現有分布式衛星快速發展的背景，利用編隊衛星對目標進行多視線測量，改善觀測的幾何條件，融合觀測數據，在保證實時性的同時，提高了對非合作目標的跟蹤精度和穩定性。此外，受空間幾何因素和光學因素的影響［7］，利用光電傳感器對目標進行觀測時會出現某段時間內測量值丟失和測量值誤差過大的現象，這些測量故障會造成跟蹤精度不高甚至跟蹤結果發散。基于此設計自適應濾波，該濾波基于新息構建自適應因子，設計增益調節矩陣，提高了濾波的魯棒性。最后，通過仿真對提出方法的有效性進行了驗證和分析。

1 空間目標狀態模型

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

圖1 衛星編隊非合作目標跟蹤系統Fig.1 Satellite formation for NCST tracking system

在僅考慮地球J2攝動時，衛星的軌道方程如下

式（1）中，f為地球非球形J2攝動下所產生的加速度。f在地心慣性坐標系（Earth Centered Inertial，ECI）下可表示為［6］

則式（2）可改寫為

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

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

式（6）離散化后可進一步寫為

利用Euler法對式（7）進行求解

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

2 測量模型

2.1 觀測方程中的坐標轉換

在典型的僅測角無源定位中，方位角和高低角的計算定義為目標與觀測器位置矢量之差。以往的文獻多是直接在ECI下建立觀測方程，公式簡單，但是觀測數據是在主星和從星各自的軌道坐標系下測量的，而雙星觀測器和目標的位置矢量都建立在ECI下，而且需要建立編隊雙星之間的坐標關系。因此，需要進行坐標轉換才更接近真實系統，Euler角的坐標轉換如圖2所示。

圖2 Euler角坐標轉換示意圖Fig.2 Schematic diagram of Euler angles coordinate transformation

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

2.2 目標運動觀測方程

設觀測主星對目標測得的高低角和方位角分別為αk1、βk1，觀測從星對目標測得的高低角和方位角分別為αk2、βk2，其定義如下

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

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

從星的測量值通過星間鏈路傳輸給主星，zk2為目標在從星軌道系下的量測值，應轉換到主星軌道系下。將zk2從從星軌道系轉換到ECI，其坐標轉換矩陣為；再從ECI轉換到主星軌道坐標系，其坐標轉換矩陣為G1。因此，為從星軌道坐標系到主星軌道坐標系的轉換矩陣。最終，雙星編隊對非合作目標的量測值為

觀測雙星對目標的觀測量均為單位矢量，故量測方程的右端也應為單位矢量。在慣性坐標系下，觀測雙星到NCST的單位觀測矢量為

式（15）建立在主星軌道坐標系下，式（16）建立在ECI下。為了方便濾波時協方差矩陣等的獲取，需進行坐標轉換，最終可寫為

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

3 跟蹤濾波算法

3.1 平方根容積Kalman濾波

天基僅測角跟蹤是非線性狀態估計，平方根容積Kalman濾波（Square Root Cubature Kalman Filter，SRCKF）使用基于容積規則的數值積分方法直接計算非線性函數的均值和方差，只需計算函數值，避免求導計算。同時，SRCKF利用協方差的平方根直接參與遞推運算，確保了協方差矩陣的對稱性和正定性，能保證濾波的精度和穩定性。

SRCKF濾波的步驟如下：

（1）初始化

（2）時間更新

計算容積點，有

式（21）中，N＝2m，m為狀態的維數；S0＝chol（P0），即。記m維單位向量為e＝［1 0 … 0］T，使用符號［1］表示對e的元素進行全排列和改變元素符號所產生的點集，稱為完整全對稱點集。

（3）量測更新

3.2 自適應SRCKF

在復雜的太空環境中，由于地球遮擋等幾何因素、目標輻射強度以及目標背景較亮等影響光電傳感器工作的光學因素，難免會出現量測信息不準確或者信息傳輸的暫時中斷，導致不能提供正確的量測信息。此時，傳統的Kalman濾波精度會大大降低，甚至導致濾波發散。因此，需要設計自適應濾波來提高其魯棒性。自適應濾波一般是基于新息和基于殘差設計的［9-10］，本文采用的是基于新息的方法來設計自適應平方根容積Kalman濾波（Adaptive Square Root Cubature Kalman Filter，ASRCKF）。

自適應濾波通過采集過去一定時間的新息，利用采集到的新息實時調整濾波器的誤差矩陣［11-12］。具體調整策略是通過比較實時估計的矩陣值和設定的矩陣值，當相差大于一定的閾值時就進行調整。

真實測量值和一步預測值計算的殘差公式如下

式（37）中，Zk為k時刻的真實測量值，為k時刻的量測估計值。

若系統量測的真實誤差統計特性與濾波遞推誤差模型一致時，以下等式成立

式（38）中，n為滑動窗口寬度，表示歷元新息值的采集個數。

當量測異常時，量測噪聲與濾波遞推模型不符，需要引入一個由比例因子組成的自適應矩陣Sk到算法中來調整測量協方差矩陣和讓理論值接近真實值，關系式如下

其中，式（39）左邊為實際的新息協方差，右邊為經調整后的理論新息協方差，則Sk可推導如下

由以上公式可知，當測量噪聲異常時，Sk是一個單位矩陣，對濾波無影響。當量測出現異常時，由比例因子組成的自適應矩陣對應的項將變得相對更大，會消除測量值異常的影響。

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

（Sk）jj為Sk的第j個主對角線元素，當量測噪聲突變時，會對濾波增益產生影響，式（34）變為

4 仿真驗證

觀測雙星編隊平臺處于高軌道，非合作目標衛星處于低軌道，編隊衛星和目標衛星的初始軌道參數如表1所示。其中，a為軌道半常軸，e為偏心率，i為軌道傾角，Ω為升交點赤經，ω為近地點幅角，f為真近點角。

表1 編隊衛星和空間非合作目標初始軌道參數Table 1 Initial orbit parameters of satellite formation and NCST

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

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

圖3 空間非合作目標狀態估計誤差Fig.3 State estimation error of NCST

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

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

4.1 量測故障

為了驗證所設計的自適應濾波算法的性能，跟蹤方式為雙星跟蹤，分別給出如下兩個算例。

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

圖4 算例1下的ASRCKF狀態估計誤差Fig.4 State estimation errors of ASRCKF in case1

圖5 算例1下的SRCKF狀態估計誤差Fig.5 State estimation errors of SRCKF in case1

算例2：模擬測量數據不準確時的場景。在歷元2500s～2510s 高低角和方位角的角度測量基礎上加入0.0004rad的常值粗差，ASRCKF和SRCKF的仿真結果如圖6、圖7 所示。

圖6 算例2下的ASRCKF狀態估計誤差Fig.6 State estimation errors of ASRCKF in case2

圖7 算例2下SRCKF狀態估計誤差Fig.7 State estimation errors of SRCKF in case2

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

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

對比表2和表3中算例1的數據可以看到，當出現量測缺失擾動時，自適應濾波結果幾乎無變化，依然能夠不受量測異常擾動的影響。濾波結果趨于平穩，而常規濾波則直接發散，濾波結果不正確。

對比表2和表3中算例2的數據可以看到，ASRCKF面對量測值不準確擾動時，濾波結果出現小幅度范圍的波動，但是和正常量測相比，濾波結果變化不大且很快重新收斂；而SRCKF則出現了較大的波動，濾波精度較差，雖然最終依然能夠收斂，但是需要較長的時間。

5 結論

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