CN110376582B - 自适应gm-phd的机动目标跟踪方法 - Google Patents

Abstract

本发明公开了一种自适应GM‑PHD的机动目标跟踪算法，主要解决现有算法对多个机动目标突发性机动跟踪精度不高、计算量和存储量大的问题。本发明实现的方法是：先读取存活目标和新生目标的状态参数和量测值。再对多个机动目标的状态值、存活概率和状态误差协方差矩阵进行预测，并通过反应机动大小的修正因子对协方差矩阵进行修正。然后利用机动误差变化率来自适应调整量测门限，并对量测值进行裁剪。然后利用状态预测值和裁剪后的量测值对多个机动目标的状态进行更新。然后通过阈值对更新后的高斯分量进行裁剪与合并。最后提取合并后机动目标的状态参数。

Description

自适应GM-PHD的机动目标跟踪方法

技术领域

本发明属于雷达技术领域，更进一步涉及目标跟踪技术领域中的一种自适应高斯混合概率假设密度滤波GM-PHD(Gaussian mixture probability hypothesis density)的机动目标跟踪方法。本发明可用于雷达对未知数目的具有突发性机动的多个地面目标进行精确跟踪。

背景技术

多目标跟踪技术广泛应用于军事和民用领域，如战场监视(地面车辆与军事设施)，海洋监视、交通导航以及智能驾驶等。传统的机动多目标跟踪方法主要包括数据关联和状态估计两个核心步骤。在稀疏杂波环境中，对非机动的目标跟踪时，此类算法结构简单，易于工程实现，跟踪效果较好。但在密集回波环境下，数据关联过程存在“组合爆炸”的问题，计算复杂度非常高而且存储量异常大。如果目标距离比较近或者轨迹一段时间内发生小角度交叉时，就会引起航迹“聚集”的问题，致使对目标跟踪出现漏跟、误跟以及航迹合并的现象。所以在对实时性要求高的系统中难以应用。另外，对于状态的估计，传统的机动模型都是依靠先验信息建立的，如果目标的运动状态能准确匹配相应的模型，则能进行准确跟踪。但是目标的机动时刻和机动水平往往是未知的，会导致模型失配滤波误差变大。近年来出现的一种基于随机集理论RST(Random Set Theory)。的多目标跟踪算法避免了数据关联的问题，直接在贝叶斯框递推框架下进行预测和更新，在降低计算量的同时达到了估计精度高、容易实现的目的。

中国石油天然气股份有限公司在其申请的专利文献“复杂环境下多机动目标跟踪方法”(专利申请号：2015109647280，申请公布号：CN106910211A)中公开了一种利用雷达跟踪多机动目标的方法。该方法的步骤是：首先，获取多个运动目标先验信息，根据所述多个运动目标信息构造估计模型。其次，根据联合概率数据关联算法对多个运动目标进行点迹和航迹关联，并获得关联后的跟踪模型。最后，根据所述模型对多个运动目标进行状态估计滤波跟踪。该方法存在的不足之处有两点：其一，该方法是利用目标的位置、速度、运动模型等先验信息对多机动目标跟踪时，存在模型不匹配导致运动目标跟踪误差大的问题；其二，由于该方法使用的联合概率数据关联算法在量测回波密集的情况下，存在计算量和存储量大的问题，难以满足运动目标跟踪中实时性的要求。

Peng Dong，Zhongliang Jing，Minzhe Li，Han Pan在其发表的论文“Thevariable structure multiple model GM-PHD filter based on likely-model setalgorithm”(201619th International Conference on Information Fusion，FUSION，04August 2016)中公开了一种基于变结构交互式多模型GM-PHD的多目标跟踪方法。该方法的具体步骤如下，首先，将目标的状态和量测分别用随机集进行近似，并建立多个机动目标的运动模型集。其次，在贝叶斯递推框架下进行状态预测和状态更新。再次，进行高斯分量的修剪与合并。最后，进行目标状态提取、个数估计、目标运动模型集更新。该方法存在的不足之处是，如果目标运动状态从匀加速、转弯等机动状态变为匀速直线运动的非机动状态，该方法所使用的模型集只能从下一时刻开始调整不能在当前时刻及时进行切换，则会导致运动目标跟踪效果变差。而且该方法在建立模型后直接进行状态预测和状态更新，如果量测集中元素过多，会产生大量的杂波高斯分量，使得该方法在复杂环境下对运动目标的跟踪依然存在计算量大的问题，难以满足目标跟踪过程中实时性的要求。

发明内容

