外辐射源雷达中一种基于随机有限集的多目标跟踪方法

摘要

本发明属于外辐射源雷达多目标跟踪技术领域，特别涉及外辐射源雷达中一种基于随机有限集的多目标跟踪方法。该外辐射源雷达中一种基于随机有限集的多目标跟踪方法包括以下步骤：利用外辐射源雷达接收数字电视信号，对外辐射源雷达接收的数字电视信号依次进行解调、重编码，重构出各个发射站的直达波信号；利用每个发射站的直达波信号对外辐射源雷达接收的数字电视信号进行自适应杂波相消处理，得到杂波相消后信号；根据所述杂波相消后信号、以及重构出的各个发射站的直达波信号，得到各个目标的观测量；采用概率假设密度粒子滤波方法对外辐射源雷达矩形观测区域内的每个目标实施定位跟踪，每个目标的定位跟踪方法为到达时间定位方法。

说明书

技术领域

本发明属于外辐射源雷达多目标跟踪技术领域，特别涉及外辐射源雷达 中一种基于随机有限集的多目标跟踪方法。本发明可在杂波背景下，实现对 多个目标的精确跟踪，且运算耗时较短。

背景技术

随着现代电子技术的发展，传统的有源雷达在电子干扰、反辐射导弹、 超低空突防及隐身技术等威胁中面临严重的生存危机。为了应对现代电子战 争中的四大威胁，外辐射源雷达成为雷达领域的研究热点。外辐射源雷达本 身不发射电磁波，而是利用已存在的民用信号(如FM调频信号、GMS(Global  System for Mobile communication，全球移动通讯系统)、电视信号等)作为第 三方照射源，对目标进行探测和定位，具有生存能力强、反隐身、抗低空突 防及成本低等优势。但是由于民用信号的发射功率较低，目标回波通常被直 达波和多径覆盖，因此实现目标探测和跟踪前需要先对回波信号进行处理。

同时对多个目标进行跟踪是雷达的重要职能之一，利用外辐射源雷达实 现多目标跟踪又称为被动多传感器多目标跟踪，相较于常规多目标跟踪，其 面临的问题更为复杂。一方面要在包含杂波的观测中配对目标与观测量；另 一方面，还要区分不同发射站(传感器)所对应的观测量；同时目标的数目 和观测量的数目是不确定而且时变的。传统的处理方法是利用数据关联算法 实现对目标的跟踪，主要算法有联合概率数据关联方法(JPDA)和多假设跟 踪方法(MHT)等。但联合概率数据关联方法只能处理目标个数固定且已知 的情况，多假设跟踪方法中的假设会随时间呈现指数型增长，出现组合爆炸。

20世纪九十年代，基于随机有限集(RFS)理论，Mahler提出了不需要 进行数据关联的概率假设密度(PHD)滤波，而概率假设密度滤波方法通过粒 子滤波(PF)实现，将概率假设密度滤波方法由理论层面带入工程实践中。 在外辐射源雷达中，利用(PF‐PHD)概率假设密度粒子滤波算法实现多目标 跟踪，虽然可以避免数据关联，但需要大量粒子以保证跟踪性能，计算量太 大。

发明内容

本发明的目的在于提出外辐射源雷达中一种基于随机有限集的多目标跟 踪方法。本发明以数字电视信号作为外辐射源雷达的第三方照射源，利用数 字电视信号的特性，获取纯净直达波，从而提高杂波相消性能；在概率假设 密度粒子滤波方法中，通过改进存活粒子状态的预测，优化新生粒子的分布， 实现用少量粒子对多个目标进行高精度跟踪。

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

外辐射源雷达中一种基于随机有限集的多目标跟踪方法包括以下步骤：

步骤1，利用外辐射源雷达接收数字电视信号，对外辐射源雷达接收的数 字电视信号依次进行解调、重编码，重构出各个发射站的直达波信号；

步骤2，利用每个发射站的直达波信号对外辐射源雷达接收的数字电视信 号进行自适应杂波相消处理，得到杂波相消后信号；根据所述杂波相消后信 号、以及重构出的各个发射站的直达波信号，得到各个目标的观测量；

步骤3，采用概率假设密度粒子滤波方法对外辐射源雷达矩形观测区域内 的每个目标实施定位跟踪，每个目标的定位跟踪方法为到达时间定位方法。

本发明的特点和进一步改进在于：

所述步骤2的具体子步骤为：

(2.1)利用每个发射站的直达波信号对外辐射源雷达接收的数字电视信 号进行自适应杂波相消处理，得到杂波相消后信号；

(2.2)对所述杂波相消后信号、以及重构出的各个发射站的直达波信号 进行距离—多普勒二维相关运算，得到各个目标的观测量；每个目标的观测 量包括：对应目标的双基地距离和以及对应目标的速度。

所述步骤3具体包括以下子步骤：

