基于幅度信息的海上多目標穩健跟蹤方法

柳超，張志國，孫進平

（1. 北京航空航天大學電子信息工程學院，北京 100191；2. 海軍92853部隊，遼寧 葫蘆島 125106）

1 引言

多目標跟蹤是當前信息融合和計算機視覺領域的熱點問題。檢測的不確定性、量測的不確定性及關聯的不確定性使多目標跟蹤非常復雜[1]，尤其在雷達對海探測環境下，目標速度慢、雜波剩余多等因素使目標連續跟蹤異常困難。目前，常見的多目標跟蹤方法主要包括聯合概率數據關聯（JPDA,joint probabilistic data association）、多假設跟蹤（MHT, multiple hypotheses tracking）、隨機有限集（RFS, random finite set）等[2]。JPDA和MHT的實現過程是先進行數據關聯，再進行單目標濾波。當量測數目增加時，數據關聯運算量呈指數級增長。近年來，基于 RFS的多目標跟蹤方法受到了很多關注，該方法將多目標狀態和傳感器量測分別建模為隨機有限集，通過貝葉斯多目標濾波公式實現多目標后驗概率密度的遞歸估計。與JPDA和MHT等傳統跟蹤方法不同的是，許多RFS方法不考慮數據關聯也可以得到準確的多目標濾波結果，如概率假設密度（PHD, probability hypothesis density）濾波器[3]、勢概率假設密度（CPHD, cardinalized PHD）濾波器[4]、多伯努利（MB, multi-Bernoulli）濾波器[5]。這3種方法由于實現過程簡單、計算效率高等優點獲得了大量應用，尤其是MB濾波器，其計算復雜度和 PHD相似，其序貫蒙特卡羅（SMC,sequential Monte Carlo）實現不需要聚類便能提取目標狀態，因而在多目標濾波中更具優勢。

在基于RFS的多目標濾波器中，雜波強度（即雜波率和雜波概率密度之積）和檢測概率需要根據應用環境設定，但在雷達對海探測中，實時獲取這2個參數非常困難。如果設定的參數值與真實值不一致，濾波器的性能便會下降。近年來，雜波強度和檢測概率未知情況下的目標跟蹤問題得到了深入研究。針對檢測概率未知的情況，Guiseppe等[6]提出2種自適應跟蹤方法；針對雜波強度未知的情況，學者們提出了多種雜波估計算法[7-10]；針對雜波強度和檢測概率均未知的情況，Mahler等[11]提出一種在濾波的同時學習雜波率和檢測概率的自適應CPHD濾波器，隨后，Beard等[12]在此基礎上提出一種自舉濾波器[12]。文獻[11-12]所提方法采用SMC實現時，需要進行聚類以提取目標狀態，當目標的勢估計誤差較大時，狀態估計性能也會變差。針對這一問題，Vo等[13]提出了穩健多伯努利濾波器（RMB, robust MB），因其SMC實現不需要聚類便能提取目標狀態，該方法已應用于傳感器選擇[14]、生物學研究[15]、視頻跟蹤[16]等領域。文獻[17]將穩健濾波思想[13]引入標簽多伯努利濾波器，但沒有給出雜波率的準確估計。

RMB濾波器在計算量測似然時，只利用目標和雜波的運動狀態（如位置、速度等）信息，當目標和雜波相距較近時，多目標濾波效果欠佳。由于雷達回波不僅包含目標的運動狀態信息，還包含幅度信息（AI, amplitude information），利用目標的幅度信息來改進量測似然，提升對目標和雜波的區分能力，實現更加準確的多目標濾波是一個可行的方法。

利用幅度信息的目標跟蹤方法主要包括檢測前跟蹤（TBD, track-before-detect）[18-19]和幅度信息輔助目標跟蹤（AIAOT, AI aided object tracking）[20-29]。前者對連續多幀回波進行非相參積累以提取可能的目標航跡，對低信雜比目標（SCR, signal to clutter ratio）跟蹤效果好，但由于每一幀都要對掃描區域進行遍歷搜索，因此運算量很大；后者是在傳統跟蹤方法的基礎上增加幅度似然的計算，能夠顯著提升算法的跟蹤效果，并且運算量小。

目前，AIAOT領域已取得許多成果。Lerro等[20-21]首先提出了 AIAOT的基本理論，隨后Ehrman等[22]和 Brekke等[23-24]對 AIAOT 進行了深入研究。文獻[20-24]只針對單目標跟蹤場景。對于多目標跟蹤的情況，Clark等[25-26]將Rayleigh雜波的幅度信息引入PHD濾波器，隨后一些學者進行了更加深入的研究[27-28]。針對未知雜波強度的情況，Yuan[29]將 Rayleigh雜波和幅度恒定目標的AI引入RMB濾波器。