本发明的目的在于针对现有技术的不足，提出一种自适应斯混合概率假设密度滤波GM-PHD的机动目标跟踪方法，以解决现有模型滤波方法不能准确的同时对多个状态突变的目标或者运动形式多样化的目标进行跟踪的问题，并通过门限对高斯分量的产生进行了合理的限制，以解决多目标跟踪方法实时性不高的问题。

实现本发明的思路是：采用自适应的高斯混合概率假设密度率波方法实现对多个机动目标的跟踪。首先，通过对目标状态参数进行预测得到预测的状态参数值，然后，通过修正因子和机动变化率自适应改变滤波增益、机动变化率和量测裁剪门限，然后，通过量测门限对量测值进行合理的裁剪，然后，结合状态预测值预测值、裁剪后的量测值和裁剪后量测预测值更新机动目标的状态，然后，裁剪与合并更新后的高斯分量、剔除虚假高斯分量，最后，提取真实的目标状态。

本发明的具体步骤如下：

(1)读取机动目标的状态值：

(1a)读取当前时刻仍然存活的所有机动目标在上一时刻的状态参数和量测值；

(1b)读取当前时刻新生机动目标的状态参数和量测值；

(2)预测当前时刻每个机动目标的状态参数：

(2a)按照下式，预测每个存活机动目标在当前时刻的状态值：

m_i＝F_im′_C+U_ia_i

其中，m_i表示当前时刻第i个存活的机动目标的状态预测值，F_i表示由雷达扫描周期和机动目标的机动频率获得的当前时刻第i个存活的机动目标状态转移矩阵，m′_C表示当前时刻的上一时刻第C个存活的机动目标的状态值，U_i表示由雷达扫描周期和机动目标的机动频率获得的当前时刻第i个存活的机动目标传输控制矩阵，a_i表示当前时刻第i个存活的机动目标加速度，其中，i＝1,2,...,J,C＝1,2,...,M，J表示当前时刻存活的机动目标总数，M表示上一时刻存活的机动目标总数，且i和C的取值对应相等，J和M的取值相等；

(2b)将每个存活目标当前时刻的上一时刻的存活概率乘以0.99，得到当前时刻该机动目标的存活概率预测值；

(2c)将当前时刻每个新生机动目标状态参数中的状态值、存活概率，分别作为该新生机动目标的状态预测值和存活概率预测值；

(2d)将所有当前时刻的存活目标和所有新生目标的存活概率预测值之和，作为当前时刻机动目标的预测目标总数；

(2e)按照下式，计算当前时刻每个机动目标的过程噪声矩阵；

其中，Q_i表示当前时刻第i个存活的机动目标的过程噪声矩阵，*表示相乘操作，α_i表示当前时刻第i个存活的机动目标的机动频率，

表示当前时刻第i个存活的机动目标的机动频率的倒数的方差，q表示当前时刻测量环境导致的随机扰动误差；

(2f)使用修正因子的计算方法，计算当前时刻修正的每个存活的机动目标预测误差修正因子；

(2g)利用修正的状态误差协方差矩阵公式，计算当前时刻每个存活的机动目标修正的预测误差协方差矩阵；

(2h)利用状态误差协方差矩阵公式，设定公式中的修正因子的值为1，计算当前时刻每个存活的机动目标实际的预测误差协方差矩阵；

(3)更新机动目标的状态参数：

(3a)利用修正的量测距离误差公式，计算当前时刻每个机动目标的量测距离误差值；

(3b)用实际的状态误差协方差替换修正的状态误差协方差矩阵，计算当前时刻每个机动目标的实际的量测距离误差；

(3c)按照下式，计算修正的当前时刻每个机动目标量测裁剪门限：

其中，D′_n表示修正的当前时刻第n个机动目标的量测裁剪门限，g_n表示当前时刻第n个机动目标的实际的量测距离误差，g′_n表示当前时刻第n个机动目标的量测距离误差值，D_n表示当前时刻第n个机动目标的量测裁剪门限，其中，n＝1,2,...,N，N＝J+W，N表示当前时刻的机动目标总数，W表示新生目标总数；

(3d)按照下式，计算修正的当前时刻每个机动目标机动频率：

其中，α′_n表示修正的当前时刻第n个机动目标的机动频率，α_n表示当前时刻第n个机动目标的机动频率；

(3e)剔除量测距离误差值大于修正的量测裁剪门限的目标；

(3f)利用增益矩阵公式，计算裁剪后当前时刻每个机动目标的增益矩阵：

