馬 嬌，董勇偉?，李 原，李凌霄，楊杰芳

(1 中國科學院電子學研究所 微波成像技術(shù)國家重點實驗室， 北京 100190； 2 中國科學院大學， 北京 100049)

近年來，小型無人機變得越來越普遍，在軍事領(lǐng)域和民用領(lǐng)域都得到廣泛應(yīng)用。但與此同時，由于其價格低廉、操作簡單，也被濫用于不安全、甚至犯罪的行為，對國家經(jīng)濟發(fā)展和國民安全構(gòu)成實際威脅[1]。因此，針對小型無人機目標的識別具有重要的應(yīng)用價值。無人機的獨特之處在于，幾乎所有無人機都有一個及以上數(shù)目的螺旋槳。旋翼的旋轉(zhuǎn)葉片會對雷達的后向散射信號產(chǎn)生周期性的頻率調(diào)制，在無人機自身多普勒譜的兩邊產(chǎn)生邊帶，這種附加的多普勒頻率調(diào)制被稱為微多普勒效應(yīng)[2-3]。

利用微多普勒特征識別小型無人機已成為研究熱點。文獻[3-5]研究不同型號無人機的微多普勒特征，提出基于頻譜圖的無人機分類方法，并通過實驗證明微多普勒特征在識別小型無人機時可以作為多普勒特征的一個補充特征，但文章未涉及無人機微多普勒參數(shù)的提取方法。文獻[6-7]研究多旋翼無人機的微多普勒模型，并對散射點積分回波模型和物理光學眼面回波模型進行比較，分析極化特性對微多普勒的影響，但沒有討論葉片數(shù)、轉(zhuǎn)速、初始相位等對微多普勒的影響。文獻[8-10]研究旋翼數(shù)、轉(zhuǎn)動頻率、葉片數(shù)及葉片長度的提取方法，但在實際場景中利用旋翼調(diào)制譜線估計葉片數(shù)較難，而采用倒頻譜估計旋翼數(shù)及轉(zhuǎn)速需要較多的有效脈沖回波數(shù)，運算速度較慢且精度不夠高。

四旋翼無人機是目前最常見的無人機類型，本文以四旋翼無人機為研究對象，對小型無人機的微多普勒特性與特征提取進行研究。首先建立多旋轉(zhuǎn)旋翼的回波模型，并通過仿真實驗分析葉片數(shù)目、旋翼轉(zhuǎn)速和初始相位等參數(shù)對微多普勒特征的影響，接著提出一種基于Gabor變換的瞬時頻率估計與FFT相結(jié)合的微動特征提取算法。該算法利用Gabor變換得到時頻特征，在此基礎(chǔ)上通過瞬時頻率極大值法提取微多普勒頻率展寬，并對瞬時頻率采用FFT提取旋翼數(shù)和轉(zhuǎn)動頻率，由此得到葉片長度估計值。大疆精靈3S無人機的實測數(shù)據(jù)驗證了算法的有效性，能較為準確地提取無人機的微動參數(shù)。

1 多旋翼無人機回波建模

美國海軍實驗室的Victor C. Chen教授建立了旋轉(zhuǎn)旋翼葉片回波的積分模型。對于多旋翼的無人機，假設(shè)各旋翼葉片的RCS相同且均為1，在直升機旋翼模型[2]基礎(chǔ)上，本文構(gòu)建出多旋轉(zhuǎn)旋翼的回波模型如下：

exp{-jΦm,k(t)},

(1)

其中

(k=0,1,2,…,N-1;m=1,2,…,M).

(2)

式中:M為旋翼總數(shù)目，N為單個旋翼總的葉片數(shù)目，L表示旋翼葉片長度，R0m為雷達到第m個旋翼中心的距離，z0m表示第m個旋翼葉片的高度，βm為雷達到第m個旋翼中心的俯仰角(近似等于雷達到無人機軸中心的俯仰角，即β1=β2=…=βM=β)，Ωm為第m個旋翼的轉(zhuǎn)動角頻率，φ0m為第m個旋翼的初始旋轉(zhuǎn)角。

