馮維婷 梁 青

（西安郵電大學電子工程學院，陜西西安 710121）

1 引言

雷達目標中的風輪機［1］、直升機［2］，無人機［3］是一類有旋轉(zhuǎn)葉片的目標，此類目標的雷達回波信號中會包含旋轉(zhuǎn)葉片的長度、轉(zhuǎn)速，初始旋轉(zhuǎn)角信息，準確地提取出這些參數(shù)是識別目標類型的重要依據(jù)［4-6］。因此，對目標旋轉(zhuǎn)葉片參數(shù)的提取分析是雷達信號處理領域的研究熱點［7-8］。

旋轉(zhuǎn)葉片的運動屬于微動［9］，目標微動會在雷達回波信號中產(chǎn)生額外地頻率調(diào)制稱為微多普勒頻率［10-11］。對葉片微動參數(shù)的估計問題，大量的文獻基于時頻譜圖處理。主要有兩種處理方法：一種處理是在時頻譜圖中搜索局部峰值點得到時頻脊線［12］，其中蘊含了目標的微多普勒頻率，從最大微多普勒頻率中反演出旋轉(zhuǎn)葉片的微動參數(shù)［13-14］；一種處理是對時頻譜圖進行逆Radon 變換［15-16］、Hough變換［17］等將微動參數(shù)估計問題轉(zhuǎn)化為變換后的參數(shù)空間的峰值檢測問題。文獻［13-14］的實測數(shù)據(jù)表明，強雜波環(huán)境中的時頻譜圖模糊，微多普勒頻率有較大波動，采用第一種時頻譜圖處理方法的微動參數(shù)估值存在較大波動，影響了參數(shù)估計精度；而第二種處理方法在雷達回波信號存在強閃爍［18-19］、正弦曲線特征不明顯時方法失效，無法估計微動參數(shù)。

針對強閃爍現(xiàn)象下旋轉(zhuǎn)葉片的參數(shù)估計問題，文獻［20］利用DnCNN網(wǎng)絡結(jié)構訓練去閃爍網(wǎng)絡，消除了時頻譜圖中的“閃爍”條，得到了正弦特征增強的旋轉(zhuǎn)葉片時頻譜圖，然后基于逆Radon 變換估計出微動參數(shù)。然而實際情況中，目標所處的噪聲、干擾環(huán)境復雜，該方法性能嚴重下降甚至失效。文獻［21］利用相位補償方法將相同葉片的“閃爍”條聚焦到特定多普勒頻率，實現(xiàn)了微動參數(shù)的估計，但情況分類較多，方法過于復雜。

針對強閃爍情況下旋轉(zhuǎn)葉片微動參數(shù)估計的問題，本文從調(diào)頻連續(xù)波（frequency modulated con?tinuous wave，F(xiàn)MCW）雷達回波信號的特點出發(fā)，提出了微動補償方法完成旋轉(zhuǎn)葉片參數(shù)估計。建立了FMCW 雷達的旋轉(zhuǎn)葉片基帶回波模型，在閃爍現(xiàn)象下回波信號具有周期性脈沖特征的基礎上，構造了基于周期脈沖模型的微動參數(shù)補償算子對回波信號進行參數(shù)補償，搜索補償后信號能量峰值完成了旋轉(zhuǎn)葉片微動參數(shù)的估計。利用FMCW 毫米波小型雷達系統(tǒng)IWR6843 搭建實驗平臺進行實測數(shù)據(jù)采集，仿真和實測數(shù)據(jù)的實驗結(jié)果表明，所提方法在強閃爍、干擾和噪聲存在的實際情況下可得到較高精度的參數(shù)估計值。

2 旋轉(zhuǎn)葉片雷達回波分析

FMCW 雷達系統(tǒng)通過獲得旋轉(zhuǎn)葉片回波可進行葉片相關參數(shù)估計。葉片上某散射點P與雷達的空間幾何關系如圖1所示。

圖1 雷達與葉片幾何關系示意圖Fig.1 Diagram of geometrical relationship between radar and blades

雷達波束中心O與旋轉(zhuǎn)葉片中心O′之間距離為R0，散射點P到O′的距離為rP，散射點P所在葉片與OO′的夾角為φ(t)，則散射點P到雷達的瞬時距離為RP(t)=圖1 中雷達坐標系XYZ與目標坐標系xyz的各坐標軸均平行，葉片在xz平面旋轉(zhuǎn)，散射點P與x軸的初始夾角為θ0，t時刻的夾角θ(t)=2πfrt+θ0，其中fr為旋轉(zhuǎn)頻率；圖中α是方位角，β是俯仰角。

