CN109031277A - 一种穿墙雷达多目标图像域稳健跟踪方法 - Google Patents

Abstract

本发明公开了一种穿墙雷达多目标图像域稳健跟踪方法，涉及穿墙雷达技术，特别涉及穿墙雷达的多目标检测跟踪技术。本发明提供适用于穿墙雷达成像后的多目标的图像域跟踪算法。在一定允许误差内快速成像，得到较高质量成像结果；基于PCF加权处理的图像增强算法能有效地抑制栅瓣的影响，增强目标；基于双轴投影的目标模板创建算法能有效地提取目标的特征信息，为跟踪做准备；基于均值漂移和卡尔曼滤波的图像域跟踪算法能有效地应对目标图像扩展和形变的问题，也能轻松应对目标遮挡问题，跟踪精度高。因此，本发明具有实时性好，跟踪精度高的优势，可以直接应用到穿墙雷达装备中。

Description

一种穿墙雷达多目标图像域稳健跟踪方法

技术领域

本发明涉及穿墙雷达技术，特别涉及穿墙雷达的多目标检测跟踪技术。

背景技术

穿墙雷达是利用发射特定频段的电磁波穿透建筑物，接收建筑物后人体等目标的散射回波数据实现对隐蔽目标成像检测跟踪的特种装备，近年来已经在反恐、巷战、灾难救援等军事及民用方面得到越来越广泛的应用。由于建筑墙体会改变雷达所发射电磁波的传播路径及速度，其引起的回波时延误差会造成成像后目标图像散焦、位置偏移和多径幻象等现象。上述目标图像在后续的检测跟踪算法中会进一步影响目标提取以及目标数据关联。此外，墙体后运动目标的高机动性会导致跟踪性能下降。因此，实际应用中，从目标图像中准确提取目标位置信息、提高目标数据关联正确率、跟踪算法适配高机动目标是穿墙雷达成像后跟踪中的关键问题。

对于穿墙雷达隐蔽目标检测跟踪的研究，国内外研究机构已提出诸多解决方法。例如，通过最小均方误差和低复杂度非线性两步跟踪滤波器来得到更加平滑的航迹；利用alpha-beta滤波及卡尔曼滤波来对直线运动目标进行跟踪。但是这些方法目标运动模型都较为简单，没有考虑目标的高机动性问题。另外，检测定位算法获取目标位置一般采用椭圆交叉定位方法，这种方法下目标位置精度很大程度上依赖于天线的数量及其摆放位置。文献“A new measurement method for through-the-wall detection and tracking ofmoving targets.Measurement,vol.46,pp.1834-1848,2013.”中用到了成像后跟踪的思想，但是成像算法运算复杂度较高，跟踪算法较为简单，不适合实际应用。从公开发表的文献资料来看，目前还没有快速实时地实现成像后对墙体后多个高机动目标进行跟踪的方法。因此，研究一种可实时的多目标图像域跟踪算法在穿墙雷达目标检测跟踪中具有重要的实际应用价值。

发明内容

本发明所要解决的技术问题是，提供一种切实可行的，实时性好的，适用于穿墙雷达的多目标图像域跟踪方法。

本发明所采用的技术方案是，本发明基于多发多八收阵列配置的穿墙雷达，使用一种穿墙雷达多目标图像域稳健跟踪方法，包括以下步骤：

步骤1：快速成像步骤；

利用M个发射天线和N个接收天线配置的阵列对一层墙体后运动目标进行探测，天线阵列平行于墙体放置，与前墙表面的距离为d₃，根据探测场景的墙体厚度d_w以及介电常数ε，对需要探测的区域运用后向投影方法获得补偿墙体传播延迟下的图像I_BP(x)；

步骤2：图像增强步骤

根据像素点处各通道的相位差异计算相位相干因子(phase coherence factor,PCF)，对步骤1得到的图像I_BP(x)进行PCF加权，经过PCF加权后的序列图像表示为{I_k(x)},k＝1,2,...,N_f，N_f为序列图像总数；

