盧哲俊，胡衛(wèi)東

國防科學(xué)技術(shù)大學(xué) 電子科學(xué)與工程學(xué)院 ATR國家重點實驗室，長沙 410073

基于隨機有限集的空間碎片群運動狀態(tài)估計

盧哲俊，胡衛(wèi)東*

國防科學(xué)技術(shù)大學(xué) 電子科學(xué)與工程學(xué)院 ATR國家重點實驗室，長沙 410073

傳統(tǒng)的空間目標(biāo)監(jiān)測是建立在單目標(biāo)狀態(tài)估計基礎(chǔ)之上，在面對突發(fā)產(chǎn)生的大量空間碎片時，由于碎片尺寸小，且密集分布以“群”的方式出現(xiàn)，傳統(tǒng)單目標(biāo)處理方法很難奏效。以“群”整體作為處理對象，基于隨機有限集(RFS)技術(shù)，對“群”的狀態(tài)特征進行估計。為了解決漏檢目標(biāo)密度分配問題和軌跡關(guān)聯(lián)問題，提出一種面向量測的改進集勢概率假設(shè)密度(CPHD)濾波器，并結(jié)合濾波后的信息處理過程，完成了對低軌空間碎片群的目標(biāo)密度分布、群內(nèi)目標(biāo)數(shù)以及群內(nèi)顯著目標(biāo)的狀態(tài)估計。在仿真實驗中，提出的濾波器表現(xiàn)明顯優(yōu)于傳統(tǒng)濾波器和標(biāo)準(zhǔn)CPHD濾波器，且在某些傳統(tǒng)濾波器和標(biāo)準(zhǔn)CPHD濾波器已失效的情況下，所提技術(shù)仍能有效工作。

空間碎片；群目標(biāo)；狀態(tài)估計；隨機有限集；改進CPHD濾波器

隨著空間碎片數(shù)量的不斷增加，對空間活動安全造成嚴重威脅，對空間碎片的監(jiān)測就顯得尤為重要。但是空間碎片數(shù)量大、尺寸小導(dǎo)致監(jiān)測困難，尤其是由于解體或碰撞產(chǎn)生的高密度空間碎片集群成為空間監(jiān)測和編目的難點[1-2]。例如2007年1月風(fēng)云1C衛(wèi)星產(chǎn)生之后6個月已編目的碎片數(shù)量達到了1 967個[3]，2009年2月俄羅斯廢棄衛(wèi)星Cosmos2251和美國商業(yè)衛(wèi)星Iridium33碰撞10個月后產(chǎn)生的已編目空間碎片數(shù)量達到1 632個[4]?？臻g碎片在剛產(chǎn)生之后的短時間里相距很近，形成了高密度的空間碎片“云”，且這些碎片中的絕大部分尺寸較小，觀測困難。這時候，地基雷達系統(tǒng)由于檢測概率低，目標(biāo)密集無法分辨，以及目標(biāo)數(shù)量大且未知等因素導(dǎo)致無法對空間碎片進行有效監(jiān)測和編目。

面對這種集群目標(biāo)，傳統(tǒng)方式是先試圖估計每個目標(biāo)的狀態(tài)，然后才對目標(biāo)群的狀態(tài)有個認識。因為傳統(tǒng)的多目標(biāo)濾波算法是先對量測和目標(biāo)進行關(guān)聯(lián)，然后使用單目標(biāo)濾波器對每個目標(biāo)的狀態(tài)進行估計。最常用的兩種方法是多假設(shè)跟蹤(Multiple Hypothesis Tracking，MHT)濾波器[5]和聯(lián)合概率數(shù)據(jù)關(guān)聯(lián)(Joint Probabilistic Data Association，JPDA)濾波器[6]。它們在面對大量密集目標(biāo)且目標(biāo)數(shù)未知的情況時，濾波表現(xiàn)會大大下降，因為目標(biāo)關(guān)聯(lián)過程會變得復(fù)雜且不穩(wěn)定，尤其當(dāng)檢測概率較低時。而且這兩種傳統(tǒng)濾波器無法對區(qū)域中的目標(biāo)數(shù)進行估計，只有通過航跡起始和航跡終止過程，依據(jù)維持的航跡數(shù)來判斷目標(biāo)數(shù)，但是航跡起始和終止過程由于檢測概率低變得不穩(wěn)定。因此，在面對這種密集空間碎片群的時候，相比對各目標(biāo)單獨處理，一種有效的策略就是對目標(biāo)“群”進行整體處理[7]。先對整個空間碎片群的狀態(tài)進行估計，在信息積累足夠或是目標(biāo)逐漸分開之后再對單個目標(biāo)的運動狀態(tài)進行估計，這樣即避免了單個目標(biāo)處理的困難，又能夠有效維持對空間碎片的監(jiān)測。“群監(jiān)測”的過程中，群特征可以用來對群進行標(biāo)識和區(qū)分，群特征包括群幾何(群大小、結(jié)構(gòu)等)、群內(nèi)目標(biāo)數(shù)量等。如果無法對群中所有目標(biāo)狀態(tài)進行估計，可以選擇群中檢測概率較高目標(biāo)(本文稱之為顯著目標(biāo))的狀態(tài)進行估計，用顯著目標(biāo)運動狀態(tài)來描述群的運動狀態(tài)。因為群內(nèi)目標(biāo)運動狀態(tài)具有相似性，以群內(nèi)顯著目標(biāo)運動狀態(tài)代表群運動狀態(tài)是合理的。

