考慮微多普勒特征的鳥類目標(biāo)多假設(shè)跟蹤算法?

（海軍工程大學(xué)電子工程學(xué)院 武漢 430033）

1 引言

鳥擊是指航空器飛行和起降過(guò)程中和鳥類、蝙蝠等飛行物相撞的事件。鳥擊給軍隊(duì)和民航部門帶來(lái)巨大的財(cái)產(chǎn)損失，同時(shí)也給機(jī)上工作人員和乘客的生命安全造成了巨大威脅［1］。防范鳥擊事件的前提是對(duì)鳥類目標(biāo)進(jìn)行實(shí)時(shí)探測(cè)跟蹤，雷達(dá)探鳥利用自身輻射的電磁波來(lái)探測(cè)鳥類，具有自動(dòng)化程度高、連續(xù)工作時(shí)間長(zhǎng)、數(shù)據(jù)便于分析存儲(chǔ)、受天氣和光線條件影響小等諸多優(yōu)點(diǎn)，故使用雷達(dá)來(lái)探測(cè)跟蹤鳥類目標(biāo)具有天然的優(yōu)勢(shì)［2～4］。

鳥類目標(biāo)具有飛行高度低、目標(biāo)小、速度慢等特點(diǎn)，其雷達(dá)回波信號(hào)弱，容易被雜波和噪聲信號(hào)掩蓋，目標(biāo)出現(xiàn)漏檢的概率較大；同時(shí)，鳥類經(jīng)常聚群飛行，數(shù)量多，機(jī)動(dòng)性強(qiáng)，空間運(yùn)動(dòng)軌跡復(fù)雜，航跡交叉與分離頻繁，雷達(dá)航跡管理難度大［3］。因此，在使用雷達(dá)探測(cè)跟蹤鳥類目標(biāo)時(shí)，對(duì)跟蹤算法要求非常高，而在大量目標(biāo)跟蹤算法中，多假設(shè)跟蹤（MHT）算法是能夠滿足此類目標(biāo)探測(cè)跟蹤需求的算法之一［4～5］。

鳥類翅膀的撲動(dòng)會(huì)產(chǎn)生微多普勒效應(yīng)，且不同的鳥，在不同的姿態(tài)和雷達(dá)視角下產(chǎn)生的微多普勒特征也不相同［6～7］，因而可利用鳥類產(chǎn)生的微多普勒特征進(jìn)行鳥類目標(biāo)探測(cè)跟蹤。本文將在MHT算法中引入微多普勒特征參量，以提高M(jìn)HT算法對(duì)鳥類目標(biāo)探測(cè)跟蹤的準(zhǔn)確性。

2 利用鳥類目標(biāo)微多普勒特征進(jìn)行參數(shù)估計(jì)

建立鳥類目標(biāo)運(yùn)動(dòng)模型如圖1所示［6］。忽略鳥身的寬度，以翅膀與鳥身的連接處為坐標(biāo)原點(diǎn)O，鳥飛行方向?yàn)閄軸正方向，建立如圖所示目標(biāo)坐標(biāo)系（X，Y，Z）。為簡(jiǎn)化分析，假設(shè)翅膀僅在YOZ平面內(nèi)作撲翼運(yùn)動(dòng)，且撲翼過(guò)程中翅膀不彎折，撲翼頻率為F0，半翼展為L(zhǎng)，撲翼角度ψ（即翼面與XOY平面的夾角）為隨時(shí)間變化的函數(shù)：

式中，A（0＜A＜π）為撲翼幅度，ω0=2πF0表示撲翼角頻率，ψ0表示初始時(shí)刻撲翼角度的相位，意義不大，為簡(jiǎn)化分析，令ψ0=0。雷達(dá)坐標(biāo)系（U，V，W）與目標(biāo)坐標(biāo)系平行，雷達(dá)位于原點(diǎn)S，到O點(diǎn)的距離為R0。雷達(dá)視線方向在USV平面上的投影與U軸的夾角為雷達(dá)方位角α，與USV平面的夾角為雷達(dá)俯仰角β。

圖1 鳥類目標(biāo)運(yùn)動(dòng)模型與雷達(dá)的位置

通過(guò)推導(dǎo)，鳥整體的回波信號(hào)s（t）為鳥身和左右翅膀回波信號(hào)之和，s（t）的微多普勒特征為鳥身回波的頻移fd和左右翅膀回波微多普勒特征fmD1（t）和fmD2（t）的總和，各分量表達(dá)式如式（2）～（4）：

