基于時頻特征的旋翼目標微動參數估計方法

宋村夫，曲智國，胡旭超

(1.武警上海總隊通信大隊，上海 200050；2.空軍預警學院，湖北 武漢 430014)

1 引言

在現代信息化戰爭中，旋翼目標由于體積小、受地形限制小、機動性強等優勢在軍事競爭中扮演著越來越重要的角色[1]。同時在民用領域，風力發電也正在蓬勃發展，風力渦輪機本質也為一種旋翼微動目標[2]。旋翼目標所產生的微動特征包含著目標的結構、運動特點等相關信息，通過其提取微動參數，對旋翼目標的探測、分類和識別具有重要的意義，因此在軍用和民用領域具有廣闊的研究場景[3][4]。2006年V.C.Chen提出將目標及其它組件除整體質心運動之外的旋轉、振動、以及翻滾等細微的運動形式稱為微動，同時提出目標微動會對雷達回波信號產生頻率調制，這種由微動引起的多普勒效應被稱為微多普勒效應[5]。

目前相關領域專家學者都提出了一些有關微動特征參數提取的技術，其主要包變換域的微動參數提取、基于圖像域的微動參提取、基于壓縮感知的微動參數提取等[6]。文獻[7]提出了一種基于變換域的短時迭代自適應-逆Radon變換(STIAA-IRT)的微多普勒特征提取方法，但是該方法步驟繁瑣，計算量大；文獻[8]提出了一種基于相位信息的微動參數提取方法，但是需要先解決相位纏繞問題；文獻[9]提出利用自相關函數峰值與回波間的對應關系實現了對微動的旋轉頻率的提取，但是該方法只用于單一周期的回波，在實際應用中局限性大，適應性不高。

針對以上問題，為了進一步有效地提取旋翼目標的微動參數，本文提出了一種基于時頻特征的旋翼目標微動參數的提取方法，對葉片旋轉頻率、葉片初始相位和葉片長度實現了快速提取。首先，構建了旋翼目標雷達回波模型，在此基礎上利用分析了旋翼目標回波的時頻特征，并利用該時頻特征分析了“時頻閃爍”與微動參數間的對應關系，最后在時頻域利用圖像特征實現了微動參數的提取。通過仿真表明，該方法不僅可以實現對旋翼目標葉片的微動參數的快速、高精度的提取，同時該方法還可適用于多個旋翼目標場景下的微動參數提取，且具有一定的抗噪性能。

2 旋翼目標雷達回波建模與分析

以旋翼目標為三葉片構造為例，為簡化分析，假設旋翼目標和雷達處在同一平面。以旋翼目標的葉片旋轉中心為原點，建立如圖所示坐標系。其中，B為雷達站，r為雷達視線(LOS)，θ葉片旋轉角，R為雷達中心距離旋翼目標旋轉中心距離，li為葉片上散射點距離旋轉中心的長度，li?R。

圖1 旋翼目標回波模型

假設雷達發射信號p(t)=exp(-j2πfct)[10]，fc為載頻，則旋翼目標單個散射點的基帶回波信號可表示為

(1)

其中，σn，i為第n個葉片上第i個散射點的散射系數，rn，i(t)為第n個葉片上第i個散射點到雷達的距離。

由于li?R，則葉片散射點與雷達間的距離rn，i可由菲涅爾近似寫為[11]

(2)

其中

(3)

式中，θ0為該葉片初始旋轉角，frot為葉片旋轉頻率，單位為r/s。

假設旋翼目標上有N個葉片，每個葉片上有K個散射點，則其回波信號可表示為

(4)

其中，exp(-j4πlicosθ(t)/λ)為多普勒項，包含了旋翼目標的多普勒信息。

在均勻散射前提下，每個散射點的散射系數可看作是相等的，散射點間隔d=L/K-1，則li=d(i-1)，式(4)可作如下變換：

(5)

當散射點間隔d趨近于0時，由近似可得sin(2πdcosθ(t)/λ)≈sin(2πdcosθ(t)/λ)，則式(5)可近似為

(6)

由式(6)可以看出，旋翼目標的回波幅度受辛克(sinc)函數調制，當cos(θ)=0時，旋翼目標回波才達到sinc函數的峰值，而此時旋翼目標葉片與雷達視線LOS的夾角θ為±π/2，所以只有當葉片與雷達視線LOS垂直時，葉片回波才達到最大值，這被稱作時域“閃爍”現象[12]，該時刻被稱作“閃爍”時刻。單個葉片在一個旋轉周期內可與LOS垂直兩次，因此會產生兩次“閃爍”現象。

