基于交互式多模型的雷达单目标跟踪算法

杨贵福，刘鲁涛

哈尔滨工程大学 信息与通信工程学院，黑龙江 哈尔滨 150001

随着科技的发展及人们生活水平的提高，民用小型无人机由于获取便利、使用便捷而呈现爆发式增长。该类无人机在航拍、娱乐等领域已被广泛应用，在带来产业升级的同时，也带来了在机场等禁飞区“黑飞”的各种安全威胁[1]。雷达系统是对该类无人机的管控与反制的常用手段之一，目标跟踪作为其中非常重要的一环，成为了一个备受关注的研究课题。

雷达目标跟踪是对目标的位置以及运动状态等参数进行估计，形成目标的运动轨迹。对目标进行跟踪过程中，噪声、目标的机动情况等会影响跟踪的结果。在雷达目标跟踪研究的初期，通常采用卡尔曼滤波算法对目标的运动状态进行估计和预测，建立起目标的运动轨迹。对于民用小型无人机这类的飞行器，由于受到人为操纵控制，并且自身运动性能不断提高，其机动性越来越强。民用小型无人机运动状态的复杂性和不确定性逐渐提高，这对雷达的跟踪性能有了更高的要求。卡尔曼滤波算法在目标运动状态比较单一时有很好的跟踪效果，当目标实施机动时(突然转弯或加、减速等)，卡尔曼滤波算法对目标的跟踪性能会下降[2]。为了提高对目标运动状态描述的精确度，Blom 等[3]提出交互式多模型(interacting multiple model，IMM)结构算法，该方法使用2 个或更多的运动模型来描述目标可能的运动状态，最后通过有效地加权融合对目标运动状态估计，很好地克服了单模型滤波算法估计误差较大的问题。交互式多模型算法是一种具有马尔可夫转移概率的结构自适应算法[4]，目前在机动目标跟踪领域得到了广泛的应用。本文主要在IMM算法的基础上，对其在单目标雷达数据处理系统中的应用进行了分析和仿真。

1 交互式多模型算法

1.1 目标运动模型

在构造目标的状态方程之前，需要建立目标运动模型。目标运动模型是对目标运动特性的描述，是目标跟踪的基本要素之一[5]。目标跟踪性能的优劣与目标运动模型建立的符合度有直接关系。目标运动模型的构造需要尽量还原目标实际的运动情况，同时应尽量简单以便进行数据处理。

假设在雷达探测范围内出现某种民用小型无人机在空中飞行，雷达系统将该无人机作为运动目标进行跟踪探测。当目标在空中做匀速直线运动时，可选用匀速(constant velocity，CV)运动模型描述目标的这种运动状态，这是目标跟踪中最基本的一种运动模型[6]。

将运动目标在t时刻的位置分量、速度分量、加速度分量分别表示为x(t)、。因为目标做匀速直线运动，则=0。但运动目标会受到环境扰动带来的随机干扰，加速度具有随机特性的扰动输入，可以用连续时间白噪声w(t)建模，并假设其服从均值为零的高斯分布，即

此时连续时间的状态方程表示为[7]

式中：W(k)为过程噪声，其为均值为零、协方差矩阵为Q的 高斯白噪声序列；F(k)为系统状态转移矩阵，其中T为采样间隔。

当目标在空中做匀加速直线运动时，可选用匀加速(constant acceleration，CA)运动模型描述目标的这种运动状态。此时加速度的一阶导数考虑随机干扰的情况，是具有随机特性的扰动输入，可以用白噪声过程建模表示为

此时连续时间的状态方程可表示为

状态方程是对运动目标的观测描述，而雷达系统对运动目标的观测过程通过量测方程表示。线性系统的量测方程表示为

式中：Z(k)为量测向量；H(k)为量测矩阵；V(k)为量测噪声序列，是对雷达在量测过程中存在误差的模拟，一般设置为均值为零、协方差矩阵为R的高斯白噪声。

1.2 算法原理

IMM 算法采用多个状态空间模型同时对目标的可能行为模式进行描述，不同的状态空间模型对应不同的卡尔曼滤波器[8]。IMM 算法的基本思想是：首先，在每一时刻通过对前一时刻所有滤波器的状态估计值混合，得出与相应模型所匹配的滤波器初始条件；然后，对每个模型按照滤波算法的步骤并行处理；最后，用模型匹配似然函数更新的模型概率作为加权值，对每个滤波器的状态估计值进行组合，得到的加权合并结果就是对目标状态的估计[9]。本文设计的IMM 算法中模型集由CV 与CA 组成，模型个数为N。IMM算法的一般步骤及计算公式如下。

1)输入交互。

对于算法模型集中的任意模型j(j=1,2,···,N)，由模型i(i=1,2,···,N)到模型j的混合概率为

