基于ISAR像序列的多旋翼無人機參數(shù)估算

任楓軒，　王忠勇

(1.河南職業(yè)技術(shù)學院電氣工程系，鄭州　450046；　2.鄭州大學信息工程學院，鄭州　450001)

0　引言

由于多旋翼無人機具有良好的空中懸停特性及易控特性，在航拍、測繪、運輸、安保等領(lǐng)域扮演著越來越

重要的角色[1-3],但由此也引發(fā)了許多安全問題。由于質(zhì)量或操作的原因極有可能在人員密集的區(qū)域墜機而引發(fā)安全事故，多旋翼攜載的可見光、紅外、電磁等探測設(shè)備也使得隱私越來越容易被窺視，曾發(fā)生過無人機黑飛影響航班安全的事件，因此我國于2013年相繼推出了民用無人機的空中管理條例，對無人機的航速、重量和飛行高度等都做出了相關(guān)要求。由于雷達電磁波受天氣影響較小，在目標探測中發(fā)揮重要作用，雷達ISAR像為目標散射源在成像平面上的投影表征，可以最直觀地反映目標的形態(tài)，在連續(xù)的ISAR像中更可以提取出目標的運動參數(shù)，為了對空中飛行的無人機進行更準確的識別，本文利用ISAR像序列對多旋翼目標進行探測[4]。目前公開發(fā)表的文獻中幾乎沒有對多旋翼無人機特征參數(shù)的研究，但在利用ISAR像獲取目標特征參數(shù)方面的研究較多。文獻[5]通過仿射配準對ISAR像進行橫向定標，再通過聚類算法提取強散射源，進而估算了目標長度特征；文獻[6]分析了目標在成像平面上不隨雷達參數(shù)變化的穩(wěn)定特征量，建立ISAR像與結(jié)構(gòu)參數(shù)的關(guān)系，提出了利用穩(wěn)定特征量反估結(jié)構(gòu)參數(shù)的方法；文獻[7]分析了空間目標運動模型及成像機理，推導了強散射源在ISAR像平面上的理論分布，并以此來估算特征參數(shù)。為此，借助雷達二維像高分辨特征，采用ISAR像序列對多旋翼無人機的旋翼片數(shù)及轉(zhuǎn)動周期進行識別，通過仿真實驗驗證了算法的有效性，且具有較低的估算誤差。

1　多旋翼無人機運動模型

多旋翼無人機一般由4片以上的旋翼構(gòu)成，每個旋翼由2片槳葉構(gòu)成，旋翼軸對稱分布且共面，為了對多旋翼無人機運動模型進行建模，首先以旋翼中心為原點，在旋翼旋轉(zhuǎn)平面內(nèi)任意方向為xrotor1軸，旋翼轉(zhuǎn)動角速度矢量方向為zrotor1軸，按右手螺旋法則確定yrotor1軸，由此建立旋翼坐標系[8]，如圖1所示。

圖1　旋翼坐標系Fig.1　Rotor coordinate system

圖1中，φrotor1_0為旋翼在旋翼坐標系中的旋轉(zhuǎn)角度，設(shè)旋翼旋轉(zhuǎn)角速度為ωrotor1，且其初始轉(zhuǎn)角為φrotor1_init_0，則φrotor1_0可表示為：φrotor1_0=φrotor1t+φrotor1_init_0。對旋翼上任意的散射源位置r0，在t時刻該散射源在旋翼坐標系下的位置r1可表示為

(1)

由于多旋翼無人機的各個旋翼共面且軸對稱分布，所以在建立機體坐標系時以多旋翼無人機的對稱軸為zbody軸，以各旋翼所在平面內(nèi)任意方向為xbody軸，再由右手螺旋定則確定zbody軸，不失一般性，設(shè)槳葉片數(shù)為4，由此建立機體坐標系如圖2所示。

圖2　機體坐標系Fig.2　Coordinate system of the aircraft body

