CN105929378A - 基于外辐射源联合时延与多普勒频率的直接跟踪方法 - Google Patents

Abstract

本发明涉及一种基于外辐射源联合时延与多普勒频率的直接跟踪方法，首先建立目标的状态转移模型以及包含外辐射源的直达波和回波的接收信号模型；通过计算接收信号的傅立叶系数将时域数据转化成频域数据，对频域接收数据构建高维最大似然估计，利用信息矩阵的最大特征值作为粒子滤波中粒子的后验概率加权，通过重采样获得权值较大的粒子，对所得粒子取平均得到对应时刻目标位置的估计。本发明相比于传统跟踪方法，综合考虑外辐射源的直达波信息与目标反射回波信息，构建包含时延与多普勒信息的多维信号模型，利用底层接收数据构建粒子后验概率加权，直接对目标位置状态进行估计，跟踪精度更高，且方法实现简单、高效，性能稳健、可靠。

Description

基于外辐射源联合时延与多普勒频率的直接跟踪方法

技术领域

本发明涉及外辐射源定位跟踪技术领域，特别涉及一种基于外辐射源联合时延与多普勒频率的直接跟踪方法。

背景技术

外辐射源无源定位与跟踪，指的是利用非合作的第三方辐射源(例如AM/FM信号、普通/数字电视信号等)作为目标的照射源，通过对直达波及反射回波进行处理，以得到时延、频率等相关参数信息，进而实现对目标的定位与跟踪。利用第三方的非合作照射源对目标进行探测，不仅能够实现对隐形目标和静默目标的探测、定位与跟踪，同样也可以用于对目标的成像与识别。是诸多军用和民用应用领域的重要组成部分，如电子对抗、雷达信号处理、空中交通管制、无线电监测、移动通信、遥测与导航等。

目前外辐射源无源定位与跟踪采用的体制是传统的无源定位与跟踪体制，其处理流程为首先进行参数估计，如到达角度、到达时间、到达时间差、多普勒频差、接收信号强度或多种参数联合估计，再通过对获取的参数进行对非线性方程的求解获得目标的位置估计。若目标运动，除了需要考虑目标的坐标信息外，还需要考虑目标的速度信息及其带来的状态变化；目标运动中的参数信息及其变化率均是状态变量的的非线性函数，必须采用非线性滤波算法，如卡尔曼滤波、扩展卡尔曼滤波、无迹卡尔曼滤波、粒子滤波算法等，否则定位系统的性能将急剧下降。寻找一种稳健且收敛速度快的跟踪滤波算法是辐射源目标跟踪中需要解决的首要问题，传统的目标跟踪处理方法为首先进行参数估计形成目标的点迹，再利用滤波算法对目标的航迹进行跟踪处理。由于参数估计与目标滤波跟踪相分离，无法保证测量的参数结果与真实目标航迹信息相匹配，同时参数估计的误差可能在后续的处理过程中被进一步放大且很难被消除从而导致整个数据处理过程中不可避免地存在信息损失，所以无法获得最优的估计性能。从信息论的角度来看，从原始接收数据到最终的处理结果，增加中间的处理环节将引入更多的不确定性，导致部分信息的损失，所以传统的目标跟踪处理方法很难取得最优的结果。为了克服传统目标跟踪方法的缺点，以色列学者A.J.Weiss和Alon Y.Sidi学者提出了单步跟踪(直接目标跟踪)方法，该方法直接从接收信号数据中获取目标的航迹信息，无需参数估计步骤，避免了误差的引入，提高了目标跟踪的精度，该方法仅考虑单根天线时的情况，未对阵列天线的情况进行研究。

发明内容

针对现有技术中的不足，本发明提供一种基于外辐射源联合时延与多普勒频率的直接跟踪方法，利用接收底层数据构造粒子后验概率权重，避免由参数估计误差引起权重与真实位置失配的问题，减少了目标跟踪信息的损失，跟踪精度明显提升，且对信噪比具有较强的鲁棒性。

按照本发明所提供的设计方案，一种基于外辐射源联合时延与多普勒频率的直接跟踪方法，具体包含如下步骤：

步骤1.外辐射源场景下利用状态转移矩阵F构建目标运动状态方程为：

x_k＝Fx_k-1+v_k-1，

其中，为状态向量，x_k,y_k与表分别示在第k次观测间隙内目标辐射源的位置坐标和相应的速度分量；T_r表示采样周期，状态转移矩阵F为：

$F=\left[\begin{array}{cccc}1& T_r& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& T_r\\ 0& 0& 0& 1\end{array}\right]$ ，

并定义初始状态与跟踪状态时的条件概率密度函数；结合粒子滤波方法构建后验概率加权迭代方程为：

$w_k^{\left(i\right)}={\tilde{w}}_k^{\left(i\right)}/\underset{i=1}{\overset{M}{\Σ}}{\tilde{w}}_k^{\left(i\right)},{\tilde{w}}_k^{\left(i\right)}=w_{k-1}^{\left(i\right)}p\left(z_k\right|x_k^{\left(i\right)})$ ，

其中，z_k表示k时刻的观测量；

步骤2.对L个观测站的双通道接收系统进行时间同步，并根据Nyquist采样定理采集外辐射源的直达波信号以及经目标反射的回波信号，获得多站接收的信号时域模型；

步骤3.针对信号时域模型，对各站双通道接收的时域数据分别计算其傅立叶系数，得到阵列信号频域模型；

步骤4.针对阵列信号频域模型，每个观测站将所获得的阵列信号频域数据传输至中心站，中心站将每个站传输的阵列信号频域数据按照观测站的顺序进行堆栈排列，构造高维阵列信号频域模型；

步骤5.针对高维阵列信号频域模型，在中心站对高维阵列信号频域数据构造高斯最大似然函数，并构造包含回波时延、多普勒信息以及直达波时延信息的数据信息矩阵；

步骤6.数据信息矩阵的最大特征值作为粒子滤波的后验概率权重；根据粒子滤波理论，设定粒子数量及初始权重，通过迭代更新粒子的后验概率权重并进行重采样，所得最终粒子的均值作为该时刻目标跟踪的结果。

上述的，步骤2中，第l个k时刻观测站的所接收到的信号时域模型为：

$r_{l,k}\left(t\right)=\left[\begin{array}{c}{b}_{d_l}s(t-{\mathrm{\&tau;}}_{d_l})\\ b_{l,k}s(t-{\mathrm{\&tau;}}_{l,k})e^{j2{\mathrm{\&pi;f}}_{l,k}t}\end{array}\right]+\left[\begin{array}{c}{w}_{d_l}\left(t\right)\\ w_{l,k}\left(t\right)\end{array}\right],0\&le;t<T;$

