多雷达系统中目标检测跟踪联合处理方法与流程

本发明涉及通信技术领域，更进一步涉及雷达信号处理技术领域中的一种多雷达系统中目标检测跟踪联合处理方法。本发明可用于实现多雷达系统对目标的检测跟踪联合处理。

背景技术：

目标跟踪在雷达及声呐系统的很多应用中扮演着重要的角色，在战场监控、空防、空中交通管制和火控等方面都起着重要作用。但是在密集杂波环境下进行目标跟踪，在每次扫描的过程中，会得到大量量测数据，但每个数据的来源是未知的(可能来源于目标，也可能来源于虚警)。

B K Habtemariam,R Tharmarasa,and T Kirubarajan等人在其发表的论文“Multiple detection probabilistic data association filter for multistatic target tracking”(Proceedings of the,International Conference on Information Fusion.IEEE,2011:1-6.)和论文“A Multiple-Detection Joint Probabilistic Data Association Filter”(IEEE Journal of Selected Topics in Signal Processing,2013,7(3):461-471)中提出了一种多探测概率数据互联(Multiple Detection Probabilistic data association,MD-PDA)方法，用组合的联合事件来处理多个量测来源于同一个目标的概率。该方法用于在多雷达系统中对单目标进行跟踪，针对在跟踪门内一个目标可能产生多个量测的情况，并计算每个联合事件的概率，然后基于得到的概率将量测与目标关联。该方法解决了概率数据互联法(Probabilistic Data Association,PDA)在多雷达系统中进行跟踪时的不足。在多雷达系统中存在一种场景，互联门内同一个目标会产生多个量测，而PDA方法假设在每次扫描中，对于目标最多有一个量测以其为源。但是，该方法存在的不足之处是，MD-PDA假设一个目标可能会产生多个量测，但没有考虑一个雷达从一个目标至多获得一个量测的限制，因此需要处理大量组合联合事件；而且MD-PDA方法的检测器依据的是奈曼皮尔逊(Neyman-Pearson,NP)准则，采用了固定门限。在这种情况下，信息只能从检测器传到跟踪器，无法根据目标运动位置自适应地调整检测门限。

西安电子科技大学在其申请的专利“基于贝叶斯理论的多传感器检测跟踪联合处理方法”(专利申请号：201110003111.4，公布号：CN 102147468B)中公开了一种基于贝叶斯理论的多传感器检测跟踪联合处理方法。该方法是基于贝叶斯理论的检测跟踪联合处理方法,直接处理原始未阈值化的观测数据，利用目标的运动模型在帧数据间以概率的形式对信号进行积累，且能够有效利用先验信息，从而提高系统的检测性能。利用先验信息和所有观测数据计算目标状态后验概率密度函数，并在计算所得的目标状态后验概率密度函数的基础上实现对目标检测和跟踪。该方法推导出了检测阈值与给定虚警率的关系，但是，仍然存在的不足之处是，该方法依据的检测准则是NP准则，跟踪门内所有检测单元的门限相同，没有合理应用从跟踪器到检测器的反馈信息，无法根据目标跟踪的位置来自适应调整波门内各检测单元的门限。

技术实现要素：

本发明的目的在于克服上述现有方法的不足，提出了一种多雷达系统中目标检测跟踪联合处理方法(Joint Detection and Tracking Processing,JDTP)方法，以增强多雷达系统对目标的检测跟踪性能。

实现本发明的基本思想是首先建立目标的运动模型以及多雷达系统对目标量测的模型，采用贝叶斯检测器，根据从跟踪器反馈得到目标的跟踪状态，设置恒虚警检测门限，在波门内平均虚警率恒定的前提下，提升目标的平均检测概率和系统的跟踪性能；同时在MD-PDA的基础上降低组合联合事件的数量，降低计算复杂度低，更符合雷达应用真实情况。

本发明的具体步骤如下：

(1)按照下式，构建目标运动模型：

ξ_Tk＝Fξ_Tk-1+u_Tk-1

