基于多普勒與微多普勒聯合利用的弱小目標檢測與估計方法

宋志勇 許云濤

(國防科技大學電子科學學院 長沙 410073)

1 引言

近年來，無人機的運用越來越廣泛，在遙感測繪、搶險救援和打擊犯罪等領域中發揮了重要的作用，然而，在許多有益應用前景的反面，無人機在機場凈空區內非法升空等一系列威脅公共安全的事件時有發生，隨著以無人機為代表的低慢小飛行器的應用快速增加，其帶來的威脅和風險也逐漸顯露。這類目標的飛行高度低、運動速度較慢、雷達反射截面積(Radar Cross Section, RCS)小，導致其雷達回波信噪比(Signal Noise Ratio, SNR)低、且多普勒頻率較低，很難將目標和背景區分開來。無人機類目標的特點是同時存在整體徑向運動和內部微動，微動是指物體的某些部分相對于物體質心的平動存在小幅振動、旋轉和其他高階運動的現象[1–4]，在頻域表現為微多普勒特征。微多普勒現象是Chen[1]發現的一種多普勒譜調制現象，由于目標內部部件存在振動或旋轉等運動，而其部件始終與目標主體連接，導致它們在雷達視角下呈現出速度的周期變化，所以微多普勒特征會表現出交叉、重疊等周期性現象[5]。與徑向運動帶來的多普勒特征不同，微多普勒特征是對多普勒頻率的高階調制。隨著雷達分辨率的提高，聯合利用多普勒和微多普勒特征，提升低慢小目標的檢測性能成為當前目標檢測一種可能的技術方向[6]。

現有對于微多普勒參數估計的研究主要分為兩種方法[5]，將回波信號經過時頻變換后，一種方法借鑒光學圖像處理中的提取曲線的經典思路，在時頻圖中搜索局部峰值點得到頻率脊[7–10]，使用一定方法將所有時刻的頻率脊做排列，即可得到微動目標各部件的運動速度變化曲線，從而反演出微多普勒參數；另一種方法則使用逆拉登變換[11]、Hough變換[12,13]等方法將時頻圖轉換到某一參數空間中，通過在參數空間檢測峰值點，來估計微多普勒參數。使用微多普勒信息進行目標檢測的研究還比較少。

將多普勒和微多普勒統一到1個數學框架下，是有效利用二者能量進行檢測的基礎。隨機有限集(Random Finite Sets, RFS)濾波器是一種序貫貝葉斯的狀態估計方法，由Mahler[14]于21世紀初提出，Vo等人[15–21]基于Mahler的工作，提出了一系列適用于多目標跟蹤的濾波器，其觀測模型主要分為兩類[22]，即檢測器輸出測量模型[17]和強度測量模型[23]，前者需要將回波先經過一個低門限檢測器，這種方法會造成一些信息丟失，后者將回波直接輸入濾波器，基于目標的運動模型和目標回波的觀測模型，可以完整使用目標回波的所有信息。勢均衡多目標多伯努利(Cardinality Balanced Multi-target Multi-Bernoulli, CBMeMBer)濾波器[23]是隨機有限集濾波器中的一種，可以在檢測目標的同時獲取目標狀態的估計，其具備良好的多目標檢測性能，同時計算量較小。

針對現有低慢小目標檢測只使用目標徑向運動產生的多普勒信息，丟失了一部分由于目標微動而分散的能量，導致目標檢測性能下降的問題，本文將微多普勒參數估計問題與微動目標檢測問題結合起來，對多普勒和微多普勒進行聯合建模，提出一種基于CBMeMBer濾波器的無人機目標聯合檢測與估計方法。所提算法在傳統點測量目標模型的基礎上對多普勒和微多普勒進行聯合建模，利用微多普勒帶來的信息增量，提高了濾波器的檢測靈敏度，同時可以獲得目標位置以及微動狀態的估計，具有適應多無人機目標檢測與估計的潛力。仿真結果表明，本文所提算法充分利用兩類信息蘊含的特征，有效實現了目標信號能量的積累，檢測性能相較傳統方法有所提升。

