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

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

0 引言

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

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

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

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

1 改進(jìn)的GM-C-CPHD 濾波器

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

1.1 AMR 與攝動(dòng)力的影響關(guān)系

在地心慣性坐標(biāo)系中，空間目標(biāo)的攝動(dòng)力模型如下：

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

與AMR 參數(shù)相關(guān)的有大氣阻力攝動(dòng)和太陽(yáng)光壓攝動(dòng)，本文主要考慮AMR 參數(shù)對(duì)太陽(yáng)光壓攝動(dòng)力的影響。其影響關(guān)系表達(dá)式如下：

式中：S/m為面質(zhì)比；rΘ為太陽(yáng)的位置矢量，m；Cr為無(wú)量綱反射系數(shù)；Pr為地球表面處的光壓強(qiáng)度；re，Θ為地日心平均距離。

1.2 狀態(tài)增廣的系統(tǒng)方程

本文采用高斯混合的實(shí)現(xiàn)形式，并在UKF 濾波框架下傳遞概率密度函數(shù)及勢(shì)分布函數(shù)，以適應(yīng)非線(xiàn)性的系統(tǒng)方程。目標(biāo)的狀態(tài)方程和量測(cè)方程如下：

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

xk，c和zk如下：

式中：A，E分別為方位角及俯仰角，rad；，分別為兩者的角速度，rad/s；x、y、z分別為以觀測(cè)站為原點(diǎn)的三維空間坐標(biāo)，m；，，分別為對(duì)應(yīng)的速度，m/s。

1.3 高斯混合實(shí)現(xiàn)形式

1）初始化

在非線(xiàn)性跟蹤系統(tǒng)中，為了降低計(jì)算難度，利用高斯近似方法，將后驗(yàn)強(qiáng)度中的非高斯分量拆解成若干高斯分量形式。這里將滿(mǎn)足均值為m、協(xié)方差為P的高斯概率密度記為N(·；m，P)。假設(shè)k時(shí)刻新生目標(biāo)的強(qiáng)度函數(shù)是高斯混合形式：

2）預(yù)測(cè)步

記k-1 時(shí)刻真實(shí)目標(biāo)的后驗(yàn)強(qiáng)度函數(shù)為

則k時(shí)刻的預(yù)測(cè)強(qiáng)度函數(shù)可計(jì)算如下：

式中：vS，k|k-1(xc)為存活目標(biāo)的強(qiáng)度。

存活目標(biāo)的強(qiáng)度計(jì)算如下：

式中：pS，k為新生目標(biāo)存活概率。

記后驗(yàn)勢(shì)分布為pk-1(·)，則k時(shí)刻的預(yù)測(cè)勢(shì)分布函數(shù)如下：

式中：pΓ，k(·)為k時(shí)刻新生目標(biāo)的勢(shì)分布函數(shù)；pS，k為新生目標(biāo)存活概率；為組合系數(shù)；Jk-1和Jk|k-1為高斯混合分量的個(gè)數(shù)；ω為分量權(quán)重。

3）更新步

k時(shí)刻目標(biāo)的后驗(yàn)勢(shì)分布及后驗(yàn)強(qiáng)度也是高斯混合形式，且兩者的更新如下：

式中：vD(·)為檢測(cè)概率下的更新強(qiáng)度函數(shù)。

1.4 濾波器實(shí)現(xiàn)

在使用高斯混合形式描述后驗(yàn)強(qiáng)度和后驗(yàn)勢(shì)分布后，使用UKF 濾波公式傳遞各個(gè)高斯分量，更新目標(biāo)強(qiáng)度函數(shù)和勢(shì)分布，并實(shí)現(xiàn)GM-C-CPHD 濾波器。具體的濾波實(shí)現(xiàn)過(guò)程如下。

1）基于k-1 時(shí)刻的狀態(tài)估計(jì)量開(kāi)始Sigma 點(diǎn)采樣，獲取1 組Sigma 點(diǎn)及其相對(duì)應(yīng)的權(quán)值wl：

2）將Sigma 點(diǎn)代進(jìn)運(yùn)動(dòng)方程，獲取Sigma 點(diǎn)的預(yù)測(cè)情況：

3）根據(jù)權(quán)值wl及各Sigma 點(diǎn)的預(yù)測(cè)值獲取狀態(tài)變量及協(xié)方差矩陣的預(yù)測(cè)情況如下：