在雷達對海探測中，運動目標的姿態和視角相對雷達經常變化，導致目標的雷達散射截面積（RCS,radar cross section）起伏不定，從而造成目標回波幅度起伏。目前，表征目標 RCS起伏的模型包括Swerling模型和χ2分布模型等，工程上常用SwerlingI～SwerlingIV 模型，其中，Swerling I型表示RCS在掃描內的脈沖間相關，而兩次掃描相互獨立，為慢起伏，其概率密度服從瑞利分布。大量海雜波實測數據分析表明，通常海雜波相關時間為毫秒級[30]，服從慢起伏模型，且采用K分布模型能有效擬合較大入射余角范圍內的海雜波幅度分布。目前，Swerling I型已廣泛應用于描述海上目標RCS的起伏[23]。

本文提出了一種K分布海雜波環境下AI輔助的RMB濾波器（AI-RMB）。首先，將K分布海雜波及Swerling I型起伏目標的幅度似然函數引入穩健多伯努利濾波器，提升對海探測環境下雷達的多目標跟蹤性能；其次，與經典 AIAOT方法[20-29]利用幅度信息計算理論檢測概率和理論虛警概率不同，本文考慮到RMB濾波器的特點，采用貝塔分布迭代估計目標伯努利項的檢測概率和雜波伯努利項的虛警概率；最后，為提升計算效率，采用積分表方法計算K分布雜波下的幅度似然函數。仿真實驗表明，相較穩健多伯努利濾波器，本文所提方法在多目標狀態估計、勢估計以及雜波率估計方面性能更優，且運行時間沒有顯著增加。

2 基本模型

本節給出單個目標運動狀態的時間演化模型和雷達的觀測模型。

2.1 目標運動模型

假設k時刻第i個目標為近似勻速運動，其運動模型如式（1）所示。

其中，T和σv分別為采樣間隔和加速度噪聲的功率譜密度。

2.2 雷達觀測模型

k時刻第i個量測矢量由運動狀態相關的量測部分和附加的幅度部分構成。

其中，h（?）為非線性量測函數，為零均值白色高斯量測噪聲，量測函數以及量測噪聲協方差分別為

3 穩健多伯努利濾波器

為了在濾波的同時估計出環境中的雜波率，Vo等[13]提出的RMB濾波器將雜波建模為一類特殊類型的目標（又稱為雜波發生器），并分別建立了目標和雜波的狀態轉移模型和量測模型，通過貝葉斯迭代濾波實現對目標狀態和雜波率的估計。

定義 χ（Δ）= [ 0 ,1]為檢測概率的取值空間，χ~=Rnx為目標的狀態空間，nx為狀態空間的維度，{0 ,1}為雜波發生器和目標的離散標簽空間（本文用u =0表示雜波發生器，u=1表示真實目標），則目標和雜波的增廣狀態空間為目標和雜波的增廣狀態為其中，α為檢測概率，為運動狀態，u為標簽。定義在增廣狀態空間上的任意函數，有

假設在k-1時刻，多目標概率密度可表示為如式（3）所示的多伯努利隨機集。

其中，有

其中，fΔ,u和fχ~,u分別為檢測概率和運動狀態的轉移密度。

將式（4）中預測的多目標密度表示為如式（8）的多伯努利隨機集。

則k時刻后驗多目標密度可表示為漏檢多伯努利隨機集和經量測更新的多伯努利隨機集的并集。

且有

4 基于AI的穩健多伯努利濾波器

在上述RMB迭代中，多目標后驗密度更新時僅利用了與運動狀態相關的量測。本文引入目標和雜波的幅度信息，建立了目標和雜波的幅度量測似然，并將其加入更新過程，從而構成AI-RMB濾波器。由于AI-RMB和RMB的預測部分完全相同，為節省篇幅，這里只給出AI-RMB的更新部分。

4.1 幅度的似然函數

假設目標的幅度d和運動狀態相互獨立，則目標的量測似然函數 g（z|x）和雜波的量測似然函數c（z）可分別表示為[21]

其中，τ為檢測門限。則經過門限檢測后，目標和雜波的幅度似然函數分別為

4.2 AI-RMB濾波器

如果在k時刻，預測的多目標密度可表示為如式（22）所示的多伯努利隨機集。

其中，有