步骤3：目标模板创建步骤

步骤3-1：利用膨胀操作确定第k帧图像I_k(x)的局部极大值，仅保留幅值较大的和高于噪声门限的极大值点；

步骤3-2：将某个极大值点占据的像素点作为初始模板，向四个方向：x轴正、负方向，y轴正、负方向各自求累加和；

步骤3-3：判断四个方向的累加和是否超过设定的门限T_s，如果高于门限则模板向对应方向的尺寸加一，如果低于门限则尺寸暂停放大；

步骤3-4：循环步骤3-2和3-3，直到四个方向：x轴正、负方向，y轴正、负方向各自的的累加和均低于门限，此时得到的矩形的内切圆为当前极大值对应的目标模板，表示目标模板包含的像素点集合，n表示集合大小；

步骤4：将第k帧图像I_k(x)的幅度范围分成m个区间，计算目标模板的幅度直方图

步骤5：航迹更新步骤；

步骤5-1：将第k-1帧的位置估计y(k-1)作为当前帧的第一个候选目标模板，设该候选模板中心为计算候选目标模板直方图

步骤5-2：计算Bhattacharyya系数

步骤5-3：计算每个像素点的权值，对于第i个像素点，其权值为w_i；

步骤5-4：更新位置估计

步骤5-5：判断是否继续迭代；判断相邻两次迭代的位置满足或迭代次数达到预设的最大值，如果不满足判断条件，将作为新的候选模板，重复步骤5-1～步骤5-4，直到满足判断条件，此时得到的候选模板中心位置作为当前帧图像的Mean-shift位置估计其对应的Bhattacharyya系数为

步骤6：航迹质量评价与平滑；

步骤6-1：构造状态向量

步骤6-2：计算一步预测状态向量

步骤6-3：假设目标帧间最大运动速度为V_max，以为圆心，R＝V_maxT为半径设定圆形波门，其中T为每帧的时间长度；设定一个门限ρ_T判断Mean-shift位置估计是否可信；

如果且Mean-shift位置估计位于圆形波门内，则表明Mean-shift位置估计是可信的，将该结果作为卡尔曼滤波的量测输入；通过卡尔曼滤波得出第k帧目标位置估计此时，当前帧目标位置为

如果或Mean-shift位置估计位于圆形波门外，表明Mean-shift结果是不可信的，此时用一步预测状态估计作为当前帧目标位置进行输出，即

本发明的创新点：提供了一种基于图像域的多目标跟踪算法，可在多目标的场景中实现多个目标的稳健跟踪。

本发明提供适用于穿墙雷达成像后的多目标的图像域跟踪算法。在一定允许误差内快速成像，得到较高质量成像结果；基于PCF加权处理的图像增强算法能有效地抑制栅瓣的影响，增强目标；基于双轴投影的目标模板创建算法能有效地提取目标的特征信息，为跟踪做准备；基于均值漂移和卡尔曼滤波的图像域跟踪算法能有效地应对目标图像扩展和形变的问题，也能轻松应对目标遮挡问题，跟踪精度高。因此，本发明具有实时性好，跟踪精度高的优势，可以直接应用到穿墙雷达装备中。

附图说明

图1为本实施的处理流程图；

图2为实验场景；

图3为算法跟踪结果图。

具体实施方式

步骤1：快速成像步骤

探测场景如图2所示，目标位于厚度与相对介电常数分别为d_w和ε的墙体后方，其坐标为T(x_tar,y_tar)采用时分MIMO(multiple-input-multiple-output)天线阵列，M个发射天线、N个接收天线的天线阵列平行于墙体放置，与前墙表面的距离为d₃；对于第i个发射天线和第j个接收天线，A_i、B_i、C_j、D_j为电磁波传播过程中的折射点，则为电磁波传播过程中的路径；发射信号s(t)为步进频连续波信号，对于第i个发射天线和第j个接收天线，得到的回波信号为：

