CN104077498A - 一种结合目标角度的外辐射源雷达多目标跟踪方法 - Google Patents

Abstract

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

Description

一种结合目标角度的外辐射源雷达多目标跟踪方法

技术领域

本发明属于外辐射源雷达多目标跟踪技术领域，特别涉及一种结合目标角度的外辐射源雷达多目标跟踪方法。本发明利用外辐射源雷达系统中采用多发一收结构，将目标双基地距离、多普勒速度与目标角度信息相结合，同时对多个目标进行定位跟踪，本发明可在杂波背景下，实现对多个目标的准确跟踪，运算量小，跟踪精度高。

背景技术

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

利用外辐射源雷达实现多目标跟踪是雷达跟踪领域中的研究热点和难点。一方面多目标跟踪在民用和军事上都有重要应用；另一方面外辐射源雷达多目标跟踪要求能从包含杂波且数目时变的观测量中提取出目标个数不定的目标状态，为实现这一目的，需要配对目标与观测量、发射站与观测量。传统的处理方法基于数据关联，存在计算量随时间增长的问题。为避免数据关联，Mahler提出概率假设密度(PHD)滤波器，随后概率假设密度粒子滤波器(PF‐PHD)的提出将PHD算法由理论层面带入工程实践中。

虽然PF‐PHD可以避免数据关联，但需要大量粒子以保证跟踪性能，计算量太大。多种观测量的结合会提高跟踪精度，在外辐射源雷达跟踪中最常用的观测量是目标的双基地距离和多普勒速度，但只利用这两种观测量概率假设密度粒子滤波器(PF‐PHD)仍然需要大量粒子来保证跟踪性能，运算量大，实时性差。

发明内容

本发明的目的在于提出一种结合目标角度的外辐射源雷达多目标跟踪方法。本发明以数字电视信号作为外辐射源雷达的第三方照射源，从数字电视信号的信道估计参数中直接获取目标时延和多普勒频移，简化目标跟踪的前期信号处理过程；在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，并按照下式计算出所对应的权值

