CN104237879A - 一种雷达系统中的多目标跟踪方法 - Google Patents

Abstract

本发明公开了一种雷达系统中的多目标跟踪方法。它是在雷达测量值分扇区处理的前提下，通过对雷达不同扇区的高斯成份进行相应的时间预测、测量修正、融合修剪、目标状态提取，实现跨扇区多目标跟踪，最终实现非同步测量条件下雷达全平面多目标跟踪，从而解决了传统概率假设密度滤波无法直接应用于雷达系统的问题。整个过程不要求所有测量值在同一时刻测量，而且能够处理变化数目的多目标跟踪问题，并且与传统的方法(如联合概率数据关联、多假设跟踪等)相比计算负担小，另外航迹起始、维持与终结都是自然完成的，不需要单独列出。

Description

一种雷达系统中的多目标跟踪方法

技术领域

本发明涉及多目标雷达跟踪技术领域，更具体的，涉及一种雷达系统中的多目标跟踪方法。

背景技术

经过数十年的研究和发展，多目标跟踪技术在理论和应用上都取得了长足进展，已被广泛应用于军事、航天、环境检测等众多领域。

传统多目标跟踪方法通过引入数据关联技术将多目标跟踪问题简化为单目标跟踪问题。常见的数据关联技术主要有：最近邻法(NN)、联合概率数据关联法(JPDA)，以及多假设跟踪法(MHT)等。NN法计算简单、易于实现，但是跟踪性能较差，如在密集杂波、目标交叉等情况下容易跟丢或跟错目标。标准的JPDA法仅能针对固定数目的目标进行跟踪，而实际战场上随时都可能有目标出现或者消失，对于这种目标数目变化的情况，JPDA方法应用起来十分不便。MHT法随着传感器扫描次数和测量回波数目的增加容易出现组合爆炸现象，严重影响算法的实时性。由此可见数据关联的引入却又使得多目标跟踪问题更加复杂化。这使得寻求不需要数据关联的多目标跟踪算法成为一个新的研究课题。

近年来，有限集统计(Finite-Set Statistics,FISST)理论[Mahler R P S.Multi-target Bayesfiltering via frst-order multitarget moments.IEEE Transactions on Aerospace and Electronic Systems,2003,39(4):1152-1178；Mahler R P.Statistical Multisource-Multi-target Information Fusion.Norwood:Artech House,2007]的提出，其作为随机有限集(Random Finite Set,RFS)的重要研究方向，为解决多目标跟踪技术中上述相关难题提供了新的解决渠道，即概率假设密度(PHD)滤波。与传统方法相比，基于FISST的多目标跟踪方法具有坚实的数学理论基础。算法实现过程中不使用数据关联，能够处理变化数目的多目标跟踪问题，并且计算复杂度与传统的方法相比要小，另外航迹起始、维持与终结都是自然完成的，不需要单独列出。不过在实际雷达系统中，由于雷达测量的非同步性，现有的PHD滤波无法直接应用。

实际雷达系统中，由于波速扫描时间不同，在同一扇区内的的测量值也会存在时间上的非同步性。以测量周期12秒，雷达全平面分为32个扇区为例，每个扇区中测量值在时间上的最大差值大约在0.3333秒，对于280米/秒速度的目标而言，这样的测量时差会导致将近100米的模型误差。不同的测量扇区的时间差甚至可达数秒，而现有的PHD滤波算法都是针对相同时刻的测量值进行测量修正，这使得现有PHD滤波无法直接应用于实际雷达系统。

现有的概率假设密度(PHD)滤波只能针对相同时刻测量值才能实现多目标跟踪功能。

发明内容

发明目的：本发明所要解决的技术问题是针对现有技术的不足，提供一种雷达系统中的多目标跟踪方法。

为了解决上述技术问题，本发明公开了一种雷达系统中的多目标跟踪方法，从雷达正北方向开始，以k-1表示雷达系统前一个扫描周期，以k表示雷达系统当前扫描周期，记雷达系统的扫描周期为将雷达系统的扫描界面平分为N个扇区，将测量值根据方位分别放入对应扇区,然后执行以下步骤：

步骤1，接收到当前扫描周期第n扇区的测量值集合n＝1,2,…,N，N取值为自然数；

步骤2，根据当前第n个扇区前一扫描周期对应时刻目标状态强度的高斯项，预测当前时刻对应扇区目标状态强度的高斯项；根据与当前扇区的两侧相邻的两个扇区的前一扫描周期对应时刻目标状态强度的高斯项，预测当前扇区对应时刻目标状态强度的高斯项；

步骤3，根据步骤2所得到的当前扇区及其相邻扇区对应时刻目标状态强度的高斯项，计算当前扇区及其两侧相邻扇区对应时刻的预测测量值及相关预测误差协方差阵；

步骤4，基于步骤3所得到的当前扇区及其两侧相邻两个扇区预测测量值及相关预测误差协方差阵，结合当前扇区的测量值进行跨扇区测量更新；

步骤5，针对步骤4所得到的测量更新后的高斯项进行裁剪与合并，裁剪与合并后的高斯项作为当前时刻的高斯项，根据当前时刻高斯项所在位置将其放入相应的扇区缓冲区中，作为下次相应扇区滤波器递归的输入；