面對這種需求，在傳統(tǒng)方法失效的情況下，隨機有限集(Random Finite Set，RFS)理論提供了一個解決集群目標(biāo)濾波的有效工具?；赗FS理論的多目標(biāo)貝葉斯濾波器(Multi-target Bayes Filter)是一種將數(shù)據(jù)關(guān)聯(lián)、目標(biāo)數(shù)估計和狀態(tài)估計統(tǒng)一在一個概率框架下的多目標(biāo)濾波算法[8]。在很多近期的研究成果中，多目標(biāo)貝葉斯濾波器已被應(yīng)用到了空間目標(biāo)濾波之中[9-10]。因為最佳多目標(biāo)貝葉斯濾波器的計算復(fù)雜度太高，且很多情況下無法實現(xiàn)，因此有很多基于特定隨機過程假設(shè)的近似多目標(biāo)貝葉斯濾波器被提出[11]。其中，基于矩近似的概率假設(shè)密度(Probability Hypothesis Density，PHD)濾波器[12]和集勢概率假設(shè)密度(Cardinalized PHD，CPHD)濾波器[13]具有簡潔的公式和低的計算復(fù)雜度，可以滿足對密集空間碎片實時處理需求。不同于MHT和JPDA濾波器，PHD和CPHD濾波器不是直接估計目標(biāo)狀態(tài)，而是對目標(biāo)密度分布進行估計，密度高的地方認為目標(biāo)出現(xiàn)概率大，密度低的地方則出現(xiàn)概率小。通過對目標(biāo)密度分布的估計可以獲得群結(jié)構(gòu)特征，即群內(nèi)目標(biāo)分布情況和群的范圍大小[14]。而目標(biāo)密度的積分即目標(biāo)數(shù)估計，因此通過對群區(qū)域內(nèi)的PHD進行積分可獲得群目標(biāo)數(shù)的估計。但是基于一階矩近似的PHD濾波器的目標(biāo)數(shù)估計十分不穩(wěn)定，且無法正確計算漏檢目標(biāo)的PHD[15]，因此傳遞部分二階矩的CPHD濾波器被提出。CPHD濾波器在濾波過程中傳遞目標(biāo)密度分布的同時還傳遞目標(biāo)數(shù)的估計函數(shù)，因此可以獲得穩(wěn)定的目標(biāo)數(shù)估計。但是CPHD濾波器仍存在一些問題，比如CPHD濾波器雖然可以正確計算漏檢目標(biāo)的PHD，但是無法將其分配到正確航跡上[16]，且CPHD濾波器不繼承目標(biāo)標(biāo)記，在提取群內(nèi)顯著目標(biāo)狀態(tài)上出現(xiàn)困難。已經(jīng)提出了一些方法使CPHD濾波器具有航跡關(guān)聯(lián)能力[17-18]。本文選擇CPHD濾波器對空間碎片群進行濾波，并提出一種新的改進CPHD濾波器，新濾波器中的PHD以量測進行聚合并以此完成數(shù)據(jù)關(guān)聯(lián)過程，不引入其他關(guān)聯(lián)算法，同時解決了漏檢目標(biāo)的PHD分配問題。完成數(shù)據(jù)關(guān)聯(lián)后，通過對提取出的空間目標(biāo)軌跡的檢測概率進行估計，取檢測概率最大的目標(biāo)作為群顯著目標(biāo)，從而完成對群顯著目標(biāo)的運動狀態(tài)估計。

在仿真實驗中，改進的CPHD濾波器在不同的觀測條件下與JPDA濾波器進行比較，展現(xiàn)了改進CPHD濾波器的在目標(biāo)狀態(tài)和目標(biāo)數(shù)估計上的優(yōu)良性能。在檢測概率很低的情況下，JPDA濾波器已經(jīng)失效，而改進CPHD濾波器仍能有效工作。本文提出的方法為空間碎片群的監(jiān)測和編目提供了一種新的有效方法。

本文結(jié)構(gòu)如下：第1節(jié)對低軌空間碎片群進行描述；第2節(jié)介紹了空間碎片群的隨機有限集描述；第3節(jié)提出改進的CPHD濾波器;第4節(jié)對濾波結(jié)果進行數(shù)據(jù)處理過程；第5節(jié)通過仿真實驗驗證方法的有效性；第6節(jié)給出本文結(jié)論。

1 空間碎片群描述

1.1 空間碎片運動模型

xk=Fk(xk-1)+Qk

(1)

式中：Fk(·)為非線性的動態(tài)模型；Qk為模型噪聲。

由于空間目標(biāo)在太空中的受力非常復(fù)雜，因此空間碎片的動態(tài)模型是非常復(fù)雜的非線性模型[19]。本文的動態(tài)模型考慮低軌空間碎片的受力包括：①n體運動,將碎片目標(biāo)、地球、太陽和月亮作為質(zhì)點處理，空間目標(biāo)與各質(zhì)點之間由于萬有引力構(gòu)成一個n體受力系統(tǒng)。② 地球非球形攝動,由于地球的密度分布不均勻，其形狀也不是標(biāo)準(zhǔn)球形，而是相當(dāng)不規(guī)則的，其形狀和密度分布的非球形部分產(chǎn)生了對空間目標(biāo)的攝動源。這里考慮10×10階的非球形攝動位函數(shù)。③ 大氣阻力,在大氣層中飛行的空間目標(biāo)受到大氣阻力的影響。大氣阻力主要影響低軌目標(biāo)，對于高軌目標(biāo)的大氣阻力可以忽略，影響大氣阻力的因素有大氣密度，碎片與大氣的相對運動速度，目標(biāo)的阻力系數(shù)和面質(zhì)比。④ 太陽光壓,太陽光線照在碎片目標(biāo)上會產(chǎn)生輻射壓，太陽光壓與碎片目標(biāo)材料的反射系數(shù)以及碎片目標(biāo)面向太陽的面質(zhì)比有關(guān)。

1.2 典型空間碎片群特征