(3g)利用更新状态值公式，计算裁剪后当前时刻每个机动目标的更新状态值：

(3h)利用更新存在概率公式，计算当前时刻所有的机动目标的更新存活概率；

(3i)利用更新误差协方差矩阵公式，计算当前时刻所有的机动目标的更新误差协方差矩阵；

(4)对机动目标进行裁剪与合并：

(4a)剔除当前时刻所有机动目标中更新权值小于0.001的机动目标；

(4b)合并当前时刻所有机动目标中，欧氏距离小于合并门限的目标；

(5)提取当前时刻裁剪与合并之后机动目标的状态参数：

(5a)提取当前时刻裁剪与合并之后所有机动目标存活概率大于0.5的目标数目，作为当前时刻的真实目标数；

(5b)提取当前时刻裁剪与合并之后所有机动目标存活概率大于0.5的目标更新状态值，作为当前时刻真实存在的目标状态值；

(5c)提取当前时刻裁剪与合并之后所有机动目标存活概率大于0.5的目标更新存活概率，作为当前时刻真实存在的更新存活概率；

(5d)提取当前时刻裁剪与合并之后所有机动目标存活概率大于0.5的目标更新协方差矩阵，作为当前时刻真实存在的目标协方差矩阵值。

本发明与现有技术相比具有以下优点：

第一，由于本发明采用自适应的方法调整量测门限值和机动频率，将当前时刻的机动变化情况反馈到下一时刻的跟踪模型中，并在计算状态误差协方差矩阵时加入修正因子，通过修正的状态误差协方差矩阵调整当前时刻的增益矩阵，克服了现有技术模型不匹配导致运动目标跟踪误差大的问题，使得本发明满足对机动目标状态存在突变的多个机动目标跟踪的准确性要求。

第二，由于本发明采用了对量测数据设置量测门限值，通过提取量测值和状态预测值的距离不超过门限值的量测数据，有效的减少了参与状态更新的量测集中的元素个数，克服了现有技术量测集中元素过多，状态更新过程会产生大量的杂波高斯分量，计算量大的问题，使得本发明满足对多个机动目标跟踪实时性的要求。

附图说明

图1是本发明的流程图；

图2是本发明的仿真图。

具体实施方式

下面结合附图对本发明作进一步的描述。

参照图1，对本发明的具体步骤作进一步的描述。

步骤1，读取机动目标的状态值。

读取当前时刻仍然存活的所有机动目标在上一时刻的状态参数和量测值，状态参数包括位置、速度、加速度、状态误差协方差和存活概率。

读取当前时刻新生机动目标的状态参数和量测值，状态参数包括位置、速度、加速度、状态误差协方差和存活概率。

步骤2，预测当前时刻每个机动目标的状态参数。

按照下式，预测每个存活机动目标在当前时刻的状态值：

m_i＝F_im′_C+U_ia_i

其中，m_i表示当前时刻第i个存活的机动目标的状态预测值，F_i表示由雷达扫描周期和机动目标的机动频率获得的当前时刻第i个存活的机动目标状态转移矩阵，m′_C表示当前时刻的上一时刻第C个存活的机动目标的状态值，U_i表示由雷达扫描周期和机动目标的机动频率获得的当前时刻第i个存活的机动目标传输控制矩阵，a_i表示当前时刻第i个存活的机动目标加速度，其中，i＝1,2,...,J,C＝1,2,...,M，且i和C的取值对应相等，J表示当前时刻存活的机动目标总数，M表示上一时刻存活的机动目标总数，且J和M的取值相等；

所述的机动目标状态转移矩阵如下：

其中，F_i表示当前时刻第i个存活的机动目标状态转移矩阵，T表示雷达扫描周期，α_i表示当前时刻第i个存活机动目标的机动频率。

所述的机动目标传输控制矩阵如下：

其中，U_i表示传输控制矩阵。

将每个存活目标当前时刻的上一时刻的存活概率乘以0.99，得到当前时刻该机动目标的存活概率预测值。

将当前时刻每个新生机动目标状态参数中的状态值、存活概率，分别作为该新生机动目标的状态预测值和存活概率预测值。

将所有当前时刻的存活目标和所有新生目标的存活概率预测值之和，作为当前时刻机动目标的预测目标总数。

按照下式，计算当前时刻每个机动目标的过程噪声矩阵：

其中，Q_i表示当前时刻第i个存活的机动目标的过程噪声矩阵，*表示相乘操作，

表示当前时刻第i个存活的机动目标的机动频率的倒数的方差，q表示当前时刻测量环境导致的随机扰动误差。

