CN103902826A - 一种冲击噪声环境下的多移动目标跟踪方法 - Google Patents

Abstract

本发明的目的在于提供一种冲击噪声环境下的多移动目标跟踪方法。该方法先设置特殊的非均匀线阵的阵列结构，将目标锁定在一个变化的搜索范围之内，再设计量子文化雁群方法，利用其搜索机制在搜索空间搜索扩展加权信号极大似然方程的最优角度值。通过逐渐减小搜索的范围以及运用智能搜索机制，有效解决了搜索方法的计算量问题。仿真结果表明这种冲击噪声环境下的多移动目标跟踪方法能够保证所设计方法的实时性，而且具备阵列扩展的能力和较好的跟踪精度，在强冲击噪声等恶劣噪声环境下同样具有较好的性能。

Description

一种冲击噪声环境下的多移动目标跟踪方法

技术领域

本发明涉及的是一种移动目标跟踪方法。

背景技术

测向通常又称空间谱估计或DOA估计，是阵列信号处理的一个重要研究领域，在无源探测、雷达、声纳和通信等领域都有广泛的应用。经过几十年的发展，基于高斯模型假设的DOA估计理论已经逐步走向完善和成熟。但实际应用中，所遇到的许多随机信号和噪声并不是高斯分布的，例如大气雷电噪声、通信线路上的瞬间尖峰语音信号和海洋环境噪声以及多种人为噪声等，这些信号中存在显著的尖峰，都可用特征指数α不同的SaS过程来描述，用传统的基于二阶和四阶统计量的方法进行处理不能得到满意的结果。如何在冲击噪声背景下找到既能改善基于分数低阶矩测向算法的测向性能，同时又可以扩展阵列孔径，且能实现对动态目标的相干信源进行测向的DOA估计算法，是工程应用中遇到的难点。对这些问题的有效解决，才可进一步实现提高雷达和无线通信系统设计的质量和效率，才能使测向技术有更强的实用性。

经对现有技术文献的检索发现，刁鸣等在《系统工程与电子技术》(2009,Vol.29,No.12,pp.2046–2049)上发表的“一种新的基于粒子群算法的DOA跟踪方法”中使用了粒子群算法和极大似然算法进行多移动目标的动态DOA估计，在高斯噪声背景下效果较好，但在冲击噪声背景下，方法性能恶化严重导致失效。赵大勇等在《山东大学学报》(2010,Vol.40,No.1,pp.133-138)上发表的“冲击噪声背景下的动态DOA跟踪”提出了粒子群算法和分数低阶矩的极大似然方法求解弱冲击噪声环境下的动态DOA估计问题，使用分数低阶矩的动态更新能在一定程度上实现动态DOA估计问题，在信源数大于阵元数的情况下算法就失效了，而且在特征指数小于1的强冲击噪声情况下基于粒子群算法和分数低阶矩的极大似然方法就失效了。

已有文献表明，在冲击噪声等恶劣噪声环境下，动态跟踪的鲁棒性问题依旧没有解决。针对冲击噪声环境下的动态跟踪是一个难题，尤其对大于阵元数的多信源进行动态跟踪更难，故不能把高斯噪声环境下的DOA估计算法直接进行移植。在冲击噪声下的鲁棒动态DOA估计首先应能建立高性能的鲁棒动态跟踪的优化方程，冲击噪声环境下经典智能计算方法很难摆脱收敛速度和收敛性能矛盾的制约，在现有计算条件下很难在有限的时间内搜索到最优解，需要设计新的智能算法求解冲击噪声环境特别是冲击噪声环境下的鲁棒动态跟踪问题。

但使用极大似然算法进行角度估计的一个主要缺点是搜索过程是一个复杂耗时的过程，运算量巨大。

发明内容

本发明的目的在于提供在恶劣测向背景达到鲁棒跟踪的一种冲击噪声环境下的多移动目标跟踪方法。

本发明的目的是这样实现的：

本发明一种冲击噪声环境下的多移动目标跟踪方法，其特征是：

（1）根据天线个数和测向要求确定天线位置，阵列远场在θ_k(k=1,2,…,N)处有N个窄带点源以平面波入射，入射波长为λ，入射到由M个各向同性的阵元构成的非均匀特殊线阵上，阵元的放置满足如下的条件：第i个阵元相对于第一个阵元的间距为d₁=0<d₂<…<d_M，其值都是半波长的整数倍，并满足Ω={d_i-d_j|i,j=1,2,…,M}是完整的，即Ω={-d_M,…,0,…,d_M}，阵列接收的快拍数据可表示为x(k)=A(θ)s(k)+n(k)，式中x(k)=[x₁(k),x₂(k),…,x_M(k)]^T为M×1维阵列快拍数据矢量，s(t)=[s₁(t),s₂(t),…,s_N(t)]^T为N×1维非高斯信号矢量，θ=(θ₁,θ₂,…,θ_N)为信源方位矢量，A(θ)=[a(θ₁)a(θ₂)…a(θ_N)]是信号导向矢量矩阵，其中第i个导向矢量 $a_i={[1,e^{-j2\mathrm{\&pi;}{d}_1\sin \left({\mathrm{\&theta;}}_i\right)/\mathrm{\&lambda;}},\&CenterDot;\&CenterDot;\&CenterDot;,e^{-j2\mathrm{\&pi;}{d}_M\sin \left({\mathrm{\&theta;}}_i\right)/\mathrm{\&lambda;}}]}^T,$ (i=1,2,…,N)；第k次采样数据的无穷范数归一化信号 $\stackrel{\&OverBar;}{x}\left(k\right)={[{\stackrel{\&OverBar;}{x}}_1\left(k\right),{\stackrel{\&OverBar;}{x}}_2\left(k\right),\&CenterDot;\&CenterDot;\&CenterDot;,{\stackrel{\&OverBar;}{x}}_M\left(k\right)]}^T=x\left(k\right)/\underset{1\&le;j\&le;M}{\max }\left\{\right|x_j\left(k\right)\left|\right\},$ 其分数低阶协方差为C(k)，其中第i维第j列元素为 $C_{\mathrm{ij}}\left(k\right)={\stackrel{\&OverBar;}{x}}_i\left(k\right){\left|{\stackrel{\&OverBar;}{x}}_i\left(k\right)\right|}^{p-1}{\stackrel{\&OverBar;}{x}}_j^{*}\left(k\right){\left|{\stackrel{\&OverBar;}{x}}_j\left(k\right)\right|}^{p-1},$ *为共轭算子，p为分数低阶协方差参数，把非均匀线阵可虚拟成更多个阵元的均匀线阵，阵列的最大相关延迟为M_a，则虚拟出的均匀直线阵阵元个数为M+1～M_a+1，设

则阵元坐标矢量为d=[d₁,d₂,…,d_M]=ε[n₁,…,n_M]，n₁,…,n_M都为整数，则接收阵列的分数低阶协方差C(k)（M×M维）可以写作 $C\left(k\right)=\left[\begin{array}{ccccc}{C}_{11}^0\left(k\right)& C_{12}^{(d_1-d_2)}\left(k\right)& C_{13}^{(d_1-d_3)}\left(k\right)& ...& C_{1M}^{(d_1-d_M)}\left(k\right)\\ C_{21}^{(d_2-d_1)}\left(k\right)& C_{22}^0\left(k\right)& C_{23}^{(d_2-d_3)}\left(k\right)& ...& C_{2M}^{(d_2-d_M)}\left(k\right)\\ C_{31}^{d_3-d_1}\left(k\right)& C_{32}^{d_3-d_2}\left(k\right)& C_{33}^0\left(k\right)& ...& C_{3M}^{(d_3-d_M)}\left(k\right)\\ ...& ...& ...& ...& ...\\ C_{M1}^{d_M-d_1}\left(k\right)& C_{M2}^{d_M-d_2}\left(k\right)& C_{M3}^{d_M-d_3}\left(k\right)& ...& C_{\mathrm{MM}}^0\left(k\right)\end{array}\right],$ 若令 ${\stackrel{\&OverBar;}{C}}^0\left(k\right)=\frac{1}{M}\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}{C}_{\mathrm{ii}}^0\left(k\right),{\stackrel{\&OverBar;}{C}}^{\mathrm{m\&epsiv;}}\left(k\right)=E\left[C_{\mathrm{ij}}^{\mathrm{m\&epsiv;}}\left(k\right)\right],{\stackrel{\&OverBar;}{C}}^{-\mathrm{m\&epsiv;}}\left(k\right)=E\left[C_{\mathrm{ij}}^{-\mathrm{m\&epsiv;}}\left(k\right)\right],m=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;,M_a,$ E()为取均值函数，虚拟均匀线阵的扩展分数低阶协方差增量