本文關(guān)注的空間碎片群指的是低軌空間突發(fā)產(chǎn)生的大量密集空間碎片組成的集群。這里以風(fēng)云1C空間碎片群為例，對空間碎片群的特征進行描述。取6個密集分布的風(fēng)云1C碎片目標(biāo)組成的集群為研究目標(biāo)，研究數(shù)據(jù)采用北美空間防空司令部(NORAD)公布的雙行根數(shù)(Two-Line Element，TLE)，這6個目標(biāo)的衛(wèi)星NORAD編號為29727～29732。雙行根數(shù)是公開數(shù)據(jù)，這6個碎片最早公布的TLE數(shù)據(jù)時間是2007年1月18日，而碎片產(chǎn)生的時間是2007年1月11日。為了獲得碎片產(chǎn)生之后不久的數(shù)據(jù)，使用NORAD 專門針對TLE數(shù)據(jù)開發(fā)的軌道傳播模型——SGP4模型進行回溯[20]，回溯至2007年1月11日中午，產(chǎn)生一段碎片群運動弧段如圖1所示。此弧段具體運動時間從2007年1月11日12:05:46到12:08:38。這6個碎片的平均軌道根數(shù)數(shù)據(jù)如表1所示。使用軌道模型回溯的結(jié)果近似真實情況，通過這個近似結(jié)果對碎片群進行研究，考慮采用下述特征對空間碎片群進行描述。

圖1 風(fēng)云1C碎片群的一段短弧段Fig.1 A short arc of FENGYUN 1C debris group

表1 風(fēng)云1C碎片群參數(shù)Table 1 Parameters of FENGYUN 1C debris group

NORADIDSemi?majoraxis/kmEccentricityInclination/(°)297277022．560．03098．83297287023．720．03199．18297297024．550．02998．93297307025．640．03098．88297317025．930．03098．93297327025．660．03198．89

1.2.1 群幾何結(jié)構(gòu)

由圖1可以看出，空間碎片群中目標(biāo)分布密集，相互之間可能由于間距很近而產(chǎn)生遮擋或是無法分辨。但是在運動的這段時間里，6個空間碎片在運動過程中保持了一個相對穩(wěn)定的位置關(guān)系。這是因為密集分布的空間碎片群的軌道模型具有相似性，在短時間內(nèi)碎片群內(nèi)目標(biāo)的運動狀態(tài)相似，空間碎片群的整體結(jié)構(gòu)較為穩(wěn)定。

1.2.2 群內(nèi)目標(biāo)數(shù)

群中目標(biāo)的數(shù)量也是群的一個重要特征。由于群中目標(biāo)的檢測概率并不相同，尺寸大的目標(biāo)可能比尺寸小的目標(biāo)檢測概率高，且檢測概率與目標(biāo)材質(zhì)、目標(biāo)翻滾、目標(biāo)姿態(tài)等因素有關(guān)。當(dāng)群內(nèi)目標(biāo)檢測概率普遍較低時，要獲得準(zhǔn)確的群目標(biāo)數(shù)估計比較困難。

1.2.3 群內(nèi)顯著目標(biāo)

上面已經(jīng)提到群內(nèi)目標(biāo)的檢測概率有區(qū)別，而空間碎片大多尺寸較小，群內(nèi)可能大部分目標(biāo)的檢測概率都較低，想要估計各個目標(biāo)的狀態(tài)比較困難。為了對群運動狀態(tài)進行描述，可以通過群內(nèi)檢測概率較高的目標(biāo)運動狀態(tài)進行描述，這里將其稱為群內(nèi)顯著目標(biāo)。通過群內(nèi)顯著目標(biāo)的狀態(tài)來對群進行標(biāo)識。從表1中可以看到，群內(nèi)目標(biāo)的軌道根數(shù)十分相近，即群內(nèi)目標(biāo)運動狀態(tài)相似，通過一個容易估計的顯著目標(biāo)的運動狀態(tài)來代表一個群的狀態(tài)是合理的。

1.3 空間碎片群狀態(tài)估計方法

在前面敘述中已經(jīng)提到一種可行的解決方案是先嘗試對群整體進行濾波，然后在信息足夠的情況下再對單個目標(biāo)進行處理。

群幾何可以通過群內(nèi)目標(biāo)的分布進行描述?；赗FS理論的CPHD濾波器在濾波過程中傳遞整個場景的目標(biāo)密度分布，通過目標(biāo)密度可以方便地顯示群的結(jié)構(gòu)特征。由于檢測概率低等因素可能導(dǎo)致估計不準(zhǔn)確，但是空間碎片群運動具有穩(wěn)定性，通過一段時間的信息積累可以獲得群幾何特征。

在對群內(nèi)目標(biāo)數(shù)的估計上面，傳統(tǒng)方法沒有集成目標(biāo)數(shù)估計功能，只有認為存在的軌跡數(shù)即目標(biāo)數(shù)，由于不穩(wěn)定的軌跡起始和軌跡終止過程，目標(biāo)數(shù)的估計結(jié)果不穩(wěn)定。CPHD濾波器在濾波過程中傳遞的目標(biāo)數(shù)分布，可以容易地對目標(biāo)數(shù)進行估計。

最后是對群內(nèi)顯著目標(biāo)的運動狀態(tài)進行估計。CPHD濾波器沒有顯式集成數(shù)據(jù)關(guān)聯(lián)過程，因此無法提取顯著目標(biāo)的軌跡。但是由于顯著目標(biāo)檢測概率最大，因此其所在位置的目標(biāo)密度會相對持續(xù)地維持較高狀態(tài)，可以以此從整個群的目標(biāo)密度分布中提取顯著目標(biāo)狀態(tài)，但需要對濾波器進行改進。通過對各量測更新的目標(biāo)密度進行標(biāo)記，然后關(guān)注檢測概率高的部分，最后取檢測概率最高的目標(biāo)狀態(tài)作為群顯著目標(biāo)狀態(tài)。

空間碎片群狀態(tài)估計過程框圖如圖2所示，包含了濾波和數(shù)據(jù)處理兩個流程，其中的詳細過程將在后面介紹。

圖2 碎片群狀態(tài)估計過程框圖Fig.2 Block diagram of debris group state estimation

2 空間碎片群的隨機有限集表示

2.1 空間碎片群的多目標(biāo)貝葉斯濾波

在低軌空間碎片群的濾波過程中，因為群內(nèi)目標(biāo)數(shù)未知，觀測到的量測數(shù)也是未知的，因此將空間碎片群內(nèi)目標(biāo)狀態(tài)和量測建模成隨機集形式。在時刻k,碎片群目標(biāo)狀態(tài)集合和量測集合分別建模成隨機有限集Xk={x1,x2,…,xnk}和Zk={z1,z2,…,zmk},xi和zi分別為單目標(biāo)狀態(tài)和單個量測狀態(tài)。空間碎片群多目標(biāo)貝葉斯濾波的預(yù)測和更新過程[8]為