(3.1)设置参数k，k＝0,1,2...；当k＝0时，在k时刻用L₀个粒子表示目 标状态的先验分布，所述L₀个粒子表示为表 示0时刻所述L₀个粒子中第i₀个粒子的权值，表示所述L₀个粒子中第i₀个 粒子在0时刻的粒子状态；

(3.2)当k≥1时，将k-1时刻的粒子集为L_k-1为k-1时刻 的粒子数，i_k-1取1至L_k-1，为k-1时刻第i_k-1个粒子的权值，为k-1时 刻第i_k-1个粒子的粒子状态；在k-1时刻每个粒子处产生一组西格玛点，k-1时 刻每个粒子处产生的一组西格玛点包括至以及至n为粒子状态的维数；

将k-1时刻每个粒子处产生的一组西格玛点带入目标状态转移方程中，得 到k时刻的一步预测状态值j取0至2n；对至进行 加权求和，得到k时刻粒子状态的初步一步预测值计算出k时刻的 初步一步预测状态协方差矩阵

将代入目标观测方程中，得到k时刻的一步预测观测值 所述目标观测方程为z_k＝h(x_k,ν_k)，其中，ν_k为观测噪声，h(x_k,ν_k)表 示目标的观测函数；对至进行加权求和，得到k时刻观测量的初 步一步预测值

在k时刻目标的观测集Z_k中，选取出与具有最近欧式距离的对应目 标观测量z_k，k时刻目标的观测集Z_k为：由k时刻每个目标的观测量组成的 集合；用z_k对粒子状态的初步一步预测值和的初步一步预测状态协 方差矩阵做进一步修正，得到初步一步预测均值及的初步一步 预测状态协方差矩阵在均值为协方差为的高斯分布上采样， 得到k时刻第i_k-1个粒子的粒子状态的一步预测值

根据的概率密度q_k，并按照下式计算出所对应的权值

${\omega }_{k|k-1}^{\left(i_{k-1}\right)}={\omega }_{k-1}^{\left(i_{k-1}\right)}{p}_{s,k}{f}_{k|k-1}\left(x_{k|k-1}^{\left(i_{k-1}\right)}\right|x_{k|k-1}^{\left(i_{k-1}\right)})/q_k{,i}_{k-1}=\mathrm{1,2},...,L_{k-1}$

其中，p_s,k为k-1时刻的目标在k时刻仍然存在的概率,为目 标状态转移概率密度函数；

(3.3)当k≥1时，得出k时刻每个发射站对应的每个目标的双基地双基 地椭圆，在k时刻每个发射站对应的每个目标的双基地椭圆曲线上，每个点对 应的双基地距离和为一固定值，每个点对应的双基地距离和指：每个点与对 应发射站的距离、以及对应点与外辐射源雷达的距离之和；根据k时刻每个发 射站对应的每个目标的双基地椭圆曲线与外辐射源雷达的矩形观测区域边缘 的交点，分布k时刻新生粒子，并计算k时刻每个新生粒子的权值，k时刻新 生粒子的个数表示为J_k，k时刻第i_new个新生粒子的粒子状态表示为k 时刻第i_new个新生粒子的权值表示为i_new取L_k-1+1至L_k-1+J_k；

(3.4)当k≥1时，根据子步骤(3.2)和子步骤(3.3)，得到预测的k时 刻的粒子集其中，表示预测的k时刻的粒子集中第i'_k个粒子的粒子状态，表示预测的k时刻的粒子集中第i'_k个粒子的权值，i'_k取1至L_k-1+J_k；利用k时刻目标的观测集Z_k，对预测的k时刻的粒子集中第 i'_k个粒子的权值进行修正，得出预测的k时刻的粒子集中第i'_k个粒子的 修正后权值

(3.5)按照预测的k时刻的粒子集中第i'_k个粒子的修正后权值对 预测的k时刻的粒子集进行重采样，得出k时刻的粒子集；对k时刻的粒子集 进行聚类分析，得到每个目标k时刻的状态。

在子步骤(3.3)中，计算k时刻每个发射站对应的每个目标的双基地椭 圆曲线与外辐射源雷达的矩形观测区域边缘的交点的位置；根据k时刻每个发 射站对应的每个目标的双基地椭圆曲线与外辐射源雷达的矩形观测区域边缘 的交点的位置，在外辐射源雷达的矩形观测区域边缘中的上边缘、下边缘、 左边缘或右边缘中，找出任一个满足以下条件的边缘：与k时刻每个发射站对 应的每个目标的双基地椭圆曲线存在交点；将找出的满足以上条件的边缘记 为：k时刻目标的新生粒子分布边缘；

在k时刻目标的新生粒子分布边缘、以及k时刻每个发射站对应的每个目 标的双基地椭圆曲线的所有交点中，找出相距最近的l个交点，这l个交点与l 个发射站一一对应；得出所述相距最近的l个交点的位置的均值，将所述相距 最近的l个交点的位置的均值记为：k时刻共同交点；根据目标观测方程，计 算k时刻共同交点对应的目标速度；

