CN116582158A - 一种大规模mimo方阵信源数与波达方向联合估计方法 - Google Patents

Abstract

本发明公开了一种大规模MIMO方阵信源数与波达方向联合估计方法包括：建立冲击噪声环境下MassiveMIMO方阵模型；构造基于拉普拉斯核相关熵的加权信号子空间拟合方程；利用分段式思想，化简得到目标函数；初始化个体量子位置，获得全局最优量子位置；初始化量子冰晶能量值，确定临时湖中心位置；更新能量值和历史量子位置空间；更新量子位置；根据轮盘赌选择产生新一代量子冰晶，更新全局最优量子位置；判断是否达到最大迭代次数，若未达到，返回步骤五；判断该信源是否存在，若存在，返回步骤三，否则输出信源数及相应波达方向。本发明具有鲁棒性、高精度的特点和更广泛的应用范围。

Description

一种大规模MIMO方阵信源数与波达方向联合估计方法

技术领域

本发明属于阵列信号处理领域，涉及一种大规模MIMO方阵信源数与波达方向联合估计方法，特别是一种冲击噪声环境下基于量子冰晶优化机制的Massive MIMO方阵信源数与波达方向联合估计方法。

背景技术

DOA估计技术作为阵列信号处理的重点研究方向之一，在无线通讯、雷达和定位等多个领域有着非常广泛的应用。尤其是在第五代移动通信系统中，大规模MIMO技术被认为是5G关键技术之一，受到业内人士的广泛关注。在大规模MIMO系统中，二维DOA估计可以同时对俯仰角和方位角同时进行估计，能够为移动台提供更精准的位置估计，且能够使波束赋形的指向性更强，是3D-MIMO技术得以实施的关键，因此对DOA估计算法进行深入研究具有重要意义。但目前的二维DOA估计算法中，多是基于线阵的MUSIC算法和ESPRIT算法扩展而来，主要通过方位角和俯仰角二维角度域的联合搜索来获取DOA估计值。这些方法不仅需要信源个数这个先验知识，往往还需要较大的快拍数保证估计精度和准确率，在面对大规模MIMO系统的庞大天线阵列时，一些高精度谱峰搜索类方法更是具有非常高的计算复杂度，并且对于相干信源，不得不损失阵列孔径才能进行解相干操作。所以针对于大规模MIMO的面阵系统，在未知信源个数情况下，设计一种抗冲击噪声、快速、高效和低计算复杂度的二维DOA估计方法是十分有必要的。

经对现有文献的检索发现，Wu W等在《IEEE Intemational Conference onSignal,Information and DataProcessing》上发表的“A low copmlexity 2D DOAestimation algorithm for massive MIMO systems”中提出基于传统的Capon DOA估计算法，通过离散傅里叶变换先对信号的DOA初步估计，然后利用降维Capon算法在小范围内进一步估计，该算法可以适当的降低计算的复杂度，但是该算法的分辨力有限，精度还需要进一步改进提高。Zheng Z等在《IEEE Transactions on Vehicular Technology》上发表的“Efficient beam space-based algorithm for two-dimensional DOA estimation ofincoherently distributed sources in massive MIMO systems”提出利用广义阵列流型矩阵进行分布源的波束空间变换，构造波束空间的移不变结构，避免了谱峰搜索的操作，以此来减少矩阵运算的维度，在计算复杂度方面具有较大的优势，但是其精度性能一般，并且在冲击噪声的情况下算法性能恶化严重。Zhu Y等在《IEEE International Conference onCommunications》上发表“DoA estimation and capacity analysis for 2D activemassive MIMO systems”提出了一种基于MUSIC算法的二维DOA角度估计方法，由于MUSIC算法需要在全域进行谱峰搜索，因此尽管文献采用分段式的方法进行角度估计，但在大规模天线阵列系统背景下，所提的算法仍具有较高的运算复杂度，且不可避免的出现量化误差。Wang A等在《IEEE GLOBECOM Workshops》上发表“Low complexity direction of arrival(DoA)estimation for 2D massive MIMO systems”提出采用Unitary ESPRIT算法进行大规模天线阵的DOA角度估计，Unitary ESPRIT算法将方向矩阵和协方差矩阵等转换成实值巧妙地化简了虚部，从而降低了运算复杂度，但是文章对于俯仰角和方位角的联合分析并不充分，估计的精度受限且需要大的快拍数。可以看出已有的大部分研究都集中在降低方法的计算复杂度上，以牺牲估计精度为代价，然而提高测向方法在低信噪比下的准确度和成功概率也是重要的研究方向，并且上述方法都是假设在信源个数已知的情况下，但实际应用中信源个数并不是已知的，因此还需要先进行额外的信源个数估计，而且信源数估计是落后于波达方向估计的，影响了高精度波达方向估计的实用化进程，且都没有考虑冲击噪声对测向结果的影响。

此外，在已有文献中，这些测向方法大多假设背景噪声为高斯噪声，通过利用二阶或高阶累积量进行分析就可以获得理想的结果。然而实际信号传输中噪声更多不是服从高斯分布的，如海杂波噪声、大气噪声和无线信道噪声等，这些噪声可利用Alpha稳定分布建模，它与高斯噪声模型失配，使得传统的基于二阶或高阶累积量的算法失效，因此研究在冲击噪声环境下的大规模MIMO系统中对信源个数以及二维波达方向进行去相干、小快拍、高精度和高鲁棒性的联合估计具有重要的意义和价值。

发明内容

针对上述现有技术，本发明要解决的技术问题是提供一种有效性和鲁棒性更高的基于量子冰晶优化机制的大规模MIMO方阵信源数与波达方向联合估计方法，在冲击噪声、小快拍和相干源存在的情况下，实现Massive MIMO系统信源数与二维波达方向的联合估计。

