一種改進的GM-C-CPHD 空間多目標跟蹤算法

謝貝旭，張 艷，陳金濤，張任莉

（中山大學 航空航天學院，廣東 深圳 518060）

0 引言

隨著航天技術不斷進步，太空探索活動逐漸增加，新發射的航天器及衛星增多，使得近地空間逐漸擁擠［1］；同時，空間的目標解體頻發，產生的碎片數量逐漸增多，對衛星和航天器的運行造成巨大的威脅［2-3］。針對空間碎片，主要有以下4 個應對措施：監測和預警、碰撞規避和防護、離軌和棄置策略，以及主動清除［4-5］。從應對措施看，對軌道上目標進行精準識別［6］、定位與跟蹤均為應對的必備前提，以此提前掌握可能發生碰撞的信息，并對應地操控衛星躲避或提前精確地清除碎片。

多目標跟蹤（Multi-Target Tracking，MTT）［7-10］是指從含噪聲的傳感器量測中（存在虛警或雜波），估計監視區域內的目標數量和狀態。MTT 算法主要分為數據關聯和非數據關聯2 類，如圖卷積網絡（Graph Convolutional Network，GCN）、聯合概率數據關聯（Joint Probabilistic Data Association，JPDA）和多假設跟蹤（Multi-hypothesis Tracking，MHT）是基于數據關聯的MTT 算法［11-12］，需先完成量測的分配步驟，再進行目標狀態估計［13］。非數據關聯的方法，是在隨機有限集（Random Finite Set，RFS）的框架下進行，如概率假設密度（Probability Hypothesis Density，PHD）和勢概率假設密度（Cardinalized Probability Hypothesis Density，CPHD）濾波算法，可以避開關聯運算模塊，并直接進行多目標狀態估計，顯著緩解了運算的復雜程度［14-19］。

R.MAHLER 從貝葉斯濾波和有限集理論出發，提出PHD 濾波算法［20］，該算法遞歸多目標后驗Ⅰ階近似矩，而不是整體地后驗概率密度。隨后又提出CPHD濾波算法，聯合傳遞多目標概率密度的Ⅰ階近似矩及目標的勢分布，進一步提高跟蹤精度。但實現PHD、CPHD 算法存在積分的困難；VO 等［21］在此基礎上實現了高斯混合形式的PHD 及CPHD 濾波算法，降低了運算難度，并在二維仿真場景中應用。空間多目標跟蹤問題的突出難點之一是由攝動力引起的運動方程不確定性，與其相關的不確定性參數會導致狀態估計難度加大，進而導致環境對目標的影響不確定，參數識別效果變差，跟蹤性能變差［22］。為了提高對空間目標的跟蹤性能，MCCABE 和DEMARS 等［23-24］推導了一種考慮不確定性參數的高斯混合勢概率假設密度（Gaussian-mixture Cardinalized Probability Hypothesis Density，GM-C-PHD）濾波器；YANG 等［25］在無跡卡爾曼濾波框架下推導了考慮面質比（Area-to-Mass Ratio，AMR）不確定性參數（Gaussian-mixture Considering Probability Hypothesis Density，GM-CPHD）的濾波器，在空間跟蹤環境中使用。CPHD 濾波器同時傳遞多目標概率密度的Ⅰ階近似矩以及目標的勢分布，因此跟蹤精度上，尤其是對目標數目的估計，其誤差小于PHD。

針對空間目標跟蹤的模型不確定性問題，本文提出一種考慮面質比（Area-to-Mass Ratio，AMR）不確定性參數（Gaussian-mixture Considering Cardinalized Probability Hypothesis Density，GM-CCPHD）的算法。該算法在軌道動力學模型中考慮了不確定參數AMR，在無跡卡爾曼濾波（Unscented Kalman Filter，UKF）濾波框架下通過協方差傳遞的方式傳遞參數AMR 對位置、速度狀態估計的影響，以此減弱模型的不確定性影響程度，提高目標跟蹤水平。仿真分析表明，相較于原來的GMCPHD 濾波器，目標的跟蹤性能有所改善。

