一种分布式外辐射源雷达组网探测的航迹融合方法

摘要

一种分布式外辐射源雷达组网探测的航迹融合方法，包括以下步骤：(1)外辐射源雷达各分布站点分别探测当前时刻的目标信息，并解析对应的当前时刻的目标状态，外辐射源雷达各分布站点获得其对应的当前时刻的局部航迹；(2)融合中心对所有来自于同一目标的外辐射源雷达分布站点的各个当前时刻的局部航迹进行加权融合，得到当前时刻的加权融合后的航迹；(3)在融合中心设置IMM滤波器，对当前时刻的加权融合后的航迹进行IMM滤波，计算得到当前时刻IMM滤波的交互输出状态向量和交互输出协方差矩阵；(4)如果外辐射源雷达各分布站点继续跟踪目标，则令当前时刻增加1，返回步骤1；否则，IMM滤波结束。

说明书

技术领域

本发明属于雷达信号处理技术领域，涉及一种分布式外辐射源雷达组网探测的航迹融 合方法，进一步涉及一种基于交互式多模型(IMM)算法的分布式外辐射源雷达组网探测 的航迹融合方法，用于在分布式外辐射源雷达组网探测中对目标进行航迹融合。

背景技术

外辐射源雷达指的是一种本身不发射电磁波，而利用非合作照射源对目标进行探测、 定位和跟踪的新体制雷达，因为其发射源的信号波形和信号功率不受控，且辐射源的空间 分布和覆盖范围受限，这就导致单站外辐射源雷达的观测信息中包含较多的杂波和虚警， 所探测的目标信息的置信度低。伴随着近年来通信技术的发展，各雷达站点间实时通信成 为了可能，由此推进了利用多站外辐射源雷达的观测信息，通过时间和空间的信息融合， 来改善多目标定位和跟踪性能的分布式外辐射源雷达组网探测的发展。目前各种体制雷达 所采用的多站组网探测主要包括三种结构：第一种是雷达各分布站点分别探测目标信息 (距离和多普勒等)，但各站点不具备解析目标状态(目标位置和速度等)的能力，只是 把探测到的目标信息传送到融合中心，由融合中心来解析目标状态，这种架构称为集中式 雷达组网探测；第二种是雷达各分布站点不但能探测目标信息，同时各站点能够根据各自 的探测信息解析目标状态，然后将解析的目标状态传送到融合中心，在融合中心进行目标 状态的融合，这种结构称为分布式雷达组网探测；第三种是雷达各分布站点中只有部分站 点同时具备目标信息探测和目标状态解析的能力，其余站点只能探测目标信息，具备目标 状态解析能力的站点将目标状态传送到融合中心，只具备探测能力的站点将探测到的目标 信息传送到融合中心，在融合中心站点完成目标信息和目标状态的融合，这种结构称为混 合式雷达组网探测。在上述三种结构中，分布式雷达组网探测相比于集中式雷达组网探测 具有传送信息量少，通信带宽窄的优势；同时，分布式雷达组网探测相比于混合式雷达组 网探测具有结构简单，造价低廉的优势，因此，分布式雷达组网探测在外辐射源雷达组网 探测中得到了较为广泛的应用。

分布式外辐射源雷达组网探测在融合中心进行目标状态融合时，较常见的是进行目标 航迹的融合，即对各站点针对同一目标局部解析形成的航迹进行融合。目前实践中采用的 航迹融合算法通常是按照一定的融合准则对各站点传送过来的目标状态进行加权，将加权 的结果作为最终的目标航迹融合的结果。这样的加权融合，由于未考虑各站点对目标进行 局部解析时因运动模型单一引入的模型误差，从而导致该模型误差仍然保留在目标航迹融 合的结果中，造成了航迹融合精度的下降。

发明内容