y_ij(t)＝σ_Ts(t-τ_ij)+ψ_ij(t)

其中σ_T表示目标的反射系数，τ_ij为目标的回波延迟，ψ_ij(t)表示噪声和杂波；

探测区域可以划分为X×Y个像素点，根据快速后向投影(BP)成像算法，对于像素点x_h(x_h,y_h)，

其中τ_ijh表示像素点x_h(x_h,y_h)对应的回波时延。

对成像区域内所有像素点进行上述操作即可得到图像I_BP(x)。

步骤2：图像增强步骤

根据像素点处各通道的相位差异计算相位相干因子(phase coherence factor,PCF)，对于图像I_BP(X,Y)中的第q个像素点，PCF的值可表示为：

其中，∠y_k(τ_q,k)表示回波信号y_k(τ_q,k)的相位，τ_q,k表示第k个通道、第q个像素点的回波延迟，p为敏感因子且大于等于1，std(·)表示标准差运算；

经过PCF加权后的序列图像表示为{I_k(x)},k＝1,2,...,N_f，N_f为序列图像总数；

步骤3：目标模板创建步骤

步骤3-1：利用膨胀操作确定第k帧图像I_k(x)的局部极大值，仅保留幅值较大的和高于噪声门限的极大值点；

步骤3-2：将某个极大值点占据的像素点作为初始模板，向四个方向(x轴正、负方向，y轴正、负方向)求累加和；

步骤3-3：判断四个方向的累加和是否超过设定的门限T_s，如果高于门限则模板向对应方向的尺寸加一，如果低于门限则尺寸暂停放大；

步骤3-4：循环步骤3-2和3-3，直到四个方向(x轴正、负方向，y轴正、负方向)的累加和均低于门限，此时得到的矩形的内切圆为当前极大值对应的目标模板，表示目标模板包含的像素点集合，n表示集合大小。

步骤4：将第k帧图像I_k(x)的幅度范围分成m个区间，目标模板的特征可以用幅度直方图表示，可表示为

其中为像素点对应的幅度区间编号，δ为Kronecker函数；C为保证成立的归一化常数：

为像素点的权重，k(·)为核函数，本发明采用Epanechnikov核函数如下所示：

步骤5：航迹更新步骤

步骤5-1：将第k-1帧的位置估计y(k-1)作为当前帧的第一个候选目标模板，假设该候选模板中心为计算候选目标模板每个幅度区间的直方图对于每一个区间u：

其中n_h表示候选目标模板包含的像素点数量，h为候选模板相对目标模板的缩放因子，b(x_i)为像素点x_i对应的幅度区间编号，δ为Kronecker函数，C_h是保证成立归一化常量，可以得到：

步骤5-2：计算Bhattacharyya系数：

其中，为目标的幅度直方图；

步骤5-3：计算权值，对于第i个像素点，其权值为：

其中，δ为Kronecker函数，抑制背景的系数W_u＝u-1；

步骤5-4：更新位置估计：

其中，w_i为第i个像素点的权值，g(x)＝-k'(x)并设k(x)的一阶导数k′(x)在区间x∈[0,∞]除有限点外都存在；

步骤5-5：判断是否继续迭代；判断相邻两次迭代的位置满足或迭代次数达到预设的最大值，ε表示根据实际情况设定的阈值；满足判断条件，将作为新的候选模板，重复步骤5-1～步骤5-4，直到满足判断条件，此时得到的候选模板中心位置作为当前帧图像的Mean-shift位置估计其对应的Bhattacharyya系数为

步骤6：航迹质量评价与平滑；

步骤6-1：构造状态向量：

其中，速度v由前两帧位置计算得到；

步骤6-2：计算一步预测状态向量：

