冲击噪声下基于量子帝王蝶的圆阵模式空间动态测向方法与流程

本发明涉及一种基于量子帝王蝶的圆阵模式空间动态测向方法，属于阵列信号处理领域。

背景技术：

测向也称为波达方向角(directionofarrival，doa)估计一直是阵列信号处理领域的热点内容，在通信、雷达和声纳等系统中有着广泛应用。均匀圆阵(uniformcirculararray,uca)由于具有特殊的阵列结构和良好的测向性能，在许多测向应用中备受青睐。均匀圆阵的阵列结构使得其阵列流行比较复杂，不具备均匀线阵的范德蒙德矩阵形式，许多适用于线阵的测向算法无法移植应用到圆阵的测向应用中。

为了使应用于线阵的优良doa估计算法可以应用到圆阵，以及减少doa估计的运算量，模式空间变换算法估计受到了学者的广泛关注。模式空间变换即通过变换将圆阵的阵列流行变换成类似线阵的形式，实现输入信号的波达方向估计。

现有技术文献检索发现，高书彦等在《电子与信息学报》(2007,vol.29,no.12,pp.2832-2835)发表的“基于均匀圆阵的模式空间矩阵重构算法”中提出了一种不受信号相关性影响的doa估计方法，并在一定程度上减少了计算量。但是该方法损失了阵列孔径，估计精度不高，不能动态测向，且在冲击噪声环境下测向失效。赵大勇等在《山东大学学报(工学版)》(2010,vol.40,no.1,pp.133-138)发表的“冲击噪声背景下的动态doa跟踪”中对粒子群算法进行改进，研究了基于最大似然算法的动态doa估计方法，该方法在避免分数低阶矩矩阵重复分解的同时，在一定程度上降低了多维搜索的计算量。但是该方法仅能在180°角度范围内进行方向估计，且使用的低阶矩阵抗冲击噪声能力有待提升，并且该方法所使用的粒子群算法容易陷入“早熟”收敛，需要较大的种群规模和迭代次数，求解精度不高。

已有的文献表明，现有动态测向方法大多使用子空间跟踪类方法，这类方法具有计算量小，实时性好的优点，但是其在低信噪比的情况下估计性能较差，且本身并不具备解相干的能力，需要运用空间平滑等解相干技术对协方差矩阵进行预处理，步骤繁琐，且在冲击噪声背景下测向性能恶化严重，甚至失效。本发明设计了一种冲击噪声下基于量子帝王蝶的圆阵模式空间动态测向方法，在冲击噪声环境下，设计了一种指数核协方差矩阵，利用圆阵模式空间的指数核协方差矩阵极大似然动态测向方法，对信源进行动态测向，并使用量子帝王蝶搜索机制在搜索区间内搜索指数核协方差矩阵极大似然方程估计的最优角度，通过逐步缩小搜索空间进而减小计算量。本发明所设计的动态测向方法，在高斯噪声、弱冲击噪声和强冲击噪声等复杂环境下均能有效测向，在低信噪比，小快拍数和相干信源的情况下具有优越的doa估计性能。

技术实现要素：

针对上述现有技术，本发明要解决的技术问题是提供一种在高斯噪声、弱冲击噪声和强冲击噪声等复杂环境下均能有效测向，在低信噪比，小快拍数和相干信源的情况下具有优越的doa估计性能的冲击噪声下基于量子帝王蝶的圆阵模式空间动态测向方法。

为解决上述技术问题，本发明一种冲击噪声下基于量子帝王蝶的圆阵模式空间动态测向方法，包括以下步骤：

步骤一：获取阵列接收的快拍数据，定义变换矩阵构建阵列信号经模式空间变换后的指数核协方差矩阵，构造基于指数核协方差矩阵的极大似然测向方程，具体为：

