基于调频连续波雷达的物体运动状态实时检测算法研究*

屈奎 张荣福† 肖鹏程

1) (上海理工大学,光电信息与计算机工程学院,上海 200093)

2) (复旦大学,专用集成电路与系统国家重点实验室,上海 201203)

微波雷达依靠非接触、响应速度快、对自然环境的适应性强等特点在物体运动状态检测中的应用越来越广泛.常用的调频连续波雷达运动检测算法基于差拍信号频谱的峰值估计,存在计算量大,抗干扰能力差等缺点.本文通过对运动物体的差拍信号做特定频率的离散傅里叶变换,将变换后的实部和虚部在互相垂直的两个方向上进行叠加,其合成轨迹近似为椭圆,求出各轨迹点的相位即可还原物体的运动状态.该算法无需对每个调频周期的拍信号做频谱分析,时间复杂度较低.静止物体的拍信号被处理成了固定的直流信号,对运动物体的测量不造成影响,具有抗静止物体干扰的能力.在雷达中心频率为24 GHz,带宽为0.15 GHz 的条件下对算法进行了验证,位移测量精度达到0.27 mm,以500 mm 作为位移的测量范围,线性度达到0.05%.速度的测量精度为1.11 mm/s.

1 引言

位移和速度是表征物体运动状态的基本物理量,运动状态的监测可以应用在精密制造、安防、自动驾驶、健康监控等方面.特别是现代社会向智能化、万物互联快速发展的过程中,对各种物体如机器人、无人机、液压杆、传送带等的精确控制提出了越来越高的要求,对物体运动状态的高精度检测是精确控制的必要前提.

运动状态检测技术可以采用接触和非接触两种方式,接触式是传感器与被测物有实体上的连接,该方式受到很多限制,某些情况下不宜采用,比如当被测物与监测仪器距离较远或者被测物对轻量化要求较高.还有些情况无法采用接触式,比如自动驾驶中对车道上其他车辆的监测.非接触测量方式较为灵活,对被测物的运动状态基本没有影响,可以在被测物无感知的情况下进行监测,常用的方法有光学[1]、超声[2]、激光[3]和微波雷达等,光学、超声、激光方法的共同缺点是测量精度受烟、尘、雾、雨等环境影响较大,有时甚至无法测量,微波雷达受环境的影响很小,是一种实际应用价值很强的运动测量方式.

微波雷达常见的调制方式有单频[4,5]、脉冲[6,7],调频[8,9]等,其中单频雷达发射单一频率的电磁波,通过计算回波的相位可以获得很高的位移测量精度,但在获得物体的绝对距离以及多目标测量中存在缺陷,虽然文献[10,11]报道了对这些缺陷的弥补,但付出了提高硬件和计算复杂度的代价.脉冲雷达要想获得较高的测距精度要求有较窄的脉冲宽度和极高的时间分辨能力[12],这些都对硬件提出了很高的要求.调频连续波(Frequency modulated continuous wave,FMCW)雷达是对发射波的频率进行调制,最简单的一种调频方式是线性调频,也是较为常用的调频雷达制式.调频雷达相比单频雷达其回波携带的信息更加丰富,相比脉冲雷达,其信号收发同时工作,不存在距离盲区,且同样的探测距离其发射功率峰值较低.

使用FMCW 雷达进行精确位移测量,一般核心步骤是使用频谱的各种估计方法得到待测物体的拍信号对应的精确频率.大致过程是在每个调频周期内,对差拍信号的幅度谱通过峰值查找得到目标谱线的粗略位置,再用频谱细化方法进行高精度修正.齐国清在研究油罐液位的测量中通过对差拍信号做离散傅里叶变换(Discrete Fourier transform,DFT)得到相位,为解决相位随距离变化的周期性带来的距离初值的模糊性,将一个调频带宽的信号分成两段,分别求出各段的相位,两者的差值用来消除距离模糊[13,14];Zwick 研究组[15−17]在对拍信号做DFT 的基础上,采用插值法,通过离散频率的振幅比值确定拍频频率,该组还尝试了使用调频Z 变换(Chirp Z-transform,CZT)进行较为精细的频率测量;Pohl 研究组[18−20]开发了中心频率为80 GHz,带宽达到20 GHz 以上的超宽带FMCW雷达系统,大的带宽带来很高的距离分辨率,但也意味着更复杂的硬件结构.以上研究组针对的测量场景都较为简单,回避了多径干扰,把问题归结为求解单一正弦信号的频率,在具有多散射物的场景中并不适用.另外,各国以国际通信联盟制定的ISM(Industrial scientific medical)频段为基础,对无线电频段的使用进行了严格的限制,在实际应用中无法通过任意提高带宽达到提高距离分辨率的目的.

