雷達鳥類目標微多普勒貝葉斯增強算法

楊 磊, 孫衛天, 毛欣瑤, 夏亞波, 桑婧隺

(1. 中國民航大學電子信息與自動化學院, 天津 300300;2. 天津港遠航國際礦石碼頭有限公司, 天津 300452)

0 引 言

近年來,鳥類入侵事件頻繁在機場凈空保護區發生,如2009年全美航空一架代號為1549的航班在起飛不久后就撞上一群加拿大黑雁,使得兩個引擎都失去了推力而不得不滑翔飛入哈德森河;2020年5月9日,CA908航班因為其發動機葉片遭遇鳥擊而被迫返航。類似鳥擊事件不斷上演,給飛行安全帶來了極大威脅。因此,實現機場附近鳥群可靠探測與監視對避免鳥群威脅航空飛行安全具有十分重要的意義[1]。

雷達由于具有全天候、全天時、遠距離探測等優勢,被廣泛運用于探測鳥類等非合作目標場合[2]。目前相對成熟的探鳥雷達系統有美國DeTect公司開發的MERLIN(梅蘭)系列、加拿大Sicom System公司開發的Accipiter(蒼鷹)系列以及荷蘭TNO公司開發的ROBIN(羅賓)系列[3]。MERLIN系列天線掃描區域廣,效率高,但無法獲得三維信息;Accipiter系列能夠獲得鳥類目標的三維信息,但其天線掃描區域有限,效率較低;ROBIN系列具有低成本、高精度、能夠對信號和圖像進行一體化處理、可同時跟蹤上萬鳥類的優點,但是其在海雜波濾除等方面的測試尚在進行中。鳥類目標的雷達回波強度受限于兩個因素:一是鳥類目標具有比較小的雷達散射截面(radar cross section, RCS),二是機場近場強雜波的影響。因此,利用傳統恒虛警率(constant false-alarm rate, CFAR)方法進行鳥類目標探測會比較困難。此外,鳥類目標是非合作目標,其機動性較強,難以實現較高的積累響應,因此有必要研究針對鳥類目標的有效探測方法。

對于鳥類目標,其飛行時翅膀規律性的撲翼運動會在由主體平動產生的雷達多普勒頻移回波信號附近引入額外調制邊帶,即微多普勒效應。而微多普勒效應反映了多普勒頻移瞬時特性,表征了目標微動的瞬時徑向速度。同時,利用微多普勒信號可以反演出目標運動特性,例如翼型、姿態等信息。因此,盡管鳥類目標雷達回波強度較弱,但可以利用其中微多普勒效應對其進行有效的探測和識別[4-5]。

由于鳥類目標微多普勒效應隨時間變化,而傳統傅里葉變換不能描述任一時刻的信號頻率成分,無法對其進行全面分析。為了解微多普勒效應隨時間變化的關系,需要使用時頻分析方法,將時域一維信號映射到時頻二維平面,全面反映非平穩信號時頻聯合特征。常見的時頻分析方法分為兩種,即非線性時頻分析方法和線性時頻分析方法[6]。其中,非線性時頻分析方法有維格納-威爾分布(Wigner-Ville distribution, WVD)、偽維格納-威爾分布(pseudo Wigner-Ville distribution, PWVD)、平滑偽維格納-威爾分布(smoothed pseudo Wigner-Ville distribution, SPWVD)、希爾伯特黃變換(Hilbert-Huang transform, HHT)。WVD將信號能量分布于時頻二維平面內,其時間帶寬積達到了測不準原理下界,但是由于WVD是雙線性變換,因此會產生交叉項,使得時頻分布受到干擾且抑制了二次型時頻分布的推廣。PWVD是對WVD進行加窗處理,使得WVD的完全非局部性變為局部化,在一定程度上抑制了多分量信號的交叉項。SPWVD是WVD與平滑函數進行卷積的結果,可以較好地消除交叉項,但其運算復雜。HHT是一種自適應方法,可通過信號經驗模態分解,使非平穩信號平穩化,但其存在邊界效應。無論是哪種非線性時頻分析方法,交叉項的產生或運算的高復雜度,均不利于后續信號處理和特征增強,因此本文采用線性時頻分析方法。