2 目標多普勒與微多普勒聯合建模

假設一無人機目標主體存在低速運動，內部微動部件存在旋轉運動，場景如圖1所示，可以對目標雷達回波進行建模。

圖1 實驗場景

2.1 目標多普勒模型

假設目標以恒定速度遠離雷達運動，其回波會產生多普勒頻移。雷達發射信號可表示為其中，t0表示目標在慢時間維的回波時延，ts為脈沖間慢時間，λ為雷達載波波長，v為目標相對雷達的運動速度，遠離雷達時為正。

2.2 目標微多普勒模型

假設目標內部存在一個雙翼螺旋槳，槳葉關于旋轉中心呈軸對稱，其旋轉軸方向與雷達視線方向垂直，槳葉初始方向與目標運動速度方向一致。旋轉目標滿足理想散射點模型[3]，一片槳葉可由一個位于槳葉尖端的散射中心等效，將位于槳葉上的散射點投影到雷達平面如圖2所示。

圖2 雷達-旋轉散射點關系示意圖

圖2中，O表示近似的雷達波束中心點，O′表示目標的旋轉中心，D為OO′的長度，即雷達與目標之間的距離。目標的旋轉部件散射點為P，旋轉半徑為r，旋轉頻率為fr，散射點運動初始相位為θ0，t時刻的相位為θt=2πfrt+θ0，散射點P到雷達的距離為

2.3 目標多普勒與微多普勒聯合模型

式(5)中的求導將式(4)中的D作為常量忽略了，如果目標存在運動速度v，那么式(3)可以寫為

假設目標的一片槳葉可使用一個散射點近似，根據式(8)，單個槳葉的回波相位為

3 基于CBMeMBer濾波器的微弱目標聯合檢測估計

基于前述分析，對式(11)中回波慢時間維進行傅里葉變換

3.1 目標運動模型

基于強度測量模型，假設目標可由式(14)中向量表示

其中，下標k表示當前時刻，rk表示目標與雷達的距離，vk表示目標與雷達的相對速度(遠離為正)，ρk表示目標的平均RCS系數，Vk表示最大相對速度差，ωk=2πfr,k表示目標螺旋槳的旋轉角速度，θk表示微多普勒調制相位，Bk表示目標螺旋槳的槳葉個數。

向量中的Vk表示距離-多普勒圖中目標多普勒最大展寬所對應的最大速度變化值，而并非實際物理意義對應的旋轉半徑rk，這是因為目標在距離-多普勒圖中的多普勒展寬并非由rk一個變量決定，其與旋轉頻率fr,k也關聯，而旋轉頻率又與目標距離維的運動周期相關聯，選用rk作為濾波器狀態輸入容易造成濾波器無法收斂，故這里選擇了Vk和ωk作為控制微多普勒頻率調制的關鍵參數。

目標狀態轉移函數為

其中

圖3展示了目標在距離-多普勒圖中的測量點與其可能的運動軌跡，其橫軸與縱軸分別為距離軸與頻率軸。t=0 時框選的1個“ × ”與t=1/4fr時框選的3個“ ×”分別表示這兩個時刻目標回波變換得到的譜峰。其中，t=0時，目標的兩個槳葉旋轉到雷達視線所在直線上，此時，槳葉的速度方向與雷達視線垂直，其速度在雷達視角為 0，所以呈現出與目標主體相同的運動速度，假設該時刻微多普勒調制相位θ=0。t=1/4fr時，目標的兩個槳葉旋轉到雷達視線的垂直線上，此時一個槳葉向遠離雷達方向運動，一個槳葉向靠近雷達方向運動，所以其速度呈現出關于目標主體速度對稱的特征，該時刻θ=π/2。

圖3 微多普勒譜與狀態參數關系示意圖