在式（25）、式（26）、式（28）和式（29）中，目標和雜波的檢測概率pD,u,k雖然可以通過式（18）和式（19）求得理論值，但該值不能直接應用于AI-RMB的更新過程。其原因如下：一方面，與僅用于估計多目標狀態的MB濾波器不同，RMB濾波器不僅需要估計目標狀態，還要在濾波的過程中估計出雜波率，由于通常情況下虛警概率比檢測概率低得多（相差幾個數量級），如果采用式（18）和式（19）計算出的理論檢測概率，將會導致雜波伯努利項的權重過低，經過航跡修剪和合并，這些項很可能丟失，從而無法實現對雜波率的估計；另一方面，由式（18）計算出的理論檢測概率并不能準確反映RMB中每個目標伯努利項的檢測概率，這是因為只有在RMB濾波器中的目標伯努利項數與實際目標數完全一致的情況下才能準確表示目標伯努利項的檢測概率，同樣，由式（19）計算出的理論虛警概率也不能準確反映RMB中每個雜波伯努利項的檢測概率。基于上述兩方面原因，本文在計算檢測概率和虛警概率時，摒棄了常規AIAOT[20-29]中所采用的方法，而是采用貝塔分布來迭代估計每個伯努利項的檢測概率和虛警概率。

4.3 K分布海雜波及目標的幅度似然函數

幅度服從K分布的海雜波具有如式（31）所示的概率密度函數（PDF, probability density function）。

其中，G（?）表示Gamma函數，Kv表示第二類修正貝塞爾函數，v是形狀參數，b是尺度參數。則根據式（19）可得海雜波的虛警概率為

其中，τ為檢測門限。

通過K分布雜波的復合解釋[23]，可得Swerling I型起伏目標加雜波的PDF為

則根據式（18）可得目標的檢測概率為

將式（31）和式（32）代入式（21）可得

將式（33）和式（34）代入式（20）可得

其中，有

式（36）中分子和分母的積分都沒有閉合解。針對這一問題，Brekke等[23]提出了一種采用網格的數值積分方法，以對g（η）的積分為例說明如下。

首先進行變量代換，如式（37）所示。

則式（36）中分子部分的積分可以表示為

其中，有

q（u）的極值up可通過對 q（u）求一階導數得到，如式（39）所示。

當u→∞時，q（u）與u-1-2v成正比。為了估計q（u）的積分，使用一個變分辨率的積分網格，該網格包含低區和高區兩部分，分別對應 q（u）極大值附近的積分和u→∞時的積分。對低區網格，采樣點為

對高區網格，采樣點為

這一數值積分方法準確性很高，但是計算量很大，不利于實時應用。為實現式（36）中積分的快速計算，本文制作了積分表（InT, integral table）。當需要計算式（36）的積分時，只需要根據K分布雜波的參數和量測幅度到積分表中查找相應的積分值即可。由于積分表是在算法執行前離線制作的，不消耗運行時間，因而可使海雜波中目標幅度似然的計算效率大大提升。具體操作步驟如下。

步驟1根據雜波參數b、v以及可能的信雜比，設定量測幅度a的可能取值。為了保證積分表在較大的信雜比范圍內可用，a的取值范圍要盡可能大。

步驟2利用Brekke提出的網格法計算雜波參數b、v以及量測a取不同值時的積分值，制作成積分表。

步驟3根據參數b、v以及a的取值，從積分表中提取相應的積分值，完成式（36）中積分的計算。

5 仿真

為檢驗所提AI-RMB的性能，通過仿真實驗將其與MB濾波器和RMB濾波器進行對比。對海雷達觀測區域如下：方位為0°～180°，距離為0～2 000 m。在總計100 s的觀測時間內先后有10個目標進入觀測區域，每個目標都為勻速運動。目標的真實軌跡如圖1所示，其中，○表示航跡起點，△表示航跡終點。

圖1 目標真實軌跡

每一幀中雜波數服從均值為1 000的泊松分布，海雜波幅度服從參數為v=4，b=0.25的K分布，空間分布為方位-距離上的均勻分布，但在直角坐標系中為非均勻分布。所有目標的平均信雜比為

目標的幅度測量值服從Swerling I型起伏。虛警概率設為 0.02，則根據式（32）可求得該信雜比下的檢測門限。當雜波和目標的量測幅度超過檢測門限時生成量測點跡。理論上檢測后雜波的數目服從均值為20的泊松分布，且在方位-距離上服從均勻分布，但在平面直角坐標系上為非均勻分布。檢測后包含目標和雜波的量測點跡如圖2所示。

圖2 量測點跡

目標的存在概率為 pS,1,k=0.99，虛警的存在概率為 pS,0,k=0.9。目標和雜波的檢測概率的轉移密度為如下的貝塔分布。

其中，有

并給定均值 μa,1,k|k-1=ak-1，對于目標，標準差σa,1,k|k-1=0.01，而對于雜波，標準差為σa,0,k|k-1=0.07。

目標的新生項為多伯努利隨機集，參數為

新生目標貝塔分布的參數為 sΓ,1,k=85，tΓ,1,k|=15，運動狀態的概率密度為

雜波的新生項也是多伯努利隨機集，其參數為

本文通過計算每種算法的 OSPA位置誤差和OSPA勢誤差[31]來評價算法的跟蹤性能。

