一种基于动态多维分配的被动多源多目标跟踪方法与流程

本发明涉及基于动态多维分配的被动多源多目标跟踪方法。

背景技术：

在多传感器多目标跟踪领域中，确定传感器中某一观测来源于哪个目标是一个很有挑战性的问题，即数据关联问题。从上世纪60年代开始对此进行了广泛的研究，并发展出了一系列算法，这些算法在复杂度和跟踪性能上各不相同。在多被动传感器多目标跟踪应用场景下，由于传感器只能获得目标的角度观测数据，无法获得完整的目标位置信息，此时的纯角度数据关联无疑具有较高的挑战性。

近年来，解决多传感器多目标数据关联问题的方法有最近邻法(NN)，联合概率数据关联(JPDA)以及多假设跟踪(MHT)等。其中最近邻法将统计意义上与被跟踪目标预测位置最近的观测作为关联的观测，仅适用于信噪比高目标密度小的条件，在实际应用中效率低；联合数据关联通过计算量大，适合杂波少的环境；而多假设跟踪试图考虑所有关联假设达到最优跟踪的效果，但是由于关联假设数量随时间指数增长使得MHT的直接实现成为一件不可能完成的任务。除此之外，多维分配(MDA)技术是解决数据关联问题的另一种方案。传统的多维分配技术处理被动多传感器跟踪问题时首先用多维(S-D)分配解决观测-观测之间的关联，利用得到的关联组合计算目标位置的最大似然估计，用二维(2-D)分配解决观测-航迹间的关联，用关联上的量测更新航迹，该方法需考虑所有的关联假设，计算量大。

被动多源多目标跟踪系统需要对获得角度观测进行数据关联，确定相关目标的位置和速度。在Musicki D.Multi-target tracking using multiple passive bearings-only asynchronous sensors.IEEE Transactions on Aerospace and Electronic Systems,2008,44(3):1151-1160.中角度观测被距离参数化，并假定各传感器每时刻只获得一个服从高斯混合分布的角度观测，用这些观测起始和更新航迹。这种方法的弊端在于假定每时刻只有一个角度观测，大大简化了数据关联的难度。在Hanselmann T,Morelande M.Multiple target tracking with asynchronousbearings-only measurements.In Proceedings of the International Conference onInformation Fusion,Quebec,Cananda.2007:1-8.中使用距离参数化技术并利用第一次扫描所得观测初始化目标，后续的角度观测用于枚举所有可能的观测来源假设，基于这些假设给出后验目标状态。该算法的弊端在于未考虑杂波观测，且枚举所有观测来源假设会产生高的时间复杂度。在高祺.多传感器多目标跟踪的数据关联算法研究[D].西安:西安电子科技大学硕士学位论文,2013:36-42.中使用了传统的多维分配算法，该算法首先用多维分配(S-D分配)进行观测-观测之间的关联，用关联上的观测计算目标位置的最大似然估计，再利用二维分配(2-D分配)进行观测-航迹之间的关联，算法的弊端在于需要对所有的关联组合进行假设，算法的时间复杂度高，且未充分利用航迹的先验信息，数据的关联准确度不高。由此，在实际应用中还需一种兼顾数据关联准确度和时间复杂度的方案。

技术实现要素：

本发明是为了解决现有技术关联准确度低和算法时间复杂度高的问题，而提出的一种基于动态多维分配的被动多源多目标跟踪方法。

一种基于动态多维分配的被动多源多目标跟踪方法按照以下步骤实现：

步骤一：利用k-1时刻的目标的状态，建立对应于该目标的航迹p的预选波门

步骤二：利用步骤一构建的预选波门以及k-1时刻的航迹p和各传感器的观测值组合构建代价函数并依据航迹p和观测值之间一一对应的关系构建二值变量

步骤三：利用步骤二构建的代价函数和二值变量构建全局关联代价函数，得到(S+1)-D分配公式，并给出约束条件；所述S表示传感器的个数，D表示维度；

步骤四：利用拉格朗日松弛算法对步骤三的约束条件进行松弛，对(S+1)-D分配公式进行降维处理，得到二维分配公式；

步骤五：利用广义拍卖算法计算步骤四所得的二维分配公式的对偶解；

步骤六：根据步骤五对约束条件进行实施，初始化次梯度向量，利用次梯度向量对拉格朗日乘子更新；

步骤七：迭代执行步骤四至步骤六，获得(S+1)-D分配公式的代价值J_S和使分配结果的原始解后，得到航迹p和对应观测值的分配组合；