本文所提算法仅对运动物体的差拍信号做特定频率的DFT,将变换结果中的虚部和实部在互相垂直的两个方向上叠加,叠加后的轨迹近似为椭圆,求出各轨迹点的相位即可还原物体的运动状态.在雷达中心频率为24 GHz,带宽为0.15 GHz的条件下对理论进行了实验验证.

2 FMCW 雷达的差拍信号

FMCW 雷达在一个调频周期T内发射信号的波形函数为

式中,f0为起始频率,Bw为调频带宽,Φ0为初相.发射波的频率随时间的变化关系如图1 所示.

图1 发射波的频率随时间的变化Fig.1.the time-frequency domain of transmitted signal.

常用FMCW 雷达收发天线的间距远小于与被测物的间距,因此可将收发天线看成在同一个位置.与天线相距为R的物体,其回波信号相比发射信号时间延迟了td2R/C,C是电磁波在介质中的传播速度,其回波函数为

回波和发射波的混频信号为

由于td≪T,可以把时间延迟量td的平方项部分省略,得到:

(4)式即是与天线相距为R的物体对应的差拍信号,式中波函数的振幅均做了归一化处理.

3 算法原理

调频周期T一般为毫秒甚至微秒级,考虑这样一种情况,在T的时间内物体的位移很小,以致可以忽略,此情况下,可以认为物体与雷达距离R不变,其差拍信号用(4)式表示.

对一个调频周期内的差拍信号做特定频率的N点采样DFT,表达式为

仅从(5)式很难看出结果的规律性,(6)式是其积分形式,采样数越密集,(5)式与(6)式的结果越接近.(5)式结果的规律性可以由(6)式的结果来反应:

将(4)式代入(6)式,可得:

R ek和 I mk分别是F(k) 的实部和虚部.令λC/(f0+Bw/2),λ等于发射波的中心波长.(7)式、(8)式可记为

(9)式、(10)式结构相似,均明显可分为三个因子相乘,即它们的第三项随R作周期性变化,变化周期为λ/2.当物体的位移不大于第三项的周期,即|∆R|≤λ/2 时,如果第一项和第二项乘积的相对变化很小,在该运动过程中就可以近似用一个常数替代它们.

将(11)式分别代入到 R ek和 I mk的前两项乘积中,可得:

当|∆R|≤λ/2 时,所引起的(12)式和(13)式的相对变化量满足如下关系:

通过分析(14)式和(15)式右边表达式的变化规律,可以知道,当 ∆Rλ/2,k1,δ−0.5 时,(14)式右边的值最大; ∆Rλ/2,k1,δ0.5时,(15)式右边的值最大.

以在实际应用中常见的24 GHz 雷达系统为例来观察(14)式、(15)式的具体值.具体参数可以取λ1.25 cm,Bw0.15 GHz,C3×108m/s,估计出(14)式和(15)式的具体值为

(16)式表示在给定波长和带宽条件下,∆R在λ/2 的范围内变化所引起的,的相对变化量很小,如果把它们看成常数,则(9)式、(10)式可以近似写为

实际测量环境中,可能会有静止物体、杂波以及其他运动物体的回波干扰.简单起见,可以只考虑静止物体、具有固定偏移量的干扰及频率稳定的杂波,对于这些干扰的拍信号做特定频率的DFT后是一个直流信号,相当于在(17)式、(18)式的基础上增加一个固定的直流分量 R e0和 I m0.

显然,合成轨迹是一个椭圆,标准形式为

通过椭圆拟合算法,可以得到椭圆参数,也即得到 R e0,I m0,Ak,Bk的值.

为方便描述,以下将这种拍信号在特定频率下DFT 分量垂直合成轨迹的椭圆化算法(Ellipse algorithm for the vertical synthesis trajectory of the DFT component of the beat signal at a specific frequency)简记为ETBF 算法.

考虑仅有单个运动物体的情况,特定离散频率系数k的确定分两种情况讨论.第一种,待测物的初始位置R0已知,则根据k的选取原则,kround(2R0Bw/c),符号round()表示对括号里面的量做四舍五入.初始的k值确定后,可以根据待测物体的位移R′对k值进行更新.更新公式为