(2)

(3)

式中：fk|k-1(Xk)和fk(Xk)分別為多目標(biāo)先驗和后驗密度；fk|k-1(Xk|·)為多目標(biāo)狀態(tài)轉(zhuǎn)移函數(shù)；g(Zk|·)為多目標(biāo)似然函數(shù)。積分過程是集積分，定義為

(4)

式中：f({x1,x2,…,xi})為多目標(biāo)分布函數(shù)。其中其遍歷分布fn(x1,x2,…,xn)的n!全排列fn(x1,x2,…,xn)的表達式為

f({x1,x2,…,xi})=

(5)

式中：ρ(n)為勢(即目標(biāo)數(shù))分布函數(shù)；θ表示一個排列，θ的和即遍歷函數(shù)的n!個全排列。

2.2 CPHD濾波器

由于集積分在很多條件下無法實現(xiàn)，可以通過對后驗多目標(biāo)分布進行假設(shè)，使得集積分形式簡化。CPHD濾波器基于獨立同分布群過程假設(shè)，在濾波過程中傳遞PHD和勢分布函數(shù)，可以獲得穩(wěn)定的勢估計結(jié)果。但是CPHD濾波器仍存在一些缺陷，且碎片群狀態(tài)估計有其特定要求，因此需要對CPHD濾波器進行改進。下面給出CPHD濾波過程。存活概率和檢測概率分別用pS(·)和pD(·)表示，泊松雜波密度函數(shù)為κ(·)。詳細的CPHD濾波器推導(dǎo)過程見文獻[13，21]。

預(yù)測假設(shè)時刻k-1多目標(biāo)后驗的PHD為vk-1(x)。預(yù)測的多目標(biāo)PHDvk|k-1(x)由存活多目標(biāo)PHDvS,k|k-1(x)和新生多目標(biāo)PHDvB,k(x)組成：

vk|k-1(x)=vS,k|k-1(x)+vB,k(x)

(6)

vS,k|k-1(x)=〈pS(·)f(x|·),vk-1(·)〉

(7)

預(yù)測的勢分布ρk|k-1(n)為

ρk|k-1(n)=

(8)

更新如果多目標(biāo)先驗PHD形式為vk|k-1(x),則多目標(biāo)后驗PHDvk(x|Z)和更新勢分布分別為

(9)

(10)

式中：

(11)

Ξk(φ,Z)={〈φz(·),1〉:z∈Z}

(12)

(13)

(14)

并規(guī)定e0(Z)=1。

遺留PHD，即漏檢目標(biāo)的PHD為

vL,k(x)=

(15)

可以看出，CPHD濾波過程不需要顯式的數(shù)據(jù)關(guān)聯(lián)過程，將整個空間碎片群作為整體處理，假設(shè)空間碎片群滿足獨立同分布過程。正是因為這種不加區(qū)分的處理，在CPHD濾波器計算出漏檢目標(biāo)PHD后直接將其分配到整個場景。如式(15)所示，式(15)中計算的漏檢目標(biāo)PHD直接分配到整個場景的預(yù)測PHDvk|k-1(x)上。因此，漏檢的PHD減少，導(dǎo)致當(dāng)連續(xù)漏檢時，可能出現(xiàn)目標(biāo)丟失。由于空間碎片尺寸小，檢測概率低，漏檢在空間碎片群濾波時將會十分常見，這個問題必須解決。同時，由于不區(qū)分各個目標(biāo)，無法提取目標(biāo)軌跡。但是為了估計群顯著目標(biāo)狀態(tài)，需要提取顯著目標(biāo)軌跡。根據(jù)這些實際需求，本文提出改進的CPHD濾波器。

3 面向量測的改進CPHD濾波器

3.1 按量測聚合的目標(biāo)密度分布

在CPHD濾波器中，每個量測是通過所有軌跡進行聯(lián)合更新[8]，濾波過程可以認為是面向軌跡的濾波過程。將濾波過程中的PHD成分按量測進行聚合，并用一個二元組進行標(biāo)記。二元組標(biāo)識設(shè)為c=(k,m),k為量測出現(xiàn)的時間,m為此量測區(qū)別于該時刻其他量測的序號。在濾波器更新之后，按同一個量測標(biāo)記的PHD成分，即標(biāo)記相同c的PHD成分被認為屬于同一個軌跡，通過這樣的方式使得PHD成分被區(qū)分開來。PHD也由原來的v(x)擴展為v(x,c)。

3.2 面向量測的關(guān)聯(lián)過程

對某一區(qū)域PHD的積分是對該區(qū)域目標(biāo)數(shù)的估計，這里使用權(quán)重w來表述對某PHD成分的積分，即w=〈v(·),1〉。令wi,j表示從軌跡ci到量測zi∈Zk的權(quán)重貢獻，通過式(16)計算:

wi,jw(ci;zj)=〈v(·,ci;zj),1〉

(16)

量測更新權(quán)重由式(17)計算:

(17)

軌跡更新權(quán)重則由式(18)計算:

(18)

在獲得了權(quán)重貢獻矩陣之后，大部分的權(quán)重可能會集中在矩陣中的部分位置，而有許多位置的權(quán)重值可能小到可以忽略。因此，通過下面的過程可以對矩陣進行簡化，通過這種簡化也可減少后面的數(shù)據(jù)關(guān)聯(lián)過程產(chǎn)生的關(guān)聯(lián)結(jié)果。

上一時刻的軌跡通過權(quán)重貢獻和該時刻觀測到的量測進行關(guān)聯(lián)。權(quán)重貢獻越大，說明量測與軌跡之間關(guān)聯(lián)的可能性越大，反之亦然。當(dāng)權(quán)重貢獻大于一個預(yù)先設(shè)定的閾值T1時，認為該量測來自該軌跡，其他關(guān)聯(lián)就不再考慮；反之，當(dāng)權(quán)重貢獻小于一個預(yù)先設(shè)定的閾值T2時，認為該量測必定不是來自該軌跡，權(quán)重貢獻設(shè)置為0。然后，所有的關(guān)聯(lián)結(jié)果都保留下來，形成多條可能軌跡。但是其中關(guān)聯(lián)到同一個量測的軌跡是互斥軌跡，即這兩條軌跡不可能同時存在，互斥軌跡在最后的結(jié)果中只能保留一條。