步骤6，根据裁剪与合并后的高斯项，提取权重大于0.5的高斯项作为滤波器的输出，相应高斯项中的均值和方差分别为存活目标的状态估计和误差估计。

本发明中，雷达系统中目标状态表示方程为：

x_k＝F_kx_k-1+w_k-1   (1)，

其中表示k时刻目标状态向量，x_k-1表示k-1时刻目标状态向量，其中x_k,y_k对应表示目标在x轴和y轴的位置、对应表示目标在x轴和y轴的速度，系统噪声w_k服从标准正态分布N(0,Q_k)，状态转移矩阵F_k及噪声方差阵Q_k分别为：

$F_k=\left(\begin{array}{cc}\tilde{F}& 0\\ 0& \tilde{F}\end{array}\right),\tilde{F}=\left(\begin{array}{cc}1& t\\ 0& 1\end{array}\right)$

$Q_k=\left(\begin{array}{cc}\tilde{Q}& 0\\ 0& \tilde{Q}\end{array}\right),\tilde{Q}={\mathrm{\σ}}_w^2\left(\begin{array}{cc}\frac{t^4}{4}& \frac{t^3}{2}\\ \frac{t^3}{2}& t^2\end{array}\right)---\left(2\right),$

其中t为每个扇区对应的时间间隔，t的上标表示t的幂次方，为过程噪声方差值，的取值范围为区间(0,100)；