其中，X={ x1,… ,xm}和Y={ y1,… ,yn}分別表示目標的真實狀態和估計狀態，m、n∈?0={ 0 ,1,2…}。本文中設置p=1和c=300。

圖3～圖5為MB、RMB和AI-RMB在某次實驗中的濾波效果。從圖3～圖5可以看到，參數設置完全正確的MB性能最優，只產生了很少的漏檢和虛假航跡；AI-RMB次之，產生了更多的漏檢和虛假航跡；未采用幅度信息的RMB性能最差，產生了大量漏檢和虛假航跡。

圖3 單次實驗中MB濾波結果

圖4 單次實驗中RMB濾波結果

圖5 單次實驗中AI-RMB濾波結果

圖6給出了MB、RMB和AI-RMB進行30次仿真后的平均 OSPA誤差，其中，圖 6（a）是 OSPA位置誤差，圖6（b）是OSPA勢誤差。從圖6可以看到，參數配置正確的MB位置誤差最小；AI-RMB的位置估計和勢估計誤差都明顯小于未采用幅度信息的RMB。這是因為，RMB在計算量測似然時只利用了目標和雜波的運動信息，當目標和雜波距離較近時難以區分，從而影響目標粒子和雜波粒子權重的計算，造成目標狀態估計性能下降。而采用幅度信息可以更好地區分目標和雜波，從而提升目標的狀態估計性能。

圖7是MB、RMB和AI-RMB在30次仿真后的平均勢估計效果。從圖7可以看到，這3種算法都能比較準確地估計出目標數目。但是，當目標數增多時（即第60 s～80 s），未采用幅度信息的RMB產生了比較明顯的欠估計，而采用幅度信息輔助的AI-RMB以及參數配置完全正確的MB對目標數目的估計更加準確。

圖6 30次實驗中平均OSPA位置與勢誤差

圖7 30次實驗中平均勢估計結果

圖8是RMB和AI-RMB對雜波數目的估計結果，其中，圖8（a）為某次仿真中每一幀的真實雜波數目和2種算法估計的雜波數目，圖8（b）為30次仿真后2種算法對雜波數目估計的均方根誤差對比。從圖8可以看出，在開始的大約20幀中，2種算法對雜波數目的估計都存在比較大的誤差，且誤差大小一致，但誤差隨濾波時間的增加而迅速降低。在20幀之后，2種算法的誤差基本上都達到了一個穩定的狀態，容易看到采用幅度信息輔助的AI-RMB的誤差小于未采用幅度信息的RMB。

表1為通過30次仿真實驗統計出的MB、RMB、未采用積分表的AI-RMB和采用積分表的AI-RMB（InT）的平均運行時間。MB、RMB和AI-RMB的計算復雜度均與當前時刻的量測數和目標數成正比，即O（|Zk|×|Xk|）。在本文仿真中，這 3種濾波器處理的量測數是相同的，區別僅在于目標數。從表1可以看到，MB的運行時間最短，這是因為該濾波器僅用于估計真實目標，因而每個時刻的|Xk|都比較小。RMB的運行時間約為MB的6倍，這是因為該濾波器需要同時估計目標和雜波，因此每個時刻的|Xk|都比較大。從仿真設置可以看到，在每個時刻，MB的新生目標項數為4，而RMB為24，為MB的6倍，因此其計算復雜度也應為MB的6倍，這與仿真結果是一致的。理論上，AI-RMB的計算復雜度與RMB相同，但由于幅度似然函數的計算中需要對K分布求積分，而該積分沒有解析表達式，只能通過數值積分求解，因此大大消耗了運行時間。從表1還可以看到，未采用積分表的AI-RMB的運行時間遠遠大于 RMB，而采用積分表的AI-RMB（InT）的運行時間則降低到接近RMB的水平。

圖8 30次實驗中平均雜波數估計和均方誤差

表1 30次實驗中平均運行時間對比

6 結束語

穩健多伯努利濾波器適用于雜波強度和檢測概率未知的場景，因而在雷達對海探測環境下具有重要的應用價值。但該濾波器僅利用與目標運動狀態相關的量測計算似然函數，當目標和雜波距離較近時性能不佳。針對這一問題，本文基于K分布海雜波和 Swerling I型起伏目標建立了幅度似然函數，并將其引入該濾波器的更新過程，以提升算法對目標和雜波的區分能力；在配置檢測概率時，舍棄常規方法中利用目標和雜波幅度計算理論檢測概率和理論虛警概率的思路，采用貝塔分布迭代估計每個目標伯努利項的檢測概率和每個雜波伯努利項的虛警概率。仿真實驗表明，本文所提方法在目標狀態估計、勢估計及雜波數估計方面均優于穩健多伯努利濾波器。且通過使用積分表，所提方法的計算效率與穩健多伯努利濾波器相當。