圖2中，φbody為多旋翼無人機在機體坐標系中的轉(zhuǎn)動角度。設(shè)槳葉中心到zbody軸的距離為l，槳葉坐標系中的r1在機體坐標系下的坐標r2表示為

(2)

由于多旋翼機體存在俯仰、斜側(cè)和旋轉(zhuǎn)等運動，以機體中心為原點，以垂直于大地表面垂直向上為z軸，以過雷達視向與z軸構(gòu)成的平面且垂直于z軸的直線為x軸，利用右手螺旋定則確定y軸[9]，los為雷達視向，由此構(gòu)建本體坐標系如圖3所示。

圖3　本體坐標系Fig.3　Ontology coordinate system

不考慮平動運動，多旋翼的任意運動可用歐拉矩陣Rconi來表示[10]，令φ為錐旋角，θ為章動角，ψ為自旋角，則歐拉矩陣Rconi(t)可表示為

Rconi(t)=

(3)

那么旋翼上r0處的散射源在本體坐標系下的坐標r可表示為

r=Rconi·r2。

(4)

2　多旋翼無人機ISAR像成像分析

由于多旋翼無人機在空中的運動軌跡復雜，各個旋翼轉(zhuǎn)動速率較快，獲取其整個旋轉(zhuǎn)狀態(tài)下回波的時間較短，即可認為在旋翼完成一周自旋的情況下多旋翼無人機機身相對于雷達視向的位置未變化，其機身自身的微動及槳葉快速旋轉(zhuǎn)都使散射源微多普勒變動且槳葉成像轉(zhuǎn)角積累較大，從而導致利用傳統(tǒng)的距離-多普勒成像算法成像時會出現(xiàn)散焦現(xiàn)象，而時頻分析技術(shù)的出現(xiàn)可以有效地解決由于多普勒時變引起的散焦現(xiàn)象，本文采用距離-瞬時多普勒成像算法完成多旋翼無人機成像[11]。

2.1　成像機理

設(shè)雷達發(fā)射線性調(diào)頻信號，對回波做解線頻調(diào)處理，通過FFT變換得到慢時間頻率域的回波信息，忽略剩余視頻項和包絡(luò)斜置項，有[7]

Sif(tm,f)=

(5)

再對每個距離單元中的回波進行時頻變換，采用Gabor變換對式(5)進行變換[12]，有

power(t,ω,l)=

(6)

上式選取特定的成像時間t，可獲得t時刻散射源在雷達視向及多普勒軸上的散射能量分布，即ISAR像，如果選擇連續(xù)的成像時間，則可獲得連續(xù)的ISAR像分布。

2.2　成像仿真實驗

設(shè)槳葉片數(shù)為4，槳葉半徑為13 cm，槳葉中心到機體中心的距離l為0.6 m，第1片旋翼在機體坐標系中的轉(zhuǎn)動角度φbody為30°，相應的第2片到第4片旋翼

在機體坐標系中的轉(zhuǎn)動角度為120°，210°，300°，設(shè)第1片旋翼旋轉(zhuǎn)角速度ωrotor1為10π，且其初始轉(zhuǎn)角φrotor1_init_0為0°；第2片旋翼旋轉(zhuǎn)角速度ωrotor2為10π，且其初始轉(zhuǎn)角φrotor2_init_0為120°；第3片旋翼旋轉(zhuǎn)角速度ωrotor3為10π，且其初始轉(zhuǎn)角φrotor3_init_0為250°；第4片旋翼旋轉(zhuǎn)角速度ωrotor4為10π，且其初始轉(zhuǎn)角φrotor4_init_0為330°，雷達視向與z軸之間的夾角α為45°，機身歐拉角中的錐旋角φ為40°，章動角θ為20°，自旋角ψ為30°；設(shè)雷達帶寬為4 GHz，起始頻率為8 GHz，脈沖重復頻率100 kHz，脈寬640 μs，信噪比為30 dB，快時間域采樣點數(shù)為64，慢時間域采樣點數(shù)為3200，回波仿真方法采用物理光學法，對回波數(shù)據(jù)的距離域及多普勒域進行10倍插值后，提取慢時間域長度為0.2 s內(nèi)的目標回波數(shù)據(jù)，將其分成0.2/25 s的時間間隔依次提取ISAR像，對多旋翼無人機ISAR像進行仿真，結(jié)果如圖4所示。共得到25幅圖像，第1行由左向右依次是第1～5幅ISAR像，第2行由左向右依次是第6～10幅ISAR像，依次類推，且豎向為多普勒軸，相應的窗長為[-781.25 Hz，781.25 Hz]，橫向為距離向軸，相應的窗長為[-1.2 m，1.2 m]。