1 改進的GM-C-CPHD 濾波器

CPHD 濾波器在計算上存在積分困難，可通過高斯混合形式將遞歸過程中的目標強度函數及勢分布函數進行高斯疊加和處理，進而有效避免積分困難的問題。本節首先簡單分析不確定性參數AMR 對空間目標攝動力模型的影響，隨后提出一種考慮不確定性參數AMR 的GM-C-CPHD 濾波器，并闡述該算法的具體計算流程。

1.1 AMR 與攝動力的影響關系

在地心慣性坐標系中，空間目標的攝動力模型如下：

式中：r為空間目標的位置矢量，m；為位置矢量加速度，m/s2；αarp為太陽光壓加速度，m/s2；αdrag為大氣阻力加速度，m/s2；αΘ為太陽引力加速度，m/s2；αΠ為月球引力加速度，m/s2；αJ2為J2攝動加速度，m/s2。

與AMR 參數相關的有大氣阻力攝動和太陽光壓攝動，本文主要考慮AMR 參數對太陽光壓攝動力的影響。其影響關系表達式如下：

式中：S/m為面質比；rΘ為太陽的位置矢量，m；Cr為無量綱反射系數；Pr為地球表面處的光壓強度；re，Θ為地日心平均距離。

1.2 狀態增廣的系統方程

本文采用高斯混合的實現形式，并在UKF 濾波框架下傳遞概率密度函數及勢分布函數，以適應非線性的系統方程。目標的狀態方程和量測方程如下：

式中：f(·)和g(·)為非線性函數；wk和vk均為高斯白噪聲；c為面質比參數，即c=S/m；xk，c為考慮AMR 參數的增廣狀態矢量；zk為觀測矢量。

xk，c和zk如下：

式中：A，E分別為方位角及俯仰角，rad；，分別為兩者的角速度，rad/s；x、y、z分別為以觀測站為原點的三維空間坐標，m；，，分別為對應的速度，m/s。

1.3 高斯混合實現形式

1）初始化

在非線性跟蹤系統中，為了降低計算難度，利用高斯近似方法，將后驗強度中的非高斯分量拆解成若干高斯分量形式。這里將滿足均值為m、協方差為P的高斯概率密度記為N(·；m，P)。假設k時刻新生目標的強度函數是高斯混合形式：

2）預測步

記k-1 時刻真實目標的后驗強度函數為

則k時刻的預測強度函數可計算如下：

式中：vS，k|k-1(xc)為存活目標的強度。

存活目標的強度計算如下：

式中：pS，k為新生目標存活概率。

記后驗勢分布為pk-1(·)，則k時刻的預測勢分布函數如下：

式中：pΓ，k(·)為k時刻新生目標的勢分布函數；pS，k為新生目標存活概率；為組合系數；Jk-1和Jk|k-1為高斯混合分量的個數；ω為分量權重。

3）更新步

k時刻目標的后驗勢分布及后驗強度也是高斯混合形式，且兩者的更新如下：

式中：vD(·)為檢測概率下的更新強度函數。

1.4 濾波器實現

在使用高斯混合形式描述后驗強度和后驗勢分布后，使用UKF 濾波公式傳遞各個高斯分量，更新目標強度函數和勢分布，并實現GM-C-CPHD 濾波器。具體的濾波實現過程如下。

1）基于k-1 時刻的狀態估計量開始Sigma 點采樣，獲取1 組Sigma 點及其相對應的權值wl：

2）將Sigma 點代進運動方程，獲取Sigma 點的預測情況：

3）根據權值wl及各Sigma 點的預測值獲取狀態變量及協方差矩陣的預測情況如下：

4）根據上一步計算的狀態變量預測值再次進行Sigma 點采樣獲取另一組Sigma 點及相對應的權值wl：

5）將上一步獲取的Sigma 點代入觀測方程，求得各Sigma 點的觀測預測值為

6）根據權值wl及各Sigma 的觀測預測值獲取系統觀測及協方差的預測值為

7）計算增益矩陣Kk為

8）計算狀態變量及協方差矩陣的更新值：