其中，cosφ(t)=cosαcosβcosθ(t)+sinβsinθ(t)。雷達發(fā)射LFM 信號，中心頻率fc對應的載波波長為λ，散射點P的基頻回波信號表示為［22］

假設目標已經(jīng)定位，且完成平動補償，去掉式（2）中RP(t)的第一項R0部分，剩下目標微動部分，目標微動回波信號表示為

式（3）中回波信號相位的導數(shù)就是微多普勒頻率，表示為

由式（4）可見，葉片上散射點P的微多普勒頻率fP(t)是正弦曲線形式，正弦曲線的頻率是葉片旋轉(zhuǎn)頻率fr、幅度和初相大小與目標的位置α、β角有關，還與葉片的參數(shù)rP、fr有關。葉片長度為r，在r上對式（3）進行線積分，則單葉片的回波信號表示為

根據(jù)信號處理理論，sinc(?)與rect(?)是一對傅里葉變換對，即sinc(t) ?FTrect(f)，所以，式（6）中sinc項的微多普勒頻率為周期性矩形脈沖信號。式（6）中指數(shù)項的微多普勒頻率是正弦信號；時域上指數(shù)項與sinc 項的乘積，在頻域是兩者微多普勒頻率卷積的結(jié)果，故時頻譜圖中出現(xiàn)周期性矩形脈沖信號，被稱為閃爍條［23］。

采用短時傅里葉變換（short time Fourier trans?form，STFT）［24］的時頻分析方法得到2 葉片和3 葉片的回波信號時頻譜，如圖2（a）、（b）所示。圖2（a）中偶數(shù)葉片的微多普勒頻率表現(xiàn)為正負頻率關于時間對稱的閃爍條。圖2（b）中奇數(shù)葉片的微多普勒頻率表現(xiàn)為正負頻率等間隔交錯的閃爍條。

圖2 基于STFT的時頻譜Fig.2 Time-frequency spectrum based on STFT

圖2（a）中，葉片數(shù)為2，每轉(zhuǎn)動π rad 一對葉片就都與雷達視線垂直，因此在一個轉(zhuǎn)動周期中包含2 個正負頻率對稱的閃爍條。圖2（b）中，葉片數(shù)為3，每轉(zhuǎn)動π/3 rad，就會有一個迎著或背著雷達的葉片交替著與雷達視線垂直，所以一個周期中包含6個正負頻率交錯的閃爍條。

另外，根據(jù)式（4），當俯仰角β=0時，即雷達視線在XOY平面時，葉尖處最大微多普勒頻率為fmax=cosα；當方位角α=0 時，fmax=由以上分析可知，最大微多普勒頻率值由參數(shù)r、fr、α，β和λ共同決定，蘊含了葉片的尺寸和運動信息。文獻［13-14］中的基于時頻譜圖的參數(shù)提取方法原理就是根據(jù)α，β和λ確定已知時，最大微多普勒頻率值fmax由微動參數(shù)r和fr決定，故基于時頻譜圖測得最大微多普勒頻率值fmax后再反演出旋轉(zhuǎn)葉片的微動參數(shù)r和fr。文獻［13-14］中具體的基于時頻譜圖的參數(shù)提取方法是先對時頻譜圖的微多普勒曲線做傅里葉變換，得到旋轉(zhuǎn)頻率fr的估計，再在時頻譜圖中估計最大微多普勒頻率fmax，得到fr和fmax這兩個參數(shù)之后，根據(jù)r=可以推算出旋轉(zhuǎn)長度r。

3 旋轉(zhuǎn)葉片的參數(shù)估計方法