圖4　慢時間域成像間隔為0.2/25 s的成像結(jié)果Fig.4　Results of slow time domain imaging at a interval of 0.2/25 s

從仿真結(jié)果可以發(fā)現(xiàn)，槳葉中心位置始終位于0 Hz處，而從雷達視向看去，多旋翼各個槳葉的旋轉(zhuǎn)中心并不在一條直線上，這是由于ISAR像橫向是由多普勒來決定的，由于在短時間內(nèi)假設(shè)多旋翼機體運動速度為0，所以每個槳葉中心位置處的運動速率為0，故在成像平面上槳葉的中心位置處始終位于0 Hz處；同時，不難發(fā)現(xiàn)當雷達視向與槳葉垂直時，槳葉的散射最強(對應以上各幅ISAR像中的粗豎線)，且槳葉散射最強時，各槳葉在多普勒軸上的分布長度是不同的，這是由于各個槳葉的轉(zhuǎn)動速度不同引起的多普勒分布不同，如果對各個槳葉分別進行橫向定標，那么獲取的各

個槳葉的尺寸是相同的。

3　旋翼數(shù)量及各槳葉轉(zhuǎn)速估算

當雷達視向與槳葉垂直時會出現(xiàn)強散射，由此可以提取旋翼數(shù)量及槳葉轉(zhuǎn)速。對ISAR像數(shù)據(jù)進行預處理，得到強散射在時間-距離向上的分布，再從中提取特征參數(shù)。

3.1　數(shù)據(jù)預處理

利用人眼觀測連續(xù)的ISAR像序列可以直接獲取槳葉片數(shù)，本文要實現(xiàn)槳葉片數(shù)的參數(shù)化估算，當槳葉垂直于雷達視向時，槳葉散射最強，設(shè)此時對應的時間序列為i，將對應的ISAR像強度分布矩陣Ii中的元素按距離單元依次求和，得到第i個時間序列散射強度沿距離像序列的分布向量li，即

li=sum(Ii)

(7)

不難發(fā)現(xiàn)，向量li應該在含有槳葉的那幾個距離單元中的幅值最高，如果將所有時刻的li匯總成一個矩陣LT，則

LT=(l1l2…lN)N×M

(8)

即LT為N×M維的矩陣，其中，N為時間序列的總數(shù)，M為距離向單元數(shù)。將采樣時間改為0.4 s，ISAR像成像時間間隔改為0.4/200 s，獲取相應的LT矩陣并繪于圖5中。

圖5　LT矩陣分布Fig.5　LT matrix distribution

3.2　特征提取

從圖5可看出，圖中黑色圓圈處均存在極大值點(rs,ts)，假設(shè)極大值點數(shù)為S，s取值范圍為1,…,S。如果已經(jīng)提取出黑色圓圈處的極大值點在距離向上的投影rs，那就可以確定槳葉片數(shù)，并可以提取相應的LT(n,rs/Δr)，其中，n為時間序列標號，取值范圍為1,…，N，Δr為距離分辨率，rs/Δr為相應的距離分辨單元標號。對LT(n,rs/Δr)進行自相關(guān)[13-14]處理獲得的周期的2倍，即該槳葉的轉(zhuǎn)動周期