（(M_a+1)×(M_a+1)维）可以表示为 $\stackrel{\&OverBar;}{C}=\left[\begin{array}{ccccc}{\stackrel{\&OverBar;}{C}}^0\left(k\right)& {\stackrel{\&OverBar;}{C}}^{-\mathrm{\&epsiv;}}\left(k\right)& {\stackrel{\&OverBar;}{C}}^{-2\mathrm{\&epsiv;}}\left(k\right)& ...& {\stackrel{\&OverBar;}{C}}^{-M_a\mathrm{\&epsiv;}}\left(k\right)\\ {\stackrel{\&OverBar;}{C}}^{\mathrm{\&epsiv;}}\left(k\right)& {\stackrel{\&OverBar;}{C}}^0\left(k\right)& {\stackrel{\&OverBar;}{C}}^{-\mathrm{\&epsiv;}}\left(k\right)& ...& {\stackrel{\&OverBar;}{C}}^{-(M_a-1)\mathrm{\&epsiv;}}\left(k\right)\\ {\stackrel{\&OverBar;}{C}}^{2\mathrm{\&epsiv;}}\left(k\right)& {\stackrel{\&OverBar;}{C}}^{\mathrm{\&epsiv;}}\left(k\right)& {\stackrel{\&OverBar;}{C}}^0\left(k\right)& ...& {\stackrel{\&OverBar;}{C}}^{-(M_a-2)\mathrm{\&epsiv;}}\left(k\right)\\ ...& ...& ...& ...& ...\\ {\stackrel{\&OverBar;}{C}}^{M_a\mathrm{\&epsiv;}}\left(k\right)& {\stackrel{\&OverBar;}{C}}^{(M_a-1)\mathrm{\&epsiv;}}\left(k\right)& {\stackrel{\&OverBar;}{C}}^{(M_a-2)\mathrm{\&epsiv;}}\left(k\right)& ...& {\stackrel{\&OverBar;}{C}}^0\left(k\right)\end{array}\right]=\left[\begin{array}{ccccc}{R}_{11}\left(k\right)& R_{12}\left(k\right)& R_{13}\left(k\right)& ...& R_{1(M_a+1)}\left(k\right)\\ R_{21}\left(k\right)& R_{22}\left(k\right)& R_{23}\left(k\right)& ...& R_{2(M_a+1)\left(k\right)}\\ R_{31}\left(k\right)& R_{32}\left(k\right)& R_{33}\left(k\right)& ...& R_{3(M_a+1)\left(k\right)}\\ ...& ...& ...& ...& ...\\ R_{(M_a+1)1}\left(k\right)& R_{(M_a+1)2}\left(k\right)& R_{(M_a+1)3}\left(k\right)& ...& R_{(M_a+1)(M_a+1)\left(k\right)}\end{array}\right]=R$ ，则虚拟成均匀线阵扩展了原阵列的阵列流型矢量为B(θ)=[b(θ₁),b(θ₂),…,b(θ_N)],其第i个导向矢量 $b\left({\mathrm{\&theta;}}_i\right)={[1,e^{-j2\mathrm{\&pi;\&epsiv;}\sin \left({\mathrm{\&theta;}}_i\right)/\mathrm{\&lambda;}},\&CenterDot;\&CenterDot;\&CenterDot;,e^{-j2\mathrm{\&pi;}\left(M_a\right)\mathrm{\&epsiv;}\sin \left({\mathrm{\&theta;}}_i\right)/\mathrm{\&lambda;}}]}^T,i=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;,N,$ 根据非均匀线阵虚拟出M_a>M个阵元，则可测信号数为M_a-1，第一次快拍采样的扩展分数低阶协方差

（M_a×M_a维）等于其扩展分数低阶协方差增量R(1),t代表量子文化雁群搜索机制的迭代次数，初始时设t=0；

（2）确定雁群中所有大雁初始状态，由H个大雁组成的雁群，每只大雁在N维搜索空间中运动，第i只大雁的量子位置定义如下:

y_i(t)=[y_i1(t),y_i2(t),…,y_iN(t)]表示第i只大雁量子位置的第n维量子位，(i=1,2,…,H)，0≤y_in(t)≤1(n=1,2,…,N)，把第i只大雁的量子位置映射到定义区间，就是该大雁的位置映射关系为 ${\stackrel{\&OverBar;}{y}}_{\mathrm{in}}\left(t\right)=y_{\mathrm{in}}\left(t\right)[h_n\left(k\right)-b_n\left(k\right)]+b_n\left(k\right),n=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;,N,$ h_n(k)和b_n(k)分别为第k次采样数据处理时第n维角度搜索区间的上限和下限值，第i只大雁对应的速度为v_i(t)=[v_i1(t),v_i2(t),…,v_iN(t)]，-0.2≤v_in(t)≤0.2(n=1,2,…,N)表示第i只大雁的第n维速度，量子位置在定义的量子域[0，1]随机初始化，速度在[-0.2，0.2]随机初始化，确定初始搜索空间确定初始搜索空 $Z\left(k\right)=\left[\begin{array}{cccc}{h}_1\left(k\right)& h_2\left(k\right)& ...& h_N\left(k\right)\\ b_1\left(k\right)& b_2\left(k\right)& ...& b_N\left(k\right)\end{array}\right]$ 为N个角度的搜索空间向量，最大迭代次数设置为

其中

为取整函数，T为取整倍数；

（3）第i只大雁位置的扩展分数低阶协方差极大似然方程的适应度值为

正交投影矩阵 $P_{B\left[{\stackrel{\&OverBar;}{y}}_i\left(t\right)\right]}=B\left({\stackrel{\&OverBar;}{y}}_i\left(t\right)\right){\left[B^H\left({\stackrel{\&OverBar;}{y}}_i\left(t\right)\right)B\left({\stackrel{\&OverBar;}{y}}_i\left(t\right)\right)\right]}^{-1}{B}^H\left({\stackrel{\&OverBar;}{y}}_i\left(t\right)\right),$ 上标H代表共轭的转置，tr()表示矩阵的求迹运算，根据

评价大雁位置和对应量子位置的优劣，第i只大雁到现在为止所经历的最优量子位置定义为该大雁的局部最优量子位置，记作p_i(t)=[p_i1(t),p_i2(t),…,p_iN(t)]，p_in(t)为第i只大雁到第t次迭代为止所经历第n维最优量子位，n=1,2,…,N，所有大雁到现在所经历的最优量子位置记作全局最优量子位置，也就是适应度值最大的量子位置，记作g(t)=[g₁(t),g₂(t),…,g_N(t)]，g_n(t)为所有大雁到第t次迭代为止所经历第n维最优量子位，n=1,2,…,N，第n维量子形势知识用I_n=[l_n,u_n]表示，下限l_n和上限u_n根据问题所给定的量子位置定义域的边界进行初始化；L_n表示第n维量子位置的下限l_n所对应的适应度值，U_n表示第n维量子位置的上限u_n所对应的适应度值，在最大值优化问题中均可初始化为-∞；

（4）更新每只大雁的速度和量子位置，H只大雁速度先根据群行为进行更新，则第i只大雁的第n维速度更新为v_in(t+1)=c₁r₁[g_n(t)-y_in(t)]+c₂r₂[p_(i-_1)n(t)-y_in(t)]+w_tv_in(t)，w_t从初次迭代时为0.9线形递减到迭代次数最大时为0.5，i=1,2,…,H，n=1,2,…,N，r₁和r₂都是[0，1]之间的均匀随机数，c₁和c₂加权常数，对于v_in(t+1)，若超出边界值，将其限制在边界，即v_in(t+1)>0.2，v_in(t+1)=0.2，若v_in(t+1)<-0.2，

v_in(t+1)=-0.2，第i只大雁的新量子位置为

y_i(t+1)=[y_i1(t+1),y_i2(t+1),…,y_iN(t+1)]，i=1,2,…,H，

$y_{\mathrm{in}}(t+1)=\mathrm{abs}[y_{\mathrm{in}}\left(t\right)\cos (-v_{\mathrm{in}}(t+1))-\sqrt{1-{\left(y_{\mathrm{in}}\left(t\right)\right)}^2}\sin (-v_{\mathrm{in}}(t+1))],$

n=1,2,…,N，abs()为量子位取绝对值函数；

（5）使用量子规范知识指导更新量子位置为y_i(t+1)，i=H+1,H+2,…,2H，第i只大雁的第n维的量子旋转角为

n=1,2,…,N，G(0,1)是均值为0方差为1的高斯随机数，η为加权常数；

（6）将第i只大雁的量子位置y_i(t+1)影射到定义区间

计算适应度值，适应度函数为 $f\left[{\stackrel{\&OverBar;}{y}}_i(t+1)\right]=\mathrm{tr}\left[P_{B\left[{\stackrel{\&OverBar;}{y}}_i(t+1)\right]}\stackrel{\&OverBar;}{R}\left(k\right)\right],i=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;,2H,$ 对于第i只大雁，若

则令y_i(t+1)=y_i+N(t+1)，若p_i(t+1)=y_i(t+1)，否则，p_i(t+1)=p_i(t)，把所有大雁的局部最优量子位置按照适应度值由大到小进行排序，把最优即适应度值最大的局部最优量子位置设置成全局最优量子位置为g(t+1)=[g₁(t+1),g₂(t+1),…,g_N(t+1)]；