3.3 漏檢目標(biāo)的PHD處理

在CPHD濾波器中，遺留PHD的分配結(jié)果不正確。在通過上述對PHD成分進行標(biāo)記之后，可以通過標(biāo)記對遺留PHD進行正確分配。

先考慮一個單獨軌跡。一條軌跡的消失權(quán)重wD(c)加上漏檢權(quán)重wL(c),再加上存活權(quán)重(即檢測更新權(quán)重)wU(c),等于1，即wD(c)+wL(c)+wU(c)=1。設(shè)該軌跡的預(yù)測權(quán)重為wk|k-1(c),則預(yù)測的消失權(quán)重為1-wk|k-1(τ)。設(shè)一個常數(shù)的檢測概率為pD,則該目標(biāo)漏檢的權(quán)重為存活權(quán)重wk|k-1(c)與漏檢概率1-pD的乘積，即(1-pD)wk|k-1(c)。

因此，在該目標(biāo)未檢測到的情況下，目標(biāo)消失和目標(biāo)漏檢權(quán)重的比值滿足：

(19)

則漏檢權(quán)重可以表示為

(20)

而軌跡的概率分布可以通過式(21)計算:

(21)

整個場景中遺留PHD的積分就是漏檢目標(biāo)數(shù)估計NL,k=〈vL,k(·),1〉。則最終通過將NL,k按各軌跡漏檢權(quán)重進行按比例分配，得到修正后的各軌跡的遺留PHD，計算公式為

(22)

3.4 新出現(xiàn)碎片模型

通過這樣的模型，來自兩個時刻的所有量測都可以建模為新生目標(biāo)狀態(tài)，但是其中只有少量相距較近的兩個量測建立的新生目標(biāo)狀態(tài)是合理的。這里通過軌道能量模型進行約束，剔除不合理的新生目標(biāo)狀態(tài)。

先假設(shè)碎片軌道是圓軌道，因為一般碎片運動的軌道偏心率都很小。則速度大小可以通過式(23)計算:

(23)

式中：μ為引力常數(shù)和地球質(zhì)量的乘積。在考慮了量測誤差的情況下，如果初始的速度大小‖v1‖大于兩倍的‖v*‖,就認為此速度是不合理的，這個新生目標(biāo)狀態(tài)被剔除。

同時，一個量測如果來自已存在目標(biāo)，則不是新目標(biāo)。因此量測來自新目標(biāo)的概率與其更新權(quán)重w(z)成反比，此概率寫為1-w(z)。還可以設(shè)定一個閾值，如果小于閾值則該量測不進行新目標(biāo)起始。

3.5 狀態(tài)估計過程

4 群狀態(tài)估計數(shù)據(jù)處理過程

在完成了實時的空間碎片群濾波之后，需要對濾波后的數(shù)據(jù)進行處理，從中獲取所需要的群狀態(tài)信息。

4.1 軌跡提取

在改進的CPHD濾波過程中，面向量測的聚合PHD過程使得各個時刻的量測可以關(guān)聯(lián)起來。在濾波結(jié)束之后，可以嘗試從各時刻濾波結(jié)果中提取各空間碎片的軌跡。當(dāng)檢測概率很低時，可能只有部分目標(biāo)能提取軌跡，比如顯著目標(biāo)，通過提取的顯著目標(biāo)軌跡可以對其進行狀態(tài)估計和定軌，用來描述整個碎片群狀態(tài)。

在碎片目標(biāo)分布十分密集的情況下，關(guān)聯(lián)結(jié)果會出現(xiàn)多種可能，最后產(chǎn)生多條軌跡。這里將每條可能軌跡的所有權(quán)重貢獻的乘積作為軌跡的權(quán)重，如出現(xiàn)漏檢則乘以漏檢權(quán)重wL(c)。最后，取某碎片的多條軌跡中權(quán)重最大的軌跡作為碎片軌跡。其中如果存在互斥軌跡，則互斥軌跡中選擇權(quán)重最大的軌跡，其他軌跡舍棄。軌跡長度須滿足一定的連續(xù)時間步長，而軌跡數(shù)可能小于目標(biāo)數(shù)估計。

4.2 顯著目標(biāo)狀態(tài)估計

由于檢測概率低等因素導(dǎo)致提取的軌跡數(shù)可能小于目標(biāo)數(shù)估計，在這種情況下，對檢測概率較高的目標(biāo)進行狀態(tài)估計的結(jié)果會相對穩(wěn)定。因此，提取檢測概率較高目標(biāo)的狀態(tài)估計結(jié)果作為群顯著目標(biāo)狀態(tài)。

在進行了面向量測的關(guān)聯(lián)過程之后，可以從濾波結(jié)果中提取軌跡。提取的航跡長度不能太短，太短沒有實際意義。可以根據(jù)實際情況設(shè)定長度閾值，大于閾值的航跡才有參考價值。航跡檢測概率可以通過式(24)進行計算:

(24)

式中：nU為檢測更新次數(shù)；nL為漏檢次數(shù)。

然后取其中檢測概率最高的目標(biāo)作為群的顯著目標(biāo)，通過其軌跡可以對其進行狀態(tài)估計和定軌，結(jié)果作為群的一個特征。在對空間碎片群進行識別的過程中，顯著目標(biāo)的狀態(tài)(如軌道根數(shù))就可以作為一個識別群的依據(jù)。

5 仿真實驗