T=2*corr(LT(n,rs/Δr))

(9)

式中，corr表示做自相關(guān)處理獲取周期。

圖6　向量L的歸一化分布Fig.6　Normalized distribution of vector L

從圖6中可以看出，L是一個多峰值分布的向量，其峰值個數(shù)及峰值位置分別是槳葉片數(shù)及槳葉在雷達視向上的位置。為了實現(xiàn)參數(shù)化提取，引入基本的蟻群算法[15]提取極大值點位置。但向量L中除了槳葉中心的位置處出現(xiàn)峰值外，在其他峰谷區(qū)域內(nèi)存在許多局部極大值點，這是由于噪聲引起的，而基本的蟻群算法無法對此進行區(qū)分，可通過設(shè)定閾值的方式來進行篩選，不失一般性，令η·max((rs,ts),s=1,…,S)為閾值，其中，η為篩選系數(shù)，取值范圍為0～1，由此可以確定槳葉片數(shù)及槳葉中心在雷達視向上的投影位置，再依據(jù)前面的方法可以獲得槳葉轉(zhuǎn)動周期。

4　仿真結(jié)果分析

對旋翼數(shù)目及旋翼轉(zhuǎn)動周期進行估算，需設(shè)置的仿真參數(shù)有修正系數(shù)η取0.8，基本蟻群算法的仿真參數(shù)設(shè)置如下：螞蟻數(shù)量為1000、螞蟻移動次數(shù)為50、信息素揮發(fā)系數(shù)為0.8、轉(zhuǎn)移概率常數(shù)為0.2、篩選系數(shù)為0.8，蟻群算法搜索獲取的極大值[16]結(jié)果如圖7所示。

圖7　蟻群算法搜索獲取的極大值結(jié)果Fig.7　Maximum results of ant colony algorithm search

圖7滿足條件的極大值共有4個，相應的雷達視向投影分別為-0.057 m，-0.007 m，0.007 m，0.057 m，提取矩陣LT中相應距離單元上的向量，分別列于圖8中。對圖8中的4個向量分別用式(9)中的方法進行周期估算，估得周期分別為0.212 s，0.216 s，0.200 s，0.204 s。相應的誤差分別為6%，8%，0%，2%。再考慮雷達視角α分別取30°，45°及60°，且信噪比分別取30 dB，20 dB，10 dB，0 dB的情況，將槳葉片數(shù)及槳葉轉(zhuǎn)動周期估算誤差列于表1。

圖8　提取的矩陣LT中相應距離單元上的向量分布Fig.8　Vector distribution on the corresponding distance unit in the extracted matrix LT

雷達視角轉(zhuǎn)動周期信噪比/dB-30-20-10010203030°旋翼1周期/%～～24222旋翼2周期/%～～42200旋翼3周期/%～～44468旋翼4周期/%～～26668旋翼片數(shù)564444445°旋翼1周期/%～～44226旋翼2周期/%～～22228旋翼3周期/%～～44680旋翼4周期/%～～44662旋翼片數(shù)854444460°旋翼1周期/%～～44222旋翼2周期/%～～22220旋翼3周期/%～～24868旋翼4周期/%～～44666旋翼片數(shù)7244444

由仿真結(jié)果可以看出，在信噪比大于-10 dB的情況下，旋翼數(shù)目及各旋翼轉(zhuǎn)動周期均可有效估算，并且旋翼轉(zhuǎn)動周期估算誤差小于8%；而在信噪比低于-20 dB時，由于旋翼數(shù)估算錯誤，致使提取出的周期與真實值差別較大，所以并未列出此時的周期估算誤差。另外，發(fā)現(xiàn)表中所有的周期估算誤差的百分數(shù)都是2的整數(shù)倍，這是由于ISAR像采樣間隔導致的，本文ISAR像采樣間隔為0.4/200=0.002 s，那么對提取的時間序列進行自相關(guān)處理后的時間分辨率也為0.002 s，由式(9)可知估算的周期分辨率為0.004 s，除以真實周期0.2 s，可得周期估算的百分數(shù)分辨率為2。