设雷达的位置为S⁰＝{x⁰,y⁰}，x⁰,y⁰分别为雷达位置的沿x轴和y轴坐标，测量方程为： $z_k=h(x_k,v_k)=\left\{\begin{array}{c}\sqrt{{(x_k-x^0)}^2+{(y_k-y^0)}^2}+v_k^1\\ \arctan (x_k-x^0,y_k-y^0)+v_k^2\end{array}\right.---\left(3\right),$

其中z_k为k时刻的测量值，h(x_k,v_k)表示雷达测量方程，v_k是k时刻雷达测量噪声向量，径向距离测量噪声服从标准正态分布其中为径向距离误差测量方差，取值范围是(0,1000000)，方位角测量噪声服从标准正态分布其中为方位误差测量方差，取值范围是(0,4)。

本发明中，所述雷达系统初始化包括：获得目标初始时刻的状态，包含目标初始时刻的位置和速度，得到第0个周期的目标状态估计集合的高斯项为其中J₀表示式中高斯项的总数，i为自然数，其取值范围为{1,2….,J₀}，代表第i个高斯项对应的权重，代表第i个高斯项对应的均值，代表第i个高斯项对应的协方差矩阵。

本发明中，步骤2中，记k-1周期第n个扇区的高斯成份集合表示为其中表示式中高斯项的总数，j为自然数，其取值范围为 代表第j个高斯项对应的权重，代表第j个高斯项对应的均值，代表第j个高斯项对应的协方差矩阵；利用雷达系统方程将k-1周期第n个扇区及其两侧相邻两个扇区的高斯成份集合进行预测，得到k周期第n个扇区所对应的预测高斯成份集合其中表示预测高斯成份集合高斯项的总数，l为自然数，其取值范围为下标k|k-1表示预测代表第l个预测高斯项对应的权重，代表第l个预测高斯项对应的均值，代表第l个预测高斯项对应的协方差矩阵；该预测高斯成份集合中包括第n个扇区及其相邻两个扇区的存活高斯成份预测项 $w_{k|k-1}^{\left(l\right)}=P_s{w}_{k-1}^{\left(l\right)},m_{k|k-1}^{\left(l\right)}=F_{k-1}{m}_{k-1}^{\left(l\right)},P_{k|k-1}^{\left(l\right)}=Q_k+F_{k-1}{P}_{k-1}^{\left(l\right)}{F}_{k-1}^T,$ 其中P_s是目标存活概率，取值范围为区间(0,1)，F_k-1为k-1时刻状态转移矩阵，上标T为转置。

本发明中，步骤3中，针对k周期第n个扇区所对应的预测高斯成份集合 ${\{w_{k|k-1}^{\left(l\right)},m_{k|k-1}^{\left(l\right)},P_{k|k-1}^{\left(l\right)},n\}}_{l=1}^{J_{k|k-1}^n},$ 分别计算：

预测测量值 ${\mathrm{\η}}_{k|k-1}^{\left(l\right)}=h(m_{k|k-1}^{\left(l\right)},0),$

预测测量误差协方差阵 $S_k^{\left(l\right)}=R_k+H_k^{\left(l\right)}{P}_{k|k-1}^{\left(l\right)}{\left[H_k^{\left(l\right)}\right]}^T,$

增益矩阵 $K_k^{\left(l\right)}=P_{k|k-1}^{\left(l\right)}{\left[H_k^{\left(l\right)}\right]}^T{\left(S_k^{\left(l\right)}\right)}^{-1},$

估计误差协方差阵 $P_{k|k}^{\left(l\right)}=(I-K_k^{\left(l\right)}{H}_k^{\left(l\right)})P_{k|k-1}^{\left(l\right)},$

其中为k时刻雷达系统测量误差协方差阵，I为对应维数的单位矩阵，为测量方程h(x_k,0)在处线性化所得到的矩阵，即：

$H_k^{\left(l\right)}=\frac{\∂h(x_k,0)}{\∂x_k}{|}_{x_k=m_{k|k-1}^{\left(l\right)}}---\left(4\right).$

本发明中，步骤4中，针对第n个扇区的每一个测量值对预测高斯成份进行测量修正得到多目标状态集合矩计算公式为：

$V_k\left(X\right|Z_{1:k}^n)=\underset{l=1}{\overset{J_{k|k-1}^n}{\mathrm{\Σ}}}{w}_{k|k}^{\left(l\right)}\·N(x,m_{k|k-1}^{\left(l\right)},P_{k|k-1}^{\left(l\right)})+\underset{z_j\∈Z_k^n}{\mathrm{\Σ}}\underset{l=1}{\overset{J_{k|k-1}^n}{\mathrm{\Σ}}}{w}_{k|k}^{(l,j)}\·N(x,m_{k|k}^{(l,j)},P_{k|k}^{(l,j)})---\left(5\right),$

其中：

表示均值为协方差阵为的标准正态分布；

表示均值为协方差阵为的标准正态分布；

为直接预测所对应的高斯项权重，为第l个预测高斯项对应的检测概率，取值范围为区间[0.4,1)；

$w_{k|k}^{(l,j)}=\frac{P_D^{\left(l\right)}{w}_{k|k-1}^{\left(l\right)}N(z_j;{\mathrm{\η}}_{k|k-1}^{\left(l\right)},S_k^{\left(l\right)})}{K_k\left(z_j\right)+\underset{q=1}{\overset{J_{k|k-1}^n}{\mathrm{\Σ}}}{P}_D^{\left(q\right)}{w}_{k|k-1}^{\left(q\right)}\·N(z_j;{\mathrm{\η}}_{k|k-1}^{\left(q\right)},S_k^{\left(q\right)})},$

为第l个高斯项和第j个测量值对应的测量修正高斯项权重，式中求和公式中的上标q为自然数，其取值范围为 表示均值为协方差阵为的标准正态分布，表示均值为协方差阵为的标准正态分布，K_k(z_j)为杂波强度函数；

为经测量值z_j修正后的高斯项对应的均值；

为经测量值z_j修正后的高斯项对应的协方差阵，

得到测量修正后的高斯成份集合其中为测量修正后所得高斯项的总数，f为自然数，其取值范围为上标f表示对应第f个测量修正后所得高斯项。

本发明中，步骤5中，删除测量修正后的高斯成份集合中权重的高斯项，其中为裁减门限，取然后将距离充分小的高斯成份合并成一个，即将的高斯项合并成一个，其中为合并门限，其取值范围为区间[4,6]。距离d_i,j定义为： $d_{i,j}={(m_{k|k}^{\left(i\right)}-m_{k|k}^{\left(j\right)})}^T{(P_{k|k}^{\left(i\right)}+P_{k|k}^{\left(j\right)})}^{-1}(m_{k|k}^{\left(i\right)}-m_{k|k}^{\left(j\right)}),$

高斯成份合并方法如下：

合并后对应高斯项的权重

合并后对应高斯项的均值 ${\tilde{m}}_{k|k}^{\left(l\right)}=\frac{1}{{\tilde{w}}_{k|k}^{\left(l\right)}}\underset{i\∈L}{\mathrm{\Σ}}{w}_{k|k}^{\left(i\right)}{m}_{k|k}^{\left(i\right)},$

合并后对应高斯项的协方差阵 ${\tilde{P}}_{k|k}^{\left(l\right)}=\frac{1}{{\tilde{w}}_{k|k}^{\left(l\right)}}\underset{i\∈L}{\mathrm{\Σ}}{w}_{k|k}^{\left(i\right)}(P_{k|k}^{\left(i\right)}+({\tilde{m}}_{k|k}^{\left(l\right)}-m_{k|k}^{\left(i\right)}){({\tilde{m}}_{k|k}^{\left(l\right)}-m_{k|k}^{\left(i\right)})}^T),$

其中L为所有可合并高斯成份的指标集合；将融合修剪后的高斯成份根据其方位分别放入对应扇区及其相邻的两个扇区，于是得到k周期第n个扇区高斯成份集合： ${\{w_k^{\left(i\right)},m_k^{\left(i\right)},P_n^{\left(i\right)},n\}}_{i=1}^{J_k^n}.$

本发明中，步骤6中，将权重大于0.5的高斯成份均值提取出来作为第n个扇区前一扇区目标状态估计值。

本发明中，当k周期第n个扇区后一扇区的测量值集合发送到数据处理器以后，重复步骤2，实现雷达系统中跨扇区全平面目标状态估计。

本发明的核心原理是：1)预测当前对应扇区目标状态的高斯项，使用当前和相邻扇区上一扫描周期对应时刻目标状态强度高斯项的并集；2)测量更新时，使用当前和相邻扇区高斯项的并集的预测测量值；3)延迟一个扇区，处理测量更新后的高斯项进行裁剪与合并，并且使用该扇区和其相邻扇区更新后的高斯项合并进行裁剪与合并，处理完毕再根据其位置分别放入相应扇区；4)滤波器的输出延迟两个扇区。

