一種變步長的微動目標參數高精度提取方法

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

0 引言

目標中存在旋轉、振動、進動和章動等微動形式時，雷達回波將產生除目標剛體部分引起的多普勒頻率外的附加頻率調制，即微多普勒效應(Micro-Doppler Effect)[1]。自2000年 Chen提出微多普勒概念以來，對微動目標微多普勒現象的研究迅速發展起來，其中，旋翼旋轉作為典型的微動形式，旋翼微動信息中包含著旋翼的轉速和尺寸等信息，利用旋翼轉速和尺寸可以對直升機一類目標進行識別，因此受到國內外學者的廣泛研究[2-7]。當前窄帶雷達對旋翼轉速估計的主要方法有時頻分析方法[8]、正交匹配追蹤算法[8]、Hough變換法[8]、高階矩函數分析法[9-10]、函數構建法[11]和自相關法[12]。其中，正交匹配追蹤算法、高階矩函數分析法需要用到回波的相位信息；Hough變換計算量大，無法快速進行轉速估計；函數構建法雖然計算量較小，但是在旋翼散射點均勻分布時會失效；自相關法必須在葉片個數已知的情況下才能準確估計旋翼轉速，并且對于形狀規則的直升機旋翼轉速估計時，真實轉速與估計轉速存在著整數倍誤差，該誤差倍數等于旋翼對稱軸個數。因此，旋翼目標的參數估計方法仍需進一步研究。

針對上述問題，本文提出了一種變步長的微動目標參數高精度提取方法，首先建立了旋翼目標的散射點模型，然后通過Gabor變換以及時頻圖骨架提取估計出旋翼最大瞬時微多普勒頻率，接著利用此先驗信息，對構建的搜索矩陣進行降維處理，并對轉速進行變步長處理，即進行粗搜索和精搜索。然后對搜索矩陣與回波矩陣進行能量積累，可得到積累效果最好的轉速值，此時所對應的轉速就是估計的最佳轉速，隨后利用轉速估計結果進一步估計出旋翼目標的葉片尺寸。通過本文方法，可以準確估計出旋翼轉速和尺寸，并且相對于遍歷搜索方法，本文方法明顯減少了算法的運算量。

1 旋翼微動目標信號建模分析

基于窄帶雷達，在觀測時間較短的情況下，旋翼目標不會出現距離單元走動，若觀測時間較長，旋翼目標出現距離單元走動現象，此時對旋翼目標進行運動補償后，旋翼目標的運動狀態可等效為懸停狀態[13]。為簡化分析，假設目標與雷達處于同一個平面，如圖1所示，雷達到旋翼目標旋轉中心的距離為RC，假設旋翼微動目標為散射點模型，旋翼某一個葉片上的散射點P到旋翼中心的距離為r(0≤r≤l，l為旋翼葉片長度)，P點到雷達的距離為RP，同時P點以角速度ω繞旋轉中心C旋轉，初始時刻P點的初始旋轉角為θ。

圖1 旋翼目標模型

設雷達發射信號為窄帶線性調頻信號[2]：

(1)

旋翼葉片上散射點的回波信號經脈沖壓縮后可表示為

(2)

式中，K為葉片數目，i為旋翼上的第i個葉片，N為目標單個葉片上的散射點數目，j為葉片上的第j個散射點，σij為散射點的散射系數，B為發射線性調頻信號帶寬，λ為發射信號波長，Rij(tm)為旋翼上散射點到雷達的距離[13]。

Rij(tm)=RC+rijcos(ωtm+θij)

(3)

由式(2)可以看出，旋翼目標微多普勒相位為

(4)

由式(4)可知，散射點的微多普勒呈現余弦函數形式，而余弦函數的周期與旋翼旋轉周期是一致的。且不同葉片上散射點引起的微多普勒曲線有著不同的初相，而同一葉片散射點初相是一致的，差異在于多普勒頻率峰值大小的不同。對應式(4)，可知旋翼目標瞬時微多普勒頻率為

(5)

2 旋翼目標微動特征提取

本文提出的變步長旋翼微動特征提取方法共分3步：第一步是通過Gabor變換和圖像骨架[5]提取的方法估計出目標的最大瞬時多普勒頻率；第二步是在第一步的基礎上通過降維處理以及將搜索過程進行變步長處理，即分為粗搜索和精搜索；第三步是利用能量積累的思想，在時頻域對回波數據矩陣與降維后的搜索矩陣進行能量積累，通過積累的效果來估計旋翼目標的轉速，最終通過旋翼轉速與最大瞬時多普勒頻率和旋翼尺寸的嚴格數學關系來估計旋翼尺寸。

2.1 基于Gabor變換的時頻分析及骨架提取方法