步骤八：根据步骤七得到的观测值的分配组合对目标当前时刻的位置利用似然函数进行最大似然估计；

步骤九：利用步骤八得到的位置的最大似然估计，根据卡尔曼滤波方法估计目标状态，用状态估计值更新航迹，实现多目标跟踪。

当维度大于等于3时为多维，大于等于3个数据来源为多源，大于等于2个目标为多目标。

发明效果：

本发明针对被动多源多目标跟踪问题，提出了一种基于动态多维分配的被动多源多目标跟踪方法。该技术使用航迹的先验航迹信息产生预选波门，对航迹和各传感器落入航迹波门内的观测值产生关联假设，利用航迹和传感器观测值组合构建(S+1)-D分配公式，对(S+1)-D分配公式进行约束松弛，求得航迹和观测值组合的关联结果，用关联上的观测值组合更新当前航迹，避免了二次求解二维分配问题。和现有方法相比，本发明利用航迹先验信息提高目标跟踪性能，采用的预选波门减少了关联假设的产生，降低算法的时间复杂度，数据关联正确率与基于多维分配的S-D/2-D方法相比提高了约5％，实现了多目标的有效跟踪。

附图说明

图1为多目标跟踪航迹图；

图2为多目标位置RMSE图；

图3为虚假航迹数示意图；

图4为成功跟踪目标数示意图；

图5为(S+1)-D分配算法的流程图。

具体实施方式

具体实施方式一：如图5所示，一种基于动态多维分配的被动多源多目标跟踪方法包括以下步骤：

步骤一：利用k-1时刻的目标的状态，建立对应于该目标的航迹p的预选波门

步骤二：利用步骤一构建的预选波门以及k-1时刻的航迹p和各传感器的观测值组合构建代价函数并依据航迹p和观测值之间一一对应的关系构建二值变量

步骤三：利用步骤二构建的代价函数和二值变量构建全局关联代价函数，得到(S+1)-D分配公式，并给出约束条件；所述S表示传感器的个数，D表示维度；

步骤四：利用拉格朗日松弛算法对步骤三的约束条件进行松弛，对(S+1)-D分配公式进行降维处理，得到二维分配公式；

步骤五：利用广义拍卖算法计算步骤四所得的二维分配公式的对偶解；

步骤六：根据步骤五对约束条件进行实施，初始化次梯度向量，利用次梯度向量对拉格朗日乘子更新；

步骤七：迭代执行步骤四至步骤六，获得(S+1)-D分配公式的代价值J_S和使分配结果的原始解后，得到航迹p和对应观测值的分配组合；

步骤八：根据步骤七得到的观测值的分配组合对目标当前时刻的位置利用似然函数进行最大似然估计；

步骤九：利用步骤八得到的位置的最大似然估计，根据卡尔曼滤波方法估计目标状态，用状态估计值更新航迹，实现多目标跟踪。

本文发明所研究问题主要涉及航迹-量测之间的关联。在使用S个被动传感器的情况下，采用动态多维分配的方案实现数据关联以及多目标跟踪。其中首先使用航迹的先验航迹信息产生预选波门，对航迹和各传感器落入航迹波门内的观测值产生关联假设，利用航迹和传感器观测值组合构建(S+1)-D分配公式，对该(S+1)-D分配公式依次进行约束松弛和约束实施，找出对应于当前航迹的观测值组合，对得到的观测值组合利用极大似然估计算法获得完整的目标位置数据，用位置的最大似然估计更新当前航迹。在传感器数目S≥3时，多维分配问题被判为NP-hard问题，解决此问题核心是合理利用拉格朗日松弛算法对多维分配问题进行降维并寻找满足实时性要求的次优解。

具体实施方式二：本实施方式与具体实施方式一不同的是：所述步骤一中利用k-1时刻的目标的状态，建立对应于该目标的航迹p的预选波门的具体过程为：

其中表示k-1时刻目标位置的估计值，表示k-1时刻目标速度的估计值，h_s(X)是传感器s的观测方程，是h_s(X)关于X的偏导数，是k时刻状态转移函数，和分别是k-1时刻目标状态和状态协方差矩阵估计值，和分别是k时刻目标状态和状态协方差矩阵的预测值，通过对k-1时刻目标状态进行预测可得到k时刻目标状态的预测值是k-1时刻过程噪声协方差矩阵，和分别是观测和观测协方差矩阵的预测值，是目标在k时刻的预测位置，(x_s,y_s)是传感器s的位置，h_s(X)是传感器s的量测方程，z是传感器s的k时刻的观测值，是观测噪声方差；τ是预先设定的阈值，用于决定观测值落入预选波门概率的大小；集合中不等式用于判断传感器s的观测是否落入航迹p的预选波门，所有满足预选波门条件的观测值构成集合

