本发明公开了一种基于时差量测箱粒子PHD的被动跟踪多目标方法，主要解决在被动跟踪多目标中现有粒子滤波技术运算效率低下的问题，本发明实现的步骤如下：(1)获取起始跟踪时刻箱粒子；(2)获取新生箱粒子(3)合并新生箱粒子和箱粒子；(4)预测合并后的箱粒子；(5)更新预测后的箱粒子；(6)重采样更新后的箱粒子；(7)获取多目标状态；(8)判断每个无源基站所获得时差量测个数是否均为0，若是，则结束跟踪，若否，则返回步骤(2)继续跟踪。本发明通过利用时差约束传播函数更新箱粒子，使得能够用较少的箱粒子实现被动跟踪量测为时差的多目标，降低了算法的时间复杂度，提高了被动跟踪多目标的效率。
1.一种基于时差量测箱粒子PHD的被动跟踪多目标方法，其特征在于，利用箱粒子进行PHD滤波，利用时差约束传播函数更新预测后的箱粒子；该方法的具体步骤如下：(1)获取起始跟踪时刻多目标状态分布的每个箱粒子：(1a)按照下式，获取起始跟踪时刻多目标状态分布的每个点粒子：其中， 表示起始跟踪时刻多目标状态分布的第i个点粒子，N表示由跟踪场景的复杂度所确定的随机采样总数，Ψ(·)表示随机采样函数，A表示由起始时刻n0个目标状态组成的矩阵，n0表示起始跟踪时刻多目标的总数，P0表示一个对角线元素为[40,1,40,1]的对角矩阵；
(1b)利用区间扩张法，将起始跟踪时刻多目标状态分布的每个点粒子扩张成区间大小均为[240,6,240,6]T的箱粒子，其中，T表示转置操作；
其中，φ(c)表示第c个新生点粒子，γ表示由跟踪场景的复杂度所确定的新生箱粒子随机采样总数，B表示由当前跟踪时刻的前一时刻无源基站所获得的所有时差量测代入Chan定位方程得到的矩阵；
(2b)利用区间扩张法，将每个新生点粒子扩张成区间大小均为[240,6,240,6]T的新生箱粒子；
(2c)将 的值作为每个新生箱粒子的权值，其中，PB表示根据待跟踪目标所处区域新生目标数目的统计知识设定的新生目标概率；
(3)合并所有的新生箱粒子与所有的箱粒子，得到合并后的所有的箱粒子；
(4a)利用状态转移公式，预测合并后每个箱粒子的状态；
(4b)利用权值预测公式，预测合并后每个箱粒子的权值；
(5)利用时差约束传播函数对预测后的每个箱粒子进行更新：(5a)按照下式，计算预测后每个箱粒子的时差量测点粒子中每个元素值：其中， 表示预测后第y个箱粒子的时差量测点粒子中的第m个元素值，||·||表示取(y)模值操作，p 表示预测后的第y个箱粒子的中心点，Sλ表示第λ个无源基站的坐标，m＝λ-1，C表示电磁波在真空中的传播速度，则预测后第y个箱粒子的预测点粒子为 M表示无源基站的总数；
(5b)按照下式，将预测后的每个箱粒子的时差量测点粒子扩张成箱粒子：[t(y)]＝t(y)+H×([p(y)]-p(y))其中，[t(y)]表示第y个预测后箱粒子的时差量测箱粒子，H表示t(y)在p(y)处的雅克比矩阵，[p(y)]表示第y个预测后箱粒子；
(5c)利用区间扩张法，将当前跟踪时刻无源基站获得的每个时差量测扩张成区间为ε的时差量测箱粒子，其中，ε为依据无源基站的性能参数得到的最大测量误差；
(5d)按照下式，用每个时差量测箱粒子约束每个箱粒子：其中， 表示第j个时差量测箱粒子对第y个预测后箱粒子进行约束后得到的箱粒子，H-1表示矩阵H的逆矩阵，∩表示求交集操作，[zj]表示当前跟踪时刻第j个时差量测箱粒子；
(5e)将约束后箱粒子 中的最大箱粒子作为更新后的第 个箱粒子， 和y的值对应相等，其中，n′表示当前跟踪时刻时差量测箱粒子的总数；
(5f)按照下式，对每个预测后的箱粒子的权值进行更新，得到更新后的每个箱粒子的权值：其中， 表示更新后的第 个箱粒子的权值，PD表示依据无源基站的性能参数设定的检测概率，∑表示求和操作，|·|表示取箱粒子面积的操作，κ表示根据目标所处跟踪区域杂波分布统计设定的杂波强度，α(y)表示预测后第y个箱粒子的权值；
将更新后的每个箱粒子带入随机子划分重采样方程，得到重采样后的箱粒子；
(7a)对重采样后所有的箱粒子的权值求和，得到当前跟踪时刻多目标的数目估计值；
(7b)利用Kmeans聚类中的欧式距离准则，对重采样后的所有的箱粒子聚出多目标的状态；
(8)判断每个无源基站所获得时差量测个数是否均为0，若是，执行步骤(9)，否则，执行步骤(2)；
2.根据权利要求1所述的基于时差量测箱粒子PHD的被动跟踪多目标方法，其特征在于，步骤(1b)，步骤(2b)，步骤(5c)中所述的区间扩张法是指，用每个点粒子减去其对应的区间长度，将差值作为箱粒子的下限，用每个点粒子加上其对应的区间长度，将差值作为箱粒子的上限。
3.根据权利要求1所述的基于时差量测箱粒子PHD的被动跟踪多目标方法，其特征在于，步骤(4a)中所述的状态转移公式为：[p(y)]＝f([x(s)],ω)
其中，[p(y)]表示预测后的第y个箱粒子，f(·)表示一个线性函数或者非线性函数，目标的运动方式为匀速直线运动时，f(·)表示一个线性函数，目标的运动方式为除匀速直线运动之外的其他任意运动方式时，f(·)表示一个非线性函数，[x(s)]表示合并后的第s个箱粒子，y和s的值对应相等，ω表示当前跟踪时刻的过程噪声。
4.根据权利要求1所述的基于时差量测箱粒子PHD的被动跟踪多目标方法，其特征在于，步骤(4b)中所述的权值预测公式为：(y) (s)
其中，α(y)表示预测后第y个箱粒子的权值，PS表示依据目标的属性设定的存活概率，w(s)表示合并后的第s个箱粒子的权值，y与s的值对应相等。
基于时差量测箱粒子PHD的被动跟踪多目标方法
技术领域
[0001] 本发明属于制导技术领域，更进一步涉及目标跟踪技术领域中的一种基于时差量测箱粒子概率假设密度PHD(Probability Hypothesis Density)滤波的被动跟踪多目标方法。本发明能够通过无源基站获得的时差量测，对可辐射电磁波的多目标进行实时的被动跟踪。本发明可用于卫星导航或军事目标被动跟踪等领域。
背景技术
[0002] 对目标的定位跟踪常用方法是利用雷达、声呐等有源设备进行，然而有源设备需要发射信号，容易暴露自己的位置。无源(被动)定位跟踪技术利用目标本身辐射的电磁波进行定位，具有隐蔽性好，生存能力强等优点，在目标定位跟踪领域具有重要的意义。对于单个移动辐射源跟踪问题可以视为经典的非线性跟踪问题，其基本框架是基于贝叶斯理论的非线性滤波。在多目标场景中，传统的多目标跟踪算法需要解决数据关联问题，计算复杂，适用性较差。近年来，Mahler提出的基于随机有限集的多目标跟踪方法避免了多目标跟踪中复杂的数据关联问题。Vo在此基础上提出了粒子概率假设密度PHD滤波，有效的解决了基于随机集理论中多目标跟踪算法的应用问题。
[0003] 西安电子科技大学在其申请的专利文献“基于粒子滤波的被动传感器多目标跟踪方法”(专利申请号201010507225.8，公开号CN101975575A)公开了一种基于被动测角粒子滤波多目标跟踪方法，该方法实现的具体步骤是，(1)根据待跟踪目标初始分布抽取目标样本；(2)计算预测联合样本；(3)优化粒子群；(4)对联合样本权值更新及归一化；(5)估计各跟踪目标状态；(6)重采样；(7)重复步骤(2)继续跟踪。该方法能根据测角被动跟踪多目标，但是，该方法仍然存在的不足之处是，跟踪过程需要大量粒子参与运算，跟踪效率较低，在很大程度上影响其跟踪效果。
[0004] 宋骊平在其发表的论文文献“基于箱粒子的多扩展目标PHD滤波”公开了一种箱粒子概率假设密度多扩展目标跟踪算法，该方法实现的具体步骤是，(1)初始化箱粒子；(2)合并新生箱粒子和存活箱粒子；(3)根据状态转移函数产生预测箱粒子；(4)利用CP准则收缩箱粒子更新权值；(5)计算箱粒子权值和得到目标估计数目；(6)重采样；(7)提取目标状态。该方法用较少的箱粒子代替较多的点粒子，在一定程度上解决了粒子滤波运算效率较低的问题，但是，该方法仍存在的不足之处是，由于量测和状态是线性的关系，即目标的状态为坐标，量测同样为坐标，量测和状态为同一物理量且它们的维度相同，在利用CP准则约束箱粒子更新权值时，只需比较两者的重叠部分。然而在时差量测被动跟踪中，量测和状态为高度非线性的关系，即目标的状态坐标，量测为时差值，状态和量测为不同的物理量且维度也不同，无法约束箱粒子更新权值，因此该方法无法适用于被动跟踪量测为时差的多目标。
发明内容
[0005] 本发明的目的是针对上述现有技术的不足，提出一种基于时差量测箱粒子PHD的被动跟踪多目标方法。
[0006] 实现本发明目的的思路是，采用箱粒子进行概率假设密度PHD滤波，将当前跟踪时刻的时差量测扩张成时差箱粒子，然后利用约束传播函数更新预测后的箱粒子，最后获取目标的状态。
[0007] 本发明的具体步骤如下：
[0008] (1)获取起始跟踪时刻多目标状态分布的每个箱粒子：
[0009] (1a)按照下式，获取起始跟踪时刻多目标状态分布的每个点粒子：
[0010]
[0011] 其中， 表示起始跟踪时刻多目标状态分布的第i个点粒子，N表示由跟踪场景的复杂度所确定的随机采样总数，Ψ(·)表示随机采样函数，A表示由起始时刻n0个目标状态组成的矩阵，n0表示起始跟踪时刻多目标的总数，P0表示一个对角线元素为[40,1,40,1]的对角矩阵；
[0012] (1b)利用区间扩张法，将起始跟踪时刻多目标状态分布的每个点粒子扩张成区间大小均为[240,6,240,6]T的箱粒子，其中，T表示转置操作；
[0013] (1c)将1/N的值作为每个箱粒子的权值；
[0014] (2)获取新生箱粒子：
[0015] (2a)按照下式，计算每个新生点粒子：
[0016]
[0017] 其中，φ(c)表示第c个新生点粒子，γ表示由跟踪场景的复杂度所确定的新生箱粒子随机采样总数，B表示由当前跟踪时刻的前一时刻无源基站所获得的所有时差量测代入Chan定位方程得到的矩阵；
[0018] (2b)利用区间扩张法，将每个新生点粒子扩张成区间大小均为[240,6,240,6]T的新生箱粒子；
[0019] (2c)将 的值作为每个新生箱粒子的权值，其中，PB表示根据待跟踪目标所处区域新生目标数目的统计知识设定的新生目标概率；
[0020] (3)合并所有的新生箱粒子与所有的箱粒子，得到合并后的所有的箱粒子；
[0021] (4)对合并后的每个箱粒子进行预测：
[0022] (4a)利用状态转移公式，预测合并后每个箱粒子的状态；
[0023] (4b)利用权值预测公式，预测合并后每个箱粒子的权值；
[0024] (5)利用时差约束传播函数对预测后的每个箱粒子进行更新：
[0025] (5a)按照下式，计算预测后每个箱粒子的时差量测点粒子中每个元素值：
[0026]
[0027] 其中， 表示预测后第y个箱粒子的时差量测点粒子中的第m个元素值，||·||表示取模值操作，p(y)表示预测后的第y个箱粒子的中心点，Sλ表示第λ个无源基站的坐标，m＝λ-1，C表示电磁波在真空中的传播速度，则预测后第y个箱粒子的预测点粒子为M表示无源基站的总数；
[0028] (5b)按照下式，将预测后的每个箱粒子的时差量测点粒子扩张成箱粒子：
[0029] [t(y)]＝t(y)+H×([p(y)]-p(y))
[0030] 其中，[t(y)]表示第y个预测后箱粒子的时差量测箱粒子，H表示t(y)在p(y)处的雅克比矩阵，[p(y)]表示第y个预测后箱粒子；
[0031] (5c)利用区间扩张法，将当前跟踪时刻无源基站获得的每个时差量测扩张成区间为ε的时差量测箱粒子，其中，ε为依据无源基站的性能参数得到的最大测量误差；
[0032] (5d)按照下式，用每个时差量测箱粒子约束每个箱粒子：
[0033]
[0034] 其中， 表示第j个时差量测箱粒子对第y个预测后箱粒子进行约束后得到的箱粒子，H-1表示矩阵H的逆矩阵，∩表示求交集操作，[zj]表示当前跟踪时刻第j个时差量测箱粒子；
[0035] (5e)将约束后箱粒子 中的最大箱粒子作为更新后的第 个箱粒子， 和y的值对应相等，其中，n′表示当前跟踪时刻时差量测箱粒子的总数；
[0036] (5f)按照下式，对每个预测后的箱粒子的权值进行更新，得到更新后的每个箱粒子的权值：
[0037]
[0038] 其中， 表示更新后的第 个箱粒子的权值，PD表示依据无源基站的性能参数设定的检测概率，∑表示求和操作，|·|表示取箱粒子面积的操作，κ表示根据目标所处跟踪区域杂波分布统计设定的杂波强度，α(y)表示预测后第y个箱粒子的权值；
[0039] (6)对更新后的每个箱粒子进行重采样：
[0040] 将更新后的每个箱粒子带入随机子划分重采样方程，得到重采样后的箱粒子；
[0041] (7)获得当前跟踪时刻多目标的状态：
[0042] (7a)对重采样后所有的箱粒子的权值求和，得到当前跟踪时刻多目标的数目估计值；
[0043] (7b)利用Kmeans聚类中的欧式距离准则，对重采样后的所有的箱粒子聚出多目标的状态；
[0044] (8)判断每个无源基站所获得时差量测个数是否均为0，若是，执行步骤(9)，否则，执行步骤(2)；
[0045] (9)结束目标跟踪。
[0046] 本发明与现有技术相比具有以下优点：
[0047] 第一，由于本发明采用箱粒子进行预测和更新，克服了现有技术中需要采用大量点粒子进行运算导致跟踪效率低的问题，使得本发明可以用较少的箱粒子代替现有技术中的较多点粒子，降低了预测和更新的运算复杂度，从而提高了被动跟踪多目标的效率。
[0048] 第二，由于本发明采用时差约束传播函数更新箱粒子，克服了现有箱粒子滤波技术在时差量测下无法更新预测箱粒子，无法适用于时差量测被动跟踪多目标，导致采用粒子滤波被动跟踪技术带来的跟踪效率较低的问题，从而使得本发明能够实现根据时差量测箱粒子快速被动跟踪多目标。
附图说明
[0049] 图1为本发明的流程图；
[0050] 图2为本发明的仿真图。
具体实施方式
[0051] 下面结合附图对本发明做进一步的描述。
[0052] 结合附图1中，对本发明的具体步骤做进一步的描述。
[0053] 步骤1，获取起始跟踪时刻多目标状态分布的每个箱粒子。
[0054] 按照下式，获取起始跟踪时刻多目标状态分布的每个点粒子：
[0055]
[0056] 其中， 表示起始跟踪时刻多目标状态分布的第i个点粒子，N表示由跟踪场景的复杂度所确定的随机采样总数，Ψ(·)表示随机采样函数，A表示由起始时刻n0个目标状态组成的矩阵，n0表示起始跟踪时刻多目标的总数，P0表示一个对角线元素为[40,1,40,1]的对角矩阵；
[0057] 利用区间扩张法，将起始跟踪时刻多目标状态分布的每个点粒子扩张成区间大小均为[240,6,240,6]T的箱粒子，其中，T表示转置操作。
[0058] 所述的区间扩张法是指用每个点粒子减去其对应的区间长度，将差值作为箱粒子的下限，用每个点粒子加上其对应的区间长度，将差值作为箱粒子的上限。
[0059] 对跟踪起始时刻多目标状态分布的每个箱粒子均赋予权值1/N。
[0060] 步骤2，获取新生箱粒子。
[0061] 按照下式，计算每个新生点粒子：
[0062]
[0063] 其中，φ(c)表示第c个新生点粒子，γ表示由跟踪场景的复杂度所确定的新生箱粒子随机采样总数，B表示由当前跟踪时刻的前一时刻无源基站所获得的所有时差量测代入Chan定位方程得到的矩阵。
[0064] 利用区间扩张法，将每个新生点粒子扩张成区间大小均为[240,6,240,6]T的新生箱粒子。
[0065] 所述的区间扩张法是指用每个点粒子减去其对应的区间长度，将差值作为箱粒子的下限，用每个点粒子加上其对应的区间长度，将差值作为箱粒子的上限。
[0066] 将 的值作为每个新生箱粒子的权值，其中，PB表示根据待跟踪目标所处区域新生目标数目的统计知识设定的新生目标概率。
[0067] 步骤3，将所有的新生箱粒子和所有的箱粒子合并，得到合并后的所有的箱粒子。
[0068] 步骤4，对合并后的每个箱粒子进行预测。
[0069] 按照状态转移公式，预测合并后每个箱粒子的状态。
[0070] 所述的状态转移公式为：
[0071] [p(y)]＝f([x(s)],[ω])
[0072] 其中，[p(y)]表示预测后的第y个箱粒子，f(·)表示一个线性函数或者非线性函数，目标的运动方式为匀速直线运动时，f(·)表示一个线性函数，目标的运动方式为除匀速直线运动之外的其他任意运动方式时，f(·)表示一个非线性函数，[x(s)]表示合并后的第s个箱粒子，y和s的值对应相等，[ω]表示当前跟踪时刻的过程噪声。
[0073] 按照权值预测公式，预测合并后每个箱粒子的权值；
[0074] 所述的权值预测公式为：
[0075] α(y)＝PS×w(s)
[0076] 其中，α(y)表示预测后第y个箱粒子的权值，PS表示依据目标的属性设定的存活概率，w(s)表示合并后的第s个箱粒子的权值，y与s的值对应相等。
[0077] 步骤5，利用时差约束传播函数对预测后的每个箱粒子进行更新。
[0078] 按照下式，计算预测后每个箱粒子的时差量测点粒子中每个元素值：
[0079]
[0080] 其中， 表示预测后第y个箱粒子的时差量测点粒子中的第m个元素值，||·||表示取模值操作，p(y)表示预测后的第y个箱粒子的中心点，Sλ表示第λ个无源基站的坐标，m＝λ-1，C表示电磁波在真空中的传播速度，则预测后第y个箱粒子的预测点粒子为M表示无源基站的总数。
[0081] 按照下式，将预测后的每个箱粒子的时差量测点粒子扩张成箱粒子：
[0082] [t(y)]＝t(y)+H×([p(y)]-p(y))
[0083] 其中，[t(y)]表示第y个预测后箱粒子的时差量测箱粒子，H表示t(y)在p(y)处的雅克比矩阵。
[0084] 利用区间扩张法，将当前跟踪时刻无源基站获得的每个时差量测扩张成区间为ε的时差量测箱粒子，其中，ε表示依据无源基站的性能参数得到的最大测量误差。
[0085] 按照下式，每个时差箱粒子对每个箱粒子进行约束：
[0086]
[0087] 其中， 表示第j个时差量测箱粒子对第y个预测后箱粒子进行约束后得到的箱粒子，H-1表示矩阵H的逆矩阵，∩表示求交集操作，[zj]表示当前跟踪时刻第j个时差量测箱粒子。
[0088] 将约束后箱粒子 中的最大箱粒子作为更新后的第 个箱粒子， 和y的值对应相等，其中，n′表示当前跟踪时刻时差量测箱粒子的总数。
[0089] 按照下式，对每个预测后的箱粒子的权值进行更新，得到更新后的每个箱粒子的权值：
[0090]
[0091] 其中， 表示更新后的第 个箱粒子的权值，PD表示依据无源基站的性能参数设定的检测概率，∑表示求和操作，|·|表示取箱粒子面积的操作，κ表示根据目标所处跟踪区域杂波分布统计设定的杂波强度，α(y)表示预测后第y个箱粒子的权值。
[0092] 步骤6，对更新后的每个箱粒子进行重采样。
[0093] 将更新后的每个箱粒子带入随机子划分重采样方程，得到重采样后的箱粒子。
[0094] 步骤7，获得当前跟踪时刻多目标的状态。
[0095] 对重采样后所有的箱粒子的权值求和，得到当前跟踪时刻多目标的数目估计值。
[0096] 利用Kmeans聚类中的欧式距离准则，对重采样后的所有的箱粒子聚出多目标的状态。
[0097] 步骤8，判断每个无源基站所获得时差量测个数是否均为0，若是，执行步骤9，否则，执行步骤2。
[0098] 步骤9，结束目标跟踪。
[0099] 下面结合仿真实验对本发明的效果做进一步说明。
[0100] 1.仿真实验条件：
[0101] 本发明仿真实验的硬件测试平台是：处理器Intel Core i3-7100CPU，主频为3.90GHz，内存8GB；软件平台为：Windows 7旗舰版，64位操作系统，MATLAB R2010a。
[0102] 2.仿真内容：
[0103] 本发明的仿真场景是，在二维跟踪场景被动跟踪数目随时间发生变化的多目标，目标所处跟踪区域杂波的个数服从均值为5的泊松分布。起始跟踪时刻多目标总数为2，多目标状态矩阵 随机采样总数N＝200。在第6个跟踪时刻新出现一个目标，在第26个跟踪时刻消失了一个目标，目标存活概率PS＝0.99，目标新生概率PB＝0.01，所有目标均为匀速直线运动，因此状态转移函数[p(y)]＝f([x(s)],ω)为一个线性函数，可以写成矩阵形式：
[0104] [p(y)]＝F×[x(s)]+G×ω
[0105] 其中，F表示状态转移矩阵，G表示过程噪声矩阵，ω＝[σ,σ]T,σ为一个服从标准正态分布的随机数。
[0106]
[0107] 本发明的仿真实验中的无源基站的总数为4，坐标分别为S1＝[0m,0m]，S2＝[0,-6000m]， 检测概率PD＝0.98，量测噪声服
从均值为0，方差为4ns的高斯分布，按照拉依达准则，设定ε＝12ns。
[0108] 3.仿真结果分析：
[0109] 附图2为本发明的仿真图，其中，图2(a)为本发明仿真实验所跟踪的真实目标轨迹图，图2(a)中用圆表示目标的新生位置，用正方形表示目标的消失位置，圆与正方形之间的连线表示目标的运动轨迹，图2(a)的横坐标表示正东方向，单位为千米，纵坐标表示正北方向，单位为千米。
[0110] 图2(b)为本发明仿真实验对场景中目标总数的估计图，图2(b)中以实线标示的曲线表示将真实目标数目连接起来的曲线，以正方形标示的点表示本发明方法估计出的目标数，以圆表示的点表示采样总数为2000的粒子滤波方法估计出的目标数，图2(b)的横坐标表示跟踪时刻，单位为秒，纵坐标表示目标个数，单位为个。
[0111] 图2(c)为本发明仿真实验目标位置误差结果图，图2(c)中以正方形标示的连线表示本发明方法得到的以OSPA距离为准则的误差连线，以圆标示的连线表示采样总数为2000的粒子滤波方法得到的以OSPA距离为准则的误差连线，图2(c)的横坐标表示当前跟踪时刻，单位为秒，纵坐标表示OSPA距离，单位为米，在第6跟踪时刻误差出现尖峰的原因是，新生箱粒子的获取是根据当前跟踪时刻的前一时刻无源基站的时差量测得到的，仅会在目标新生的时刻出现漏估，因此OSPA距离会出现一个尖峰。
[0112] 由图2(b)和图2(c)可以看出，本发明方法采用200个箱粒子所得到的跟踪效果与现有技术采用2000个粒子进行粒子滤波所得到的跟踪效果大抵相当。综上所述，在实现跟踪精度相同的条件下，本发明可以降低被动跟踪多目标的计算复杂度，提高了跟踪效率。