使用修正因子的计算方法，计算当前时刻修正的每个存活的机动目标预测误差修正因子，来修正状态误差协方差，以便运动目标状态突变时能及时的调整滤波增益矩阵的值。

所述的修正因子的计算方法如下：

第1步，按照下式，根据当前时刻每个存活的机动目标的量测值，计算当前时刻每个存活的机动目标的量测矩阵：

其中，H_n表示当前时刻第n个的机动目标量测矩阵，x_n表示当前时刻第n个的机动目标量测值的横坐标值，y_n表示当前时刻第n个机动目标量测值的纵坐标值，

表示求平方根操作，其中，n＝1,2,...,N，N＝J+W，N表示当前时刻的机动目标总数，W表示新生目标总数；

第2步，将当前时刻每个存活的机动目标的量测矩阵与当前时刻该存活的机动目标的状态预测值相乘，得到当前时刻该存活的机动目标的量测预测值：

第3步，求当前时刻每个存活的机动目标的量测值和当前时刻该存活的机动目标的量测预测值的差值，并对该差值求期望值，得到当前时刻该存活的机动目标的平均残差；

第4步，按照下式，计算当前时刻每个存活的机动目标的状态误差协方差矩阵的中间修正因子：

其中，λ_i表示当前时刻第i个存活的机动目标的状态误差协方差矩阵的中间修正因子，tra表示矩阵求迹操作，V_i表示当前时刻第i个机动目标的平均残差，R表示根据雷达的测量精度得到的量测噪声，H_i表示当前时刻第i个存活的机动目标的量测矩阵，T表示矩阵转置操作，P′_C表示当前时刻的上一时刻存活的第C个机动目标的预测误差协方差矩阵；

第5步，取中间修正因子与1两者中的最大值，作为当前时刻每个存活的机动目标的状态误差协方差矩阵的修正因子。

利用状态误差协方差矩阵公式，计算当前时刻每个存活的机动目标的预测误差协方差矩阵。

所述修正的状态误差协方差矩阵公式如下：

P_i ^M＝Q_i+λ′_iF_iP′_CF_i ^T

其中，P_i ^M表示修正的当前时刻存活的第i个机动目标状态的预测误差协方差矩阵，λ′_i表示当前时刻存活的第i个机动目标的状态预测误差方差的修正因子。

利用状态误差协方差矩阵公式，设定公式中的修正因子的值为1，计算当前时刻每个存活的机动目标实际的预测误差协方差矩阵。

步骤3，更新机动目标的状态参数。

利用修正的量测距离误差公式，计算当前时刻每个机动目标的量测距离误差值。

所述的量测距离误差公式如下：

其中，g′_n表示当前时刻第n个机动目标的量测距离误差值，z_n表示当前时刻第n个机动目标的量测值，d_n表示当前时刻第n个机动目标的量测预测值，-1表示矩阵求逆操作。

将修正的量测距离误差公式中的修正的状态误差协方差矩阵替换成实际的状态误差协方差，计算当前时刻每个机动目标的实际的量测距离误差。

按照下式，计算修正的当前时刻每个机动目标量测裁剪门限：

其中，D′_n表示修正的当前时刻第n个机动目标的量测裁剪门限，g_n表示当前时刻第n个机动目标的实际的量测距离误差，g′_n表示当前时刻第n个机动目标的量测距离误差值，D_n表示当前时刻第n个机动目标的量测裁剪门限，其中，n＝1,2,...,N，N＝J+W，N表示当前时刻的机动目标总数，W表示新生目标总数。

按照下式，计算修正的当前时刻每个机动目标机动频率：

其中，α′_n表示修正的当前时刻第n个机动目标的机动频率，α_n表示当前时刻第n个机动目标的机动频率。

剔除量测距离误差值大于修正的量测裁剪门限的目标。

利用增益矩阵公式，计算裁剪后当前时刻每个机动目标的增益矩阵。

所述的增益矩阵公式如下：

其中，K_k表示当前时刻第k个机动目标的增益矩阵，k＝1,2,...,Y，Y表示量测裁剪后机动目标的总数。

利用更新状态值公式，计算裁剪后当前时刻每个机动目标的更新状态值。

所述的更新状态值公式如下：

其中，

表示当前时刻第k个机动目标的更新状态值。

利用更新存在概率公式，计算当前时刻所有的机动目标的更新存活概率。

所述的更新存活概率方法如下：

第1步，按照下式，计算当前时刻每个机动目标存活概率的中间值：