其中，ξ_Tk表示Tk时刻目标的运动状态，F表示匀速运动情况下目标状态的转移矩阵，ξ_Tk-1表示Tk-1时刻目标的运动状态，u_Tk-1表示分布服从均值为零、协方差为Q_Tk-1的高斯分布的运动过程噪声；

(2)按照下式，构建目标量测模型：

其中，表示t时刻多雷达系统中第n个雷达对目标测量得到过门限值的第j个含噪声的量测值，h_n,t(ξ_t)表示t时刻多雷达系统中第n个雷达对目标状态ξ_t的量测值，w_n,t表示t时刻多雷达系统中第n个雷达分布服从均值为零、协方差为Σ_n,t的高斯分布的量测噪声，v_n,t表示在t时刻在多雷达系统中第n个雷达跟踪波门内服从均匀分布的虚假量测值；

(3)初始化：

(3a)将目标状态的先验信息作为待跟踪目标初始时刻的跟踪状态；

(3b)按照下式，计算目标运动状态初始协方差：

C₀＝CRLB(y₀)

其中，C₀表示目标运动状态初始协方差，CRLB(y₀)表示初始时刻目标运动的跟踪状态y₀的克拉美罗限；

(4)选择目标检测的雷达：

从多雷达系统中，任选一个雷达作为目标检测的雷达；

(5)确定恒虚警检测门限值：

(5a)按照下式，计算目标检测雷达的检测器在恒虚警率条件时的工作参数：

其中，η_k表示k时刻目标检测雷达的检测器在恒虚警率条件时的工作参数，V_k表示k时刻目标检测雷达跟踪波门的大小，μ_k表示k时刻目标检测的雷达得到回波的平均信噪比，P表示目标检测的雷达所要求的虚警率；n_z表示目标检测雷达对目标量测的维数，Σ_k表示k时刻目标检测雷达的量测噪声协方差；

(5b)按照下式，计算目标检测雷达的恒虚警检测门限值：

其中，γ_k表示k时刻目标检测雷达的检测门限值，μ_k表示k时刻目标检测的雷达得到回波的平均信噪比，ln(·)表示取自然对数操作，η_k表示k时刻目标检测雷达的检测器在恒虚警率条件时的工作参数，V_k表示k时刻目标检测的雷达跟踪波门的大小，表示服从均值为方差为的正态分布，表示k时刻目标检测的雷达在跟踪门中的第l个检测单元处的量测值，表示k时刻目标检测的雷达对目标预测状态的量测值，表示坐标变换后k时刻目标检测的雷达对目标预测状态的量测协方差值；

(6)确定有效量测：

将目标检测雷达量测值中所有大于检测门限值的量测作为有效量测；

(7)判断多雷达系统中所有的目标检测雷达是否选择完，若是，则执行步骤(8)，否则，执行步骤(4)；

(8)按照下式，计算组合联合事件的数量：

其中，N′_e表示e时刻组合联合事件的数量，S表示多雷达系统中雷达的总数，∏表示连乘操作，s表示多雷达系统中各个雷达的标号，表示从M_s,r个量测中选择1个量测来源于目标的排列组合操作，表示r时刻多雷达系统中第s个雷达所有M_s,r个量测都来源于虚警，M_s,r表示r时刻多雷达系统中第s个雷达的有效量测总数，r的取值与e相同；

(9)按照下式，更新当前目标跟踪状态：

其中，y_d表示d时刻更新后的目标跟踪状态，∑表示求和操作，y_J,d表示d时刻组合联合事件中第J个事件更新后的状态，J表示在所有组合联合事件中的标号，表示d时刻所有组合联合事件中第J个事件发生的概率；

(10)按照下式，更新目标跟踪状态的协方差：

其中，C_d表示更新后的目标状态协方差，∑表示求和操作，J表示在所有组合联合事件中的标号，表示c时刻所有组合联合事件中第J个事件发生的概率，C_J,d表示y_J,d的协方差值，y_J,d表示d时刻组合联合事件中第J个事件更新后的状态，(·)^T表示矩阵转置操作，y_d表示d时刻更新后的目标跟踪状态；