第二种情况,待测物初始位置未知.对第一个调频周期内的拍信号做幅度谱的峰值检测,在这些峰值对应的离散频率下将拍信号用ETBF 算法处理,静止物体的信号被ETBF 算法处理后是一个固定值,运动物体的信号在ETBF 算法处理下会形成椭圆化的合成轨迹,而且当离散频率最接近被测物的差拍信号频率时,由于此时(6)式积分中的两个乘积项具有最强的相关性,得到的椭圆幅度将会最大,所以幅度最大的椭圆对应的离散频率序数即为初始的k值,再用频谱细化方法得到较为精确的频率finitial,根据拍信号频率和距离的关系可以得到初始位置R0CTfinitial/(2Bw).同样,采用(22)式对k值进行更新.此情况仅在初始位置的确定中使用了常用的频谱估计算法.

4 相位解缠(phase unwrapping)算法和运动状态的计算

得到拟合椭圆的参数后,即可通过(19)式和(20)两式计算各轨迹点对应相位φ,即,易知:

根据(24)式计算反正切得到的相位φ′被限制在[ 0,2π).

分析(19)式、(20)式可知,当R增大时,Rek和 I mk的合成轨迹点做逆时针运动,反之做顺时针运动.相位φ与R成正比,所以,当合成轨迹点做逆时针运动,意味着相位增大,反之减小.考察轨迹点逆时针越过图2 所示的横坐标轴的正半轴时的情况,此时实际相位增大,但直接根据(24)式计算出的相位减小;当轨迹点做顺时针运动越过横坐标的正半轴时,实际相位减小,但根据(24)式计算出的相位增大.所以,直接根据(24)式的计算结果无法反应真实的相位变化,需要通过相位解缠算法解决这一问题.而且由于初始时刻的相位值无法获得,最终得到的是相对初始时刻的相对相位.

图2 变量R ek和I mk的合成轨迹示意图.红色实心点处在横坐标轴的正半轴上,代表相位φ=2jπ,j∈Z 的位置Fig.2.The synthetic ellipse trajectory diagram of variables Rek and I mk.The red solid point is on the positive half axis of abscissa,which represents the position of φ=2jπ,j∈Z.

将第n个由(24)式计算出的相位采样值记为φn,解缠后的相对相位记为,则:

jn可以称为相对相位级数,n=1 时,j10,φ1是n=1 时的φn.

解缠绕的过程就是确定jn的过程.定义 ∆φnφn −φn−1,根据∆φn的值可以判断φn相对φn−1是否越过了横坐标的正半轴.

1) −π<∆φn <π

后一采样点相对于前一采样点在椭圆轨迹上没有越过横坐标的正半轴,当前相对相位级数不变

2) π<∆φn <2π

此情况代表轨迹点做顺时针运动且越过了横坐标的正半轴,相位减小了 2π.所以当前相对相位级数

3) −2π<∆φn <−π

此情况代表轨迹点做逆时针运动且越过了横坐标的正半轴,相位增大了 2π.所以当前相对相位级数

由相对相位可求出其对应的位移为

位移的符号上加了撇号以示和绝对距离相区分.

当位移的测量比较准确时,可以直接对时间求导得到速度

对于离散值,(30)式的数值计算表达式为

实际是用 ∆t时间内的平均速度代替瞬时速度,具体计算时,可根据实际情况选取合适的 ∆t值.

5 与其他算法的比较

时间复杂度是对算法的消耗时间进行度量,由于算法消耗时间与工作量成正比,所以也等价于对工作量进行度量.通过大O符号表示法[23]给出计算工作量的渐进函数,可以观察不同算法在计算工作量上的本质差异.

采用表1 中其他文献的算法时,杂波会对待测物的测量造成干扰,杂波越强烈,干扰越大,测量误差越大,如果杂波强度接近或超过待测物的信号强度时,甚至导致待测物无法检出,所以文献采用的测量场景都较为简单,回避了其他散射物的杂波干扰,实际测量环境很难满足这种严格的条件,算法的适应性受到很大的限制.

表1 与其他算法的比较Table 1.Comparison of this work with other methods.