本文提出的改進CPHD濾波器使用高斯混合(Gaussian Mixture，GM)實現(xiàn)[22]，由于目標(biāo)運動模型是復(fù)雜的非線性模型，使用不敏卡爾曼濾波器(Unscented Kalman Filter，UKF)進行濾波[23]。同時使用經(jīng)典CPHD濾波器和JPDA濾波器作為比較。在仿真實驗中，考慮一個包含4個碎片的低軌空間碎片群目標(biāo)，設(shè)這4個目標(biāo)是在2007年1月11日07:00:00產(chǎn)生。剛產(chǎn)生時的軌道根數(shù)數(shù)據(jù)見表2，其中顯著目標(biāo)是目標(biāo)1。碎片的運動模型已經(jīng)在1.1節(jié)中介紹，而量測模型為

zk=Hk(x)+Rk

(25)

式中：Rk為量測誤差；Hk(·)為非線性的量測方程。雷達可觀測到目標(biāo)相對雷達站址的方位角、俯仰角和距離。雷達的站址設(shè)為

rR=[-1 787.37 5 500.97 2 679.10] km

雷達的觀測區(qū)域為方位角75°～105°、俯仰角0°～90°范圍。在一天時間里，空間碎片群兩次經(jīng)過雷達觀測區(qū)域，此過程示意見圖3。圖3表示雷達在先觀測到一次碎片群后，隨地球自轉(zhuǎn)再次觀測到碎片群的過程。第1次碎片群通過雷達觀測區(qū)域的時間為8:43:04到8:44:05，第2次通過雷達觀測區(qū)域的時間為20:36:47到 20:39:32。

采樣間隔設(shè)置為1 s。雷達的量測誤差考慮[0.1° 0.1° 50 m]和[0.2° 0.2° 100 m]兩種情況，而每次量測的平均雜波數(shù)取nC=10或nC=50。雜波平均分布區(qū)域為

V=[0.4π 0.6π]rad×[0.05π 0.25π]rad×[300 2 300] km

因此平均雜波密度為λC=1.37×10-2(rad2·km)-1或λC=6.85×10-2(rad2·km)-1。檢測概率考慮3種情況，除顯著目標(biāo)外的其他目標(biāo)的檢測概率取pD=0.3，0.6，0.9,顯著目標(biāo)檢測概率則分別取pD=0.5,0.8，0.95。目標(biāo)的存活概率為pS=0.99。CPHD和改進CPHD濾波需要通過修剪和合并過程減少高斯混合分量數(shù)，修剪的閾值為Tp=10-4,合并的閾值為Tm=4,且規(guī)定高斯分量數(shù)不能超過100個。新生目標(biāo)權(quán)重取0.05。另外，改進CPHD濾波過程中的兩個閾值分別取T1=0.9和T2=10-4。由于目標(biāo)檢測概率不同，而濾波器只能設(shè)置固定的檢測概率，改進CPHD對目標(biāo)數(shù)估計可能出現(xiàn)不穩(wěn)定，這里對每個更新量測權(quán)重設(shè)置必須滿足w(z)>0.3才進行狀態(tài)估計。

表2 仿真實驗碎片群參數(shù)Table 2 Parameters of simulated debris group

圖3 碎片群的兩段觀測弧段(紅色)和兩次觀測的 雷達位置(黑點)示意圖 Fig.3 Sketch of two observed arcs of debris group (red line) and two observing position of radar (black dot)

此外，在JPDA濾波器中使用與CPHD和改進CPHD濾波器相似的面向量測的軌跡起始過程，連續(xù)關(guān)聯(lián)兩次就起始一條軌跡；如果一條軌跡連續(xù)3次在確認門限中沒有量測關(guān)聯(lián)則終止此軌跡。JPDA的目標(biāo)數(shù)估計就是存活的軌跡數(shù)。

在改變觀測條件進行仿真時，通過200次蒙特卡羅仿真實驗的平均結(jié)果來檢驗改進的CPHD濾波器的性能，每次蒙特卡羅實驗中碎片的軌道不變，隨機產(chǎn)生空間碎片的量測數(shù)據(jù)。圖4～圖9分別顯示了改進CPHD濾波器對群幾何、群內(nèi)目標(biāo)數(shù)和群內(nèi)顯著目標(biāo)狀態(tài)估計的結(jié)果。隨著檢測概率降低，量測誤差增大，以及雜波密度上升，改進的CPHD濾波器表現(xiàn)有所下降，但是都能有效估計出群內(nèi)目標(biāo)分布、目標(biāo)數(shù)，并提取出顯著目標(biāo)狀態(tài)。證實了在觀測條件比較惡劣的情況下，所提方法對空間碎片群運動狀態(tài)估計取得較好的效果。而且，在相同條件下，CPHD濾波器和JPDA濾波器也進行了仿真，通過比較來觀察改進CPHD的效果。下面詳細介紹各仿真實驗結(jié)果。

在對群幾何結(jié)構(gòu)的估計中， 使用最優(yōu)子模式分配(Optimal Sub-Pattern Assignment，OSPA)[24]距離來評估其估計結(jié)果。OSPA距離是用來評價估計結(jié)果的，它在考慮目標(biāo)分布準(zhǔn)確性的同時，也考慮了目標(biāo)數(shù)估計的準(zhǔn)確性，是合理的誤差距離準(zhǔn)則。200次蒙特卡羅實驗的平均OSPA距離如圖4和圖5所示。

圖4 弧段1的200次蒙特卡羅平均OSPA距離Fig.4 Mean OSPA distance over 200 Monte Carlo runs of the first arc

圖5 弧段2的200次蒙特卡羅平均OSPA距離Fig.5 Mean OSPA distance over 200 Monte Carlo runs of the second arc

從圖4和圖5中可以看出，觀測條件越惡劣，CPHD濾波器、改進CPHD濾波器和JPDA濾波器的濾波表現(xiàn)越差，但是總體上改進CPHD濾波器的表現(xiàn)都比CPHD和JPDA濾波器好很多，說明改進CPHD濾波器對群內(nèi)目標(biāo)分布估計的精度和穩(wěn)定性比CPHD和JPDA濾波器好。在檢測概率為0.9的情況下，3個濾波器都能比較好地工作。在檢測概率降低到0.6時，3個濾波器的表現(xiàn)都出現(xiàn)大幅度下降。由于CPHD濾波器無法正確分配漏檢目標(biāo)的PHD，CPHD濾波器表現(xiàn)下降很快。當(dāng)檢測概率進一步降低到0.3時，CPHD和JPDA濾波器已經(jīng)無法工作，而改進CPHD濾波器仍能工作，展示了改進CPHD濾波器的穩(wěn)定表現(xiàn)，也證實了其對漏檢目標(biāo)PHD的正確處理。

