一种基于连续量子鸽群的360度全方位动态测向方法与流程

本发明属于阵列信号处理技术领域，具体涉及一种基于连续量子鸽群的360度全方位动态测向方法。

背景技术：

测向又称为空间谱估计或者波达方向(directionofarrival,doa)估计，是阵列信号处理的一个重要研究领域，在卫星和移动通信系统、雷达、地震学等方面有着广泛的应用。传统的doa估计往往针对的是固定信源，然而在实际情况中，信源的角度是随时变化的，因此，对入射角度随时间变化的信号源进行动态doa跟踪是doa估计理论应用中的一个重要课题，如何在非相干信源情况下以及相干信源情况下实现360度全方位动态doa估计更是工程应用中遇到的难点之一。

经对已有文献的检索发现，刁鸣等在《系统与电子技术》(2009,vol.29,no.12,pp.2046-2049)

上发表的“一种新的基于粒子群算法的doa跟踪方法”中利用均匀线阵进行多移动目标的动态doa估计，然而均匀线阵只能提供180度的方位角测向。为了实现360度全方位的动态测向，本发明选择均匀圆阵的阵列结构，均匀圆阵以其特殊的阵列结构和良好的测向性能，在测向应用中备受青睐，但是由于圆阵的阵列流形矩阵是非vandermonde矩阵，在数学处理上比较困难，所以要对均匀圆形阵列的输出信号进行相位模式激励，使其成为模式空间内的虚拟线阵，其阵列流形矩阵才具有vandermonde结构，一些只适用于等距均匀线阵的算法才可以运用于等距均匀圆阵。现有的圆阵模式空间处理方法没有应用到动态目标的doa跟踪研究方向上。

此外，已有文献还表明测向方法大多使用基于子空间的测向方法，这些方法不仅在低信噪比时性能较差，而且在跟踪相干信源的情况下还需要解相干处理，步骤繁琐。本发明设计的是圆阵模式空间下基于极大似然的动态测向方法，不仅在信噪比低，快拍数小以及相干信源的情况下具有优越的doa估计性能，还能直接处理阵列数据协方差矩阵，避免了子空间类测向方法不断进行协方差矩阵分解的问题，实现了非相干源和相干源360度的全方位动态测向。但是，极大似然方法的一个主要缺点是它涉及一个多维非线性优化问题，计算量巨大，耗时长。

因此，本发明设计了一种圆阵模式空间下基于连续量子鸽群的360度全方位动态测向方法。该方法利用连续量子鸽群机制在搜索区间内搜索极大似然方程估计的最优角度，通过逐步缩小搜索区间减少运算量，同时依据量子编码和模拟量子演化方程设计的鸽群计算方法还可以加快算法的收敛速度，增强全局搜索能力，提高跟踪精度。

技术实现要素：

本发明的目的在于提供可减少搜索时间，提高搜索速度和跟踪速度，具有良好的实时性的一种基于连续量子鸽群的360度全方位动态测向方法。

本发明的目的通过如下技术方案来实现：

一种基于连续量子鸽群的360度全方位动态测向方法，包括如下步骤：

步骤一：获取阵列接收的快拍数据，定义变换矩阵t，计算阵列信号经模式变换后的协方差矩阵。

假设n个窄带远场信源以方位角θn，俯仰角方向入射到一个由m个阵元构成的半径为r的均匀圆阵上，n＝1,2,…,n，入射波长为λ，那么阵列接收的第k次快拍数据的数学模型为式中，k为最大快拍数，x(k)＝[x1(k),x2(k),…,xm(k)]^t为m×1维的阵列接收的快拍数据矢量，为阵列m×n维的流形矩阵，θ＝[θ1,θ2,…,θn]和为信源的方向矢量，是流形矩阵的第n个导向矢量，其中s(k)为n×1维的信号矢量，n(k)为m×1维的噪声矢量，阵列噪声假定为空间和时间均独立的高斯白噪声。这里只讨论所有信源都与阵列共面的情形，即所有的俯仰角均为90°。

均匀圆阵可激发的最大模式为向下取整函数。变换矩阵t定义为式中j＝diag{j^-lj-l(-β),…,j^ljl(-β)}，其中jl(·)为l阶第一类贝塞尔函数，l＝-l,…,0,…,l。式中f^h＝[w-l,w-l+1,…,wl]，其中用t左乘x(k)可得y(k)＝tx(k)＝bs(k)+tn(k)，其中阵列流形其具有vandermonde矩阵结构，此时均匀圆阵已经转换为虚拟线阵。第1次快拍采样数据经过模式空间变换后的协方差矩阵为