利用微多普勒特征進(jìn)行參數(shù)估計(jì)的步驟如下：

1）計(jì)算回波信號(hào)的自相關(guān)函數(shù)，估計(jì)目標(biāo)撲翼頻率ω0。

可知散射點(diǎn)P′和Q′到雷達(dá)的距離由初始距離、平動(dòng)和微動(dòng)三部分構(gòu)成，其中平動(dòng)速度根據(jù)回波頻譜容易進(jìn)行估計(jì)，微動(dòng)部分是以T0=2π/ω0為周期的周期函數(shù)［6，8］。由復(fù)合函數(shù)的性質(zhì)可知，平動(dòng)補(bǔ)償后的目標(biāo)回波s（t）也是周期為T0的時(shí)間序列?；夭ㄐ盘?hào)的自相關(guān)函數(shù)和平均幅度差函數(shù)均可用于估計(jì)信號(hào)周期，假設(shè)平動(dòng)補(bǔ)償后的雷達(dá)回波為長(zhǎng)度為N的周期序列x（n）（n=1，2，…，N），x（n）的周期延拓為，其自相關(guān)函數(shù)Ф1（k）和平均幅度差函數(shù)Ф2（k）為

當(dāng)信號(hào)長(zhǎng)度超過(guò)一個(gè)周期時(shí)，其自相關(guān)函數(shù)Ф1（k）將在0，±T0，±2T0，…處出現(xiàn)極大值，平均幅度差函數(shù)Ф2（k）將0，T0，2T0，…在處出現(xiàn)極小值。目標(biāo)的撲翼頻率和回波信號(hào)的頻率相等，因此可以通過(guò)估計(jì)回波信號(hào)的周期得到撲翼頻率的估計(jì)值。

2）對(duì)目標(biāo)回波信號(hào)進(jìn)行時(shí)頻分析，提取目標(biāo)微多普勒頻率。

3）對(duì)目標(biāo)微多普勒頻率進(jìn)行快速傅里葉變換，得到各諧波分量的幅度絕對(duì)值c1～c4。

4）根據(jù)雷達(dá)視角，選擇使用式（6）或者式（7）求得到撲翼幅度估計(jì)值。

5）根據(jù)式（3），令

則式（3）可以重寫為

將撲翼幅度的估計(jì)值A(chǔ)和撲翼頻率的估計(jì)值ω0代入式（9），可以得到g1（t）的估計(jì)值。定義微多普勒譜寬BmD1和函數(shù)g1（t）的譜寬Bg1為

由式（10）可知：

式（11）中，BmD1由目標(biāo)回波獲取，撲翼幅度的估計(jì)值A(chǔ)和撲翼頻率的估計(jì)值ω0已知，再代入其他參數(shù)，即可計(jì)算得到半翼展的估計(jì)值。

同理，可由式（11）計(jì)算得到半翼展的估計(jì)值，因此，鳥類目標(biāo)半翼展的估計(jì)值確為和的均值：

3 面向航跡的MHT算法

面向航跡的MHT算法主要分為：1）對(duì)現(xiàn)有的航跡進(jìn)行預(yù)測(cè)，將新的量測(cè)和預(yù)測(cè)波門進(jìn)行關(guān)聯(lián)，生成假設(shè)得到新的航跡；2）計(jì)算新航跡得分，并基于得分進(jìn)行航跡刪除和確認(rèn)；3）對(duì)航跡刪除和確認(rèn)后幸存航跡進(jìn)行聚類；4）在每個(gè)聚類中獨(dú)立生成全局假設(shè)，并尋找最優(yōu)全局假設(shè)；5）基于最優(yōu)全局假設(shè)進(jìn)行N掃描回溯剪枝，并輸出最優(yōu)航跡［9］。

3.1 航跡生成