回波信號的瞬時多普勒頻率可由信號的相位函數(shù)求時間導數(shù)得到，對式(2)求時間導數(shù)得到第m個旋翼的第k個葉片的等效瞬時微多普勒頻率為

(3)

而葉片上的一個散射點P的瞬時多普勒頻率是

(4)

式中:lP為散射點P到旋翼旋轉(zhuǎn)中心的距離，且0≤lP≤L。

式(4)表明，葉片上散射點的微多普勒頻率是正弦變化的曲線，曲線條數(shù)即葉片數(shù)目，且正弦曲線頻率ω與葉片轉(zhuǎn)動角頻率相同，則有

ωm,k=Ωm.

(5)

葉片葉尖處的多普勒頻率幅度值最大，則最大微多普勒頻率為

(6)

目標平動速度為0時，葉片的微多普勒頻率的最大展寬為

(7)

由此推導出旋翼的葉片長度為

(8)

通過估計的旋翼數(shù)、葉片數(shù)、葉片長度、旋翼轉(zhuǎn)速等參數(shù)可判斷無人機型號和運動狀態(tài)，達到識別多旋翼無人機的目的。

2 多旋翼無人機微多普勒特性分析

本文以四旋翼無人機為典型目標進行微多普勒特性分析和特征提取，四旋翼無人機結(jié)構(gòu)如圖1所示，其中旋翼1、3逆時針方向旋轉(zhuǎn)，旋翼2、4順時針方向旋轉(zhuǎn)。

圖1 四旋翼飛行器的結(jié)構(gòu)示意圖Fig.1 Diagram of four-rotor UAV

由式(3)可知，多旋翼無人機微多普勒特征是由M×N條正弦形式的曲線組成，且受到載頻、旋翼數(shù)目、旋翼轉(zhuǎn)速、葉片數(shù)、葉片長度、初始相位和俯仰角的影響。其中，載頻、葉片長度和俯仰角僅與微多普勒頻率幅度有關(guān)，而旋翼數(shù)目、旋翼轉(zhuǎn)速、葉片數(shù)和初始相位將影響微多普勒特征曲線的幅度和相位，下面對主要參數(shù)的影響進行仿真分析。

2.1 不同葉片數(shù)目對微多普勒特性的影響

由式(3)可知，多旋翼無人機微多普勒特征曲線條數(shù)及相位與葉片數(shù)目有關(guān)。

設(shè)雷達微波為Ka波段，載波頻率為34.6 GHz，帶寬為1.2 GHz，脈沖重復頻率PRF為125 kHz。雷達俯仰角β為10°，方位角α為0°，雷達到無人機軸中心距離是100 m，無人機運動速度為0，旋翼數(shù)目為4個，各旋翼葉片的初始旋轉(zhuǎn)角φ0為0°，旋翼轉(zhuǎn)動速度均為30 r/s，各旋翼旋轉(zhuǎn)中心至無人機軸中心距離15 cm，葉片長度為12 cm，觀測時長0.1 s。采用式(1)仿真多旋翼散射點積分的回波模型，采用Gabor變換對不同葉片數(shù)目的旋翼微多普勒信號進行仿真，結(jié)果如圖2所示。

圖2 不同葉片數(shù)的四旋翼無人機微多普勒特征Fig.2 Micro-Doppler characteristics of four-rotor UAVs with two and three blades

分析式(1)和圖2(a)、2(b)，在時域回波出現(xiàn)sinc峰值的時刻，對應(yīng)的時頻域會有一條頻率帶，是時頻域閃爍現(xiàn)象，在葉片與雷達視線垂直時刻產(chǎn)生[11-12]；由圖2知，四旋翼無人機回波的時域閃爍和對應(yīng)的時頻域閃爍次數(shù)相同，閃爍頻率為fT=kNfr(fr為轉(zhuǎn)動頻率，N為葉片數(shù)，N為奇數(shù)時，k取2；N為偶數(shù)時，k取1)[13]，2葉片0.1 s內(nèi)出現(xiàn)6次閃爍，3葉片0.1 s內(nèi)出現(xiàn)18次閃爍，可估計出葉片轉(zhuǎn)動頻率30 Hz；從圖2(b)、2(d)發(fā)現(xiàn)，四旋翼無人機葉片的微多普勒關(guān)于平均多普勒頻率對稱，原因是四旋翼無人機有兩組旋轉(zhuǎn)方向相反的旋翼，產(chǎn)生的微多普勒頻率相反，而在葉片數(shù)為3時有明顯的零頻帶；圖2(b)提取的微多普勒頻率最大展寬值為fW-max=5 259×2=10 518 Hz，再根據(jù)式(8)計算得到葉片長度為L=12.28 cm，葉長估計值的相對誤差為2.33%，驗證了理論分析的正確性。