在線性時頻分析方法中,線性時頻分析方法有短時傅里葉變換(short-time Fourier transform, STFT)、小波變換(wavelet transform, WT)以及拉普拉斯變換(Laplace transform, S)。STFT需要在時域對信號進行加窗處理,這意味著STFT無法同時達到高時間分辨率與高頻率分辨率,即當窗函數變窄時,時間分辨率高而頻率分辨率低;當窗函數變寬時,時間分辨率低而頻率分辨率高。WT對時頻二維平面是一種機械式的劃分,且所獲結果不一定是真正的時頻譜。S變換相對于STFT運算較為復雜。因此有必要研究在時頻域上對微多普勒特征進行特征增強的方法。

由于鳥類目標微多普勒特征在時頻二維平面呈現出明顯的稀疏性,根據壓縮感知(compressed sensing, CS)理論可知,能夠利用特征增強的方法恢復微多普勒特征[7]。基于CS進行特征增強的方法主要包括貪婪類算法、凸優化類算法和貝葉斯類算法[8-9]。貪婪類算法通過迭代搜索均方根誤差最小的局部稀疏得到最優解,并在解的稀疏度小于某一經驗值時停止迭代[10-11]。常見的貪婪類算法有匹配追蹤(matching pursuit, MP)算法、正交匹配追蹤(orthogonal matching pursuit, OMP)算法等。雖然OMP算法運算效率快,但其稀疏度需預先自行設定,這使得稀疏度無法自適應地滿足實際目標特征的需求,難以適應現實先驗特征建模的需要,導致特征恢復精度有限。凸優化類算法通過最小化目標的1范數并結合基追蹤抑噪進行稀疏恢復[12-14],例如凸優化算法,其基于不斷優化設計最終得到最優解,但在高維數據情況下,凸優化算法運行性能低,耗費時間長。貝葉斯類算法通過引入目標特征先驗分布,利用分層模型進行推理,進而得到后驗概率密度函數估計[15-18],常見算法如稀疏貝葉斯學習(sparse Bayesian learning, SBL),其對目標特征以先驗概率表征,通過分層可以保證獲得解析解,求解精度高且無需多次獨立同分布實驗,僅需單次測量結果即可。通常,先驗分布選擇拉普拉斯分布、高斯分布。拉普拉斯先驗在降采樣、低信噪比情況下會造成求解過稀疏,即目標特征恢復結果損失一部分連續性。高斯先驗抑制噪聲能力不強,會造成目標特征恢復結果過于平滑。

針對以上問題,本文提出利用廣義高斯分布(generalized Gaussian distribution, GGD)在時頻域對鳥類目標微多普勒特征進行先驗自適應建模。該分布可以根據不同場景,靈活選擇先驗模型,從而避免了目標特征恢復結果因過于稀疏而丟失連續性、因過于平滑而保留過多噪聲的問題。然而,GGD先驗與所假設的高斯似然函數由于不共軛而無法得到閉合解。針對該問題,本文利用貝葉斯分層模型對鳥類目標微多普勒特征及相關參數進行層級建模推理,并建立參數依賴的后驗分布。考慮到動態可變GGD先驗造成目標特征后驗分布非平滑、求解運算效率低等困難,本文在可以快速高效求解后驗分布的未調整朗之萬算法(unadjusted Langevin algorithm, ULA)基礎上引入了近端算子,以高效求解目標特征非平滑后驗分布的解析解。最后,通過模擬強雜波噪聲等現實場景,將本文所提貝葉斯增強算法與STFT算法、凸優化算法特征恢復結果作對比,分析了本文所提貝葉斯增強算法在微多普勒特征增強方面的優勢,并利用實測數據驗證了該貝葉斯增強算法的有效性。

1 微動目標信號模型

自然界中,鳥類的飛行模式主要有兩種,一種是鳥類達到一定高度或速度后,伸開雙翅不動在空中滑行前進。此飛行模式稱為滑翔,多見于高空。另一種是撲翼飛行,這是所有鳥類最常見、最基本的飛行模式。鳥類在近地面以及低空飛行時大多采用該飛行模式。此時,鳥類飛行高度正好與起飛或者降落的飛行器所處高度重合,常常容易引發鳥擊事件。因此,研究鳥類低空撲翼飛行模式對于避免鳥類撞擊飛行器是十分有必要的。鑒于鳥類飛行時翅膀規律性撲翼運動是產生微多普勒效應的主要原因,下面將研究鳥類撲翼時的微動信號模型。