其中，y_k表示第k个机动目标的存活概率的中间值，e表示以自然常数为底的指数操作；

第2步，按照下式，计算当前时刻每个机动目标的更新存活概率：

其中，w_k表示当前时刻第k个机动目标的更新存活概率，w′_k表示当前时刻第k个机动目标的预测存活概率，pd表示雷达检测到目标可能性的检测概率，η表示当前时刻的虚假量测个数，Σ表示求和操作。

利用更新误差协方差矩阵公式，计算当前时刻所有的机动目标的更新误差协方差矩阵。

所述的更新误差协方差矩阵的计算方法如下：

其中，P_k表示第k个机动目标的更新误差协方差矩阵，I表示单位矩阵，K_k表示当前时刻第k个机动目标的增益矩阵，H_k表示第k个机动目标的量测矩阵，

表示第k个机动目标的状态预测协方差矩阵其中。

步骤4，对机动目标进行裁剪与合并。

剔除当前时刻所有机动目标中更新权值小于0.001的机动目标。

合并当前时刻所有机动目标中，欧氏距离小于合并门限的目标。

步骤5，提取当前时刻裁剪与合并之后机动目标的状态参数。

提取当前时刻裁剪与合并之后所有机动目标存活概率大于0.5的目标数目，作为当前时刻的真实目标数。

提取当前时刻裁剪与合并之后所有机动目标存活概率大于0.5的目标更新状态值，作为当前时刻真实存在的目标状态值。

提取当前时刻裁剪与合并之后所有机动目标存活概率大于0.5的目标更新存活概率，作为当前时刻真实存在的更新存活概率。

提取当前时刻裁剪与合并之后所有机动目标存活概率大于0.5的目标更新协方差矩阵，作为当前时刻真实存在的目标协方差矩阵值。

步骤6，使用平均状态误差方法，得到当前时刻的所有机动目标的平均状态误差。

所述的平均状态误差方法如下：

第1步，按照下式，计算当前时刻的状态误差值中间值：

其中，

表示当前时刻第k个机动目标的状态预测值和当前时刻第k个机动目标的更新状态值的距离误差中间值，

表示当前时刻第k个机动目标更新状态值；

第2步，按照下式，计算当前时刻的所有机动目标的平均状态误差：

其中，s表示当前时刻的所有机动目标的平均状态误差，min表示求最小值操作，c表示根据雷达测量误差确定的常数，p表示雷达测量的机动目标状态向量的维数。

以下结合仿真实验，对本发明技术效果进行进一步说明：

1.仿真条件

本发明在Intel(R)Pentium(R)2CPU@3.00GHz 3.00GHz处理器的电脑上，采用MATLAB R2014a软件完成仿真。

仿真场景设置：为了验证本发明提出的自适应GM-PHD方法能精确跟踪多个机动目标，本发明的仿真实验场景为在的二维空间[0,6×10⁴]×[0,2.5×10⁴]内有4个做变速运动的目标，相继出现和消失在不同的时间整个时间持续120s，目标1存在时间是t＝1s到t＝120s，目标2存在时间是从t＝10s到t＝110s，目标3存在时间是从t＝15s到t＝120s，目标4存在时间是从t＝1s到t＝100s。附图2(a)为目标真实轨迹，其中横轴表示目标的横坐标，纵轴表示目标的纵坐标，*表示目标的实际位置。

定义目标状态向量

分别表示目标x方向距离、x方向速度、x方向加速度，y方向距离、y方向速度、y方向加速度,给定其中两个存活目标初始状态为：

m₁＝[47000m,-426m/s,0,12000m,0,0]^T，

m₂＝[5000m,100m/s,0,12000m,0,0]^T，

状态协方差矩阵P₁＝diag[100,2,1,100,2,1]，P₂＝diag[100,2,1,100,2,1]，存活概率w₁和w₂都为0.99。

第10s出现的新生目标的初始状态为：

bm1＝[40000m,100m/s,-5m/s²,13000m,100,15m/s²]^T

协方差矩阵bP₁＝diag[100,2,1,100,2,1]，存活概率bw₁为0.2。

第15s出现的新生目标的初始状态为：

bm2＝[6000m,300m/s,0,20000m,-100,0]^T

协方差矩阵bP₂＝diag[100,2,1,100,2,1]，存活概率bw₂为0.2。