2.2 不同旋翼轉(zhuǎn)速對微多普勒特性的影響

由式(3)可知，多旋翼無人機旋翼的轉(zhuǎn)速會影響微多普勒頻率的幅度大小及微多普勒特征曲線的頻率。

在2.1節(jié)的仿真參數(shù)下，假定各旋翼葉片數(shù)設(shè)為2，觀測時間設(shè)為0.2 s，設(shè)置不同的旋翼轉(zhuǎn)速，仿真多旋翼散射點積分的回波模型，用Gabor變換進行微多普勒仿真。旋翼1、3轉(zhuǎn)速相等(30 r/s)，旋翼2、4轉(zhuǎn)速相等(35 r/s)時微多普勒仿真結(jié)果為圖3。

4個旋翼轉(zhuǎn)速相同時，各旋翼微多普勒特征重合，且只有唯一的時頻閃爍周期，如圖2(b)所示，而圖3出現(xiàn)不同的閃爍周期，表明旋翼數(shù)不只一個；圖3中明顯看出兩種不同頻率的正弦曲線，其中頻率較小的時頻閃爍次數(shù)為m1=12，對應(yīng)的轉(zhuǎn)動頻率為fr1=30 Hz，微多普勒頻率的最大值為fmdMax1=5 259 Hz，頻率較大的時頻閃爍次數(shù)為m2=14，對應(yīng)的轉(zhuǎn)動頻率為fr2=35 Hz,微多普勒頻率的最大值為fmdMax=6 238 Hz,有m1/m2=fr1/fr2=0.857 1，fmdMax1/fmdMax2=0.843 1≈fr1/fr2，則微多普勒特征的時頻閃爍頻率之比等于旋翼轉(zhuǎn)動頻率之比，微多普勒峰值頻率與葉片轉(zhuǎn)動頻率成正比。綜合以上分析，旋翼轉(zhuǎn)動頻率會影響微多普勒頻率的幅度、時頻閃爍的周期及次數(shù)。

圖3 不同轉(zhuǎn)速的四旋翼無人機微多普勒特征Fig.3 Micro-Doppler characteristics of four-rotor UAV with different rotation rates

2.3 不同初始相位對微多普勒特性的影響

由式(3)知，多旋翼無人機各旋翼葉片的初始相位影響微多普勒特征曲線的相位。

在2.1節(jié)的仿真參數(shù)下，假定各旋翼葉片數(shù)設(shè)為2，設(shè)置不同的旋翼初始相位，仿真多旋翼散射點積分的回波模型，采用Gabor變換提取旋翼的微多普勒信息。旋翼2初始相位為60°，旋翼1、3、4初始相位相同(0°)時微多普勒仿真結(jié)果為圖4(a)；各旋翼初始相位各不相同時微多普勒仿真結(jié)果為圖4(b)。

圖4 不同初始相位的四旋翼無人機微多普勒特征Fig.4 Micro-Doppler characteristics of four-rotor UAV with different initial phases

圖2(b)中4個旋翼初始相位相同時，微多普勒特征曲線重合；對比圖4(a)、4(b)，當旋翼初始相位不同時，微多普勒特征曲線的初始相位不同，且時頻閃爍時刻出現(xiàn)差異，使得微多普勒特征曲線周期發(fā)生變化。圖2(b)第2次時頻閃爍時刻是在0.022 61 s，而圖4(a)僅旋翼2初相位改變時出現(xiàn)新的時頻閃爍，其第2次時頻閃爍時刻是在0.016 94 s，時間間隔為Δt=0.005 7 s，根據(jù)時頻閃爍頻率fT=30 Hz可得兩次閃爍的初始相位差為Δφ0=2πfTΔt=61.236°，接近仿真值60°。從上述分析可知，旋翼葉片不同的初始相位會改變微多普勒峰值出現(xiàn)的時刻及閃爍次數(shù)，不改變閃爍周期，可通過對無人機不同角度的照射得到有相位差異的微多普勒特征曲線，進而對葉片數(shù)的判斷提供參考。

