基于连续波雷达的人体目标微多普勒频率估计方法

丁一鹏，柳润金，许雪梅

(中南大学物理与电子学院，湖南长沙，410083)

连续波多普勒雷达在人体目标探测方面[1-2]具有独特优势，它可以抑制静止杂波，穿透墙壁等障碍物，且检测条件不受天气、光线影响，因此被广泛应用于反恐、边境安防、灾难救援等重要领域。微多普勒频率作为雷达目标所具有的独特特征，可用于目标运动状态的识别[3-4]。因此，微多普勒频率的精确估计对提高雷达探测和目标识别的准确性至关重要。为了更精确地识别感兴趣目标的散射部位(例如手臂和腿)的运动，需要先提取这些部位的基本回波分量。然而，人体柔性关节和复杂的运动模式使得其雷达回波动态范围大，时频谱带相互交叠，人体各部位的瞬时微多普勒频率很难被精确估计。此外，人体各部位微多普勒回波为非平稳随机信号，如何准确拟合回波局部特性并精确估计目标微多普勒频率，是检测算法设计中需要考虑的问题。

学者们针对微多普勒频率估计技术进行了大量研究。最常用的分析方法是线性时频分析技术，如短时傅里叶变换[5-6]和S变换[7]，这些技术在分析处理多分量的信号时可免受交叉项的干扰且易于实现，但时频分辨率的折中处理使其在雷达目标探测领域的应用受到限制。传统的非参数信号分离算法如经验模式分解(EMD)[8-9]、主成分分析(PCA)[10]和Hilbert-Huang 变换[11]等都是基于回波频率特性而不是目标散射特性，这会给分离后的频率估计结果带来很大的偏差。参数信号分离算法[12]的参数搜索过程则通常涉及大量计算。稀疏分解算法[13]近年来在语音信号、地震信号和雷达信号的分离和去噪研究中也取得了较大进展，然而，算法所需原子库的完备性和原子匹配会对算法的收敛性造成影响。借鉴图论中的边缘检测算法思想，利用维特比算法[14]可以对瞬时频率轨迹进行跟踪，但该算法只能准确估计单分量信号频率。改进的维特比算法[15]增加了局部斜率差惩罚函数，然而，在多散射点条件下该算法仍表现出稳健性不足、估计精度较低的特点。结合航迹关联思想，通过设置隶属度惩罚函数[16]可以提高改进型维特比算法对多分量信号频率估计的精度，但在人体目标特定散射部位的微多普勒频率估计中会出现路径分叉问题。

本文作者基于三次动态指数平滑预测方法获取时频分布中最佳路径的非线性变化趋势，将其表示为新型惩罚函数的形式并对传统维特比算法进行改进，从而增强最佳路径的识别能力；为了更好地拟合人体散射部位的局部特性，借鉴Hough变换的频率估计思想[17]，用动态平滑系数代替传统指数平滑法的静态平滑系数；在人体微多普勒频率估计的应用场景中，该方法以历史路径为参考，并通过动态确定候选路径的搜索范围以提高最佳路径的寻找效率。

1 维特比瞬时频率估计算法

维特比瞬时频率估计算法的基本思想来源于图论和数字图像处理中的边缘跟踪问题[14]。对信号进行时频分析后可以得到二维时频分布，利用维特比算法从时频分布中搜索最佳路径从而提取信号瞬时频率主要基于以下2个假设：1)每个时刻的瞬时频率对应的时频点的幅值要尽可能大；2)两相邻时刻的瞬时频率变化不会太大，瞬时频率曲线比较平滑。

对信号进行时频分析后，可得到N行M列的时频分布，将时频分布数据中时刻m的时频幅值T排列成非递增序列：

式中：j= 1,2,…,N；m= 1,2,…,M；kj为m时刻的第j个时频点。

采用如下算式定义惩罚函数f(·)：

考虑分别位于相邻时刻m和m+ 1的2个时频点(m,kl)和(m,kj)， 其 中，kl= 1,2,…,N且kj=1,2,…,N，则惩罚函数g(·)可定义为