由式(2)可以看出，旋翼轉動的微多普勒是一個非平穩時變信號，因此需要采用局部分析方法即時頻分析來對信號進行處理。由于在窄帶雷達中，散射點的微動幅度通常在一個距離分辨單元內，其一維距離像表現為一條平行于方位向的直線，只需要取出該距離門的數據進行時頻分析即可。常用的時頻分析方法有短時傅里葉變換(Short Time Fourier Transform, STFT)、Wigner-Ville分布(Wigner-Ville Distribution, WVD)、平滑偽WVD (Smoothed Pseudo Wigner-Ville Distribution, SP-WVD)、S方法(S-Method, SM)等[14]。其中，STFT是最常用的基于匹配濾波的時頻分析方法，當STFT選取的窗函數為高斯窗時，稱之為Gabor變換，根據測不準原理，在所有可能的窗函數中，高斯窗函數能得到最好的時頻效果。

對式(2)進行Gabor變換可以得到微動目標信號的時頻圖，其表達式[9]為

(6)

由式(5)可以看出，在時頻圖上，旋翼目標的微多普勒曲線是一個三參數的正弦曲線，信號的能量集中在曲線上。圖2給出了一個三旋翼目標，且每個旋翼上僅有葉尖處有散射點時的時頻特征曲線示意圖。

圖2 時頻特征曲線示意圖

由于時頻分析受到其分辨率的影響，直接提取參數的誤差較大，為了減少后續處理的運算量以及提高參數提取的準確度，從圖像域的角度出發，對時頻圖圖像進行閾值處理、平滑處理和二值化處理，而后提取出微動目標正弦曲線的骨架，此時的圖像矩陣中只包含“0”和“1”兩個值，其中“1”出現的位置即為特征曲線的骨架。這樣的骨架同樣表征了信號的頻率隨時間的變化規律，提取骨架的表達式為

(7)

由式(7)可知，此時可以通過時頻圖骨架準確地找出曲線的極值位置，從而提取出微動目標的最大瞬時微動頻率：

(8)

從式(8)可以看出，對于旋翼目標的轉速和尺寸參數估計，可以先估計其中一個參數，然后利用嚴格的數學關系估計出另外一個參數。下面本文提出一種基于該原理的旋翼目標轉速和尺寸估計方法。

2.2 基于降維后的變步長微動參數高精度估計方法

目前除了傳統的轉速估計方法外，還有許多利用微多普勒正弦調制曲線為先驗知識的微動參數提取方法，如復經驗模式分解(CEMD)、匹配傅里葉變換、擴展Hough變換等，但是上述方法均計算復雜，運算量大。

針對上述問題，本文提出一種降維后的變步長微動參數估計方法，該方法將對微動目標尺寸參數和轉速參數的二維聯合估計，通過將提取到的微動目標的最大瞬時微動頻率作為先驗信息，對其轉化為分別估計。在估計轉速時，首先構建搜索函數：

(9)

為了方便分析，這里假設目標初相θ已知，由分析知道旋翼目標的微多普勒曲線是一條正弦曲線，其表達式可以寫為

fm-D=fm-Dmaxsin(ωtm+θ)

(10)

由式(5)、式(10)可知，當r為最外側散射點時，瞬時微多普勒頻率取到最大值，此時

(11)

將搜索函數改寫為

(12)

那么搜索矩陣可表示為

(13)

由于構建的搜索矩陣僅利用了一個葉片最外側散射點的最大瞬時微多普勒頻率，因此可以在不知道葉片個數的情況下進行轉速的搜索。對HQ×M每行進行Gabor變換，可以得到q個時頻數據矩陣V：

V=[V1,V2,…,Vq-1,Vq]

(14)

由于時頻分析方法是利用時間變量和頻率變量來表征信號的能量分布密度，因此可利用回波信號經時頻分析后的數據矩陣J與搜索函數的時頻分析矩陣V進行點乘來進行能量的積累，能量積累的效果僅與搜索函數中的ω有關。

F=J·V(q)

(15)

(16)

由分析可知，當搜索函數中搜索步進Δω越小時，搜索結果越精確，但是由于本算法需要對每個轉速可能值進行時頻分析，會導致運算量急劇增大，因此本文在此前提下對搜索過程進行了粗估計和精估計的處理，即變步長處理，利用粗估計找到轉速真實值的大致范圍，此時估計精度較差，必須利用精估計對粗估計的轉速值附近范圍內的轉速值進行精搜索，最后得到轉速估計值。精搜索的次數越多，其估計的精度越高，但是同時會增加算法的運算量，因此本文在進行粗估計后僅進行兩次精估計搜索。

在轉速搜索范圍[ω1,ω2]內，假設一次Gabor變換的運算量為Nr，設置搜索步進為Δ進行粗搜索，所需的運算量為(ω2-ω1)Nr/Δ，而兩次精估計的運算量均為10Nr，故變步長的轉速估計方法總運算量為

(17)