通過以上對多旋翼無人機微多普勒特征的建模和仿真分析，弄清微多普勒特征的影響因素，為后文的微多普勒特征提取及參數(shù)估計提供理論基礎(chǔ)。

3 多旋翼無人機微多普勒特征提取

微多普勒信號是非平穩(wěn)信號，而傳統(tǒng)的傅里葉變換不能得到與時間有關(guān)的頻率信息，聯(lián)合時頻分析提供了一種分析信號瞬時頻率隨時間變化規(guī)律的方法[11]。時頻分析方法主要包括線性變換法和雙線性變換法。線性時頻分析法主要有如短時傅里葉變換(STFT)、Gabor變換、分數(shù)階傅里葉變換(FRFT)等，雙線性時頻分析法有Wigner-Ville分布(WVD)、偽WVD分布、平滑偽WVD分布、廣義S變換等。這些方法都有各自的優(yōu)缺點。例如STFT是最常用的時頻分析法，其實現(xiàn)簡單且無交叉項，但不能同時獲得高的時間和頻率分辨率；FRFT為參數(shù)搜索類方法，無交叉項的干擾，但進行匹配搜索運算量較大，且參數(shù)估計精度受搜索步長的限制，多用于分析chirp信號[14]；平滑偽WVD分布是在WVD基礎(chǔ)上產(chǎn)生的方法，對交叉項進行了較好的抑制，但計算量增加且時頻分辨率下降[15]；廣義S變換時頻分辨率較好且能有效抑制交叉項，但窗函數(shù)的雙參數(shù)調(diào)節(jié)增加運算的復雜度。

多分量信號的時頻分析中分辨率與交叉項存在矛盾，Gabor變換是一種特殊的STFT變換，不會產(chǎn)生交叉項，且運算簡單，其采用高斯窗獲得最小的時間和頻率分辨率的乘積，時頻特征明顯。對于多旋翼無人機目標，其微多普勒信號為多分量信號，對實時性要求高，且時頻特征要突出，因此考慮采用Gabor變換進行時頻特征分析。

3.1 瞬時頻率提取

信號的時頻分布代表信號的瞬時能量，其在時頻面上的投影就是信號的瞬時頻率[16]。多分量信號通過Gabor變換后，各分量信號的峰值在時頻面上的投影等價于時頻面的極大值，提取時頻圖的極大值即為各分量信號的瞬時頻率。

設(shè)信號的離散Gabor分布Gabor(m,n)是M×N的矩陣，為了減小噪聲的影響，將Gabor(m,n)通過設(shè)定的門限值，然后按列求得極大值，由此得到所有分量信號的瞬時頻率[17]，瞬時頻率的表達式近似為

(9)

式中：fk(i)表示第i時刻第k個分量的瞬時多普勒頻率，GaborM×N(:,i)表示M×N陣列的第i列，arg(·)表示取宗量運算，peak[·]為求極大值。若目標的徑向平動速度引起的多普勒頻移為fu，則瞬時微多普勒頻率為：

fm,k(i)=fk(i)-fu, 1≤i≤N.

(10)

采用瞬時頻率極大值法可以得到微多普勒頻率的最大展寬值fW-max。

3.2 旋翼數(shù)和轉(zhuǎn)動頻率提取

多旋翼無人機旋翼旋轉(zhuǎn)產(chǎn)生的微多普勒信號是多個正弦函數(shù)組成，對其做傅里葉變換，則可得到各正弦函數(shù)的頻率[16]，即各旋翼的轉(zhuǎn)動頻率。當旋翼轉(zhuǎn)速不同時，會得到不同的頻率成分，可判斷旋翼數(shù)目。因此，對無人機旋翼產(chǎn)生的微多普勒進行時間維的傅里葉變換得到的頻譜圖，可估計微多普勒信號的正弦頻率及不同頻率成分的個數(shù)。