4）根據(jù)上一步計(jì)算的狀態(tài)變量預(yù)測(cè)值再次進(jìn)行Sigma 點(diǎn)采樣獲取另一組Sigma 點(diǎn)及相對(duì)應(yīng)的權(quán)值wl：

5）將上一步獲取的Sigma 點(diǎn)代入觀測(cè)方程，求得各Sigma 點(diǎn)的觀測(cè)預(yù)測(cè)值為

6）根據(jù)權(quán)值wl及各Sigma 的觀測(cè)預(yù)測(cè)值獲取系統(tǒng)觀測(cè)及協(xié)方差的預(yù)測(cè)值為

7）計(jì)算增益矩陣Kk為

8）計(jì)算狀態(tài)變量及協(xié)方差矩陣的更新值：

2 仿真分析

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

表1 仿真實(shí)驗(yàn)參數(shù)設(shè)置Tab.1 Parameter settings of the simulation tests

表2 地球同步軌道目標(biāo)的初始軌道根數(shù)Tab.2 Initial Keplerian elements of the geosynchronous orbit objects

表3 地心地固坐標(biāo)系下測(cè)站參數(shù)說(shuō)明Tab.3 Description of the ground station parameters in the coordinate system

2.1 仿真實(shí)驗(yàn)1

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

圖1 OSPA 誤差和目標(biāo)數(shù)目估計(jì)Fig.1 Estimation of the OSPA errors and target numbers

表4 第30 時(shí)刻位置、速度OSPA 誤差和目標(biāo)數(shù)目錯(cuò)誤估計(jì)時(shí)刻數(shù)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 時(shí)刻不同分位值的OSPA 誤差Tab.5 OSPA errors regarding different quantile values at the 30th moment

2.2 仿真實(shí)驗(yàn)2

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

圖2 GM-CPHD 和GM-C-CPHD 的OSPA 誤差和目標(biāo)數(shù)目估計(jì)Fig.2 OSPA errors and estimated target number

表6 第30 時(shí)刻位置、速度OSPA 誤差和目標(biāo)數(shù)目錯(cuò)誤估計(jì)時(shí)刻數(shù)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 時(shí)刻不同分位值的OSPA 誤差Tab.7 OSPA errors regarding different quantile values at the 30th moment

2.3 仿真實(shí)驗(yàn)3

在仿真2 的基礎(chǔ)上，進(jìn)一步考慮雜波的情況，雜波的產(chǎn)生服從泊松分布，泊松分布參數(shù)λ=1，其他參數(shù)設(shè)置與仿真實(shí)驗(yàn)2 保持一致。GM-C-CPHD 和GMCPHD 濾波的仿真結(jié)果如圖3、表8 和表9 所示。

圖3 GM-CPHD 和GM-C-CPHD 的OSPA 誤差和目標(biāo)數(shù)目估計(jì)Fig.3 OSPA errors and the estimated target number by GM-CPHD and GM-C-CPHD

表8 第30 時(shí)刻位置、速度OSPA 誤差和目標(biāo)數(shù)目錯(cuò)誤估計(jì)時(shí)刻數(shù)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 時(shí)刻不同分位值的OSPA 誤差Tab.9 OSPA errors regarding different quantile values at the 30th moment

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

2.4 其他指標(biāo)比較及分析

除了使用OSPA 指標(biāo)評(píng)價(jià)濾波器的性能外，還統(tǒng)計(jì)了2 個(gè)濾波器在所有時(shí)刻的OSPA 誤差均值、標(biāo)準(zhǔn)差及均方根誤差，統(tǒng)計(jì)結(jié)果見(jiàn)表10。由表10 可知，GM-C-CPHD 濾波器的OSPA 誤差標(biāo)準(zhǔn)差比GM-CPHD 濾波器的標(biāo)準(zhǔn)差大，說(shuō)明后期誤差減小幅度更大，而平均值及均方根誤差均比GM-CPHD濾波器小，表明了GM-C-CPHD 濾波器對(duì)目標(biāo)狀態(tài)估計(jì)的整體誤差水平更低。

表10 OSPA 誤差統(tǒng)計(jì)值Tab.10 Statistical OSPA errors for all the tests

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

3 結(jié)束語(yǔ)

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