其中，表示外辐射源相对于第l个观测站的直达路径时延；τ_l,k＝(||p_e-p_0,k||+||p_l-p_0,k||)/c表示外辐射源照射k时刻的目标并反射至观测站产生的时延，c表示信号传播速度，||·||表示2范数；p_e为外辐射源位置，发射信号带宽W，p_0,k为目标在k时刻时的位置，速度为v_k＝[v_x,k,v_y,k]^T；与w_l,k(t)分别表示均值为0，方差为与直达波通道与k时刻回波通道的加性平稳复高斯白噪声；b_l,k与分别表示k时刻回波信道与直达波信道的衰减系数；f_l,k表示目标与观测站之间的多普勒频率，其包含两部分，一部分为外辐射源照射至目标时信号的多普勒频率，另一部分是k时刻反射回波到达观测站时产生的多普勒频率，故f_l,k表示为：

$f_{l,k}\stackrel{\mathrm{\&Delta;}}{=}\frac{f_c}{c}(\frac{v_k^T(p_{0,k}-p_l)}{\left|\right|p_{0,k}-p_l\left|\right|}+\frac{v_k^T(p_{0,k}-p_e)}{\left|\right|p_{0,k}-p_e\left|\right|})$ ，

其中，p_l＝[x_l,y_l]^T(l＝1,2,…,L)为第l个观测站的位置坐标。

上述的，步骤3中第l个观测站k时刻的所接收到的信号频域模型为：

${\tilde{r}}_{l,k}\left(f_n\right)=\left[\begin{array}{c}{b}_{d_l}\tilde{s}\left(f_n\right)e^{-j2{\mathrm{\&pi;f}}_n{\mathrm{\&tau;}}_{d_l}}\\ b_{l,k}\tilde{s}(f_n-f_{l,k})e^{-j2{\mathrm{\&pi;f}}_n{\mathrm{\&tau;}}_{l,k}}\end{array}\right]+\left[\begin{array}{c}{\tilde{w}}_{d_l}\left(f_n\right)\\ {\tilde{w}}_{l,k}\left(f_n\right)\end{array}\right]$ ，

其中，与分别表示接收信号、发射信号以及噪声项的傅立叶系数。

上述的，步骤4中，中心站所获得的高维信号频域模型为：

${\tilde{r}}_{l,k}=\left[\begin{array}{c}{b}_{d_l}{\tilde{B}}_l\\ b_{l,k}{\tilde{A}}_{l,k}{\tilde{F}}_{l,k}\end{array}\right]\tilde{s}+\left[\begin{array}{c}{\tilde{w}}_{d_l}\\ {\tilde{w}}_{l,k}\end{array}\right]$ ，

其中，