${\mathrm{\ω}}_{k|k-1}^{\left(i_{k-1}\right)}={\mathrm{\ω}}_{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-1时刻每个目标相对于外辐射源雷达的角度、以及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-1时刻目标的观测集Z_k-1，获取k-1时刻每个目标相对于外辐射源雷达的角度，将k-1时每个目标相对于外辐射源雷达的角度组成k-1时刻目标观测量角度集合θ_k-1，k-1时刻目标观测量角度集合θ_k-1中元素的个数表示为<θ_k-1>，k-1时刻目标观测量角度集合θ_k-1中第α_k-1个元素表示为α_k-1取1至<θ_k-1>；根据k时刻目标的观测集Z_k，获取k时刻每个目标相对于外辐射源雷达的角度，将k时刻每个目标相对于外辐射源雷达的角度组成k时刻目标观测量角度集合θ_k，k时刻目标观测量角度集合θ_k中元素的个数表示为<θ_k>，k时刻目标观测量角度集合θ_k中第α_k个元素表示为α_k取1至<θ_k>；

判断k时刻是否产生新生粒子，如果k时刻不产生新生粒子，则k时刻新生粒子的个数J_k为零；如果k时刻产生新生粒子，则得出k时刻产生的每个新生粒子的分布区域中心点方位角，k时刻产生的每个新生粒子的分布区域中心点方位角指：k时刻产生的对应新生粒子的位置的纵坐标与横坐标的比值的反正切值；

判断k时刻是否产生新生粒子的过程为：将α_k-1从1至<θ_k-1>进行遍历，α_k-1每取一个数值时，在k时刻目标观测量角度集合θ_k中找出满足第一新生粒子存在条件的元素，所述第一新生粒子存在条件为：k时刻目标观测量角度集合θ_k中的元素与的差值的绝对值不小于γ_θ，γ_θ为设定的角度门限；在将α_k-1从1至<θ_k-1>进行遍历的过程中，若k时刻目标观测量角度集合θ_k中不存在满足第一新生粒子存在条件的元素，则k时刻不产生新生粒子；反之，k时刻产生新生粒子；

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

得出k时刻目标的新生粒子分布边缘；得出k时刻目标的新生粒子分布边缘的过程为：当k时刻目标的候选新生粒子分布边缘的个数等于1时，k时刻目标的新生粒子分布边缘为k时刻目标的候选新生粒子分布边缘；当k时刻目标的候选新生粒子分布边缘的个数大于1时，在k时刻目标的所有候选新生粒子分布边缘中找出k时刻目标的新生粒子分布边缘；

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

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

当k时刻目标的候选新生粒子分布边缘的个数大于1时，找出k时刻目标的新生粒子分布边缘的过程为：判断k时刻是否产生新生粒子，如果k时刻不产生新生粒子，则没有k时刻目标的新生粒子分布边缘；如果k时刻产生新生粒子，针对k时刻的每个新生粒子，判断每个候选新生粒子分布边缘的k时刻共同交点的方位角是否满足第二新生粒子存在条件，所述第二新生粒子存在条件为：每个候选新生粒子分布边缘的k时刻共同交点的方位角与k时刻目标观测量角度集合θ_k中满足第一新生粒子存在条件的元素的差值的绝对值不大于设定的方位角门限值；如果每个候选新生粒子分布边缘的k时刻共同交点的方位角不满足第二新生粒子存在条件，则令k时刻新生粒子的个数J_k取0；如果至少一个候选新生粒子分布边缘满足第二新生粒子存在条件，则在满足第二新生粒子存在条件的候选新生粒子分布边缘中，选择其k时刻共同交点的方位角与k时刻的对应新生粒子的分布区域中心点方位角最接近的候选新生粒子分布边缘，则k时刻目标的新生粒子分布边缘为：在满足第二新生粒子存在条件的候选新生粒子分布边缘中选择的候选新生粒子分布边缘，此时，k时刻共同交点为：k时刻目标的新生粒子分布边缘的k时刻共同交点，k时刻共同交点对应的目标速度为：k时刻目标的新生粒子分布边缘的k时刻共同交点对应的目标速度；

在k时刻共同交点处采样J_k个点，得出k时刻J_k个新生粒子；k时刻每个新生粒子的粒子状态包括：k时刻对应新生粒子的方位状态、k时刻对应新生粒子的速度状态；

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

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

其中，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{\&CenterDot;}{x}}=2,Q_y=Q_{\stackrel{\&CenterDot;}{y}}=1.$

在子步骤(3.4)中，当k≥1时，当k≥1时，对预测的k时刻的粒子集中第i'_k个粒子的权值进行修正的过程包括以下子步骤：

a)设置迭代参数j，j＝1,2,...当j＝1时，执行子步骤b)；

b)根据以下公式计算第j个发射站在处的检测概率

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

其中，表示已知的第j个发射站对应的处的信噪比，p_FA为设定的虚警概率， $Q[a,b]=\underset{b}{\overset{\text{\&infin;}}{\text{\&Integral;}}}\mathrm{xexp}(-\frac{x^2+a^2}{2})I_0\left(\mathrm{ax}\right)\text{dx,}$ I₀表示零阶贝塞尔函数；

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

${\mathrm{\&lambda;}}_k=\frac{L_1}{l_1}\&CenterDot;\frac{L_2}{l_2}\&CenterDot;p_{\mathrm{FA}},$ c_k＝1/(L₁L₂)

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

d)在k时刻目标的观测集Z_k中，利用与第j个发射站对应的k时刻每个目标的观测量，组成k时刻第j个发射站的目标的观测集然后利用中的观测量计算似然函数