在k时刻共同交点处采样J_k个点，得出k时刻J_k个新生粒子、k时刻J_k个新生粒子的粒子状态；按下式计算k时刻第i_new个新生粒子的权值

$\left(\begin{array}{c}{\omega }_{k|k-1}^{\left(i_{\mathrm{new}}\right)}=Q_x·q_x\left(x_k^{\left(i_{\mathrm{new}}\right)}\right|z_k)·Q_y·q_y\left(y_k^{\left(i_{\mathrm{new}}\right)}\right|z_k)·\\ Q_{\stackrel{·}{x}}·q_{\stackrel{·}{x}}\left({\stackrel{·}{x}}_k^{\left(i_{\mathrm{new}}\right)}\right|z_k)·Q_{\stackrel{·}{y}}·q_{\stackrel{·}{y}}\left({\stackrel{·}{y}}_k^{\left(i_{\mathrm{new}}\right)}\right|z_k)\end{array}\right)$

其中，i_new取L_k-1+1至L_k-1+J_k，为k时刻第i_new个新生粒子的位置 在x方向上的概率密度函数，为k时刻第i_new个新生粒子的位置在 y方向上的概率密度函数，为k时刻第i_new个新生粒子的速度在x 方向上的概率密度函数，为k时刻第i_new个新生粒子的速度在y方 向上的概率密度函数；当k时刻目标的新生粒子分布边缘为外辐射源雷达的矩 形观测区域边缘中的上边缘或下边缘时，当k时 刻目标的新生粒子分布边缘为外辐射源雷达的矩形观测区域边缘中的左边缘 或右边缘时，$Q_x=Q_{\stackrel{·}{x}}=2,Q_y=Q_{\stackrel{·}{y}}=1.$

在子步骤(3.4)中，首先由下式计算处的检测概率

$p_{d,k}\left(x_{k|k-1}^{\left({i^{\prime }}_k\right)}\right)=Q[\sqrt{2\mathrm{SNR}\left(x_{k|k-1}^{\left({i^{\prime }}_k\right)}\right)},\sqrt{2\ln (1/p_{\mathrm{FA}})}]$

其中，表示处的信噪比，p_FA为设定的虚警概率，Q[·]表示 marcum Q函数，$Q[a,b]=\underset{b}{\overset{\infty }{\int }}\mathrm{xexp}(-\frac{x^2+a^2}{2})I_0\left(\mathrm{ax}\right)\mathrm{dx},$I₀表示零阶贝塞尔函数；

然后，按照下式计算k时刻杂波平均数λ_k和杂波概率密度c_k：

${\lambda }_k=\frac{L_1}{l_1}·\frac{L_2}{l_2}·p_{\mathrm{FA}},c_k=1/\left(L_1{L}_2\right)$

其中，L₁为外辐射源雷达接收信号的距离范围，L₂为外辐射源雷达接收信号 的多普勒频率范围，l₁外辐射源雷达接收信号的距离分辨率，l₂为外辐射源雷 达接收信号的多普勒分辨率；

利用z_k计算似然函数

$g\left(z_k\right|x_{k|k-1}^{\left({i^{\prime }}_k\right)})=\frac{1}{2\pi \sqrt{R_v}}\exp (-\frac{1}{2}(z_k-h(x_{k|k-1}^{\left({i^{\prime }}_k\right)},v_k))R_v^{-1}{(z_k-h(x_{k|k-1}^{\left({i^{\prime }}_k\right)},v_k))}^T)$

其中，R_ν是观测噪声ν_k的协方差矩阵，上标‐1表示矩阵的逆，上标T表示矩 阵或向量的转置；

按照下式计算出预测的k时刻的粒子集中第i'_k个粒子的修正后权值

${\omega }_k^{\left({i^{\prime }}_k\right)}={\omega }_{k|k-1}^{\left({i^{\prime }}_k\right)}·[1-p_{d,k}\left(x_{k|k-1}^{\left({i^{\prime }}_k\right)}\right)+\underset{z_k\in Z_k}{\Sigma }\frac{{\psi }_{k,z}\left(x_{k|k-1}^{\left({i^{\prime }}_k\right)}\right)}{{\lambda }_k{c}_k+<{\omega }_{k|k-1},{\psi }_{k,z}>}]$

${\psi }_{k,z}\left(x_{k|k-1}^{\left({i^{\prime }}_k\right)}\right)=p_{d,k}\left(x_{k|k-1}^{\left({i^{\prime }}_k\right)}\right)·g\left(z_k\right|x_{k|k-1}^{\left({i^{\prime }}_k\right)})$

$<{\omega }_{k|k-1},{\psi }_{k,z}>=\underset{{i^{\prime }}_k=1}{\overset{L_{k-1}+J_k}{\Sigma }}{\psi }_{k,z}\left(x_{k|k-1}^{\left({i^{\prime }}_k\right)}\right)·{\omega }_{k|k-1}^{\left({i^{\prime }}_k\right)}$