从所有大雁局部最优位置中选取最优的前εH个局部最优量子位置去更新量子规范知识，设第j只大雁影响第i维量子规范知识的下限，第m只大雁影响第i维量子规范知识的上限，通过下面的量子规范知识更新函数对第n维量子规范知识进行更新：

其中，

分别表示第t代第n维量子变量的下限

和上限

所对应适应度；

分别为第j只大雁和第m只大雁的局部最优量子位置p_j和p_k映射到定义域区间的潜在解向量；

（7）判断是否达到最大迭代次数，若是，记录最优位置，执行步骤（8）；否则t=t+1，返回步骤（4）进行迭代；

（8）快拍采样新数据x(k+1)=[x₁(k+1),x₂(k+1),…,x_M(k+1)]^T的无穷范数归一化信号为 $\stackrel{\&OverBar;}{x}(k+1)={[{\stackrel{\&OverBar;}{x}}_1(k+1),{\stackrel{\&OverBar;}{x}}_2(k+1),\&CenterDot;\&CenterDot;\&CenterDot;,{\stackrel{\&OverBar;}{x}}_M(k+1)]}^T=x(k+1)/\underset{1\&le;j\&le;M}{\max }\left\{\right|x_j(k+1)\left|\right\},$ 则 $C_{\mathrm{ij}}(k+1)={\stackrel{\&OverBar;}{x}}_i(k+1){\left|{\stackrel{\&OverBar;}{x}}_i(k+1)\right|}^{p-1}{\stackrel{\&OverBar;}{x}}_j^{*}(k+1){\left|{\stackrel{\&OverBar;}{x}}_j(k+1)\right|}^{p-1}$ 是新增加的第k+1个采样数据的分数低阶协方差增量矩阵，依照步骤（1）的方法，根据分数低阶协方差增量矩阵构造出虚拟均匀线阵的扩展分数低阶协方差增量

（(M_a+1)×(M_a+1)维），且可以表示为 $\stackrel{\&OverBar;}{C}=\left[\begin{array}{ccccc}{\stackrel{\&OverBar;}{C}}^0\left(k\right)& {\stackrel{\&OverBar;}{C}}^{-\mathrm{\&epsiv;}}\left(k\right)& {\stackrel{\&OverBar;}{C}}^{-2\mathrm{\&epsiv;}}\left(k\right)& ...& {\stackrel{\&OverBar;}{C}}^{-M_a\mathrm{\&epsiv;}}\left(k\right)\\ {\stackrel{\&OverBar;}{C}}^{\mathrm{\&epsiv;}}\left(k\right)& {\stackrel{\&OverBar;}{C}}^0\left(k\right)& {\stackrel{\&OverBar;}{C}}^{-\mathrm{\&epsiv;}}\left(k\right)& ...& {\stackrel{\&OverBar;}{C}}^{-(M_a-1)\mathrm{\&epsiv;}}\left(k\right)\\ {\stackrel{\&OverBar;}{C}}^{2\mathrm{\&epsiv;}}\left(k\right)& {\stackrel{\&OverBar;}{C}}^{\mathrm{\&epsiv;}}\left(k\right)& {\stackrel{\&OverBar;}{C}}^0\left(k\right)& ...& {\stackrel{\&OverBar;}{C}}^{-(M_a-2)\mathrm{\&epsiv;}}\left(k\right)\\ ...& ...& ...& ...& ...\\ {\stackrel{\&OverBar;}{C}}^{M_a\mathrm{\&epsiv;}}\left(k\right)& {\stackrel{\&OverBar;}{C}}^{(M_a-1)\mathrm{\&epsiv;}}\left(k\right)& {\stackrel{\&OverBar;}{C}}^{(M_a-2)\mathrm{\&epsiv;}}\left(k\right)& ...& {\stackrel{\&OverBar;}{C}}^0\left(k\right)\end{array}\right]=\left[\begin{array}{ccccc}{R}_{11}\left(k\right)& R_{12}\left(k\right)& R_{13}\left(k\right)& ...& R_{1(M_a+1)}\left(k\right)\\ R_{21}\left(k\right)& R_{22}\left(k\right)& R_{23}\left(k\right)& ...& R_{2(M_a+1)\left(k\right)}\\ R_{31}\left(k\right)& R_{32}\left(k\right)& R_{33}\left(k\right)& ...& R_{3(M_a+1)\left(k\right)}\\ ...& ...& ...& ...& ...\\ R_{(M_a+1)1}\left(k\right)& R_{(M_a+1)2}\left(k\right)& R_{(M_a+1)3}\left(k\right)& ...& R_{(M_a+1)(M_a+1)\left(k\right)}\end{array}\right]=R$ ，通过下面的方式更新当前采样的扩展分数低阶协方差：其中R_ij(k+1)新增加的第k+1个采样数据的扩展分数低阶协方差的第i行第j列增量，更新搜索空间为 $Z(k+1)=\left[\begin{array}{cccc}{h}_1(k+1)& h_2(k+1)& ...& h_N(k+1)\\ b_1(k+1)& b_2(k+1)& ...& b_N(k+1)\end{array}\right],$ $b_n(k+1)={\stackrel{\&OverBar;}{\mathrm{\&theta;}}}_{\mathrm{zn}}\left(k\right)-{\mathrm{\&beta;}}^k|b_n\left(k\right)-{\stackrel{\&OverBar;}{\mathrm{\&theta;}}}_n\left(k\right)|-r,h_n(k+1)={\stackrel{\&OverBar;}{\mathrm{\&theta;}}}_{\mathrm{zn}}\left(k\right)+{\mathrm{\&beta;}}^k|h_n\left(k\right)-{\stackrel{\&OverBar;}{\mathrm{\&theta;}}}_n\left(k\right)|+r,$ n=1,2,…,N，β为收敛因子；常数r为搜索空间在锁定状态下的搜索半径；

为第n个方向在第k次采样时的估计值；

为第n个方向在第k次采样搜索空间的中心值，其更新公式为

其中δ为遗传因子；

（9）如果达到最大快拍采样数，执行步骤（10）；否则设k=k+1，迭代次数t=0，返回步骤（2）继续估计动态目标下一个时刻的方向；

（10）得到所有快拍采样下的一系列全局最优量子位置，其映射得到的全局最优位置就是跟踪到的动态目标方向值，输出动态跟踪结果。

本发明还可以包括：

1、加权常数η的取值为η=0.06。

本发明的优势在于：

(1)本发明适合冲击噪声环境，同时也适应高斯噪声和强冲击噪声环境，可以改善相同阵元数均匀线阵列的动态波达方向角估计性能。

(2)较好地解决了冲击噪声环境测向的阵列扩展问题，在天线数小于阵元数时也可进行有效的跟踪。

(3)根据具体的应用环境通过调整天线位置就可以满足测更多的信源或获得更好的解相干能力，在位置受限的情况下，可以避开一些不可放置天线的点，也可以取得较好的效果。

(4)使用所设计的量子文化雁群搜索方法可快速对扩展分数低阶协方差极大似然方程进行高精度求解。

附图说明

图1为本发明的流程图；

图2为本发明的搜索机制结构示意图；

图3为特征指数α=1.2时对于3个独立信源使用粒子群算法和分数低阶矩极大似然算法的角度跟踪情况；

图4为特征指数α=1.2时对于3个独立信源使用所提的量子文化雁群动态跟踪方法的角度跟踪情况；

图5为特征指数α=1.6时对于4个信源使用粒子群算法和分数低阶矩极大似然算法的角度跟踪情况；

图6为特征指数α=1.6时对于4个信源使用所提的量子文化雁群动态跟踪方法的角度跟踪情况；

图7为特征指数α=1.8时所提量子文化雁群动态跟踪方法对5个独立信源进行动态跟踪情况；

图8为特征指数α=1.9时所提量子文化雁群动态跟踪方法对6个独立信源进行动态跟踪情况；

图9为特征指数α=0.8时使用粒子群算法和分数低阶矩极大似然算法对2个独立信源的进行动态跟踪情况；

图10为特征指数α=0.8时所提量子文化雁群动态跟踪方法对2个独立信源的动态跟踪情况；

图11为量子文化雁群动态跟踪方法的总体结构图。

具体实施方式

下面结合附图举例对本发明做更详细地描述：

结合图1～10，本发明主要包括以下步骤：