步骤二：初始化搜索区间。

在第k次快拍，n个角度的搜索区间定义为其中un(k)和gn(k)分别为第k次快拍第n维角度搜索区间的上限和下限，其初始值分别取搜索区间定义域的上限和下限，n＝1,2,…,n。连续量子鸽群的最大迭代次数取第k次快拍所有维搜索区间上限与下限之差整倍数的最大值，即其中取整数，为向下取整函数。

步骤三：初始化鸽群每只鸽子的量子位置和速度，计算适应度值，更新全局最优量子位置。

鸽群中共有p只鸽子，每只鸽子的量子位置在量子域[0,1]内随机初始化，速度在[-vmax,vmax]内随机初始化。第i只鸽子量子位置定义为其中t为当前迭代次数，为第i只鸽子的第n维量子位置，且第i只鸽子的位置为量子位置的相应映射态，映射关系为第i只鸽子的速度为且

评估第i只鸽子位置状态的适应度函数为其中为的映射矩阵，trace()为求矩阵的迹。利用此适应度函数计算每只鸽子位置的适应度值，选择适应度值最大的解的量子位置作为全局最优量子位置，记作其中为到第t次迭代为止鸽群最优量子位置的第n维，n＝1,2,…,n，其相应映射态为全局最优鸽子位置

步骤四：根据地图和指南针算子更新每只鸽子的速度和量子位置及其相应映射态，计算适应度值，更新全局最优量子位置。

第i只鸽子速度的更新方程为其中i＝1,2,…,p，n＝1,2,…,n，r是地图和指南针因数，取值范围为(0,1)，是(0,1)之间的随机数。如果速度超出边界值，则将其限制在边界，即若若

第i只鸽子量子位置更新方程为同时按照映射方程获得量子位置的相应映射态。

利用适应度函数计算第i只鸽子位置的适应度值，更新全局最优量子位置，即对于第i只鸽子，若否则，

步骤五：判断地图和指南针算子操作的循环次数是否达到所要求的迭代次数：若未达到，令t＝t+1，返回步骤四继续循环；否则，令t＝t+1，停止地图和指南针算子操作，进入地标算子操作。

步骤六：进入地标算子操作，选择适应度值较优的前一半作为当前鸽群并计算当前鸽群的量子位置中心。

按照适应度值从大到小的规则对鸽子进行排序，当前代鸽子数量为p^t+1＝p^t/2，选择适应度值较优的前一半作为当前鸽群，计算当前鸽群的量子位置中心的第n维，即并将其作为参考方向，i＝1,2,…,p，n＝1,2,…,n。

步骤七：更新每只鸽子的量子位置及其相应映射态，计算适应度值，更新全局最优量子位置。

按照方程更新每只鸽子的量子位置，其中i＝1,2,…,p，n＝1,2,…,n，是(0,1)之间的随机数。按照映射方程获得量子位置的相应映射态。利用适应度函数计算第i只鸽子位置的适应度值，更新全局最优量子位置，即对于第i只鸽子，若否则，

步骤八：判断地标算子操作的循环次数是否达到要求的迭代次数：若未达到，令t＝t+1，返回步骤六；否则，输出鸽群全局最优位置b^t+1。

步骤九：更新n个角度的搜索区间。

在第(k+1)次快拍，n个角度搜索区间更新为其中为第(k+1)个快拍第n维的中心值，即δ为遗传因子，ω为搜索区间中影响收敛速度的收敛因子，为第k个快拍第n维的估计值，d为搜索空间的搜索半径，n＝1,2,…,n。

步骤十：判断是否达到最大快拍数：若未达到，继续获取下一个快拍采样经模式空间变换后的数据y(k+1)＝[y1(k+1),…,ym(k+1)]^t，进而协方差矩阵的更新方程为其中令k＝k+1，返回步骤三；否则，根据得到的所有快拍采样下动态目标的估计值，输出动态测向的结果。

本发明的效果在于：

基于圆阵模式空间实现了适用于非相干信源和相干信源的360度全方位动态测向，具有更广泛的适用性。

设计了连续量子鸽群计算方法进而能对阵列的极大似然方程进行高精度求解，具有收敛速度快，收敛精度高的优点。