为解决上述技术问题，本发明的一种大规模MIMO方阵信源数与波达方向联合估计方法，包括：

步骤一：建立冲击噪声环境下的Massive MIMO方阵模型，获取阵列接收的快拍数据；

步骤二：利用接收数据构造基于拉普拉斯核相关熵的低阶协方差矩阵，得到基于拉普拉斯核相关熵的通用加权信号子空间拟合方程；

步骤三：利用分段式逐个求解方法，得到待搜索来波最终的目标函数；

步骤四：初始化估计第l个来波的量子冰晶机制的冰晶个体数及每个量子冰晶的量子位置和相应的映射态位置并获得全局最优量子位置；

步骤五：初始化量子冰晶能量值，确定一个临时湖中心位置，与外界进行能量交换；

步骤六：更新量子冰晶能量，能量值最低的N_q个量子冰晶开始沉淀结冰，加入到已有的冰晶壳中，并更新历史量子位置空间

步骤七：根据量子演化规则，使用模拟量子旋转门更新未沉淀的量子冰晶的量子位置和相应的映射态位置；

步骤八：计算每个量子冰晶的适应度值，根据轮盘赌选择策略选择量子冰晶来产生新一代量子冰晶的最终位置，然后更新全局最优量子位置；

步骤九：判断是否达到其最大迭代次数T_max，若未达到，令t=t+1，返回步骤五继续迭代；若达到最大迭代次数，则选择当前的全局最优量子位置的映射态位置作为最后的结果，并输出此时的映射态位置

步骤十：判断信源l是否存在：设置阈值若/>则判断空间中已无未知信源，然后输出空间中的已测得的未知信源的个数以及相应的二维波达方向的估计值，信源数为l-1；若/>则选择当前的映射态位置/>作为该信源二维波达方向的估计值，再令l=l+1，返回步骤三，并将得到的波达方向估计值作为求得的先验信息带入下一次求解的目标函数化简中，然后再一次搜索。

进一步的，步骤一所述阵列接收的快拍数据满足：

假设空间中存在B个窄带远场信源，分别以方位角θ_b和俯仰角方向入射到由M×N个阵元构成的Massive MIMO水平均匀方阵上，b=1,2,...,B，阵元间距为d，入射波长为λ，则阵列接收的第k次快拍采样数据的数学模型为：

其中k=1,2,...,K，K为最大快拍数，z(k)=[z₁(k),z₂(k),...,z_MN(k)]^T为MN×1维的阵列接收的第k次快拍数据矢量，上标T表示转置，为MN×B维的流型矩阵，其中θ＝[θ₁,θ₂,...,θ_B]和/>分别为信源的方位角矢量与俯仰角矢量，/>表示第b个信源的导向矢量，其中为流型矩阵关于x轴方向上的第b个信源的导向矢量，/>为流型矩阵关于y轴方向上的第b个信源的导向矢量，/>为Kronecker积，s(k)为B×1维信号矢量，n(k)为MN×1维服从SαS稳定分布的冲击噪声矢量。

进一步的，步骤二所述基于拉普拉斯核相关熵的低阶协方差矩阵和基于拉普拉斯核相关熵的通用加权信号子空间拟合方程具体为：

基于拉普拉斯核相关熵的低阶协方差矩阵R的第行/>列元素/>表示为：

式中，和/>分别表示接收的第k次快拍信号数据矢量的第/>维和第/>维，(·)^*表示共轭运算，η为核函数的核长；

已知空间中的未知信源个数B，然后对低阶协方差矩阵R进行特征分解：

其中，U_S是由B个较大特征值对应的特征向量张成的信号子空间，Σ_S是由B个较大特征值构成的对角阵，U_N是由剩余较小特征值对应的特征向量张成的噪声子空间，Σ_N是由剩余较小特征值构成的对角阵，(·)^H表示共轭转置运算，均匀方阵的正交投影矩阵为(·)^-1表示求逆运算；

则基于拉普拉斯核相关熵的加权信号子空间拟合方程通式为：

其中，为信号子空间的最优权矩阵，其中σ²表示噪声功率，tr(·)为矩阵求迹函数。

进一步的，步骤三所述待搜索来波最终的目标函数包括：

其中，空间中未知的第l个来波标号为l，初始时未知来波标号l=1，为定义的归一化矢量，满足：

其中， 和/>分别表示已经搜索出的l-1维方位角参数矢量和俯仰角参数矢量，低阶协方差矩阵R分解为：/>其中/>由前l个较大特征值对应的特征向量张成的信号子空间，/> 是由前l个较大特征值构成的对角阵，σ²表示噪声功率。

进一步的，步骤四所述初始化估计第l个来波的量子冰晶机制的冰晶个体数及每个量子冰晶的量子位置和相应的映射态位置并获得全局最优量子位置包括：

设置量子冰晶的个体数为N_p，每个量子冰晶所具有的空间维数为最大的迭代次数为T_max，t代表迭代次数，对N_p个量子冰晶的量子位置在量子域[0,1]内随机初始化，则估计第l个信源时第t代第i个量子冰晶的量子位置定义为/>其中表示第t代第i个量子冰晶量子位置的第/>维，相应的映射态位置为其映射关系为/>其中 和/>分别表示第/>维变量的下边界和上边界，i=1,2,...,N_p，计算每个量子冰晶映射态位置的适应度值，在估计第l个信源时将第t代第i个量子冰晶映射态位置/>代入适应度函数表达式中，得到记录搜索到的指导第t代的量子冰晶中适应度值最大的量子位置记为全局最优量子冰晶位置/>其所对应的映射态位置记为/>

进一步的，步骤五所述初始化量子冰晶能量值，确定一个临时湖中心位置，与外界进行能量交换包括：