步骤一，根据天线个数和测向要求确定天线位置。阵列远场在θ_k(k=1,2,…,N)处有N个窄带点源以平面波入射，入射波长为λ。入射到由M个各向同性的阵元构成的非均匀特殊线阵上，阵元的放置满足如下的条件：第i个阵元相对于第一个阵元的间距为d₁=0<d₂<…<d_M,其值都是半波长的整数倍，并满足Ω={d_i-d_j|i,j=1,2,…,M}是完整的，即Ω={-d_M,…,0,…,d_M}。阵列接收的快拍数据可表示为x(k)=A(θ)s(k)+n(k)，式中，x(k)=[x₁(k),x₂(k),…,x_M(k)]^T为M×1维阵列快拍数据矢量。s(t)=[s₁(t),s₂(t),…,s_N(t)]^T为N×1维非高斯信号矢量，θ=(θ₁,θ₂,…,θ_N)为信源方位矢量，A(θ)=[a(θ₁)a(θ₂)…a(θ_N)]是信号导向矢量矩阵，其中第i个导向矢量 $a_i={[1,e^{-j2\mathrm{\&pi;}{d}_1\sin \left({\mathrm{\&theta;}}_i\right)/\mathrm{\&lambda;}},\&CenterDot;\&CenterDot;\&CenterDot;,e^{-j2\mathrm{\&pi;}{d}_M\sin \left({\mathrm{\&theta;}}_i\right)/\mathrm{\&lambda;}}]}^T,$ 第k次采样数据的无穷范数归一化信号 $\stackrel{\&OverBar;}{x}\left(k\right)={[{\stackrel{\&OverBar;}{x}}_1\left(k\right),{\stackrel{\&OverBar;}{x}}_2\left(k\right),\&CenterDot;\&CenterDot;\&CenterDot;,{\stackrel{\&OverBar;}{x}}_M\left(k\right)]}^T=x\left(k\right)/\underset{1\&le;j\&le;M}{\max }\left\{\right|x_j\left(k\right)\left|\right\},$ 其分数低阶协方差为C(k)，其中第i维第j列元素为 $C_{\mathrm{ij}}\left(k\right)={\stackrel{\&OverBar;}{x}}_i\left(k\right){\left|{\stackrel{\&OverBar;}{x}}_i\left(k\right)\right|}^{p-1}{\stackrel{\&OverBar;}{x}}_j^{*}\left(k\right){\left|{\stackrel{\&OverBar;}{x}}_j\left(k\right)\right|}^{p-1},$ 式中，*为共轭算子，p为分数低阶协方差参数。根据特殊线阵的特点，把非均匀线阵可虚拟成更多个阵元的均匀线阵，阵列的最大相关延迟为M_a，则虚拟出的均匀直线阵阵元个数为M+1～M_a+1。设

则阵元坐标矢量为d=[d₁,d₂,…,d_M]=ε[n₁,…,n_M]，n₁,…,n_M都为整数，为了更好表述接收的分数低阶协方差和虚拟出均匀线阵的扩展分数低阶协方差的关系，则接收阵列的分数低阶协方差C(k)（M×M维）可以写作 $C\left(k\right)=\left[\begin{array}{ccccc}{C}_{11}^0\left(k\right)& C_{12}^{(d_1-d_2)}\left(k\right)& C_{13}^{(d_1-d_3)}\left(k\right)& ...& C_{1M}^{(d_1-d_M)}\left(k\right)\\ C_{21}^{(d_2-d_1)}\left(k\right)& C_{22}^0\left(k\right)& C_{23}^{(d_2-d_3)}\left(k\right)& ...& C_{2M}^{(d_2-d_M)}\left(k\right)\\ C_{31}^{d_3-d_1}\left(k\right)& C_{32}^{d_3-d_2}\left(k\right)& C_{33}^0\left(k\right)& ...& C_{3M}^{(d_3-d_M)}\left(k\right)\\ ...& ...& ...& ...& ...\\ C_{M1}^{d_M-d_1}\left(k\right)& C_{M2}^{d_M-d_2}\left(k\right)& C_{M3}^{d_M-d_3}\left(k\right)& ...& C_{\mathrm{MM}}^0\left(k\right)\end{array}\right],$ 若令 ${\stackrel{\&OverBar;}{C}}^0\left(k\right)=\frac{1}{M}\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}{C}_{\mathrm{ii}}^0\left(k\right),{\stackrel{\&OverBar;}{C}}^{\mathrm{m\&epsiv;}}\left(k\right)=E\left[C_{\mathrm{ij}}^{\mathrm{m\&epsiv;}}\left(k\right)\right],{\stackrel{\&OverBar;}{C}}^{-\mathrm{m\&epsiv;}}\left(k\right)=E\left[C_{\mathrm{ij}}^{-\mathrm{m\&epsiv;}}\left(k\right)\right],m=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;,M_a,$ E()为取均值函数。虚拟均匀线阵的扩展分数低阶协方差增量

（(M_a+1)×(M_a+1)维）可以表示为 $\stackrel{\&OverBar;}{C}=\left[\begin{array}{ccccc}{\stackrel{\&OverBar;}{C}}^0\left(k\right)& {\stackrel{\&OverBar;}{C}}^{-\mathrm{\&epsiv;}}\left(k\right)& {\stackrel{\&OverBar;}{C}}^{-2\mathrm{\&epsiv;}}\left(k\right)& ...& {\stackrel{\&OverBar;}{C}}^{-M_a\mathrm{\&epsiv;}}\left(k\right)\\ {\stackrel{\&OverBar;}{C}}^{\mathrm{\&epsiv;}}\left(k\right)& {\stackrel{\&OverBar;}{C}}^0\left(k\right)& {\stackrel{\&OverBar;}{C}}^{-\mathrm{\&epsiv;}}\left(k\right)& ...& {\stackrel{\&OverBar;}{C}}^{-(M_a-1)\mathrm{\&epsiv;}}\left(k\right)\\ {\stackrel{\&OverBar;}{C}}^{2\mathrm{\&epsiv;}}\left(k\right)& {\stackrel{\&OverBar;}{C}}^{\mathrm{\&epsiv;}}\left(k\right)& {\stackrel{\&OverBar;}{C}}^0\left(k\right)& ...& {\stackrel{\&OverBar;}{C}}^{-(M_a-2)\mathrm{\&epsiv;}}\left(k\right)\\ ...& ...& ...& ...& ...\\ {\stackrel{\&OverBar;}{C}}^{M_a\mathrm{\&epsiv;}}\left(k\right)& {\stackrel{\&OverBar;}{C}}^{(M_a-1)\mathrm{\&epsiv;}}\left(k\right)& {\stackrel{\&OverBar;}{C}}^{(M_a-2)\mathrm{\&epsiv;}}\left(k\right)& ...& {\stackrel{\&OverBar;}{C}}^0\left(k\right)\end{array}\right]=\left[\begin{array}{ccccc}{R}_{11}\left(k\right)& R_{12}\left(k\right)& R_{13}\left(k\right)& ...& R_{1(M_a+1)}\left(k\right)\\ R_{21}\left(k\right)& R_{22}\left(k\right)& R_{23}\left(k\right)& ...& R_{2(M_a+1)\left(k\right)}\\ R_{31}\left(k\right)& R_{32}\left(k\right)& R_{33}\left(k\right)& ...& R_{3(M_a+1)\left(k\right)}\\ ...& ...& ...& ...& ...\\ R_{(M_a+1)1}\left(k\right)& R_{(M_a+1)2}\left(k\right)& R_{(M_a+1)3}\left(k\right)& ...& R_{(M_a+1)(M_a+1)\left(k\right)}\end{array}\right]=R$ ，则虚拟成均匀线阵扩展了原阵列的阵列流型矢量为B(θ)=[b(θ₁),b(θ₂),…,b(θ_N)],其第i个导向矢量 $b\left({\mathrm{\&theta;}}_i\right)={[1,e^{-j2\mathrm{\&pi;\&epsiv;}\sin \left({\mathrm{\&theta;}}_i\right)/\mathrm{\&lambda;}},\&CenterDot;\&CenterDot;\&CenterDot;,e^{-j2\mathrm{\&pi;}\left(M_a\right)\mathrm{\&epsiv;}\sin \left({\mathrm{\&theta;}}_i\right)/\mathrm{\&lambda;}}]}^T,i=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;,N,$ 根据非均匀线阵虚拟出M_a>M个阵元，则可测信号数为M_a-1。第一次快拍采样的扩展分数低阶协方差

（M_a×M_a维）等于其扩展分数低阶协方差增量R(1),t代表量子文化雁群搜索机制的迭代次数，初始时设t=0。

步骤二，确定雁群中所有大雁初始状态。考虑由H个大雁组成的雁群，每只大雁在N维搜索空间中运动，大雁的位置代表优化问题的解。第i只大雁的量子位置定义如下:y_i(t)=[y_i1(t),y_i2(t),…,y_iN(t)]，(i=1,2,…,H)，其中:0≤y_in(t)≤1(n=1,2,…,N)表示第i只大雁量子位置的第n维量子位，把第i只大雁的量子位置映射到定义区间，就是该大雁的位置

映射关系为 ${\stackrel{\&OverBar;}{y}}_{\mathrm{in}}\left(t\right)=y_{\mathrm{in}}\left(t\right)[h_n\left(k\right)-b_n\left(k\right)]+b_n\left(k\right),n=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;,N,$ h_n(k)和b_n(k)分别为第k次采样数据处理时第n维角度搜索区间的上限和下限值，第i只大雁对应的速度为v_i(t)=[v_i1(t),v_i2(t),…,v_iN(t)]，-0.2≤v_in(t)≤0.2(n=1,2,…,N)表示第i只大雁的第n维速度。为了使初始位置具有一定的分散性和均布性，量子位置在定义的量子域[0，1]随机初始化，速度在[-0.2，0.2]随机初始化。确定初始搜索空间确定初始搜索空 $Z\left(k\right)=\left[\begin{array}{cccc}{h}_1\left(k\right)& h_2\left(k\right)& ...& h_N\left(k\right)\\ b_1\left(k\right)& b_2\left(k\right)& ...& b_N\left(k\right)\end{array}\right]$ 为N个角度的搜索空间向量，最大迭代次数设置为