對式(4)進行FFT變換，得到第m個旋翼的第k個葉片的傅里葉變換結(jié)果為

[ej(φ0m+k2π/N)δ(ω-Ωm)-e-j(φ0m+k2π/N)δ(ω+Ωm)].

(11)

從式(11)可知，第m個旋翼的微多普勒傅里葉變換后的譜線應(yīng)該出現(xiàn)在其轉(zhuǎn)動頻率處。

3.3 算法流程

本文提出的多旋翼無人機微動特征提取算法的處理流程如圖5所示。

圖5 本文提出的微動特征提取算法流程圖Fig.5 Flow chart of the proposed feature parameter extraction algorithm

算法的具體實現(xiàn)過程為：將接收的Ka波段調(diào)頻連續(xù)波雷達回波脈沖按列存放為二維矩陣，對快時間維進行傅里葉變換(即脈沖壓縮)得到距離維信息，慢時間維運用高通濾波抑制靜止雜波，取二維矩陣中目標所在的距離單元沿距離維累加形成沿時間維的一維信號，對該一維信號進行Gabor變換則可得到目標的微多普勒信號。由微多普勒譜提取瞬時微多普勒頻率來估計微多普勒頻率展寬值，對提取的瞬時頻率沿多普勒頻率維疊加并進行時間維的傅里葉變換，由此可估計旋翼轉(zhuǎn)動頻率及旋翼數(shù)目，最后由公式(8)計算出葉片長度值。

3.4 性能分析與實測數(shù)據(jù)驗證

3.4.1 性能分析

在2.1節(jié)的仿真參數(shù)下，設(shè)葉片數(shù)為2，在未做積累的回波信號中分別加入信噪比為10、5、0、-5、-10、-11 dB的高斯白噪聲，表1是不同信噪比條件下微多普勒參數(shù)的提取結(jié)果，由20次仿真實驗取得的估計均值。根據(jù)式(8)可得各旋翼葉片長度的估計值為

(12)

將各旋翼葉片長度估計值的均值作為無人機葉片長度的最終估計值：

(13)

葉長L估計值的相對誤差為

(14)

式中：Le表示葉片長度的真實值。

從表1可以看出，在信號信噪比低至-11 dB時，該算法仍能準確地提取出旋翼數(shù)目和各旋翼轉(zhuǎn)速，具有一定的抗噪性；信噪比為-10 dB以上時，隨著信噪比增大，葉片長度值的估計誤差越來越小，且在-10 dB時估計誤差也能達到14.08%，而在-11 dB及以下信噪比條件下，由于噪聲較強，時頻圖中的微多普勒特征很模糊，難以估計最大微多普勒頻率的展寬值，故不能推導出葉片長度估計值。分析結(jié)果表明，本文提出的算法在低信噪比情況下，能實現(xiàn)無人機微多普勒參數(shù)的高精度提取。

表1 不同信噪比條件下微多普勒參數(shù)提取結(jié)果Table 1 Results of micro-Doppler parameters extraction under different SNR conditions

3.4.2 四旋翼無人機懸停試驗

針對旋翼數(shù)目、轉(zhuǎn)速、葉片長度的估計，進行檢測懸停的大疆精靈3S無人機微多普勒特征的試驗。試驗使用Ka波段調(diào)頻連續(xù)波雷達，載頻為34.6 GHz，帶寬1.2 GHz，PRF為125 kHz，雷達雷達定向照射，俯仰角約20°，雷達到無人機的距離約20 m左右。圖6為四旋翼無人機懸停試驗數(shù)據(jù)處理結(jié)果。

圖6 四旋翼無人機懸停試驗結(jié)果Fig.6 Test results for a hovering four-rotor UAV