1.1 鳥類目標撲翼模型

鳥類翅膀結構如圖1(a)所示,若類比并借用人體手臂術語,鳥類翅膀包括了翼弦、半翼展、上臂、前臂以及手,能夠完成撲翼、掃翼及扭翼等運動形式。撲翼是指上臂或翼弦圍繞某一個關節以一定撲翼角度做上下運動。掃翼是指肩部做前向延伸或后向回收運動。扭翼是指翅膀圍繞鳥自身主軸做旋轉運動,該運動可以引發翅膀前緣下降以及后緣上升。接下來,主要側重分析鳥類的撲翼運動。

圖1 鳥類翅膀結構及運動學簡化模型Fig.1 Structure and simplified kinematic model of bird wing

在以上假設基礎上,上臂撲打角度的調和時變函數為

φ1(t)=A1cos(2πfflapt)+φ10

(1)

前臂撲打角度的調和時變函數為

φ2(t)=A2cos(2πfflapt)+φ20

(2)

前臂扭轉角度的調和時變函數為

(3)

(4)

鳥翼型肘關節位置為

(5)

綜上,借用鳥類翅膀運動學簡化模型,便可產生鳥類目標撲翼模型的雷達回波信號。

1.2 鳥類目標微多普勒模型

忽略鳥身寬度且假設鳥類翅膀在運動過程中翼面不彎曲,建立如圖2所示的空間坐標系。具體地,以雷達S為原點建立坐標系(X,Y,Z)和以鳥兩側翅膀根部連線中點O(以下稱為對稱中心)為原點建立鳥類目標撲翼坐標系(x,y,z),其中y軸方向為鳥首方向。雷達坐標系與鳥類目標撲翼坐標系平行。雷達S到O點的距離為R0,方位角為α,俯仰角為?。為簡便計算,假定鳥類目標以撲翼頻率Fbird沿y軸正半軸飛行,其翅膀上任意點的移動可以認為是平飛和翅膀撲翼兩種運動的合成。若平飛速度為v,撲翼角度φ(即翼面與xoy平面的夾角)為隨時間變化的函數,即

φ(t)=θmaxsin(2πFbirdt+φ0)=θmaxsin(ω0t+φ0)

(6)

式中:θmax為撲翼角度幅值,即翼面與xoy平面的最大夾角;ω0=2πFbird為撲翼角頻率;φ0表示初始時刻鳥類目標撲翼角度相位。為簡化分析,令φ0=0。

圖2 雷達與鳥類目標幾何關系圖Fig.2 Geometric relationship diagram of radar and bird target

R(t)=|R0+OO′+O′Q|2=

{(R0sin ?cosα+lcos[φ(t)])2+(R0sin ?sinα+vt)2+

(7)

R(t)≈R0+vtsin ?sinα+

lcosαsin ?cos[φ(t)]+lcos ?sin[φ(t)]

(8)

(9)

式中:-4πR(t)/λ為t時刻Q點相位。將式(8)代入式(9)并對其求導,可以得到Q點回波多普勒頻率為

[cos ?cosφ-cosαsin ?sinφ]=fd+fmicro(t)

(10)

式中:第一項fd表示由平動引起的多普勒頻率;第二項fmicro(t)表示由撲翼運動引起的微多普勒頻率。若觀測初始雷達徑向對準該鳥類目標,即α=0,fmicro(t)可進一步化簡為

cos[θmaxsin(ω0t)+?-ω0t]}

(11)

從式(11)可以看出,鳥類目標微多普勒頻率至少以2π為周期進行周期性變化。

2 雷達回波信號時頻分析