对于第t代第i个量子冰晶其初始能量为然后确定一个临时湖中心，采用角度联合距离的方式，动态调整湖水中心位置，定义第t代第i个量子冰晶的多样性度量方程为/>表示第t代第i个量子冰晶/>与第j个晶体之间的的最小夹角，i=1,2,...,N_p，j=1,2,...,N_p且i≠j，/>越大表示该晶体与其他晶体分隔越明显，其多样性也就越好，其中/>||·||₂表示向量的二阶范式，(·)^T表示转置；然后定义第t代第i个量子冰晶角度联合距离动态调整方程为/>其中定义Z(ρ)为角度影响因子，τ为控制系数，/>为第t代第i个量子冰晶的能量值，最后获得临时的湖水中心位置为/>

将湖水结冰过程总结为析出和冰冻两个阶段，首先判断其迭代次数t是否小于等于T_max/3，若是，则量子冰晶进入析出阶段，更新每个晶体能量变化值；否则，量子冰晶进入冰冻阶段，更新每个晶体能量变化值；

对于第t代第i个量子冰晶，定义其能量变化方程为其中主要分为三部分，第一部分/>为第t代第i个量子冰晶从湖水中心获取到的能量，/>为第t代第i个量子冰晶到湖水中心的距离，q为一个能量定值；第二部分/>表示第t代第i个量子冰晶与已完成结冰的冰晶之间进行的能量交换，O为一个能量定值；第三部分F_i ^t=C表示第t代第i个量子冰晶被风带走的能量，C为一个能量定值，λ₁、λ₂和λ₃均为比例系数。通过调整λ₁、λ₂和λ₃三个比例系数来区别处于析出阶段或是冰冻阶段。

进一步的，步骤六所述更新量子冰晶能量，能量值最低的N_q个量子冰晶开始沉淀结冰，加入到已有的冰晶壳中，并更新历史量子位置空间包括：

更新量子冰晶能量，第t+1代第i个量子冰晶的能量为能量更新后，能量值最低的N_q个量子冰晶沉淀结冰，同时记录这N_q个冰晶中适应度值最大的前/>个映射位置所对应的量子位置，组成历史量子位置空间/> 其中其相应的映射态位置所代表的空间为其中/>随着迭代的进行，不断地更新替换历史量子位置空间/>中的量子位置。

进一步的，步骤七所述根据量子演化规则，使用模拟量子旋转门更新未沉淀的量子冰晶的量子位置和相应的映射态位置包括：

定义第t+1代第i个量子冰晶的第维量子旋转角为/> 其中为转移因子，/>为第t+1代第i个晶体的能量值，μ为比例系数，/>为在[0,1]之间服从均匀分布的随机数，/>ν^t为第t代惯性权重系数，ν_max为ν^t的最大值，v_min为v^t的最小值，/>用于查找第t代/>中最接近/>的量子位置的第/>维，最后对每个量子冰晶的量子位置通过简化的模拟量子旋转门进行更新，则第t+1代第i个量子冰晶的第/>维量子位置更新为同时通过映射方程获得第t+1代第i个量子冰晶的第/>维映射态位置为/>

进一步的，步骤八：计算每个量子冰晶的适应度值，根据轮盘赌选择策略选择量子冰晶来产生新一代量子冰晶的最终位置，然后更新全局最优量子位置包括：

首先计算第i个量子冰晶个体被选中的概率对于第t+1代第i个量子冰晶，其被选中的概率为/> 表示第t+1代第i个量子冰晶的适应度值，然后确定第i个个体的累计概率/>对于第i个冰晶，其累积概率/>接着生成一个[0,1]之间服从均匀分布的随机数/>若/>大于等于累积概率/>而小于/>则第i个量子冰晶被选中来指导生成一个新的量子冰晶，其第/>维量子位置生成公式为 需要约束在[0,1]之间，超过边界值，取到边界，/>为权重因子，/>为在[0,1]之间服从均匀分布的随机数，然后计算其相应的映射态位置以此重复生成N_q个量子冰晶，补充更新过程中沉淀结冰的量子冰晶，保证量子冰晶个体数不变；最后将生成冰晶的映射态位置带入到适应度值公式中，计算对应的适应度值，并与原有更新后的量子冰晶的适应度值放在一起，选择其中的最大值与上一代的全局最优量子位置的适应度值进行比较，最终选择更大适应度值所对应的量子位置作为指导第t+1代量子冰晶的全局最优量子位置，记为/>其相应的映射态位置记为/>

进一步的，所述阈值的设置包括：

设置其中δ是通过计算迭代过程中所有映射位置的平均适应度，/>是设定的倍数，同时确定一个阈值下限/>取/>和/>中的较大值作为阈值。

本发明的有益效果：本发明将量子优化理论和冰晶优化机制结合进而得到量子冰晶优化机制，进一步提升冰晶算法的全局收敛性能。其中利用Massive MIMO天线阵接收的信息数据构造基于拉普拉斯核相关熵的低阶协方差矩阵，消除了二阶及以上矩抗冲击噪声能力弱的不足，设计出了基于拉普拉斯核相关熵的加权信号子空间拟合方程。而且加权信号子空间拟合方法在小快拍和低信噪比情况下测向性能优异，也不需要额外的解相干处理即可处理相干信源，唯一不足就是作为一个非线性多维搜索问题，具有极高的复杂度，尤其是应用在Massive MIMO系统中。因此通过利用分段式逐个求解的思想，每次搜索只对一个未知信源进行求解，而其他已求得信源参数保持不变，这样不仅将搜索区域从高维转换到了低维，克服大规模天线阵列计算量大的问题，极大地降低了计算的复杂度，而且还可以在信源个数未知的情况下，进行未知信源个数与波达方向的联合估计。并且基于量子编码和模拟量子旋转门设计了量子旋转角和量子演化方程，进而设计了有更快收敛速度的量子冰晶优化机制，进一步提高了已有DOA估计方法在搜索区间内搜索效率、搜索精度和抗干扰的能力，可以更快速地获得信源个数和二维波达方向的全局最优解。与现有技术相比，本发明具有以下特点：

