基于時頻特征的旋翼目標(biāo)微動參數(shù)估計方法

宋村夫，曲智國，胡旭超

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

1 引言

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

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

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

2 旋翼目標(biāo)雷達(dá)回波建模與分析

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

圖1 旋翼目標(biāo)回波模型

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

(1)

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

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

(2)

其中

(3)

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

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

(4)

其中，exp(-j4πl(wèi)icosθ(t)/λ)為多普勒項，包含了旋翼目標(biāo)的多普勒信息。

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

(5)

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

(6)

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

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

(7)

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

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

(8)

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

(9)

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

圖2 多普勒頻率正負(fù)示意圖

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

3.1 旋翼目標(biāo)時頻特征分析

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

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

(10)

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

圖3 STFT變換示意圖

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

綜上，旋翼目標(biāo)回波時頻特征具有以下特點：

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

2) 時頻“閃爍”與旋翼目標(biāo)旋轉(zhuǎn)周期具有對應(yīng)關(guān)系。在一個葉片旋轉(zhuǎn)周期內(nèi)時頻域“閃爍”兩次，且只在葉片與LOS垂直時發(fā)生。其中，當(dāng)葉片由背離方向轉(zhuǎn)向雷達(dá)方向時，為“正閃爍”，其頻率區(qū)間為[0，fmax]，由靠近轉(zhuǎn)向背離雷達(dá)時，為“負(fù)閃爍”，其頻率區(qū)間為[-fmax，0]。

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

3.2 旋翼目標(biāo)葉片旋轉(zhuǎn)頻率提取

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

(11)

圖4 葉片轉(zhuǎn)動與閃爍時刻對應(yīng)關(guān)系示意圖

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

旋翼目標(biāo)回波時頻“閃爍”頻率帶受STFT分辨率的影響具有展寬效應(yīng)，其閃爍頻帶往往具有一定的“寬度”，因此在橫向占據(jù)多個時頻單元，該展寬效應(yīng)可能會造成Δt估計誤差的產(chǎn)生，所以對每個頻帶寬度都選取該橫向?qū)挾鹊闹虚g時頻單元為計算單元，這樣有利于消除展寬效應(yīng)，提高葉片旋轉(zhuǎn)頻率的估計精度。

3.3 葉片初始相位提取

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

(12)

圖5 葉片旋轉(zhuǎn)角度關(guān)系示意圖

(13)

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

3.4 旋翼葉片長度提取

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

(14)

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

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

4 仿真與分析

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

仿真1：單風(fēng)輪機(jī)回波微動參數(shù)提取

假設(shè)某風(fēng)輪機(jī)葉片旋轉(zhuǎn)頻率frot=0.4976Hz，葉片初始相位為θ0=1.2129rad，回波信噪比SNR為30dB。其回波微動特征如圖6所示。

圖6 單風(fēng)輪機(jī)回波

由圖6可以看出旋翼目標(biāo)回波呈現(xiàn)辛格函數(shù)調(diào)制，回波幅度周期性變化，同時其時頻域特征由零頻帶、時頻閃爍和正弦包絡(luò)三部分組成。因此通過回波時域特征和時頻域特征可以直觀看出風(fēng)輪機(jī)回波的周期性，但是難以直接精準(zhǔn)的提取其葉片旋轉(zhuǎn)頻率和葉片初始相位。在時頻域?qū)夭ㄟM(jìn)行降噪、平滑處理和二值化處理，其時頻閃爍提取結(jié)果如圖7所示，其中二值化能量門限η=20dB。

圖7 時頻“閃爍”提取結(jié)果

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

仿真2：多風(fēng)輪機(jī)回波微動參數(shù)提取

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

圖8 多風(fēng)輪機(jī)回波

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

圖9 多風(fēng)輪機(jī)回波時頻閃爍提取結(jié)果

通過時頻提取結(jié)果結(jié)果，可以區(qū)分出風(fēng)輪機(jī)1的負(fù)閃爍時刻為0.4380s、1.0695s、1.7000s，風(fēng)輪機(jī)2的負(fù)閃爍時刻為0.0685s、0.7320s、1.3955s，風(fēng)輪機(jī)3的負(fù)閃爍時刻為0.2700s、0.9390s、1.6125s。多風(fēng)輪機(jī)微動參數(shù)提取結(jié)果如表1所示。

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

表1 多風(fēng)輪機(jī)微動參數(shù)提取結(jié)果

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

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

對于估計的旋轉(zhuǎn)頻率，其RMSE的表達(dá)式為

(15)

本文取M=100，設(shè)置SNR的范圍為2～40dB，步進(jìn)2dB，圖10為不同SNR場景下估計結(jié)果的RMSE。

圖10 估計誤差結(jié)果

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

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

5 結(jié)論

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