基于四元数的人运动雷达回波计算方法

李 飞,黄斌科,冯恩信,汪文秉

(西安交通大学电信学院微波所,陕西西安 710049)

0 引言

真实环境中的运动目标不仅存在整体的平动,目标或者目标部件还往往伴随有加速、振动、旋转等非刚体性质的微运动。目标微运动对于雷达散射信号会产生调制作用,在目标多普勒频率附近产生边带,这种现象称为微多普勒现象[1]。对于目标微运动信号的获取和处理正在成为一个新的技术热点[2-4]。人在行走过程中由于手臂和腿等部分的摆动,存在明显的微运动特点,可以利用这些特征进行人体的探测与识别。因此研究人体运动的微多普勒特征,在城市反恐、自然灾害救援等方面有着实际意义,利用雷达对人运动特征进行分析近来得到了广泛关注[5-7]。

在微多普勒效应分析中需要得到人运动情况下的雷达回波。2002年Geisheimer应用运动捕捉(Motion capture)技术对于人行走的雷达回波进行了模拟,并分别讨论了每个部位的回波特征[5]。V.Chen基于欧拉角方法,对人运动回波的合成和微多普勒特征提取进行了研究[6]。由于条件限制或者出于研究一种不受具体设备影响的通用模型的考虑,在研究人体运动的微多普勒特征以及应用时首先需要解决的问题是人运动情况下的雷达回波模拟。在已有的工作中,对回波的仿真主要是基于运动捕捉技术和欧拉角方法,在仿真过程中需要复杂的矩阵运算,这种三维坐标不便于用来计算旋转关节的位置,还可能出现万向节死锁问题[8],因此本文针对一个常用的运动模型特点,提出利用四元数方法简化计算过程。

1 人模型与四元数概念

1.1 人体典型模型

本文应用的人体模型如图1所示,该模型由12个不同的部件组成,可以比较准确地描述人的运动情况。其中头部是一个球体,而其它部件用圆柱或者椭球体表示。

图1 人体组合模型Fig.1 Human model

在人行走过程中的任意时刻,各个部分的运动都可以被分解为相对于父关节的转动和相对于人体参考点的平动。仿真中各部位运动幅度(转动角度和平移大小)来源于文献[9],下文工作中应用这些数据计算了人运动时各个部件的实时位置。

1.2 四元数概念

四元数(Quaternion)理论是数学家Hamilton于1843年首先提出的。其形式为q=w+x i+y j+z k,因此四元数概念是复数的一个推广。其中w为实部或者标量,i、j、k为虚数部分。假设有两个四元数:

那么它们共轭和相乘的定义如下:

2 应用四元数方法计算旋转位置

在旋转位置计算中,一般应用欧拉角的方法,即令目标依次绕三个坐标轴旋转后得到旋转结果,这一方法比较直观,但存在被称为万向节死锁的问题。在处理人运动时,这一问题可能会导致出现角度插值后运动不自然甚至反向运动的情况。另外欧拉角方法需要复杂的矩阵运算。因此本文利用了另外一种旋转算法——四元数方法。和欧拉角方法相比,四元数方法解决旋转问题非常有效[8],尤其是在分析人运动时由于获得的数据量不足,需要进行插值的情况下。

2.1 用四元数表示旋转

四元数具有一个基本性质:假设经过原点的向量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)就可以对人运动的实时位置进行计算。

2.2 行人组合部件位置计算

在对人行走时各组合部件位置的计算中,首先需要考虑坐标系的转换。全局坐标系是以人的腰部为原点建立的,并以该点为计算基准点。而各部件旋转位置的计算是在以该部分父关节点为原点的局部坐标系下进行,例如小腿部分的旋转计算需要在膝关节为原点的坐标系下进行。因此计算中首先需要获得父关节在运动时刻的位置,即位置的计算是由基准点开始向远端部位顺序进行的。然后由式(3)、(4)计算运动过程中各部件的位置。在计算中各部件的旋转轴所对应的式(3)中的四元数V是由经过父关节的垂直轴来构造的,例如在图1的情况下,小腿部分的旋转轴是x轴,旋转对应的四元数R则由旋转部分的尺寸来确定。

图3是应用四元数旋转式(1)～式(4)计算各部件位置数据后得到的人行走时脚踝关节速度曲线,该结果与文献[6]经欧拉角方法得到的结果一致。

图3 人行走时脚踝关节的运动速度曲线Fig.3 Velocity of ankle when human walks

3 人行走回波信号仿真与分析

前面工作给出了一个常用人体运动模型,并且利用四元数方法计算了不同时刻行人各部件的位置。本节将基于这一结果进行人行走回波信号的仿真。

3.1 人体组合部件形状选择

在雷达方程里用雷达散射截面(RCS)来描述人体各个部分的散射,RCS大小是与目标形状、入射角度等相关的。图4给出了三种不同形状的物体的RCS曲线,可以看到和椭球体相比,圆柱体在电磁波入射角度变化时其RCS变化要剧烈的多。因此当接近或者刚好为后向入射时,圆柱体结构的RCS会产生激烈的变化,从而导致仿真雷达回波在该入射角度以及附近产生峰值,反映在仿真结果的时频图中会出现类似于冲激函数造成的频域扩散,这与实际情况不符合;因此在本文计算中人体各部件的建模使用椭球体,文献[11]指出应用该结构仿真结果更接近实际测量值。

图4 球体、圆柱体、椭球体RCS随入射角度变化曲线Fig.4 RCS of sphere,ellipsoid and cylinder

当椭球体横截面为圆时,其RCS可以由式(5)来描述[6]:

3.2 人行走回波仿真

电磁波经过人体部件时会发生散射。本文对于行人目标的回波计算采用走停模型来处理,这一模型假设在同一个脉冲周期内目标静止,在相邻脉冲之间目标运动,可以比较精确地描述低速运动目标的散射回波。

在计算人体运动的散射回波时,可将人体等效为如图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

3.3 回波数据处理

由于合成的行人雷达回波是非平稳信号,从中难以直接得到有价值的信息,因此根据非平稳信号的特点对其进行时频分析处理。本文采用自适应核函数的方法[12],时间窗为长度为31的 Hamming窗,时频谱图如图6(a)所示。从图6(a)可以看出人行走中身体不同部分所对应的信号,其中散射信号最强的位置对应于躯干及相邻的大腿等部件的运动,而多普勒频移最大的位置对应于脚的运动,这和实际情况是一致的:实际运动中人体各部件的速度是不一致的,躯干速度相对稳定,而越到身体边缘处的运动速度变化越剧烈,人行走时脚的运动速度可能远大于人体躯干的速度。为了验证本文合成信号的可靠性,图6(b)给出了利用运动捕捉技术合成雷达回波并进行时频分析处理后的结果[5],可见两者之间具有较好的一致性。

图6 人行走雷达回波的时频谱Fig.6 Spectrum of human walking

4 结论

本文提出利用四元数来计算人体运动中伴随有微运动的各组成部件位置,并用计算结果进行人运动雷达回波合成的方法。这一方法可以有效地减少位置计算中的矩阵运算,避免出现万向节死锁等情况。应用该方法仿真获得的人体运动散射回波与国外实测数据对比结果表明了这一方法的有效性。

对人运动的多普勒谱图进行分析,可从中识别出人体组成部件的微多普勒效应。下一步工作将对合成运动回波的微多普勒现象进行具体分析,提取合理的微多普勒特征应用于人目标身份等方面。

[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.