同時，旋翼目標散射點的多普勒頻率可以表示為

(7)

由于各散射點的旋轉半徑不同，因此各散射點具有不同的線速度vi=ωili，所以散射點的多普勒頻率也是不同的。由式(7)可知旋翼目標單個葉片的回波頻率由各散射點的多普勒頻率組成

fd1(t)=[fd11(t)，fd12(t)，…，fd1N(t)]

(8)

此外，旋翼目標各散射點的回波多普勒頻率是一個隨時間變化的量，隨著旋轉角度的變化，單個散射點的多普勒頻率呈現出正弦變化趨勢，其周期與葉片旋轉周期保持一致。而對于某個t時刻而言，旋翼目標葉尖的散射點具有最大的多普勒頻率

(9)

而葉片旋轉中心附近的散射點多普勒頻率近乎為零。因此旋翼目標葉片回波的多普勒頻率實際上是一個較寬的頻率帶，其范圍為[-4πfrotL/λ，4πfrotL/λ]，當葉片轉向雷達方向時，其頻率為正，當葉片背離雷達方向時，其多普勒頻率為負，其葉片轉向與頻率正負關系如圖2所示。

圖2 多普勒頻率正負示意圖

3 基于時頻特征的微動參數提取原理

3.1 旋翼目標時頻特征分析

由式(8)可知，旋翼目標的多普勒頻率是隨時間變化的量，因此其回波實際是一種時變非平穩信號。而常規雷達信號處理是建立在平穩信號的基礎上的，其頻域分析通常常采用傅里葉變換，只能從整體上表示信號的頻率構成，無法反映出信號的頻譜特性隨時間變化的情況。所以頻域分析方法不再適用于旋翼葉片回波這一非平穩時變信號，需要利用時頻域聯合分析處理的方法來對回波信號進行分析[13]，本文采用Garber-短時間傅里葉變換(Garber-STFT)作為旋翼回波時頻聯合分析方法。

短時傅里葉變換是通過對信號加滑動的時間窗，在時間窗內的信號可看作是平穩信號，因此將長的非平穩信號等價于一系列短平穩信號的疊加，通過時間窗內的信號進行頻域變換表示該時刻內信號的頻譜特征，由此得到了非平穩信號的時頻聯合分布，其表達式為：

(10)

式中g(t)為窗函數，其中u為窗函數的長度.當g(t)為高斯窗函數時，該短時傅里葉變換就是Garber變換，STFT時頻分析過程可以如圖3所示。

圖3 STFT變換示意圖

當高斯時間窗滑動到以時域閃爍為中心的信號時，回波信號顯示為一個相對完整的辛格峰值包絡形式，根據傅里葉變換關系，此時該段信號的頻域為矩形包絡rect(f)形式，即旋翼目標回波在“閃爍”時刻頻率為矩形頻帶，這被稱作旋翼葉片回波在時頻域的“閃爍效應”，與時域“閃爍”具有一致性。而當高斯時間窗滑動在非閃爍時刻時，此時回波主要是由葉尖和旋轉中心的散射點反射的強電磁散射回波組成，所以在這些時刻的微多普勒特征主要表現為由葉尖引起的正弦曲線包絡形式和旋轉中心散射點引起的零頻帶組成。

綜上，旋翼目標回波時頻特征具有以下特點：

1) 旋翼目標回波時頻特征主要由三部分組成：由旋轉中心引起的零頻帶、葉尖引起的正弦包絡和以及回波峰值產生的時頻“閃爍”。

2) 時頻“閃爍”與旋翼目標旋轉周期具有對應關系。在一個葉片旋轉周期內時頻域“閃爍”兩次，且只在葉片與LOS垂直時發生。其中，當葉片由背離方向轉向雷達方向時，為“正閃爍”，其頻率區間為[0，fmax]，由靠近轉向背離雷達時，為“負閃爍”，其頻率區間為[-fmax，0]。

3) 由于旋翼目標回波峰值積累了大量回波能量，相應地，其時頻“閃爍”也具有很高的能量。

3.2 旋翼目標葉片旋轉頻率提取

由于旋翼目標回波的“時頻”閃爍現象與其旋轉周期具有對應關系，因此可以通過計算相鄰閃爍頻帶間的時間間隔來估計旋翼目標的旋轉頻率。以三葉片的旋翼目標為例，由于相鄰兩個葉片間隔120o，所以對于時頻“閃爍”而言，相鄰兩個同向“閃爍”間的時間間隔為1/3個旋轉周期，假設同向閃爍時刻間的間隔Δt，那么旋翼目標葉片的轉速估計值為

(11)

圖4 葉片轉動與閃爍時刻對應關系示意圖