旋轉(zhuǎn)葉片的參數(shù)主要是：旋轉(zhuǎn)長度r、旋轉(zhuǎn)頻率fr和初始旋轉(zhuǎn)角θ0。對于旋轉(zhuǎn)頻率fr的估計，常用的算法有自相關、奇異值分解方法，以及對微多普勒曲線進行頻譜分析方法，實際的雷達回波信號存在目標能量微弱或正負頻譜能量不對稱的情況，這些方法不適用于實際復雜情況。對于旋轉(zhuǎn)長度r的估計，目前常采用的方法是基于時頻譜圖直接估計最大微多普勒頻率fmax，然后推算出旋轉(zhuǎn)長度。然而，時頻譜圖直接估計最大微多普勒頻率的誤差較大，這是由于旋轉(zhuǎn)葉片形狀、雷達的位置，以及各種干擾和噪聲都會影響雷達回波時頻譜圖的清晰度，導致這種方法估計得到的r誤差較大?？紤]到實際情況，為了有效進行微動參數(shù)的估計，結(jié)合雷達回波信號模型特征的基礎上，本節(jié)采用微動參數(shù)補償?shù)姆椒▽崿F(xiàn)旋轉(zhuǎn)葉片微動參數(shù)的估計。

從圖2 的旋轉(zhuǎn)葉片時頻譜圖可看出，回波信號的能量分布在閃爍條和正弦曲線上。如果構造補償信號進行微動補償，當補償完全時，補償后信號能量集中于0 頻處、達到最大；當補償不完全時，補償后信號能量仍分散在時頻譜圖的不同頻率處，小于完全補償時的能量。因此，可進行補償后信號能量峰值搜索，最大值處對應完全補償，基于此來估計旋轉(zhuǎn)葉片的微動參數(shù)。

現(xiàn)對微動補償原理進行說明，為簡化表述，假設α=0、β=0，此時將式（5）的單葉片回波信號重寫為

其中，A為信號振幅。根據(jù)式（7）的回波信號模型，構造補償算子為

其中，Pi={ri，fri，θ0i}為待估計的微動參數(shù)集合。將式（7）與式（8）相乘得到補償后信號

設定參數(shù)集合中各參數(shù)取值范圍，搜索補償后信號能量和的最大值，即

當完全補償時，補償后信號sb(t)=A，式（10）達到最大值。此時的參數(shù)集就是待估計參數(shù)的最優(yōu)估計值。故通過式（10）求得的便是微動參數(shù)的最優(yōu)估計值。

式（8）的補償算子中有除法運算，當分母取值較小時受外在干擾影響補償后的信號中會出現(xiàn)尖峰脈沖噪聲，會影響參數(shù)估計精度。為了抑制此類噪聲，將式（9）修正為

式（11）中進行了非線性變換，對原補償后的信號除以自身的幅值，將信號中的尖峰脈沖噪聲幅值衰減到1，可抑制脈沖噪聲。

總回波信號模型是多目標多葉片的多分量信號時，可以采用迭代方法得到各分量信號的參數(shù)。先利用所提方法估計出最強分量的參數(shù)，接著構造出最強分量信號，從總回波中剔除掉該分量信號，重復該過程直到估計出所有信號分量參數(shù)。微動補償參數(shù)估計算法具體步驟如下：

步驟1設置參數(shù)集Pi={ri，fri，θ0i}中各參數(shù)的取值范圍、搜索步長；多分量信號迭代進行初始化：k=1。

步驟2對雷達基帶回波信號sint(t)乘以補償算子，見式（8）；再根據(jù)式（11）修正補償后信號。

步驟3在參數(shù)設定范圍內(nèi)計算補償后信號的能量，在結(jié)果中搜索峰值得到參數(shù)估計值，見式（10）。

步驟4利用所估計的參數(shù)重構單分量信號，并從原回波信號sint(t)中剔除該分量信號。

步驟5判斷是否滿足條件k≥Q，不滿足時迭代：k=k+1，重復步驟2～4；直到滿足條件停止。

旋轉(zhuǎn)葉片微動參數(shù)估計算法的處理流程如圖3所示。

圖3 算法流程圖Fig.3 Flow chart of the proposed algorithm

4 仿真和實測實驗結(jié)果

4.1 仿真分析

設置雷達波長為0.005 m，脈沖重復頻率PRF為5 kHz。兩個旋轉(zhuǎn)目標，俯仰角β都為0，方位角α分別為0、π/6 rad，目標主體平動已經(jīng)完全補償，第1個目標2葉片，相鄰葉片相隔π rad；第2個目標3葉片，相鄰葉片相隔2π/3 rad；采用式（6）的基帶回波模型，觀測時長0.5 s。目標參數(shù)具體設置如表1所示。

表1 仿真目標的參數(shù)Tab.1 Parameters of the simulated targets

旋轉(zhuǎn)目標基于STFT 的時頻譜如圖4 所示。圖中不同目標的微多普勒曲線具有不同的最大多普勒頻率值，部分閃爍條重合；目標1 的2 葉閃爍條是正負頻率成對出現(xiàn)，而目標2 的3 葉閃爍條是正負頻率等間隔交錯出現(xiàn)。