假设第n个窄带远场信源从方位角θn，俯仰角入射到一个由m个阵元构成的半径为r的均匀圆阵上，n＝1,2,...,n，入射波长为λ，阵列接收的第k次快拍数据的数学模型为x(k)＝[x1(k),x2(k),...,xm(k)]^t为m×1维阵列接收快拍数据矢量，为m×n维的阵列流形矩阵，θ＝[θ1,θ2,...,θn]和分别为n个信源的方位角和俯仰角矢量，是阵列流形矩阵的第n个导向矢量，其中，n＝1,2,...,n；j为虚数单位；t表示转置，s(k)为n×1维的信号矢量，n(k)为m×1维服从特征指数为α的标准sαs分布的冲击噪声矢量，全部信源均与阵列共面，即俯仰角

均匀圆阵可激发的最大模式为向下取整函数，β＝2πr/λ；定义变换矩阵为式中j＝diag{j^-lj-l(-β),...,j^ljl(-β)}，其中，jl(·)为l阶第一类贝塞尔函数；diag{θ}表示由矢量θ中的元素构成的对角矩阵，式中f＝[w-l,w-l+1,...,wl]，其中，l＝-l,...,0,...,l；h表示共轭转置；用左乘x(k)得到y(k)＝[y1(k),y2(k),...,y2l+1(k)]^t，其中阵列流形具有范德蒙矩阵结构，则均匀圆阵转换为虚拟均匀线阵；第1次快拍采样数据经过模式空间变换后的指数核协方差矩阵为矩阵中第i行第d列个元素其中，i＝1,2,...,2l+1；d＝1,2,...,2l+1；η为指数核常数；*表示共轭；构造基于指数核协方差矩阵的极大似然测向方程为其中pb(θ)＝b(θ)[b^h(θ)b(θ)^-1]b^h(θ)为正交投影矩阵，为第k次快拍采样的指数核协方差矩阵，argmax()表示寻找具有最大函数值的变量，trace表示矩阵的求迹；

步骤二：初始化帝王蝶的搜索区间，具体为：

在第k次快拍，n个方位角度的搜索区间定义为其中，un(k)和gn(k)分别为第k次快拍第n维方位角度搜索区间的上限和下限，第1次快拍搜索区间的初始值分别取搜索区间定义域的上限和下限；n＝1,2,...,n；

步骤三：初始化帝王蝶种群中每只帝王蝶的量子位置、临时量子位置和位置，构造适应度函数，计算每只帝王蝶所在位置的适应度，确定帝王蝶种群的全局最优量子位置，确定量子帝王蝶搜索机制的终止迭代次数，具体为：

帝王蝶的种群规模为通过随机方法在量子域内初始化只帝王蝶的量子位置和临时量子位置，第t次迭代第h只帝王蝶的量子位置为其中，第t次迭代第h只帝王蝶的临时量子位置为其中，第t次迭代第h只帝王蝶的位置为是第t次迭代第h只帝王蝶量子位置的映射态，映射方式为其中，t为迭代次数；定义第t次迭代第h只帝王蝶的量子旋转角为其中，对于第k次快拍的角度估计，量子帝王蝶搜索机制的终止迭代次数取第k次快拍所有维搜索区间上限与下限之差整倍数的最大值，即其中，取整数；为向下取整函数，初始时迭代次数设为1；

第t次迭代第h只帝王蝶所在位置的适应度函数为根据适应度函数计算每只帝王蝶所在位置的适应度函数值，将帝王蝶种群至第t代为止所搜索到的适应度值最大的位置记录为全局最优位置对应的全局最优量子位置为

步骤四：将帝王蝶种群分为两个子群，使用两种不同的策略更新子群1和子群2中帝王蝶的临时量子位置，具体为：

选择只帝王蝶组成子群1，其中，p1为常数，表示子群1中帝王蝶所占比例，余下的只帝王蝶构成子群2；

在子群1中，在第t次迭代，如果第h1只帝王蝶第n维临时量子位置的更新方式为否则，第h1只帝王蝶第n维临时量子位置的更新方式为其中，和均为[0,1]之间均匀分布的随机数；p2为常数，表示迁移周期；为在子群1中随机选择的帝王蝶标号；为在子群2中随机选择的帝王蝶标号；