式中：c1为惩罚因子，可以反映g(·)惩罚函数值对总惩罚值的贡献度；Δ1为2个连续点之间瞬时频率变化的最大期望值。路径的总惩罚值计算方法如下：

式中：p为总惩罚值。最终的最优路径即为使惩罚值最小的时频曲线Kmin：

为了对多分量信号瞬时频率进行估计且增强算法稳健性，改进的维特比算法[16]设置了如下隶属度惩罚函数：

式中：f′k,i为通过多点多步预测的瞬时频率；fk,i为需要进行隶属度判别的瞬时频率；u为常数。则式(5)所示路径总惩罚值可表示为

传统维特比算法在单分量信号和频谱混叠不严重的多分量信号瞬时频率估计方面表现良好，然而，传统维特比算法是基于时频分布中脊线的幅值定义f(·)惩罚函数，其所定义的g(·)惩罚函数仅通过限制2个连续采样点绝对频率变化值来提取路径最短的平滑脊线，而不考虑脊线的变化趋势。由于采样间隔太小，惩罚函数的约束是较弱的局部约束[18]，在面临时频分量交叠处的路径搜索时，会产生路径分叉问题。传统维特比算法路径分叉示意图如图1 所示，其中Im(t)和In(t)分别为时频域交叠的2 个分量信号，gm(k4,k5)为由Im(t)分量上的点4 和5 的频率变化而引起的g(·)惩罚值；gn(k4,k8)为由In(t)分量上的点4和8的频率变化而引起的g(·)惩罚值。对于待估分量Im(t)，点1至点7为最佳路径。 然而， 当gm(k4,k5)≈gn(k4,k8)且点8的幅值比点5的幅值更大时，传统维特比算法会将干扰分量In(t)上的点8 误选为待估分量Im(t)的路径，后续的路径搜索可能完全偏离待估分量，从而将点1至点4和8至点10误搜索为最佳路径，产生路径分叉问题。在隶属度惩罚函数[14]的设置过程中，所采用的预测方法是选择局部若干时刻的瞬时频率对当前时刻瞬时频率进行线性预测而不是基于全部历史数据，故该方法不适用于非平稳信号序列的预测。此外，在整个时频平面范围内搜索候选路径会使算法计算量增大。

为解决上述问题，提高微多普勒频率估计精度和最佳路径的搜索效率，本文采用一种基于三次动态指数平滑预测的具有高搜索效率的新型维特比算法用于估计人体目标微多普勒频率。

2 基于三次动态指数平滑预测的新型维特比算法用于人体微多普勒频率估计

2.1 人体步态微多普勒回波模型

为了保持实验模型所固有的动态特性，采用基于简单运动学方法设计的经典Boulic人体模型[19]进行人体微多普勒回波分析。当连续波多普勒雷达发射天线发射电磁载波时，运动的目标对载波信号产生频率调制，接收天线接收到的人体基带回波模型r(t)可用式(9)表示，Boulic 模型将人体目标分为若干散射部分，回波信号r(t)即人体各散射部位(如头部、躯干和腿部)的回波之和。

式中：Nh为人体散射部位的数量；c和fc分别为光速和雷达载波频率；θi为目标运动速度与雷达视线方向的夹角；Ri(0)为第i个散射中心与雷达的初始视距；vi(t)为第i个散射部位运动速度。第i个散射部位的微多普勒频率fdi(t)可以表示为

微多普勒频率是由人体运动模式和柔性关节的复杂相互作用对雷达回波产生的非线性频率，在时频谱上表现为整体多普勒频率周围的旁瓣和边带效应。

2.2 微多普勒频率估计算法

人体正常行走时具有大惯量、低速度的特点，且其微多普勒频率曲线为类正弦曲线。基于此，本文考虑将可靠性和实时性高的三次指数平滑预测技术应用于人体目标微多普勒频率估计。下面介绍基于动态指数平滑预测模型的新型维特比算法的原理和构建步骤。