其中，和分别代表从第k-1帧一步预测第k帧的状态向量和位置估计，F为状态转移矩阵；

步骤6-3：假设目标帧间最大运动速度为V_max，以为圆心，R＝V_maxT为半径设定圆形波门，其中T为每帧的时间长度；根据实际情况设定一个门限ρ_T判断Mean-shift位置估计是否可信；

如果且Mean-shift位置估计位于圆形波门内，则表明Mean-shift位置估计是可信的，将该结果作为卡尔曼滤波的量测输入；卡尔曼滤波的步骤具体如下：

P(k|k-1)＝FP(k-1)F^T+Q_w

K(k)＝P(k|k-1)C^T[CP(k|k-1)C^T+Q_v]^-1

P(k)＝[I-K(k)C]P(k|k-1)

其中，P(k|k-1)为一步预测协方差矩阵，Q_w和Q_v分别为系统噪声协方差矩阵和量测噪声协方差矩阵，K(k)为第k帧卡尔曼增益矩阵，C为观测矩阵，P(k)为第k帧状态估计的误差自相关矩阵，为第k帧目标位置估计，具体可表示为：

此时，当前帧目标位置为

如果或Mean-shift位置估计位于圆形波门外，表明Mean-shift结果是不可信的，此时用一步预测状态估计作为当前帧目标位置进行输出，即

本发明的效果通过以下实验进一步说明：

实验场景：

探测场景如图2所示，利用超宽带穿墙雷达系统探测墙后两个运动人体目标。两个发射天线分别放在(-1.2m,0m)和(1.2m,0m)，八个接收天线放置在(-1.05m,0m)到(1.05m,0m)，间距为0.3m。发射天线发射1GHz～2GHz的步进频率连续波信号，步进间隔为2MHz。探测区域大小为8m×8m。墙体厚度为0.25m，天线与墙体平行，离墙5m。两个人体目标在墙后交叉走动，目标A和B分别从(3.0m,12.0m)走到(-1.0m,6.0m)和从(-2.0m,11.0m)走到(2.0m,7.0m)，两者速度几乎一致。

对于上述所示的穿墙雷达系统探测得到的回波数据，采用本发明的处理方案，得到的跟踪结果如图3所示，其中虚线表示两个目标的真实航迹，椭圆表示当前帧的目标模板轮廓。由图3(a)-(c)可以看出，目标A和B均出现扩展和散焦现象，而且随着时间推移，目标形状和轮廓都发生改变，本发明的算法依然能够用一个椭圆包裹住目标的大部分能量，实现目标的整体跟踪；由图3(b)可以看出，当目标发生遮挡时，本发明的跟踪算法能够发挥其预测功能，有效应对多目标环境下出现的目标遮挡问题；如图3(c)所示，本文的跟踪算法能从图像域中对目标进行稳健的跟踪，航迹整体比较平滑，而且非常接近真实航迹。

通过本发明的具体实施可以看出，即使目标图像扩展成多个像素点、发生散焦现象，本发明能在图像域中对目标进行稳健的跟踪。

Claims (6)

1.一种穿墙雷达多目标图像域稳健跟踪方法，包括以下步骤：

步骤1：快速成像步骤；

利用M个发射天线和N个接收天线配置的阵列对一层墙体后运动目标进行探测，天线阵列平行于墙体放置，与前墙表面的距离为d₃，根据探测场景的墙体厚度d_w以及介电常数ε，对需要探测的区域运用后向投影方法获得补偿墙体传播延迟下的图像I_BP(x)；

步骤2：图像增强步骤；

根据像素点处各通道的相位差异计算相位相干因子，对步骤1得到的图像I_BP(x)进行PCF加权，经过PCF加权后的序列图像表示为{I_k(x)},k＝1,2,...,N_f，N_f为序列图像总数；

步骤3：目标模板创建；

步骤3-1：利用膨胀操作确定第k帧图像I_k(x)的局部极大值，仅保留幅值较大的和高于噪声门限的极大值点；

