低信噪比條件下直升機微動特征提取方法

張朝偉, 夏賽強, 楊 軍, 劉建衛, 徐穎鑫

(空軍預警學院， 湖北武漢 430019)

0 引 言

微多普勒概念由Chen引入到雷達觀測領域，開創了基于微多普勒特征的雷達目標檢測與識別新領域[1-4]。

帶有螺旋槳部件的空中復雜運動目標，自身發生平動的同時伴隨著螺旋槳的旋轉運動，具有典型微動形式，被國內外學者廣泛研究并不斷深入，取得了一系列研究成果。時頻分析方法[5-7]是一種最直觀最常用的微多普勒特征分析方法，但其時頻分辨率受窗函數的約束，存在分辨率和交叉項的矛盾問題。文獻[8]聯合時頻分析方法和圖像處理方法，提取出時頻圖中曲線特征來獲得目標微動特征，但是在低信噪比條件下失效。文獻[9]通過求取循環相關系數，實現了低信噪比條件下的微動特征提取，但是在閃爍情況下性能不佳。文獻[10]研究了一種基于自相關處理的閃爍情況下微動特征提取方法，該方法運算量較大，信號重構能力弱，不易于工程實現。文獻[11]研究了經驗模態分解方法，將復雜的微多普勒時頻曲線快速分解為多個模態分量，從分量中可以獲得目標的運動參數，在提取進動目標微多普勒特征時具有優勢，但是在發生模態混疊現象時，估計誤差較大。文獻[12-13]研究了高階矩函數方法，通過對回波數據進行處理并計算其高階矩函數來得到旋轉頻率和旋轉半徑的估計，具有計算復雜度小、運算速度快的優點，但抗噪性性能較差。文獻[14-15]研究了正交匹配追蹤方法，通過將信號在過完備原子庫上稀疏分解來得到微動參數的估計，具有更為優良的魯棒性，抗噪性能好，但隨著估計參數的增加，運算量會急劇增加。

為解決信噪比適應性與估計精度和運算量的矛盾問題，本文以空中旋翼目標為研究對象，對其微多普勒效應進行建模仿真，深入分析微動特征參數的內在聯系，并使用時頻分析聯合正交匹配追蹤方法對微多普勒特征進行提取。仿真和實驗結果表明，本方法能夠實現低信噪比條件下旋翼目標微動特征的快速高精度提取。

1 旋翼目標回波模型

不失一般性，假設目標回波已完成運動補償，建立雷達坐標系O-XYZ，如圖1所示。僅考慮直升機主旋翼，該旋翼具有N個葉片，以轉動角頻率ω旋轉，旋翼葉片長度為l，雷達到目標旋翼中心距離為R0，方位角為α，俯仰角為β。

圖1 散射點模型示意圖

將葉片看作積分模型，進行單基地直升機旋翼建模，其運動補償后回波[16]可表示為

(1)

式中，

(2)

(3)

σ為目標散射系數，λ為雷達波長，fr為目標旋翼轉速，φ0為第一個葉片初相，整數k表示第k個葉片(1≤k≤N)。在式(1)可以看出，時域信號幅度受到sinc(·)函數調制，第k個葉片引起的瞬時多普勒頻率為

(4)

(5)

將式(5)代入式(3)可知該葉片的相位可表示為

(6)

兩相鄰閃爍之間的時間間隔滿足

(7)

其中，%為取余數運算，當出現閃爍時，瞬時多普勒頻率取到最大值：

(8)

2 目標閃爍參數估計

直升機旋翼回波的微多普勒呈現非線性變化，對于這類信號通常采用時頻分析方法來揭示信號頻率的時變特性。Gabor變換是一種常用時頻分析方法，分辨率較好，且不存在交叉項[17]。將旋翼直升機回波信號變換到時頻域：

(9)

其中，w(·)是高斯窗函數。旋翼微多普勒效應特征曲線在時頻域表現為余弦調制曲線，在時域出現閃爍現象的時刻，在時頻域會出現垂直于時間橫軸的頻率帶，即時頻域閃爍現象。

圖2和圖3分別是直升機葉片數為奇數情況和偶數情況下的閃爍現象示意圖。從圖2(a)和圖3(a)可以看出，無論奇葉片數還是偶葉片數，在時域波形上并不能進行有效區分，但是在時頻域，當直升機葉片數為奇數時，時頻圖中閃爍帶交替出現，如圖2(b)所示，當直升機葉片數為偶數時，時頻圖中的閃爍帶則是同時出現，如圖3(b)所示。

(a) 時域閃爍

(a) 時域閃爍

同理，將時頻圖數據幅值按頻率進行累加，通過尋找峰值包絡的寬度，可以計算得到目標微動的最大瞬時多普勒頻率max。