本文所提出的算法回避了频谱的峰值估计,只需计算一个特定频率的DFT,时间复杂度显著低于其他文献中的算法.由于静止物体的拍信号被处理成了固定的直流信号,对运动物体的测量不造成影响,所以具有抗静止物体干扰的能力.由于计算的是相对相位,限制了测量内容为相对运动.

6 实 验

6.1 硬件系统

得益于集成电路制造技术的飞速发展,雷达的关键部件已实现芯片化,大大简化了小型雷达系统的设计和制造难度.实验所用雷达系统主要由ADI 公司的ADF4158 芯片和Infineon 公司的BGT24MTR 芯片组成.ADF4158 芯片产生调制信号,控制BGT24MTR 芯片中的压控振荡器产生24 GHz 调频波形,信号经过功率放大器(Power amplifier,PA)放大后,采用微带天线发送和接收反射回来的信号.接收到的信号经过低噪放大器(Low noise amplifier,LNA)放大滤波后进入混频器与发射信号的一部分进行混频,产生的差拍信号由NI 公司的PXI-4461 采集卡进行数据采集,然后输入到计算机通过MATLAB 进行信号处理.硬件采用了单通道输入和输出.雷达前端如图3 所示,实验系统的结构如图4 所示.

图3 雷达前端 (a)正面为收发天线;(b)反面为电路板Fig.3.Radar front end:(a) The front side is transceiver antenna;(b) the reverse side is circuit board.

图4 实验系统结构图Fig.4.Block diagram of the FMCW radar system.

实验中调频信号的参数为f024 GHz,Bw0.15 GHz,T=4 ms.

在校准过的步进电机上固定一块金属板作为反射板,步进电机的运动方向垂直雷达的天线面,装置如图5 所示.

图5 测量装置与载有金属板的步进电机Fig.5.Measuring set and stepper motor with metal plate.

6.2 R ek 和 I mk 采样点的合成运动轨迹

由理论公式(19)和(20)可知,λ/2 的运动距离刚好可以使 R ek和 I mk的合成运动轨迹形成一个完整的椭圆.图6 给出了用计算机仿真的待测物在距离天线800 mm,1200 mm,1600 mm,2000 mm四个起始位置上运动λ/2 时所形成的运动轨迹.仿真采用的调频信号参数和实验相同,将待测物拍信号的振幅设置为单位1,在4 个起始距离上设置拍信号振幅为0.5 的干扰物体,并添加功率为–20 dBW 的高斯白噪声.图中星号代表采样点,实线是根据这些采样点拟合出的椭圆,可以看出,这些采样点形成的轨迹和拟合椭圆符合得很好.

图6 计算机模拟在4 种不同起始距离下做 λ/2 位移时的采样点轨迹变化,起始距离分别为 (a) 800 mm;(b) 1200 mm;(c) 1600 mm;(d) 2000 mmFig.6.Computer simulation of trajectory change by these sampling points at four different starting distances,the displacement is λ/2,the starting distances are:(a) 800 mm;(b) 1200 mm;(c) 1600 mm;(d) 2000 mm.

控制步进电机带动金属板分别在距离天线800 mm,1200 mm,1600 mm,2000 mm 的距离上运动λ/2.图7 给出了实验测量出的采样点和根据这些采样点拟合出的椭圆,采样点基本落在拟合椭圆上,仿真和实验结果都与理论分析一致,说明把运动轨迹看成椭圆具有合理性.由于仿真和实验中拍信号的振幅并不一致,所以图6 和图7 的坐标范围不一致.

图7 实验测量在4 种不同起始距离下做 λ/2 位移时的采样点轨迹变化,起始距离分别为 (a) 800 mm;(b) 1200 mm;(c) 1600 mm;(d) 2000 mmFig.7.Experimental measurement of trajectory changes by these sampling points at four different starting distances,the displacement is λ/2,the starting distances are:(a) 800 mm;(b) 1200 mm;(c) 1600 mm;(d) 2000 mm.

6.3 位移的测量精度

让金属板初始时刻距离雷达天线1500 mm,控制步进电机带动金属板以10 mm/s,20 mm/s,30 mm/s 三种速度分别匀速运动300 mm,400 mm,500 mm.重复测量20 次,表2 给出了测量的统计结果.平均值的偏差处于(–0.1 mm,0.1 mm)之间,测量的最大偏差为0.27 mm.测量结果的标准差小于0.15 mm.

表2 金属板运动的测量结果Table 2.Measurement results of metal plate movement.