其它步骤及参数与具体实施方式一相同。

具体实施方式三：本实施方式与具体实施方式一或二不同的是：所述步骤二中利用步骤一构建的预选波门以及k-1时刻的航迹p和各传感器的观测值组合构建代价函数并依据航迹p和观测值之间一一对应的关系构建二值变量的具体过程为：

其中代价函数表示观测组合中的观测值来源于目的代价，表示观测组合来源于目标的概率，表示观测组合来源于伪信号源的概率，表示空集，X_p为k时刻目标状态的真实值，为k时刻目标状态的估计值，用k时刻目标状态的预测值代替是观测向量基于目标状态估计值的条件均值，是对应的协方差矩阵，H是S个传感器组合的堆叠观测函数，由单个传感器的观测方程组成，H_X是堆叠观测函数H的偏导数，是观测向量的条件概率密度函数，假定服从高斯分布，是传感器s的探测概率，u(i_s)是二值函数，ψ_s是传感器s监督区域的体积，则代价函数具体表示为：

依据航迹p和观测值间一一对应的关系，二值变量的具体公式为：

其它步骤及参数与具体实施方式一或二相同。

具体实施方式四：本实施方式与具体实施方式一至三之一不同的是：所述步骤三中利用步骤二构建的代价函数和二值变量构建全局关联代价函数，得到(S+1)-D分配公式，并给出约束条件的具体过程为：

其它步骤及参数与具体实施方式一至三之一相同。

具体实施方式五：本实施方式与具体实施方式一至四之一不同的是：所述步骤四中利用拉格朗日松弛算法对步骤三的约束条件进行松弛，对(S+1)-D分配公式进行降维处理，得到二维分配公式的具体过程为：

其中u_r,r＝S+1,S,S-1,S-2,...,3为无约束的拉格朗日乘子，利用拉格朗日乘子将S+1个约束条件中的后S+1-r个约束松弛，得到r维分配公式，满足前r个约束条件；和为对应于该公式的代价函数和二值变量，r维分配公式的具体公式为：

对约束条件依次松弛，直到r＝2时获得二维分配公式。

其它步骤及参数与具体实施方式一至四之一相同。

具体实施方式六：本实施方式与具体实施方式一至五之一不同的是：所述步骤五中利用广义拍卖算法计算步骤四所得的二维分配公式的对偶解的具体过程为：

利用拍卖算法完成对二维分配公式的处理后，给出下述定义：R＝2；γ_k←{k}；m＝R，n₁表示传感器1观测值的数目，J₂表示二维分配公式单次迭代的对偶解；若则表示γ_k和i_m成功关联，γ_k表示利用拍卖算法得到的分配组合，i_m表示第m个传感器观测值的序号，i_m＝1,…,n_m，n_m为传感器m观测值的数目(共有S个传感器，m为传感器1到传感器S中的任意一个传感器)；表示二者无关；初始化对于分配结果的集合，进行操作依次把取值为1的分配组合加入新的集合存储。

按拍卖算法的理论，将γ_k作为投标人，i_r作为物品，投标人与物品匹配的价值是算法的目的就是找出投标人和物品之间的最佳组合，使得总的价值最大，定义分配给第γ_k个投标人的物品集合为A(γ_k)，所有分配给第γ_k的物品与第γ_k个投标人的二元组集合为B，

定义分配T为人与物品的二元对集合，满足：对任意(γ_k,i_r)∈T都有i_r∈A(γ_k)，对每个投标人γ_k至多存在一对(γ_k,i_r)∈T，对每个物品i_r至多存在一对(γ_k,i_r)∈T；给定分配T，若存在一对(γ_k,i_r)∈T，则称第γ_k个投标人被分配，否则称未被分配；若每个投标人和每个物品都被分配，则就称该分配为可行分配或者完全分配，否则称该分配为偏分配；

定义价值向量以及松弛条件参数ε，用于限定投标人γ_k和物品i_r在何种程度上是最佳分配对，若

则称物品i_r在ε范围内是投标人γ_k的最佳物品，若对任意(γ_k,i_r)∈T都有

则称T和满足ε补松弛条件；拍卖算法不断迭代，直至获得完全分配后终止；迭代从满足ε补松弛条件的偏分配和价值向量开始；

先进行投标阶段；设分配T中未被分配的投标人构成集合I，对任意i∈I，寻找最佳物品使

和对应值