式中：r为状态空间模型个数；pij为模型i到模型j的转移概率；μi(k)为当前时刻每个滤波器的模型概率。

模型j的混合状态估计：

模型j的混合协方差估计：

2)对于模型j，进行卡尔曼滤波。

状态预测：

预测误差协方差：

残差：

残差协方差：

卡尔曼增益：

状态更新：

预测误差协方差更新：

3)模型概率更新。

采用似然函数更新模型概率μj(k)，模型j的似然函数为

模型j的概率更新为

式中c=

4)输出融合

总的状态估计：

总的协方差估计：

2 坐标转换误差分析

在许多雷达跟踪系统中，目标的量测信息通常在极坐标系下获得，但在笛卡尔直角坐标系中对目标运动状态进行建模与跟踪滤波处理，需要利用坐标变换的方法将极坐标下的量测值转换成直角坐标形式。由于雷达在极坐标系下量测的数据中含有噪声，转换至直角坐标系下的数据有误差存在，这成为雷达跟踪系统必须考虑的问题[10]。

设二坐标雷达在极坐标系下对目标进行量测，得到的测量距离为rm，测量方位角为 θm，且有

式中：r为目标真实的距离；θ为目标真实的方位角；为相应的测量误差，且其是不相关的零均值高斯分布，方差分别为

通过极-直坐标转换，可得直角坐标系下目标位置的测量值为

由于式(1)与式(2)中包含了非线性转换，如果直接将转换的数据作为跟踪算法的目标观测数据进行滤波，则得到的滤波估计结果为有偏估计[11]。因此需要对这种转换的误差进行分析，然后对误差进行补偿，得到目标观测状态的无偏转换数据。

设xe、ye是直角坐标系的误差，则有

按三角函数的和角公式展开得

同理可得：

因此，转换后的直角坐标误差(xe,ye)在真值(r,θ)条件下的均值为

以上是在真值(r,θ)条件下求取的均值，在工程实践中无法直接应用。为了获取能在实践中使用的结果，在量测数据(rm,θm)的条件下求取xe和ye的均值。

在 μr中，用代换r，用代换θ，根据三角函数和差化积公式，可求得坐标转换的转换误差为

因此，雷达量测数据坐标转换的无偏转换公式为

即直角坐标(xz,yz)是极坐标(rm,θm)的无偏转换，μμ为坐标转换的转换误差。

转换误差的协方差矩阵为[12]

式中：

3 单目标雷达数据处理系统仿真

假定使用二坐标雷达对低空飞行的民用小型无人机进行观测。目标的起始点为(500 m，2 000 m)，在t=0～50s 沿y轴以运动速度为-20 m/s 做匀速直线运动；在t=50～60 s 向x轴方向做90°的慢转弯运动，加速度ax=ay=2 m/s2，完成慢转弯后加速度将降为零；在t=60～66 s 沿x轴做匀速运动；从t=66s 开始做 90°的快转弯运动，加速度为5 m/s2；在t=70s 结束转弯，加速度降为零；在t=70～110 s沿y轴做匀速直线运动。雷达扫描周期T=2 s，目标距离与方位角测量噪声的标准差分别为σr=10 m，σθ=1°。转换误差的协 方差矩阵为

IMM 滤波算法采用3 个模型对目标的运动状态进行描述：第一个模型是非机动的CV 模型，过程噪声的方差为零；第二、三个模型均为CA 机动模型，但过程噪声的方差不同，第二个模型的过程噪声方差为Q=0.001I2×2，第三个模型的过程噪声方差为Q=0.011 4I2×2。控制模型转换的马尔可夫链的转移概率矩阵为

设定各模型的模型概率初始值分别为μ1=0.5、μ2=0.25、μ3=0.25。

定义滤波误差的均值为

定义滤波误差的标准差为

式中：M为蒙特卡洛模拟次数；k=1,2,···,N，N为采样次数。

进行1 000 次蒙特卡洛仿真实验，得到运动目标在极坐标下的真实轨迹、观测轨迹及一次滤波估计结果曲线，如图1—图3 所示。其中，图3 滤波估计结果是在无偏坐标转换的条件下得出的。

图1 目标真实轨迹

图2 目标观测轨迹

图3 目标估计轨迹(一次滤波曲线)

从图3 中可以看出，无论目标在非机动还是机动状态，滤波轨迹相对于观测轨迹更加接近于目标真实轨迹，具有较好的滤波效果，说明了IMM算法能够实现对雷达单目标进行有效的跟踪。

IMM 算法的滤波过程包括了对目标状态原始数据初始化、目标状态预测及目标状态更新。在极坐标下对目标距离、方位角滤波过程的误差分析曲线如图4 所示。

图4 极坐标下滤波误差分析曲线