對同向閃爍時刻間的間隔Δt的求解，可以利用圖像處理技術手段[15]對時頻圖進行預處理以提取Δt。由于回波中除旋翼目標外往往還包含大量的噪聲，因此首先需要對時頻圖像進行降噪處理。其次，對降噪后的時頻圖進行二值化處理，提取時頻“閃爍”特征，由于時頻“閃爍”回波能量較強，因此通過設置能量門限η實現對時頻圖的二值化處理，剔除與時頻“閃爍”無關的圖像信息。同時由于二值化結果較為粗糙，還需要進一步作平滑處理。

旋翼目標回波時頻“閃爍”頻率帶受STFT分辨率的影響具有展寬效應，其閃爍頻帶往往具有一定的“寬度”，因此在橫向占據多個時頻單元，該展寬效應可能會造成Δt估計誤差的產生，所以對每個頻帶寬度都選取該橫向寬度的中間時頻單元為計算單元，這樣有利于消除展寬效應，提高葉片旋轉頻率的估計精度。

3.3 葉片初始相位提取

基于時頻特征的旋翼目標微動參數提取方法不僅可以提取葉片的旋轉頻率，同時還可以根據閃爍頻帶出現的位置估計葉片的初始相位。假設初始時刻t0的葉片初始相位為θ0，其到第一“閃爍時刻”為t1，此時葉片轉過了θ1，其中θ1=2πfrott1，其角度關系如圖5所示。因此以“閃爍”作為參照點，利用幾何知識可得葉片的初始相位估計值為

(12)

圖5 葉片旋轉角度關系示意圖

(13)

由式(8)可知，旋翼目標同一葉片上的散射點具有相同的初始相位，不同葉片上的散射點初始相位不同，相鄰連個葉片間的初始相位相差2π/3，其它葉片的初始相位可在一個葉片的初始相位估計值0基礎上加上相應的葉片角度間隔即可。

3.4 旋翼葉片長度提取

由式(9)回波頻率的最大值可知，其回波最大頻率與旋翼的轉速和葉片長度有關，因此旋翼葉片的長度可由式提取

(14)

fdmax可由時頻圖中頻率范圍直接獲取，利用式(11)可以得到旋翼轉速的估計值，

同時，基于時頻特征的旋翼，目標微動參數提取方法也可適用于多個旋翼目標的場景，如對風電場中多個旋翼目標微動參數的提取。由于各旋翼目標間旋轉頻率均不相同，因此其所有散射點中的最大的多普勒頻率也不盡相同，也就是說其時頻“閃爍”的頻帶范圍不同，所以可以通過頻帶范圍區分不同的旋翼目標，再根據每個旋翼目標的同向閃爍時刻間的間隔即可計算其葉片旋轉頻率。

4 仿真與分析

為了驗證本文方法對旋翼目標回波的微動參數提取方法的有效性，以三葉片風輪機為旋翼目標進行仿真，雷達仿真參數為：脈沖重復頻率為1000Hz；載波頻率fc=1GHz；觀測時間為2s；雷達波束與風輪機葉片旋轉中心距離為20km；風輪機葉片長度均為L=20m。同時雷達波束在風輪機葉片上發生均勻散射，且葉片上散射點間隔均為d=λ/10，所加噪聲為高斯白噪聲。

仿真1：單風輪機回波微動參數提取

假設某風輪機葉片旋轉頻率frot=0.4976Hz，葉片初始相位為θ0=1.2129rad，回波信噪比SNR為30dB。其回波微動特征如圖6所示。

圖6 單風輪機回波

由圖6可以看出旋翼目標回波呈現辛格函數調制，回波幅度周期性變化，同時其時頻域特征由零頻帶、時頻閃爍和正弦包絡三部分組成。因此通過回波時域特征和時頻域特征可以直觀看出風輪機回波的周期性，但是難以直接精準的提取其葉片旋轉頻率和葉片初始相位。在時頻域對回波進行降噪、平滑處理和二值化處理，其時頻閃爍提取結果如圖7所示，其中二值化能量門限η=20dB。

圖7 時頻“閃爍”提取結果

利用圖7中閃爍提取結果，通過尋峰函數可以得到時頻“閃爍”負向頻帶的閃爍時刻分別為t1=0.1140s，t2=0.7835s，t3=1.4535s，則葉片旋轉頻率估計值為rot=0.4978Hz，與仿真設置值的誤差為0.04%，同時，由式(12)可以求得葉片初始相位0=1.2143rad，與仿真設置值的誤差為0.12%，相應地，由時頻圖可知，fdmax=416.87Hz，代入式(14)可知其葉片長度估計值為=19.9920m，其誤差為0.4%。由以上結果可知，該微動參數估計結果具有較高的精度，誤差均在可接受范圍內。