并寻找物品以外其他物品提供的最佳值

计算第γ_k个投标人的标值

进行分配阶段；每个物品i_r可能收到多个投标人的投标，记投标人的集合为P(i_r)，若P(i_r)非空，记最高投标为即

从分配T中去掉(γ_k,i_r)对，加上对，其中是P(i_r)中取得上述最大值的投标人；重复上述投标阶段和分配阶段，直到每个投标人和每个物品都被分配，获得完全分配，得到(γ_k,i_r)的分配组合。

其它步骤及参数与具体实施方式一至五之一相同。

具体实施方式七：本实施方式与具体实施方式一至六之一不同的是：所述步骤六中根据步骤五对约束条件进行实施，初始化次梯度向量，利用次梯度向量对拉格朗日乘子更新的具体过程为：

初始化次梯度向量其具体计算公式如下：

g_(r+1)j＝g_(r+1)j-1

其中为通过二维分配得到的组合；次梯度向量的具体意义为：当取值为0时，表明r+1维分配问题中的i_r+1没有违反约束，取值为1表明i_r+1没有被分配，小于0表明i_r+1进行了多次分配。在r≥3时更新r维分配公式的拉格朗日乘子，具体公式如下：

其中是每次迭代中r维分配公式的最优代价值，f_dual是各次迭代的最优对偶值，即若R<S+1，令R＝R+1,r＝R，重复步骤五至步骤六直到R＝S+1。

其它步骤及参数与具体实施方式一至六之一相同。

具体实施方式八：本实施方式与具体实施方式一至七之一不同的是：所述步骤七中重复执行步骤四至步骤六，获得(S+1)-D分配公式的代价值J_S和使分配结果的原始解后，得到航迹p和对应观测值的分配组合的具体过程为：

获得(S+1)-D分配公式的代价值J_S和使分配结果的原始解后，对此次迭代结果是够满足要求进行判断；求最优原始值f_primal＝min(f_primal,J_S)及其与最优对偶值之间的差值gap＝(f_primal-f_dual)/|f_primal|，若gap<mingap，则迭代终止，其中mingap为预先设置的阈值，取值为0.01至0.05；最优原始解即为各传感器观测数据的关联结果，即为最终的分配结果，即航迹p和与之对应的观测值的组合，迭代终止；否则继续迭代步骤四至步骤七的过程，直至迭代终止。

其它步骤及参数与具体实施方式一至七之一相同。

具体实施方式九：本实施方式与具体实施方式一至八之一不同的是：所述步骤八中根据步骤七得到的观测值的分配组合对目标当前时刻的位置利用似然函数进行最大似然估计的具体过程为：

其中满足高斯分布，使似然函数取值最大的X即为目标状态的最大似然估计。

其它步骤及参数与具体实施方式一至八之一相同。

具体实施方式十：本实施方式与具体实施方式一至九之一不同的是：所述步骤九中利用步骤八得到的位置的最大似然估计，根据卡尔曼滤波方法估计目标状态，用状态估计值更新航迹，实现多目标跟踪的具体过程为：

P_k|k-1＝FP_k-1|k-1F^T+ΓQ_k-1Γ^T

S_k＝ΨP_k|k-1Ψ^T+R_k

所述Ψ为卡尔曼滤波中使用的观测矩阵，为k时刻目标观测预测值，P_k|k-1为k时刻目标状态预测误差协方差矩阵，Q_k-1为k-1时刻目标过程噪声协方差矩阵，Γ为过程噪声分布矩阵，S_k为新协方差矩阵，R_k为观测噪声协方差矩阵，K_k为卡尔曼滤波增益，Z_k为k时刻目标观测值，为k时刻目标状态估计值，P_k|k为k时刻目标状态误差协方差矩阵估计值。

其它步骤及参数与具体实施方式一至九之一相同。

实施例一：

由图1看出，利用本发明对多目标的跟踪效果好，且产生的虚假航迹数少。由图2可看出，目标位置的RMSE呈现出收敛的趋势，且RMSE的值维持在1之内，说明本发明有较高的目标状态估计精度。由图3可知，每时刻的虚假航迹数维持在0.4附近，这说明本发明有较好的的虚假航迹剔除能力，能将虚假航迹保持在较低水平。由图4可知，每时刻成功跟踪的目标数逐渐维持在4个，说明本发明在目标跟踪发现方面，能实现对全部目标的有效跟踪。在未损失目标跟踪发现能力的条件下，维持了较低的虚假航迹数，具有较强的分辨目标与虚警的能力。因此应用本发明能实现多被动传感器多目标的有效跟踪。