其中，i'_k取1至L_k-1+J_k。

在子步骤(3.5)中，首先根据预测的k时刻的粒子集中第i'_k个粒子的修 正后权值计算k时刻的估计目标数

针对每个目标固定采样m个粒子，对预测的k时刻的粒子集进行重采样， 得到k时刻的粒子集i_k取1至L_k，表示k时刻 的粒子集中第i_k个粒子的粒子状态，表示k时刻的粒子集中第i_k个粒子权 值。

在步骤3之后，针对每个目标0时刻至k时刻的状态，进行航迹关联， 剔除虚假航迹，得出每个目标的航迹。

本发明的有益效果为：

第一，本发明通过对接收到的数字电视信号进行解调、重编码重构出各 个发射站的直达波信号，与借助辅助天线及滤波算法所获取的直达波相比， 其纯净度更高。

第二，本发明通过数字电视信号中的信息传输参数序列(TPS)区分接收 信号中对应于各个发射站的回波，用每个发射站的直达波与其相应的回波做 距离—多普勒二维相关，得到对应每个发射站的观测量，从而不用在后续处 理中进行观测量与发射站的配对。

第三，本发明在预测每个时刻存活粒子状态时，利用观测量对每一个存 活粒子进行无迹卡尔曼滤波，得到每个时刻粒子的一步预测状态和权值。由 于利用了观测值，预测的粒子状态更接近于真值，从而减轻了粒子退化现象。

第四，本发明在预测每个时刻新生粒子的状态时，将新生粒子分布在目 标、接收站与各个发射站构成的双基地椭圆与观测区域边缘的共同交点处， 缩小了新生粒子分布范围，从而减少了跟踪所需粒子数。

附图说明

图1为本发明的使用场景示意图；

图2为本发明应用的多目标到达时间定位原理示意图；

图3为本发明中对每个目标实施定位跟踪的步骤框图；

图4a为仿真实验中每个发射站对应的对应目标的双基地椭圆曲线与外辐 射源雷达的矩形观测区域边缘的交点示意图；

图4b为图4a的局部放大图；

图5a为改进预测粒子状态方法之前每时刻估计目标数、本发明得出的每 时刻估计目标数与真实目标数的对比图；

图5b为改进预测粒子状态方法之前与本发明得出的方位估计误差的对比 图；

图5c为改进预测粒子状态方法之前与本发明得出的速度估计误差的对比 图。

具体实施方式

下面结合附图对本发明作进一步说明：

步骤1，参照图1，为本发明的使用场景示意图。本发明实施例中，利用 外辐射源雷达(即图1中的接收站)接收数字电视信号，所述数字电视信号 包括各个发射站的直达波信号(图1中，3个发射站直接发送到接收站的数字 电视信号)、各个发射站的多径信号以及各个发射站的目标回波信号。本发明 实施例中，发射站的个数为l。

对外辐射源雷达接收的数字电视信号依次进行解调、重编码，重构出各 个发射站的直达波信号。

步骤2，利用每个发射站的直达波信号对外辐射源雷达接收的数字电视信 号进行自适应杂波相消处理，得到杂波相消后信号；根据所述杂波相消后信 号、以及重构出的各个发射站的直达波信号，得到各个目标的观测量。

其具体子步骤为：

(2.1)利用每个发射站的直达波信号对外辐射源雷达接收的数字电视信 号进行自适应杂波相消处理，消除外辐射源雷达接收的数字电视信号中各个 发射站的直达波信号和各个发射站的多径信号，得到杂波相消后信号。

(2.2)根据所述杂波相消后信号、以及重构出的各个发射站的直达波信 号，得到各个目标的观测量。具体地说，对所述杂波相消后信号、以及重构 出的各个发射站的直达波信号进行距离—多普勒二维相关运算，得到各个目 标的观测量。

本发明实施例中，每个目标的观测量包括：对应目标的双基地距离和以 及对应目标的速度(多普勒速度)。对应目标的双基地距离和包括：对应目标 相对于每个发射站的双基地距离和，对应目标相对于第l'个发射站的双基地距 离和为：对应目标与第l'个发射站的距离、以及对应目标与外辐射源雷达的距 离的和，l'取1至l。

步骤3，采用概率假设密度粒子滤波器(PF‐PHD)对外辐射源雷达矩形观测 区域内的每个目标实施定位跟踪，定位方法采用到达时间(TOA)定位方法； 参照图2，为本发明应用的多目标到达时间定位原理示意图。参照图3，为本 发明中对每个目标实施定位跟踪的步骤框图。

具体地，步骤3具体包括以下子步骤：

(3.1)初始化概率假设密度粒子滤波器(PF‐PHD)中的粒子状态和权值。 本发明实施例中，粒子状态由粒子的方位状态和速度状态组成。