设置仿真参数：采样时间间隔T＝1s,测量时间N＝120s,高斯修剪阈值T_prun＝10^-4，量测裁剪门限γ＝16，高斯合并门限U_merg＝4，最大高斯分量数J_max＝40，目标检测概率pd＝0.98，杂波平均数rc＝20，机动频率a＝0.05，目标加速度上限a_max＝20m/s²，量测噪声矩阵R＝diag[20,20]^T。

为了证明仿真效果，进行100次的蒙特卡罗仿真，计算估计的目标数目均值和平均状态误差OSPA，OSPA的参数为c＝200，p＝2。附图2(b)是目标量测结果图，横轴表示量测点的横坐标，纵轴表示量测点的纵坐标，*表示目标和杂波共同组成的量测结果。附图2(c)是采用EK-GMPHD跟踪方法的跟踪结果和实际轨迹，图2(c)中的横轴表示目标的横坐标，纵轴表示目标的纵坐标，黑点表示EK-GMPHD跟踪方法的跟踪结果，实线表示实际轨迹。附图2(d)是采用本发明跟踪方法的跟踪结果，图2(d)中的横轴表示跟踪结果的横坐标，纵轴表示跟踪结果的纵坐标，轨迹中的黑点表示本发明跟踪方法的跟踪结果，实线表示实际轨迹。附图2(e)是EK-GMPHD跟踪方法平均估计目标数目，图2(e)中的横轴表示时间，纵轴表示目标数目，黑圈o表示应用EK-GMPHD跟踪方法平均估计目标数目，实线表示实际目标个数。附图2(f)是本发明跟踪方法平均估计目标数目，图2(f)中的横轴表示时间，纵轴表示目标数目，黑圈o表示应用本发明跟踪方法平均估计目标数目，实线表示实际目标个数。附图2(g)是采用本发明跟踪方法和EK-GMPHD跟踪方法得到的目标跟踪误差，图2(g)中的横轴表示时间，纵轴表示平均OSPA误差，加*曲线表示本发明跟踪方法的目标跟踪误差，另一条曲线表示EK-GMPHD跟踪方法得到的目标跟踪误差。

2.仿真内容与结果分析：

附图2(b)给出了本仿真实验的量测结果图，由图2(b)可知，本仿真实验设置的量测回波个数多，也就是虚假量测个数较多，符合真实复杂环境。附图2(c)给出了本仿真实验中EK-GMPHD跟踪方法的跟踪结果，由图2(c)可知，EK-GMPHD跟踪方法的在目标状态发生突变时跟踪结果与实际轨迹相差较大。附图2(d)给出了本仿真实验中本发明的跟踪方法的跟踪结果，由图2(d)可知，本发明的跟踪方法对目标的状态是处于机动状态还是非机动状态还是出现状态突变，都能实现多目标的准确跟踪。附图2(e)给出了EK-GMPHD跟踪方法目标数的估计均值，由图2(e)可知，该方法在估计多目标个数时，由于不能及时应对目标状态突变或出现新生目标，有多个时刻存在个数估计误差。附图2(f)给出了本发明跟踪方法目标数的估计均值，由图2(f)可知，该方法能较准确的估计新生目标和状态突变的目标个数，只有三个时刻存在多估计或少估计的情况。附图2(g)给出了本发明跟踪方法和EK-GMPHD跟踪方法状态误差结果，由图2(g)可知，当有新生目标或者目标消失时，EK-GMPHD跟踪方法误差出现较大的起伏，但是本发明的跟踪方法虽然在估计个数存在误差的地方也表现出起伏，但是无论出现新生目标还是目标状态复杂多变的情况，相比于EK-GMPH跟踪方法误差都明显减小，并在目标全部目标出现时状态误差趋于稳定。

综上所述，从仿真效果图的分析可知，本发明提出的一种自适应EK-GMPHD跟踪方法，实现了对多个机动目标的准确跟踪。目标数估计准确度高，目标跟踪误差小，跟踪性能良好，同时由于采用了量测门限对量测值进行裁剪，使得目标状态更新后高斯分量数大大减少，降低了计算的复杂度和储存空间，并且通过自适应的方法改变门限值和机动频率的大小，能适应各种复杂的运动形式和量测环境，使该发明在现实工程应用中更有优势。

Claims (8)

1.一种自适应高斯混合概率假设密度滤波GM-PHD的机动目标跟踪方法，其特征在于，利用修正因子调整当前时刻的滤波增益，生成机动变化率和量测门限，利用机动变化率自适应地调整机动频率和量测门限值，该方法的步骤包括如下：

(1)读取机动目标的状态值：

(1a)读取当前时刻仍然存活的所有机动目标在上一时刻的状态参数和量测值；