三次指数平滑预测是一种序列分析和预测方法，能将非线性时间数据序列的数量差异抽象化，对序列进行修匀从而消除不规则和随机扰动，显示出预测对象变动的基本趋势。利用传统维特比算法在短时傅里叶变换后的时频分布中搜索初始路径，并计算得到三次指数平滑预测的历史路径平滑值：

获取历史路径的平滑值后，可得三次指数平滑法的预测系数λ，γ和η的计算公式[20]如下：

为增强传统维特比算法惩罚函数的最佳路径选择能力，减少频率模糊区域的微多普勒瞬时频率分叉问题，将新型维特比算法的惩罚函数μ(·)定义如下：

式中：ε为三次指数平滑基于全部历史数据平滑值的短期预测值，且ε=λ+γ+η；Δ2为设定的频率阈值，在人体步态微多普勒频率估计的应用场景中，Δ2一般设置为3～5 Hz；c2为惩罚因子，用于表征μ(·)惩罚函数对总惩罚值的贡献度，当预测值与当前估计值距离较小时，惩罚值较小，该预测值所在路径属于最佳路径的可能性较大；若路径在频率模糊区域向其他路径分叉，远远偏离预测路径，则将产生较大的μ(·)惩罚值，使之被选为最佳路径的可能性减小。因此，这样定义的罚函数既考虑了时频脊线的变化趋势，也加强了确定最佳路径需要的整体约束，可对当前搜索路径进行修正。

将惩罚函数μ(mi,kj)代入式(5)，有

以p(mi+1,kj)为总惩罚值对时频分布进行搜索，则以αml为模型参数的微多普勒频率估计模型F(αml,t)表达式如下：

下面根据Hough变换频率估计的思想[17]确定动态平滑系数αml。调整平滑系数αml，以F(αml,t)为频率估计模型对雷达回波信号r(t)进行解调，并将傅里叶变换应用于解调后的信号，得到如下表达式：

在理想情况下，若待估计分量的微多普勒频率fp(t)≈F(αml,t)，则解调后该分量将成为恒定信号，其能量将以脉冲的形式在频域中累积。对于估计模型以外的其他分量，其能量则不会在频域中收敛，因此，当S(f)的信号能量收敛到基带的最大值时，可以得到感兴趣的人体散射部位的微多普勒频率fp'(t)≈F(α′ml,t)，α′ml即为最佳模型参数。基于指数平滑预测模型的新型维特比算法构建流程如图2所示。

将窄带滤波器应用于具有最佳模型参数的S(α′ml,f)，将当前所估计分量从原始信号中分离出来后，可对其他部位的微多普勒频率进行估计：

式中：s(t)为被解调分量；FB(f)为滤波器的频率响应；IFFT为逆傅里叶算子；sp(t)为由F(α′ml,t)分离的分量。

在路径搜索过程中，并不需要搜索整个时频平面上的所有点，只有在能量收敛区域附近的搜索才有价值，其他搜索点会耗费更多处理时间。由于正常人体行走时的运动为惯性运动，其雷达回波具有局部静态特性，且两相邻时刻的瞬时频率变化不会非常剧烈，瞬时频率曲线应该比较平滑。因此，通过式(6)确定某一时刻的历史路径Kpi后，以[Kpi-δ,Kpi+δ]mi作为参考间隔，[Kpi-δ,Kpi+δ]mi+1为搜索区间来进行候选路径的动态搜索，参考区间和搜索区间长度均为2δ+ 1，若时频分布中的极差为R，则δ可估计为R/2。在一般情况下，人体正常行走时微多普勒频率估计的δ不超过5，则在确定从mi到mi+1的路径时，此搜索方法只需要不超过(2δ+ 1)2即121次搜索。