各個槳葉雖然存在與主體不同的速度，但其仍然是目標上的一部分，所以其速度會在t=1/4fr之后漸漸又接近主體速度，并呈現周期變化，其最大頻譜展寬為 ?f=-4V/λ。假設目標運動速度不變，則目標螺旋槳旋轉一個周期時目標相對雷達的運動距離為 ?r=v/fr。

假設距離-多普勒圖的分辨率為Nr×Nv，其距離維分辨率為Nr，速度維分辨率為Nv，k時刻測量數據可表示為Yk={y1,k,...,ym,k}Nr×Nv，一個狀態矢量xk會影響集合T(xk)中所包含的所有測量，即yl,k,l ∈T(xk)，yl,k表示一個像素，l為該像素的絕對索引，T表示狀態空間向觀測空間點的映射函數。像素絕對索引與該像素所對應的距離和速度之間存在映射L(r,v)=l，其逆映射為L-1(l)=r,v。

本文使用點測量假設，一個狀態對應B+1個測量點，分為目標主體的一個測量點，以及目標螺旋槳產生的B個測量點，各散射點的強度會因為雷達視線角不同而變化，數量上一般不會變化。

假設狀態xk共對應Bk+1個測量點，對于某個測量點集合ti ?T(xk),i ∈[1,Bk+1] ， 其距離為rk，速度為

一個狀態所影響的集合T(xk) 為所有Bk+1個測量點子集的并集，即

3.2 目標觀測模型

3.3 目標測量似然函數

對于某個狀態xk，觀測似然比[23]可表示為

3.4 目標聯合檢測估計濾波器構建

假設一個定義在狀態空間χ上的伯努利RFS為空的概率為 1-r，存在單個狀態x的概率為r，x服從概率密度函數p所定義的分布。

一個定義在χ上的多伯努利RFS是固定個數為M的獨立伯努利RFSX(i)的集合。其概率生成泛函(Probability Generating Functional, PGFl)定義為

其中，r(i)和p(i)為X(i)的單項存在概率和狀態的概率密度函數。

則多伯努利RFS的概率密度可表示為

根據文獻[23]的推導，CBMeMBer濾波可分為預測和更新兩個過程，其序貫蒙特卡羅(Sequential Monte Carlo, SMC)實現可用算法1表示，其中，所有帶有k或k-1 下標的變量分別表示k時刻或者k-1時刻的該變量值。

在使用算法1對測量進行濾波后，需要計算當前時刻的估計，可以使用

4 仿真實驗

為了驗證本文所提算法的有效性，仿真使用高重頻脈沖多普勒(Pulse Doppler, PD)雷達對遠距離微動目標進行檢測估計。為了簡化計算量，仿真以PD雷達距離-多普勒圖作為算法輸入，其參數設置如表1所示。總共40個時刻，目標自0時刻出現，一直存在至最后一個時刻，其部分參數如表2所示。

表1 仿真雷達參數

表2 仿真目標參數

仿真場景如圖4所示，雷達靜止不動，目標主體存在徑向運動，同時目標內部存在旋轉運動。

圖4 仿真場景示意圖

根據仿真設置，可以生成各時刻目標的距離-多普勒圖，第8幀與第16幀如圖5所示。從圖5可以清晰地看到由目標所產生的3個多普勒譜峰，這表明在同一距離上出現了3個不同速度的散射點。對比圖5(a)和圖5(b)兩幅子圖可以觀察到，3個峰值點中，有1個峰值點的速度始終不變，兩個峰值點的速度變化呈現出周期性調制的特性。由此可以得到與式(13)一致的結論，即目標螺旋槳的回波呈現出周期性的多普勒頻率調制，且當其旋轉軸與雷達視線垂直時，該調制函數為正弦函數。

圖5 第8幀與第16幀距離-多普勒圖

采用本文的方法對無人機目標進行檢測估計，濾波得到的目標狀態估計結果如圖6所示，從中可以看出，濾波器可以準確估計出各時刻目標的距離、速度信息，并且能夠同時準確估計目標的微多普勒信息。