圖6 弧段1的200次蒙特卡羅目標(biāo)數(shù)估計Fig.6 Mean cardinality estimate over 200 Monte Carlo runs of the first arc

圖7 弧段2的200次蒙特卡羅目標(biāo)數(shù)估計Fig.7 Mean cardinality estimate over 200 Monte Carlo runs of the second arc

200次蒙特卡羅實驗的平均目標(biāo)數(shù)估計如圖6和圖7所示，其中不同顏色虛線為對應(yīng)的平均估計標(biāo)準(zhǔn)差范圍，兩種濾波器在觀測條件不斷惡化的時候標(biāo)準(zhǔn)差都在不斷增大，即估計結(jié)果變得不穩(wěn)定。從圖6和圖7中可以看出，3個濾波器在檢測概率為0.9時能估計出群中包含了4個目標(biāo)，其中CPHD在弧段2后半段出現(xiàn)估計偏差。當(dāng)量測誤差增大時，CPHD和JPDA濾波器的表現(xiàn)有所下降。當(dāng)檢測概率降低到0.6時，CPHD和JPDA濾波器對目標(biāo)數(shù)的估計出現(xiàn)偏差，已經(jīng)無法準(zhǔn)確估計目標(biāo)數(shù)，而改進CPHD濾波器仍能估計出4個目標(biāo)；當(dāng)檢測概率降低到0.3時，JPDA濾波器的估計已經(jīng)無效，CPHD濾波器的表現(xiàn)優(yōu)于JPDA濾波器，但也無法準(zhǔn)確估計目標(biāo)數(shù)，而改進CPHD濾波器仍能在一段時間后估計出有4個目標(biāo)，說明其在檢測概率很低的情況下仍能工作，只是估計結(jié)果的標(biāo)準(zhǔn)差(起伏)變大。

圖8 弧段1的200次蒙特卡羅群顯著目標(biāo)估計 均方根誤差Fig.8 RMSE of conspicuous debris object over 200 Monte Carlo runs of the first arc

圖9 弧段2的200次蒙特卡羅群顯著目標(biāo)估計 均方根誤差(改進CPHD)Fig.9 RMSE of conspicuous debris object over 200 Monte Carlo runs of the second arc (improved CPHD)

200次蒙特卡羅實驗對群顯著目標(biāo)估計的均方根誤差(RMSE)如圖8和圖9所示。從圖8和圖9中可以看出，隨著量測誤差的增大，對顯著目標(biāo)RMSE的估計也越來越大；而檢測概率越低，RMSE誤差減小得越慢，且估計精度也缺差。但是，即使是在檢測概率很低的情況下，改進CPHD濾波器仍能估計出較為準(zhǔn)確的結(jié)果，可以對群內(nèi)的顯著目標(biāo)狀態(tài)給出較為準(zhǔn)確的估計，從而通過顯著目標(biāo)狀態(tài)來了解群狀態(tài)就更加可信，也可以通過此狀態(tài)來對群進行標(biāo)識和分辨。

6 結(jié) 論

1) 在對面空間碎片群時，由于低檢測概率、群內(nèi)目標(biāo)密集分布且目標(biāo)數(shù)未知，傳統(tǒng)單目標(biāo)處理方法已經(jīng)無法有效工作。

2) 改進的CPHD濾波器解決了漏檢目標(biāo)PHD分配問題，且集成軌跡標(biāo)識過程，使其可以提取群內(nèi)顯著目標(biāo)狀態(tài)。

3) 仿真實驗通過多種不同的觀測條件，驗證了本文提出的方法是一種對空間碎片群運動狀態(tài)估計的有效方法。

因為實際的空間碎片群觀測更加復(fù)雜，未來的工作需要將濾波過程推廣到檢測概率未知的條件下去，同時通過定軌和軌道改進等方式對空間碎片群以及群內(nèi)目標(biāo)不同弧段的觀測進行弧段關(guān)聯(lián)。

[1] 劉靜, 王榮蘭, 張宏博, 等. 空間碎片碰撞預(yù)警研究[J]. 空間科學(xué)學(xué)報, 2004, 24(6): 462-469.

LIU J, WANG R L, ZHANG H B, et al. Space debris collision prediction research[J]. Chinese Journal of Space Science, 2004, 24(6): 462-469 (in Chinese).

[2] 李春來, 歐陽自遠, 都亨. 空間碎片與空間環(huán)境[J]. 第四紀研究, 2002, 22(6): 540-551.

LI C L, OUYANG Z Y, DU H. Space debris and space environment[J]. Quaternary Sciences, 2002, 22(6): 540-551 (in Chinese).

[3] JOHNSON N L, STANSBERY E, LIOU J C, et al. The characteristics and consequences of the break-up of the Fengyun-1C spacecraft[J]. Acta Astronautica, 2008, 63(1-4): 128-135.

[4] PARDINI C, ANSELMO L. Physical properties and long-term evolution of the debris clouds produced by two catastrophic collisions in Earth orbit[J]. Advances in Space Research, 2011, 48(3): 557-569.

[5] BLACKMAN S S. Multiple hypothesis tracking for multiple target tracking[J]. IEEE Aerospace & Electronic Systems Magazine, 2004, 19(1): 5-18.

[6] FORTMANN T E, BAR-SHALOM Y, SCHEFFE M. Sonar tracking of multiple targets using joint probabilistic data association[J]. IEEE Journal of Oceanic Engineering, 1983,8(3): 173-184.

[7] HUANG J, HU W D. MCMC-particle-based group tracking of space objects within Bayesian framework[J]. Advances in Space Research, 2014, 53(2): 280-294.

