李 飛,黃斌科,馮恩信,汪文秉
(西安交通大學電信學院微波所,陜西西安 710049)
真實環境中的運動目標不僅存在整體的平動,目標或者目標部件還往往伴隨有加速、振動、旋轉等非剛體性質的微運動。目標微運動對于雷達散射信號會產生調制作用,在目標多普勒頻率附近產生邊帶,這種現象稱為微多普勒現象[1]。對于目標微運動信號的獲取和處理正在成為一個新的技術熱點[2-4]。人在行走過程中由于手臂和腿等部分的擺動,存在明顯的微運動特點,可以利用這些特征進行人體的探測與識別。因此研究人體運動的微多普勒特征,在城市反恐、自然災害救援等方面有著實際意義,利用雷達對人運動特征進行分析近來得到了廣泛關注[5-7]。
在微多普勒效應分析中需要得到人運動情況下的雷達回波。2002年Geisheimer應用運動捕捉(Motion capture)技術對于人行走的雷達回波進行了模擬,并分別討論了每個部位的回波特征[5]。V.Chen基于歐拉角方法,對人運動回波的合成和微多普勒特征提取進行了研究[6]。由于條件限制或者出于研究一種不受具體設備影響的通用模型的考慮,在研究人體運動的微多普勒特征以及應用時首先需要解決的問題是人運動情況下的雷達回波模擬。在已有的工作中,對回波的仿真主要是基于運動捕捉技術和歐拉角方法,在仿真過程中需要復雜的矩陣運算,這種三維坐標不便于用來計算旋轉關節的位置,還可能出現萬向節死鎖問題[8],因此本文針對一個常用的運動模型特點,提出利用四元數方法簡化計算過程。
本文應用的人體模型如圖1所示,該模型由12個不同的部件組成,可以比較準確地描述人的運動情況。其中頭部是一個球體,而其它部件用圓柱或者橢球體表示。

圖1 人體組合模型Fig.1 Human model
在人行走過程中的任意時刻,各個部分的運動都可以被分解為相對于父關節的轉動和相對于人體參考點的平動。仿真中各部位運動幅度(轉動角度和平移大小)來源于文獻[9],下文工作中應用這些數據計算了人運動時各個部件的實時位置。
四元數(Quaternion)理論是數學家Hamilton于1843年首先提出的。其形式為q=w+x i+y j+z k,因此四元數概念是復數的一個推廣。其中w為實部或者標量,i、j、k為虛數部分。假設有兩個四元數:

那么它們共軛和相乘的定義如下:

在旋轉位置計算中,一般應用歐拉角的方法,即令目標依次繞三個坐標軸旋轉后得到旋轉結果,這一方法比較直觀,但存在被稱為萬向節死鎖的問題。在處理人運動時,這一問題可能會導致出現角度插值后運動不自然甚至反向運動的情況。另外歐拉角方法需要復雜的矩陣運算。因此本文利用了另外一種旋轉算法——四元數方法。和歐拉角方法相比,四元數方法解決旋轉問題非常有效[8],尤其是在分析人運動時由于獲得的數據量不足,需要進行插值的情況下。
四元數具有一個基本性質:假設經過原點的向量A繞向量B轉動θ角度,旋轉結果可以表示為[8]:

該式中四元數R=[cos(θ/2),B*sin(θ/2)],R*是R的共軛,由式(1)給出,V=[0,A]。計算時R需要預處理為單位四元數,即使其模值為1。
如圖2所示,有時需要將式(3)推廣到旋轉軸為過任意定點p的向量pa的情況。

圖2 向量繞過任意點的軸旋轉Fig.2 Vector rotation around axis through any point
當向量 pv繞向量pa旋轉θ角度到向量 pv1時,旋轉后向量ov1所對應的四元數公式如式(4)[10]:

該式中四元數P=[0,op],其他四元數變量與式(3)中意義是相同的,按相同方法構造。由式(3)、式(4)就可以對人運動的實時位置進行計算。
在對人行走時各組合部件位置的計算中,首先需要考慮坐標系的轉換。全局坐標系是以人的腰部為原點建立的,并以該點為計算基準點。而各部件旋轉位置的計算是在以該部分父關節點為原點的局部坐標系下進行,例如小腿部分的旋轉計算需要在膝關節為原點的坐標系下進行。因此計算中首先需要獲得父關節在運動時刻的位置,即位置的計算是由基準點開始向遠端部位順序進行的。然后由式(3)、(4)計算運動過程中各部件的位置。在計算中各部件的旋轉軸所對應的式(3)中的四元數V是由經過父關節的垂直軸來構造的,例如在圖1的情況下,小腿部分的旋轉軸是x軸,旋轉對應的四元數R則由旋轉部分的尺寸來確定。
圖3是應用四元數旋轉式(1)~式(4)計算各部件位置數據后得到的人行走時腳踝關節速度曲線,該結果與文獻[6]經歐拉角方法得到的結果一致。

圖3 人行走時腳踝關節的運動速度曲線Fig.3 Velocity of ankle when human walks
前面工作給出了一個常用人體運動模型,并且利用四元數方法計算了不同時刻行人各部件的位置。本節將基于這一結果進行人行走回波信號的仿真。
在雷達方程里用雷達散射截面(RCS)來描述人體各個部分的散射,RCS大小是與目標形狀、入射角度等相關的。圖4給出了三種不同形狀的物體的RCS曲線,可以看到和橢球體相比,圓柱體在電磁波入射角度變化時其RCS變化要劇烈的多。因此當接近或者剛好為后向入射時,圓柱體結構的RCS會產生激烈的變化,從而導致仿真雷達回波在該入射角度以及附近產生峰值,反映在仿真結果的時頻圖中會出現類似于沖激函數造成的頻域擴散,這與實際情況不符合;因此在本文計算中人體各部件的建模使用橢球體,文獻[11]指出應用該結構仿真結果更接近實際測量值。

圖4 球體、圓柱體、橢球體RCS隨入射角度變化曲線Fig.4 RCS of sphere,ellipsoid and cylinder
當橢球體橫截面為圓時,其RCS可以由式(5)來描述[6]:

電磁波經過人體部件時會發生散射。本文對于行人目標的回波計算采用走停模型來處理,這一模型假設在同一個脈沖周期內目標靜止,在相鄰脈沖之間目標運動,可以比較精確地描述低速運動目標的散射回波。
在計算人體運動的散射回波時,可將人體等效為如圖1所示的12個部件組成,這一模型可以比較精確地描述人體的運動。其散射信號由位于各部件幾何中心的散射點單獨散射信號的矢量合成來表示[6]。用該方法合成的結果比較接近真實觀測的信號,還可以通過插入更多的散射點數目來提高建模精度,但是在回波仿真中計算量會增加[5]。本文中,當雷達工作頻率為 f時,得到的回波可表示為:

式中,td,m=2r d,m/c代表第m個散射點的回波時延,rd,m對應于該部分散射中心到雷達的距離,c為光速,總的散射點數目M=12,各散射點的相對散射強度Am由雷達方程式(7)給出:

這里幅度值Am是接收的信號功率(W),P avg是平均發射功率(W),G是天線增益,λ是雷達工作波長(m),L是損耗,σm是該部分的RCS(m2),按照式(5)進行計算,橢球體尺寸參數由人體該部分的典型規格[11]給出,RCS公式中的夾角θ則是由入射信號方向矢量與該時刻橢球體長軸方向矢量夾角來表示。根據式(6)、式(7),可以得到人行走時的雷達回波數據。仿真中雷達位于坐標原點,人從20 m遠處沿著雷達視線方向朝向原點運動,移動速度1.3 m/s,雷達工作頻率10 GHz,脈沖重復頻率(PRF)為1 000。仿真結果如圖5所示。