如圖2所示，假定在某個(gè)區(qū)域1～k時(shí)刻收到了多批量測(cè)。k-1時(shí)刻，已經(jīng)生成了多個(gè)航跡，可能存在的航跡如圖中實(shí)線，并且通過(guò)卡爾曼預(yù)測(cè)，得到了量測(cè)②的跟蹤波門和量測(cè)③的跟蹤波門。在k時(shí)刻，收到了五個(gè)量測(cè)，其中量測(cè)①、②和④落入波門中，量測(cè)①和③落入波門中，量測(cè)⑤沒有落入任何波門（如圖3所示）。因此，將k時(shí)刻的量測(cè)和k-1時(shí)刻的航跡關(guān)聯(lián)，可得到如圖2中虛線所示假設(shè)：1）k時(shí)刻五個(gè)量測(cè)都是新目標(biāo)；2）k時(shí)刻五個(gè)量測(cè)都是虛警；3）k時(shí)刻量測(cè)①、②和④是k-1時(shí)刻量測(cè)③的繼續(xù)，k時(shí)刻量測(cè)①和③是k-1時(shí)刻量測(cè)②的繼續(xù)，其中，k-1時(shí)刻量測(cè)②可能是k-2時(shí)刻兩個(gè)量測(cè)的繼續(xù)，因此從該量測(cè)到k時(shí)刻量測(cè)①和③都分別有兩個(gè)假設(shè)。同時(shí)，可以得到1～k時(shí)刻的航跡樹，如圖4所示。

圖2 k-1時(shí)刻航跡

圖3 k時(shí)刻量測(cè)關(guān)聯(lián)

圖4 k時(shí)刻航跡樹

3.2 航跡得分計(jì)算

航跡得分由遞歸累積產(chǎn)生。每一個(gè)航跡的得分等于它的上一次的值加上一個(gè)得分增量Δlk，即：

其中，Zk是目標(biāo)測(cè)量值，Xk|k-1是目標(biāo)預(yù)測(cè)位置，S是殘差協(xié)方差矩陣；Pd為目標(biāo)的探測(cè)概率；βf為虛警的空間密度；Pf為虛警概率；M為量測(cè)的維數(shù)。且航跡初始得分為l1=ln[(1-βf)/βf]。

通過(guò)航跡得分，可以進(jìn)行航跡的刪除和確認(rèn)，即lk≤ TL，則刪除航跡；lk≥ TU，則確認(rèn)航跡；TL≤ lk≤TU，等待更多的數(shù)據(jù)更新航跡，其中TL定義為航跡刪除閾值值，TU為航跡確認(rèn)閾值。這樣很多得分低的航跡在此過(guò)程中被刪除，大大減少了后面生成假設(shè)的數(shù)量［10～11］。

3.3 聚類

由于某些目標(biāo)在時(shí)空中的可區(qū)分性很高，因此可將所有的假設(shè)航跡劃分成子集，在每個(gè)子集中獨(dú)立地進(jìn)行假設(shè)生成、全局級(jí)航跡刪減和航跡更新等操作，即可將大規(guī)模跟蹤問(wèn)題化解為若干個(gè)小規(guī)模跟蹤問(wèn)題來(lái)處理，降低了算法復(fù)雜度和計(jì)算量，假設(shè)航跡子集即為航跡類。

圖5 k時(shí)刻聚類及全局假設(shè)生成

在新掃描周期內(nèi)分別將新接收到的量測(cè)與以前的類進(jìn)行互聯(lián)，若上一周期處理中兩個(gè)獨(dú)立的類與同一個(gè)量測(cè)相關(guān)，則這兩個(gè)類合并成一個(gè)類，不與任何舊聚類相關(guān)的量測(cè)形成新的類。將圖4中的航跡樹進(jìn)行聚類，由于航跡1和航跡2有共同量測(cè)，必須合并成一個(gè)類，航跡3為一個(gè)新類，最終得到兩個(gè)類，如圖5所示。

3.4 全局假設(shè)生成

在每個(gè)類中，利用圖論算法，將假設(shè)航跡映射為圖論中圖的頂點(diǎn)，航跡共享量測(cè)則對(duì)應(yīng)頂點(diǎn)是相鄰的（即頂點(diǎn)間存在邊）。全局假設(shè)生成即為頂點(diǎn)集的劃分操作，根據(jù)聚類原則，相鄰頂點(diǎn)必在同一集合中（如圖5），而尋找全局最優(yōu)假設(shè)則為尋求極大連通子集的過(guò)程，此過(guò)程可由深度優(yōu)先搜索算法實(shí)現(xiàn)。圖5中，深色部分即為全局最優(yōu)假設(shè)［12］。

3.5 N-scan剪枝

N-scan剪枝法通過(guò)限制航跡樹深度來(lái)控制假設(shè)數(shù)量，強(qiáng)制在k-N時(shí)刻產(chǎn)生的不確定性在k時(shí)刻延遲解決。當(dāng)軌跡樹的深度大于N時(shí)，N-scan剪枝法將保留全局最優(yōu)假設(shè)中的根節(jié)點(diǎn)分枝，刪除其余分枝［10］。如圖6所示，N=2時(shí)在k時(shí)刻航跡樹執(zhí)行N-scan剪枝操作，上一步已得到最優(yōu)假設(shè)為132、223、005，保留該分枝，而刪除其他分枝，最終得到k時(shí)刻的航跡（如圖7所示）。