仿真2：多風輪機回波微動參數提取

當多個風輪機分別位于不同距離單元時，可分別對不同距離單元逐個進行微動參數提取，其本質與單個風輪機回波微動參數提取無異。本仿真主要考慮的是在同一距離單元內具有多個風輪機回波，風輪機間的參數均不相同。因此假設距離雷達20km附近有一風電場，場內有三臺型號相同風輪機，其葉片轉速分別為fr1=0.5282Hz，fr2=0.5025Hz，fr3=0.4976Hz，其葉片初相對應分別為θ01=0.1148rad，θ02=1.352rad，θ03=0.7253rad，其它參數與仿真1保持不變，不考慮風輪機葉間的相互影響，其回波時域與時頻域微動特征如圖8所示。

圖8 多風輪機回波

通過圖8可知，由于三臺風輪機葉片旋轉頻率各不相同，多風輪機回波時域和時頻域特征不同于單個風輪機，其回波周期特性難以通過時域直接體現出來，文獻[9]所提利用相關函數提取方法失效。但由于不同風輪機葉片旋轉頻率不同，因此其時頻域“閃爍”頻帶范圍也不相同，通過閃爍頻帶范圍(時頻圖中時頻閃爍的長度)即可區分不同風輪機的回波。對多風輪機回波時頻圖提取時頻“閃爍”結果如圖9所示。

圖9 多風輪機回波時頻閃爍提取結果

通過時頻提取結果結果，可以區分出風輪機1的負閃爍時刻為0.4380s、1.0695s、1.7000s，風輪機2的負閃爍時刻為0.0685s、0.7320s、1.3955s，風輪機3的負閃爍時刻為0.2700s、0.9390s、1.6125s。多風輪機微動參數提取結果如表1所示。

由表1提取結果可以看出，本文所提基于時頻特征的微動參數提取方法對多個風輪機的多個旋轉頻率和初始相位實現有效提取，其中旋轉頻率誤差小于0.5%，初始相位誤差小于2%，葉片初相估計誤差明顯大于旋轉頻率誤差，這是因為葉片初相需要利用旋轉頻率估計值求解而造成了誤差積累，但是該誤差在可接受范圍內。因此本文所提方法可適用于多個風輪機場景下的微動參數提取。

表1 多風輪機微動參數提取結果

仿真3：不同SNR下估計誤差分析

為驗證在不同信噪比(SNR)下本文方法提取微動特征的魯棒性，本文進行M次蒙特卡羅實驗，并進行均方根誤差(Root Mean Square Error，RMSE)分析。

對于估計的旋轉頻率，其RMSE的表達式為

(15)

本文取M=100，設置SNR的范圍為2～40dB，步進2dB，圖10為不同SNR場景下估計結果的RMSE。

圖10 估計誤差結果

由圖10(a)可知，估計旋轉頻率的RMSE隨著SNR的提高而降低，當SNR小于6dB時，該方法估計的葉片旋轉頻率誤差較大，這是由于噪聲功率過大無法通過時頻圖像提取有效信息。當SNR達到6dB以上時，RMSE的值小于1%，說明該方法在高信雜比條件下具有很高的精度，此時RMSE趨于穩定，說明該算法具有一定的穩健性。

圖10(b)為葉片初始相位估計值的RMSE結果，由圖可知，初始相位的估計值同樣也隨著SNR的提高而降低，在6dB以下誤差精度較大，得不到初始相位的有效信息，此時葉片初始相位RMSE結果大于旋轉頻率的RMSE，也是由于噪聲功率過高，在時頻域無法提取到有效信息，同時，葉片旋轉頻率誤差過大在初始相位估計時也造成了誤差積累。當SNR在6dB以上時，初始相位估計誤差區域穩定且小于2%，此時算法具有較好的參數估計精度。

5 結論

旋翼目標的微動參數在目標探測、識別領域具有重要意義。本文提出了一種基于時頻特征的雷達旋翼目標微動參數估計方法。該方法主要通過分析旋翼目標回波時頻特征，利用時頻域“閃爍”特征與微動參數間的對應關系，在時頻圖像域實現了對微動參數的有效估計。仿真結果表明該估計方法不僅可以實現對單個旋翼目標的葉片旋轉頻率、葉片初始相位和葉片長度的高效、快速提取，而且在多旋翼目標場景下也同樣適用。同時該估計方法在高信雜比條件下，估計結果具有良好的精度且算法具有一定的穩定性。