(11)判断目标跟踪航迹是否发散，若是，则执行步骤(12)，否则，当前时刻加1，将跟踪器当前的跟踪状态反馈到检测中心后执行步骤(4)；

(12)目标跟踪结束。

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

第一，由于本发明通过将跟踪器当前的跟踪状态反馈到检测中心，计算目标检测雷达的恒虚警检测门限值，克服现有技术无法根据目标跟踪的位置来自适应调整波门内各检测单元的门限，使得本发明设计的在波门内平均虚警率恒定的前提下，提升目标的平均检测概率和系统的跟踪性能。

第二，由于本发明根据计算多雷达检测跟踪系统中组合联合事件的数量，更新当前目标跟踪状态，克服了现有技术需要处理大量组合联合事件的不足，使得本发明具有降低计算复杂度低，更符合雷达应用真实情况的优点。

附图说明

图1为本发明的流程图；

图2为目标量测模型的示意图；

图3为多雷达系统的雷达与目标方位示意图；

图4为多雷达系统每个雷达采用不同方法检测器的检测概率对比曲线；

图5为NP检测器和贝叶斯检测器门限的示意图；

图6为多雷达系统采用不同方法对目标跟踪的精度对比图。

具体实施方式

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

参照图1，本发明的具体实施步骤如下：

步骤1，构建目标运动模型。

按照下式，构建目标运动模型：

ξ_Tk＝Fξ_Tk-1+u_Tk-1

其中，ξ_Tk表示Tk时刻目标的运动状态，F表示匀速运动情况下目标状态的转移矩阵，ξ_Tk-1表示Tk-1时刻目标的运动状态，u_Tk-1表示分布服从均值为零、协方差为Q_Tk-1的高斯分布的运动过程噪声。

匀速运动情况下目标状态的转移矩阵F由下式给出：

其中，表示矩阵直积操作，T₀表示采样间隔。

运动噪声协方差矩阵Q_k-1由下式给出：

其中q₁表示控制目标动态模型的过程噪声强度的参数。

步骤2，构建目标量测模型。

按照下式，构建目标量测模型：

其中，表示t时刻多雷达系统中第n个雷达对目标测量得到过门限值的第j个含噪声的量测值，h_n,t(ξ_t)表示t时刻多雷达系统中第n个雷达对目标状态ξ_t的量测值，w_n,t表示t时刻多雷达系统中第n个雷达分布服从均值为零、协方差为Σ_n,t的高斯分布的量测噪声，v_n,t表示在t时刻在多雷达系统中第n个雷达跟踪波门内服从均匀分布的虚假量测值。

参照附图2所示的两个雷达观测目标示意图，其中，方框表示雷达1的的量测来源于目标，三角表示雷达1的量测来源于虚警，圆点表示雷达2的的量测来源于目标，菱形表示雷达2的量测来源于虚警，星型表示目标预测位置，V_1,t和V_2,t分别表示雷达1和雷达2的跟踪波门。

步骤3，初始化。

将目标状态的先验信息作为待跟踪目标初始时刻的跟踪状态。

本发明的实施例中将先验信息作为目标航迹起始的目标状态。

按照下式，计算目标运动状态初始协方差：

C₀＝CRLB(y₀)

其中，C₀表示目标运动状态初始协方差，CRLB(y₀)表示初始时刻目标运动的跟踪状态y₀的克拉美罗限。

步骤4，选择目标检测的雷达。

从多雷达系统中，任选一个雷达作为目标检测的雷达。

步骤5，确定恒虚警检测门限值。

按照下式，计算目标检测雷达的检测器在恒虚警率条件时的工作参数：

其中，η_k表示k时刻目标检测雷达的检测器在恒虚警率条件时的工作参数，V_k表示k时刻目标检测雷达跟踪波门的大小，μ_k表示k时刻目标检测的雷达得到回波的平均信噪比，P表示目标检测的雷达所要求的虚警率；n_z表示目标检测雷达对目标量测的维数，Σ_k表示k时刻目标检测雷达的量测噪声协方差。