第1.2節給出了鳥類目標微多普勒模型,可基于此進行微多普勒效應仿真模擬。由式(11)可知,鳥類目標微多普勒頻率是周期時變信號,若采用常規傅里葉變換對整個信號進行分析,只能得到信號在頻域上的特性,不能獲取信號在時域上的信息。因此,需要對鳥類目標微多普勒模型進行時頻分析。本文采用物理意義明確的STFT時頻分析方法,其原理是在時域內選取一個合適的窗函數,并對該窗函數進行采樣。同時,對雷達回波信號y(t)進行采樣,得到回波采樣序列y。接著,利用窗函數對該采樣序列進行截取,并假定截取序列是平穩的。通過不斷移動窗函數得到不同時刻截取序列的傅里葉變換,這些傅里葉變換的集合就是對整段信號進行時頻分析的結果[20]。本文選擇高斯窗作為STFT窗函數。

xi=FHyi

(12)

式中:yi∈CL×1;H∈CL×L是由L個窗函數采樣點構成的對角矩陣;F∈CN×L是傅里葉變換矩陣;xi∈CN×1是對應yi的待恢復微多普勒特征。進一步,可將上述公式改寫為

(13)

式中:A′∈CL×N是STFT矩陣FH的偽逆。在實際探測鳥類目標時,由于周圍環境存在大量噪聲,雷達回波往往易受噪聲及近場強雜波的影響。因此,將上式修正為

(14)

式中:ni∈CL×1表示噪聲以及近場強雜波。由于觀測噪聲的影響,若直接求解,可能會出現頻域分辨率下降等問題。由于微多普勒特征在時頻二維平面具有稀疏特性,利用CS理論可知,可以把提取微多普勒特征轉化為基于CS的特征增強問題[7-9]。

2017年，五建進行了區域化改革。區域化運行在促進“資源優化、責任落實”的同時，也盡可能地實現了“安居樂業、人文關懷”。“讓員工可以在相對穩定的區域內工作是我們的一大初衷。這樣可以對自己的職業進行有針對性地規劃，比如在哪里工作，在哪里安家，孩子在哪里上學，如何照顧老人等等，解決了干事創業的后顧之憂。”蔣德軍說。

為了方便解釋后續參數更新過程,將全部M段截取信號排列成一個列向量,則有

Y=AX+N

(15)

式中:Y∈CLM×1;A∈CLM×NM是由A′構成的分塊對角矩陣;X∈CNM×1;N∈CLM×1。

3 微多普勒特征貝葉斯學習

貝葉斯學習將關于未知參數的先驗信息與樣本數據相結合,基于貝葉斯原理得出未知參數的后驗概率密度函數,并選擇合適的貝葉斯推理方式求解,進而推斷未知參數。

3.1 廣義高斯先驗概率模型

假設待恢復的微多普勒特征的每個元素均獨立同分布。為了實現動態模擬散射場先驗統計特性,自由選擇先驗模型,本文采用靈活自適應的GGD[21-22],其表示為

(16)

圖3 一維GGD概率密度函數Fig.3 Probability density function of one-dimensional GGD

3.2 回波數據似然函數模型

通常,與模型相應的似然函數可寫為

f(Y|X,β)=CN(Y|AX,βI0)

(17)

式中:CN(·)代表復高斯分布;β表示噪聲方差;I0表示單位矩陣。

3.3 貝葉斯分層模型

由于鳥類目標微多普勒場景GGD先驗與式(17)高斯似然函數不共軛,無法求得后驗概率密度函數的閉合解析解[23]。針對此問題,采用貝葉斯分層模型對先驗分布和似然函數分層以保證得到閉合解。由于GGD等價于p范數,為了能夠更好地實現先驗自適應,減少預固定先驗影響,在貝葉斯分層模式下將形狀參數p建模為均勻分布,即

(18)

式中:L[0,2]表示在區間[0,2]上的指示函數,當p在此區間時,L[0,2]=1,否則,L[0,2]=0。

為了減少系統誤差,將尺度參數γ建模為無信息Jeffrey’s先驗

(19)

式中:當γ>0時,LR+=1,否則為0。

通常,噪聲服從加性高斯白噪聲(additive white Gaussian noise, AWGN)。此時,可將噪聲方差β建模為共軛逆Gamma分布

f(β|a,b)=TG(β|a,b)=

(20)

式中:TG(·)表示逆Gamma分布;a和b分別代表形狀參數和比例參數,調整其數值可以改變Gamma分布的特性,通常令a=10-3,b=10-3,此時可以確保先驗沒有信息。