(1)采用传统经典的方法，像MUSIC算法、ESPRIT算法或其他由此扩展算法，当面对实际应用中复杂的电磁环境，尤其在低信噪比、小快拍和冲击噪声背景下进行MassiveMIMO方阵的二维波达方向估计时，其鲁棒性与准确性多难以满足实际应用要求。本发明通过采用加权信号子空间拟合算法，设计的一种基于拉普拉斯核相关熵的加权信号子空间拟合方程，可以有效的消除冲击噪声的影响，并且在低信噪比、小快拍和相干源的情况下依旧进行有效的波达方向估计。

(2)采用传统方法进行Massive MIMO系统的波达方向估计时，不可避免会面临一个由于大规模的天线阵列所带来的巨大的计算量和非常高的计算复杂度的问题，由此往往会以牺牲一部分算法估计的精度为代价，来降低算法的计算复杂度，然而提高算法的估计准确度也是十分重要的，因此本发明中利用分段式逐个求解的思想，将加权信号子空间拟合方程进一步的转化，把高维搜索问题通过低维搜索方式进行了巧妙地解决，极大的降低计算的复杂度，同时还保持了加权信号子空间拟合算法具有的高精度的特点，并且本发明中将量子优化理论与冰晶优化机制相结合，基于量子编码和模拟量子旋转门设计的量子演化方程和量子冰晶优化机制具有更快的收敛速度，可以对所设计的目标函数进行快速高精度求解，能够再进一步降低运算量，提高求解的精度和鲁棒性，使波达方向估计结果具有更好的鲁棒性与收敛性。

(3)传统方法在进行波达方向估计时往往需要信源个数作为先验知识，这是算法实施的前提条件，因此往往需要先进行额外的未知信源个数估计。但本发明通过利用分段式逐个求解的思想，不仅不需要未知信源个数这个先验知识，并且还可以完成对未知信源个数及二维波达方向的联合估计，因此本发明具有更小的限制和更广泛的适用范围。

附图说明

图1是一种大规模MIMO方阵信源数与波达方向联合估计方法流程图

图2是信源数与二维波达方向联合估计示意图

图3是冲击噪声下两个信源30次波达方向估计示意图

图4是两种波达方向估计方法均方根误差与快拍数关系曲线

图5是两种波达方向估计方法均方根误差与广义信噪比关系曲线

具体实施方式

下面结合说明书附图和实施例对本发明做进一步说明。

本发明的基于量子冰晶优化机制的Massive MIMO方阵信源数与波达方向联合估计方法总体流程如图1所示，本发明包括如下步骤：

步骤一：建立冲击噪声环境下的Massive MIMO方阵模型，获取阵列接收的快拍数据。

假设空间中存在B个窄带远场信源，分别以方位角θ_b和俯仰角方向入射到由M×N个阵元构成的Massive MIMO水平均匀方阵上，b=1,2,...,B，阵元间距为d，入射波长为λ，那么阵列接收的第k次快拍采样数据的数学模型为/>其中k=1,2,...,K，K为最大快拍数，z(k)=[z₁(k),z₂(k),…,z_MN(k)]^T为MN×1维的阵列接收的第k次快拍数据矢量，上标T表示转置，/>为MN×B维的流型矩阵，其中θ＝[θ₁,θ₂,...,θ_B]和/>分别为信源的方位角矢量与俯仰角矢量，表示第b个信源的导向矢量，其中为流型矩阵关于x轴方向上的第b个信源的导向矢量，/>为流型矩阵关于y轴方向上的第b个信源的导向矢量，/>为Kronecker积，s(k)为B×1维信号矢量，n(k)为MN×1维服从SαS稳定分布的冲击噪声矢量。

步骤二：利用接收数据构造基于拉普拉斯核相关熵的低阶协方差矩阵，得到基于拉普拉斯核相关熵的通用加权信号子空间拟合方程。

利用接收的数据构造基于拉普拉斯核相关熵的低阶协方差矩阵R，其第行/>列元素/>表示为，/>式中/>和/>分别表示接收的第k次快拍信号数据矢量的第/>维和第/>维，/>(·)^*表示共轭运算，η为核函数的核长。

对于采用传统的加权信号子空间拟合方法进行波达方向估计，首先需要已知空间中的未知信源个数B，然后对由接收数据构造的低阶协方差矩阵R进行特征分解，其中U_S是由B个较大特征值对应的特征向量张成的信号子空间，Σ_S是由B个较大特征值构成的对角阵，U_N是由剩余较小特征值对应的特征向量张成的噪声子空间，Σ_N是由剩余较小特征值构成的对角阵，(·)^H表示共轭转置运算，均匀方阵的正交投影矩阵为/>(·)^-1表示求逆运算，由此再根据MIMO方阵的二维波达方向估计问题，设计出基于拉普拉斯核相关熵的加权信号子空间拟合方程通式为/> 为信号子空间的最优权矩阵，其中σ²表示噪声功率，tr(·)为矩阵求迹函数。

步骤三：利用分段式逐个求解思想，得到待搜索来波最终的目标函数。

本发明在加权信号子空间拟合方法的基础上，利用分段式逐个求解的思想，每次搜索只对一个未知信源进行求解，而其他已求得信源参数保持不变，并且不需要关于信源数的先验信息。