(1b)读取当前时刻新生机动目标的状态参数和量测值；

(2)预测当前时刻每个机动目标的状态参数：

(2a)按照下式，预测每个存活机动目标在当前时刻的状态值：

m_i＝F_im′_C+U_ia_i

其中，m_i表示当前时刻第i个存活的机动目标的状态预测值，F_i表示由雷达扫描周期和机动目标的机动频率获得的当前时刻第i个存活的机动目标状态转移矩阵，m′_C表示当前时刻的上一时刻第C个存活的机动目标的状态值，U_i表示由雷达扫描周期和机动目标的机动频率获得的当前时刻第i个存活的机动目标传输控制矩阵，a_i表示当前时刻第i个存活的机动目标加速度，其中，i＝1,2,...,J,C＝1,2,...,M，J表示当前时刻存活的机动目标总数，M表示上一时刻存活的机动目标总数，且i和C的取值对应相等，J和M的取值相等；

(2b)将每个存活目标当前时刻的上一时刻的存活概率乘以0.99，得到当前时刻该机动目标的存活概率预测值；

(2c)将当前时刻每个新生机动目标状态参数中的状态值、存活概率，分别作为该新生机动目标的状态预测值和存活概率预测值；

(2d)将所有当前时刻的存活目标和所有新生目标的存活概率预测值之和，作为当前时刻机动目标的预测目标总数；

(2e)按照下式，计算当前时刻每个机动目标的过程噪声矩阵；

其中，Q_i表示当前时刻第i个存活的机动目标的过程噪声矩阵，*表示相乘操作，α_i表示当前时刻第i个存活的机动目标的机动频率，

表示当前时刻第i个存活的机动目标的机动频率的倒数的方差，q表示当前时刻测量环境导致的随机扰动误差；

(2f)使用修正因子的计算方法，计算当前时刻修正的每个存活的机动目标预测误差修正因子；

(2g)利用修正的状态误差协方差矩阵公式，计算当前时刻每个存活的机动目标修正的预测误差协方差矩阵；

(2h)利用状态误差协方差矩阵公式，设定公式中的修正因子的值为1，计算当前时刻每个存活的机动目标实际的预测误差协方差矩阵；

(3)更新机动目标的状态参数：

(3a)利用修正的量测距离误差公式，计算当前时刻每个机动目标的量测距离误差值；

(3b)用实际的状态误差协方差替换修正的状态误差协方差矩阵，计算当前时刻每个机动目标的实际的量测距离误差；

(3c)按照下式，计算修正的当前时刻每个机动目标量测裁剪门限：

其中，D_n′表示修正的当前时刻第n个机动目标的量测裁剪门限，g_n表示当前时刻第n个机动目标的实际的量测距离误差，g_n′表示当前时刻第n个机动目标的量测距离误差值，D_n表示当前时刻第n个机动目标的量测裁剪门限，其中，n＝1,2,...,N，N＝J+W，N表示当前时刻的机动目标总数，W表示新生目标总数；

(3d)按照下式，计算修正的当前时刻每个机动目标机动频率：

其中，α_n′表示修正的当前时刻第n个机动目标的机动频率，α_n表示当前时刻第n个机动目标的机动频率；

(3e)剔除量测距离误差值大于修正的量测裁剪门限的目标；

(3f)利用增益矩阵公式，计算裁剪后当前时刻每个机动目标的增益矩阵：

(3g)利用更新状态值公式，计算裁剪后当前时刻每个机动目标的更新状态值：

(3h)利用更新存在概率公式，计算当前时刻所有的机动目标的更新存活概率；

(3i)利用更新误差协方差矩阵公式，计算当前时刻所有的机动目标的更新误差协方差矩阵；

(4)对机动目标进行裁剪与合并：

(4a)剔除当前时刻所有机动目标中更新权值小于0.001的机动目标；

(4b)合并当前时刻所有机动目标中，欧氏距离小于合并门限的目标；

(5)提取当前时刻裁剪与合并之后机动目标的状态参数：

(5a)提取当前时刻裁剪与合并之后所有机动目标存活概率大于0.5的目标数目，作为当前时刻的真实目标数；

(5b)提取当前时刻裁剪与合并之后所有机动目标存活概率大于0.5的目标更新状态值，作为当前时刻真实存在的目标状态值；

(5c)提取当前时刻裁剪与合并之后所有机动目标存活概率大于0.5的目标更新存活概率，作为当前时刻真实存在的更新存活概率；