针对上述已有技术的不足，本发明的目的在于提出一种分布式外辐射源雷达组网探测 的航迹融合方法，该方法提出基于IMM算法对加权融合后的航迹进行滤波，即通过多模型 交互转换更加精确地得到目标的运动模型；本发明方法能够克服分布式外辐射源雷达组网 探测中融合中心未考虑分布站点的模型误差而造成的融合航迹精度差的缺点，因而能够提 高对目标航迹融合的精度。

为达到上述技术目的，本发明采用以下技术方案予以实现。

一种分布式外辐射源雷达组网探测的航迹融合方法，其特征在于，包括以下步骤：

步骤1，外辐射源雷达各分布站点分别探测当前时刻的目标信息，并分别解析得到各 分布站点对应的当前时刻的目标状态，外辐射源雷达各分布站点获得其对应的当前时刻的 局部航迹，该局部航迹包含状态向量和协方差矩阵；外辐射源雷达各分布站点将获得的各 个当前时刻的局部航迹传送到融合中心；

步骤2，融合中心根据外辐射源雷达每个分布站点传送过来的对应的当前时刻的局部 航迹中所包含的状态向量和协方差矩阵，判断外辐射源雷达各分布站点的各个当前时刻的 局部航迹是否来自于同一目标，对所有来自于同一目标的外辐射源雷达分布站点的当前时 刻的局部航迹进行加权融合，得到当前时刻的加权融合后的航迹，该当前时刻的加权融合 后的航迹包含当前时刻的融合状态向量和当前时刻的融合协方差矩阵；

步骤3，在融合中心设置IMM滤波器，对步骤2所述当前时刻的加权融合后的航迹进 行IMM滤波，计算得到当前时刻IMM滤波的交互输出状态向量和交互输出协方差矩阵；

步骤4，如果外辐射源雷达各分布站点继续跟踪目标，则令当前时刻增加1，返回步 骤1；否则，IMM滤波结束，保存步骤3所得当前时刻IMM滤波的交互输出状态向量和交 互输出协方差矩阵。

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

第一，本发明在融合中心基于IMM算法对加权融合后的航迹进行滤波，在滤波过程中 将目标运动模型设定为多种运动模型的叠加，从而减小了外辐射源雷达各分布站点在局部 滤波时因为使用单一目标运动模型而引入的模型误差，提高了对目标航迹融合的精度；

第二，本发明在融合中心设置IMM滤波器，避免了在外辐射源雷达各分布站点同时使 用IMM滤波器，减少了滤波器的个数，从而减少了运算量，并加快了运算速度，节省了存 储空间；同时本发明可根据实际需要在融合中心增加或减少目标运动模型的个数，而不需 要改变融合中心的IMM滤波器的结构，可移植性强；

第三，本发明方法使得分布式外辐射源雷达组网探测的整体结构更为简单，造价低廉， 应用面广。

附图说明

下面结合附图说明和具体实施方式对本发明作进一步详细说明。

图1为本发明的IMM滤波的数据流向图；

图2为两种方法所得的目标航迹沿X方向的航迹精度对比图；

图3为两种方法所得的目标航迹沿X方向的航迹精度对比图；

图4为实测实验的航迹误差图。

具体实施方式

参照图1，本发明的一种分布式外辐射源雷达组网探测的航迹融合方法，包括以下具 体步骤：

步骤1，外辐射源雷达各分布站点分别探测当前时刻的目标信息，并分别解析得到各 分布站点对应的当前时刻的目标状态，外辐射源雷达各分布站点获得其对应的当前时刻的 局部航迹，该局部航迹包含状态向量和协方差矩阵；外辐射源雷达各分布站点将获得的各 个当前时刻的局部航迹传送到融合中心。

步骤2，融合中心根据外辐射源雷达每个分布站点传送过来的对应的当前时刻的局部 航迹中所包含的状态向量和协方差矩阵，判断外辐射源雷达各分布站点的各个当前时刻的 局部航迹是否来自于同一目标，对所有来自于同一目标的外辐射源雷达分布站点的当前时 刻的局部航迹进行加权融合，得到当前时刻的加权融合后的航迹，该当前时刻的加权融合 后的航迹包含当前时刻的融合状态向量和当前时刻的融合协方差矩阵；