本发明通过分扇区处理雷达测量数据，解决雷达系统中测量扇区时间非同步问题；通过跨扇区处理相邻扇区目标状态强度的高斯成分，解决雷达全平面目标跨扇区连续跟踪问题。

本发明为了使得PHD滤波可以针对雷达系统不同扇区测量时间作相应的时间预测和测量更新，本发明要基于有限集统计(Finite-Set Statistics,FISST)理论，重新设计PHD滤波中的时间预测步和测量修正步，从而实现雷达系统中分扇区测量条件下的多目标跟踪功能。

有益效果：本发明实现过程中不使用数据关联，能够处理变化数目的多目标跟踪问题，并且与传统的方法(如联合概率数据关联、多假设跟踪等)相比计算负担小，另外航迹起始、维持与终结都是自然完成的，不需要单独列出。值得一提的是，本发明可以直接应用于分扇区雷达系统中测量时间非同步条件下的多目标跟踪，而现有的PHD滤波却无法实现该功能。

附图说明

下面结合附图和具体实施方式对本发明做更进一步的具体说明，本发明的上述和/或其他方面的优点将会变得更加清楚。

图1是本发明的多目标跟踪方法原理框图。

图2是雷达显示界面示意图。

图3是本发明实施例目标编队飞行的轨迹图。

图4是本发明方法确定的目标运动轨迹。

具体实施方式

雷达系统中基于FISST的测量扇区时间不同步多目标跟踪方法。具体实施步骤是：假设从雷达正北方向开始，将雷达扫描界面平分为N(例如N＝36)个扇区，如图2所示，将测量值根据方位分别放入不同扇区；

步骤1：接收到当前第n(n＝1,2,…,N)个扇区的测量值集合；

步骤2：根据当前第n(n＝1,2,…,N)个扇区前一扫描周期对应时刻目标状态强度的高斯项，预测当前时刻对应扇区目标状态强度的高斯项；根据与当前扇区相邻的两个扇区的前一扫描周期对应时刻目标状态强度的高斯项，预测当前扇区对应时刻目标状态强度的高斯项；

步骤3：根据步骤2所得到的当前扇区(及其两侧相邻扇区)对应时刻目标状态强度的高斯项，计算当前扇区(及其两侧相邻扇区)对应时刻的预测测量值及相关预测误差协方差阵等；

步骤4：基于步骤3所得到的(当前扇区、及其两侧相邻两个扇区)预测测量值及相关预测误差协方差阵，结合当前扇区的测量值进行跨扇区测量更新；

步骤5：针对步骤4所得到的测量更新后的高斯项进行裁减与合并，裁减与合并后的高斯项作为当前时刻的高斯项，并根据其位置分别放入相应的扇区，当前时刻的高斯成份作为下一次相应扇区滤波器递归的输入；

步骤6：根据裁减与合并后的高斯项，提取权重大于0.5的高斯项作为滤波器的输出，相应高斯项中的均值和方差分别为存活目标的状态估计和误差估计。

本发明的雷达系统多目标跟踪方法能够作为滤波器输出的高斯项，在滤波器预测、更新和裁剪合并时，滤波器输出目标状态。

本实施例中雷达系统中目标状态方程为：

x_k＝F_kx_k-1+w_k-1   (1)，

其中表示k时刻目标状态向量，x_k-1表示k-1时刻目标状态向量，其中x_k,y_k对应表示目标在x轴和y轴的位置、对应表示目标在x轴和y轴的速度，系统噪声w_k服从标准正态分布N(0,Q_k)，状态转移矩阵F_k及噪声方差阵Q_k分别为：

$F_k=\left(\begin{array}{cc}\tilde{F}& 0\\ 0& \tilde{F}\end{array}\right),\tilde{F}=\left(\begin{array}{cc}1& t\\ 0& 1\end{array}\right)$

$Q_k=\left(\begin{array}{cc}\tilde{Q}& 0\\ 0& \tilde{Q}\end{array}\right),\tilde{Q}={\mathrm{\σ}}_w^2\left(\begin{array}{cc}\frac{t^4}{4}& \frac{t^3}{2}\\ \frac{t^3}{2}& t^2\end{array}\right)---\left(2\right),$

其中t为每个扇区对应的时间间隔，t的上标表示t的幂次方，为过程噪声方差值，的取值范围为区间(0,100)；

设雷达的位置为S⁰＝{x⁰,y⁰}，x⁰,y⁰分别为雷达位置的沿x轴和y轴坐标，测量方程为：