通过更新搜索中心、搜索区间和迭代次数，可减少搜索时间，提高搜索速度和跟踪速度，具有良好的实时性。

附图说明

图1为基于连续量子鸽群的360度全方位动态测向方法的示意图；

图2为连续量子鸽群计算方法的流程图；

图3为均匀圆阵示意图；

图4为信噪比snr＝10db时两个独立信源的动态测向结果；

图5为信噪比snr＝10db时两个相干信源的动态测向结果；

图6为信噪比snr＝10db时三个独立信源的动态测向结果；

图7为信噪比snr＝10db时两个相干信源一个独立信源的动态测向结果；

图8为信噪比snr＝10db时三个相干信源的动态测向结果。

具体实施方式

下面结合附图对本发明的具体实施方式作进一步说明：

图1、图2分别为圆阵模式空间下基于连续量子鸽群的360度全方位动态测向方法的示意图和连续量子鸽群计算方法的流程图。本发明能够完成高斯噪声环境下的全方位360度测向，设计了连续量子鸽群计算方法来计算阵列的极大似然方程，通过不断更新搜索中心、搜索区间和迭代次数，减少搜索时间，提高搜索速度和跟踪速度，以满足性能要求。本发明解决问题所采用的方案步骤如下：

步骤一：获取阵列接收的快拍数据，定义变换矩阵t，计算阵列信号经模式变换后的协方差矩阵。

假设n个窄带远场信源以方位角θn，俯仰角方向入射到一个由m个阵元构成的半径为r的均匀圆阵上，n＝1,2,…,n，入射波长为λ，那么阵列接收的第k次快拍数据的数学模型为式中，k为最大快拍数，x(k)＝[x1(k),x2(k),…,xm(k)]^t为m×1维的阵列接收的快拍数据矢量，为阵列m×n维的流形矩阵，θ＝[θ1,θ2,…,θn]和为信源的方向矢量，是流形矩阵的第n个导向矢量，其中s(k)为n×1维的信号矢量，n(k)为m×1维的噪声矢量，阵列噪声假定为空间和时间均独立的高斯白噪声。这里只讨论所有信源都与阵列共面的情形，即所有的俯仰角均为90°。

均匀圆阵可激发的最大模式为向下取整函数。变换矩阵t定义为式中j＝diag{j^-lj-l(-β),…,j^ljl(-β)}，其中jl(·)为l阶第一类贝塞尔函数，l＝-l,…,0,…,l。式中f^h＝[w-l,w-l+1,…,wl]，其中m＝1,2,…,m-1。用t左乘x(k)可得y(k)＝tx(k)＝bs(k)+tn(k)，其中阵列流形其具有vandermonde矩阵结构，此时均匀圆阵已经转换为虚拟线阵。第1次快拍采样数据经过模式空间变换后的协方差矩阵为

步骤二：初始化搜索区间。

在第k次快拍，n个角度的搜索区间定义为其中un(k)和gn(k)分别为第k次快拍第n维角度搜索区间的上限和下限，其初始值分别取搜索区间定义域的上限和下限，n＝1,2,…,n。连续量子鸽群的最大迭代次数取第k次快拍所有维搜索区间上限与下限之差整倍数的最大值，即其中取整数，为向下取整函数。

步骤三：初始化鸽群每只鸽子的量子位置和速度，计算适应度值，更新全局最优量子位置。

鸽群中共有p只鸽子，每只鸽子的量子位置在量子域[0,1]内随机初始化，速度在[-vmax,vmax]内随机初始化。第i只鸽子量子位置定义为其中t为当前迭代次数，为第i只鸽子的第n维量子位置，且第i只鸽子的位置为量子位置的相应映射态，映射关系为第i只鸽子的速度为且

评估第i只鸽子位置状态的适应度函数为其中为的映射矩阵，trace()为求矩阵的迹。利用此适应度函数计算每只鸽子位置的适应度值，选择适应度值最大的解的量子位置作为全局最优量子位置，记作其中为到第t次迭代为止鸽群最优量子位置的第n维，n＝1,2,…,n，其相应映射态为全局最优鸽子位置

步骤四：根据地图和指南针算子更新每只鸽子的速度和量子位置及其相应映射态，计算适应度值，更新全局最优量子位置。

第i只鸽子速度的更新方程为其中i＝1,2,…,p，n＝1,2,…,n，r是地图和指南针因数，取值范围为(0,1)，是(0,1)之间的随机数。如果速度超出边界值，则将其限制在边界，即若若