步骤2的具体子步骤为：

2.1按照Singer-Kanyuck算法判断外辐射源雷达各分布站点的各个当前时刻的局部航 迹是否来自于同一目标；

2.2对步骤2.1中所有判断为来自于同一目标的外辐射源雷达分布站点的各个当前时 刻的局部航迹按照Bar-Shalom-Campo算法进行加权融合，得到当前时刻的加权融合后的 航迹，该当前时刻的加权融合后的航迹包括：当前时刻的融合状态向量、当前时刻的融合 协方差矩阵；将当前时刻记为k时刻，当前时刻的融合状态向量记为当前时刻 的融合协方差矩阵记为

步骤3，在融合中心设置IMM滤波器，对步骤2所述当前时刻的加权融合后的航迹进 行IMM滤波，计算得到当前时刻IMM滤波的交互输出状态向量和交互输出协方差矩阵；

步骤3的具体子步骤为：

3.1在融合中心设置IMM滤波器，并设定IMM滤波初始时刻，如果当前时刻是IMM滤 波初始时刻，初始化目标运动模型集，该目标运动模型集包含常速度模型和常加速度模型， 设定IMM滤波初始时刻目标运动模型为常速度模型的概率μ₁(0)＝0.5，设定IMM滤波初始时 刻目标运动模型为常加速度模型的概率μ₂(0)＝0.5；

构建常速度模型和常加速度模型的马尔科夫转移概率矩阵p，p的维数为2×2，p中 的元素p_ij表示模型i向模型j转移的概率，i∈{1，2}，j∈{1，2}，1表示常速度模型，2表示 常加速度模型，p中的同模型间转移概率，即p中的元素p₁₁和p₂₂在区间[0.8，1]中选取；p 中的不同模型间转移概率，即p中的元素p₁₂和p₂₁在区间[0，0.2]中选取；

根据IMM滤波初始时刻，目标运动模型分别为常速度模型和常加速度模型的概率μ₁(0) 和μ₂(0)，以及所构建的常速度模型和常加速度模型的马尔科夫转移概率矩阵p，计算IMM 滤波初始时刻的归一化马尔科夫转移概率矩阵的维数为2×2，中的元素表示为 i∈{1，2}，j∈{1，2}，1表示常速度模型，2表示常加速度模型；

将当前时刻记为k时刻；将步骤2.2中所述当前时刻的融合状态向量分别 作为IMM滤波初始时刻常速度模型和常加速度模型的卡尔曼滤波输出状态向量和 即${\hat{x}}_1\left(k\right|k)={\hat{x}}^f\left(k\right|k),{\hat{x}}_2\left(k\right|k)={\hat{x}}^f\left(k\right|k);$将步骤2.2中所述当前时刻的融合 协方差矩阵分别作为IMM滤波初始时刻常速度模型和常加速度模型的卡尔曼滤 波输出协方差矩阵和即${\hat{p}}_1\left(k\right|k)={\hat{p}}^f\left(k\right|k),{\hat{p}}_2\left(k\right|k)={\hat{p}}^f\left(k\right|k);$并 令当前时刻增加1，即令k时刻增加1，返回步骤1；

如果当前时刻不是IMM滤波初始时刻，则进行下一步操作；

3.2分别计算当前时刻常速度模型和常加速度模型的交互输入条件：即根据上一时刻 常速度模型和常加速度模型的卡尔曼滤波输出状态向量、上一时刻常速度模型和常加速度 模型的卡尔曼滤波输出协方差矩阵、上一时刻目标运动模型分别为常速度模型和常加速度 模型的概率以及上一时刻的归一化马尔科夫转移概率矩阵，分别计算当前时刻常速度模型 和常加速度模型的交互输入状态向量，以及当前时刻常速度模型和常加速度模型的交互输 入协方差矩阵；

具体地，将当前时刻记为k时刻；首先，计算k-1时刻模型j的归一化常数其中，j∈{1，2}，1表示常速度模型，2表示常加速度模型，的表达式为：