$z_k=h(x_k,v_k)=\left\{\begin{array}{c}\sqrt{{(x_k-x^0)}^2+{(y_k-y^0)}^2}+v_k^1\\ \arctan (x_k-x^0,y_k-y^0)+v_k^2\end{array}\right.---\left(3\right),$

其中z_k为k时刻的测量值，h(x_k,v_k)表示雷达测量方程，v_k是k时刻雷达测量噪声向量，径向距离测量噪声服从标准正态分布其中为径向距离误差测量方差，取值范围是(0,1000000)，方位角测量噪声服从标准正态分布其中为方位误差测量方差，取值范围是(0,4)。

如图1所示，本发明主要包括：初始化模块、测量值分区模块、高斯成份预测模块、预测测量值构建模块、高斯成份测量修正模块、高斯成份融合修剪模块、目标状态提取模块。结合流程图说明具体实现步骤为：

步骤1：初始化，通过传感器获得目标初始时刻的状态，包含目标的位置和速度信息,得到第0个周期的目标状态估计集合的高斯项为

步骤2：测量值分区，以k-1表示雷达系统前一扫描周期，以k表示当前雷达扫描周期,记雷达扫描周期为雷达全平面被平分为N个扇区，将测量值根据方位不同分置于不同扇区，当前第n(n＝1,2,…,N)个扇区的测量值集合记为

步骤3：高斯成份预测，记k-1周期第n(n＝1,2,…,N)个扇区的高斯成份集合表示为利用系统方程将k-1周期第n(n＝1,2,…,N)个扇区及其相邻两个扇区的高斯成份集合进行预测，得到k周期第n(n＝1,2,…,N)个扇区所对应的预测高斯成份集合其中表示预测高斯成份集合高斯项的总数，l为自然数，其取值范围为下标k|k-1表示预测代表第l个预测高斯项对应的权重，代表第l个预测高斯项对应的均值，代表第l个预测高斯项对应的协方差矩阵；该预测高斯成份集合中包括第n个扇区及其相邻两个扇区的存活高斯成份预测项 $w_{k|k-1}^{\left(l\right)}=P_s{w}_{k-1}^{\left(l\right)},m_{k|k-1}^{\left(l\right)}=F_{k-1}{m}_{k-1}^{\left(l\right)},P_{k|k-1}^{\left(l\right)}=Q_k+F_{k-1}{P}_{k-1}^{\left(l\right)}{F}_{k-1}^T,$ 需要指出的是：不同扇区所采用的状态转移矩阵中的时间间隔t并不相同。对应地，分别取t₁，t₂，t₃。第n个扇区高斯成份预测时，取第n个扇区前一扇区高斯成份预测时，取第n个扇区后一扇区高斯成份预测时，取)，以及新出生高斯成份、衍生高斯成份。

步骤4：预测测量值构建，针对步骤3所得到的k周期第n个扇区所对应的预测高斯成份集合 ${\{w_{k|k-1}^{\left(l\right)},m_{k|k-1}^{\left(l\right)},P_{k|k-1}^{\left(l\right)},n\}}_{l=1}^{J_{k|k-1}^n},$ 分别计算：

预测测量值 ${\mathrm{\η}}_{k|k-1}^{\left(l\right)}=h(m_{k|k-1}^{\left(l\right)},0),$

预测测量误差协方差阵 $S_k^{\left(l\right)}=R_k+H_k^{\left(l\right)}{P}_{k|k-1}^{\left(l\right)}{\left[H_k^{\left(l\right)}\right]}^T,$

增益矩阵 $K_k^{\left(l\right)}=P_{k|k-1}^{\left(l\right)}{\left[H_k^{\left(l\right)}\right]}^T{\left(S_k^{\left(l\right)}\right)}^{-1},$

估计误差协方差阵 $P_{k|k}^{\left(l\right)}=(I-K_k^{\left(l\right)}{H}_k^{\left(l\right)})P_{k|k-1}^{\left(l\right)},$

其中为k时刻雷达系统测量误差协方差阵，I为对应维数的单位矩阵，为测量方程h(x_k,0)在处线性化所得到的矩阵，即：

$H_k^{\left(l\right)}=\frac{\∂h(x_k,0)}{\∂x_k}{|}_{x_k=m_{k|k-1}^{\left(l\right)}}---\left(4\right)$ 。

步骤5：高斯成份测量修正，针对第n个扇区的每一个测量值对预测高斯成份进行测量修正得到多目标状态集合矩计算公式为：

$V_k\left(X\right|Z_{1:k}^n)=\underset{l=1}{\overset{J_{k|k-1}^n}{\mathrm{\Σ}}}{w}_{k|k}^{\left(l\right)}\·N(x,m_{k|k-1}^{\left(l\right)},P_{k|k-1}^{\left(l\right)})+\underset{z_j\∈Z_k^n}{\mathrm{\Σ}}\underset{l=1}{\overset{J_{k|k-1}^n}{\mathrm{\Σ}}}{w}_{k|k}^{(l,j)}\·N(x,m_{k|k}^{(l,j)},P_{k|k}^{(l,j)})---\left(5\right),$

其中：

表示均值为协方差阵为的标准正态分布；

表示均值为协方差阵为的标准正态分布；

为直接预测所对应的高斯项权重，为第l个预测高斯项对应的检测概率，取值范围为区间[0.4,1)；

$w_{k|k}^{(l,j)}=\frac{P_D^{\left(l\right)}{w}_{k|k-1}^{\left(l\right)}N(z_j;{\mathrm{\η}}_{k|k-1}^{\left(l\right)},S_k^{\left(l\right)})}{K_k\left(z_j\right)+\underset{q=1}{\overset{J_{k|k-1}^n}{\mathrm{\Σ}}}{P}_D^{\left(q\right)}{w}_{k|k-1}^{\left(q\right)}\·N(z_j;{\mathrm{\η}}_{k|k-1}^{\left(q\right)},S_k^{\left(q\right)})},$