由此,通過對形狀參數p和尺度參數γ以及噪聲方差β進行層級建模,構建出了似然函數和先驗分布之間的聯系,為下節求解后驗概率密度函數奠定了基礎。根據式(16)～式(20),可得X,p,γ,β的聯合后驗分布為

f(X,p,γ,β|Y)∝

f(Y|X,β)f(X|p,γ)×

f(p)f(γ)f(β|a,b)∝

CN(Y|AX,βI0)GGD(X|p,γ)×

(21)

至此,鳥類目標微多普勒特征貝葉斯分層模型已構建完成。根據上述分層信息,可以使用有向無環圖(directed acyclic graph, DAG)表示各個隨機變量之間的關系,如圖4所示。

圖4 概率圖模型Fig.4 Probabilistic graphical model

圖中,Y為顯隨機變量,表示雷達回波數據;X,p,γ,β為隱隨機變量,X表示待恢復的微多普勒特征,p,γ,β為貝葉斯分層建模引入的3個超參數,分別為GGD先驗的形狀參數、尺度參數以及目標場景內加性噪聲方差;a和b為固定常量參量。

3.4 后驗概率密度函數求解

求解后驗分布的方法可分為兩類:近似確定性算法和統計采樣算法[24-27]。其中,Gibbs采樣作為一種典型的統計采樣算法,基于馬爾可夫鏈蒙特卡羅(Markov chain Monte Carlo,MCMC)[24]方法建立,能夠在確定后驗概率分布的條件下,通過循環采樣獲取參數估計值。變分貝葉斯期望最大化(variational Bayesian expectation maximization,VB-EM)是典型的近似確定性算法[25-27],可以實現隱隨機變量推理,但其運算涉及矩陣求逆,求解復雜且效率低。為了降低運算復雜度,本文選擇統計采樣算法求解各個參數服從的條件后驗分布。

根據式(17)和式(20),針對噪聲方差進行邊緣后驗推導,結果如下:

f(β|Y,X,p,γ)∝

(22)

式(22)為常規逆Gamma分布,可直接根據該后驗分布式進行Gibbs采樣[15]。對于可變GGD先驗形狀參數p和尺度參數γ,由式(16)、式(18)和式(13)可得邊緣后驗分布為

f(p|Y,X,γ,β)∝

(23)

(24)

f(X|Y,p,γ,β)∝

(25)

ULA算法源自分子動力學,可用Langevin方程來表示隨機粒子在流體中的布朗運動[28],其表達式為

(26)

式中:X代表粒子某一時刻的位置;m為粒子質量;-λ(dX/dt)表示流體分子粘滯力;η(X)表示分子熱運動帶來的熱漲落力,通常表示為零均值高斯分布。求解式(26)即可得到粒子位置X分布。因此,模擬Langevin分子動力學,可以得到Langevin采樣估計算法,該算法可以對以下形式對數凸分布進行采樣

f(X)∝e-U(X)

(27)

(28)

式中:ε表示更新步長;X(r-1)、X(r)分別為X第r-1次和第r次迭代后的采樣值。考慮勢能函數U(X)中含有p范數,此時勢能函數U(X)很可能為非平滑函數,用ULA算法可能無法求解。因此,引入近端算子來解決非平滑勢能函數不可微問題,其標準定義形式為

(29)

式中:U(X)為非平滑不可微勢能函數;proxU(v)表示U(X)在v處的近端算子,近端算子可理解為某種求梯度的步驟,其中ε表示步長。由式(28)及近端算子與梯度之間的關系[29-30],近端ULA(proximal unadjusted langevin algorithm, P-ULA)中X的更新公式為

(30)

式中:proxU(X(r-1))為X在第r-1次更新時的近端算子,可用對偶前向后向(dual forward-backward, DFB)算法求解[29-31]。

4 貝葉斯增強算法流程

經過上述分析與討論,所提雷達鳥類目標微多普勒貝葉斯增強算法流程如圖5所示,主要包括初始化參數、超參數采樣與微多普勒特征采樣3個步驟。

圖5 雷達鳥類目標微多普勒特征貝葉斯增強算法流程圖Fig.5 Flowchart of radar Bayesian enhanced algorithm of micro-Doppler feature of bird target