本发明的实施例中雷达接收到的回波幅度的概率密度可表示为：

其中，H₀表示该位置没有目标，这时接收的数据仅来自于噪声。H₁表示该位置存在目标，得到的数据将是目标的回波和噪声的叠加。表示k时刻检测雷达在第l个分辨单元的回波幅度(在目标起伏模型为Swerling I型的前提下，I、Q两路的平方和输出为指数分布)，μ_k表示k时刻检测雷达得到回波的平均信噪比。

按照下式，给出目标是否存在的先验概率：

其中，表示目标不存在的先验概率，表示目标存在的先验概率，V_k表示k时刻检测雷达跟踪波门的大小，表示k时刻检测雷达对目标预测量测和表示k时刻雷达对目标预测协方差。h_k(ξ_k)表示k时刻检测雷达对目标状态的量测值，表示服从均值为方差为的正态分布。

根据在跟踪过程反馈得到的先验信息，将判别式写成如下的贝叶斯形式：

其中和表示检测雷达在第l个分辨单元处量测的预测分布，H₀表示该位置没有目标的情况，H₁表示该位置存在目标的情况，η_k表示k时刻检测雷达贝叶斯检测器在恒虚警率条件时的常数。

计算平均虚警率由下式计算得到：

其中，P表示平均虚警率，V_k表示k时刻检测雷达跟踪波门的大小，exp(·)表示指数操作，μ_k表示k时刻检测雷达得到回波的平均信噪比，η_k表示k时刻检测雷达贝叶斯检测器在恒虚警率条件时的常数，P表示检测雷达所要求的虚警率；n_z表示雷达对目标量测的维数，Σ_k表示k时刻检测雷达的量测噪声协方差。

所以，给定平均虚警率，可以得到贝叶斯检测器在恒虚警率条件时的常数η_k。

按照下式，计算目标检测雷达的恒虚警检测门限值。

其中，γ_k表示k时刻目标检测雷达的检测门限值，μ_k表示k时刻目标检测的雷达得到回波的平均信噪比，ln(·)表示取自然对数操作，η_k表示k时刻目标检测雷达的检测器在恒虚警率条件时的工作参数，V_k表示k时刻目标检测的雷达跟踪波门的大小，表示服从均值为方差为的正态分布，表示k时刻目标检测的雷达在跟踪门中的第l个检测单元处的量测值，表示k时刻目标检测的雷达对目标预测状态的量测值，表示坐标变换后k时刻目标检测的雷达对目标预测状态的量测协方差值；多雷达系统的检测门限通过将贝叶斯判别式进行变换可以得到：

其中，表示k时刻检测雷达在第l个分辨单元的回波幅度，将其与检测门限γ_k比较，可以判断量测来源于目标还是虚警。从上式可以看出，若越靠近则越大，检测门限γ_k越低。换句话说，贝叶斯定理检测器的检测概率在不同的位置处是不同的。

步骤6，确定有效量测。

将目标检测雷达量测值中所有大于检测门限值的量测作为有效量测。

步骤7，判断多雷达系统中所有的目标检测雷达是否选择完，若是，则执行步骤8，否则，执行步骤4。

步骤8，按照下式，计算组合联合事件的数量。

其中，N′_e表示e时刻组合联合事件的数量，S表示多雷达系统中雷达的总数，∏表示连乘操作，s表示多雷达系统中各个雷达的标号，表示从M_s,r个量测中选择1个量测来源于目标的排列组合操作，表示r时刻多雷达系统中第s个雷达所有M_s,r个量测都来源于虚警，M_s,r表示r时刻多雷达系统中第s个雷达的有效量测总数，r的取值与e相同。

本发明的实施例中组合联合事件为目标量测模型量测来源于目标或者虚警的所有可能事件。

通常用于解决在多雷达系统中单目标的跟踪问题的方法为MD-PDA方法，联合事件的数量定义为N_q。