由式(6)可知，葉片的相位與轉速和葉片數密切相關。由式(7)知葉片的轉速與葉片數和閃爍間隔密切相關。由式(8)知葉片長度與轉速和最大多普勒頻率密切相關。此時聯立式(2)、式(6)、式(7) 和式(8)可得

3 旋翼目標參數估計

由式(1)知時域回波信號可分解為

(11)

式中，cm為原子系數，gm為第m個原子，Λ為待估計的參量，D為字典矩陣,D=[g1g2g3…gM]∈CNt×M，M為原子個數，Nt為采樣點數，α∈CM為系數矢量，是稀疏的。將參數估計問題轉化為最優gm范數問題進行稀疏向量求解，正交匹配追蹤算法常用來解決該問題，通過構建字典矩陣，將信號在過完備原子庫上進行稀疏分解，不斷選取與信號最匹配的原子進行稀疏逼近，可實現低信噪比條件下的多微動參數估計。

由式(1)知直升機旋翼的回波信號由參數(fr,N,φ0,l)確定，同時對4個參數進行估計將會導致運算量大幅提高，且估計精度也會受到影響，若將已估計得到的目標閃爍參數作為先驗信息，那么時域回波僅由N確定，葉片數通常在3～8個，且為整數值，大大減少了估計難度，運算量減少了Cfr·CN·Cφ0倍。由OMP算法原理可知，字典中原子可通過待分解信號的內在特性來構造，目標微多普勒信號是一個正弦調頻信號，因此選擇SFM原子來構造過完備原子庫，第m個原子可表示為

(12)

式中，m=pNt。對原子集里每個原子進行歸一化處理，即

(13)

算法實現流程[18-19]如下：

輸入： 測試樣本r0=s，訓練字典D=[g1g2g3…gM]，稀疏度K。

初始化： 初始化殘差余量r0=s，匹配原子集賦空值Λ0=?，匹配原子記錄矩陣賦空值D0=?，迭代初始值k=1。

2) 將尋找到的最相關字典元素的索引λk加入匹配原子集Λk=Λk-1∪{λk}；

4 仿真分析

為了更清晰地分析本文方法對旋翼目標參數的估計性能和可行性，參考當前窄帶雷達系統的相關參數確定雷達仿真參數，并選取兩種常見直升機的真實參數作為目標仿真參數，脈沖壓縮后信噪比為10 dB，仿真雷達參數設置如表1所示，仿真目標參數如表2所示。

表1 雷達參數

表2 目標參數

利用SA341直升機真實參數進行仿真，其仿真結果如圖4所示。從圖4(a)可以清晰觀察到正負頻率部分的時頻圖幅值縱向積累曲線峰值交替出現，即可判斷出目標葉片數目為奇數，且曲線相對于時域閃爍曲線更加平滑，有利于尋找峰值點位置，此時第一個閃爍時刻為0.014 s，為負頻率閃爍，第十八個閃爍時刻為0.486 s，為正頻率閃爍，通過這兩個時刻可以估計出閃爍間隔時間為(0.486-0.014)/17=0.027 8 s，通過首尾兩個閃爍時間差求平均值的方式求取閃爍間隔相比于直接計算相鄰兩個閃爍時間差作為閃爍間隔的方式更為精確。圖4(b)是時頻圖幅值橫向積累曲線，由于時頻圖中曲線具有一定寬度以及一定的邊緣效應，導致圖中曲線具有明顯的上升沿和下降沿，且持續時間較長，不利于尋找最大瞬時多普勒頻率數值，對于這類信號，通常采樣3 dB寬度來衡量其包絡的寬度，此時估計出目標最大瞬時多普勒頻率為1 315 Hz，與理論值1 319.5 Hz基本一致，誤差為0.34%。將估計得到的第一個閃爍時刻、閃爍間隔和最大瞬時多普勒頻率作為先驗信息，利用OMP方法估計葉片數，結果如圖4(c)所示，可以看出在葉片數為3時取到最大值，與預設值一致。利用式(7)可估計出目標葉片轉動頻率為5.995 2 Hz，與預設值誤差為0.08%，利用式(8)可估計出目標旋翼半徑為5.236 4 m，與預設值誤差為0.26%。

(a) 時頻圖幅值縱向積累曲線

AH-64直升機仿真結果如圖5所示。從圖5(a)觀察到正負頻率部分的時頻圖幅值縱向積累曲線峰值同時出現，即判斷目標葉片數目為偶數，此時第一個閃爍時刻為0 s，第十個閃爍時刻為0.468 6 s，估計出閃爍間隔為(0.468 6-0)/9=0.052 s，從圖5(b)中估計出目標最大瞬時多普勒頻率為1 468 Hz，與理論值1 469.8 Hz基本一致，誤差為0.12%。圖5(c) 是利用OMP方法對葉片數的估計結果，在葉片數為4時取到最大值，與預設值一致。此時，目標葉片轉動頻率估計值為4.807 7 Hz，與預設值誤差為0.16%，目標旋翼半徑估計值為7.289 5 m，與預設值誤差為0.28%。