$g\left(z_k^{\left(j\right)}\right|x_{k|k-1}^{\left({i^{\&prime;}}_k\right)})=\frac{1}{2\mathrm{\&pi;}\sqrt{R_v}}\exp (-\frac{1}{2}(z_k^{\left(j\right)}-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_v是观测噪声ν_k的协方差矩阵，上标‐1表示矩阵的逆，上标T表示矩阵或向量的转置；

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

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

${\mathrm{\&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)\&CenterDot;g\left(z_k^{\left(j\right)}\right|x_{k|k-1}^{\left({i^{\&prime;}}_k\right)})$

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

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

然后，令j值自增1；

f)判断j与发射站的个数为l的大小关系，如果j<l，则令然后返回至子步骤b)，如果j＝l，则迭代终止，预测的k时刻的粒子集中第i'_k个粒子的修正后权值为

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

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

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

本发明的有益效果为：

第一，本发明通过数字电视信号的信道估计参数直接得到所需目标状态参数(包括时延和多普勒频率)，不需要用每个发射站的直达波与其相应的回波做距离—多普勒二维相关处理，从而简化了外辐射源雷达信号处理流程。

第二，本发明除利用目标的双基地距离和多普勒速度外，还结合目标的角度信息实现目标跟踪。在分布新生粒子前，先利用角度信息判断当前时刻是否有新目标出现，而不是在每一时刻都产生新生粒子，从而降低了运算量。

第三，本发明在分布新生粒子时，利用角度信息剔除了观测区域边缘上虚假的共同交点，缩小了新生粒子分布范围，进一步降低跟踪所需粒子数。

第四，本发明利用相对于每个发射站的观测量依次对预测权值进行更新，有利于剔除对应单个发射站的观测量中可能存在的杂波，从而提高跟踪精度。

附图说明

图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。对应目标相对于外辐射源雷达(接收站)的角度指：对应目标的外辐射源雷达连线方向与水平面的夹角(是否正确)，对应目标的外辐射源雷达连线方向为：对应目标与外辐射源雷达的连线方向。

下面说明对应目标相对于每个发射站的双基地距离和的计算过程、以及对应目标的速度的计算过程。首先利用数字电视信号的信道估计参数直接得到对应目标相对于每个发射站的时延、以及对应目标的多普勒频率。对于其中一个发射站来说，如果对应目标相对于对应发射站的时延为τ，对应目标的多普勒频率为f_d，每个发射站发射信号的波长为γ，则对应目标相对于发射站的时延为bτ/2，对应目标的多普勒速度为f_dγ/2。

步骤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时， ${\mathrm{\&chi;}}_{k-1}^{(i_{k-1},j)}=x_{k-1}^{\left(i_{k-1}\right)},{\mathrm{\&omega;}}_{k-1}^{(i_{k-1},j)}=\mathrm{\&lambda;}/(n+\mathrm{\&lambda;}),$ λ为尺度因子，

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

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

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

当j取n+1至2n时， ${\mathrm{\&chi;}}_{k-1}^{(i_{k-1},j)}=x_{k-1}^{\left(i_{k-1}\right)}-\sqrt{(n+\mathrm{\&lambda;})P_{k-1}},{\mathrm{\&omega;}}_{k-1}^{(i_{k-1},j)}=\mathrm{\&lambda;}/2(n+\mathrm{\&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时刻的初步一步预测状态协方差矩阵具体操作过程如下：

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

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

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

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

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

${\stackrel{\&OverBar;}{z}}_{k|k-1}^{(i_{k-1},j)}=h({\mathrm{\&chi;}}_{k|k-1}^{(i_{k-1},j)},v_k),{\stackrel{\&OverBar;}{\mathrm{\&eta;}}}_{k|k-1}^{\left(i_{k-1}\right)}=\underset{j=0}{\overset{2n}{\mathrm{\&Sigma;}}}{\mathrm{\&omega;}}_{k-1}^{(i_{k-1},j)}{\stackrel{\&OverBar;}{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)}={\stackrel{\&OverBar;}{x}}_{k|k-1}^{\left(i_{k-1}\right)}+K_k^{\left(i_{k-1}\right)}(z_k-{\mathrm{\&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}{\mathrm{\&Sigma;}}}{\mathrm{\&omega;}}_{k-1}^{(i_{k-1},j)}({\stackrel{\&OverBar;}{z}}_{k|k-1}^{(i_{k-1},j)}-{\stackrel{\&OverBar;}{\mathrm{\&eta;}}}_{k|k-1}^{\left(i_{k-1}\right)}){({\stackrel{\&OverBar;}{z}}_{k|k-1}^{(i_{k-1},j)}-{\stackrel{\&OverBar;}{\mathrm{\&eta;}}}_{k|k-1}^{\left(i_{k-1}\right)})}^T$

$G_k^{\left(i_{k-1}\right)}=\underset{j=0}{\overset{2n}{\mathrm{\&Sigma;}}}{\mathrm{\&omega;}}_{k-1}^{(i_{k-1},j)}({\mathrm{\&chi;}}_{k-1}^{(i_{k-1},j)}-{\stackrel{\&OverBar;}{x}}_{k|k-1}^{\left(i_{k-1}\right)}){({\stackrel{\&OverBar;}{z}}_{k|k-1}^{(i_{k-1},j)}-{\stackrel{\&OverBar;}{\mathrm{\&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，按照下式计算出一步预测状态为的粒子的权值

${\mathrm{\&omega;}}_{k|k-1}^{\left(i_{k-1}\right)}={\mathrm{\&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},\&CenterDot;\&CenterDot;\&CenterDot;,L_{k-1}$

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

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

根据k-1时刻每个目标相对于外辐射源雷达的角度、以及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.3)中，首先，根据k-1时刻目标的观测集Z_k-1，获取k-1时刻每个目标相对于外辐射源雷达的角度，将k-1时每个目标相对于外辐射源雷达的角度组成k-1时刻目标观测量角度集合θ_k-1，k-1时刻目标观测量角度集合θ_k-1中元素的个数表示为<θ_k-1>，k-1时刻目标观测量角度集合θ_k-1中第α_k-1个元素表示为α_k-1取1至<θ_k-1>。根据k时刻目标的观测集Z_k，获取k时刻每个目标相对于外辐射源雷达的角度，将k时刻每个目标相对于外辐射源雷达的角度组成k时刻目标观测量角度集合θ_k，k时刻目标观测量角度集合θ_k中元素的个数表示为<θ_k>，k时刻目标观测量角度集合θ_k中第α_k个元素表示为α_k取1至<θ_k>。

判断k时刻是否产生新生粒子，如果k时刻不产生新生粒子，则k时刻新生粒子的个数J_k为零。如果k时刻产生新生粒子，则得出k时刻产生的每个新生粒子的分布区域中心点方位角，k时刻产生的每个新生粒子的分布区域中心点方位角指：k时刻产生的对应新生粒子的位置的纵坐标与横坐标的比值的反正切值。

判断k时刻是否产生新生粒子的过程为：将α_k-1从1至<θ_k-1>进行遍历，α_k-1每取一个数值时，在k时刻目标观测量角度集合θ_k中找出满足第一新生粒子存在条件的元素，所述第一新生粒子存在条件为：k时刻目标观测量角度集合θ_k中的元素与的差值的绝对值不小于γ_θ，γ_θ为设定的角度门限；在将α_k-1从1至<θ_k-1>进行遍历的过程中，若k时刻目标观测量角度集合θ_k中不存在满足第一新生粒子存在条件的元素，则k时刻不产生新生粒子；反之，k时刻产生新生粒子。

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

然后得出k时刻目标的新生粒子分布边缘。得出k时刻目标的新生粒子分布边缘的过程为：当k时刻目标的候选新生粒子分布边缘的个数等于1时，k时刻目标的新生粒子分布边缘为k时刻目标的候选新生粒子分布边缘。当k时刻目标的候选新生粒子分布边缘的个数大于1时，在k时刻目标的所有候选新生粒子分布边缘中找出k时刻目标的新生粒子分布边缘。

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

当k时刻目标的候选新生粒子分布边缘的个数大于1时，在k时刻目标的每个候选新生粒子分布边缘、以及k时刻每个发射站对应的每个目标的双基地椭圆曲线的所有交点中，找出相距最近的l个交点，这l个交点与l个发射站一一对应。得出所述相距最近的l个交点的位置的均值，将所述相距最近的l个交点的位置的均值为：对应候选新生粒子分布边缘的k时刻共同交点；根据目标观测方程，针对每个候选新生粒子分布边缘的k时刻共同交点，得出对应候选新生粒子分布边缘的k时刻共同交点对应的目标速度。例如，k时刻目标的候选新生粒子分布边缘的个数表示为W_k，W_k为小于或等于4的自然数。在k时刻目标的第w_k个候选新生粒子分布边缘、以及k时刻每个发射站对应的每个目标的双基地椭圆曲线的所有交点中，找出相距最近的l个交点，w_k取1至W_k。得出所述相距最近的l个交点的位置的均值，将所述相距最近的l个交点的位置的均值为：第w_k个候选新生粒子分布边缘的k时刻共同交点。

当k时刻目标的候选新生粒子分布边缘的个数大于1时，找出k时刻目标的新生粒子分布边缘的过程为：判断k时刻是否产生新生粒子，如果k时刻不产生新生粒子，则没有k时刻目标的新生粒子分布边缘。如果k时刻产生新生粒子，针对k时刻的每个新生粒子，判断每个候选新生粒子分布边缘的k时刻共同交点的方位角是否满足第二新生粒子存在条件，所述第二新生粒子存在条件为：每个候选新生粒子分布边缘的k时刻共同交点的方位角(每个候选新生粒子分布边缘的k时刻共同交点的位置的纵坐标与横坐标的比值的反正切值)与k时刻目标观测量角度集合θ_k中满足第一新生粒子存在条件的元素的差值的绝对值不大于设定的方位角门限值；如果每个候选新生粒子分布边缘的k时刻共同交点的方位角不满足第二新生粒子存在条件，则令k时刻新生粒子的个数J_k取0。如果至少一个候选新生粒子分布边缘满足第二新生粒子存在条件，则在满足第二新生粒子存在条件的候选新生粒子分布边缘中，选择其k时刻共同交点的方位角与k时刻的对应新生粒子的分布区域中心点方位角最接近的候选新生粒子分布边缘，则k时刻目标的新生粒子分布边缘为：在满足第二新生粒子存在条件的候选新生粒子分布边缘中选择的候选新生粒子分布边缘，此时，k时刻共同交点为：k时刻目标的新生粒子分布边缘的k时刻共同交点，k时刻共同交点对应的目标速度为：k时刻目标的新生粒子分布边缘的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个新生粒子的权值

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

其中，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{\&CenterDot;}{x}}=2,Q_y=Q_{\stackrel{\&CenterDot;}{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个粒子的修正后权值

当k≥1时，对预测的k时刻的粒子集中第i'_k个粒子的权值进行修正的过程包括以下子步骤：

a)设置迭代参数j，j＝1,2,...当j＝1时，执行子步骤b)。

b)根据以下公式计算第j个发射站在处的检测概率

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

$\mathrm{SNR}(x_{k|k-1}^{\left({i^{\&prime;}}_k\right)},j)=\frac{K_j}{R_{T,j}^2{R}_R^2}$

$K_j=\frac{P_{T,j}{G}_{T,j}{G}_R{\mathrm{\&lambda;}}_{f,j}^2{\mathrm{\&sigma;}}_{\mathrm{rcs}}{F}_{T,j}^2{F}_R^2}{{\left(4\mathrm{\&pi;}\right)}^3{k}_0{T}_0(1/\mathrm{CPI})\left(\mathrm{NF}\right)}$

其中，表示第j个发射站对应的处的信噪比(本发明实施例中，为已知数值)，p_FA为设定的虚警概率，Q[·]表示marcum Q函数， $Q[a,b]=\underset{b}{\overset{\text{\&infin;}}{\text{\&Integral;}}}\mathrm{xexp}(-\frac{x^2+a^2}{2})I_0\left(\mathrm{ax}\right)\text{dx,}$ I₀表示零阶贝塞尔函数；R_T,j为预测的k时刻的粒子集中第i'_k个粒子对应的目标到第j个发射站的距离，R_R为预测的k时刻的粒子集中第i'_k个粒子对应的目标到外辐射源雷达(接收站)的距离；P_T,j为第j个发射站的信号发射功率，G_T,j为第j个发射站的发射天线增益，G_R为外辐射源雷达的接收天线增益，λ_f,j为第j个发射站的发射信号波长，σ_rcs为预测的k时刻的粒子集中第i'_k个粒子对应的目标的双基雷达(指第j个发射站与外辐射源雷达)反射面积，F_T,j为第j个发射站的发射信道中的信号传播衰减因子，F_R为外辐射源雷达的接收信道中的信号传播衰减因子。k₀为玻尔兹曼常数，T₀为环境温度，CPI为外辐射源雷达接收信号时的脉冲积累时间，NF指环境噪声系数(此处包含干扰)。

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

${\mathrm{\&lambda;}}_k=\frac{L_1}{l_1}\&CenterDot;\frac{L_2}{l_2}\&CenterDot;p_{\mathrm{FA}},$ c_k＝1/(L₁L₂)

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

d)在k时刻目标的观测集Z_k中，利用与第j个发射站对应的k时刻每个目标的观测量，组成k时刻第j个发射站的目标的观测集然后利用中的观测量计算似然函数

$g\left(z_k^{\left(j\right)}\right|x_{k|k-1}^{\left({i^{\&prime;}}_k\right)})=\frac{1}{2\mathrm{\&pi;}\sqrt{R_v}}\exp (-\frac{1}{2}(z_k^{\left(j\right)}-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_v是观测噪声ν_k的协方差矩阵，上标‐1表示矩阵的逆，上标T表示矩阵或向量的转置。

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

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

${\mathrm{\&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)\&CenterDot;g\left(z_k^{\left(j\right)}\right|x_{k|k-1}^{\left({i^{\&prime;}}_k\right)})$

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

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

然后，令j值自增1。

f)判断j与发射站的个数为l的大小关系，如果j<l，则令然后返回至子步骤b)，如果j＝l，则迭代终止，预测的k时刻的粒子集中第i'_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个粒子权值， ${\mathrm{\&omega;}}_k^{\left(i_k\right)}={\stackrel{\&OverBar;}{N}}_k/L_k.$

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

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

1)仿真条件

本发明的仿真实验中，设置三个发射站与一个接收站(外辐射源雷达)，本发明的仿真实验软件平台为MATLAB，操作系统为Win XP系统。每个发射站的信号发射功率为1kW,每个发射站的信号发射频率为780MHz，每个发射站的信号发射带宽为8MHz，每个发射站的发射天线增益为8dB，接收站的接收天线增益为13dB，信号在传播中无衰减。设定的虚警概率p_FA为10^-4，环境噪声系数(包含干扰)为30dB。目标在观测区域内匀速运动，且每一帧最多有一个新生目标出现，每个目标对每个发射站的RCS(双基雷达反射面积)都相同，均为10dBm²。目标存活概率为0.98，每个目标固定分配30个粒子，每个新生目标分配20个粒子。

2)仿真结果分析

参照图4a，为仿真实验中k时刻每个发射站对应的对应目标的双基地椭圆曲线与外辐射源雷达的矩形观测区域边缘的交点示意图。参照图4b，为图4a的局部放大图。图4b中，新生粒子分布在A点和B点的附近区域内。如图4a和图4b所示，A、B两点都为k时刻共同交点，若只利用目标的距离和信息和多普勒信息，无法进一步选择。本发明方法利用当前时刻观测量中角度信息判断新目标更可能出现的位置，从而剔除虚假公共交点。

参照图5a，为增加目标角度信息之前和之后的每时刻估计目标数和真实目标数的对比图，图5a中，横轴表示时间，纵轴表示目标数，无标示直线代表真实目标数，以三角形标示的直线代表增加目标角度信息之前的估计目标数，以菱形标示的增加目标角度信息之后(对应本发明)的估计目标数。

对增加目标角度信息之前和本发明(增加目标角度信息之后)得出的方位估计误差和速度估计误差进行对比，采用最优子模式分配(OptimalSub‐pattern Assignment,OSPA)作为多目标跟踪性能评估标准，其中距离敏感性参数选择2，关联敏感性参数选择5。参照图5b，为增加目标角度信息之前和之后得出的目标方位估计误差的对比图。图5b中，横轴表示时间，纵轴表示目标方位估计误差，单位为m，图5b中，以三角形标示的直线代表增加目标角度信息之前的目标方位估计误差，以菱形标示的直线代表增加目标角度信息之后(对应本发明)的目标方位估计误差。参照图5c，为增加目标角度信息之前和之后得出的目标速度估计误差的对比图。图5c中，横轴表示时间，纵轴表示目标速度估计误差，单位为m，图5c中，以三角形标示的直线代表增加目标角度信息之前的目标速度估计误差，以菱形标示的直线代表增加目标角度信息之后(对应本发明)的目标速度估计误差。从图5a至图5c可以看出，采用本发明之后，目标数目的估计精度虽然有所下降，但与增加目标角度信息之前相比，性能相差不大；但是增加目标角度信息后，目标的方位和速度估计精度得到了显著提高。此外，由于并非每个时刻都产生新生粒子，且在分布新生粒子时有效剔除了虚假点，增加目标角度信息后，降低了运算耗时。

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

Claims (7)

1.一种结合目标角度的外辐射源雷达多目标跟踪方法，其特征在于，包括以下步骤：

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

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

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

2.如权利要求1所述的一种结合目标角度的外辐射源雷达多目标跟踪方法，其特征在于，所述步骤2的具体子步骤为：

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

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

3.如权利要求1所述的一种结合目标角度的外辐射源雷达多目标跟踪方法，其特征在于，所述步骤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，并按照下式计算出所对应的权值

${\mathrm{\&omega;}}_{k|k-1}^{\left(i_{k-1}\right)}={\mathrm{\&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},\&CenterDot;\&CenterDot;\&CenterDot;,L_{k-1}$

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

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

根据k-1时刻每个目标相对于外辐射源雷达的角度、以及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时刻的状态。

4.如权利要求3所述的一种结合目标角度的外辐射源雷达多目标跟踪方法，其特征在于，在子步骤(3.3)中，首先，根据k-1时刻目标的观测集Z_k-1，获取k-1时刻每个目标相对于外辐射源雷达的角度，将k-1时每个目标相对于外辐射源雷达的角度组成k-1时刻目标观测量角度集合θ_k-1，k-1时刻目标观测量角度集合θ_k-1中元素的个数表示为<θ_k-1>，k-1时刻目标观测量角度集合θ_k-1中第α_k-1个元素表示为α_k-1取1至<θ_k-1>；根据k时刻目标的观测集Z_k，获取k时刻每个目标相对于外辐射源雷达的角度，将k时刻每个目标相对于外辐射源雷达的角度组成k时刻目标观测量角度集合θ_k，k时刻目标观测量角度集合θ_k中元素的个数表示为<θ_k>，k时刻目标观测量角度集合θ_k中第α_k个元素表示为α_k取1至<θ_k>；

判断k时刻是否产生新生粒子，如果k时刻不产生新生粒子，则k时刻新生粒子的个数J_k为零；如果k时刻产生新生粒子，则得出k时刻产生的每个新生粒子的分布区域中心点方位角，k时刻产生的每个新生粒子的分布区域中心点方位角指：k时刻产生的对应新生粒子的位置的纵坐标与横坐标的比值的反正切值；

判断k时刻是否产生新生粒子的过程为：将α_k-1从1至<θ_k-1>进行遍历，α_k-1每取一个数值时，在k时刻目标观测量角度集合θ_k中找出满足第一新生粒子存在条件的元素，所述第一新生粒子存在条件为：k时刻目标观测量角度集合θ_k中的元素与的差值的绝对值不小于γ_θ，γ_θ为设定的角度门限；在将α_k-1从1至<θ_k-1>进行遍历的过程中，若k时刻目标观测量角度集合θ_k中不存在满足第一新生粒子存在条件的元素，则k时刻不产生新生粒子；反之，k时刻产生新生粒子；

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

得出k时刻目标的新生粒子分布边缘；得出k时刻目标的新生粒子分布边缘的过程为：当k时刻目标的候选新生粒子分布边缘的个数等于1时，k时刻目标的新生粒子分布边缘为k时刻目标的候选新生粒子分布边缘；当k时刻目标的候选新生粒子分布边缘的个数大于1时，在k时刻目标的所有候选新生粒子分布边缘中找出k时刻目标的新生粒子分布边缘；

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

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

当k时刻目标的候选新生粒子分布边缘的个数大于1时，找出k时刻目标的新生粒子分布边缘的过程为：判断k时刻是否产生新生粒子，如果k时刻不产生新生粒子，则没有k时刻目标的新生粒子分布边缘；如果k时刻产生新生粒子，针对k时刻的每个新生粒子，判断每个候选新生粒子分布边缘的k时刻共同交点的方位角是否满足第二新生粒子存在条件，所述第二新生粒子存在条件为：每个候选新生粒子分布边缘的k时刻共同交点的方位角与k时刻目标观测量角度集合θ_k中满足第一新生粒子存在条件的元素的差值的绝对值不大于设定的方位角门限值；如果每个候选新生粒子分布边缘的k时刻共同交点的方位角不满足第二新生粒子存在条件，则令k时刻新生粒子的个数J_k取0；如果至少一个候选新生粒子分布边缘满足第二新生粒子存在条件，则在满足第二新生粒子存在条件的候选新生粒子分布边缘中，选择其k时刻共同交点的方位角与k时刻的对应新生粒子的分布区域中心点方位角最接近的候选新生粒子分布边缘，则k时刻目标的新生粒子分布边缘为：在满足第二新生粒子存在条件的候选新生粒子分布边缘中选择的候选新生粒子分布边缘，此时，k时刻共同交点为：k时刻目标的新生粒子分布边缘的k时刻共同交点，k时刻共同交点对应的目标速度为：k时刻目标的新生粒子分布边缘的k时刻共同交点对应的目标速度；

在k时刻共同交点处采样J_k个点，得出k时刻J_k个新生粒子；k时刻每个新生粒子的粒子状态包括：k时刻对应新生粒子的方位状态、k时刻对应新生粒子的速度状态；

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

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

其中，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{\&CenterDot;}{x}}=2,Q_y=Q_{\stackrel{\&CenterDot;}{y}}=1.$

5.如权利要求3所述的一种结合目标角度的外辐射源雷达多目标跟踪方法，其特征在于，在子步骤(3.4)中，当k≥1时，当k≥1时，对预测的k时刻的粒子集中第i'_k个粒子的权值进行修正的过程包括以下子步骤：

a)设置迭代参数j，j＝1,2,...当j＝1时，执行子步骤b)；

b)根据以下公式计算第j个发射站在处的检测概率

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

其中，表示已知的第j个发射站对应的处的信噪比，p_FA为设定的虚警概率， $Q[a,b]=\underset{b}{\overset{\text{\&infin;}}{\text{\&Integral;}}}\mathrm{xexp}(-\frac{x^2+a^2}{2})I_0\left(\mathrm{ax}\right)\text{dx,}$ I₀表示零阶贝塞尔函数；

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

${\mathrm{\&lambda;}}_k=\frac{L_1}{l_1}\&CenterDot;\frac{L_2}{l_2}\&CenterDot;p_{\mathrm{FA}},$ c_k＝1/(L₁L₂)

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

d)在k时刻目标的观测集Z_k中，利用与第j个发射站对应的k时刻每个目标的观测量，组成k时刻第j个发射站的目标的观测集然后利用中的观测量计算似然函数

$g\left(z_k^{\left(j\right)}\right|x_{k|k-1}^{\left({i^{\&prime;}}_k\right)})=\frac{1}{2\mathrm{\&pi;}\sqrt{R_v}}\exp (-\frac{1}{2}(z_k^{\left(j\right)}-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_v是观测噪声ν_k的协方差矩阵，上标‐1表示矩阵的逆，上标T表示矩阵或向量的转置；

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

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

${\mathrm{\&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)\&CenterDot;g\left(z_k^{\left(j\right)}\right|x_{k|k-1}^{\left({i^{\&prime;}}_k\right)})$

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

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

然后，令j值自增1；

f)判断j与发射站的个数为l的大小关系，如果j<l，则令然后返回至子步骤b)，如果j＝l，则迭代终止，预测的k时刻的粒子集中第i'_k个粒子的修正后权值为

6.如权利要求3所述的一种结合目标角度的外辐射源雷达多目标跟踪方法，其特征在于，在子步骤(3.5)中，首先根据预测的k时刻的粒子集中第i'_k个粒子的修正后权值计算k时刻的估计目标数

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

7.如权利要求3所述的一种结合目标角度的外辐射源雷达多目标跟踪方法，其特征在于，在步骤3之后，针对每个目标0时刻至k时刻的状态，进行航迹关联，剔除虚假航迹，得出每个目标的航迹。