首先令空间中未知的第l个来波标号为l，初始时未知来波标号l=1，然后对加权信号子空间拟合方程的通式进行分解转换，那么对于转换后的第l个未知信源估计值的加权信号子空间拟合方程可以表示为其中/>分别表示已经搜索出的l-1维方位角参数矢量和俯仰角参数矢量，/>其中/>由前l个较大特征值对应的特征向量张成的信号子空间，/>是由前l个较大特征值构成的对角阵，/>是由剩余较小特征值对应的特征向量张成的噪声子空间，/>是由剩余较小特征值构成的对角阵，/>然后将拟合方程中的投影矩阵进行分解化简得到/>其中第一项与第l个信源的俯仰角和方位角估计无关，所以忽略，进而得到再次化简后的第l个信源的加权信号子空间拟合方程其中最后定义一个归一化矢量得到所以接下来对空间中假设存在的第l个未知来波采用量子冰晶优化机制搜索的加权信号子空间拟合目标方程表达为

步骤四：初始化估计第l个来波的量子冰晶机制的冰晶个体数及每个量子冰晶的量子位置和相应的映射态位置并获得全局最优量子位置。

设置量子冰晶的个体数为N_p，每个量子冰晶所具有的空间维数为最大的迭代次数为T_max，t代表迭代次数。对N_p个量子冰晶的量子位置在量子域[0,1]内随机初始化，则估计第l个信源时第t代第i个量子冰晶的量子位置定义为/>其中表示第t代第i个量子冰晶量子位置的第/>维，相应的映射态位置为其映射关系为/>其中 和/>分别表示第/>维变量的下边界和上边界，i=1,2,...,N_p，计算每个量子冰晶映射态位置的适应度值，在估计第l个信源时将第t代第i个量子冰晶映射态位置/>代入适应度函数表达式中，得到记录搜索到的指导第t代的量子冰晶中适应度值最大的量子位置记为全局最优量子冰晶位置/>其所对应的映射态位置记为/>

步骤五：初始化量子冰晶能量值，确定一个临时湖中心位置，与外界进行能量交换。

对于第t代第i个量子冰晶其初始能量为然后确定一个临时湖中心，为了避免冰晶在结冰过程中受到外部干扰，陷入局部最优，同时加快收敛速度，采用角度联合距离的方式，动态调整湖水中心位置。定义第t代第i个量子冰晶的多样性度量方程为表示第t代第i个量子冰晶/>与第j个晶体之间的的最小夹角，i=1,2,...,N_p，j=1,2,...,N_p且i≠j，/>越大表示该晶体与其他晶体分隔越明显，其多样性也就越好，其中/>||·||₂表示向量的二阶范式，(·)^T表示转置。然后定义第t代第i个量子冰晶角度联合距离动态调整方程为/>其中定义Z(ρ)为角度影响因子，τ为控制系数，/>为第t代第i个量子冰晶的能量值，最后获得临时的湖水中心位置为/>

自然界中湖水结冰是一个复杂的动态的能量变化过程，其中冰晶结冰主要受湖水中心、已结冰冰晶和大气三部分影响。为了简化模型，本搜索机制将这一过程总结为析出和冰冻两个阶段，其主要原因是由于在湖水结冰初期，湖中心对外散发的能量占主导地位，影响结冰进程，而到了结冰后期，已结成的冰晶壳与大气影响占据主导地位，吸收能量，因此在该搜索机制中，析出和冰冻两个阶段的区别可进一步简化表达为在不同阶段中其影响因素的能量占比不同。所以首先判断其迭代次数t是否小于等于T_max/3，若是，则量子冰晶进入析出阶段，更新每个晶体能量变化值；否则，量子冰晶进入冰冻阶段，更新每个晶体能量变化值。

对于第t代第i个量子冰晶，定义其能量变化方程为其中主要分为三部分，第一部分/>为第t代第i个量子冰晶从湖水中心获取到的能量，/>为第t代第i个量子冰晶到湖水中心的距离，q为一个能量定值；第二部分/>表示第t代第i个量子冰晶与已完成结冰的冰晶之间进行的能量交换，O为一个能量定值；第三部分F_i ^t=C表示第t代第i个量子冰晶被风带走的能量，C为一个能量定值，λ₁、λ₂和λ₃均为比例系数。通过调整λ₁、λ₂和λ₃三个比例系数来区别处于析出阶段或是冰冻阶段。

步骤六：更新量子冰晶能量，能量值最低的N_q个量子冰晶开始沉淀结冰，加入到已有的冰晶壳中，并更新历史量子位置空间

更新量子冰晶能量，第t+1代第i个量子冰晶的能量为能量更新后，能量值最低的N_q个量子冰晶沉淀结冰，同时记录这N_q个冰晶中适应度值最大的前/>个映射位置所对应的量子位置，组成历史量子位置空间/>其中/>其相应的映射态位置所代表的空间为其中/>随着迭代的进行，不断地更新替换历史量子位置空间/>中的量子位置。

步骤七：根据量子演化规则，使用模拟量子旋转门更新未沉淀的量子冰晶的量子位置和相应的映射态位置。

定义第t+1代第i个量子冰晶的第维量子旋转角为/> 其中为转移因子，/>为第t+1代第i个晶体的能量值，μ为比例系数，/>为在[0,1]之间服从均匀分布的随机数，/>ν^t为第t代惯性权重系数，ν_max为ν^t的最大值，v_min为v^t的最小值，/>用于查找第t代/>中最接近/>的量子位置的第/>维。最后对每个量子冰晶的量子位置通过简化的模拟量子旋转门进行更新，则第t+1代第i个量子冰晶的第/>维量子位置更新为同时通过映射方程获得第t+1代第i个量子冰晶的第/>维映射态位置为/>

步骤八：计算每个量子冰晶的适应度值，根据轮盘赌选择策略选择量子冰晶来产生新一代量子冰晶的最终位置，然后更新全局最优量子位置。