其中，N为多雷达系统中雷达的总数，min(·)表示取最小值操作，M_s表示多雷达系统接收到的量测总数，a表示量测来源于目标的量测，∑表示求和操作，C表示排列组合操作，(·)！表示阶乘操作。然而，在实际应用中，更真实的情况为一个雷达从一个目标至多能收到一个量测。在这种情况下，联合事件的数量定义为N′_s：

其中，N′_e表示e时刻组合联合事件的数量，S表示多雷达系统中雷达的总数，∏表示连乘操作，s表示多雷达系统中各个雷达的标号，表示从M_s,r个量测中选择1个量测来源于目标的排列组合操作，表示r时刻多雷达系统中第s个雷达所有M_s,r个量测都来源于虚警，M_s,r表示r时刻多雷达系统中第s个雷达的有效量测总数。

步骤9，按照下式，更新当前目标跟踪状态。

其中，y_d表示d时刻更新后的目标跟踪状态，∑表示求和操作，y_J,d表示d时刻组合联合事件中第J个事件更新后的状态，J表示在所有组合联合事件中的标号，表示d时刻所有组合联合事件中第J个事件发生的概率。

本发明的实施例中d时刻组合联合事件中第J个事件更新后的状态y_J,d由下式得到：

其中，y_J,d表示d时刻组合联合事件中第J个事件更新后的状态，表示d时刻多雷达系统对目标的预测状态，表示d时刻第b个目标检测雷达的第J(b)个量测的卡尔曼增益，表示d时刻第b个雷达第J(b)个量测的相关新息。

本发明的实施例中d时刻所有组合联合事件中第J个事件发生的概率是由下式得到的：

其中，表示d时刻所有组合联合事件中第J个事件发生的概率，∏表示连乘操作，i表示多雷达系统中的雷达标号，N表示多雷达系统中雷达的总数，p{·}表示求概率操作，表示m时刻多雷达系统中第i个雷达的第J(i)个量测来源于目标的事件，表示第i个雷达的测量全都来源于虚警。

步骤10，按照下式，更新目标跟踪状态的协方差：

其中，C_d表示更新后的目标状态协方差，∑表示求和操作，J表示在所有组合联合事件中的标号，表示c时刻所有组合联合事件中第J个事件发生的概率，C_J,d表示y_J,d的协方差值，y_J,d表示d时刻组合联合事件中第J个事件更新后的状态，(·)^T表示矩阵转置操作，y_d表示d时刻更新后的目标跟踪状态。

步骤11，判断目标跟踪航迹是否发散，若是，则执行步骤12，否则，当前时刻加1，将跟踪器当前的跟踪状态反馈到检测中心后执行步骤4。

在多雷达系统中，当目标不在多雷达系统的检测范围内时，目标跟踪的航迹就会发散。

步骤12，目标跟踪结束。

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

1.仿真条件：

本发明的仿真运行系统为Intel(R)Core(TM)i7-4590CPU@3.30GHz，32位Windows操作系统，仿真软件采用MATLAB(R2014b)。

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

本发明的仿真实验设定一个由在6个不同地点的雷达组成的雷达网，每个雷达的信号有效带宽和有效持续时间分别为B＝2MHz和T＝1ms，仿真序列数据为60帧，观测间隔T₀＝2s，门系数g＝8；目标初始位置在(0,0)km，并以速度(-30,0)m/s作匀速飞行，当R₀＝10km时信噪比μ₀＝8dB。

参照图3多雷达系统的雷达与目标方位示意图，本发明仿真实验针对一种场景：一个目标接近部署在边界上的雷达网。

图3中方框表示多雷达系统中的雷达，虚线表示目标运动的轨迹，箭头表示目标运动方向。本发明的实施例是针对下面两种情况进行仿真：(1)第一种情况的虚警率P_i＝10^-2；(2)第二种情况的虚警P_i＝10^-3，i＝1,2,...,N。