具體解釋如下:首先,對待更新的參數集進行初始化,同時設置算法總迭代次數、Burn-in期。Burn-in期指在參數采樣估計過程中,樣本未收斂至目標分布前的階段,所以在更新參數時應舍去該階段樣本數。

接著,進行超參數采樣求解,如圖5中左側虛線方框所示。根據式(22),利用Gibbs采樣器對噪聲方差進行迭代采樣更新。根據式(23)和式(24)分別利用MH和Gibbs采樣算法迭代更新GGD先驗形狀參數p和尺度參數γ。

5 實驗驗證

根據上節所述概率模型,本文所提貝葉斯增強算法可以實現不同場景下鳥類目標微多普勒特征的增強。需要明確的是,本算法強調在時頻域進行特征增強,不涉及時頻分析前校正距離徙動等內容。換言之,本文所提算法僅僅考慮在進行時頻分析時,在時頻域上直接進行鳥類目標微多普勒特征增強。本章由仿真數據實驗與實測數據實驗組成。其中,仿真數據實驗選擇強雜波噪聲、窄窗超分辨兩種仿真場景;實測數據實驗將以X波段雷達探測飛行鴿子數據為例,通過與STFT算法、凸優化算法[14]進行對比,驗證了本文所提雷達鳥類目標微多普勒貝葉斯增強算法的有效性和優越性。

仿真數據實驗中雷達工作在X波段,發射帶寬為600 MHz,脈沖重復頻率約為820 Hz,脈沖數為8 192,目標觀測時間約為10 s。鳥類目標翅膀撲翼頻率為1 Hz,正向平移速度為1 m/s,翅膀長度為1 m,上臂和前臂長度均為0.5 m,上臂撲打角度幅值為40°,上臂撲打角度滯后15°,前臂撲打角度幅值為30°,前臂撲打角度滯后40°,前臂掃角幅值為20°。在本文所提貝葉斯增強算法中,近端算子的步長為自適應調節,設置Gibbs采樣老化期為40,設置總迭代次數為50。根據第1.1節給出的鳥類目標撲翼模型及以上參數,建立如圖6所示的鳥類目標運動仿真簡化模型[19]。其中,圖6(a)為鳥類目標翅膀撲翼圖,圖6(b)為鳥類目標翼尖及肘關節飛行軌跡。圖7為雷達追蹤鳥類目標仿真及鳥類目標撲翼距離像。其中,圖7(a)為雷達追蹤鳥類目標仿真,雷達位置為(20,0,-10),圖7(b)為鳥類目標撲翼距離像。

圖6 鳥類目標運動仿真簡化模型Fig.6 Simplified model of simulating motion of bird target

圖7 雷達追蹤鳥類目標仿真及鳥類目標撲翼距離像Fig.7 Simulation of radar tracking bird target and range profile of flapping bird target

5.1 仿真數據實驗

本節由強雜波噪聲與窄窗超分辨仿真實驗組成,并定量、定性分析每組仿真實驗的結果。

5.1.1 強雜波噪聲仿真實驗

為驗證本文所提算法在抑制噪聲的同時保留微多普勒特征完整性方面的有效性,本節以信噪比(signal-noise ratio, SNR)分別為0 dB、5 dB場景為例,將本文算法與STFT算法、凸優化算法作對比。圖8和圖9給出了SNR=0 dB、SNR=5 dB下運用STFT算法、凸優化算法以及本文所提算法得到的鳥類目標微多普勒特征結果。當SNR=0 dB時,圖8(a)為采用STFT算法得到的結果,不難看出背景噪聲斑點較多,鳥類目標微多普勒特征淹沒在背景噪聲中,難以有效識別鳥類目標。圖8(b)為運用凸優化算法進行稀疏特征增強后的結果。從該鳥類目標微多普勒特征圖中可以看出,盡管經過稀疏特征增強處理后,大量背景噪聲被抑制,但鳥類目標微多普勒特征的連續性被破壞。圖8(c)為運用本文所提算法提取鳥類目標微多普勒特征的結果,相較于圖8(a)、圖8(b),圖8(c)不僅抑制背景噪聲效果明顯,而且保留了微多普勒特征更多的連續性,其連續性特征更加明顯。

圖8 SNR=0 dB時仿真數據實驗結果對比Fig.8 Comparison of experimental results of simulation data with SNR=0 dB