圖6 k時(shí)刻剪枝

圖7 k時(shí)刻航跡

4 引入微多普勒特征的MHT算法

考慮在MHT算法的第二步，即計(jì)算航跡得分時(shí)引入微多普勒頻率。利用雷達(dá)歷史回波得到目標(biāo)的參數(shù)估計(jì)，根據(jù)MHT算法中與目標(biāo)相關(guān)聯(lián)量測(cè)的坐標(biāo)得到目標(biāo)的雷達(dá)視角，將目標(biāo)的參數(shù)估計(jì)和雷達(dá)視角代入式（2）～（4），計(jì)算目標(biāo)的多普勒特征理論值ft。根據(jù)當(dāng)前回波提取目標(biāo)多普勒特征觀測(cè)值fa，并與理論值ft做對(duì)比，將兩者差值的均方根作為目標(biāo)多普勒特征得分lf。給式（13）中MHT算法航跡得分lk分配權(quán)重ak，目標(biāo)多普勒特征得分lf分配權(quán)重af，則假設(shè)航跡最終得分為

其中，權(quán)重ak、af與雷達(dá)性能、目標(biāo)特征、雷達(dá)工作環(huán)境有關(guān)，需經(jīng)過(guò)多次實(shí)驗(yàn)確定。

5 仿真計(jì)算

仿真中，模擬產(chǎn)生兩只不同種類的鳥目標(biāo)，其半翼展、撲翼幅度和撲翼頻率均不相同，目標(biāo)1做直線運(yùn)動(dòng)，目標(biāo)2做曲線運(yùn)動(dòng)，兩個(gè)目標(biāo)在空間交匯兩次（如圖8所示）。雜波點(diǎn)個(gè)數(shù)服從泊松分別，雜波密度3×10-5，N-scan剪枝深度為3，量測(cè)噪聲服從高斯分布，標(biāo)準(zhǔn)差為100m。分別用本文算法和傳統(tǒng)MHT算法進(jìn)行100次Monte Carlo仿真計(jì)算，得到兩種算法失跟率如表1所示，兩種算法對(duì)兩個(gè)目標(biāo)跟蹤的距離均方根誤差曲線分別如圖9、10所示。

圖8 目標(biāo)軌跡

表1 算法失跟率統(tǒng)計(jì)表

圖9 目標(biāo)1跟蹤誤差RMSE曲線

從表1可以看出，在目標(biāo)跟蹤穩(wěn)定跟蹤率上，本文算法要優(yōu)于傳統(tǒng)MHT算法。而從圖9、10可以看出，在位置的跟蹤精度上，兩種算法差距不大。另外，將圖8放大，可以看到在兩個(gè)目標(biāo)交匯處，點(diǎn)跡距離很近，航跡關(guān)聯(lián)難度更大，對(duì)算法的性能提出了挑戰(zhàn)。通過(guò)統(tǒng)計(jì)，傳統(tǒng)MHT算法將交匯處點(diǎn)跡錯(cuò)誤關(guān)聯(lián)的概率高達(dá)30%，而本文算法錯(cuò)誤關(guān)聯(lián)的概率僅1%。因此，在航跡精確關(guān)聯(lián)方面，本文算法性能要優(yōu)于傳統(tǒng)MHT算法。

圖10 目標(biāo)2跟蹤誤差RMSE曲線

圖11 目標(biāo)軌跡局部圖

6 結(jié)語(yǔ)

本文介紹了根據(jù)微多普勒特征進(jìn)行鳥類目標(biāo)參數(shù)估計(jì)的方法和面向?qū)ο蟮腗HT算法基本流程，提出利用參數(shù)估計(jì)的結(jié)果計(jì)算鳥類目標(biāo)微多普勒特征得分，并將該得分引入MHT算法，作為MHT算法中航跡得分的組成部分，從而實(shí)現(xiàn)了考慮微多普勒特征的鳥類目標(biāo)MHT算法。通過(guò)仿真計(jì)算表明，該算法與傳統(tǒng)MHT算法相比，失跟率更低，且在處理交叉航跡時(shí)，航跡關(guān)聯(lián)精度更高。