$\left\{\begin{array}{c}{\tilde{r}}_{l,k}={\&lsqb;{\tilde{r}}_{l,k}\left(f_{-2N}\right),{\tilde{r}}_{l,k}\left(f_{-2N+1}\right),\mathrm{...},{\tilde{r}}_{l,k}\left(f_{2N}\right)\&rsqb;}^T\\ \tilde{s}={\&lsqb;\tilde{s}\left(f_{-N}\right),\tilde{s}\left(f_{-N+1}\right),\mathrm{...},\tilde{s}\left(f_N\right)\&rsqb;}^T\\ {\tilde{B}}_l=diag\{e^{-j2{\mathrm{\&pi;f}}_{-N}{\mathrm{\&tau;}}_{d_l}},e^{-j2{\mathrm{\&pi;f}}_{-N+1}{\mathrm{\&tau;}}_{d_l}},\mathrm{...},e^{-j2{\mathrm{\&pi;f}}_N{\mathrm{\&tau;}}_{d_l}}\}\\ {\tilde{A}}_{l,k}=diag\{e^{-j2{\mathrm{\&pi;f}}_{-N}{\mathrm{\&tau;}}_{l,k}},e^{-j2{\mathrm{\&pi;f}}_{-N+1}{\mathrm{\&tau;}}_{l,k}},\mathrm{...},e^{-j2{\mathrm{\&pi;f}}_N{\mathrm{\&tau;}}_{l,k}}\}\\ {\tilde{w}}_{l,k}={\&lsqb;{\tilde{w}}_{l,k}\left(f_{-N}\right),{\tilde{w}}_{l,k}\left(f_{-N+1}\right),\mathrm{...},{\tilde{w}}_{l,k}\left(f_N\right)\&rsqb;}^T\\ {\tilde{w}}_{d_l}={\&lsqb;{\tilde{w}}_{d_l}\left(f_{-N}\right),{\tilde{w}}_{d_l}\left(f_{-N+1}\right),\mathrm{...},{\tilde{w}}_{d_l}\left(f_N\right)\&rsqb;}^T\end{array}\right.$ ，

式中，表示循环移位矩阵，其形式表示为：

，

即将单位阵循环向下移位floor(Tf_l,k)+1行；即将矩阵按行循环向下移位floor(Tf_l,k)，floor(·)表示向下取整。

上述的，步骤5中，构造的高斯最大似然函数为：

$p\left({\tilde{r}}_k\right|x_k^{\left(i\right)})\&Proportional;\exp \left(-\frac{1}{\det \left({\tilde{C}}_w\right)}\underset{l=1}{\overset{L}{\&Sigma;}}\left|\right|{\tilde{r}}_{l,k}-b_l\left(\begin{array}{c}{\tilde{A}}_{l,k}^{\left(i\right)}{\tilde{F}}_{l,k}^{\left(i\right)}\\ {\tilde{B}}_l\end{array}\right)\tilde{s}|{|}^2\right)$ ，

其中，det(·)表示求行列式； 为了估计目标的位置坐标，应使上式右端最小化，此时第i个粒子似然表示如下：

$p\left(r_k\right|x_k^{\left(i\right)})\&Proportional;\exp (-\frac{1}{\det \left({\tilde{C}}_w\right)}\underset{p_{0,k}}{\arg \min }{C}_1(x_k^{\left(i\right)},p_{0,k}\left)\right)$ ，

其中，代价函数为：

$C_1(x_k^{\left(i\right)},p_{0,k})=\underset{l=1}{\overset{L}{\&Sigma;}}\left|\right|{\tilde{r}}_{l,k}-b_{l,k}\left(\begin{array}{c}{\tilde{A}}_{l,k}^{\left(i\right)}{\tilde{F}}_{l,k}^{\left(i\right)}\\ {\tilde{B}}_l\end{array}\right)\tilde{s}|{|}^2$ ，

利用b_l,k使得上式最小化，可得：

$b_{l,k}=\&lsqb;\left(\begin{array}{c}{\tilde{A}}_{l,k}^{\left(i\right)}{\tilde{F}}_{l,k}^{\left(i\right)}\\ {\tilde{B}}_{l,k}^{\left(i\right)}\end{array}\right){\tilde{s}}^H\tilde{s}\left(\begin{array}{c}{\tilde{A}}_{l,k}^{\left(i\right)}{\tilde{F}}_{l,k}^{\left(i\right)}\\ {\tilde{B}}_{l,k}^{\left(i\right)}\end{array}\right)\&rsqb;{\left({\tilde{Q}}_{l,k}^{\left(i\right)}\tilde{s}\right)}^H{\tilde{r}}_{l,k}$ ，

其中，不失一般性，假设结合矩阵的特殊结构，可得：

$b_{l,k}=2{\left({\tilde{Q}}_{l,k}^{\left(i\right)}\tilde{s}\right)}^H{\tilde{r}}_{l,k}$ ，

经过数学推导简化，得到如下形式：

$C_2(x_k^{\left(i\right)},p_{0,k})=\underset{l=1}{\overset{L}{\&Sigma;}}\left|\right|{\tilde{r}}_{l,k}-2{\left({\tilde{Q}}_{l,k}^{\left(i\right)}\tilde{s}\right)}^H{\tilde{r}}_{l,k}{\tilde{Q}}_{l,k}^{\left(i\right)}\tilde{s}|{|}^2=\underset{l=1}{\overset{L}{\&Sigma;}}|\left|{\tilde{r}}_{l,k}\right|{|}^2-\sqrt{2}|{\left({\tilde{Q}}_{l,k}^{\left(i\right)}\tilde{s}\right)}^H{\tilde{r}}_{l,k}{|}^2$ ，

上式的最小化等价于式子右端的最大化，即构造代价函数为：

$C_3(x_k^{\left(i\right)},p_{0,k})=\underset{l=1}{\overset{L}{\&Sigma;}}|{\left({\tilde{Q}}_{l,k}^{\left(i\right)}\tilde{s}\right)}^H{\tilde{r}}_{l,k}{|}^2={\tilde{s}}^H{Q}_{c,k}^{\left(i\right)}\tilde{s}$ ，

其中，

$Q_{c,k}^{\left(i\right)}=V_k^{\left(i\right)}{V}_k^{\left(i\right)H}$ ，

数据信息矩阵表示为：

$V_k^{\left(i\right)}=\&lsqb;({\tilde{F}}_{1,k}^{\left(i\right)H}{\tilde{A}}_{1,k}^{\left(i\right)H},{\tilde{B}}_1^H){\tilde{r}}_{1,k}({\tilde{F}}_{2,k}^{\left(i\right)H}{\tilde{A}}_{2,k}^{\left(i\right)H},{\tilde{B}}_2^H){\tilde{r}}_{2,k},\mathrm{...},({\tilde{F}}_{L,k}^{\left(i\right)H}{\tilde{A}}_{L,k}^{\left(i\right)H},{\tilde{B}}_L^H){\tilde{r}}_{L,k}\&rsqb;.$

上述的，步骤6中具体包含如下：取粒子对应数据信息矩阵对应的最大特征值为：

$C_3(x_k^{\left(i\right)},p_{0,k})={\mathrm{\&lambda;}}_{\max }\left(Q_{c,k}^{\left(i\right)}\right)$ ，

其中，λ_max(·)表示求最大特征值；的维度为(2N+1)×(2N+1)，考虑给定矩阵X，X^HX与XX^H的非零特征值是一致的，故将变换为

$C_3(x_k^{\left(i\right)},p_{0,k})={\mathrm{\&lambda;}}_{\max }\left({\stackrel{\&OverBar;}{Q}}_{c,k}^{\left(i\right)}\right)$ ，

式中， 的维度为L×L维且仅与观测站数量有关，利用该代价函数计算粒子滤波加权，即：

${\tilde{w}}_k\&Proportional;\exp \left(\frac{1}{\det \left({\tilde{C}}_w\right)}{\mathrm{\&lambda;}}_{max}\right({\stackrel{\&OverBar;}{Q}}_{c,k}^{\left(i\right)}\left)\right)$ ，

通过迭代更新粒子后验概率权重并重采样，所得最终粒子的均值即为该时刻目标跟踪的结果。

本发明的有益效果：

1、本发明结合外辐射源无源定位，为计算粒子后验概率权重，首先通过计算回波信号与直达波信号的傅立叶系数将时域数据转化成频域数据，然后对频域数据构建高维数据模型，结合高斯最大似然估计器得到估计代价函数；最后利用信息矩阵的最大特征值作为粒子后验概率权重，实现基于直接数据域的目标跟踪；与传统跟踪算法相比，本发明直接利用接收底层数据构造粒子后验概率权重，避免了由参数估计误差引起权重与真实位置信息失配的问题，减少了目标跟踪信息的损失，跟踪精度明显提升，且对信噪比具有较强的鲁棒性。

2、本发明相比于传统跟踪方法中先对目标的参数进行估计再进行滤波跟踪，本发明综合考虑了外辐射源的直达波信息与目标反射回波信息，构建了包含时延与多普勒信息的多维信号模型，通过利用底层接收数据构建粒子的后验概率加权，并直接对目标位置状态进行估计，获取了更高的跟踪精度，且方法实现简单、高效，性能稳健、可靠，实现高精度定位跟踪。

附图说明：

图1为本发明的流程示意图；

图2为本发明的实例场景示意图；

图3为相对于传统基于时差与频差跟踪方法的跟踪性能对比图；

图4为相对于传统基于时差与频差跟踪方法的不同时刻下跟踪误差对比图；

图5为时刻k＝15时算法定位误差对比图。

具体实施方式：

下面结合附图和技术方案对本发明作进一步详细的说明，并通过优选的实施例详细说明本发明的实施方式，但本发明的实施方式并不限于此。

实施例一，参见图1所示，一种基于外辐射源联合时延与多普勒频率的直接跟踪方法，具体包含如下步骤：

步骤1.外辐射源场景下利用状态转移矩阵F构建目标运动状态方程为：

x_k＝Fx_k-1+v_k-1，

其中，为状态向量，x_k,y_k与表分别示在第k次观测间隙内目标辐射源的位置坐标和相应的速度分量；T_r表示采样周期，状态转移矩阵F为：

$F=\left[\begin{array}{cccc}1& T_r& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& T_r\\ 0& 0& 0& 1\end{array}\right]$ ，

并定义初始状态与跟踪状态时的条件概率密度函数；结合粒子滤波方法构建后验概率加权迭代方程为：

$w_k^{\left(i\right)}={\tilde{w}}_k^{\left(i\right)}/\underset{i=1}{\overset{M}{\&Sigma;}}{\tilde{w}}_k^{\left(i\right)},{\tilde{w}}_k^{\left(i\right)}=w_{k-1}^{\left(i\right)}p\left(z_k\right|x_k^{\left(i\right)})$ ，

其中，z_k表示k时刻的观测量；

步骤2.对L个观测站的双通道接收系统进行时间同步，并根据Nyquist采样定理采集外辐射源的直达波信号以及经目标反射的回波信号，获得多站接收的信号时域模型；

步骤3.针对信号时域模型，对各站双通道接收的时域数据分别计算其傅立叶系数，得到阵列信号频域模型；

步骤4.针对阵列信号频域模型，每个观测站将所获得的阵列信号频域数据传输至中心站，中心站将每个站传输的阵列信号频域数据按照观测站的顺序进行堆栈排列，构造高维阵列信号频域模型；

步骤5.针对高维阵列信号频域模型，在中心站对高维阵列信号频域数据构造高斯最大似然函数，并构造包含回波时延、多普勒信息以及直达波时延信息的数据信息矩阵；

步骤6.数据信息矩阵的最大特征值作为粒子滤波的后验概率权重；根据粒子滤波理论，设定粒子数量及初始权重，通过迭代更新粒子的后验概率权重并进行重采样，所得最终粒子的均值作为该时刻目标跟踪的结果。

实施例二，参见图1～4，一种基于外辐射源联合时延与多普勒频率的直接跟踪方法，具体包含如下步骤：

步骤1.外辐射源场景下利用状态转移矩阵F构建目标运动状态方程为：

x_k＝Fx_k-1+v_k-1，

其中，为状态向量，x_k,y_k与表分别示在第k次观测间隙内目标辐射源的位置坐标和相应的速度分量；T_r表示采样周期，状态转移矩阵F为：

$F=\left[\begin{array}{cccc}1& T_r& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& T_r\\ 0& 0& 0& 1\end{array}\right]$ ，

并定义初始状态与跟踪状态时的条件概率密度函数，若假设状态噪声v_k-1为零均值白高斯过程，其协方差矩阵为：

${\mathrm{\&Lambda;}}_v=E\&lsqb;v_k{v}_k^T\&rsqb;=q\left[\begin{array}{cccc}{T}_r^3/3& T_r^2/2& 0& 0\\ T_r^2/2& T_r& 0& 0\\ 0& 0& T_r^3/3& T_r^2/2\\ 0& 0& T_r^2/2& T_r\end{array}\right]$ ，

其中，q为系统定义参数；目标的条件概率密度函数p(x_k+1|x_k)同样为高斯过程，即p(x_k+1|x_k)～N(Fx_k,Λ_v)；对于目标初始状态x₁，假设为已知的先验高斯概率密度，即结合粒子滤波方法构建后验概率加权迭代方程为：

$w_k^{\left(i\right)}={\tilde{w}}_k^{\left(i\right)}/\underset{i=1}{\overset{M}{\&Sigma;}}{\tilde{w}}_k^{\left(i\right)},{\tilde{w}}_k^{\left(i\right)}=w_{k-1}^{\left(i\right)}p\left(z_k\right|x_k^{\left(i\right)})$ ，

其中，z_k表示k时刻的观测量，考虑粒子滤波中权重更新迭代：粒子滤波的基本思想是在给定一组观测量Z_k＝[z₁,z₂,…,z_k]的基础上估计器后验概率密度p(x_k|Z_k)，一般通过若干离散样本(粒子)来估计：

$p\left(x_k\right|Z_k)\&ap;\hat{p}\left(x_k\right|Z_k)=\underset{i=1}{\overset{M}{\&Sigma;}}{w}_k^{\left(i\right)}\mathrm{\&delta;}(x_k-x_k^{\left(i\right)})$ ，

其中，表示k时刻第i个粒子的权值，且当粒子数目足够多时，上式的估计越接近真实的p(x_k|Z_k)；一般由下式给出

${\tilde{w}}_k^{\left(i\right)}=w_{k-1}^{\left(i\right)}\frac{p\left(z_k\right|x_k^{\left(i\right)})p\left(x_k^{\left(i\right)}\right|x_{k-1}^{\left(i\right)})}{q\left(x_k^{\left(i\right)}\right|x_{k-1}^{\left(i\right)},z_k)}$

$w_k^{\left(i\right)}={\tilde{w}}_k^{\left(i\right)}/\underset{i=1}{\overset{M}{\&Sigma;}}{\tilde{w}}_k^{\left(i\right)}$ ，

式中，表示重要性函数，一般选择为则粒子权值可写为：为了避免粒子滤波中粒子退化的问题，最后需要对粒子进行重采样；

步骤2.对L个观测站的双通道接收系统进行时间同步，并根据Nyquist采样定理采集外辐射源的直达波信号以及经目标反射的回波信号，获得多站接收的信号时域模型；第l个k时刻观测站的所接收到的信号时域模型为：

$r_{l,k}\left(t\right)=\left[\begin{array}{c}{b}_{d_l}s(t-{\mathrm{\&tau;}}_{d_l})\\ b_{l,k}s(t-{\mathrm{\&tau;}}_{l,k})e^{j2{\mathrm{\&pi;f}}_{l,k}t}\end{array}\right]+\left[\begin{array}{c}{w}_{d_l}\left(t\right)\\ w_{l,k}\left(t\right)\end{array}\right],0\&le;t<T;$

其中，表示外辐射源相对于第l个观测站的直达路径时延；τ_l,k＝(||p_e-p_0,k||+||p_l-p_0,k||)/c表示外辐射源照射k时刻的目标并反射至观测站产生的时延，c表示信号传播速度，||·||表示2范数；p_e为外辐射源位置，发射信号带宽W，p_0,k为目标在k时刻时的位置，速度为v_k＝[v_x,k,v_y,k]^T；与w_l,k(t)分别表示均值为0，方差为与直达波通道与k时刻回波通道的加性平稳复高斯白噪声；b_l,k与分别表示k时刻回波信道与直达波信道的衰减系数；f_l,k表示目标与观测站之间的多普勒频率，其包含两部分，一部分为外辐射源照射至目标时信号的多普勒频率，另一部分是k时刻反射回波到达观测站时产生的多普勒频率，故f_l,k表示为：

$f_{l,k}=\frac{f_c}{c}(\frac{v_k^T(p_{0,k}-p_l)}{\left|\right|p_{0,k}-p_l\left|\right|}+\frac{v_k^T(p_{0,k}-p_e)}{\left|\right|p_{0,k}-p_e\left|\right|})$ ，

其中，p_l＝[x_l,y_l]^T(l＝1,2,…,L)为第l个观测站的位置坐标；

步骤3.针对信号时域模型，对各站双通道接收的时域数据分别计算其傅立叶系数，得到阵列信号频域模型；第l个观测站k时刻的所接收到的信号频域模型为：

${\tilde{r}}_{l,k}\left(f_n\right)=\left[\begin{array}{c}{b}_{d_l}\tilde{s}\left(f_n\right)e^{-j2{\mathrm{\&pi;f}}_n{\mathrm{\&tau;}}_{d_l}}\\ b_{l,k}\tilde{s}(f_n-f_{l,k})e^{-j2{\mathrm{\&pi;f}}_n{\mathrm{\&tau;}}_{l,k}}\end{array}\right]+\left[\begin{array}{c}{\tilde{w}}_{d_l}\left(f_n\right)\\ {\tilde{w}}_{l,k}\left(f_n\right)\end{array}\right]$ ，

其中，与分别表示接收信号、发射信号以及噪声项的傅立叶系数；

步骤4.针对阵列信号频域模型，每个观测站将所获得的阵列信号频域数据传输至中心站，中心站将每个站传输的阵列信号频域数据按照观测站的顺序进行堆栈排列，构造高维阵列信号频域模型；中心站所获得的高维信号频域模型为：

${\tilde{r}}_{l,k}=\left[\begin{array}{c}{b}_{d_l}{\tilde{B}}_l\\ b_{l,k}{\tilde{A}}_{l,k}{\tilde{F}}_{l,k}\end{array}\right]\tilde{s}+\left[\begin{array}{c}{\tilde{w}}_{d_l}\\ {\tilde{w}}_{l,k}\end{array}\right]$ ，

其中，

$\left\{\begin{array}{c}{\tilde{r}}_{l,k}={\&lsqb;{\tilde{r}}_{l,k}\left(f_{-2N}\right),{\tilde{r}}_{l,k}\left(f_{-2N+1}\right),\mathrm{...},{\tilde{r}}_{l,k}\left(f_{2N}\right)\&rsqb;}^T\\ \tilde{s}={\&lsqb;\tilde{s}\left(f_{-N}\right),\tilde{s}\left(f_{-N+1}\right),\mathrm{...},\tilde{s}\left(f_N\right)\&rsqb;}^T\\ {\tilde{B}}_l=diag\{e^{-j2{\mathrm{\&pi;f}}_{-N}{\mathrm{\&tau;}}_{d_l}},e^{-j2{\mathrm{\&pi;f}}_{-N+1}{\mathrm{\&tau;}}_{d_l}},\mathrm{...},e^{-j2{\mathrm{\&pi;f}}_N{\mathrm{\&tau;}}_{d_l}}\}\\ {\tilde{A}}_{l,k}=diag\{e^{-j2{\mathrm{\&pi;f}}_{-N}{\mathrm{\&tau;}}_{l,k}},e^{-j2{\mathrm{\&pi;f}}_{-N+1}{\mathrm{\&tau;}}_{l,k}},\mathrm{...},e^{-j2{\mathrm{\&pi;f}}_N{\mathrm{\&tau;}}_{l,k}}\}\\ {\tilde{w}}_{l,k}={\&lsqb;{\tilde{w}}_{l,k}\left(f_{-N}\right),{\tilde{w}}_{l,k}\left(f_{-N+1}\right),\mathrm{...},{\tilde{w}}_{l,k}\left(f_N\right)\&rsqb;}^T\\ {\tilde{w}}_{d_l}={\&lsqb;{\tilde{w}}_{d_l}\left(f_{-N}\right),{\tilde{w}}_{d_l}\left(f_{-N+1}\right),\mathrm{...},{\tilde{w}}_{d_l}\left(f_N\right)\&rsqb;}^T\end{array}\right.$ ，

式中，表示循环移位矩阵，其形式表示为：

，

即将单位阵循环向下移位floor(Tf_l,k)+1行；即将矩阵按行循环向下移位floor(Tf_l,k)，floor(·)表示向下取整

步骤5.针对高维阵列信号频域模型，在中心站对高维阵列信号频域数据构造高斯最大似然函数，并构造包含回波时延、多普勒信息以及直达波时延信息的数据信息矩阵；构造的高斯最大似然函数为：

$p\left(\tilde{r}\right|x_k^{\left(i\right)})\&Proportional;\exp \left(-\frac{1}{\det \left({\tilde{C}}_w\right)}\underset{l=1}{\overset{L}{\&Sigma;}}\left|\right|{\tilde{r}}_{l,k}-b_l\left(\begin{array}{c}{\tilde{A}}_{l,k_l}^{\left(i\right)}{\tilde{F}}_{l,k}^{\left(i\right)}\\ {\tilde{B}}_l\end{array}\right)\tilde{s}|{|}^2\right)$ ，

其中，det(·)表示求行列式； 为了估计目标的位置坐标，应使上式右端最小化，此时第i个粒子似然表示如下：

$p\left(r_k\right|x_k^{\left(i\right)})\&Proportional;\exp (-\frac{1}{\det \left({\tilde{C}}_w\right)}\underset{p_{0,k}}{\arg \min }{C}_1(x_k^{\left(i\right)},p_{0,k}\left)\right)$ ,

其中，代价函数为：

$C_1(x_k^{\left(i\right)},p_{0,k})=\underset{l=1}{\overset{L}{\&Sigma;}}\left|\right|{\tilde{r}}_{l,k}-b_{l,k}\left(\begin{array}{c}{\tilde{A}}_{l,k}^{\left(i\right)}{\tilde{F}}_{l,k}^{\left(i\right)}\\ {\tilde{B}}_l\end{array}\right)\tilde{s}|{|}^2$

利用b_l,k使得上式最小化，可得：

$b_{l,k}=\&lsqb;{\left(\begin{array}{c}{\tilde{A}}_{l,k}^{\left(i\right)}{\tilde{F}}_{l,k}^{\left(i\right)}\\ {\tilde{B}}_{l,k}^{\left(i\right)}\end{array}\right)}^H{\tilde{s}}^H\tilde{s}\left(\begin{array}{c}{\tilde{A}}_{l,k}^{\left(i\right)}{\tilde{F}}_{l,k}^{\left(i\right)}\\ {\tilde{B}}_{l,k}^{\left(i\right)}\end{array}\right)\&rsqb;{\left({\tilde{Q}}_{l,k}^{\left(i\right)}\tilde{s}\right)}^H{\tilde{r}}_{l,k}$ ,

其中，不失一般性，假设结合矩阵和的特殊结构，可得：

$b_{l,k}=2{\left({\tilde{Q}}_{l,k}^{\left(i\right)}\tilde{s}\right)}^H{\tilde{r}}_{l,k}$ ,

经过数学推导简化，得到如下形式：

$C_2(x_k^{\left(i\right)},p_{0,k})=\underset{l=1}{\overset{L}{\&Sigma;}}\left|\right|{\tilde{r}}_{l,k}-2{\left({\tilde{Q}}_{l,k}^{\left(i\right)}\tilde{s}\right)}^H{\tilde{r}}_{l,k}{\tilde{Q}}_{l,k}^{\left(i\right)}\tilde{s}|{|}^2=\underset{l=1}{\overset{L}{\&Sigma;}}|\left|{\tilde{r}}_{l,k}\right|{|}^2-\sqrt{2}|{\left({\tilde{Q}}_{l,k}^{\left(i\right)}\tilde{s}\right)}^H{\tilde{r}}_{l,k}{|}^2$ ,

上式的最小化等价于式子右端的最大化，即构造代价函数为：

$C_3(x_k^{\left(i\right)},p_{0,k})=\underset{l=1}{\overset{L}{\&Sigma;}}|{\left({\tilde{Q}}_{l,k}^{\left(i\right)}\tilde{s}\right)}^H{\tilde{r}}_{l,k}{|}^2={\tilde{s}}^H{Q}_{c,k}^{\left(i\right)}\tilde{s}$ ，

其中，

$Q_{c,k}^{\left(i\right)}=V_k^{\left(i\right)}{V}_k^{\left(i\right)H}$ ，

数据信息矩阵表示为：

$V_k^{\left(i\right)}=\&lsqb;({\tilde{F}}_{1,k}^{\left(i\right)H}{\tilde{A}}_{1,k}^{\left(i\right)H},{\tilde{B}}_1^H){\tilde{r}}_{1,k},({\tilde{F}}_{2,k}^{\left(i\right)H}{\tilde{A}}_{2,k}^{\left(i\right)H},{\tilde{B}}_2^H){\tilde{r}}_{2,k},\mathrm{...},({\tilde{F}}_{L,k}^{\left(i\right)H}{\tilde{A}}_{L,k}^{\left(i\right)H},{\tilde{B}}_L^H){\tilde{r}}_{L,k}\&rsqb;;$

步骤6.数据信息矩阵的最大特征值作为粒子滤波的后验概率权重；根据粒子滤波理论，设定粒子数量及初始权重，通过迭代更新粒子的后验概率权重并进行重采样，所得最终粒子的均值作为该时刻目标跟踪的结果；取粒子对应数据信息矩阵对应的最大特征值为：

$C_3(x_k^{\left(i\right)},p_{0,k})={\mathrm{\&lambda;}}_{max}\left(Q_{c,k}^{\left(i\right)}\right)$ ，

其中，λ_max(·)表示求最大特征值；的维度为(2N+1)×(2N+1)，考虑给定矩阵X，X^HX与XX^H的非零特征值是一致的，故将变换为，

式中， 的维度为L×L维且仅与观测站数量有关，利用该代价函数计算粒子滤波加权，即：

${\tilde{w}}_k\&Proportional;\exp \left(\frac{1}{\det \left({\tilde{C}}_w\right)}{\mathrm{\&lambda;}}_{max}\right({\stackrel{\&OverBar;}{Q}}_{c,k}^{\left(i\right)}\left)\right)$ ，

通过迭代更新粒子后验概率权重并重采样，所得最终粒子的均值即为该时刻目标跟踪的结果。

参见图3～5所示，结合具体的试验数据对本发明做进一步解释说明：

如图3所示，目标初始位置p_0,1＝[1,1]^Tkm，速度v_k＝[0.1,0.05]^Tkm/s，外辐射源位置p_e＝[-5,4]^Tkm，6个观测站的位置分别为[10,9]^Tkm，[7,2]^Tkm，[-6,-8]^Tkm，[6,10]^Tkm，[-2,6]^Tkm，[0,-4]^Tkm。目标初始位置，外辐射源及观测站的地理几何关系如图所示。发射信号载频为f_c＝0.5MHz，带宽为100kHz的高斯信号，其传播速度为c＝3×10⁵km/s；观测站采样频率f_c＝2MHz，跟踪模型中T_r＝1，观测次数K＝30，每次观测时间为32us，系统参数q＝0.032，产生粒子数1000；各站接收回波信噪比10dB，接收直达波信噪比70dB。。

下面将本专利公开的多目标直接定位方法与传统的两步定位方法，以及Weiss-Amar方法进行性能比较，这里的两步定位方法是指利用多重信号分类估计算法(即经典MUSIC算法)进行到达角度估计，然后基于Taylor级数迭代定位算法估计目标位置。

图3和图4分别为本发明基于直接数据域的跟踪方法和传统的基于时差与频差信息的跟踪方法跟踪效果图与误差对比图，其中，基于时差与频差信息方法为，首先采用最大似然法，以第1观测站为参考站估计目标的时差与频差信息，再利用参数信息构建后验概率加权，通过粒子滤波实现目标跟踪。

由图可以看出，采用本发明直接数据方法的估计精度高于传统的方法；在目标运动过程中，直接数据方法能够较好的对目标进行跟踪保持，不同时刻下的误差均较小，而先进行参数估计再利用滤波的传统方法随着目标的运动估计结果起伏较大，其原因就是由于参数误差对跟踪性能的影响。

为了进一步说明直接数据域算法较之传统算法的优势，选择观测间隙k＝15时，对比两种算法在不同信噪比下对目标的定位性能，不同信噪比条件下进行200次蒙特卡洛仿真，结果如图5所示。由图可以看出，在k＝15时刻基于直接数据域的方法定位精度远高于传统两步定位方法，且对信噪比较为稳健。相比于传统的两步定位方法，其非线性门限也高于本专利所提出的直接数据域方法，也是由于其精度同时受到参数估计方法和定位方法的约束；本专利在低信噪比条件下的优势更加明显。

本发明并不局限于上述具体实施方式，本领域技术人员还可据此做出多种变化，但任何与本发明等同或者类似的变化都应涵盖在本发明权利要求的范围内。

Claims (6)

1.一种基于外辐射源联合时延与多普勒频率的直接跟踪方法，具体包含如下步骤：

步骤1.外辐射源场景下利用状态转移矩阵F构建目标运动状态方程为：

x_k＝Fx_k-1+v_k-1，

其中，为状态向量，x_k,y_k与表分别示在第k次观测间隙内目标辐射源的位置坐标和相应的速度分量；T_r表示采样周期，状态转移矩阵F为：

$F=\left[\begin{array}{cccc}1& T_r& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& T_r\\ 0& 0& 0& 1\end{array}\right]$ ，

并定义初始状态与跟踪状态时的条件概率密度函数；结合粒子滤波方法构建后验概率加权迭代方程为：

$w_k^{\left(i\right)}={\tilde{w}}_k^{\left(i\right)}/\underset{i=1}{\overset{M}{\&Sigma;}}{\tilde{w}}_k^{\left(i\right)},{\tilde{w}}_k^{\left(i\right)}=w_{k-1}^{\left(i\right)}p\left(z_k\right|x_k^{\left(i\right)})$ ，

其中，z_k表示k时刻的观测量；

步骤2.对L个观测站的双通道接收系统进行时间同步，并根据Nyquist采样定理采集外辐射源的直达波信号以及经目标反射的回波信号，获得多站接收的信号时域模型；

步骤3.针对信号时域模型，对各站双通道接收的时域数据分别计算其傅立叶系数，得到阵列信号频域模型；

步骤4.针对阵列信号频域模型，每个观测站将所获得的阵列信号频域数据传输至中心站，中心站将每个站传输的阵列信号频域数据按照观测站的顺序进行堆栈排列，构造高维阵列信号频域模型；

步骤5.针对高维阵列信号频域模型，在中心站对高维阵列信号频域数据构造高斯最大似然函数，并构造包含回波时延、多普勒信息以及直达波时延信息的数据信息矩阵；

步骤6.数据信息矩阵的最大特征值作为粒子滤波的后验概率权重；根据粒子滤波理论，设定粒子数量及初始权重，通过迭代更新粒子的后验概率权重并进行重采样，所得最终粒子的均值作为该时刻目标跟踪的结果。

2.根据权利要求1所述的基于外辐射源联合时延与多普勒频率的直接跟踪方法，其特征在于：步骤2中，第l个k时刻观测站的所接收到的信号时域模型为：

$r_{l,k}\left(t\right)=\left[\begin{array}{c}{b}_{d_l}s(t-{\mathrm{\&tau;}}_{d_l})\\ b_{l,k}s(t-{\mathrm{\&tau;}}_{l,k})e^{j2{\mathrm{\&pi;f}}_{l,k}t}\end{array}\right]+\left[\begin{array}{c}{w}_{d_l}\left(t\right)\\ w_{l,k}\left(t\right)\end{array}\right],0\&le;t<T;$

其中，表示外辐射源相对于第l个观测站的直达路径时延；

τ_l,k＝(||p_e-p_0,k||+||p_l-p_0,k||)/c表示外辐射源照射k时刻的目标并反射至观测站产生的时延，c表示信号传播速度，||·||表示2范数；p_e为外辐射源位置，发射信号带宽W，p_0,k为目标在k时刻时的位置，速度为v_k＝[v_x,k,v_y,k]^T；与w_l,k(t)分别表示均值为0，方差为与直达波通道与k时刻回波通道的加性平稳复高斯白噪声；b_l,k与分别表示k时刻回波信道与直达波信道的衰减系数；f_l,k表示目标与观测站之间的多普勒频率，其包含两部分，一部分为外辐射源照射至目标时信号的多普勒频率，另一部分是k时刻反射回波到达观测站时产生的多普勒频率，故f_l,k表示为：

$f_{l,k}\stackrel{\mathrm{\&Delta;}}{=}\frac{f_c}{c}(\frac{v_k^T(p_{0,k}-p_l)}{\left|\right|p_{0,k}-p_l\left|\right|}+\frac{v_k^T(p_{0,k}-p_e)}{\left|\right|p_{0,k}-p_e\left|\right|})$ ，

其中，p_l＝[x_l,y_l]^T(l＝1,2,…,L)为第l个观测站的位置坐标。

3.根据权利要求2所述的基于外辐射源联合时延与多普勒频率的直接跟踪方法，其特征在于：步骤3中第l个观测站k时刻的所接收到的信号频域模型为：

${\tilde{r}}_{l,k}\left(f_n\right)=\left[\begin{array}{c}{b}_{d_l}\tilde{s}\left(f_n\right)e^{-j2{\mathrm{\&pi;f}}_n{\mathrm{\&tau;}}_{d_l}}\\ b_{l,k}\tilde{s}(f_n-f_{l,k})e^{-j2{\mathrm{\&pi;f}}_n{\mathrm{\&tau;}}_{l,k}}\end{array}\right]+\left[\begin{array}{c}{\tilde{w}}_{d_l}\left(f_n\right)\\ {\tilde{w}}_{l,k}\left(f_n\right)\end{array}\right]$ ，

其中，与分别表示接收信号、发射信号以及噪声项的傅立叶系数。

4.根据权利要求3所述的基于外辐射源联合时延与多普勒频率的直接跟踪方法，其特征在于：步骤4中，中心站所获得的高维信号频域模型为：

${\tilde{r}}_{l,k}=\left[\begin{array}{c}{b}_{d_l}{\tilde{B}}_l\\ b_{l,k}{\tilde{A}}_{l,k}{\tilde{F}}_{l,k}\end{array}\right]\tilde{s}+\left[\begin{array}{c}{\tilde{w}}_{d_l}\\ {\tilde{w}}_{l,k}\end{array}\right]$ ，

其中，

$\left\{\begin{array}{c}{\tilde{r}}_{l,k}={\&lsqb;{\tilde{r}}_{l,k}\left(f_{-2N}\right),{\tilde{r}}_{l,k}\left(f_{-2N+1}\right),\mathrm{...},{\tilde{r}}_{l,k}\left(f_{2N}\right)\&rsqb;}^T\\ \tilde{s}={\&lsqb;\tilde{s}\left(f_{-N}\right),\tilde{s}\left(f_{-N+1}\right),\mathrm{...},\tilde{s}\left(f_N\right)\&rsqb;}^T\\ {\tilde{B}}_l=diag\{e^{-j2{\mathrm{\&pi;f}}_{-N}{\mathrm{\&tau;}}_{d_l}},e^{-j2{\mathrm{\&pi;f}}_{-N+1}{\mathrm{\&tau;}}_{d_l}},\mathrm{...},e^{-j2{\mathrm{\&pi;f}}_N{\mathrm{\&tau;}}_{d_l}}\}\\ {\tilde{A}}_{l,k}=diag\{e^{-j2{\mathrm{\&pi;f}}_{-N}{\mathrm{\&tau;}}_{l,k}},e^{-j2{\mathrm{\&pi;f}}_{-N+1}{\mathrm{\&tau;}}_{l,k}},\mathrm{...},e^{-j2{\mathrm{\&pi;f}}_N{\mathrm{\&tau;}}_{l,k}}\}\\ {\tilde{w}}_{l,k}={\&lsqb;{\tilde{w}}_{l,k}\left(f_{-N}\right),{\tilde{w}}_{l,k}\left(f_{-N+1}\right),\mathrm{...},{\tilde{w}}_{l,k}\left(f_N\right)\&rsqb;}^T\\ {\tilde{w}}_{d_l}={\&lsqb;{\tilde{w}}_{d_l}\left(f_{-N}\right),{\tilde{w}}_{d_l}\left(f_{-N+1}\right),\mathrm{...},{\tilde{w}}_{d_l}\left(f_N\right)\&rsqb;}^T\end{array}\right.$ ，

式中，表示循环移位矩阵，其形式表示为：

，

即将单位阵循环向下移位floor(Tf_l,k)+1行；即将矩阵按行循环向下移位floor(Tf_l,k)，floor(·)表示向下取整。

5.根据权利要求4所述的基于外辐射源联合时延与多普勒频率的直接跟踪方法，其特征在于：步骤5中，构造的高斯最大似然函数为：

$p\left({\tilde{r}}_k\right|x_k^{\left(i\right)})\&Proportional;\exp \left(-\frac{1}{\det \left({\tilde{C}}_w\right)}\underset{l=1}{\overset{L}{\&Sigma;}}\left|\right|{\tilde{r}}_{l,k}-b_l\left(\begin{array}{c}{\tilde{A}}_{l,k}^{\left(i\right)}{\tilde{F}}_{l,k}^{\left(i\right)}\\ {\tilde{B}}_l\end{array}\right)\tilde{s}|{|}^2\right)$ ，

其中，det(·)表示求行列式； 为了估计目标的位置坐标，应使上式右端最小化，此时第i个粒子似然表示如下：

$p\left(r_k\right|x_k^{\left(i\right)})\&Proportional;\exp (-{\frac{1}{\det \left({\tilde{C}}_w\right)}}_p\underset{p_{0,k}}{\arg \min }{C}_1(x_k^{\left(i\right)},p_{0,k}\left)\right)$ ，

其中，代价函数为：

$C_1(x_k^{\left(i\right)},p_{0,k})=\underset{l=1}{\overset{L}{\&Sigma;}}\left|\right|{\tilde{r}}_{l,k}-b_{l,k}\left(\begin{array}{c}{\tilde{A}}_{l,k}^{\left(i\right)}{\tilde{F}}_{l,k}^{\left(i\right)}\\ {\tilde{B}}_l\end{array}\right)\tilde{s}|{|}^2$

利用b_l,k使得上式最小化，可得：

$b_{l,k}=\&lsqb;{\left(\begin{array}{c}{\tilde{A}}_{l,k}^{\left(i\right)}{\tilde{F}}_{l,k}^{\left(i\right)}\\ {\tilde{B}}_{l,k}^{\left(i\right)}\end{array}\right)}^H{\tilde{s}}^H\tilde{s}\left(\begin{array}{c}{\tilde{A}}_{l,k}^{\left(i\right)}{\tilde{F}}_{l,k}^{\left(i\right)}\\ {\tilde{B}}_{l,k}^{\left(i\right)}\end{array}\right)\&rsqb;{\left({\tilde{Q}}_{l,k}^{\left(i\right)}\tilde{s}\right)}^H{\tilde{r}}_{l,k}$ ，

其中，不失一般性，假设结合矩阵和的特殊结构，可得：

$b_{l,k}=2{\left({\tilde{Q}}_{l,k}^{\left(i\right)}\tilde{s}\right)}^H{\tilde{r}}_{l,k}$ ，

经过数学推导简化，得到如下形式：

$C_2(x_k^{\left(i\right)},p_{0,k})=\underset{l=1}{\overset{L}{\&Sigma;}}\left|\right|{\tilde{r}}_{l,k}-2{\left({\tilde{Q}}_{l,k}^{\left(i\right)}\tilde{s}\right)}^H{\tilde{r}}_{l,k}{\tilde{Q}}_{l,k}^{\left(i\right)}\tilde{s}|{|}^2=\underset{l=1}{\overset{L}{\&Sigma;}}|\left|{\tilde{r}}_{l,k}\right|{|}^2-\sqrt{2}|{\left({\tilde{Q}}_{l,k}^{\left(i\right)}\tilde{s}\right)}^H{\tilde{r}}_{l,k}{|}^2$ ，

上式的最小化等价于式子右端的最大化，即构造代价函数为：

$C_3(x_k^{\left(i\right)},p_{0,k})=\underset{l=1}{\overset{L}{\&Sigma;}}|{\left({\tilde{Q}}_{l,k}^{\left(i\right)}\tilde{s}\right)}^H{\tilde{r}}_{l,k}{|}^2={\tilde{s}}^H{Q}_{c,k}^{\left(i\right)}\tilde{s}$ ，

其中，

$Q_{c,k}^{\left(i\right)}=V_k^{\left(i\right)}{V}_k^{\left(i\right)H}$ ，

数据信息矩阵表示为：

$V_k^{\left(i\right)}=\&lsqb;({\tilde{F}}_{1,k}^{\left(i\right)H}{\tilde{A}}_{1,k}^{\left(i\right)H},{\tilde{B}}_1^H){\tilde{r}}_{1,k},({\tilde{F}}_{2,k}^{\left(i\right)H}{\tilde{A}}_{2,k}^{\left(i\right)H},{\tilde{B}}_2^H){\tilde{r}}_{2,k},...,({\tilde{F}}_{L,k}^{\left(i\right)H}{\tilde{A}}_{L,k}^{\left(i\right)H},{\tilde{B}}_L^H){\tilde{r}}_{L,k}\&rsqb;.$

6.根据权利要求5所述的基于外辐射源联合时延与多普勒频率的直接跟踪方法，其特征在于：步骤6中具体包含如下：取粒子对应数据信息矩阵对应的最大特征值为：

$C_3(x_k^{\left(i\right)},p_{0,k})={\mathrm{\&lambda;}}_{\max }\left(Q_{c,k}^{\left(i\right)}\right)$ ，

其中，λ_max(·)表示求最大特征值；的维度为(2N+1)×(2N+1)，考虑给定矩阵X，X^HX与XX^H的非零特征值是一致的，故将变换为，

式中， 的维度为L×L维且仅与观测站数量有关，利用该代价函数计算粒子滤波加权，即：

${\tilde{w}}_k\&Proportional;\exp \left(\frac{1}{\det \left({\tilde{C}}_w\right)}{\mathrm{\&lambda;}}_{max}\right({\stackrel{\&OverBar;}{Q}}_{c,k}^{\left(i\right)}\left)\right)$ ，

通过迭代更新粒子后验概率权重并重采样，所得最终粒子的均值即为该时刻目标跟踪的结果。