具体地说，设置时刻参数k，k＝0,1,2...；

当k＝0时，在k时刻用L₀个粒子表示目标状态的先验分布，所述L₀个粒子 (即0时刻粒子集)表示为L₀为大于1的自然数，i₀＝1,2,…,L₀； 表示0时刻所述L₀个粒子中第i₀个粒子的权值，与估计目标数和粒子 数目有关，表示0时刻设定的目标的个数；表示所述L₀个粒子中第i₀个粒子在0时刻的粒子状态。

(3.2)当k≥1时，基于贝叶斯滤波原理，利用k时刻目标的观测集Z_k， 通过对k时刻存活粒子在k-1时刻的粒子状态进行无迹卡尔曼滤波，得到k时 刻存活粒子的一步预测状态和权值。k-1时刻的粒子集为L_k-1为k-1时刻的粒子数，i_k-1取1至L_k-1，为k-1时刻第i_k-1个粒子的权值，为k-1时刻第i_k-1个粒子的粒子状态。k时刻目标的观测集Z_k为：由k时刻每个 目标的观测量组成的集合。

具体地说，在子步骤(3.2)中，首先在k-1时刻每个粒子处产生一组西 格玛点(Sigma Point)n为粒子状态的维数，本发明实 施例中，n＝4。

当j＝0时，${\chi }_{k-1}^{(i_{k-1},j)}=x_{k-1}^{\left(i_{k-1}\right)},{\omega }_{k-1}^{(i_{k-1},j)}=\lambda /(n+\lambda ),$λ为尺度因子，

λ＝α²(n+ε₀)-n,(α<1,ε₀＝0)

其中，α为小于1的常数，ε₀＝0。

当j取1至n时，${\chi }_{k-1}^{(i_{k-1},j)}=x_{k-1}^{\left(i_{k-1}\right)}+\sqrt{(n+\lambda )P_{k-1}},{\omega }_{k-1}^{(i_{k-1},j)}=\lambda /2(n+\lambda ),$P_k-1为 k-1时刻第i_k-1个粒子的状态协方差矩阵。

当j取n+1至2n时，${\chi }_{k-1}^{(i_{k-1},j)}=x_{k-1}^{\left(i_{k-1}\right)}-\sqrt{(n+\lambda )P_{k-1}},{\omega }_{k-1}^{(i_{k-1},j)}=\lambda /2(n+\lambda ).$

然后，将k-1时刻每个粒子处产生的一组西格玛点带入目标状态转移方程 中，得到k时刻的一步预测状态值所述目标状态转移方程为： x_k＝f(x_k-1,u_k)，x_k-1为k-1时刻目标的状态，x_k为k时刻目标的状态，u_k为过 程噪声，f(x_k-1,u_k)表示目标的状态转移函数；对至进行加权求 和，得到k时刻粒子状态的初步一步预测值计算出k时刻的初步一 步预测状态协方差矩阵具体操作过程如下：

${\chi }_{k|k-1}^{(i_{k-1},j)}=f({\chi }_{k-1}^{(i_{k-1},j)},u_k)$

${\bar{x}}_{k|k-1}^{\left(i_{k-1}\right)}=\underset{j=0}{\overset{2n}{\Sigma }}{\omega }_{k-1}^{(i_{k-1},j)}{\chi }_{k|k-1}^{(i_{k-1},j)}$

$P_{k|k-1}^{\left(i_{k-1}\right)}=\underset{j=0}{\overset{2n}{\Sigma }}{\omega }_{k-1}^{(i_{k-1},j)}[{\chi }_{k-1}^{(i_{k-1},j)}-{\bar{x}}_{k|k-1}^{\left(i_{k-1}\right)}]{[{\chi }_{k-1}^{(i_{k-1},j)}-{\bar{x}}_{k|k-1}^{\left(i_{k-1}\right)}]}^T$

其中，上标T表示向量或矩阵的转置。

然后，将代入目标观测方程中，得到k时刻的一步预测观测值 所述目标观测方程为z_k＝h(x_k,ν_k)，其中，ν_k为观测噪声，h(x_k,ν_k)表 示目标的观测函数。对至进行加权求和，得到k时刻观测量的初 步一步预测值具体操作过程如下：

${\bar{z}}_{k|k-1}^{(i_{k-1},j)}=h({\chi }_{k|k-1}^{(i_{k-1},j)},v_k),{\bar{\eta }}_{k|k-1}^{\left(i_{k-1}\right)}=\underset{j=0}{\overset{2n}{\Sigma }}{\omega }_{k-1}^{(i_{k-1},j)}{\bar{z}}_{k|k-1}^{(i_{k-1},j)}$

然后，在k时刻目标的观测集Z_k中，选取出与具有最近欧式距离的 目标观测量z_k，用z_k对粒子状态的初步一步预测值和的初步一步预 测状态协方差矩阵做进一步修正，得到初步一步预测均值及的 初步一步预测状态协方差矩阵在均值为协方差为的高斯分 布上采样，得到k时刻第i_k-1个粒子的粒子状态的一步预测值具体操作 过程如下：