新型维特比算法用于人体特定散射部位微多普勒频率估计的步骤如下：在三次动态指数平滑预测的基础上获取最佳路径的非线性变化趋势，并依据g(·)惩罚函数形式定义新惩罚函数μ(·)，对当前搜索的最佳路径进行修正从而提高传统惩罚函数的最佳路径识别能力，抑制传统维特比算法在频率模糊区域的路径分叉，同时动态调整候选路径的搜索范围，提高最优路径的寻找效率。

3 仿真实验与分析

为了验证所提出的基于三次动态指数平滑预测的维特比算法用于人体目标微多普勒频率估计的有效性和优越性，采用经典的Boulic人体模型作为分析标准，通过MATLAB 软件进行仿真实验，对比不同算法的相应估计结果，并进行分析。所应用的Boulic 模型中的人体不同部位的尺寸见表1，这些数据是通过对1 个身高约1.78 m 的男性目标进行实验测量获得的。雷达探测场景示意图如图3(a)所示，连续波雷达所使用的载波频率为2.4 GHz，发射单元和接收单元间隔距离约为6.5 cm。目标以0.5 m/s 的恒定径向速度向雷达走去，每0.005 s采集1次回波信号，接收回波的频谱如图3(b)所示。

表1 人体目标各部位尺寸Table 1 Dimensions of different parts for human targets

在基于微多普勒频率对目标进行识别和分类的实际应用中，目标四肢的微多普勒信息是区分人类和其他动物的典型特征。由于四肢运动的动态范围大且其微多普勒频率形状不规则，最难精确估计，下面分别采用基于短时傅里叶变换的峰值检测法、传统维特比算法和本文所提出算法依次对人体左腿和左臂的微多普勒频率进行估计，结果分别如图4和图5所示。在传统算法的频率估计过程中，其他部位的雷达回波会对估计结果产生较大干扰，产生路径分叉问题，故在图4 和图5中用细实线标示出其他干扰散射部位的实际微多普勒频率以用于对比。

由图4(a)和图5(a)可见，无论是否处于频谱交叠区，峰值检测法所估计的左腿微多普勒瞬时频率曲线都显示出来自其他分量的干扰。这是因为峰值检测法按能量脊线提取频率，而干扰分量的能量大于感兴趣部位分量的能量，使得所估计的频率曲线往其他分量的频率曲线跳变。

由图4(b)和图5(b)可见：与峰值检测法相比，传统维特比算法在一定程度上减小了微多普勒频率估计误差。然而，该方法在频率模糊区域仍存在路径分叉问题。这是因为在这些区域能量分布密集，传统维特比算法的f(·)惩罚函数值相近，而g(·)通过仅限制2个连续采样点绝对频率变化值来提取路径最短的平滑脊线，不考虑脊线的变化趋势，导致维特比算法搜索的路径分叉至干扰分量频率曲线上惩罚值最小的点。

由图4(c)和图5(c)可见：新型维特比算法所得结果与实际目标微多普勒频率几乎相同。不同算法所得手臂微多普勒频率的均方根误差和处理时间对比见表2。由表2 可见，与峰值检测法和传统维特比算法相比，本文所提方法的均方根误差显著降低。

表2 不同算法所得微多普勒频率的均方根误差及处理时间对比Table 2 Comparison of micro-Doppler frequency root mean square error and processing time of different algorithms

4 结论

1)基于三次指数平滑预测获取传统维特比算法历史路径的变化趋势，并定义新型维特比算法的惩罚函数μ(·)，可增强传统维特比算法惩罚函数的最佳路径选择能力，抑制路径分叉，提高连续波雷达频率模糊区域的微多普勒瞬时频率估计精度。

2)采用Hough变换频率估计思想，以F(αml,t)为频率估计模型对雷达回波信号进行解调，依据雷达回波的局部特性动态调整平滑系数αml；同时，动态调整维特比算法候选路径的搜索范围，提高了最优路径的寻找效率。

3)与基于短时傅里叶变换的峰值检测算法和传统维特比算法相比，基于三次指数平滑预测的新型维特比算法得到的人体正常行走时左臂微多普勒频率估计值的均方根误差显著降低。