其中为取整函数，T为取整倍数。

步骤三，第i只大雁位置

的扩展分数低阶协方差极大似然方程的适应度值为

式中，正交投影矩阵 $P_{B\left[{\stackrel{\&OverBar;}{y}}_i\left(t\right)\right]}=B\left({\stackrel{\&OverBar;}{y}}_i\left(t\right)\right){\left[B^H\left({\stackrel{\&OverBar;}{y}}_i\left(t\right)\right)B\left({\stackrel{\&OverBar;}{y}}_i\left(t\right)\right)\right]}^{-1}{B}^H\left({\stackrel{\&OverBar;}{y}}_i\left(t\right)\right),$ 上标H代表共轭的转置，tr()表示矩阵的求迹运算，适应度值越大，量子位置和位置越优秀，估计的角度越准确。根据

评价大雁位置和对应量子位置的优劣，适应度值越大，位置越优。第i只大雁到现在为止所经历的最优量子位置定义为该大雁的局部最优量子位置，记作p_i(t)=[p_i1(t),p_i2(t),…,p_iN(t)]，p_in(t)为第i只大雁到第t次迭代为止所经历第n维最优量子位，n=1,2,…,N。所有大雁到现在所经历的最优量子位置记作全局最优量子位置，也就是适应度值最大的量子位置，记作g(t)=[g₁(t),g₂(t),…,g_N(t)]，g_n(t)为所有大雁到第t次迭代为止所经历第n维最优量子位，n=1,2,…,N。第n维量子形势知识用I_n=[l_n,u_n]表示，下限l_n和上限u_n根据问题所给定的量子位置定义域的边界进行初始化；L_n表示第n维量子位置的下限l_n所对应的适应度值，U_n表示第n维量子位置的上限u_n所对应的适应度值，在最大值优化问题中均可初始化为-∞。

步骤四，更新每只大雁的速度和量子位置。H只大雁速度先根据群行为进行更新，则第i只大雁的第n维速度更新为v_in(t+1)=c₁r₁[g_n(t)-y_in(t)]+c₂r₂[p_(i-1)n(t)-y_in(t)]+w_tv_in(t)，w_t从初次迭代时为0.9线形递减到迭代次数最大时为0.5，i=1,2,…,H，n=1,2,…,N，r₁和r₂都是[0，1]之间的均匀随机数，c₁和c₂加权常数，；对于v_in(t+1)，若超出边界值，将其限制在边界，即v_in(t+1)>0.2，v_in(t+1)=0.2，若v_in(t+1)<-0.2，

v_in(t+1)=-0.2。第i只大雁的新量子位置为

y_i(t+1)=[y_i1(t+1),y_i2(t+1),…,y_iN(t+1)]，i=1,2,…,H，其中

$y_{\mathrm{in}}(t+1)=\mathrm{abs}[y_{\mathrm{in}}\left(t\right)\cos (-v_{\mathrm{in}}(t+1))-\sqrt{1-{\left(y_{\mathrm{in}}\left(t\right)\right)}^2}\sin (-v_{\mathrm{in}}(t+1))],$

n=1,2,…,N，abs()为量子位取绝对值函数。

步骤五，使用量子规范知识指导更新量子位置为y_i(t+1)，i=H+1,H+2,…,2H，第i只大雁的第n维的量子旋转角为

第i只大雁的量子位置的第n维更新为

n=1,2,…,N，G(0,1)是均值为0方差为1的高斯随机数，η为加权常数，一般取η=0.06。

步骤六，将第i只大雁的量子位置y_i(t+1)影射到定义区间计算适应度值，适应度函数为 $f\left[{\stackrel{\&OverBar;}{y}}_i(t+1)\right]=\mathrm{tr}\left[P_{B\left[{\stackrel{\&OverBar;}{y}}_i(t+1)\right]}\stackrel{\&OverBar;}{R}\left(k\right)\right],i=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;,2H,$ 对于第i（i=1,2,…,H）只大雁，若

则令y_i(t+1)=y_i+N(t+1)。若p_i(t+1)=y_i(t+1)，否则，p_i(t+1)=p_i(t)。把所有大雁的局部最优量子位置按照适应度值由大到小进行排序，把最优（适应度值最大）的局部最优量子位置设置成全局最优量子位置为g(t+1)=[g₁(t+1),g₂(t+1),…,g_N(t+1)]。

从所有大雁局部最优位置中选取最优秀的前εH个局部最优量子位置去更新量子规范知识，设第j只大雁影响第i维量子规范知识的下限，第m只大雁影响第i维量子规范知识的上限，通过下面的量子规范知识更新函数对第n维量子规范知识进行更新：

其中，

分别表示第t代第n维量子变量的下限和上限所对应适应度；

分别为第j只大雁和第m只大雁的局部最优量子位置p_j和p_k映射到定义域区间的潜在解向量。

步骤七，判断是否达到最大迭代次数，若是，记录最优位置，执行步骤八；否则，t=t+1，返回步骤四。

步骤八，快拍采样新数据x(k+1)=[x₁(k+1),x₂(k+1),…,x_M(k+1)]^T的无穷范数归一化信号为 $\stackrel{\&OverBar;}{x}(k+1)={[{\stackrel{\&OverBar;}{x}}_1(k+1),{\stackrel{\&OverBar;}{x}}_2(k+1),\&CenterDot;\&CenterDot;\&CenterDot;,{\stackrel{\&OverBar;}{x}}_M(k+1)]}^T=x(k+1)/\underset{1\&le;j\&le;M}{\max }\left\{\right|x_j(k+1)\left|\right\},$ 则 $C_{\mathrm{ij}}(k+1)={\stackrel{\&OverBar;}{x}}_i(k+1){\left|{\stackrel{\&OverBar;}{x}}_i(k+1)\right|}^{p-1}{\stackrel{\&OverBar;}{x}}_j^{*}(k+1){\left|{\stackrel{\&OverBar;}{x}}_j(k+1)\right|}^{p-1}$ 是新增加的第k+1个采样数据的分数低阶协方差增量矩阵。依照步骤一的方法，根据分数低阶协方差增量矩阵构造出虚拟均匀线阵的扩展分数低阶协方差增量

维），且可以表示为 $\stackrel{\&OverBar;}{C}=\left[\begin{array}{ccccc}{\stackrel{\&OverBar;}{C}}^0\left(k\right)& {\stackrel{\&OverBar;}{C}}^{-\mathrm{\&epsiv;}}\left(k\right)& {\stackrel{\&OverBar;}{C}}^{-2\mathrm{\&epsiv;}}\left(k\right)& ...& {\stackrel{\&OverBar;}{C}}^{-M_a\mathrm{\&epsiv;}}\left(k\right)\\ {\stackrel{\&OverBar;}{C}}^{\mathrm{\&epsiv;}}\left(k\right)& {\stackrel{\&OverBar;}{C}}^0\left(k\right)& {\stackrel{\&OverBar;}{C}}^{-\mathrm{\&epsiv;}}\left(k\right)& ...& {\stackrel{\&OverBar;}{C}}^{-(M_a-1)\mathrm{\&epsiv;}}\left(k\right)\\ {\stackrel{\&OverBar;}{C}}^{2\mathrm{\&epsiv;}}\left(k\right)& {\stackrel{\&OverBar;}{C}}^{\mathrm{\&epsiv;}}\left(k\right)& {\stackrel{\&OverBar;}{C}}^0\left(k\right)& ...& {\stackrel{\&OverBar;}{C}}^{-(M_a-2)\mathrm{\&epsiv;}}\left(k\right)\\ ...& ...& ...& ...& ...\\ {\stackrel{\&OverBar;}{C}}^{M_a\mathrm{\&epsiv;}}\left(k\right)& {\stackrel{\&OverBar;}{C}}^{(M_a-1)\mathrm{\&epsiv;}}\left(k\right)& {\stackrel{\&OverBar;}{C}}^{(M_a-2)\mathrm{\&epsiv;}}\left(k\right)& ...& {\stackrel{\&OverBar;}{C}}^0\left(k\right)\end{array}\right]=\left[\begin{array}{ccccc}{R}_{11}\left(k\right)& R_{12}\left(k\right)& R_{13}\left(k\right)& ...& R_{1(M_a+1)}\left(k\right)\\ R_{21}\left(k\right)& R_{22}\left(k\right)& R_{23}\left(k\right)& ...& R_{2(M_a+1)\left(k\right)}\\ R_{31}\left(k\right)& R_{32}\left(k\right)& R_{33}\left(k\right)& ...& R_{3(M_a+1)\left(k\right)}\\ ...& ...& ...& ...& ...\\ R_{(M_a+1)1}\left(k\right)& R_{(M_a+1)2}\left(k\right)& R_{(M_a+1)3}\left(k\right)& ...& R_{(M_a+1)(M_a+1)\left(k\right)}\end{array}\right]=R$ 。通过下面的方式更新当前采样的扩展分数低阶协方差：

