閃爍現(xiàn)象下基于微動補(bǔ)償?shù)男D(zhuǎn)葉片參數(shù)估計方法

馮維婷 梁 青

（西安郵電大學(xué)電子工程學(xué)院，陜西西安 710121）

1 引言

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

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

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

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

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

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

圖1 雷達(dá)與葉片幾何關(guān)系示意圖Fig.1 Diagram of geometrical relationship between radar and blades

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

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

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

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

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

根據(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)為正負(fù)頻率關(guān)于時間對稱的閃爍條。圖2（b）中奇數(shù)葉片的微多普勒頻率表現(xiàn)為正負(fù)頻率等間隔交錯的閃爍條。

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

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

另外，根據(jù)式（4），當(dāng)俯仰角β=0時，即雷達(dá)視線在XOY平面時，葉尖處最大微多普勒頻率為fmax=cosα；當(dāng)方位角α=0 時，fmax=由以上分析可知，最大微多普勒頻率值由參數(shù)r、fr、α，β和λ共同決定，蘊(yùn)含了葉片的尺寸和運(yùn)動信息。文獻(xiàn)［13-14］中的基于時頻譜圖的參數(shù)提取方法原理就是根據(jù)α，β和λ確定已知時，最大微多普勒頻率值fmax由微動參數(shù)r和fr決定，故基于時頻譜圖測得最大微多普勒頻率值fmax后再反演出旋轉(zhuǎn)葉片的微動參數(shù)r和fr。文獻(xiàn)［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的估計，常用的算法有自相關(guān)、奇異值分解方法，以及對微多普勒曲線進(jìn)行頻譜分析方法，實際的雷達(dá)回波信號存在目標(biāo)能量微弱或正負(fù)頻譜能量不對稱的情況，這些方法不適用于實際復(fù)雜情況。對于旋轉(zhuǎn)長度r的估計，目前常采用的方法是基于時頻譜圖直接估計最大微多普勒頻率fmax，然后推算出旋轉(zhuǎn)長度。然而，時頻譜圖直接估計最大微多普勒頻率的誤差較大，這是由于旋轉(zhuǎn)葉片形狀、雷達(dá)的位置，以及各種干擾和噪聲都會影響雷達(dá)回波時頻譜圖的清晰度，導(dǎo)致這種方法估計得到的r誤差較大?？紤]到實際情況，為了有效進(jìn)行微動參數(shù)的估計，結(jié)合雷達(dá)回波信號模型特征的基礎(chǔ)上，本節(jié)采用微動參數(shù)補(bǔ)償?shù)姆椒▽崿F(xiàn)旋轉(zhuǎn)葉片微動參數(shù)的估計。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4.1 仿真分析

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

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

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

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

利用所提的方法估計微動參數(shù)，設(shè)置各參數(shù)搜索范圍，構(gòu)造微動補(bǔ)償算子，頻率搜索范圍為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）中兩處峰值對應(yīng)的旋轉(zhuǎn)頻率估值分別是4 Hz、3 Hz，對應(yīng)于目標(biāo)1、2；圖5（b）是目標(biāo)1 的旋轉(zhuǎn)長度r和初相θ0估計結(jié)果m，兩葉片的初相分別為=0 rad、=3.14 rad，與目標(biāo)1的參數(shù)值相同；圖5（c）中峰值處對應(yīng)目標(biāo)2的r和θ0估計結(jié)果22 m，三葉片的初相分別為1.57 rad、3.66 rad5.76 rad，與目標(biāo)2的參數(shù)值相同。由圖5可見，各目標(biāo)的參數(shù)估計值與真值一致，表明了所提方法的有效性。

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

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

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

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

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

表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 雷達(dá)是德州儀器公司的IWR6843 毫米波雷達(dá)系統(tǒng)，采用1 發(fā)4 收模式采集信號。雷達(dá)接收到的回波信號與發(fā)射信號經(jīng)過混頻處理和低通濾波后得到中頻信號。

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

實驗中設(shè)置雷達(dá)參數(shù)如下：雷達(dá)波長λ=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ù)，脈沖重復(fù)頻率為6.4 kHz，一個調(diào)頻周期內(nèi)的采樣點數(shù)為256點，采樣頻率為4 MHz。方位角α=30°，俯仰角β=0。

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

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

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

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

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

采用所提微動補(bǔ)償方法進(jìn)行旋轉(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 真實目標(biāo)參數(shù)估計結(jié)果Fig.10 The real target’s parameter estimation results

旋轉(zhuǎn)葉片的實際初始相位未知，但是三葉相鄰葉片相位差是確定的，3 葉片分量的初始相位和相位差估計值見表5。表5 中相位差的取值范圍設(shè)置在0～π 之間，相位差的估值與真值基本符合。因此，所提方法能有效估計實際雷達(dá)旋轉(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ù)提取對目標(biāo)類型的識別有重要意義，然而閃爍現(xiàn)象下的微動參數(shù)提取方法近幾年才得到關(guān)注研究。本文推導(dǎo)了強(qiáng)閃爍情況下旋轉(zhuǎn)葉片的回波信號模型，基于該信號模型提出了微動補(bǔ)償方法，最終實現(xiàn)了閃爍下的旋轉(zhuǎn)葉片微動參數(shù)估計。仿真分析和FMCW 雷達(dá)系統(tǒng)實測實驗結(jié)果表明，所提方法能有效提取旋轉(zhuǎn)葉片旋轉(zhuǎn)頻率、長度和初始相位值，參數(shù)估計精度高，為復(fù)雜背景中雷達(dá)旋轉(zhuǎn)目標(biāo)微動參數(shù)獲取提供了一種途徑。