通过轮盘赌选择策略，以一定的概率从完成模拟量子旋转门更新后的量子冰晶中选择若干个个体，来指导生成新的量子冰晶。因此首先计算第i个量子冰晶个体被选中的概率对于第t+1代第i个量子冰晶，其被选中的概率为/> 表示第t+1代第i个量子冰晶的适应度值，然后确定第i个个体的累计概率/>对于第i个冰晶，其累积概率/>接着生成一个[0,1]之间服从均匀分布的随机数/>若/>大于等于累积概率/>而小于/>则第i个量子冰晶被选中来指导生成一个新的量子冰晶，其第维量子位置生成公式为/> 需要约束在[0,1]之间，超过边界值，取到边界，/>为权重因子，/>为在[0,1]之间服从均匀分布的随机数，然后计算其相应的映射态位置/>以此重复生成N_q个量子冰晶，补充更新过程中沉淀结冰的量子冰晶，保证量子冰晶个体数不变。最后将生成冰晶的映射态位置带入到适应度值公式中，计算对应的适应度值，并与原有更新后的量子冰晶的适应度值放在一起，选择其中的最大值与上一代的全局最优量子位置的适应度值进行比较，最终选择更大适应度值所对应的量子位置作为指导第t+1代量子冰晶的全局最优量子位置，记为/>其相应的映射态位置记为/>

步骤九：判断是否达到其最大迭代次数T_max。若未达到，令t=t+1，返回步骤五继续迭代；若达到最大迭代次数，则选择当前的全局最优量子位置的映射态位置作为最后的结果，并输出此时的映射态位置

步骤十：判断信源l是否存在。量子冰晶优化机制属于最大值搜索并且根据搜索的目标函数可知，若存在未知信源，那么输出映射态位置的适应度值相较于迭代过程中的其他位置应处于绝对优势地位，因此通过设置一个阈值判断输出映射态位置的适应度值是否远远优于其他位置，进而判断是否为未知信源。若/>则认为输出的映射态位置的未达到远远优于的状态，因此判断空间中已无未知信源，然后输出空间中的已测得的未知信源的个数以及相应的二维波达方向的估计值，信源数为l-1。若/>则选择当前的映射态位置/>作为该信源二维波达方向的估计值，再令l=l+1，返回步骤三，并将得到的波达方向估计值作为求得的先验信息带入下一次求解的目标函数化简中，然后再一次搜索，直至输出的适应度值小于阈值，空间中无未知信源为止。

对于阈值的设置，通过计算迭代过程中所有映射位置的平均适应度δ，确定一个倍数/>设置/>同时根据实际问题情况确定一个阈值下限/>来保证尽可能消除搜索机制的偶然性或者误判情况，借此可以通过该阈值判断出输出位置的适应度值是否远远优于其他位置，进而判断输出的映射态位置是否为一个未知信源的二维波达方向估计值。

为便于叙述，将基于量子冰晶机制和拉普拉斯核函数相关熵的加权信号子空间拟合方法简记为QCEO-LLOC-WSSF，对比基于无穷范数归一化的MUSIC方阵波达方向估计方法记为IN-MUSIC。其中仿真实验参数设计如下：对于Massive MIMO方阵，其x轴方向阵元数M=16，y轴方向阵元数N=16，x方向上和y方向上阵元间距均为d＝λ/2，拉普拉斯核函数的核长η＝129，最大迭代次数设置为T_max=150， 量子冰晶个体总数N_p=80，空间维数/>控制系数τ＝2.5，冰晶壳中量子位置个数/>能量定值q=1，O=0.6，C=0.6，对于比例系数λ₁、λ₂和λ₃，当算法初期，处于析出阶段时λ₁=1，λ₂=0，λ₃＝-0.5，当处于沉淀结冰阶段时λ₁=0.1，λ₂＝-0.85，λ₃＝-0.5，沉淀个数N_q=20，权重因子比例系数μ＝0.005，倍数/> ν_max=2，ν_min=0，冲击噪声设置α＝1.5。对于IN-MUSIC方法，采用扫描间隔为1°的谱峰搜索方式。

图2中，设置信源个数B=5，分别从{60°,20°}、{100°,50°}、{160°,30°}、{160°,50°}、{300°,75°}方向入射，在冲击噪声的背景下、广义信噪比为10dB、快拍数K=20时，信源个数与波达方向联合估计示意图。从仿真图2中可以看出，本发明可以准确的测得信源个数，精度较高，测角范围全覆盖。

图3中，设置信源个数B=2，两个信源分别从{60°,20°}和{250°,70°}方向入射，冲击噪声特征指数α＝1.5，广义信噪比为10dB，快拍数K=20，Monte Carlo实验次数为30次。从仿真图3中可以看出，测向效果比较理想，估计值与真实值之间偏差较小，只有极个别点存在一定的偏差。因此可以看出所设计的QCEO-LLOC-WSSF可以在冲击噪声环境下进行有效的DOA估计。

图4中，设置信源个数B=2，两个信源分别从{250°,70°}和{100.5°,30°}方向入射，采用本发明中的QCEO-LLOC-WSSF与经典的IN-MUSIC方法同时进行处理，比较快拍数对两种方法测向效果的影响。冲击噪声特征指数为1.5，广义信噪比为10dB，Monte Carlo试验次数30次。从仿真图4的实验结果可以看出，在冲击噪声且快拍数较小的背景下，所设计的QCEO-LLOC-WSSF方法具有更大的优越性，可以在恶劣的环境下较好的完成波达方向估计，精度更高。而反观IN-MUSIC算法，则无法适应此种环境。所以也展示了所发明的QCEO-LLOC-WSSF方法在小快拍和冲击噪声等较差条件下依旧可以进行有效的DOA估计。

图5中，设置信源个数B=2，两个信源分别从{250°,70°}和{100.5°,30°}方向入射，采用本发明中的QCEO-LLOC-WSSF与经典的IN-MUSIC方法同时进行处理，比较广义信噪比对两种方法测向效果的影响。冲击噪声特征指数α＝1.5，快拍数K=20，Monte Carlo试验次数30次。从仿真图5的实验结果可以看出，所设计的QCEO-LLOC-WSSF方法还是具有很好的优越性，能够看出在小信噪比时，在通信质量较差的情况下依旧可以完成信源数与波达方向的联合估计。