(a) 閃爍參數

從以上兩組仿真可以看出，本文方法能夠較為準確地估計出目標的旋翼參數，在信噪比為10 dB情況下，其估計誤差小于0.3%。以SA341仿真為例，圖6給出了不同信噪比情況下其目標閃爍參數估計精度和旋翼參數估計精度。從圖6(a)可以看出，閃爍間隔被準確估計，而最大瞬時多普勒頻率估計準確度隨信噪比的降低整體上呈現下降趨勢，且在信噪比小于8 dB時下降較為明顯。在圖6(b)中，葉片長度估計準確度曲線與圖6(a)中最大瞬時多普勒頻率估計準確度基本一致，這是由于在估計葉片長度時使用了最大瞬時多普勒頻率估計值，導致估計誤差向下傳遞。同理，目標轉動頻率估計值與葉片數估計值和閃爍間隔估計值息息相關，而這兩個參數估計準確度在不同信噪比條件下均能保持很高的準確度，因此目標轉動頻率估計較為準確。實際上，OMP算法對網格匹配程度要求較高，一旦網格失配就會導致葉片數估計不準確，從而導致轉動頻率和葉片長度估計失敗，通過仿真發現OMP算法對閃爍間隔估計準確度的要求要遠高于瞬時多普勒頻率估計準確度，即式(10)中sin(·)函數項相比其系數對OMP算法性能的影響更大。

5 實驗分析

為驗證本文方法在實測數據中的有效性，利用一組實測數據進行微動特征提取實驗，輻射源為某型窄帶雷達，目標為某型直升機，其公開參數為葉片數5個，葉片長度8.6 m，旋翼轉動頻率4 Hz。

在實測數據中，目標直升機背站飛行，且外場環境存在強地物雜波，因此需要對回波進行雜波抑制和運動補償，其結果如圖7所示。圖7(a)為雜波抑制后的時域回波，能夠觀察到幅度不一的時域閃爍現象，其中幅度相對較強的閃爍共計出現七次。對時域回波進行運動補償并變換到時頻域，結果如圖7(b)所示，時頻閃爍現象明顯，在[-300 Hz,-600 Hz]范圍內還存在覆蓋范圍較廣的調制頻帶，這是由于葉轂跟隨葉片一起旋轉以及尾槳旋轉導致的，但是葉轂和尾槳相對主旋翼而言具有尺寸小，轉動半徑小的特點，其微多普勒頻率較低，且沒發生閃爍現象，從而不會影響本文方法對直升機葉片的特征提取。

圖8是本文方法的特征提取結果，從圖8(a)可以觀察到正負頻率部分的時頻圖幅值縱向積累曲線峰值交替出現，即判斷目標葉片數目為奇數，此時第一個閃爍時刻為0.161 s，最后一個閃爍時刻為0.315 s，估計出閃爍間隔為0.025 7 s，在圖8(b)中正負頻率上的最大值并不嚴格對稱，這是由于目標姿態導致的，此時對其取均值估計出目標最大瞬時多普勒頻率為1 396 Hz，圖8(c) 是利用OMP方法對葉片數的估計結果，葉片數估計值為5。將上述三個參數估計值帶入公式進行計算，目標葉片轉動頻率估計值為3.891 Hz，與目標公開參數誤差為2.72%，目標旋翼半徑估計值為8.565 m，與目標公開參數誤差僅為0.41%。

(a) 時域波形

(a) 時頻圖幅值縱向積累曲線

6 結束語

本文從目標微動特性出發，提出了一種低信噪比條件下的直升機微動特征提取方法。通過對直升機回波進行Gabor變換，獲取其時頻信息，采用時頻圖幅值積累法較好地獲取旋翼回波閃爍參數，如葉片奇偶數、閃爍周期和最大瞬時多普勒頻率。并將其作為先驗信息，對稀疏字典中原子進行降維處理，將多參數搜索優化為單參數搜索，極大地降低了正交匹配追蹤算法的運算量，葉片數估計準確度在網格不失配情況下達到100%，其目標葉片數和轉動頻率估計準確度達到99.9%，葉片長度估計準確度受最大多普勒頻率估計誤差傳遞的影響，能夠保持在95%以上，且在低信噪比條件下保持優良的魯棒性。實驗結果表明本文方法具有較強的實用性，為后續直升機目標識別奠定了基礎，具有一定的工程應用價值。