步骤3-2：将某个极大值点占据的像素点作为初始模板，向四个方向：x轴正、负方向，y轴正、负方向各自求累加和；

步骤3-3：判断四个方向的累加和是否超过设定的门限T_s，如果高于门限则模板向对应方向的尺寸加1，如果低于门限则尺寸暂停放大；

步骤3-4：循环步骤3-2和3-3，直到四个方向：x轴正、负方向，y轴正、负方向各自的的累加和均低于门限，此时得到的矩形的内切圆为当前极大值对应的目标模板，表示目标模板包含的像素点集合，n表示集合大小；

步骤4：将第k帧图像I_k(x)的幅度范围分成m个区间，计算目标模板的幅度直方图

步骤5：航迹更新步骤；

步骤5-1：将第k-1帧的位置估计y(k-1)作为当前帧的第一个候选目标模板，设该候选模板中心为计算候选目标模板直方图

步骤5-2：计算Bhattacharyya系数

步骤5-3：计算每个像素点的权值，对于第i个像素点，其权值为w_i；

步骤5-4：更新位置估计

步骤5-5：判断是否继续迭代；判断相邻两次迭代的位置满足或迭代次数达到预设的最大值，如果不满足判断条件，将作为新的候选模板，重复步骤5-1～步骤5-4，直到满足判断条件，此时得到的候选模板中心位置作为当前帧图像的Mean-shift位置估计其对应的Bhattacharyya系数为

步骤6：航迹质量评价与平滑；

步骤6-1：构造状态向量

步骤6-2：计算一步预测状态向量

步骤6-3：假设目标帧间最大运动速度为V_max，以为圆心，R＝V_maxT为半径设定圆形波门，其中T为每帧的时间长度；设定一个门限ρ_T判断Mean-shift位置估计是否可信；

如果且Mean-shift位置估计位于圆形波门内，则表明Mean-shift位置估计是可信的，将该结果作为卡尔曼滤波的量测输入；通过卡尔曼滤波得出第k帧目标位置估计此时，当前帧目标位置为

如果或Mean-shift位置估计位于圆形波门外，表明Mean-shift结果是不可信的，此时用一步预测状态估计作为当前帧目标位置进行输出，即

2.如权利要求1所述的一种穿墙雷达多目标图像域稳健跟踪方法，其特征在于所述步骤1的具体方法为：

设目标位于厚度与相对介电常数分别为d_w和ε的墙体后方，其坐标为T(x_tar,y_tar)采用时分MIMO(multiple-input-multiple-output)天线阵列，M个发射天线、N个接收天线的天线阵列平行于墙体放置，与前墙表面的距离为d₃；对于第i个发射天线和第j个接收天线，A_i、B_i、C_j、D_j为电磁波传播过程中的折射点，则为电磁波传播过程中的路径；发射信号s(t)为步进频连续波信号，对于第i个发射天线和第j个接收天线，得到的回波信号为：

y_ij(t)＝σ_Ts(t-τ_ij)+ψ_ij(t)

其中σ_T表示目标的反射系数，τ_ij为目标的回波延迟，ψ_ij(t)表示噪声和杂波；

探测区域可以划分为X×Y个像素点，根据快速后向投影(BP)成像算法，对于像素点x_h(x_h,y_h)，

其中τ_ijh表示像素点x_h(x_h,y_h)对应的回波时延；

对成像区域内所有像素点进行上述操作即可得到图像I_BP(x)。

3.如权利要求2所述的一种穿墙雷达多目标图像域稳健跟踪方法，其特征在于所述步骤2的具体方法为：根据像素点处各通道的相位差异计算相位相干因子(phase coherencefactor,PCF)，对于图像I_BP(X,Y)中的第q个像素点，PCF的值可表示为：

∠Y_q＝{∠y_k(τ_q,k),k＝1,2,...,MN},