圖5 人行走雷達回波仿真結果Fig.5 Radar echo of human walking
由于合成的行人雷達回波是非平穩信號,從中難以直接得到有價值的信息,因此根據非平穩信號的特點對其進行時頻分析處理。本文采用自適應核函數的方法[12],時間窗為長度為31的 Hamming窗,時頻譜圖如圖6(a)所示。從圖6(a)可以看出人行走中身體不同部分所對應的信號,其中散射信號最強的位置對應于軀干及相鄰的大腿等部件的運動,而多普勒頻移最大的位置對應于腳的運動,這和實際情況是一致的:實際運動中人體各部件的速度是不一致的,軀干速度相對穩定,而越到身體邊緣處的運動速度變化越劇烈,人行走時腳的運動速度可能遠大于人體軀干的速度。為了驗證本文合成信號的可靠性,圖6(b)給出了利用運動捕捉技術合成雷達回波并進行時頻分析處理后的結果[5],可見兩者之間具有較好的一致性。


圖6 人行走雷達回波的時頻譜Fig.6 Spectrum of human walking
本文提出利用四元數來計算人體運動中伴隨有微運動的各組成部件位置,并用計算結果進行人運動雷達回波合成的方法。這一方法可以有效地減少位置計算中的矩陣運算,避免出現萬向節死鎖等情況。應用該方法仿真獲得的人體運動散射回波與國外實測數據對比結果表明了這一方法的有效性。
對人運動的多普勒譜圖進行分析,可從中識別出人體組成部件的微多普勒效應。下一步工作將對合成運動回波的微多普勒現象進行具體分析,提取合理的微多普勒特征應用于人目標身份等方面。
[1]Chen V C,Li F,Ho S,et al.Micro-doppler effect in radar phenomenon,model and simulation study[J].IEEE Trans.On AES,2006,42(1):2-21.
[2]莊釗文,劉永祥,黎湘.目標微動特性研究進展[J].電子學報,2007,35(3):520-525.ZHUANG Zhaowen,LIU Yongxiang,LI Xiang.The achievements of target characteristic with micro-motionActa[J].Electronica Sinica,2007,35(3):520-525.
[3]高紅衛,謝良貴,文樹梁,等.基于微多普勒分析的彈道導彈目標進動特性研究[J].系統工程與電子技術,2008,30(1):50-52.GAO Hongwei,XIE Lianggui,WEN Shuliang,et al.Research on precession of ballistic missile warhead based on micro-Doppler analysis[J].Systems Engineering and E-lectronics,2008,30(1):50-52.
[4]李寶柱,袁起,何佩琨,等.目標自旋引起微多普勒的補償新方法[J].現代雷達,2008,30(10):49-51.LI Baozhu,YUAN Qi,HE Peikun,et al.New method for compensation of micro-doppler caused by target spin[J].Modern Radar,2008,30(10):49-51.
[5]Jonathan Geisheimer.A high-resolution doppler model of human gait[J].Proceedings of SPIE,2002:8-18.
[6]Chen V C.Doppler signatures of radar back-scattering from objects with micro-motions[J].IET Signal Process,2008,2(3):291-300.
[7]Ram SS,Hao Ling.Microdoppler signature simulation of computer animated human and animal motions[C]//Antennas and Propagation Society International Symposium,SAN Diego:IEEE,2008:679-699.
[8]鄧恩,帕貝利.3d數學基礎:圖形與游戲開發[M].史銀雪,譯.北京:清華大學出版社,2005:135-171.
[9]Ronan Boulic,N Thalmann,D Thalmann.A global human walking model with real-time kinematic personification[J].Visual Computer,1990(6):344-358.
[10]李亞平,黃崇超.實四元數組與旋轉[J].武漢水利水電大學學報,1995,28(6):607-612.LI Yaping,HUANG Chongchao.Using real quaternions to represent rotation in three dimensions[J].Journal of Wuhan University of Hydraulic and Electric Engineering,1995,28(6):607-612.
[11]Dorp van P,Groen F C A.Human walking estimation with radar[J].IEEProceedings on Radar,Sonar and Navigation,2003,150(5):356-365.
[12]Deming Zhang.Adaptive time frequency analysis[CP/OL].[2009-04-11]http://www.math works.com/matlabcentral/fileexchange/11551.