${\bar{c}}_j(k-1)=\underset{i=1}{\overset{2}{\Sigma }}{p}_{ij}{\mu }_i(k-1)$

其中，p_ij是常速度模型和常加速度模型的马尔科夫转移概率矩阵p的元素，p_ij表示模 型i向模型j转移的概率，μ_i(k-1)表示k-1时刻目标运动模型为模型i的概率；i∈{1，2}， 1表示常速度模型，2表示常加速度模型；

然后，根据k-1时刻模型j的归一化常数计算k-1时刻的归一化马尔科夫转 移概率矩阵中的元素的表达式为：

${\bar{\mu }}_{ij}^{k-1}=\frac{1}{{\bar{c}}_j(k-1)}{p}_{ij}{\mu }_i(k-1)$

最后，根据k-1时刻的归一化马尔科夫转移概率矩阵计算k时刻模型j的交互 输入状态向量以及k时刻模型j的的交互输入协方差矩阵其中， j∈{1，2}，1表示常速度模型，2表示常加速度模型，和的表达式分别为：

$x_j^{*}\left(k\right|k)=\underset{i=1}{\overset{2}{\Sigma }}{\hat{x}}_i(k-1|k-1){\bar{\mu }}_{ij}^{k-1}$

$p_j^{*}\left(k\right|k)=\underset{i=1}{\overset{2}{\Sigma }}{\bar{\mu }}_{ij}^{k-1}\{{\hat{p}}_i(k-1|k-1)+[{\hat{x}}_i(k-1|k-1)-x_j^{*}\left(k\right|k)]*{[{\hat{x}}_i(k-1|k-1)-x_j^{*}\left(k\right|k)]}^T\}$

其中，表示k-1时刻模型i的卡尔曼滤波输出状态向量，表 示k-1时刻模型i的卡尔曼滤波输出协方差矩阵，i∈{1，2}，1表示常速度模型，2表示常 加速度模型，上标T表示矩阵转置；

3.3将步骤2.2中所述当前时刻的融合状态向量分别作为当前时刻常速度 模型和常加速度模型的输入观测向量x₁(k|k)和x₂(k|k)，即将步骤2.2中所述当前时刻的融合协方差矩阵分别作为当 前时刻常速度模型和常加速度模型的输入观测方差矩阵p₁(k|k)和p₂(k|k)，即 $p_1\left(k\right|k)={\hat{p}}^f\left(k\right|k),p_2\left(k\right|k)={\hat{p}}^f\left(k\right|k);$结合步骤3.2所述当前时刻常速度模型和常加 速度模型的交互输入条件，对常速度模型和常加速度模型分别进行卡尔曼滤波，即分别计 算得到当前时刻常速度模型和常加速度模型的卡尔曼滤波输出状态向量和 以及当前时刻常速度模型和常加速度模型的卡尔曼滤波输出协方差矩阵和

3.4根据当前时刻常速度模型和常加速度模型的输入观测向量x₁(k|k)和x₂(k|k)、当 前时刻常速度模型和常加速度模型的输入观测方差矩阵p₁(k|k)和p₂(k|k)、步骤3.1所 述常速度模型和常加速度模型的马尔科夫转移概率矩阵p，以及步骤3.2所述当前时刻常 速度模型和常加速度模型的交互输入条件，分别计算当前时刻常速度模型和常加速度模型 的似然函数Λ₁(k)和Λ₂(k)；再根据当前时刻常速度模型和常加速度模型的似然函数Λ₁(k) 和Λ₂(k)，计算当前时刻目标运动模型分别为常速度模型和常加速度模型的概率μ₁(k)和 μ₂(k)；

具体地，将当前时刻记为k时刻；首先，计算k时刻模型j的似然函数Λ_j(k)，其中， j∈{1，2}，1表示常速度模型，2表示常加速度模型，Λ_j(k)的表达式为：