圖6 2dB信噪比下濾波器估計結果與真值

為了驗證本文算法在不同信噪比下的檢測性能，可以調整ks使目標在相同距離下的信噪比變化，在不同信噪比下進行100次蒙特卡羅實驗，統計目標存在情況，可以得到目標檢測概率隨信噪比變化曲線，如圖7中圈號實線所示。可以看到，隨著信噪比增加，濾波器對目標的檢測概率逐漸趨近于1。

圖7 兩種算法的檢測概率隨信噪比變化曲線

為了對比本文方法的性能，同時給出僅使用多普勒特征進行目標檢測估計的檢測性能，如圖7中加號虛線所示。通過對比可以看到，由于微多普勒特征的加入，算法可以使用的信息增加了，在回波信噪比為1 dB以下時的檢測概率相較僅使用多普勒信息的情況下更高，在檢測概率為0.92時，檢測靈敏度提高了2 dB左右。

為了分析各信噪比下濾波器的估計性能，可以計算濾波器估計結果與真值之間的最優子模式分配(Optimal Sub Pattern Assignment, OSPA)距離[14]，對所有時刻的OSPA距離取平均，可以得到如圖8所示平均OSPA距離隨信噪比變化曲線。可以看到，隨著信噪比增加，平均OSPA距離不斷減小，在0.5 dB開始趨于穩定，這說明算法在該信噪比下趨于收斂。

圖8 本文算法的平均OSPA誤差隨信噪比變化曲線

通過改變假設的槳葉個數B，可以得到如圖9所示不同槳葉個數下檢測概率隨信噪比變化曲線，可以觀察到，隨著槳葉數量的增加，在低信噪比下的目標檢測概率整體上呈現增加的趨勢，這說明槳葉個數的增加會帶來額外的信息增量。

圖10所示為不同旋轉速度ω下檢測概率隨信噪比變化曲線，可以看到，雖然旋轉速度ω有所變化，但是各信噪比下的檢測概率整體上呈現出近似的變化。這說明旋轉速度在一定范圍內變化時，不會對檢測概率帶來非常大的影響。

為了對比本文算法的性能，嘗試將粒子濾波器中粒子數分別設置為700, 1 000, 1 500，運行100幀，計算每一幀平均消耗時間，結果如表3所示。在不同粒子數設置情況下，對比不同信噪比下算法檢測概率，如圖11所示。計算機使用Windows 10操作系統，CPU為AMD R4750U，內存為48 GB，運行MATLAB R2022b實驗軟件。

結合表3與圖11，可以看到，本文算法在3種粒子數設置下的檢測概率差別不大，在檢測概率0.9時，檢測靈敏度提高了2 dB左右。在粒子數設置為700時，計算性能優于僅使用多普勒信息的情況。

5 結束語

本文提出聯合多普勒和微多普勒的目標檢測前跟蹤方法，通過對微動目標所產生的多普勒和微多普勒現象進行分別分析，建立了兩者的聯合信號模型，通過選擇合適的參數，計算了濾波器所需要的運動方程，并分析了由目標狀態得到目標測量的觀測模型；對本文所提濾波器進行了100次蒙特卡羅實驗，給出了檢測概率性能曲線，并對濾波結果做了誤差分析，對比了在不同槳葉數、旋轉速度下濾波器的檢測性能，驗證了濾波器的估計準確性。

與僅使用多普勒信息的傳統算法相比，本文算法得益于充分利用了多普勒和微多普勒兩個特征，在檢測概率為0.9時的檢測靈敏度提高了2 dB左右，算法在信噪比為0.5 dB時能夠穩定實現目標檢測，同時OSPA誤差趨近于0.1，估計精度較高。后續工作可以根據實測數據驗證本文算法的有效性，針對實際可能遇到的回波信號強閃爍問題，修改觀測模型。