在子群2中，如果第h2只帝王蝶第n维临时量子位置的更新方式为否则，第h2只帝王蝶第n维临时量子位置的更新方式为其中，为[0,1]之间均匀分布的随机数；p3为常数，表示迁移率；为在子群2中随机选择的帝王蝶标号；为帝王蝶种群全局最优量子位置的第n维；

步骤五：重组两个子群为一个种群，更新每只帝王蝶的量子位置和位置，计算每只帝王蝶所在新位置的适应度值，更新全局最优量子位置和全局最优位置，具体为：

产生一个[0,1]之间均匀分布的随机数如果p4为帝王蝶调整概率，使用模拟量子旋转门更新第h只帝王蝶第n维量子位置：否则，其中，为第h只帝王蝶第n维量子旋转角；为[0,1]之间均匀分布的随机数；根据映射关系获得第h只帝王蝶第n维位置为

根据适应度函数计算所有帝王蝶所在新位置的适应度，根据适应度函数值更新帝王蝶种群的全局最优位置和全局最优量子位置，从个量子位置即和中选择个优秀的量子位置作为下一次迭代帝王蝶的演进初始量子位置，且记作

步骤六：判断是否达到所要求的最大迭代次数：若未达到，令迭代次数加1，返回步骤四继续循环；否则，输出帝王蝶种群全局最优位置和全局最优量子位置并进入下一步；

步骤七：更新下一次快拍n个方位角度的搜索区间，判断是否达到最大快拍数：若未达到，继续获取下一次快拍采样经模式空间变换后的数据，更新指数核协方差矩阵，令快拍数加1，返回步骤三；否则，根据得到的所有快拍采样下动态目标的估计值，输出动态测向结果，具体为：

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

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

本发明有益效果：与现有技术相比本发明具有以下优点：

(1)本发明在冲击噪声环境下，经过模式空间变换将均匀圆阵转换为虚拟均匀线阵，并使用基于指数核协方差矩阵的极大似然动态测向方法对动态目标进行角度估计，所设计方法能够对独立源、相干源和混合源进行动态测向，动态测向性能优，且可进行冲击噪声环境下动态目标的360°全方位跟踪。

(2)本发明所设计的动态测向方法，设计了指数核协方差矩阵，在低信噪比、小快拍数、高斯噪声、弱冲击噪声和强冲击噪声情况下均能有效动态测向，应用范围广泛。

(3)设计了量子帝王蝶搜索方法进而能对阵列的指数核协方差矩阵极大似然方程进行高精度求解，具有收敛速度快，收敛精度高的优点。

本发明搜索速度快，跟踪精度高，可在高斯噪声、弱冲击噪声和强冲击噪声下动态测向，应用范围广泛。

附图说明

图1为冲击噪声下基于量子帝王蝶的圆阵模式空间动态测向方法示意图；

图2为特征指数为0.95，广义信噪比snr＝15db时两个独立信源的动态测向结果；

图3为特征指数为1.8，广义信噪比snr＝15db时三个独立信源的动态测向结果；

图4为特征指数为1.8，广义信噪比snr＝15db时三个相干信源的动态测向结果；

图5为特征指数为1.8，广义信噪比snr＝15db时三个混合信源的动态测向结果；

图6为高斯噪声下即特征指数为2，信噪比snr＝10db时三个独立信源的动态测向结果。

具体实施方式

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

如图1所示，本发明技术方案包括如下步骤：

步骤一：获取阵列接收的快拍数据，定义变换矩阵构建阵列信号经模式空间变换后的指数核协方差矩阵，构造基于指数核协方差矩阵的极大似然测向方程。

假设第n个窄带远场信源从方位角θn，俯仰角入射到一个由m个阵元构成的半径为r的均匀圆阵上，n＝1,2,...,n，入射波长为λ，那么阵列接收的第k次快拍数据的数学模型为x(k)＝[x1(k),x2(k),...,xm(k)]^t为m×1维阵列接收快拍数据矢量，为m×n维的阵列流形矩阵，θ＝[θ1,θ2,...,θn]和分别为n个信源的方位角和俯仰角矢量，是阵列流形矩阵的第n个导向矢量。其中，n＝1,2,...,n；j为虚数单位；t表示转置。s(k)为n×1维的信号矢量，n(k)为m×1维服从特征指数为α的标准sαs分布的冲击噪声矢量。这里讨论的是所有信源都与阵列共面的情况，即俯仰角