${\Lambda }_j\left(k\right)=p\left\{x_j\left(k\right|k)\right|M_j\left(k\right),x_j^{*}\left(k\right|k),p_j^{*}\left(k\right|k)\}$

其中，M_j(k)表示k时刻模型j有效，x_j(k|k)表示k时刻模型j的输入观测向量， 表示k时刻模型j的交互输入状态向量，表示k时刻模型j的交互输入协 方差矩阵；

然后，根据k时刻模型j的似然函数Λ_j(k)，计算k时刻似然函数的归一化常数c(k)， c(k)的表达式为：

$c\left(k\right)=\underset{j=1}{\overset{2}{\Sigma }}{\Lambda }_j\left(k\right){\bar{c}}_j(k-1)$

其中，表示k-1时刻模型j的归一化常数；

最后，根据k时刻似然函数的归一化常数c(k)，计算k时刻目标运动模型为模型j的 概率μ_j(k)，μ_j(k)的表达式为：

${\mu }_j\left(k\right)=\frac{1}{c\left(k\right)}{\Lambda }_j\left(k\right){\bar{c}}_j(k-1);$

3.5根据步骤3.3所述当前时刻常速度模型和常加速度模型的卡尔曼滤波输出状态向 量和以及当前时刻常速度模型和常加速度模型的卡尔曼滤波输出协方差 矩阵和并结合步骤3.4所述当前时刻目标运动模型分别为常速度模型 和常加速度模型的概率μ₁(k)和μ₂(k)，加权求得当前时刻IMM滤波的交互输出结果，即加 权求得当前时刻IMM滤波的交互输出状态向量和交互输出协方差矩阵

具体地，将当前时刻记为k时刻；首先，计算k时刻IMM滤波的交互输出状态向量 的表达式为：

${\hat{x}}^{*}\left(k\right|k)=\underset{j=1}{\overset{2}{\Sigma }}{\hat{x}}_j\left(k\right|k){\mu }_j\left(k\right)$

其中，表示k时刻模型j的卡尔曼滤波输出状态向量，μ_j(k)表示k时刻目标 运动模型为模型j的概率；

然后，计算k时刻IMM滤波的交互输出状态向量的表达式为：

${\hat{p}}^{*}\left(k\right|k)=\underset{j=1}{\overset{2}{\Sigma }}{\mu }_j\left(k\right)\{{\hat{p}}_j\left(k\right|k)+[{\hat{x}}_j\left(k\right|k)-{\hat{x}}^{*}\left(k\right|k)]*{[{\hat{x}}_j\left(k\right|k)-{\hat{x}}^{*}\left(k\right|k)]}^T\}$

其中，表示k时刻模型j的卡尔曼滤波输出协方差矩阵。

步骤4，如果外辐射源雷达各分布站点继续跟踪目标，则令当前时刻增加1，返回步 骤1；否则，IMM滤波结束，保存步骤3所得当前时刻IMM滤波的交互输出状态向量和交 互输出协方差矩阵。

本发明的效果可由以下仿真实验和实测实验作进一步说明：

1)仿真实验参数

仿真实验的数据是在总采样帧数为200帧的采样时间内由多种运动模型叠加产生，采 样周期T＝1s，目标初始位置在原点，初始速度为(200m/s，200m/s)，在1-50帧的采样时 间内采用的运动模型是匀速直线运动叠加标准方差为0.1m/s²的高斯白噪声驱动，在 51-150帧的采样时间内采用的运动模型是匀速圆周运动叠加标准方差为0.1m/s²的高斯白 噪声驱动，在151-200帧的采样时间内采用的运动模型是加速度为(5m/s²，5m/s²)的匀 加速直线运动叠加标准方差为0.1m/s²的高斯白噪声驱动，外辐射源雷达各分布站点的观 测噪声是标准方差为0.1m/s²的高斯白噪声。

2)仿真实验内容及结果分析：