5　結(jié)束語

借助雷達二維像高分辨特征，利用ISAR像序列提取旋翼強散射特征進行參數(shù)估算，通過對多旋翼無人機運動模型建模，對多旋翼無人機ISAR像成像機理分析，再利用強散射在時間及雷達視向上的分布，結(jié)合蟻群算法及自相關(guān)對多旋翼無人機的旋翼片數(shù)及其轉(zhuǎn)動周期進行估計。通過仿真實驗結(jié)果表明：當信噪比大于-10 dB時，各個參數(shù)估算誤差不會超過8%，這說明提出的算法可有效提取相應參數(shù)，對飛行的旋翼無人機的探測識別提供了有力的數(shù)據(jù)支持，為低空安全管理提供保障。

[1]戴冬,王果,王磊.博弈論在固定翼無人機地面目標跟蹤控制中的應用[J].計算機工程,2016,42(7):287-292,298.

[2]高扉扉,陳念年,范勇,等.一種旋翼式無人機的視覺著陸位姿估計方法[J].電光與控制,2017,24(2):35-38,80.

[3]王慶賀,萬剛,曹雪峰,等.基于數(shù)據(jù)融合的多旋翼無人機定位控制器設(shè)計[J].系統(tǒng)仿真學報,2016,28 (10):2593-2599.

[4]張凌曉,王寶順,賀思三,等.基于圖像旋轉(zhuǎn)匹配的組網(wǎng)雷達ISAR圖像橫向定標[J].計算機工程與科學,2015,37(4):796-801.

[5]JIN G H,GAO X Z,DONG Z.Two-dimensional length extraction of ballistic target from ISAR images using a new scaling method by affine registration[J].Defence Science Journal,2014,64(5):458- 463.

[6]徐少坤,劉記紅,袁翔宇,等.基于ISAR圖像的中段目標二維幾何特征反演方法[J].電子與信息學報,2015, 37(2):339-345.

[7]束長勇,陳世春,吳洪騫,等.基于ISAR像序列的錐體目標進動及結(jié)構(gòu)參數(shù)估計[J].電子與信息學報,2015,37(5):1078-1084.

[8]陳麗城,李春濤,張孝偉,等.無人機地面動力學建模及分析[J].計算機仿真,2016,33(6):13-18,35.

[9]劉東輝,奚樂樂,孫曉云.矢量拉力垂直起降無人機姿態(tài)縱向控制研究[J].計算機工程與應用,2017,53(1):260-264.

[10]XU S K,LIU J H,WEI X Z,et al.Wideband electromagnetic characteristics modeling and analysis of missile targets in ballistic midcourse[J].Science China Technolo-gical Sciences,2012,55(6):1655-1666.

[11]芮力,錢廣紅,張國慶,等.基于自適應最優(yōu)核時頻分布理論的ISAR成像方法[J].電光與控制,2014,21(7):46-50,102.

[12]?ZDEMIR C.Inverse synthetic aperture radar imaging with MATLAB algorithms[M].New York:John Wiley & Sons,2012:274-287.

[13]周磊磊,羅炬鋒,付耀先,等.低信噪比下基于自相關(guān)函數(shù)的頻率估計方法[J].華中科技大學學報:自然科學版,2014,42(4):45- 49.

[14]陳剛,王鵬飛,李金玲.基于自相關(guān)函數(shù)的模糊時間序列優(yōu)化算法[J].控制與決策,2015(10):1797-1802.

[15]DORIGO M,DI CARO G.Ant colony optimization:a new meta-heuristic[J].Proceedings of the Congress on Evolutionary Computation,1999.doi:10.1109/CEC.1999.782657.

[16]武光輝,童創(chuàng)明,李西敏,等.逆合成孔徑雷達目標成像識別優(yōu)化仿真[J].計算機仿真,2016,33(6):9-12,339.