为第l个高斯项和第j个测量值对应的测量修正高斯项权重，式中求和公式中的上标q为自然数，其取值范围为 表示均值为协方差阵为的标准正态分布，表示均值为协方差阵为的标准正态分布，K_k(z_j)为杂波强度函数；

为经测量值z_j修正后的高斯项对应的均值；

为经测量值z_j修正后的高斯项对应的协方差阵，

得到测量修正后的高斯成份集合其中为测量修正后所得高斯项的总数，f为自然数，其取值范围为上标f表示对应第f个测量修正后所得高斯项。

步骤6：高斯成份融合修剪，删除测量修正后的高斯成份集合中权重的高斯项，其中为裁减门限，取然后将距离充分小的高斯成份合并成一个，即将的高斯项合并成一个，其中为合并门限，其取值范围为区间[4,6]。距离d_i,j定义为： $d_{i,j}={(m_{k|k}^{\left(i\right)}-m_{k|k}^{\left(j\right)})}^T{(P_{k|k}^{\left(i\right)}+P_{k|k}^{\left(j\right)})}^{-1}(m_{k|k}^{\left(i\right)}-m_{k|k}^{\left(j\right)}),$

高斯成份合并方法如下：

合并后对应高斯项的权重

合并后对应高斯项的均值 ${\tilde{m}}_{k|k}^{\left(l\right)}=\frac{1}{{\tilde{w}}_{k|k}^{\left(l\right)}}\underset{i\∈L}{\mathrm{\Σ}}{w}_{k|k}^{\left(i\right)}{m}_{k|k}^{\left(i\right)},$

合并后对应高斯项的协方差阵 ${\tilde{P}}_{k|k}^{\left(l\right)}=\frac{1}{{\tilde{w}}_{k|k}^{\left(l\right)}}\underset{i\∈L}{\mathrm{\Σ}}{w}_{k|k}^{\left(i\right)}(P_{k|k}^{\left(i\right)}+({\tilde{m}}_{k|k}^{\left(l\right)}-m_{k|k}^{\left(i\right)}){({\tilde{m}}_{k|k}^{\left(l\right)}-m_{k|k}^{\left(i\right)})}^T),$

其中L为所有可合并高斯成份的指标集合；将融合修剪后的高斯成份根据其方位分别放入对应扇区及其相邻的两个扇区，于是得到k周期第n个扇区高斯成份集合： ${\{w_k^{\left(i\right)},m_k^{\left(i\right)},P_n^{\left(i\right)},n\}}_{i=1}^{J_k^n}.$

步骤7：目标状态提取，完成步骤六以后，k周期第n个扇区前一扇区中的高斯成份将不再发生变化，故可以将权重大于0.5的高斯成份均值提取出来作为第n个扇区前一扇区目标状态估计值。

进一步地，当k周期第n个扇区后一扇区的测量值集合发送到数据处理器以后，重复步骤3，于是可实现雷达系统中跨扇区全平面目标状态估计。

本实施例在既有新目标出现和已有目标消失的情况下，对图3所示的仿真数据进行处理，本实施例得到的跟踪结果如图4所示。从图4中可以看出，本实施例提出的跟踪方法能探测出观测空间的全部6批目标并能进行有效跟踪。

和现有的多目标跟踪方法相比，本实施例特点是：通过对雷达不同扇区的高斯成份进行相应的时间预测、测量修正、融合修剪，实现跨扇区多目标跟踪，最终实现非同步测量条件下雷达全平面多目标跟踪，从而解决了传统概率假设密度滤波无法直接应用于雷达系统的问题，提高了多目标跟踪方法的工程使用价值。

本申请的技术内容，还承担了基金项目：中国博士后科学基金资助项目(2013M541643)，国家自然科学基金资助项目(61403352)。

本发明提供了一种雷达系统中的多目标跟踪方法，具体实现该技术方案的方法和途径很多，以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。本实施例中未明确的各组成部分均可用现有技术加以实现。

Claims (9)

1.一种雷达系统中的多目标跟踪方法，其特征是，从雷达正北方向开始，以k-1表示雷达系统前一个扫描周期，以k表示雷达系统当前扫描周期，记雷达系统的扫描周期为将雷达系统的扫描界面平分为N个扇区，将测量值根据方位分别放入对应扇区,然后执行以下步骤：

步骤1，接收到当前扫描周期第n扇区的测量值集合n＝1,2,…,N，N取值为自然数；

步骤2，根据当前第n个扇区前一扫描周期对应时刻目标状态强度的高斯项，预测当前时刻对应扇区目标状态强度的高斯项；根据与当前扇区的两侧相邻的两个扇区的前一扫描周期对应时刻目标状态强度的高斯项，预测当前扇区对应时刻目标状态强度的高斯项；

步骤3，根据步骤2所得到的当前扇区及其相邻扇区对应时刻目标状态强度的高斯项，计算当前扇区及其两侧相邻扇区对应时刻的预测测量值及相关预测误差协方差阵；