图4 曲线表示了滤波过程中的误差。目标观测状态是由雷达系统测量得到的，通过初始误差曲线可以看出，在各个时刻目标量测误差都较大。预测误差曲线是由IMM 算法对目标状态的预测所形成的，当目标处于机动状态时，IMM 算法的预测结果会出现较大的误差。在滤波过程中IMM 算法利用观测数据对预测状态进行修正。对比预测误差曲线与更新误差曲线可以看出：在目标机动时滤波得到的目标状态估计结果比预测结果的误差要小；当目标状态由机动转换至非机动时，目标状态估计误差会趋近于零。

在极坐标下，方位角滤波误差标准差曲线和距离误差标准差曲线如图5 和图6 所示。

图5 极坐标下方位角误差标准差

图6 极坐标下距离误差标准差

图5 是在极坐标下对方位角误差标准差的分析曲线。分别将方位角在测量、有偏估计和无偏估计的条件下进行滤波误差标准差分析。方位角测量值的误差标准差约为1°，与仿真条件中设置的方位角量测噪声标准差相符。对量测值进行极-直坐标转换，通过交互式多模型算法对转换结果进行滤波处理，得到的方位角估计误差标准差约为0.7°。由于量测值直接进行坐标转换存在转换误差，因此对转换误差补偿之后再用滤波算法处理，得到的方位角无偏估计误差标准差约为0.55°，这比测量值与有偏估计值更加接近于目标运动状态的真实值。

同样，由图6 可以看出，在极坐标下距离测量误差标准差约为10 m，量测值通过坐标转换及滤波处理后，距离的有偏估计误差标准差约为6.5 m，而对坐标转换误差补偿之后滤波得到的距离误差标准差稍小，说明无偏估计结果更精确一些。目标在第二次转弯时要比第一次转弯更急，因此在误差标准差曲线上产生了更大的波动；在滤波进入稳定状态之后，误差标准差又减小，这体现出IMM 算法受目标机动的影响很小。

为了比较交互式多模型算法与卡尔曼滤波算法对雷达目标的跟踪性能，分别将雷达目标的量测值使用这2 种算法进行滤波处理，得到目标位置估计结果在极坐标下的估计误差如图7 所示。

图7 极坐标下估计误差曲线

由图7 的估计误差曲线可以看出：雷达目标在初始阶段的匀速运动状态下，交互式多模型算法与卡尔曼滤波算法的估计误差都比较小；但目标进行机动运动之后，卡尔曼滤波算法对目标机动的适应性较低，因此估计误差逐渐增大，使其失去对目标跟踪的能力，而交互式多模型算法对机动目标有很好的适应性，得到估计结果的误差较小，能够实现对目标进行有效跟踪。

为了分析交互式多模型算法在实际应用时的可行性，在外场环境实验中，将低空飞行的大疆公司小型无人机作为雷达目标，利用线性调频连续波雷达对该目标进行观测，无人机在一个运动过程通过实测数据形成的目标运动轨迹如图8 所示。

图8 实测数据形成的目标运动轨迹

无人机在这个运动过程中通过飞行记录所显示的目标真实运动轨迹如图9 所示。

图9 飞行记录显示的无人机运动轨迹

通过卡尔曼滤波算法与交互式多模型算法对雷达目标实测数据进行跟踪滤波处理，得到的滤波估计结果所形成的目标运动轨迹如图10 所示。由图10 可以看出，卡尔曼滤波算法的滤波估计轨迹出现了很大的偏差，交互式多模型算法的滤波估计轨迹与目标飞行记录显示的运动轨迹比较符合。这是由于卡尔曼滤波算法在滤波过程中只通过单一模型对目标的运动状态进行估计，但是无人机的真实运动状态会出现机动情况，只通过单一模型不能概括目标所有的运动状态，使得滤波结果会出现较大的误差，还可能出现滤波发散的现象。交互式多模型算法能够采用多个目标运动模型对目标的运动状态进行滤波估计，因此在目标出现机动运动状态时该滤波算法能实现对目标运动状态的估计，并且滤波误差较小。

图10 目标的滤波估计轨迹

4 结论

本文对交互式多模型算法在雷达单目标跟踪中的应用情况进行了描述。针对民用小型无人机的跟踪问题，对目标的运动模型进行介绍，分析了量测数据坐标转换的转换误差，并给出无偏坐标转换公式。通过蒙特卡洛仿真实验对交互式多模型算法在机动目标跟踪时的跟踪性能进行分析，仿真结果表明，滤波估计轨迹相对于观测轨迹更加接近于目标真实轨迹，并且不会有大的偏差，实现了对机动目标的跟踪能力，证明了算法估计的有效性。

通过外场环境实验所得到的实测数据对卡尔曼滤波算法与交互式多模型算法的目标跟踪性能进行比较，实测数据分析结果表明交互式多模型算法能够实现对雷达单目标进行有效的跟踪，说明了该算法在雷达目标跟踪过程的实用性，体现出一定的工程应用价值。