2 仿真分析

為了比較GM-C-CPHD 濾波器與GM-CPHD濾波器在不同空間目標跟蹤環境下的跟蹤效果，共設置3 組仿真實驗，其變量控制見表1，3 個目標初始軌道根數見表2。其中目標1 存在時間為第1至第30 時刻；目標2 存在時間為第1 至第21 時刻；目標3 存在時間為第5 至第30 時刻。目標的濾波初始狀態標準差設置為：x、y、z方向的位置標準差為σx=σy=σz=1 km；x、y、z方向的速度標準差為σvx=σvy=σvz=1 m/s。地心地固坐標系（Earthcentered，Earth-fixed，CEF）測站參數設置見表3。使用最優子模式分配（Optimal Subpattern Assignment，OSPA）一致性度量進行性能評估，其中懲罰參數c=50 km，階次參數p=2。

表1 仿真實驗參數設置Tab.1 Parameter settings of the simulation tests

表2 地球同步軌道目標的初始軌道根數Tab.2 Initial Keplerian elements of the geosynchronous orbit objects

表3 地心地固坐標系下測站參數說明Tab.3 Description of the ground station parameters in the coordinate system

2.1 仿真實驗1

GM-C-CPHD 和GM-CPHD 濾波的位置、速度OSPA 誤差及目標數目估計結果如圖1 所示。其中，第30 時刻的位置、速度OSPA 誤差統計值及目標數目錯誤估計的時刻數見表4。可以看到，隨著時間的迭代，2 個濾波器最后均能收斂，但GM-CCPHD 濾波器的位置及速度OSPA 誤差更小。GM-C-CPHD 濾波器的位置和速度OSPA 誤差的最后收斂數值分別為130.311 和0.019，而GMCPHD 濾波器則為3951.62 和0.74。在目標數目估計方面，圖1 顯示2 個濾波器在大部分時間均可以準確估計，其中GM-C-CPHD 在第15 時刻目標數目估計錯誤，估計錯誤時刻數為1；而GM-CPHD濾波器在第15 和18 時刻均估計錯誤，錯誤估計時刻數為2。2 個濾波器在第30 時刻不同分位值的OSPA 誤差見表5。由表5 可知，GM-C-CPHD 濾波器在第30 時刻的估計精度更高，顯示出所提方法的有效性。

圖1 OSPA 誤差和目標數目估計Fig.1 Estimation of the OSPA errors and target numbers

表4 第30 時刻位置、速度OSPA 誤差和目標數目錯誤估計時刻數Tab.4 OSPA errors of the position and velocity at the 30th moment and the total moment number of false estimation regarding the target number

表5 第30 時刻不同分位值的OSPA 誤差Tab.5 OSPA errors regarding different quantile values at the 30th moment

2.2 仿真實驗2

仿真實驗2 進一步降低了檢測概率，設置為0.95，其他實驗參數設置與仿真實驗1 一致。GM-C-CPHD 和GM-CPHD 濾波的仿真結果如圖2、表6 和表7 所示。由圖2 可知，2 個濾波器的跟蹤誤差值最后同樣收斂，GM-C-CPHD 濾波器的位置及速度OSPA 誤差收斂于更小的數值。由表6 可知，本文方法收斂數值分別為189.751 和0.025 3，而GM-CPHD 濾波器則為7 548.05 和1.463 8。在目標的數目估計方面，由于檢測概率減小，2 個濾波器的估計效果相比于仿真1 有所降低。其中GM-CCPHD 濾波器在第11、15 和第22 時刻估計不準確；而GM-CPHD 濾波器在第11、12、15 和第22 時刻估計不準確。此外，表7 表明在5%分位值情況下，2 個濾波器的OSPA 誤差較為接近，但其他分位值情況下，GM-C-CPHD 濾波器的位置及速度估計誤差小于GM-CPHD 濾波器。

圖2 GM-CPHD 和GM-C-CPHD 的OSPA 誤差和目標數目估計Fig.2 OSPA errors and estimated target number