步骤4，基于步骤3所得到的当前扇区及其两侧相邻两个扇区预测测量值及相关预测误差协方差阵，结合当前扇区的测量值进行跨扇区测量更新；

步骤5，针对步骤4所得到的测量更新后的高斯项进行裁剪与合并，裁剪与合并后的高斯项作为当前时刻的高斯项，根据当前时刻高斯项所在位置将其放入相应的扇区缓冲区中，作为下次相应扇区滤波器递归的输入；

步骤6，根据裁剪与合并后的高斯项，提取权重大于0.5的高斯项作为滤波器的输出，相应高斯项中的均值和方差分别为存活目标的状态估计和误差估计。

2.根据权利要求1所述的雷达系统的多目标跟踪方法，其特征在于：所述雷达系统中目标状态方程为：

x_k＝F_kx_k-1+w_k-1   (1)，

其中表示k时刻目标状态向量，x_k-1表示k-1时刻目标状态向量，其中x_k,y_k对应表示目标在x轴和y轴的位置、对应表示目标在x轴和y轴的速度，系统噪声w_k服从标准正态分布N(0,Q_k)，状态转移矩阵F_k及噪声方差阵Q_k分别为：

$F_k=\left(\begin{array}{cc}\tilde{F}& 0\\ 0& \tilde{F}\end{array}\right),\tilde{F}=\left(\begin{array}{cc}1& t\\ 0& 1\end{array}\right)$

$Q_k=\left(\begin{array}{cc}\tilde{Q}& 0\\ 0& \tilde{Q}\end{array}\right),\tilde{Q}={\mathrm{\σ}}_w^2\left(\begin{array}{cc}\frac{t^4}{4}& \frac{t^3}{2}\\ \frac{t^3}{2}& t^2\end{array}\right)---\left(2\right),$

其中t为每个扇区对应的时间间隔，t的上标表示t的幂次方，为过程噪声方差值，的取值范围为区间(0,100)；

设雷达的位置为S⁰＝{x⁰,y⁰}，x⁰,y⁰分别为雷达位置沿x轴和y轴的坐标，测量方程为：