如果直接采用遍歷式的搜索方法，且要達到相同的估計精度，則需要將搜索步進設置為0.01Δ，此時需要的運算量為

(18)

相對于遍歷式的搜索方法，變步長的搜索方法的運算量減少了

(19)

在估計出旋翼目標轉速后，通過式(8)、式(16)可以估計出目標的尺寸：

(20)

本文提出的基于Gabor變換的變步長微動目標參數高精度提取方法的處理流程如圖3所示。

圖3 微動目標參數快速提取方法流程

3 仿真分析

為了更加清晰得到本文方法對微動目標轉動頻率(frot=ω/2π)和尺寸的提取，仿真選擇3個旋翼葉片共10個散射點，模型如圖4所示，葉片長度為l=6 m，且本文定義的信噪比為脈壓后的信噪比。具體仿真參數如表1所示。

表1 仿真參數

圖4 旋翼葉片模型

圖5是SNR=5 dB，frot=1.41 Hz時，利用本文方法對旋翼目標轉速進行估計的結果；圖6是SNR=5 dB，frot=3.14 Hz時，旋翼目標轉速的估計結果；圖5(a)和圖6(a)是旋翼目標回波所在距離單元的Gabor變換時頻分析結果；圖5(b)和圖6(b)是對時頻圖進行骨架提取的結果；圖5(c)和圖6(c)是設置搜索步進為0.005 Hz對轉速進行遍歷搜索的搜索結果；圖5(d)和圖6(d)是將遍歷搜索方法改為變步長搜索方法后，設置搜索步進為0.5 Hz對轉速進行粗搜索的搜索結果；圖5(e)和圖6(e)是在粗搜索得到的最佳搜索結果附近，將搜索步進設置為0.05 Hz對轉速進行第一次精搜索的搜索結果；圖5(f)和圖6(f)是在第一次精搜索的最佳搜索結果附近，將搜索步進設置為0.005 Hz對轉速進行第二次精搜索的搜索結果。

在圖5(a)和圖6(a)中，可以看出此時目標的微動頻率是受正弦調制的，利用時頻圖可以提取出圖5(b)和圖6(b)中所示的時頻圖骨架，并且可以利用骨架估計出最大微多普勒瞬時頻率。對轉速進行遍歷搜索時，由于已經將搜索步長設置得足夠小，因此在圖5(c)和圖6(c)中可以看到其搜索精度很高。另一方面，正是由于其搜索步長很小，需要對搜索范圍內800個轉動頻率點進行搜 索，導致遍歷搜索過程很慢，因此需要提出一種既能保證搜索精度，又能提高搜索速度的算法，即變步長的搜索方法。將搜索方法改為變步長的搜索方式后，此時僅需要對30個轉動頻率點進行搜索，搜索量大幅下降。但是在進行粗搜索時，此時轉動頻率與搜索步進存在失配情況，本文方法只能夠估計出目標轉動頻率的大致范圍，估計誤差很大，因此還需要通過改變搜索步長對轉動頻率進行估計，其估計結果逐漸向真實轉速逼近，經過第二次精估計后，本文方法能準確地估計出旋翼目標的轉動頻率，如果旋翼目標的轉動頻率還存在多位小數的情況，但是本文在準確估計出小數前幾位的情況下，其估計誤差帶來的影響很小。

圖5 本文方法對frot=1.41 Hz的估計結果

圖6 本文方法對frot=3.14 Hz的估計結果

表2給出了本文方法在不同信噪比條件下的轉動頻率和尺寸估計誤差。

通過表2可以發現，本文提出的變步長微動目標參數提取方法對旋翼轉速和尺寸的估計精度和對信噪比的適應程度都具有較好的性能。但是在SNR<5 dB時，本文方法的估計誤差急劇增大，這是由于此時時頻分析方法已經無法得到有效時頻圖以及相應的骨架，從而無法估計出旋翼目標最大瞬時微多普勒頻率，缺失了這一先驗信息，本文方法將會失效。

4 結束語

本文提出了一種變步長的微動目標參數高精度提取方法，通過利用旋翼最大瞬時微多普勒頻率的先驗信息，對構建的搜索矩陣進行降維處理，然后從時頻域能量積累的角度出發，對轉速進行了粗估計和兩次精估計，積累效果最好的轉速值即為估計的最佳轉速，并通過轉速估計值和最大瞬時微多普勒頻率估計值估計出旋翼目標的葉片尺寸。通過仿真結果可以看出，本文方法可以較為快速地估計出旋翼轉速和尺寸，并且具有很高的估計精度，對信噪比的適應程度也較好。但是，該方法在進行估計時是基于骨架提取結果完成的，在低信噪比條件下該方法的估計精度將急劇下降甚至估計失效，下一步還需要從提高時頻分辨率和降噪處理的角度出發，來解決本文方法在低信噪比條件下估計精度不足的缺陷。