其中，∠y_k(τ_q,k)表示回波信号y_k(τ_q,k)的相位，τ_q,k表示第k个通道、第q个像素点的回波延迟，p为敏感因子且大于等于1，std(·)表示标准差运算；

经过PCF加权后的序列图像表示为{I_k(x)},k＝1,2,...,N_f，N_f为序列图像总数。

4.如权利要求2所述的一种穿墙雷达多目标图像域稳健跟踪方法，其特征在于所述步骤4的具体方法为：

步骤4：将第k帧图像I_k(x)的幅度范围分成m个区间，目标模板的特征可以用幅度直方图表示，可表示为

其中为像素点对应的幅度区间编号，δ为Kronecker函数；C为保证成立的归一化常数：

为像素点的权重，k(·)为核函数，本发明采用Epanechnikov核函数如下所示：

5.如权利要求2所述的一种穿墙雷达多目标图像域稳健跟踪方法，其特征在于所述步骤5的具体方法为：

步骤5-1：将第k-1帧的位置估计y(k-1)作为当前帧的第一个候选目标模板，假设该候选模板中心为计算候选目标模板每个幅度区间的直方图对于每一个区间u：

其中n_h表示候选目标模板包含的像素点数量，h为候选模板相对目标模板的缩放因子，b(x_i)为像素点x_i对应的幅度区间编号，δ为Kronecker函数，C_h是保证成立归一化常量，可以得到：

步骤5-2：计算Bhattacharyya系数：

其中，为目标的幅度直方图；

步骤5-3：计算权值，对于第i个像素点，其权值为：

其中，δ为Kronecker函数，抑制背景的系数W_u＝u-1；

步骤5-4：更新位置估计：

其中，w_i为第i个像素点的权值，g(x)＝-k'(x)并设k(x)的一阶导数k′(x)在区间x∈[0,∞]除有限点外都存在；

步骤5-5：判断是否继续迭代；判断相邻两次迭代的位置满足或迭代次数达到预设的最大值，ε表示根据实际情况设定的阈值；满足判断条件，将作为新的候选模板，重复步骤5-1～步骤5-4，直到满足判断条件，此时得到的候选模板中心位置作为当前帧图像的Mean-shift位置估计其对应的Bhattacharyya系数为

6.如权利要求2所述的一种穿墙雷达多目标图像域稳健跟踪方法，其特征在于所述步骤6的具体方法为：

步骤6-1：构造状态向量：

其中，速度v由前两帧位置计算得到；

步骤6-2：计算一步预测状态向量：

其中，和分别代表从第k-1帧一步预测第k帧的状态向量和位置估计，F为状态转移矩阵；

步骤6-3：假设目标帧间最大运动速度为V_max，以为圆心，R＝V_maxT为半径设定圆形波门，其中T为每帧的时间长度；根据实际情况设定一个门限ρ_T判断Mean-shift位置估计是否可信；

如果且Mean-shift位置估计位于圆形波门内，则表明Mean-shift位置估计是可信的，将该结果作为卡尔曼滤波的量测输入；卡尔曼滤波的步骤具体如下：

P(k|k-1)＝FP(k-1)F^T+Q_w

K(k)＝P(k|k-1)C^T[CP(k|k-1)C^T+Q_v]^-1

P(k)＝[I-K(k)C]P(k|k-1)

其中，P(k|k-1)为一步预测协方差矩阵，Q_w和Q_v分别为系统噪声协方差矩阵和量测噪声协方差矩阵，K(k)为第k帧卡尔曼增益矩阵，C为观测矩阵，P(k)为第k帧状态估计的误差自相关矩阵，为第k帧目标位置估计，具体可表示为：

此时，当前帧目标位置为

如果或Mean-shift位置估计位于圆形波门外，表明Mean-shift结果是不可信的，此时用一步预测状态估计作为当前帧目标位置进行输出，即