${\tilde{x}}_{k|k-1}^{\left(i_{k-1}\right)}={\bar{x}}_{k|k-1}^{\left(i_{k-1}\right)}+K_k^{\left(i_{k-1}\right)}(z_k-{\eta }_{k|k-1}^{\left(i_{k-1}\right)}),K_k^{\left(i_{k-1}\right)}=G_k^{\left(i_{k-1}\right)}{\left[S_k^{\left(i_{k-1}\right)}\right]}^{-1}$

$P_k^{\left(i_{k-1}\right)}=P_{k|k-1}^{\left(i_{k-1}\right)}-G_k^{\left(i_{k-1}\right)}{\left[S_k^{\left(i_{k-1}\right)}\right]}^{-1}{\left(G_k^{\left(i_{k-1}\right)}\right)}^T$

$S_k^{\left(i_{k-1}\right)}=\underset{j=0}{\overset{2n}{\Sigma }}{\omega }_{k-1}^{(i_{k-1},j)}({\bar{z}}_{k-1}^{(i_{k-1},j)}-{\bar{\eta }}_{k|k-1}^{\left(i_{k-1}\right)}){({\bar{z}}_{k|k-1}^{(i_{k-1},j)}-{\bar{\eta }}_{k|k-1}^{\left(i_{k-1}\right)})}^T$

$G_k^{\left(i_{k-1}\right)}=\underset{j=0}{\overset{2n}{\Sigma }}{\omega }_{k-1}^{(i_{k-1},j)}({\chi }_{k-1}^{(i_{k-1},j)}-{\bar{x}}_{k|k-1}^{\left(i_{k-1}\right)}){({\bar{z}}_{k|k-1}^{(i_{k-1},j)}-{\bar{\eta }}_{k|k-1}^{\left(i_{k-1}\right)})}^T$

$x_{k|k-1}^{\left(i_{k-1}\right)}~N({\tilde{x}}_{k|k-1}^{\left(i_{k-1}\right)},P_k^{\left(i_{k-1}\right)})$

其中，上标-1表示矩阵的逆，上标T表示矩阵或向量的转置，表 示均值为协方差为的高斯分布。

然后，根据的概率密度q_k，并按照下式计算出所对应的权值

${\omega }_{k|k-1}^{\left(i_{k-1}\right)}={\omega }_{k-1}^{\left(i_{k-1}\right)}{p}_{s,k}{f}_{k|k-1}\left(x_{k|k-1}^{\left(i_{k-1}\right)}\right|x_{k-1}^{\left(i_{k-1}\right)})/q_k{,i}_{k-1}=\mathrm{1,2},...,L_{k-1}$

其中，p_s,k为k-1时刻的目标在k时刻仍然存在的概率,为目 标状态转移概率密度函数。

(3.3)当k≥1时，得出k时刻每个发射站对应的每个目标的双基地双基 地椭圆，在k时刻每个发射站对应的每个目标的双基地椭圆曲线上，每个点对 应的双基地距离和为一固定值，每个点对应的双基地距离和指：每个点与对 应发射站的距离、以及对应点与外辐射源雷达的距离之和。

根据k时刻l个发射站对应的每个目标的双基地椭圆曲线与外辐射源雷达 的矩形观测区域边缘的交点，分布k时刻新生粒子，并计算k时刻每个新生粒 子的权值。k时刻新生粒子的个数表示为J_k，k时刻第i_new个新生粒子的粒子 状态表示为k时刻第i_new个新生粒子的权值表示为i_new取 L_k-1+1至L_k-1+J_k。

具体地说，在子步骤(3.3)中，首先，计算k时刻每个发射站对应的每 个目标的双基地椭圆曲线与外辐射源雷达的矩形观测区域边缘的交点的位 置。

然后，根据k时刻每个发射站对应的每个目标的双基地椭圆曲线与外辐射 源雷达的矩形观测区域边缘的交点的位置，在外辐射源雷达的矩形观测区域 边缘中的上边缘、下边缘、左边缘或右边缘中，找出任一个满足以下条件的 边缘：与k时刻每个发射站对应的每个目标的双基地椭圆曲线存在交点；将找 出的满足以上条件的边缘记为：k时刻目标的新生粒子分布边缘。

在k时刻目标的新生粒子分布边缘、以及k时刻每个发射站对应的每个目 标的双基地椭圆曲线的所有交点中，找出相距最近的l个交点，这l个交点与l 个发射站一一对应。得出所述相距最近的l个交点的位置的均值，将所述相距 最近的l个交点的位置的均值记为：k时刻共同交点；根据目标观测方程，计 算k时刻共同交点对应的目标速度(k时刻共同交点位置代替目标观测方程的 x_k)。