圖4 基于STFT的時頻譜Fig.4 Time-frequency spectrum based on STFT

利用所提的方法估計微動參數(shù)，設置各參數(shù)搜索范圍，構造微動補償算子，頻率搜索范圍為1.0～6.0 Hz，步長為0.1 Hz；旋轉(zhuǎn)長度范圍0～0.30 m，步長0.01 m；初相角范圍0～6.27 rad，步長0.01 rad。所提方法的參數(shù)估計結(jié)果如圖5所示。圖5（a）中兩處峰值對應的旋轉(zhuǎn)頻率估值分別是4 Hz、3 Hz，對應于目標1、2；圖5（b）是目標1 的旋轉(zhuǎn)長度r和初相θ0估計結(jié)果m，兩葉片的初相分別為=0 rad、=3.14 rad，與目標1的參數(shù)值相同；圖5（c）中峰值處對應目標2的r和θ0估計結(jié)果22 m，三葉片的初相分別為1.57 rad、3.66 rad5.76 rad，與目標2的參數(shù)值相同。由圖5可見，各目標的參數(shù)估計值與真值一致，表明了所提方法的有效性。

圖5 仿真目標參數(shù)估計結(jié)果Fig.5 Parameter estimation results for the simulated targets

為說明所提方法的有效性和對不同目標參數(shù)的適應性，設置與4.1節(jié)相同的仿真條件，利用其他方法進行信號處理。文獻［13］提出的基于時頻譜提取最大多普勒頻率方法對參數(shù)估計的結(jié)果如圖6。圖6中可以看出目標1的正最大多普勒頻率約為949.4 Hz，目標2的負最大多普勒頻率為?1616 Hz，在目標1、2的旋轉(zhuǎn)頻率估值分別是4 Hz、3 Hz時，反演推算出目標旋轉(zhuǎn)長度的估值分別為0.094 m、0.214 m。該方法結(jié)果誤差大于所提方法，這是由于時頻譜圖的測量結(jié)果僅能給出粗略的最大多普勒頻率值，在噪聲、干擾環(huán)境下測量誤差更大，故基于時頻譜圖的最大多普勒頻率提取方法不利于進行微動參數(shù)精確估計。

圖6 時頻譜圖方法的結(jié)果Fig.6 The results of time-frequency spectrum method

圖7（a）、（b）是旋轉(zhuǎn)頻率估值分別為4 Hz、3 Hz時閃爍現(xiàn)象下采用文獻［20］逆Radon變換方法的結(jié)果，可以看出，整個圖中無明顯的聚焦點，因此閃爍條件下使用逆Radon 變換是無法估計出目標微動參數(shù)。圖7（c）～（e）是理想的去除閃爍后時頻譜圖及其逆Radon 變換的結(jié)果，圖7（d）、（e）中有明顯的聚焦點，檢測出峰值點即可得到微動參數(shù)估值。然而實測環(huán)境中強噪聲、干擾存在時，去除閃爍的處理能力十分有限，難以完全去除閃爍，此情況下如同圖7（a）、（b），此時逆Radon變換目標微動的參數(shù)方法失效。

采用本文所提方法在信噪比（signal to noise ratio，SNR）分別為0 dB、3 dB、6 dB，9 dB 時各進行100 點蒙特卡羅實驗，得到旋轉(zhuǎn)目標的參數(shù)平均估計值見表2 到表4。由表可以看出，旋轉(zhuǎn)長度r和頻率fr在SNR 為0 dB 時仍有較高估計精度，初始相位估值在SNR 為0 dB時誤差較大；但隨著SNR 的增加各參數(shù)估計誤差都減小，當SNR 為9 dB 時，初始相位估值也達到較高估計精度。

表2 旋轉(zhuǎn)長度估計值Tab.2 Parameter estimation of rotating length

表3 旋轉(zhuǎn)頻率估計值Tab.3 Parameter estimation of rotating frequency

表4 初始相位估計值Tab.4 Parameter estimation of initial phase

4.2 實測實驗

實測實驗使用的FMCW 雷達是德州儀器公司的IWR6843 毫米波雷達系統(tǒng)，采用1 發(fā)4 收模式采集信號。雷達接收到的回波信號與發(fā)射信號經(jīng)過混頻處理和低通濾波后得到中頻信號。