第i只鸽子量子位置更新方程为同时按照映射方程获得量子位置的相应映射态。

利用适应度函数计算第i只鸽子位置的适应度值，更新全局最优量子位置，即对于第i只鸽子，若否则，

步骤五：判断地图和指南针算子操作的循环次数是否达到所要求的迭代次数：若未达到，令t＝t+1，返回步骤四继续循环；否则，令t＝t+1，停止地图和指南针算子操作，进入地标算子操作。

步骤六：进入地标算子操作，选择适应度值较优的前一半作为当前鸽群并计算当前鸽群的量子位置中心。

按照适应度值从大到小的规则对鸽子进行排序，当前代鸽子数量为p^t+1＝p^t/2，选择适应度值较优的前一半作为当前鸽群，计算当前鸽群的量子位置中心的第n维，即并将其作为参考方向，i＝1,2,…,p，n＝1,2,…,n。

步骤七：更新每只鸽子的量子位置及其相应映射态，计算适应度值，更新全局最优量子位置。

按照方程更新每只鸽子的量子位置，其中i＝1,2,…,p，n＝1,2,…,n，是(0,1)之间的随机数。按照映射方程获得量子位置的相应映射态。利用适应度函数计算第i只鸽子位置的适应度值，更新全局最优量子位置，即对于第i只鸽子，若否则，

步骤八：判断地标算子操作的循环次数是否达到要求的迭代次数：若未达到，令t＝t+1，返回步骤六；否则，输出鸽群全局最优位置b^t+1。

步骤九：更新n个角度的搜索区间。

在第(k+1)次快拍，n个角度搜索区间更新为其中为第(k+1)个快拍第n维的中心值，即δ为遗传因子，ω为搜索区间中影响收敛速度的收敛因子，为第k个快拍第n维的估计值，d为搜索空间的搜索半径，n＝1,2,…,n。

步骤十：判断是否达到最大快拍数：若未达到，继续获取下一个快拍采样经模式空间变换后的数据y(k+1)＝[y1(k+1),…,ym(k+1)]^t，进而协方差矩阵的更新方程为其中令k＝k+1，返回步骤三；否则，根据得到的所有快拍采样下动态目标的估计值，输出动态测向的结果。

图3中，均匀圆阵的m个阵元均匀分布在半径为r的圆周上，采用球坐标系来表示来波方向，坐标原点位于均匀圆阵的圆心位置上，来波的俯仰角是来波方向与z轴的夹角且方位角θ是来波方向在xoy平面内的投影与x轴的夹角且θ∈[0,360°]。

图4至图8所示的仿真中，均匀圆阵的阵元数m＝16，阵元间距d＝0.3λ，半径为噪声为空间和时间均独立的高斯白噪声。波达方向估计时的参数如下：信噪比snr＝10db，最大快拍数k＝300，每个方向的初始搜索区间均为[0°,360°]，收敛因子ω＝0.995，搜索半径d＝3，遗传因子δ＝0.8，协方差矩阵更新方程中的μ＝0.95。连续量子鸽群计算方法中鸽群的鸽子数量p＝40，地图和指南针算子r＝0.3，鸽子速度最大值vmax＝0.15。最大迭代次数取第k次快拍所有维搜索区间上限与下限之差整数倍的最大值，倍数地图和指南针算子操作迭代次数取最大迭代次数的地标算子操作的迭代次数取最大迭代次数的假设有2个信号源分别从θ1(k)＝[210+5sin(2πk/200)]°，θ2(k)＝[100+5sin(2πk/200)]°方向入射到均匀圆阵上，这两个独立信号源和相干信号源的动态测向结果分别如图4和图5所示。假设有3个信号源分别从θ1(k)＝[210+5sin(2πk/200)]°，θ2(k)＝[150+5sin(2πk/200)]°，θ3(k)＝[100+5sin(2πk/200)]°方向入射到均匀圆阵上，这三个独立信号源，两个相干源一个独立源和三个相干信号源的动态测向结果分别如图6、图7、图8所示。通过仿真图可以看出，仿真环境下本发明所设计的方法可有效地对动态来波方向进行估计和跟踪，解决了动态目标360度全方位测向难题，同时在高斯噪声环境下无论是非相干信源还是相干信源都能够被准确的跟踪也说明了本发明所设计的方法具有广泛的适用性。

以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。