亚毫米精度的运动测量一般应用在机器人室内定位、机械臂控制、倒车监测等近距探测场景中,范围一般在5 m 以内,在这样的距离范围内,实验测量未发现有明显的精度变化,说明算法在近距测量场景中具有较好的精度稳定性.当测量距离进一步增大后,由于待测物的散射立体角减小而导致的信号强度降低越来越明显,测量精度会逐步下降.

6.4 位移的测量线性度

线性度可以反应系统对待测运动过程还原的准确性.从金属板以20 mm/s 的速度移动500 mm的20 组实验数据中任意选取一组,位移随时间变化的测量结果如图8(a)所示.匀速运动时,时间与位移应该是线性关系,图中同时给出了用最小二乘法拟合出的直线作为线性比较的标准.图8(a)中的子图是主图的部分放大图,从中可看出测量值和拟合直线之间的轻微偏差.图8(b)为测量偏差变化图,偏差的平均值为 2.7×10−13mm,标准差为0.22 mm,基本满足零均值的高斯分布.最大偏差为0.25 mm,以500 mm 作为测量范围,线性度达到0.05%.

图8 位移测量结果 (a)位移随时间的变化;(b)位移偏差随位移的变化Fig.8.Measurement of displacement of metal plate:(a) Change of displacement with time;(b) change of displacement deviation.

6.5 速度的测量精度

为检验速度测量的精确度,设计了如下的运动过程:首先从静止开始以 1 0 mm/s2的加速度运动3 s,此时速度达到30 mm/s,然后匀速运动8 s,再以 − 10 mm/s2的加速度运动3 s,此时速度为零,静止5 s,然后以 − 5 mm/s2加速度运动6 s.整个过程可分为5 个阶段,分别为匀加速、匀速、匀减速、静止及一个反向匀加速阶段.通过对步进电机编程让其带动金属板按照设定运动.

由于采用(31)式计算速度,需要先得到运动过程的位移.图9(a)同时给出了位移的设定值和实际测量值,两者的曲线变化几乎完全一致.图9(b)是测量值与设定值的差值变化,平均偏差为0.02 mm,最大偏差δR.max为0.29 mm,比线性运动的偏差稍大.偏差较大的部分出现在加速过程中,可能是步进电机的加速不稳定导致.

图9 (a)测量位移和设定位移的比较;(b)测量位移的偏差变化Fig.9.(a) Change of measured displacement and set displacement with time;(b) deviation of measured displacement.

测量速度和设定速度随时间的变化如图10(a)所示.其中计算平均速度所选取的时间间隔 ∆t为0.3 s.根据误差的传递公式,测量速度的理论偏差σ应该满足.图10(b)是测量值与设定值的差值变化,最大偏差为1.11 mm/s,标准差为0.31 mm/s,符合误差分析结果.通过适当增加时间间隔可以提高速度的测量精度,但会带来速度更新率的下降,实际应用中,可根据具体情况灵活设置.

图10 (a)速度测量值和速度设定值的比较;(b)速度测量值的偏差变化Fig.10.(a) Change of measured speed and set speed with time;(b) deviation of measured speed.

7 总结

本文所提出的算法可以在单通道接收和输出的硬件条件下实现,单通道的输出无需考虑正交解调输出时的信号失衡问题,且硬件结构更简单.对运动物体的差拍信号做特定频率的离散傅里叶变换,避免了使用复杂方法进行频率估算,减少了计算量,易于实时测量.测量中可以隔离其他静止物体的干扰,直接提取运动物体的信号,降低了对环境的要求.由于算法得到是相对相位,限制了测量内容为相对运动.在一些需要测量相对位移,比如机械的振动频率,生命信号探测,机械手操控等领域有很大的应用潜力.

在雷达中心频率为24 GHz,带宽为0.15 GHz的条件下对算法进行了验证.位移测量精度达到0.27 mm,以500 mm 作为位移的测量范围,线性度达到0.05%.速度的测量精度为1.11 mm/s.实验采用了符合ISM 要求的频段,更符合实际应用场景,但中心频率和带宽低于文献[16,19−22]中的设置,高的频率和带宽可以显著提高位移的测量精度,由于硬件条件不一致,不能把本实验的测量精度简单的与上述文献比较.研发出更高频段的雷达硬件对本算法进行验证及实现绝对距离的测量是下一步要做的工作.