根据上述仿真实验参数，进行500次蒙特卡洛实验，在单个外辐射源雷达站点运用卡 尔曼滤波计算得到目标航迹，简称卡尔曼滤波方法；同时，在多个外辐射源雷达站点运用 本发明方法计算得到目标航迹，简称本发明方法；将上述两种方法所得的目标航迹分别与 未叠加高斯白噪声驱动的目标的真实航迹进行比较，绘制两种方法所得的目标航迹沿X方 向的航迹精度对比图，如图2所示，图2中，横坐标为采样时刻，单位为帧，纵坐标为X 方向误差，单位为米(m)，虚线表示在单个外辐射源雷达站点运用卡尔曼滤波计算得到的 目标航迹与目标真实航迹沿X方向的误差，实线表示在多个外辐射源雷达站点运用本发明 方法计算得到的目标航迹与目标真实航迹沿X方向的误差；绘制两种方法所得的目标航迹 沿Y方向的航迹精度对比图，如图3所示，图3中，横坐标为采样时刻，单位为帧，纵坐 标为Y方向误差，单位为米(m)，虚线表示在单个外辐射源雷达站点运用卡尔曼滤波计算 得到的目标航迹与目标真实航迹沿Y方向的误差，实线表示在多个外辐射源雷达站点运用 本发明方法计算得到的目标航迹与目标真实航迹沿Y方向的误差；各个采样时刻的上述误 差取500次蒙特卡洛实验结果的几何平均值。

从图2和图3中可以看出：在匀速直线运动的采样时间内，两种方法所得的目标航迹 沿X方向和沿Y方向的误差比较一致，都能比较准确的跟踪目标；在匀速圆周运动的采样 时间内，由于匀速圆周运动在X方向和Y方向的速度和加速度都不是均匀变化的，所以两 种方法所得的目标航迹沿X方向和Y方向的航迹精度都开始变差，但本发明方法相比较于 卡尔曼滤波方法可以得到更精确的目标航迹；在匀加速直线运动的采样时间内，本发明方 法的优势更能体现出来，这是由于相比较于卡尔曼滤波方法只考虑单一目标运动模型，本 发明方法将目标运动模型考虑为多种运动模型的叠加，再通过IMM滤波获得多模型交互输 出的结果，减小了卡尔曼滤波方法存在的模型匹配误差，提高了对目标的跟踪精度。

3)实测实验参数：

实测实验数据在湖北省荆州市江北农场采集，两个外辐射源雷达站点分别位于东经 112.4°、北纬30.1252°和东经112.414°、北纬30.1251°，两个外辐射源雷达站点相距1.3465 公里，所探测的目标为正在飞行的民航飞机。两个外辐射源雷达站点的接收机具有同样的 软件、硬件架构，采集速率均为1帧/秒，总采样帧数为200帧。在实测实验中，将融合中 心已关联航迹设定为所选航迹，并设定所选航迹不存在错误关联和数据采集错误的情况， 与所选目标航迹对应的目标的真实航迹，称为标准航迹，该标准航迹包含所探测的目标(即 正在飞行的民航飞机)的真实状态信息，该标准航迹由民航ADS-B接收机获得；

4)实测实验内容及结果分析：

在实测实验中，首先，在所述两个外辐射源雷达站点中的一个外辐射源雷达站点运用 卡尔曼滤波得到目标航迹；然后，在所述两个外辐射源雷达站点运用本发明方法得到目标 航迹。图4为实测实验的航迹误差图，横轴是采样时刻，单位为帧，纵轴是距离误差，单 位是千米(Km)，实线表示在单个外辐射源雷达站点运动卡尔曼滤波所得到的目标航迹与 标准航迹的距离误差，虚线表示在两个外辐射源雷达站点运用本发明方法所得到的目标航 迹与标准航迹的距离误差。

从图4中可以看出，相比较于在一个外辐射源雷达站点运用卡尔曼滤波所得到的目标 航迹，在两个外辐射源雷达站点运用本发明方法所得到的目标航迹具有更小的距离误差， 所以本发明方法在实际应用中具有较好的优越性和实用性。

显然，本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和 范围。这样，倘若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围之内， 则本发明也意图包含这些改动和变型在内。