然后，在k时刻共同交点处采样J_k个点，得出k时刻J_k个新生粒子、k时 刻J_k个新生粒子的粒子状态。k时刻每个新生粒子的粒子状态包括：k时刻对 应新生粒子的方位状态、k时刻对应新生粒子的速度状态。具体地，在k时刻 共同交点的位置高斯分布中采样J_k个点，得出k时刻J_k个新生粒子的方位状 态，k时刻共同交点的位置高斯分布指：均值为k时刻共同交点的位置坐标(包 括横坐标和纵坐标)，标准差为双基地距离单元(外辐射源雷达的距离分辨率) 的高斯分布。在k时刻共同交点的速度高斯分布采样J_k个点，得出k时刻J_k个 新生粒子的速度状态，k时刻共同交点的速度高斯分布指：均值为k时刻共同 交点的速度(包括速度的横坐标以及速度的纵坐标)，标准差为外辐射源雷达 的多普勒分辨率的高斯分布。

然后，按下式计算k时刻第i_new个新生粒子的权值

$\left(\begin{array}{c}{\omega }_{k|k-1}^{\left(i_{\mathrm{new}}\right)}=Q_x·q_x\left(x_k^{\left(i_{\mathrm{new}}\right)}\right|z_k)·Q_y·q_y\left(y_k^{\left(i_{\mathrm{new}}\right)}\right|z_k)·\\ Q_{\stackrel{·}{x}}·q_{\stackrel{·}{x}}\left({\stackrel{·}{x}}_k^{\left(i_{\mathrm{new}}\right)}\right|z_k)·Q_{\stackrel{·}{y}}·q_{\stackrel{·}{y}}\left({\stackrel{·}{y}}_k^{\left(i_{\mathrm{new}}\right)}\right|z_k)\end{array}\right)$

其中，i_new取L_k-1+1至L_k-1+J_k，为k时刻第i_new个新生粒子的位置 在x方向上的概率密度函数，为k时刻第i_new个新生粒子的位置在 y方向上的概率密度函数，为k时刻第i_new个新生粒子的速度在x 方向上的概率密度函数，为k时刻第i_new个新生粒子的速度在y方 向上的概率密度函数。当k时刻目标的新生粒子分布边缘为外辐射源雷达的矩 形观测区域边缘中的上边缘或下边缘时，当k时 刻目标的新生粒子分布边缘为外辐射源雷达的矩形观测区域边缘中的左边缘 或右边缘时，$Q_x=Q_{\stackrel{·}{x}}=2,Q_y=Q_{\stackrel{·}{y}}=1.$

(3.4)当k≥1时，根据子步骤(3.2)和子步骤(3.3)，得到预测的k时 刻的粒子集其中，表示预测的k时刻的粒子集中第i'_k个粒子的粒子状态，表示预测的k时刻的粒子集中第i'_k个粒子的权值，i'_k取1至L_k-1+J_k。利用k时刻目标的观测集Z_k，对预测的k时刻的粒子集中第 i'_k个粒子的权值进行修正，得出预测的k时刻的粒子集中第i'_k个粒子的 修正后权值

在已知处的信噪比的条件下，首先由下式计算处的检测概率

$p_{d,k}\left(x_{k|k-1}^{\left({i^{\prime }}_k\right)}\right)=Q[\sqrt{2\mathrm{SNR}\left(x_{k|k-1}^{\left({i^{\prime }}_k\right)}\right)},\sqrt{2\ln (1/p_{\mathrm{FA}})}]$

其中表示处的信噪比(每个发射站对应的处的信噪比，本 发明实施例中，每个发射站对应的处的信噪比相同)p_FA为设定的虚警概 率，Q[·]表示marcum Q函数，$Q[a,b]=\underset{b}{\overset{\infty }{\int }}\mathrm{xexp}(-\frac{x^2+a^2}{2})I_0\left(\mathrm{ax}\right)\mathrm{dx},$I₀表示零 阶贝塞尔函数。

然后，按照下式计算k时刻杂波平均数λ_k和杂波概率密度c_k：

${\lambda }_k=\frac{L_1}{l_1}·\frac{L_2}{l_2}·p_{\mathrm{FA}},c_k=1/\left(L_1{L}_2\right)$

其中，L₁为外辐射源雷达接收信号的距离范围，L₂为外辐射源雷达接收信号 的多普勒频率范围，l₁外辐射源雷达接收信号的距离分辨率，l₂为外辐射源雷 达接收信号的多普勒分辨率。

然后，利用z_k，计算似然函数z_k∈Z_k。

$g\left(z_k\right|x_{k|k-1}^{\left({i^{\prime }}_k\right)})=\frac{1}{2\pi \sqrt{R_v}}\exp (-\frac{1}{2}(z_k-h(x_{k|k-1}^{\left({i^{\prime }}_k\right)},v_k))R_v^{-1}{(z_k-h(x_{k|k-1}^{\left({i^{\prime }}_k\right)},v_k))}^T)$