為定量描述不同算法恢復微多普勒特征的性能,引入歸一化均方誤差(normalized mean square error, NMSE),其數學公式[15]為

(31)

當SNR=5 dB時,不同算法提取的鳥類目標微多普勒特征結果如圖9所示。從圖9(a)、圖9(b)、圖9(c)可以看出,隨著SNR的增大,所有算法恢復的鳥類目標微多普勒特征完整性均比圖8有所提升。特別地,對于利用本文所提貝葉斯增強算法恢復的圖9(c)結果,在有效抑制背景噪聲的同時,保留了更多微多普勒特征的連續性。經過MCM仿真計算,圖9(a)、圖9(b)、圖9(c)對應算法的NMSE分別為0.490 7、0.505 4、0.234 7,由此可見本文所提貝葉斯增強算法微多普勒特征恢復精度高。觀察以上兩組仿真實驗結果,鳥類目標微多普勒特征由多組重疊類正弦波多普勒圖形組成,每組多普勒圖形中多普勒頻偏最大的類正弦波圖形為鳥類翅膀翅尖撲打所產生,由此可推算出該鳥類目標翅膀總長度為1 m,符合預設。而重疊在內部的類正弦波圖形由鳥類翅膀前后臂間的節點撲打所產生,由此可推算出該鳥類目標翅膀前臂長度為0.5 m,符合預設。

圖10給出了在一定SNR范圍內,STFT算法、凸優化算法、貝葉斯增強算法恢復鳥類目標微多普勒特征性能曲線的對比示意圖。繼續沿用MCM仿真,采用NMSE評價不同算法在不同SNR條件下恢復鳥類目標微多普勒特征的性能。由圖10可知,隨著SNR不斷增加,STFT算法的NMSE呈直線式下降,表明STFT算法易受SNR影響,對噪聲敏感,算法恢復穩定性差。而凸優化算法雖然在低SNR條件下恢復性能優于STFT算法,對噪聲不敏感且性能曲線平穩,但是在高SNR條件下,其僅能進行稀疏特征增強的凸優化算法,不能夠滿足恢復精度的要求。相比前兩種算法,本文所提貝葉斯增強算法,不僅性能曲線相對平緩,不易受SNR影響,而且在NMSE數值上明顯優于其他兩種算法。綜上所述,在不同SNR條件下,本文所提貝葉斯增強算法不僅能夠以較高精度有效提取鳥類目標微多普勒特征,而且具有一定的性能穩定性。

圖9 SNR=5 dB時仿真數據實驗結果對比Fig.9 Comparison of experimental results of simulation data with SNR=5 dB

圖10 不同SNR條件下不同算法的性能曲線對比Fig.10 Comparison of performance curves of different algorithms under different SNR conditions

5.1.2 窄窗超分辨仿真實驗

為驗證所提算法在超分辨方面的有效性,在SNR為5 dB前提下,將強雜波噪聲仿真實驗中采用的短時窗函數長度縮短至原來長度的50%,再將本文所提算法與STFT算法、凸優化算法作對比,考察各個算法恢復鳥類目標微多普勒特征的性能。圖11(a)為運用STFT算法提取鳥類目標微多普勒特征的結果,由圖11(a)可見,短時窗函數變窄使得微多普勒特征變得模糊,即頻率分辨率降低,同時STFT算法無法抑制背景噪聲,難以通過該微多普勒特征圖準確獲取與鳥類目標相關的特征信息。圖11(b)為運用凸優化算法進行稀疏特征增強后的結果,此算法雖然能夠有效抑制背景噪聲,但同時也抑制了微多普勒特征的連續性,使得微多普勒特征不夠連續,超分辨效果一般。圖11(c)為運用本文所提貝葉斯增強算法進行特征增強后的結果。從圖11中可以看出,該算法既能有效抑制背景噪聲,同時又盡可能保留了微多普勒特征中更多的連續性,相較于圖11(a)、圖11(b),本文所提貝葉斯增強算法所恢復的微多普勒特征連續性更加明顯,超分辨效果較好。經過MCM仿真計算,圖11(a)、圖11(b)、圖11(c)對應算法的NMSE分別為0.674 1、0.553 3、0.510 8,驗證了本文所提貝葉斯增強算法在超分辨方面的有效性。