其中R_ij(k+1)新增加的第k+1个采样数据的扩展分数低阶协方差的第i行第j列增量。更新搜索空间为 $Z(k+1)=\left[\begin{array}{cccc}{h}_1(k+1)& h_2(k+1)& ...& h_N(k+1)\\ b_1(k+1)& b_2(k+1)& ...& b_N(k+1)\end{array}\right],$ $b_n(k+1)={\stackrel{\&OverBar;}{\mathrm{\&theta;}}}_{\mathrm{zn}}\left(k\right)-{\mathrm{\&beta;}}^k|b_n\left(k\right)-{\stackrel{\&OverBar;}{\mathrm{\&theta;}}}_n\left(k\right)|-r,h_n(k+1)={\stackrel{\&OverBar;}{\mathrm{\&theta;}}}_{\mathrm{zn}}\left(k\right)+{\mathrm{\&beta;}}^k|h_n\left(k\right)-{\stackrel{\&OverBar;}{\mathrm{\&theta;}}}_n\left(k\right)|+r,$ n=1,2,…,N，β为收敛因子，决定了搜索空间的收敛速度；常数r为搜索空间在锁定状态下的搜索半径；

为第n个方向在第k次采样时的估计值；

为第n个方向在第k次采样搜索空间的中心值，其更新公式为

其中δ为遗传因子。

步骤九，如果达到最大快拍采样数，执行步骤十；否则，设k=k+1，迭代次数t=0，返回步骤二继续估计动态目标下一个时刻的方向。

步骤十，得到所有快拍采样下的一系列全局最优量子位置，其映射得到的全局最优位置就是跟踪到的动态目标方向值，输出动态跟踪结果。

本发明考虑到冲击噪声环境下完成动态跟踪的过程中能够同时考虑收敛精度，收敛速度和阵列扩展能力，使用量子文化雁群搜索机制去求解扩展分数低阶协方差极大似然方程获得最优估计值。所设计的动态跟踪系统还可以根据时间要求和性能要求自动确定迭代次数，从而使所设计的量子文化雁群动态跟踪方法满足更高性能要求。

图2是量子文化雁群搜索机制的结构示意图，在实际动态跟踪中还会遇到一种情况，就是试图放置阵元的位置已被占用不能再放置阵元，这个时候，可以根据冗余的理论，寻找一些次优的位置，同样获得阵列扩展和优良的解相干性能，冗余量大有助于解相干。现在举例说明，假设在距线阵原点位置3.5λ和4.5λ处已有其它用处，不可放置阵元，还可以选择其它阵列扩展形式的阵元放置位置。同样是5个阵元，阵元放置位置为0.5λ[0 1 4 5 8]，但该方法放置的阵元扩展出虚拟阵元数为9。所提的量子文化雁群动态跟踪方法，扩展分数低阶协方差更新公式中的μ=0.95，方法中收敛因子取β=0.995，遗传因子δ=0.8，收敛半径取r=3，ε=0.2，c₁=c₂=0.78。在实验仿真过程中使用的用于比较的动态跟踪方法是基于粒子群算法和分数低阶矩极大似然算法（PSO-FLOM）的动态DOA估计(赵大勇等在《山东大学学报》(2010,Vol.40,No.1,pp.133-138)上发表的“冲击噪声背景下的动态DOA跟踪”)，实验中使用阵元间隔为0.5倍波长的等距均匀线阵，阵元数M=5。为了考察两种动态跟踪方法的收敛速度和性能，搜索区间设为[-90°,90°]，量子文化雁群动态跟踪方法和粒子群算法的种群规模（所含个体数）都设为30，最大迭代次数设置方法相同，都取该次快拍第一维搜索边界区域差的整倍数，倍数值为T=4。假设冲击噪声环境下的广义信噪比为15dB，仿真中所有信源设置为等功率信源，若信源的初值为θ_i(0)，最大采样数为K，当前采样数为k，则动态目标的角度的运动轨迹为 ${\mathrm{\&theta;}}_i\left(k\right)={\mathrm{\&theta;}}_i\left(0\right)+5\sin \left(2\mathrm{\&pi;}\frac{k}{500}\right),i=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;,N,$ N为信源个数。

从图3、图4、图5和图6可以看出，所提量子文化雁群多移动目标跟踪方法在冲击噪声背景下对独立信源的跟踪精度要好于基于粒子群算法和分数低阶矩极大似然算法的动态DOA估计方法，有更好的收敛精度。

从图7和图8可以看出，所提量子文化雁群动态跟踪方法在信源数大于等于阵元数的情况下可对多个目标可进行较精确的动态跟踪，基于粒子群算法和分数低阶矩极大似然算法的动态DOA方法失效，因为基于粒子群算法和分数低阶矩极大似然算法的动态DOA跟踪仅适用阵元数大于信源数的情况。

从图9和图10可以看出，所提的量子文化雁群动态跟踪方法在特征指数小于1的强冲击噪声情况下可对2个目标可进行较精确的动态跟踪，而基于粒子群算法和分数低阶矩极大似然算法的动态DOA方法失效，因为基于粒子群算法和分数低阶矩极大似然算法的动态DOA跟踪仅适用特征指数大于1的情况。

以图11的鲁棒动态跟踪装置为例，其它情况可以依此类推，在动态DOA估计系统存在M个天线，使用特殊阵列结构放置阵元，例如可选择最小冗余线阵、最大连续延迟线阵和最小间隙阵列等特殊的非均匀线阵，对N个移动目标所发射的信号对快拍采样信号进行无穷范数归一化后求其分数低阶协方差，根据特殊阵列结构形式把分数低阶协方差构造出扩展分数低阶协方差。使用量子文化雁群动态目标跟踪方法搜索波达方向的估计值。以当前采样的波达方向估计值为中心确定下一个采样的搜索范围，对接收到的当前快拍采样数据的无穷范数归一化数据的扩展分数低阶协方差增量和前一快拍采样信号扩展分数低阶协方差用于更新扩展分数低阶信号协方差，使用量子文化雁群优化搜索机制去优化扩展分数低阶协方差极大似然方程，使用量子群行为和知识策略更新大雁的速度和量子位置，根据局部最优量子位置更新量子规范知识，迭代求得大雁的最优量子位置，根据映射态求得多移动目标方向值，每次快拍采样后都循环迭代上述波达方向估计过程，完成波达方向的动态跟踪过程。

Claims (2)

1.一种冲击噪声环境下的多移动目标跟踪方法，其特征是：