均匀圆阵可激发的最大模式为向下取整函数；β＝2πr/λ。定义变换矩阵为式中j＝diag{j^-lj-l(-β),...,j^ljl(-β)}，其中，jl(·)为l阶第一类贝塞尔函数；diag{θ}表示由矢量θ中的元素构成的对角矩阵，式中f＝[w-l,w-l+1,...,wl]，其中，l＝-l,...,0,...,l；h表示共轭转置。用左乘x(k)可得y(k)＝[y1(k),y2(k),...,y2l+1(k)]^t，其中阵列流形具有vandermonde矩阵结构，此时均匀圆阵已经转换为虚拟均匀线阵。第1次快拍采样数据经过模式空间变换后的指数核协方差矩阵为矩阵中第i行第d列个元素其中，i＝1,2,...,2l+1；d＝1,2,...,2l+1；η为指数核常数；*表示共轭。构造基于指数核协方差矩阵的极大似然测向方程为其中pb(θ)＝b(θ)[b^h(θ)b(θ)^-1]b^h(θ)为正交投影矩阵，为第k次快拍采样的指数核协方差矩阵，argmax()表示寻找具有最大函数值的变量，trace表示矩阵的求迹。

步骤二：初始化帝王蝶的搜索区间。

在第k次快拍，n个方位角度的搜索区间定义为其中，un(k)和gn(k)分别为第k次快拍第n维方位角度搜索区间的上限和下限，第1次快拍搜索区间的初始值分别取搜索区间定义域的上限和下限；n＝1,2,...,n。

步骤三：初始化帝王蝶种群中每只帝王蝶的量子位置、临时量子位置和位置，构造适应度函数，计算每只帝王蝶所在位置的适应度，确定帝王蝶种群的全局最优量子位置，确定量子帝王蝶搜索机制的终止迭代次数。

帝王蝶的种群规模为通过随机方法在量子域内初始化只帝王蝶的量子位置和临时量子位置，第t次迭代第h只帝王蝶的量子位置为其中，第t次迭代第h只帝王蝶的临时量子位置为其中，第t次迭代第h只帝王蝶的位置为是第t次迭代第h只帝王蝶量子位置的映射态，映射方式为其中，t为迭代次数。定义第t次迭代第h只帝王蝶的量子旋转角为其中，对于第k次快拍的角度估计，量子帝王蝶搜索机制的终止迭代次数取第k次快拍所有维搜索区间上限与下限之差整倍数的最大值，即其中，取整数；为向下取整函数。初始时迭代次数设为1。

第t次迭代第h只帝王蝶所在位置的适应度函数为根据适应度函数计算每只帝王蝶所在位置的适应度函数值，将帝王蝶种群至第t代为止所搜索到的适应度值最大的位置记录为全局最优位置对应的全局最优量子位置为

步骤四：将帝王蝶种群分为两个子群，使用两种不同的策略更新子群1和子群2中帝王蝶的临时量子位置。

选择只帝王蝶组成子群1，其中，p1为常数，表示子群1中帝王蝶所占比例。余下的只帝王蝶构成子群2。

在子群1中，在第t次迭代，如果第h1只帝王蝶第n维临时量子位置的更新方式为否则，第h1只帝王蝶第n维临时量子位置的更新方式为其中，和均为[0,1]之间均匀分布的随机数；p2为常数，表示迁移周期；为在子群1中随机选择的帝王蝶标号；为在子群2中随机选择的帝王蝶标号。

在子群2中，如果第h2只帝王蝶第n维临时量子位置的更新方式为否则，第h2只帝王蝶第n维临时量子位置的更新方式为其中，为[0,1]之间均匀分布的随机数；p3为常数，表示迁移率；为在子群2中随机选择的帝王蝶标号；为帝王蝶种群全局最优量子位置的第n维。