實驗對象為一旋轉(zhuǎn)的3 葉風扇，雷達與風扇之間的距離約2.4 m，風扇轉(zhuǎn)速為330 r/min，即旋轉(zhuǎn)頻率為5.5 Hz，葉片長度為0.21 m。

實驗中設置雷達參數(shù)如下：雷達波長λ=4.9 mm，中心頻率fc=61.56 GHz，發(fā)射LFM 信號的調(diào)頻斜率K=40.99 MHz/us，實際帶寬B=2.62 GHz；一幀數(shù)據(jù)包含128個調(diào)頻周期數(shù)，脈沖重復頻率為6.4 kHz，一個調(diào)頻周期內(nèi)的采樣點數(shù)為256點，采樣頻率為4 MHz。方位角α=30°，俯仰角β=0。

為了獲取旋轉(zhuǎn)葉片參數(shù)，首先對回波數(shù)據(jù)中每個調(diào)頻周期內(nèi)的快時間信號進行脈沖壓縮得到一維距離像，再對每個距離單元的慢時間信號進行傅里葉變換就得到了距離-多普勒像。一幀實測數(shù)據(jù)的距離-多普勒像如圖8 所示。目標在第45 距離單元，在該距離單元有多普勒展寬現(xiàn)象，目標主體的能量遠大于旋轉(zhuǎn)葉片的能量，同時較小距離單元附近存在雷達系統(tǒng)本身帶來的強靜止干擾。

圖8 一幀數(shù)據(jù)距離-多普勒像Fig.8 Range-Dopple imaging of one frame’s data

脈沖壓縮后的信號，僅提取目標所在距離單元的信號，得到了包含目標多普勒信息的多普勒信號。圖9是多普勒信號進行STFT變換后的時頻譜圖。

圖9 STFT時頻譜Fig.9 Time-frequency spectrum based on STFT

圖9（a）是包含了目標主體分量的時頻譜，零多普勒頻率對應的是目標主體分量和靜止干擾，由于主體分量和干擾的能量遠遠大于旋轉(zhuǎn)葉片分量的能量，圖中無法觀察到旋轉(zhuǎn)葉片分量的微多普勒頻率。對原多普勒信號濾波處理，圖9（b）是目標主體分量和靜止干擾濾除后的時頻譜，可觀察出正負頻率交錯且等間隔周期性的閃爍條，這與三葉片時頻特性的理論分析結(jié)果相符合。

采用所提微動補償方法進行旋轉(zhuǎn)葉片的參數(shù)估計結(jié)果如圖10 所示。圖10（a）中旋轉(zhuǎn)頻率估值5.6 Hz，定義參數(shù)fr估值的相對誤差為er=則頻率估值相對誤差為0.02。圖10（b）～（d）中葉片長度的估值21 m，相對誤差為0，各分量的初相估值3.39 rad、5.54 rad1.37 rad。

圖10 真實目標參數(shù)估計結(jié)果Fig.10 The real target’s parameter estimation results

旋轉(zhuǎn)葉片的實際初始相位未知，但是三葉相鄰葉片相位差是確定的，3 葉片分量的初始相位和相位差估計值見表5。表5 中相位差的取值范圍設置在0～π 之間，相位差的估值與真值基本符合。因此，所提方法能有效估計實際雷達旋轉(zhuǎn)葉片的旋轉(zhuǎn)頻率、旋轉(zhuǎn)長度和初始相位值。

表5 實測數(shù)據(jù)初相相位差估計值Tab.5 Parameter estimation of initial phase difference for experimental data

5 結(jié)論

旋轉(zhuǎn)葉片的微動參數(shù)提取對目標類型的識別有重要意義，然而閃爍現(xiàn)象下的微動參數(shù)提取方法近幾年才得到關注研究。本文推導了強閃爍情況下旋轉(zhuǎn)葉片的回波信號模型，基于該信號模型提出了微動補償方法，最終實現(xiàn)了閃爍下的旋轉(zhuǎn)葉片微動參數(shù)估計。仿真分析和FMCW 雷達系統(tǒng)實測實驗結(jié)果表明，所提方法能有效提取旋轉(zhuǎn)葉片旋轉(zhuǎn)頻率、長度和初始相位值，參數(shù)估計精度高，為復雜背景中雷達旋轉(zhuǎn)目標微動參數(shù)獲取提供了一種途徑。