参照图4所示的采用本发明在上述两种情况下，多雷达系统每个雷达采用NP检测器和贝叶斯检测器下的检测概率对比曲线。图4中的纵坐标表示检测概率，横坐标表示对目标跟踪的时刻。图4中以点虚线标示的曲线表示图3部署图中的雷达1，图4中以加号线标示的曲线表示图3部署图中的雷达2，图4中以点划线标示的曲线表示图3部署图中的雷达3，图4中以实线标示的曲线表示图3部署图中的雷达4，图4中以叉号线标示的曲线表示图3部署图中的雷达5，图4中以虚线标示的曲线表示图3部署图中的雷达6。图4(a)所示为多雷达系统在第一种情况下采用NP检测器的检测概率曲线，图4(b)所示为多雷达系统在第一种情况下采用贝叶斯检测器的检测概率曲线。图4(c)所示为多雷达系统在第二种情况下采用NP检测器的检测概率曲线，图4(d)所示为多雷达系统在第二种情况下采用贝叶斯检测器的检测概率曲线。

比较图4(a)和图4(b)可以看到，在第一种情况下，采用贝叶斯检测器的检测性能要优于NP检测器。比较图4(c)和图4(d)，同样可以看到采用贝叶斯检测器的检测性能要优于NP检测器。因此，综合图4可以看出，采用本发明的贝叶斯检测器，虚警率越高，对检测性能的提升越好。

参照图5所示，本发明取第10个时刻的图3部署图中的雷达1作为例子，比较本发明的贝叶斯检测器和NP检测器的检测门限。图5中的纵坐标表示雷达对目标的检测门限，横坐标表示有效波门内的检测单元范围。图5中的实线表示采用NP检测器的检测门限曲线，虚线表示采用本发明的贝叶斯检测器的检测门限曲线，实心矩形表示目标在波门内所处的检测单元位置。图5(a)所示为本实施例在第一种情况下采用贝叶斯检测器和NP检测器检测门限的对比图，图5(b)所示为本实施例在第二种情况下采用贝叶斯检测器和NP检测器检测门限的对比图。

综合图5可以看出，贝叶斯检测门限在反馈的预测中心附近区域较低；越远离预测中心，检测门限越高。

表1为多探测概率数据互联MD-PDA方法和本发明JDTP方法联合事件的数量对比一览表。从表1可以看到，在第一种情况下，JDTP方法相对多探测概率数据互联MD-PDA方法得到更少的联合事件；在第二种情况下，两种方法得到的联合事件数几乎相同。在高虚警率的情况下，相比于原始的多探测概率数据互联MD-PDA方法，本发明可以明显降低降低联合事件的数量；如果虚警率太低，与原始的多探测概率数据互联MD-PDA方法比较，本发明就不能明显降低联合事件的数量。

表1联合事件的数量一览表

参照图6所示的采用本发明比较了标准概率数据互联PDA方法、多探测概率数据互联方法MD-PDA以及本发明JDTP方法的跟踪精度，精度采用均方根误差(Root MSEs,RMSE)来衡量。图6中纵坐标表示均方根误差，横坐标表示对目标跟踪的时刻，图6种点虚线表示标准概率数据互联PDA方法的均方根误差RMSE，带叉号的虚线表示多探测概率数据互联MD-PDA方法的均方根误差RMSE，虚线表示多探测概率数据互联MD-PDA方法的贝叶斯克拉美罗下限BCRLB，带加号的虚线表示本发明JDTP方法的均方根误差RMSE，实线表示本发明JDTP方法的贝叶斯克拉美罗下限BCRLB。图6(a)表示多雷达系统在第一种情况下采用不同方法对目标跟踪的精度对比图，图6(b)表示多雷达系统在第二种情况下采用不同方法对目标跟踪的精度对比图。

综合图6可以看出，在高虚警率的情况下，本发明相对于其他两种方法在估计精度方面具有较高的提高；在低虚警率的情况下，采用本发明的贝叶斯检测器仅仅能将跟踪精度轻微提高。