步骤五：重组两个子群为一个种群，更新每只帝王蝶的量子位置和位置，计算每只帝王蝶所在新位置的适应度值，更新全局最优量子位置和其相应的全局最优位置。

产生一个[0,1]之间均匀分布的随机数如果p4为帝王蝶调整概率，使用模拟量子旋转门更新第h只帝王蝶第n维量子位置：否则，其中，为第h只帝王蝶第n维量子旋转角；为[0,1]之间均匀分布的随机数；根据映射关系获得第h只帝王蝶第n维位置为

根据适应度函数计算所有帝王蝶所在新位置的适应度，根据适应度函数值更新帝王蝶种群的全局最优位置和全局最优量子位置。从个量子位置即和中选择个优秀的量子位置作为下一次迭代帝王蝶的演进初始量子位置，且记作

步骤六：判断是否达到所要求的最大迭代次数：若未达到，令迭代次数加1，返回步骤四继续循环；否则，输出帝王蝶种群全局最优位置和全局最优量子位置并进入下一步。

步骤七：更新下一次快拍n个方位角度的搜索区间，判断是否达到最大快拍数：若未达到，继续获取下一次快拍采样经模式空间变换后的数据，更新指数核协方差矩阵，令快拍数加1，返回步骤三；否则，根据得到的所有快拍采样下动态目标的估计值，输出动态测向结果。

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

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

在特征指数为2即高斯噪声下snr代表信噪比，其他情况下snr代表广义信噪比。冲击噪声下基于量子帝王蝶的圆阵模式空间动态测向方法的参数设置如下：

在圆阵模式空间测向系统中，均匀圆阵的阵元数m＝16，阵元间距圆阵半径为波达方向估计参数设置情况：最大快拍数k＝500，每个方位角度的初始搜索区间均为[0°,360°]，收敛因子ω＝0.995，搜索半径遗传、因子更新因子μ＝0.95。量子帝王蝶搜索机制的参数设置情况：种群规模p1＝0.5，迁移周期p2＝1.2，迁移率帝王蝶调整概率

图2中，两个独立信源从θ1(k)＝[120+10sin(2πk/500)]°和θ2(k)＝[-120+10sin(2πk/500)]°方向入射到均匀圆阵，冲击噪声特征指数为0.95，广义信噪比为15db。从仿真图2中可以看出，在强冲击噪声下，所设计的方法在大多数情况下能有效的对动态来波方向进行跟踪。

图3中，三个独立信源从θ1(k)＝[120+10sin(2πk/500)]°，θ2(k)＝[0+10sin(2πk/500)]°和θ3(k)＝[-120+10sin(2πk/500)]°方向入射到均匀圆阵，冲击噪声特征指数为1.8，广义信噪比为15db。从仿真图3中可以看出，所设计的方法能有效的对独立信源进行动态跟踪。

图4中，三个相干信源从θ1(k)＝[120+10sin(2πk/500)]°，θ2(k)＝[0+10sin(2πk/500)]°和θ3(k)＝[-120+10sin(2πk/500)]°方向入射到均匀圆阵，冲击噪声特征指数为1.8，广义信噪比为15db。从仿真图4中可以看出，所设计的方法能有效的对相干信源进行动态跟踪。

图5中，三个混合信源从θ1(k)＝[120+10sin(2πk/500)]°，θ2(k)＝[0+10sin(2πk/500)]°和θ3(k)＝[-120+10sin(2πk/500)]°方向入射到均匀圆阵，其中θ1(k)与θ2(k)相干，θ2(k)与θ3(k)独立，冲击噪声特征指数为1.8，广义信噪比为15db。从仿真图5中可以看出，所设计的方法能有效的对混合信源进行动态跟踪。

图6中，在高斯噪声下，三个相干信源从θ1(k)＝[120+10sin(2πk/500)]°，θ2(k)＝[0+10sin(2πk/500)]°和θ3(k)＝[-120+10sin(2πk/500)]°方向入射到均匀圆阵，信噪比为10db。从仿真图6中可以看出，所设计的方法在低信噪比下能有效的对动态来波方向进行跟踪。