$z_k=h(x_k,v_k)=\left\{\begin{array}{c}\sqrt{{(x_k-x^0)}^2+{(y_k-y^0)}^2}+v_k^1\\ \arctan (x_k-x^0,y_k-y^0)+v_k^2\end{array}\right.---\left(3\right),$

其中z_k为k时刻的测量值，h(x_k,v_k)表示雷达测量方程，v_k是k时刻雷达测量噪声向量，径向距离测量噪声服从标准正态分布其中为径向距离误差测量方差，取值范围是(0,1000000)，方位角测量噪声服从标准正态分布其中为方位误差测量方差，取值范围是(0,4)。

3.根据权利要求2所述的雷达系统的多目标跟踪方法，其特征在于：所述雷达系统初始化包括：获得目标初始时刻的状态，包含目标初始时刻的位置和速度，得到第0个周期的目标状态估计集合的高斯项为其中J₀表示式中高斯项的总数，i为自然数，取值范围为1,2….,J₀，代表第i个高斯项对应的权重，代表第i个高斯项对应的均值，代表第i个高斯项对应的协方差矩阵。

4.根据权利要求3所述的雷达系统的多目标跟踪方法，其特征在于：步骤2中，记k-1周期第n个扇区的高斯成份集合表示为其中表示式中高斯项的总数，j为自然数，取值范围为1,2...., 代表第j个高斯项对应的权重，代表第j个高斯项对应的均值，代表第j个高斯项对应的协方差矩阵；利用雷达系统方程将k-1周期第n个扇区及其两侧相邻两个扇区的高斯成份集合进行预测，得到k周期第n个扇区所对应的预测高斯成份集合其中表示预测高斯成份集合高斯项的总数，l为自然数，取值范围为1,2....,下标k|k-1表示预测，代表第l个预测高斯项对应的权重，代表第l个预测高斯项对应的均值，代表第l个预测高斯项对应的协方差矩阵；该预测高斯成份集合中包括第n个扇区及其相邻两个扇区的存活高斯成份预测项 其中P_s是目标存活概率，取值范围为区间(0,1)，F_k-1为k-1时刻状态转移矩阵，上标T为转置。

5.根据权利要求4所述的雷达系统的多目标跟踪方法，其特征在于：步骤3中，针对k周期第n个扇区所对应的预测高斯成份集合分别计算：

预测测量值 ${\mathrm{\η}}_{k|k-1}^{\left(l\right)}=h(m_{k|k-1}^{\left(l\right)},0),$

预测测量误差协方差阵 $S_k^{\left(l\right)}=R_k+H_k^{\left(l\right)}{P}_{k|k-1}^{\left(l\right)}{\left[H_k^{\left(l\right)}\right]}^T,$

增益矩阵 $K_k^{\left(l\right)}=P_{k|k-1}^{\left(l\right)}{\left[H_k^{\left(l\right)}\right]}^T{\left(S_k^{\left(l\right)}\right)}^{-1},$

估计误差协方差阵 $P_{k|k}^{\left(l\right)}=(I-K_k^{\left(l\right)}{H}_k^{\left(l\right)})P_{k|k-1}^{\left(l\right)},$

其中为k时刻雷达系统测量误差协方差阵，I为对应维数的单位矩阵，为测量方程h(x_k,0)在处线性化所得到的矩阵，即：

$H_k^{\left(l\right)}=\frac{\∂h(x_k,0)}{\∂x_k}{|}_{x_k=m_{k|k-1}^{\left(l\right)}}---\left(4\right).$

6.根据权利要求5所述的雷达系统的多目标跟踪方法，其特征在于：步骤4中，针对第n个扇区的每一个测量值对预测高斯成份进行测量修正得到多目标状态集合矩计算公式为：

$V_k\left(X\right|Z_{1:k}^n)=\underset{l=1}{\overset{J_{k|k-1}^n}{\mathrm{\Σ}}}{w}_{k|k}^{\left(l\right)}\·N(x,m_{k|k-1}^{\left(l\right)},P_{k|k-1}^{\left(l\right)})+\underset{z_j\∈Z_k^n}{\mathrm{\Σ}}\underset{l=1}{\overset{J_{k|k-1}^n}{\mathrm{\Σ}}}{w}_{k|k}^{(l,j)}\·N(x,m_{k|k}^{(l,j)},P_{k|k}^{(l,j)})---\left(5\right),$

其中：

表示均值为协方差阵为的标准正态分布；

表示均值为协方差阵为的标准正态分布；

为直接预测所对应的高斯项权重，为第l个预测高斯项对应的检测概率，取值范围为区间[0.4,1)；

$w_{k|k}^{(l,j)}=\frac{P_D^{\left(l\right)}{w}_{k|k-1}^{\left(l\right)}N(z_j;{\mathrm{\η}}_{k|k-1}^{\left(l\right)},S_k^{\left(l\right)})}{K_k\left(z_j\right)+\underset{q=1}{\overset{J_{k|k-1}^n}{\mathrm{\Σ}}}{P}_D^{\left(q\right)}{w}_{k|k-1}^{\left(q\right)}\·N(z_j;{\mathrm{\η}}_{k|k-1}^{\left(q\right)},S_k^{\left(q\right)})},$

为第l个高斯项和第j个测量值对应的测量修正高斯项权重，式中求和公式中的上标q表示自然数，表示均值为协方差阵为的标准正态分布，表示均值为协方差阵为的标准正态分布，K_k(z_j)为杂波强度函数；

为经测量值z_j修正后的高斯项对应的均值；

为经测量值z_j修正后的高斯项对应的协方差阵，

得到测量修正后的高斯成份集合

其中为测量修正后所得高斯项的个数，f为自然数，取值范围为1,2....,上标f表示对应第f个测量修正后所得高斯项。

7.根据权利要求6所述的雷达系统的多目标跟踪方法，其特征在于：步骤5中，删除测量修正后的高斯成份集合中权重的高斯项，其中为裁减门限，取然后将距离充分小的高斯成份合并成一个，即将的高斯项合并成一个，其中为合并门限，其取值范围为区间[4,6]，距离d_i,j定义为： $d_{i,j}={(m_{k|k}^{\left(i\right)}-m_{k|k}^{\left(j\right)})}^T{(P_{k|k}^{\left(i\right)}+P_{k|k}^{\left(j\right)})}^{-1}(m_{k|k}^{\left(i\right)}-m_{k|k}^{\left(j\right)}),$

高斯成份合并方法如下：

合并后对应高斯项的权重

合并后对应高斯项的均值 ${\tilde{m}}_{k|k}^{\left(l\right)}=\frac{1}{{\tilde{w}}_{k|k}^{\left(l\right)}}\underset{i\∈L}{\mathrm{\Σ}}{w}_{k|k}^{\left(i\right)}{m}_{k|k}^{\left(i\right)},$

合并后对应高斯项的协方差阵 ${\tilde{P}}_{k|k}^{\left(l\right)}=\frac{1}{{\tilde{w}}_{k|k}^{\left(l\right)}}\underset{i\∈L}{\mathrm{\Σ}}{w}_{k|k}^{\left(i\right)}(P_{k|k}^{\left(i\right)}+({\tilde{m}}_{k|k}^{\left(l\right)}-m_{k|k}^{\left(i\right)}){({\tilde{m}}_{k|k}^{\left(l\right)}-m_{k|k}^{\left(i\right)})}^T),$ 其中L为所有可合并高斯成份的指标集合；将融合修剪后的高斯成份根据其方位分别放入对应扇区及其相邻的两个扇区，于是得到k周期第n个扇区高斯成份集合： ${\{w_k^{\left(i\right)},m_k^{\left(i\right)},P_n^{\left(i\right)},n\}}_{i=1}^{J_k^n}.$

8.根据权利要求7所述的雷达系统的多目标跟踪方法，其特征在于：步骤6中，将权重大于0.5的高斯成份均值提取出来作为第n个扇区前一扇区目标状态估计值。

9.根据权利要求8所述的雷达系统的多目标跟踪方法，其特征在于：当k周期第n个扇区后一扇区的测量值集合发送到数据处理器以后，重复步骤2，实现雷达系统中跨扇区全平面目标状态估计。