其中，R_ν是观测噪声ν_k的协方差矩阵，上标‐1表示矩阵的逆，上标T表示矩 阵或向量的转置。

然后，按照下式计算出预测的k时刻的粒子集中第i'_k个粒子的修正后权 值

${\omega }_k^{\left({i^{\prime }}_k\right)}={\omega }_{k|k-1}^{\left({i^{\prime }}_k\right)}·[1-p_{d,k}\left(x_{k|k-1}^{\left({i^{\prime }}_k\right)}\right)+\underset{z_k\in Z_k}{\Sigma }\frac{{\psi }_{k,z}\left(x_{k|k-1}^{\left({i^{\prime }}_k\right)}\right)}{{\lambda }_k{c}_k+<{\omega }_{k|k-1},{\psi }_{k,z}>}]$

${\psi }_{k,z}\left(x_{k|k-1}^{\left({i^{\prime }}_k\right)}\right)=p_{d,k}\left(x_{k|k-1}^{\left({i^{\prime }}_k\right)}\right)·g\left(z_k\right|x_{k|k-1}^{\left({i^{\prime }}_k\right)})$

$<{\omega }_{k|k-1},{\psi }_{k,z}>=\underset{{i^{\prime }}_k=1}{\overset{L_{k-1}+J_k}{\Sigma }}{\psi }_{k,z}\left(x_{k|k-1}^{\left({i^{\prime }}_k\right)}\right)·{\omega }_{k|k-1}^{\left({i^{\prime }}_k\right)}$

其中，i'_k取1至L_k-1+J_k。

(3.5)按照预测的k时刻的粒子集中第i'_k个粒子的修正后权值对 预测的k时刻的粒子集进行重采样，得出k时刻的粒子集；对k时刻的粒子集 进行聚类分析，得到每个目标k时刻的状态。

具体地，在子步骤(3.5)中，首先根据预测的k时刻的粒子集中第i'_k个 粒子的修正后权值计算k时刻的估计目标数

针对每个目标固定采样m个粒子，对预测的k时刻的粒子集进行重采样， 得到k时刻的粒子集i_k取1至L_k，表示k时刻 的粒子集中第i_k个粒子的粒子状态，表示k时刻的粒子集中第i_k个粒子权 值，${\omega }_k^{\left(i_k\right)}={\bar{N}}_k/L_k.$

在步骤3之后，针对每个目标0时刻至k时刻的状态，进行航迹关联，剔 除虚假航迹，得出每个目标的航迹。

本发明的效果可以通过以下仿真实验进行验证。

1)仿真条件

本发明的仿真实验中，设置三个发射站与一个接收站(外辐射源雷达)， 本发明的仿真实验软件平台为MATLAB，操作系统为Win XP系统。设定的虚 警概率p_FA为10^-4，信噪比SNR为15dB。目标在外辐射源雷达的观测区域内匀 速运动，且每一帧最多有一个新生目标出现。目标存活概率p_s＝0.98，每个目 标固定分配30个粒子，新生目标分配20个粒子。

2)仿真结果分析

参照图4a，为仿真实验中k时刻每发射站对应的对应目标的双基地椭圆 曲线与外辐射源雷达的矩形观测区域边缘的交点示意图。参照图4b，为图4a 的局部放大图。图4b中，新生粒子分布在A点和B点的附近区域内。

参照图5a，为改进预测粒子状态方法之前每时刻估计目标数、本发明得 出的每时刻估计目标数与真实目标数的对比图，图5a中，横轴表示时间，纵 轴表示目标数。无标示直线代表真实目标数，以圆形标示的直线代表算法改 进前的估计目标数，以三角形标示的直线代表改进后算法的估计目标数。

对改进预测粒子状态方法之前与本发明得出的方位估计误差和速度估计 误差进行对比，采用最优子模式分配(Optimal Sub‐pattern Assignment,OSPA)作 为多目标跟踪性能评估标准，其中距离敏感性参数选择2，关联敏感性参数选 择5。参照图5b，为改进预测粒子状态方法之前与本发明得出的方位估计误 差的对比图。图5b中，横轴表示时间，纵轴表示方位估计误差，单位为m， 图5b中，以圆形标示的直线代表改进预测粒子状态方法之前得出的方位估计 误差，以三角形标示的直线代表本发明得出的方位估计误差。

参照图5c，为改进预测粒子状态方法之前与本发明得出的速度估计误差 的对比图。图5b中，横轴表示时间，纵轴表示速度估计误差，单位为m/s， 图5b中，以圆形标示的直线代表改进预测粒子状态方法之前得出的速度估计 误差，以三角形标示的直线代表本发明得出的速度估计误差。从图5a至图5c 可以看出，采用本发明之后，目标数目的估计精度及方位和速度估计精度都 得到了很大改善。此外，由于提高了预测存活粒子状态的准确度，优化了新 生粒子分布，所以减少了跟踪所需粒子数，降低了运算耗时。

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