圖6(a)是取1 s時頻圖做傅里葉變換的結(jié)果，使得頻率分辨率達到Hz級，圖中僅得到1種頻率成分，為31 Hz，由式(11)可知其對應(yīng)的微多普勒信號頻率應(yīng)為31 Hz，而對于2個葉片的四旋翼無人機，其微多普勒特征關(guān)于平均多普勒頻率對稱，則微多普勒信號的頻率為轉(zhuǎn)動頻率的2倍，由此可知對應(yīng)旋翼的轉(zhuǎn)動頻率應(yīng)為15.5 Hz；圖6(b)的時頻圖中能明顯看出兩種不同相位的時頻閃爍且間隔不均勻，說明至少有2個轉(zhuǎn)速均為31 r/s的旋翼，而四旋翼無人機在懸停時4個旋翼轉(zhuǎn)速相等[18]，因此另2個旋翼轉(zhuǎn)速均為15.5 r/s；從圖6(c)估計出微多普勒頻率的最大展寬值為fW-max=2 568-(-2 324)=4 892 Hz，由公式(8)計算得到葉片長度估計值為11.92 cm，其真實值為12 cm，則葉長估計值的相對誤差為0.67%，試驗證明該算法能有效提取四旋翼無人機的微動參數(shù)。

3.4.3 強雜波背景下四旋翼無人機慢速運動試驗

在強雜波背景中對空中慢速運動的大疆精靈3S無人機進行檢測，無人機垂直于雷達視線運動，背景中有樹木、小型建筑物、行人等，有3～4級風。雷達使用Ka波段調(diào)頻連續(xù)波，載頻為34.6 GHz，信號帶寬為1.2 GHz，PRF為125 kHz，雷達定向照射，俯仰角約22°，雷達到無人機的距離約為25 m。根據(jù)采集的數(shù)據(jù)提取樹木雜波背景中四旋翼無人機的微多普勒特征，回波信號脈壓后的信雜比為-25 dB。在使用本文方法之前采用高通濾波對雜波進行抑制，濾波后的信雜比為-6 dB。圖7所示為強雜波背景下四旋翼無人機慢速運動試驗處理結(jié)果。

圖7(a)得到2種頻率成分，分別為46、61 Hz，對應(yīng)的轉(zhuǎn)動頻率分別為23、30.5 Hz，則判斷有2個轉(zhuǎn)速分別為23、30.5 r/s的旋翼，由于四旋翼無人機在平移運動過程中有2對轉(zhuǎn)速分別相同的旋翼[18]，則另2個旋翼轉(zhuǎn)速也分別為23、30.5 r/s；圖7(b)是由0.2 s信號的時頻圖提取的瞬時頻率，各旋翼微多普勒頻率交叉在一起，不宜直接估計微多普勒展寬頻率，則從目標回波所在的多個距離單元中選取不同距離單元分別進行分析，結(jié)果如圖7(c)、7(d)所示；由圖7(c)、7(d)提取的最大微多普勒頻率寬度分別是fW-max1=3 792-(-3 547)=7 339 Hz，fW-max2=4 281-(-5 015)=9 296 Hz，對應(yīng)的旋轉(zhuǎn)頻率分別為23和30.5 Hz，根據(jù)式(8)，估計出葉片長度分別是12.15和11.60 cm，相對誤差分別為 1.25%和3.33%，試驗表明該算法對多旋翼無人機微動參數(shù)提取有較高的準確率, 且在雜波環(huán)境中也具有可行性。

圖7 強雜波背景下四旋翼無人機慢速運動試驗結(jié)果Fig.7 Test results for a hovering four-rotor UAV in strong clutter background

4 結(jié)論

本文在直升機旋翼回波模型的基礎(chǔ)上，建立多旋翼無人機的回波模型，并基于該模型進行仿真實驗，分析不同葉片數(shù)目、不同旋翼轉(zhuǎn)速及不同初始相位對多旋翼無人機微多普勒特性的影響。針對多旋翼無人機目標，提出一種基于Gabor變換的瞬時頻率估計與FFT相結(jié)合的微動特征提取算法，并通過四旋翼無人機懸停試驗和強雜波背景下四旋翼無人機慢速運動試驗驗證了該方法的有效性，且對旋翼數(shù)、轉(zhuǎn)速、葉片長度的估計都有較高的準確性。但對于無人機運動時轉(zhuǎn)速不穩(wěn)定的情況，該方法不再適用，在后續(xù)工作中將對這些情況進行研究。