（1）根据天线个数和测向要求确定天线位置，阵列远场在θ_k(k=1,2,…,N)处有N个窄带点源以平面波入射，入射波长为λ，入射到由M个各向同性的阵元构成的非均匀特殊线阵上，阵元的放置满足如下的条件：第i个阵元相对于第一个阵元的间距为d₁=0<d₂<…<d_M，其值都是半波长的整数倍，并满足Ω={d_i-d_j|i,j=1,2,…,M}是完整的，即Ω={-d_M,…,0,…,d_M}，阵列接收的快拍数据可表示为x(k)=A(θ)s(k)+n(k)，式中x(k)=[x₁(k),x₂(k),…,x_M(k)]^T为M×1维阵列快拍数据矢量，s(t)=[s₁(t),s₂(t),…,s_N(t)]^T为N×1维非高斯信号矢量，θ=(θ₁,θ₂,…,θ_N)为信源方位矢量，A(θ)=[a(θ₁)a(θ₂)…a(θ_N)]是信号导向矢量矩阵，其中第i个导向矢量 $a_i={[1,e^{-j2\mathrm{\&pi;}{d}_1\sin \left({\mathrm{\&theta;}}_i\right)/\mathrm{\&lambda;}},\&CenterDot;\&CenterDot;\&CenterDot;,e^{-j2\mathrm{\&pi;}{d}_M\sin \left({\mathrm{\&theta;}}_i\right)/\mathrm{\&lambda;}}]}^T,$ (i=1,2,…,N)；第k次采样数据的无穷范数归一化信号 $\stackrel{\&OverBar;}{x}\left(k\right)={[{\stackrel{\&OverBar;}{x}}_1\left(k\right),{\stackrel{\&OverBar;}{x}}_2\left(k\right),\&CenterDot;\&CenterDot;\&CenterDot;,{\stackrel{\&OverBar;}{x}}_M\left(k\right)]}^T=x\left(k\right)/\underset{1\&le;j\&le;M}{\max }\left\{\right|x_j\left(k\right)\left|\right\},$ 其分数低阶协方差为C(k)，其中第i维第j列元素为 $C_{\mathrm{ij}}\left(k\right)={\stackrel{\&OverBar;}{x}}_i\left(k\right){\left|{\stackrel{\&OverBar;}{x}}_i\left(k\right)\right|}^{p-1}{\stackrel{\&OverBar;}{x}}_j^{*}\left(k\right){\left|{\stackrel{\&OverBar;}{x}}_j\left(k\right)\right|}^{p-1},$ *为共轭算子，p为分数低阶协方差参数，把非均匀线阵可虚拟成更多个阵元的均匀线阵，阵列的最大相关延迟为M_a，则虚拟出的均匀直线阵阵元个数为M+1～M_a+1，设则阵元坐标矢量为d=[d₁,d₂,…,d_M]=ε[n₁,…,n_M]，n₁,…,n_M都为整数，则接收阵列的分数低阶协方差C(k)（M×M维）可以写作 $C\left(k\right)=\left[\begin{array}{ccccc}{C}_{11}^0\left(k\right)& C_{12}^{(d_1-d_2)}\left(k\right)& C_{13}^{(d_1-d_3)}\left(k\right)& ...& C_{1M}^{(d_1-d_M)}\left(k\right)\\ C_{21}^{(d_2-d_1)}\left(k\right)& C_{22}^0\left(k\right)& C_{23}^{(d_2-d_3)}\left(k\right)& ...& C_{2M}^{(d_2-d_M)}\left(k\right)\\ C_{31}^{d_3-d_1}\left(k\right)& C_{32}^{d_3-d_2}\left(k\right)& C_{33}^0\left(k\right)& ...& C_{3M}^{(d_3-d_M)}\left(k\right)\\ ...& ...& ...& ...& ...\\ C_{M1}^{d_M-d_1}\left(k\right)& C_{M2}^{d_M-d_2}\left(k\right)& C_{M3}^{d_M-d_3}\left(k\right)& ...& C_{\mathrm{MM}}^0\left(k\right)\end{array}\right],$ 若令 ${\stackrel{\&OverBar;}{C}}^0\left(k\right)=\frac{1}{M}\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}{C}_{\mathrm{ii}}^0\left(k\right),{\stackrel{\&OverBar;}{C}}^{\mathrm{m\&epsiv;}}\left(k\right)=E\left[C_{\mathrm{ij}}^{\mathrm{m\&epsiv;}}\left(k\right)\right],{\stackrel{\&OverBar;}{C}}^{-\mathrm{m\&epsiv;}}\left(k\right)=E\left[C_{\mathrm{ij}}^{-\mathrm{m\&epsiv;}}\left(k\right)\right],m=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;,M_a,$ E()为取均值函数，虚拟均匀线阵的扩展分数低阶协方差增量（(M_a+1)×(M_a+1)维）可以表示为 $\stackrel{\&OverBar;}{C}=\left[\begin{array}{ccccc}{\stackrel{\&OverBar;}{C}}^0\left(k\right)& {\stackrel{\&OverBar;}{C}}^{-\mathrm{\&epsiv;}}\left(k\right)& {\stackrel{\&OverBar;}{C}}^{-2\mathrm{\&epsiv;}}\left(k\right)& ...& {\stackrel{\&OverBar;}{C}}^{-M_a\mathrm{\&epsiv;}}\left(k\right)\\ {\stackrel{\&OverBar;}{C}}^{\mathrm{\&epsiv;}}\left(k\right)& {\stackrel{\&OverBar;}{C}}^0\left(k\right)& {\stackrel{\&OverBar;}{C}}^{-\mathrm{\&epsiv;}}\left(k\right)& ...& {\stackrel{\&OverBar;}{C}}^{-(M_a-1)\mathrm{\&epsiv;}}\left(k\right)\\ {\stackrel{\&OverBar;}{C}}^{2\mathrm{\&epsiv;}}\left(k\right)& {\stackrel{\&OverBar;}{C}}^{\mathrm{\&epsiv;}}\left(k\right)& {\stackrel{\&OverBar;}{C}}^0\left(k\right)& ...& {\stackrel{\&OverBar;}{C}}^{-(M_a-2)\mathrm{\&epsiv;}}\left(k\right)\\ ...& ...& ...& ...& ...\\ {\stackrel{\&OverBar;}{C}}^{M_a\mathrm{\&epsiv;}}\left(k\right)& {\stackrel{\&OverBar;}{C}}^{(M_a-1)\mathrm{\&epsiv;}}\left(k\right)& {\stackrel{\&OverBar;}{C}}^{(M_a-2)\mathrm{\&epsiv;}}\left(k\right)& ...& {\stackrel{\&OverBar;}{C}}^0\left(k\right)\end{array}\right]=\left[\begin{array}{ccccc}{R}_{11}\left(k\right)& R_{12}\left(k\right)& R_{13}\left(k\right)& ...& R_{1(M_a+1)}\left(k\right)\\ R_{21}\left(k\right)& R_{22}\left(k\right)& R_{23}\left(k\right)& ...& R_{2(M_a+1)\left(k\right)}\\ R_{31}\left(k\right)& R_{32}\left(k\right)& R_{33}\left(k\right)& ...& R_{3(M_a+1)\left(k\right)}\\ ...& ...& ...& ...& ...\\ R_{(M_a+1)1}\left(k\right)& R_{(M_a+1)2}\left(k\right)& R_{(M_a+1)3}\left(k\right)& ...& R_{(M_a+1)(M_a+1)\left(k\right)}\end{array}\right]=R$ ，则虚拟成均匀线阵扩展了原阵列的阵列流型矢量为B(θ)=[b(θ₁),b(θ₂),…,b(θ_N)],其第i个导向矢量 $b\left({\mathrm{\&theta;}}_i\right)={[1,e^{-j2\mathrm{\&pi;\&epsiv;}\sin \left({\mathrm{\&theta;}}_i\right)/\mathrm{\&lambda;}},\&CenterDot;\&CenterDot;\&CenterDot;,e^{-j2\mathrm{\&pi;}\left(M_a\right)\mathrm{\&epsiv;}\sin \left({\mathrm{\&theta;}}_i\right)/\mathrm{\&lambda;}}]}^T,i=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;,N,$ 根据非均匀线阵虚拟出M_a>M个阵元，则可测信号数为M_a-1，第一次快拍采样的扩展分数低阶协方差（M_a×M_a维）等于其扩展分数低阶协方差增量R(1),t代表量子文化雁群搜索机制的迭代次数，初始时设t=0；

（2）确定雁群中所有大雁初始状态，由H个大雁组成的雁群，每只大雁在N维搜索空间中运动，第i只大雁的量子位置定义如下:

y_i(t)=[y_i1(t),y_i2(t),…,y_iN(t)]表示第i只大雁量子位置的第n维量子位，

(i=1,2,…,H)，0≤y_in(t)≤1(n=1,2,…,N)，把第i只大雁的量子位置映射到定义区间，就是该大雁的位置

映射关系为 ${\stackrel{\&OverBar;}{y}}_{\mathrm{in}}\left(t\right)=y_{\mathrm{in}}\left(t\right)[h_n\left(k\right)-b_n\left(k\right)]+b_n\left(k\right),n=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;,N,$ h_n(k)和b_n(k)分别为第k次采样数据处理时第n维角度搜索区间的上限和下限值，第i只大雁对应的速度为v_i(t)=[v_i1(t),v_i2(t),…,v_iN(t)]，-0.2≤v_in(t)≤0.2(n=1,2,…,N)表示第i只大雁的第n维速度，量子位置在定义的量子域[0，1]随机初始化，速度在[-0.2，0.2]随机初始化，确定初始搜索空间确定初始搜索空 $Z\left(k\right)=\left[\begin{array}{cccc}{h}_1\left(k\right)& h_2\left(k\right)& ...& h_N\left(k\right)\\ b_1\left(k\right)& b_2\left(k\right)& ...& b_N\left(k\right)\end{array}\right]$ 为N个角度的搜索空间向量，最大迭代次数设置为

其中为取整函数，T为取整倍数；

（3）第i只大雁位置

的扩展分数低阶协方差极大似然方程的适应度值为

正交投影矩阵 $P_{B\left[{\stackrel{\&OverBar;}{y}}_i\left(t\right)\right]}=B\left({\stackrel{\&OverBar;}{y}}_i\left(t\right)\right){\left[B^H\left({\stackrel{\&OverBar;}{y}}_i\left(t\right)\right)B\left({\stackrel{\&OverBar;}{y}}_i\left(t\right)\right)\right]}^{-1}{B}^H\left({\stackrel{\&OverBar;}{y}}_i\left(t\right)\right),$ 上标H代表共轭的转置，tr()表示矩阵的求迹运算，根据

评价大雁位置和对应量子位置的优劣，第i只大雁到现在为止所经历的最优量子位置定义为该大雁的局部最优量子位置，记作p_i(t)=[p_i1(t),p_i2(t),…,p_iN(t)]，p_in(t)为第i只大雁到第t次迭代为止所经历第n维最优量子位，n=1,2,…,N，所有大雁到现在所经历的最优量子位置记作全局最优量子位置，也就是适应度值最大的量子位置，记作g(t)=[g₁(t),g₂(t),…,g_N(t)]，g_n(t)为所有大雁到第t次迭代为止所经历第n维最优量子位，n=1,2,…,N，第n维量子形势知识用I_n=[l_n,u_n]表示，下限l_n和上限u_n根据问题所给定的量子位置定义域的边界进行初始化；L_n表示第n维量子位置的下限l_n所对应的适应度值，U_n表示第n维量子位置的上限u_n所对应的适应度值，在最大值优化问题中均可初始化为-∞；

