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

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

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

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

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

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

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

隨著空間碎片數(shù)量的不斷增加，對空間活動安全造成嚴重威脅，對空間碎片的監(jiān)測就顯得尤為重要。但是空間碎片數(shù)量大、尺寸小導致監(jiān)測困難，尤其是由于解體或碰撞產(chǎn)生的高密度空間碎片集群成為空間監(jiān)測和編目的難點[1-2]。例如2007年1月風云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]。空間碎片在剛產(chǎn)生之后的短時間里相距很近，形成了高密度的空間碎片“云”，且這些碎片中的絕大部分尺寸較小，觀測困難。這時候，地基雷達系統(tǒng)由于檢測概率低，目標密集無法分辨，以及目標數(shù)量大且未知等因素導致無法對空間碎片進行有效監(jiān)測和編目。

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

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

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

本文結(jié)構如下：第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為模型噪聲。

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

1.2 典型空間碎片群特征

本文關注的空間碎片群指的是低軌空間突發(fā)產(chǎn)生的大量密集空間碎片組成的集群。這里以風云1C空間碎片群為例，對空間碎片群的特征進行描述。取6個密集分布的風云1C碎片目標組成的集群為研究目標，研究數(shù)據(jù)采用北美空間防空司令部(NORAD)公布的雙行根數(shù)(Two-Line Element，TLE)，這6個目標的衛(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 風云1C碎片群的一段短弧段Fig.1 A short arc of FENGYUN 1C debris group

表1 風云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é)構

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

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

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

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

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

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

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

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

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

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

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

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

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

2.1 空間碎片群的多目標貝葉斯濾波

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

(2)

(3)

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

(4)

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

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

(5)

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

2.2 CPHD濾波器

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

預測假設時刻k-1多目標后驗的PHD為vk-1(x)。預測的多目標PHDvk|k-1(x)由存活多目標PHDvS,k|k-1(x)和新生多目標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)

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

ρk|k-1(n)=

(8)

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

(9)

(10)

式中：

(11)

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

(12)

(13)

(14)

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

遺留PHD，即漏檢目標的PHD為

vL,k(x)=

(15)

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

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

3.1 按量測聚合的目標密度分布

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

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

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

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

(16)

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

(17)

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

(18)

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

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

3.3 漏檢目標的PHD處理

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

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

因此，在該目標未檢測到的情況下，目標消失和目標漏檢權重的比值滿足：

(19)

則漏檢權重可以表示為

(20)

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

(21)

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

(22)

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

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

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

(23)

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

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

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

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

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

4.1 軌跡提取

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

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

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

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

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

(24)

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

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

5 仿真實驗

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

zk=Hk(x)+Rk

(25)

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

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。

采樣間隔設置為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種情況，除顯著目標外的其他目標的檢測概率取pD=0.3，0.6，0.9,顯著目標檢測概率則分別取pD=0.5,0.8，0.95。目標的存活概率為pS=0.99。CPHD和改進CPHD濾波需要通過修剪和合并過程減少高斯混合分量數(shù)，修剪的閾值為Tp=10-4,合并的閾值為Tm=4,且規(guī)定高斯分量數(shù)不能超過100個。新生目標權重取0.05。另外，改進CPHD濾波過程中的兩個閾值分別取T1=0.9和T2=10-4。由于目標檢測概率不同，而濾波器只能設置固定的檢測概率，改進CPHD對目標數(shù)估計可能出現(xiàn)不穩(wěn)定，這里對每個更新量測權重設置必須滿足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ù)關聯(lián)兩次就起始一條軌跡；如果一條軌跡連續(xù)3次在確認門限中沒有量測關聯(lián)則終止此軌跡。JPDA的目標數(shù)估計就是存活的軌跡數(shù)。

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

在對群幾何結(jié)構的估計中， 使用最優(yōu)子模式分配(Optimal Sub-Pattern Assignment，OSPA)[24]距離來評估其估計結(jié)果。OSPA距離是用來評價估計結(jié)果的，它在考慮目標分布準確性的同時，也考慮了目標數(shù)估計的準確性，是合理的誤差距離準則。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)目標分布估計的精度和穩(wěn)定性比CPHD和JPDA濾波器好。在檢測概率為0.9的情況下，3個濾波器都能比較好地工作。在檢測概率降低到0.6時，3個濾波器的表現(xiàn)都出現(xiàn)大幅度下降。由于CPHD濾波器無法正確分配漏檢目標的PHD，CPHD濾波器表現(xiàn)下降很快。當檢測概率進一步降低到0.3時，CPHD和JPDA濾波器已經(jīng)無法工作，而改進CPHD濾波器仍能工作，展示了改進CPHD濾波器的穩(wěn)定表現(xiàn)，也證實了其對漏檢目標PHD的正確處理。

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

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

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

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

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

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

6 結(jié) 論

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

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

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

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

[1] 劉靜, 王榮蘭, 張宏博, 等. 空間碎片碰撞預警研究[J]. 空間科學學報, 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多目標跟蹤算法[J]. 電子與信息學報, 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)絡出版時間

時間：2017-04-281648

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

國家自然科學基金(61372162)

.E-mailwdhu@nudt.edu.cn

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

(責任編輯：蘇磊)