[8] MAHLER R P S. Statistical multisource-multitarget information fusion[M]. Norwood, MA: Artech House, 2007: 305-538.

[9] WEI B, NENER B, LIU W. Tracking of space debris via CPHD and consensus[C]∥International Conference on Control, Automation and Information Sciences, 2015: 436-441.

[10] JONES B A, VO B N. A labeled multi-Bernoulli filter for space object tracking[C]∥AAS/AIAA Space Flight Mechanics Meeting, 2014.

[11] MAHLER R P S. Advances in statistical multisource-multitarget information fusion[M]. Norwood, MA: Artech House, 2014: 181-215,379-403

[12] MAHLER R P S. Multitarget Bayes filtering via first-order multitarget moments[J]. IEEE Transactions on Aerospace & Electronic Systems, 2003, 39(4): 1152-1178.

[13] MAHLER R. PHD filters of higher order in target number[J]. IEEE Transactions on Aerospace & Electronic Systems, 2007, 43(4): 1523-1543.

[14] MAHLER R P S, ZAJIC T. Bulk multitarget tracking using a first-order multitarget moment filter[J]. Proceedings of SPIE-The International Society for Optical Engineering, 2002, 4729: 175-186.

[15] ERDINC O, WILLETT P, BAR-SHALOM Y. Probability hypothesis density filter for multitarget multisensor tracking[C]∥International Conference on Information Fusion, 2005: 1-8.

[16] FRANKEN D, SCHMIDT M, ULMKE M. “Spooky action at a distance” in the cardinalized probability hypothesis density filter[J]. IEEE Transactions on Aerospace & Electronic Systems, 2009, 45(4): 1657-1664.

[17] PANTA K, CLARK D E, VO B N. Data association and track management for the gaussian mixture probability hypothesis density filter[J]. IEEE Transactions on Aerospace & Electronic Systems, 2009, 45(3): 1003-1016.

[18] 歐陽成, 姬紅兵, 張俊根. 一種改進的CPHD多目標(biāo)跟蹤算法[J]. 電子與信息學(xué)報, 2010, 32(9): 2112-2118.

OUYANG C, JI H B, ZHANG J G. Improved CPHD filter for multitarget tracking[J]. Journal of Electronics & Information Technology, 2010, 32(9): 2112-2118 (in Chinese).

[19] 劉林. 航天器軌道理論[M]. 北京： 國防工業(yè)出版社, 2000: 99-326.

LIU L. Orbit theory of spacecraft[M]. Beijing: National Defense Industry Press, 2000: 99-326 (in Chinese).

[20] VALLADO D, CRAWFORD P, HUJSAK R, et al. Revisiting spacetrack report #3[C]∥AIAA/AAS Astrodynamics Specialist Conference and Exhibit. Reston, VA: AIAA, 2006: 1-88.

[21] VO B T, VO B N, CANTONI A. Analytic implementations of the cardinalized probability hypothesis density filter[J]. IEEE Transactions on Signal Processing, 2007, 55(7): 3553-3567.

[22] VO B N, MA W K. The Gaussian mixture probability hypothesis density filter[J]. IEEE Transactions on Signal Processing, 2006, 54(11): 4091-4104.

[23] JULIER S J, UHLMANN J K. Unscented filtering and nonlinear estimation[J]. Proceedings of the IEEE, 2004, 92(3): 401-422.

[24] SCHUHMACHER D, VO B T, VO B N. A Consistent metric for performance evaluation of multi-object filters[J]. IEEE Transactions on Signal Processing, 2008, 56(8): 3447-3457.

Stateestimationofspacedebrisgroupbasedonrandomfiniteset

LUZhejun,HUWeidong*

ATRKeyLab,CollegeofElectronicScienceandEngineering,NationalUniversityofDefenseTechnology,Changsha410073,China

Basedonthesingletargetstateestimation,theconventionalapproachisnotabletoworkwellwhenfacedwithalargenumberofsuddenlygeneratedspacedebrisobjects,asthoseobjectsareclosely-spacedasagroupwithsmallsize.Thus,basedontheRandomFiniteSet(RFS)theory,thespacedebrisgroupistreatedastheprocessingobjectanditsstatesareestimatedinthiswork.Inordertoaddresstheissuesofmissedobjectdensitydistributionandtrajectoryassociation,animprovedmeasurement-orientedCardinalizedProbabilityHypothesisDensity(CPHD)filterisproposed.Withadataprocessingusedafterfiltering,thisfilteraccomplishestheestimationofobjectdensitydistribution,objectnumberandconspicuousobjectstateinagroup.Insimulations,theproposedfiltersignificantlyoutperformstheconventionalfilterandCPHDfilter.Itcanworkinchallengingenvironment,andmeanwhile,theconventionalfilterandCPHDfilterfail.

spacedebris;groupobject;stateestimation;randomfiniteset;improvedCPHDfilter

2017-02-28；Revised2017-03-28；Accepted2017-04-27；Publishedonline2017-04-281648

URL：http://hkxb.buaa.edu.cn/CN/html/20171124.html

NationalNaturalScienceFoundationofChina(61372162)

.E-mailwdhu@nudt.edu.cn

http://hkxb.buaa.edu.cnhkxb@buaa.edu.cn

10.7527/S1000-6893.2017.321200

V412

A

1000-6893(2017)11-321200-11

2017-02-28；退修日期2017-03-28；錄用日期2017-04-27；< class="emphasis_bold">網(wǎng)絡(luò)出版時間

時間：2017-04-281648

http://hkxb.buaa.edu.cn/CN/html/20171124.html

國家自然科學(xué)基金(61372162)

.E-mailwdhu@nudt.edu.cn

盧哲俊， 胡衛(wèi)東. 基于隨機有限集的空間碎片群運動狀態(tài)估計方法J. 航空學(xué)報，2017,38(11):321200.LUJZ,HUWD.StateestimationofspacedebrisgroupbasedonrandomfinitesetJ.ActaAeronauticaetAstronauticaSinica,2017,38(11):321200.

(責(zé)任編輯：蘇磊)