Claims (10)

1.一种大规模MIMO方阵信源数与波达方向联合估计方法，其特征在于，包括：

步骤一：建立冲击噪声环境下的Massive MIMO方阵模型，获取阵列接收的快拍数据；

步骤二：利用接收数据构造基于拉普拉斯核相关熵的低阶协方差矩阵，得到基于拉普拉斯核相关熵的通用加权信号子空间拟合方程；

步骤三：利用分段式逐个求解方法，得到待搜索来波最终的目标函数；

步骤四：初始化估计第l个来波的量子冰晶机制的冰晶个体数及每个量子冰晶的量子位置和相应的映射态位置并获得全局最优量子位置；

步骤五：初始化量子冰晶能量值，确定一个临时湖中心位置，与外界进行能量交换；

步骤六：更新量子冰晶能量，能量值最低的N_q个量子冰晶开始沉淀结冰，加入到已有的冰晶壳中，并更新历史量子位置空间

步骤七：根据量子演化规则，使用模拟量子旋转门更新未沉淀的量子冰晶的量子位置和相应的映射态位置；

步骤八：计算每个量子冰晶的适应度值，根据轮盘赌选择策略选择量子冰晶来产生新一代量子冰晶的最终位置，然后更新全局最优量子位置；

步骤九：判断是否达到其最大迭代次数T_max，若未达到，令t＝t+1，返回步骤五继续迭代；若达到最大迭代次数，则选择当前的全局最优量子位置的映射态位置作为最后的结果，并输出此时的映射态位置

步骤十：判断信源l是否存在：设置阈值若/>则判断空间中已无未知信源，然后输出空间中的已测得的未知信源的个数以及相应的二维波达方向的估计值，信源数为l-1；若/>则选择当前的映射态位置/>作为该信源二维波达方向的估计值，再令l＝l+1，返回步骤三，并将得到的波达方向估计值作为求得的先验信息带入下一次求解的目标函数化简中，然后再一次搜索。

2.根据权利要求1所述的一种大规模MIMO方阵信源数与波达方向联合估计方法，其特征在于：步骤一所述阵列接收的快拍数据满足：

假设空间中存在B个窄带远场信源，分别以方位角θ_b和俯仰角方向入射到由M×N个阵元构成的Massive MIMO水平均匀方阵上，b＝1,2,…,B，阵元间距为d，入射波长为λ，则阵列接收的第k次快拍采样数据的数学模型为：

其中k＝1,2,...,K，K为最大快拍数，z(k)＝[z₁(k),z₂(k),...,z_MN(k)]^T为MN×1维的阵列接收的第k次快拍数据矢量，上标T表示转置，为MN×B维的流型矩阵，其中θ＝[θ₁,θ₂,...,θ_B]和/>分别为信源的方位角矢量与俯仰角矢量，/>表示第b个信源的导向矢量，其中为流型矩阵关于x轴方向上的第b个信源的导向矢量，/>为流型矩阵关于y轴方向上的第b个信源的导向矢量，/>为Kronecker积，s(k)为B×1维信号矢量，n(k)为MN×1维服从SαS稳定分布的冲击噪声矢量。

3.根据权利要求2所述的一种大规模MIMO方阵信源数与波达方向联合估计方法，其特征在于：步骤二所述基于拉普拉斯核相关熵的低阶协方差矩阵和基于拉普拉斯核相关熵的通用加权信号子空间拟合方程具体为：

基于拉普拉斯核相关熵的低阶协方差矩阵R的第行/>列元素/>表示为：

式中，z_i(k)和z_j(k)分别表示接收的第k次快拍信号数据矢量的第维和第/>维，(·)^*表示共轭运算，η为核函数的核长；

已知空间中的未知信源个数B，然后对低阶协方差矩阵R进行特征分解：

其中，U_S是由B个较大特征值对应的特征向量张成的信号子空间，Σ_S是由B个较大特征值构成的对角阵，U_N是由剩余较小特征值对应的特征向量张成的噪声子空间，Σ_N是由剩余较小特征值构成的对角阵，(·)^H表示共轭转置运算，均匀方阵的正交投影矩阵为(·)^-1表示求逆运算；

则基于拉普拉斯核相关熵的加权信号子空间拟合方程通式为：

其中，为信号子空间的最优权矩阵，其中σ²表示噪声功率，tr(·)为矩阵求迹函数。

4.根据权利要求3所述的一种大规模MIMO方阵信源数与波达方向联合估计方法，其特征在于：步骤三所述待搜索来波最终的目标函数包括：

其中，空间中未知的第l个来波标号为l，初始时未知来波标号l＝1，为定义的归一化矢量，满足：

其中， 和/>分别表示已经搜索出的l-1维方位角参数矢量和俯仰角参数矢量，低阶协方差矩阵R分解为：其中/>由前l个较大特征值对应的特征向量张成的信号子空间，/> 是由前l个较大特征值构成的对角阵，σ²表示噪声功率。

5.根据权利要求4所述的一种大规模MIMO方阵信源数与波达方向联合估计方法，其特征在于：步骤四所述初始化估计第l个来波的量子冰晶机制的冰晶个体数及每个量子冰晶的量子位置和相应的映射态位置并获得全局最优量子位置包括：