表6 第30 時刻位置、速度OSPA 誤差和目標數目錯誤估計時刻數Tab.6 OSPA errors of the position and velocity at the 30th moment and the total moment number of false estimation regarding the target number

表7 第30 時刻不同分位值的OSPA 誤差Tab.7 OSPA errors regarding different quantile values at the 30th moment

2.3 仿真實驗3

在仿真2 的基礎上，進一步考慮雜波的情況，雜波的產生服從泊松分布，泊松分布參數λ=1，其他參數設置與仿真實驗2 保持一致。GM-C-CPHD 和GMCPHD 濾波的仿真結果如圖3、表8 和表9 所示。

圖3 GM-CPHD 和GM-C-CPHD 的OSPA 誤差和目標數目估計Fig.3 OSPA errors and the estimated target number by GM-CPHD and GM-C-CPHD

表8 第30 時刻位置、速度OSPA 誤差和目標數目錯誤估計時刻數Tab.8 OSPA errors of the position and velocity at the 30th moment and the total moment number of false estimation regarding the target number

表9 第30 時刻不同分位值的OSPA 誤差Tab.9 OSPA errors regarding different quantile values at the 30th moment

如圖3 所示，GM-C-CPHD 濾波器的位置及速度OSPA 誤差更小，跟蹤精度更優。由表8 可知，本文方法位置及速度OSPA 誤差收斂數值分別為154.473 0 和0.023 2，而GM-CPHD 濾波器則為13 843.3 和2.737 4。在目標數目估計方面，由于跟蹤場景存在雜波，估計效果相比于仿真2 均有所減小。其中GM-CPHD 濾波器在第11、14、15、19、20和21 時刻有估計錯誤情況，而GM-C-CPHD 濾波器只在第14 和15 時刻估計存在偏差。此外，表9顯示，2 個濾波器在第30 時刻不同分位值的OSPA誤差，可知GM-C-CPHD 濾波器在各分位值情況下的位置及速度估計誤差均小于GM-CPHD 濾波器。

2.4 其他指標比較及分析

除了使用OSPA 指標評價濾波器的性能外，還統計了2 個濾波器在所有時刻的OSPA 誤差均值、標準差及均方根誤差，統計結果見表10。由表10 可知，GM-C-CPHD 濾波器的OSPA 誤差標準差比GM-CPHD 濾波器的標準差大，說明后期誤差減小幅度更大，而平均值及均方根誤差均比GM-CPHD濾波器小，表明了GM-C-CPHD 濾波器對目標狀態估計的整體誤差水平更低。

表10 OSPA 誤差統計值Tab.10 Statistical OSPA errors for all the tests

綜合仿真結果及其他指標進行分析，由于不確定性參數AMR 對軌道動力學模型存在影響，使得濾波過程的協方差矩陣不能清晰地反映多目標跟蹤的精度。其中，GM-CPHD 濾波器在遞歸過程中將AMR 參數假設為常數，參數的影響效果就會確定化，那么濾波器估計參數的不確定性也會受到影響，進而在協方差中表現出來。而GM-C-CPHD 濾波器通過協方差增廣的方式將AMR 參數考慮到協方差傳遞過程中，提高該不確定性參數的可觀性，進而降低其對濾波估計的影響，使得目標跟蹤精度得到改善。

3 結束語

本文針對空間目標運動模型的不確定性問題，提出一種考慮面質比不確定性參數的GM-CCPHD 多目標跟蹤算法，該算法在軌道動力學模型中考慮了面質比不確定參數，并在UKF 濾波框架下，通過協方差傳遞的方式傳遞該參數對位置、速度狀態估計的影響，以此緩解空間目標運動模型的不確定性影響。通過考慮降低檢測概率、增加雜波率等條件設置3 組仿真實驗，仿真結果表明，與GMCPHD 濾波器相比，GM-C-CPHD 濾波器在OSPA誤差及其均值、均方根誤差等指標方面均有所改善。未來，主要從以下2 個方面開展工作：1）跟蹤更多數量且密集的空間多目標；2）進一步優化算法，考慮檢測概率的不確定性影響。