(5d)提取当前时刻裁剪与合并之后所有机动目标存活概率大于0.5的目标更新协方差矩阵，作为当前时刻真实存在的目标协方差矩阵值。

2.根据权利要求1所述的自适应高斯混合概率假设密度滤波GM-PHD的机动目标跟踪方法，其特征在于，步骤(2f)中所述修正因子的计算方法如下：

第1步，按照下式，根据当前时刻每个存活的机动目标的量测值，计算当前时刻每个存活的机动目标的量测矩阵：

其中，H_n表示当前时刻第n个的机动目标量测矩阵，x_n表示当前时刻第n个的机动目标量测值的横坐标值，y_n表示当前时刻第n个机动目标量测值的纵坐标值，

表示求平方根操作；

第2步，将当前时刻每个存活的机动目标的量测矩阵与当前时刻该存活的机动目标的状态预测值相乘，得到当前时刻该存活的机动目标的量测预测值：

第3步，求当前时刻每个存活的机动目标的量测值和当前时刻该存活的机动目标的量测预测值的差值，并对该差值求期望值，得到当前时刻该存活的机动目标的平均残差；

第4步，按照下式，计算当前时刻每个存活的机动目标的状态误差协方差矩阵的中间修正因子：

其中，λ_i表示当前时刻第i个存活的机动目标的状态误差协方差矩阵的中间修正因子，tra表示矩阵求迹操作，V_i表示当前时刻第i个机动目标的平均残差，R表示根据雷达的测量精度得到的量测噪声，H_i表示当前时刻第i个存活的机动目标的量测矩阵，T表示矩阵转置操作，P_C′表示当前时刻的上一时刻存活的第C个机动目标的预测误差协方差矩阵；

第5步，取中间修正因子与1两者中的最大值，作为当前时刻每个存活的机动目标的状态误差协方差矩阵的修正因子。

3.根据权利要求2所述的自适应高斯混合概率假设密度滤波GM-PHD的机动目标跟踪方法，其特征在于，步骤(2g)中所述修正的状态误差协方差矩阵公式如下：

P_i ^M＝Q+λ_i′FP_C′F^T

其中，P_i ^M表示修正的当前时刻存活的第i个机动目标状态的预测误差协方差矩阵，λ_i′表示当前时刻存活的第i个机动目标的状态预测误差方差的修正因子。

4.根据权利要求2所述的自适应高斯混合概率假设密度滤波GM-PHD的机动目标跟踪方法，其特征在于，步骤(3a)中所述的修正的量测距离误差公式如下：

其中，g_n′表示当前时刻第n个机动目标的量测距离误差值，z_n表示当前时刻第n个机动目标的量测值，d_n表示当前时刻第n个机动目标的量测预测值，-1表示矩阵求逆操作。

5.根据权利要求1所述的自适应高斯混合概率假设密度滤波GM-PHD的机动目标跟踪方法，其特征在于，步骤(3f)中所述的增益矩阵公式如下：

其中，K_k表示当前时刻第k个机动目标的增益矩阵，k＝1,2,...,Y，Y表示量测裁剪后机动目标的总数。

6.根据权利要求5所述的自适应高斯混合概率假设密度滤波GM-PHD的机动目标跟踪方法，其特征在于，步骤(3g)中所述的更新状态值公式如下：

其中，m_k表示当前时刻第k个机动目标的更新状态值。

7.根据权利要求5所述的自适应高斯混合概率假设密度滤波GM-PHD的机动目标跟踪方法，其特征在于，步骤(3h)中所述的更新存活概率方法如下：

第1步，按照下式，计算当前时刻每个机动目标存活概率的中间值：

其中，y_k表示第k个机动目标的存活概率的中间值，e表示以自然常数为底的指数操作；

第2步，按照下式，计算当前时刻每个机动目标的更新存活概率：

其中，w_k表示当前时刻第k个机动目标的更新存活概率，w′_k表示当前时刻第k个机动目标的预测存活概率，pd表示雷达检测到目标可能性的检测概率，η表示当前时刻的虚假量测个数，Σ表示求和操作。

8.根据权利要求1所述的自适应高斯混合概率假设密度滤波GM-PHD的机动目标跟踪方法，其特征在于，步骤(3i)中所述的更新误差协方差矩阵的计算方法如下：

其中，P_k表示第k个机动目标的更新误差协方差矩阵，I表示单位矩阵，K_k表示当前时刻第k个机动目标的增益矩阵，H_k表示第k个机动目标的量测矩阵，

表示第k个机动目标的状态预测协方差矩阵其中。

Images