设置量子冰晶的个体数为N_p，每个量子冰晶所具有的空间维数为最大的迭代次数为T_max，t代表迭代次数，对N_p个量子冰晶的量子位置在量子域[0,1]内随机初始化，则估计第l个信源时第t代第i个量子冰晶的量子位置定义为/>其中/>表示第t代第i个量子冰晶量子位置的第/>维，相应的映射态位置为其映射关系为/>其中和/>分别表示第/>维变量的下边界和上边界，i＝1,2,...,N_p，计算每个量子冰晶映射态位置的适应度值，在估计第l个信源时将第t代第i个量子冰晶映射态位置/>代入适应度函数表达式中，得到记录搜索到的指导第t代的量子冰晶中适应度值最大的量子位置记为全局最优量子冰晶位置/>其所对应的映射态位置记为/>

6.根据权利要求5所述的一种大规模MIMO方阵信源数与波达方向联合估计方法，其特征在于：步骤五所述初始化量子冰晶能量值，确定一个临时湖中心位置，与外界进行能量交换包括：

对于第t代第i个量子冰晶其初始能量为然后确定一个临时湖中心，采用角度联合距离的方式，动态调整湖水中心位置，定义第t代第i个量子冰晶的多样性度量方程为/>表示第t代第i个量子冰晶/>与第j个晶体之间的的最小夹角，i＝1,2,...,N_p，j＝1,2,...,N_p且i≠j，/>越大表示该晶体与其他晶体分隔越明显，其多样性也就越好，其中/>||·||₂表示向量的二阶范式，(·)^T表示转置；然后定义第t代第i个量子冰晶角度联合距离动态调整方程为/>其中定义Z(ρ)为角度影响因子，τ为控制系数，/>为第t代第i个量子冰晶的能量值，最后获得临时的湖水中心位置为/>

将湖水结冰过程总结为析出和冰冻两个阶段，首先判断其迭代次数t是否小于等于T_max/3，若是，则量子冰晶进入析出阶段，更新每个晶体能量变化值；否则，量子冰晶进入冰冻阶段，更新每个晶体能量变化值；

对于第t代第i个量子冰晶，定义其能量变化方程为其中主要分为三部分，第一部分/>为第t代第i个量子冰晶从湖水中心获取到的能量，/>为第t代第i个量子冰晶到湖水中心的距离，q为一个能量定值；第二部分/>表示第t代第i个量子冰晶与已完成结冰的冰晶之间进行的能量交换，O为一个能量定值；第三部分F_i ^t＝C表示第t代第i个量子冰晶被风带走的能量，C为一个能量定值，λ₁、λ₂和λ₃均为比例系数。通过调整λ₁、λ₂和λ₃三个比例系数来区别处于析出阶段或是冰冻阶段。

7.根据权利要求6所述的一种大规模MIMO方阵信源数与波达方向联合估计方法，其特征在于：步骤六所述更新量子冰晶能量，能量值最低的N_q个量子冰晶开始沉淀结冰，加入到已有的冰晶壳中，并更新历史量子位置空间包括：

更新量子冰晶能量，第t+1代第i个量子冰晶的能量为能量更新后，能量值最低的N_q个量子冰晶沉淀结冰，同时记录这N_q个冰晶中适应度值最大的前/>个映射位置所对应的量子位置，组成历史量子位置空间/> 其中其相应的映射态位置所代表的空间为其中/>随着迭代的进行，不断地更新替换历史量子位置空间/>中的量子位置。

8.根据权利要求7所述的一种大规模MIMO方阵信源数与波达方向联合估计方法，其特征在于：步骤七所述根据量子演化规则，使用模拟量子旋转门更新未沉淀的量子冰晶的量子位置和相应的映射态位置包括：

定义第t+1代第i个量子冰晶的第维量子旋转角为/> 其中为转移因子，/>为第t+1代第i个晶体的能量值，μ为比例系数，/>为在[0,1]之间服从均匀分布的随机数，/>ν^t为第t代惯性权重系数，ν_max为ν^t的最大值，ν_min为ν^t的最小值，/>用于查找第t代/>中最接近/>的量子位置的第/>维，最后对每个量子冰晶的量子位置通过简化的模拟量子旋转门进行更新，则第t+1代第i个量子冰晶的第/>维量子位置更新为同时通过映射方程获得第t+1代第i个量子冰晶的第/>维映射态位置为/>

9.根据权利要求8所述的一种大规模MIMO方阵信源数与波达方向联合估计方法，其特征在于：步骤八：计算每个量子冰晶的适应度值，根据轮盘赌选择策略选择量子冰晶来产生新一代量子冰晶的最终位置，然后更新全局最优量子位置包括：

首先计算第i个量子冰晶个体被选中的概率对于第t+1代第i个量子冰晶，其被选中的概率为/> 表示第t+1代第i个量子冰晶的适应度值，然后确定第i个个体的累计概率/>对于第i个冰晶，其累积概率/>接着生成一个[0,1]之间服从均匀分布的随机数/>若/>大于等于累积概率/>而小于/>则第i个量子冰晶被选中来指导生成一个新的量子冰晶，其第/>维量子位置生成公式为 需要约束在[0,1]之间，超过边界值，取到边界，/>为权重因子，/>为在[0,1]之间服从均匀分布的随机数，然后计算其相应的映射态位置以此重复生成N_q个量子冰晶，补充更新过程中沉淀结冰的量子冰晶，保证量子冰晶个体数不变；最后将生成冰晶的映射态位置带入到适应度值公式中，计算对应的适应度值，并与原有更新后的量子冰晶的适应度值放在一起，选择其中的最大值与上一代的全局最优量子位置的适应度值进行比较，最终选择更大适应度值所对应的量子位置作为指导第t+1代量子冰晶的全局最优量子位置，记为/>其相应的映射态位置记为/>

10.根据权利要求1所述的一种大规模MIMO方阵信源数与波达方向联合估计方法，其特征在于：所述阈值的设置包括：

设置其中δ是通过计算迭代过程中所有映射位置的平均适应度，/>是设定的倍数，同时确定一个阈值下限/>取/>和/>中的较大值作为阈值。