圖11 SNR=5 dB時窄窗超分辨仿真數據實驗結果對比Fig.11 Comparison of super-resolution experimental results of simulation data with SNR=5 dB and narrow window

5.2 實測數據實驗

本節將雷達探測飛行鴿子實測數據用于驗證鳥類目標微動信號模型的通用性與所提貝葉斯增強算法的有效性。該實測數據對應陣列雷達系統工作在X波段,發射信號瞬時帶寬為16 MHz,陣列單元數目為16,工作模式為單發16。圖12給出了該鴿子目標探測響應結果。其中,下面的線表示單天線目標探測響應結果,上面的線表示多天線目標探測響應結果。很明顯,通過觀察單天線目標探測響應結果,可疑距離單元445處的鴿子目標淹沒在了近場強雜波中。在多天線目標探測響應結果中,雖然該鴿子目標探測響應強度得以提高,但近場強雜波的影響依然顯著存在。因此,利用傳統CFAR檢測方法不易探測類似鴿子這種低慢小目標,而且容易導致虛警率高。接下來,考慮提取該鴿子目標微動特征來實現對其的探測。

圖13給出了飛行鴿子的微多普勒特征圖,其中圖13(a)、圖13(b)、圖13(c)分別為由STFT算法、凸優化算法和本文所提貝葉斯增強算法提取的微多普勒特征結果。為了能夠一致地評價各個算法的恢復效果,對每個算法恢復的微多普勒特征結果進行歸一化處理。在圖13(a)中,STFT算法作為線性運算,雖然運算效率快,可以實現多分量微多普勒特征提取,但是其恢復的微多普勒特征相對模糊,微多普勒特征周期性不明顯,難以有效從圖中獲取該鴿子目標的特征信息。相比STFT算法,在圖13(b)中,由凸優化算法恢復的微多普勒特征雖然在頻率分辨率方面得到一定程度提升,但是造成了微多普勒連續性特征損失,同樣不利于從其中準確提取該鴿子目標特征信息,甚至可能會遺漏重要的特征信息。相比前兩種算法,由本文所提貝葉斯增強算法提取的微多普勒特征在有效改善頻率分辨率的同時,保留了微多普勒特征更多的連續性,提高了獲取該鴿子目標相關信息的準確性。由圖13(c)所示,觀測目標總時間約為3.3 s,由此可推算出鴿子撲翼的頻率約為5 Hz,屬于鴿子飛行時正常撲翼頻率范圍。因此,通過本文所提貝葉斯增強算法可以獲得鳥類目標撲翼周期,以用于后續鳥類目標的檢測與識別。綜上所述,經過雷達探測飛行鴿子實測數據的驗證,本文所提貝葉斯增強算法在處理實測數據時仍然具有一定的有效性。

圖12 鳥類目標探測響應圖Fig.12 Figure of detection response of bird target

圖13 實測數據實驗結果對比Fig.13 Comparison of experimental results of measured data

6 結 論

針對鳥類目標RCS小、近場強雜波干擾等因素造成探測鳥類目標困難的問題,本文提出運用鳥類目標撲翼運動產生的微多普勒特征對鳥類進行探測與識別。首先建立鳥類目標撲翼模型,繼而推導出鳥類目標雷達回波及微多普勒模型,最后針對STFT線性時頻分析方法獲取的微多普勒特征分辨率低、特征不明顯以及場景復雜多變等問題,提出利用GGD對先驗進行自適應建模,借助貝葉斯推理,采用P-ULA算法高效求解后驗概率密度函數,在時頻域實現鳥類目標微多普勒特征增強。經過仿真數據實驗驗證,在低分辨率、高噪聲的情況下,本文所提算法可以得到特征明顯的微多普勒特征。同時,本文所提算法在實測數據上依然可以有效提升微多普勒特征的分辨率,保留微多普勒特征更多的連續性。本文為雷達鳥類目標識別提供了新思路,在此基礎上,可以研究更加復雜的鳥類模型以及提高時頻分辨率的算法。同時,在利用本文所提貝葉斯增強算法進行特征增強時,雖然保留了連續性特征,但是損失了一部分結構性特征,有必要研究適用于更多特征的增強算法,這將是后續研究工作的重點。