（4）更新每只大雁的速度和量子位置，H只大雁速度先根据群行为进行更新，则第i只大雁的第n维速度更新为v_in(t+1)=c₁r₁[g_n(t)-y_in(t)]+c₂r₂[p(_i-1)_n(t)-y_in(t)]+w_tv_in(t)，w_t从初次迭代时为0.9线形递减到迭代次数最大时为0.5，i=1,2,…,H，n=1,2,…,N，r₁和r₂都是[0，1]之间的均匀随机数，c₁和c₂加权常数，对于v_in(t+1)，若超出边界值，将其限制在边界，即v_in(t+1)>0.2，v_in(t+1)=0.2，若v_in(t+1)<-0.2，

v_in(t+1)=-0.2，第i只大雁的新量子位置为

y_i(t+1)=[y_i1(t+1),y_i2(t+1),…,y_iN(t+1)]，i=1,2,…,H，

$y_{\mathrm{in}}(t+1)=\mathrm{abs}[y_{\mathrm{in}}\left(t\right)\cos (-v_{\mathrm{in}}(t+1))-\sqrt{1-{\left(y_{\mathrm{in}}\left(t\right)\right)}^2}\sin (-v_{\mathrm{in}}(t+1))],$

n=1,2,…,N，abs()为量子位取绝对值函数；

（5）使用量子规范知识指导更新量子位置为y_i(t+1)，i=H+1,H+2,…,2H，第i只大雁的第n维的量子旋转角为

第i只大雁的量子位置的第n维更新为

n=1,2,…,N，G(0,1)是均值为0方差为1的高斯随机数，η为加权常数；

（6）将第i只大雁的量子位置y_i(t+1)影射到定义区间

计算适应度值，适应度函数为 $f\left[{\stackrel{\&OverBar;}{y}}_i(t+1)\right]=\mathrm{tr}\left[P_{B\left[{\stackrel{\&OverBar;}{y}}_i(t+1)\right]}\stackrel{\&OverBar;}{R}\left(k\right)\right],i=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;,2H,$ 对于第i只大雁，若

则令y_i(t+1)=y_i+N(t+1)，若

p_i(t+1)=y_i(t+1)，否则，p_i(t+1)=p_i(t)，把所有大雁的局部最优量子位置按照适应度值由大到小进行排序，把最优即适应度值最大的局部最优量子位置设置成全局最优量子位置为g(t+1)=[g₁(t+1),g₂(t+1),…,g_N(t+1)]；

从所有大雁局部最优位置中选取最优的前εH个局部最优量子位置去更新量子规范知识，设第j只大雁影响第i维量子规范知识的下限，第m只大雁影响第i维量子规范知识的上限，通过下面的量子规范知识更新函数对第n维量子规范知识进行更新：

其中，

分别表示第t代第n维量子变量的下限

和上限

所对应适应度；

分别为第j只大雁和第m只大雁的局部最优量子位置p_j和p_k映射到定义域区间的潜在解向量；

（7）判断是否达到最大迭代次数，若是，记录最优位置，执行步骤（8）；否则t=t+1，返回步骤（4）进行迭代；

（8）快拍采样新数据x(k+1)=[x₁(k+1),x₂(k+1),…,x_M(k+1)]^T的无穷范数归一化信号为 $\stackrel{\&OverBar;}{x}(k+1)={[{\stackrel{\&OverBar;}{x}}_1(k+1),{\stackrel{\&OverBar;}{x}}_2(k+1),\&CenterDot;\&CenterDot;\&CenterDot;,{\stackrel{\&OverBar;}{x}}_M(k+1)]}^T=x(k+1)/\underset{1\&le;j\&le;M}{\max }\left\{\right|x_j(k+1)\left|\right\},$ 则 $C_{\mathrm{ij}}(k+1)={\stackrel{\&OverBar;}{x}}_i(k+1){\left|{\stackrel{\&OverBar;}{x}}_i(k+1)\right|}^{p-1}{\stackrel{\&OverBar;}{x}}_j^{*}(k+1){\left|{\stackrel{\&OverBar;}{x}}_j(k+1)\right|}^{p-1}$ 是新增加的第k+1个采样数据的分数低阶协方差增量矩阵，依照步骤（1）的方法，根据分数低阶协方差增量矩阵构造出虚拟均匀线阵的扩展分数低阶协方差增量

（(M_a+1)×(M_a+1)维），且可以表示为 $\stackrel{\&OverBar;}{C}=\left[\begin{array}{ccccc}{\stackrel{\&OverBar;}{C}}^0\left(k\right)& {\stackrel{\&OverBar;}{C}}^{-\mathrm{\&epsiv;}}\left(k\right)& {\stackrel{\&OverBar;}{C}}^{-2\mathrm{\&epsiv;}}\left(k\right)& ...& {\stackrel{\&OverBar;}{C}}^{-M_a\mathrm{\&epsiv;}}\left(k\right)\\ {\stackrel{\&OverBar;}{C}}^{\mathrm{\&epsiv;}}\left(k\right)& {\stackrel{\&OverBar;}{C}}^0\left(k\right)& {\stackrel{\&OverBar;}{C}}^{-\mathrm{\&epsiv;}}\left(k\right)& ...& {\stackrel{\&OverBar;}{C}}^{-(M_a-1)\mathrm{\&epsiv;}}\left(k\right)\\ {\stackrel{\&OverBar;}{C}}^{2\mathrm{\&epsiv;}}\left(k\right)& {\stackrel{\&OverBar;}{C}}^{\mathrm{\&epsiv;}}\left(k\right)& {\stackrel{\&OverBar;}{C}}^0\left(k\right)& ...& {\stackrel{\&OverBar;}{C}}^{-(M_a-2)\mathrm{\&epsiv;}}\left(k\right)\\ ...& ...& ...& ...& ...\\ {\stackrel{\&OverBar;}{C}}^{M_a\mathrm{\&epsiv;}}\left(k\right)& {\stackrel{\&OverBar;}{C}}^{(M_a-1)\mathrm{\&epsiv;}}\left(k\right)& {\stackrel{\&OverBar;}{C}}^{(M_a-2)\mathrm{\&epsiv;}}\left(k\right)& ...& {\stackrel{\&OverBar;}{C}}^0\left(k\right)\end{array}\right]=\left[\begin{array}{ccccc}{R}_{11}\left(k\right)& R_{12}\left(k\right)& R_{13}\left(k\right)& ...& R_{1(M_a+1)}\left(k\right)\\ R_{21}\left(k\right)& R_{22}\left(k\right)& R_{23}\left(k\right)& ...& R_{2(M_a+1)\left(k\right)}\\ R_{31}\left(k\right)& R_{32}\left(k\right)& R_{33}\left(k\right)& ...& R_{3(M_a+1)\left(k\right)}\\ ...& ...& ...& ...& ...\\ R_{(M_a+1)1}\left(k\right)& R_{(M_a+1)2}\left(k\right)& R_{(M_a+1)3}\left(k\right)& ...& R_{(M_a+1)(M_a+1)\left(k\right)}\end{array}\right]=R$ ，通过下面的方式更新当前采样的扩展分数低阶协方差：

其中R_ij(k+1)新增加的第k+1个采样数据的扩展分数低阶协方差的第i行第j列增量，更新搜索空间为 $Z(k+1)=\left[\begin{array}{cccc}{h}_1(k+1)& h_2(k+1)& ...& h_N(k+1)\\ b_1(k+1)& b_2(k+1)& ...& b_N(k+1)\end{array}\right],$ $b_n(k+1)={\stackrel{\&OverBar;}{\mathrm{\&theta;}}}_{\mathrm{zn}}\left(k\right)-{\mathrm{\&beta;}}^k|b_n\left(k\right)-{\stackrel{\&OverBar;}{\mathrm{\&theta;}}}_n\left(k\right)|-r,h_n(k+1)={\stackrel{\&OverBar;}{\mathrm{\&theta;}}}_{\mathrm{zn}}\left(k\right)+{\mathrm{\&beta;}}^k|h_n\left(k\right)-{\stackrel{\&OverBar;}{\mathrm{\&theta;}}}_n\left(k\right)|+r,$ n=1,2,…,N，β为收敛因子；常数r为搜索空间在锁定状态下的搜索半径；为第n个方向在第k次采样时的估计值；

为第n个方向在第k次采样搜索空间的中心值，其更新公式为

其中δ为遗传因子；

（9）如果达到最大快拍采样数，执行步骤（10）；否则设k=k+1，迭代次数t=0，返回步骤（2）继续估计动态目标下一个时刻的方向；

（10）得到所有快拍采样下的一系列全局最优量子位置，其映射得到的全局最优位置就是跟踪到的动态目标方向值，输出动态跟踪结果。

2.根据权利要求1所述的一种冲击噪声环境下的多移动目标跟踪方法，其特征是：加权常数η的取值为η=0.06。

Images