CN114019448B - 毫米波大规模多输入多输出系统二维波达角快速估计方法 - Google Patents

Abstract

本发明提供了毫米波大规模多输入多输出系统二维波达角快速估计方法，包括如下步骤：步骤1，对多输入多输出基站的接收阵列的信号进行空间压缩采样；步骤2，进行协方差分解获得信号子空间；步骤3，进行二维波达角估计。

Description

毫米波大规模多输入多输出系统二维波达角快速估计方法

技术领域

本发明涉及毫米波大规模多输入多输出系统二维波达角快速估计方法。

背景技术

传感器阵列是5G/6G网络中最重要的基础设施之一，是传输、定位和感知的基本组成部分。通过传感器阵列进行波达角(DOA)估计是一个众所周知的非线性问题。目前，已经有许多算法被用于解决这个问题。例如，旋转时不变技术(ESPRIT)、最大似然估计(ML)、多重信号分类(MUSIC)。通常，谱搜索方法，如MUSIC和ML，往往效率低下。此外，他们很难避免Off-grid问题(Z.Yang,L.Xie and C.Zhang,"Off-grid direction of arrivalestimation using sparse Bayesian inference,"IEEE Trans.Signal Process.,vol.61,no.1,pp.38-43,Jan.1,2013.)。ESPRIT的高效率已被普遍认可，因为它可以提供参数估计的闭式解。

现有的绝大多数研究主要集中在使用标量传感器阵列对进行一维(1D)DOA估计。在实际工程中，二维(2D)DOA更具吸引力。要提供2D-DOA估计，传统标量传感器的阵列几何结构必须严格非线性，例如L形、圆形、矩形。但是，对于空间有限的平台来说，标量传感器阵列可能无法实现。此外，在参数估计时必须提供传感器的位置的先验信息，否则需要进行复杂的差误补偿。使用电磁矢量传感器(EMVS)可以避免上述缺点(A.Nehorai,E.Paldi,"Vector-sensor array processing for electromagnetic source localization,"IEEETrans.Signal Process.,vol.42,no.2,pp.376-398,Feb.1994.)一个完整的EMVS不仅能提供入射信号的2D-DOA估计，还能提供额外的极化状态估计，这可能有助于探测弱信源。此外，单个EMVS具有6个自由度(DOF)，远高于传统的标量传感器，因而具有更好的可辨识性。

与标量传感器阵列中的估计任务不同，使用EMVS阵列的测算问题往往更为复杂，因为它涉及2D-DOA(方位角和俯仰角)和二维极化状态(极化相差和辅助极化角)。在文献K.T.Wong,M.D.Zoltowski,“Closed-form direction finding and polarizationestimation with arbitrarily spaced electromagnetic vector-sensors at unknownlocations,”IEEE Trans.Antennas Propagat.,vol.48,no.5,pp.671-681,May 2000中，引入了ESPRIT算法，其首先利用阵列的旋转/位移不变特性得到归一化坡印亭矢量，再使用矢量叉乘技术获得2D-DOA；文献M.D.Zoltowski,K.T.Wong,"ESPRIT-based 2-D directionfinding with a sparse uniform array of electromagnetic vector sensors,"IEEETrans.Signal Process.,vol.48,no.8,pp.2195-2204,Aug.2000研究了基于稀疏均匀矩形阵列(URA)的2D-DOA估计框架。它首先首先利用阵列的旋转/位移不变特性得到模糊的坡印亭向量，然后利用阵列几何结构的稀疏性实现无模糊的角度估计。在文献K.T.Wong,X.Yuan,"“Vector cross-product direction-finding”with an electromagneticvector-sensor of six orthogonally oriented but spatially noncollocatingdipoles/loops,"IEEE Trans.Signal Process.,vol.59,no.1,pp.160-171,Jan.2011中，作者将矢量叉积方法推广到具有分布偶极/环的EMVS阵列；为了减少子空间方法中的计算负担，J.He,Z.Liu,“Computationally efficient 2D direction finding andpolarization estimation with arbitrarily spaced electromagnetic vectorsensors at unknown locations using the propagator method,”Digital SignalProcess.vol.19,no.3,pp.491-503,May 2009中讨论了传播子估计器，它不需要进行复杂的特征分解计算。在T.Ahmed,Z.Xiaofei,Z.Wang,“DOA estimation for coprime EMVSarrays via minimum distance criterion based on PARAFAC analysis,”IET RadarSonar Nav.,vol.13,no.1,pp.65–73,Jan.2019钟，作者提出了基于互质阵列结构的EMVS流形，并设计了一个基于张量分解的估计器来以利用EMVS阵列的多维特性。最近，在L.Wan,K.Liu,Y.C.Liang,T.Zhu,“DOA and polarization estimation for non-circularsignals in 3-D millimeter wave polarized massive MIMO systems,”IEEEWirel.Commun.,vol.20,no.5,pp.3152-3167,May 2021，作者提出了毫米波极化大规模MIMO系统的概念，并给出了一种类似MUSIC算法的方法。不幸的是，它只适合两个偶极子传感器配置的URA几何形状。与传统对应方一样，它对传感器位置误差敏感。此外，由于需要进行谱峰搜索，它的计算效率很低，并不适合大规模MIMO系统。

发明内容

发明目的：为解决背景技术中存在的技术问题，本发明提出一种毫米波大规模MIMO系统二维波达角快速估计方法，采用共点式EMVS作为MIMO系统的收发天线，并采用压缩感知技术对数据进行压缩。该框架结构灵活，适用于任意几何形状的阵列结构，且对传感器位置误差不敏感。本发明方法包括如下步骤：

步骤1，对MIMO基站的接收信号进行空间压缩采样；

步骤2，对样本协方差矩阵进行分解获得信号子空间；

步骤3，进行二维波达角估计。

步骤1包括：对MIMO基站的接收信号进行空间压缩采样。

步骤1还包括：用N表示MIMO基站的共址电磁矢量传感器EMVS接收阵元，N为整数；用K表示与MIMO基站上行链路中的远场用户源信号的个数，K为整数；第n个电磁矢量传感器EMVS的坐标r_n为其中，n＝1,2,…,N，(·)^T表示转置运算。

步骤1还包括：令θ_k，φ_k，γ_k，η_k分别表示需要估计的k个用户源信号的俯仰角、方位角、辅助极化角和极化相位差，含噪的阵列响应y(t)为：

其中，k＝1,2,…,K，t表示快拍索引；和⊙分别表示克罗内克积和KhatriRao积：表示第k个源信号对应的导向矢量，表示以自然指数为底、以-j2πτ_1,k/λ为指数的数，表示复数标识符，g_k为中间变量，具体为λ是载波波长；为第k个源信号对应的极化响应矢量，b_k(q)表示b_k的第q个元素，q＝1,2,…,6；s_k(t)表示第k个源信号，n(t)表示阵列噪声，表示阵列导向矩阵，表示阵列的极化响应矩阵，如果阵列噪声n(t)满足均值为零，方差为σ²的高斯分布，且源信号s(t)与噪声是不相关的。

步骤1还包括：定义空间下采样算子空间下采样算子将N维矢量线性映射为N₁维矢量，N₁<N，即：

其中，为压缩空间的响应矢量；具体地，对高维矢量的压缩过程通过一个降维矩阵实现：

是一个降维矩阵，M<N；定义同步压缩矩阵F：

其中，对MIMO基站的阵列接收数据进行同步压缩后得到：

其中，I₆表示6×6的单位矩阵，为中间变量，

步骤1中，设定EMVS阵列具有任意的几何形状。

步骤2包括：对权利要求2所获得的样本x(t)的协方差矩阵进行分解，获得信号子空间；

为高斯过程n(t)通过同步压缩矩阵F的结果，它也是一个高斯的过程，其均值为零，协方差矩阵为σ²(FF^H)，上标(·)^H表示共轭转置运算；设定D行是相互正交的，则存在则得到：

其中，E{·}表示数学期望；令所以x(t)的协方差R_x为：

当存在L个观测样本时，t＝1,2,…,L，R_x通过下式进行估计：

其中，表示R_x的估计值；对进行特征分解，获得信号子空间E_s，信号子空间E_s由K个最大的特征值所对应的特征向量组成；且存在一个非奇异的矩阵满足：

步骤3包括：设计估计器，估计器的二维波达角估计基于以下旋转不变性原理：

其中，以对角矩阵的形式返回B的第q行，即diag{·}表示对角化运算，q＝1,2,…,6，/表示除法。

步骤3还包括：定义选择性矩阵J_q如下：

其中，i_6,q表示为I₆的第q行，得到：

即存在公式(10)被重写为：

步骤3还包括：将公式(9)的结果代入公式(13)，得到：

其中，上标(·)^-1表示逆；等价地，公式(14)表述为：

其中，上标表示伪逆；从的特征分解，能够得到对Φ^(1,2)的估计和T的估计，T的估计为用替换T，并计算公式(15)除第一行除外的左侧部分，得出Φ^(1,3)、Φ^(1,4)、Φ^(1,5)、Φ^(1,6)的估计，从Φ^(1,3)、Φ^(1,4)、Φ^(1,5)、Φ^(1,6)的对角线元素获得对的估计

将b_k重新表示为：

得到：

其中，*表示矢量交叉乘，(·)^*表示共轭；||·||表示向量的F-范数；由于||b_k(1)||²是常数，则：

将的估计值带入(18)中，获得u_k、v_k、w_k的估计值，并分别用表示，二维波达角估计通过下式计算得到：

其中，和分别表示θ_k的估计值和φ_k的估计值。

本发明具有如下有益效果：1.本发明适用于任意的阵列流形；

2.本发明与传统的ESPRIT算法性能相差并不大；

3.本发明方法具有更快的运算速度，适合于大规模阵列天线的场景；

4.本发明方法对天线的位置误差不敏感。

附图说明

下面结合附图和具体实施方式对本发明做更进一步的具体说明，本发明的上述和/或其他方面的优点将会变得更加清楚。

图1是本发明方法的架构。

图2是本发明所提估计器的2D-DOA散点图。

图3是2D-DOA估计随SNR变化的的平均RMSE比较示意图。

图4是2D-DOA估计随M变化的平均RMSE比较示意图。

图5是不同方案平均运行时间对M变化的比较示意图。

具体实施方式

一个完整的共置电磁矢量传感器EMVS包括六个传感器：三个电偶极子(测量电场的强度)和三个磁环(测量磁场的强度)。考虑一个远场源信号入射到一个完备的电磁矢量传感器EMVS阵元上的场景。用θ、φ、γ、η分别表示源信号相对于电磁矢量传感器EMVS阵元的俯仰角、方位角、辅助极化角和极化相位差。电磁矢量传感器EMVS的极化方向矢量由下式给出：

其中表示‘定义’，上标(·)^T表示转置运算，e和m分别表示电场响应矢量和磁场响应矢量，D、p分别为：

此外，e和m之间的矢量叉乘满足：

其中*表示矢量交叉乘，(·)^*表示共轭，||·||表示向量的F-范数。

在本发明中，考虑一个大规模MIMO(multiple-in multipleout；MIMO，多进多出)基站(BS)的上行链路，设定基站BS包括N个共址电磁矢量传感器EMVS阵元。设定K个远场用户源信号通过上行链路同基站连接在一起。特别地，设定电磁矢量传感器EMVS具有任意的几何形状。设定第n个电磁矢量传感器EMVS的坐标为r_n＝[x_n,y_n,z_n]^T，令(θ_k，φ_k，γ_k，η_k)为需要估计的k个(k＝1,2,…,K)用户源信号的的角度参数。含噪的阵列响应则为：

其中，t表示快拍索引；和⊙分别表示克罗内克积和KhatriRao积(按列克罗内克积)：表示第k个源信号对应的导向矢量， λ是载波长；为第k个源信号对应的极化响应矢量。s_k(t)表示第k个源信号，表示阵列噪声，在进一步推导之前，将确定以下假设：

假设1：阵列噪声n(t)假定满足均值为零，方差为σ²的高斯分布；

假设2：源信号与噪声是不相关的；

根据假设1和假设2，数据矢量y(t)的协方差矩阵表述为：

在这里上标(·)^H表示共轭转置运算，E{·}表示数学期望；I_6N表示6N×6N的单位矩阵。当存在L个观测样本时，数据矩阵写成：

其中，实际工程中，R_y通过下式进行估计：

传统的ESPRIT算法是基于奇异值分解或特征值分解在大规模MIMO配置的情况下计算效率低下。在本发明中，提出了一个基于空间压缩采样的框架，如图1所示，具体描述如下。

空间压缩采样：

定义空间下采样算子该算子将N维矢量线性映射为N₁维(N₁<N)矢量，即：

其中，为压缩空间的响应矢量。应当强调，这种压缩是线性的。实际上，这种下采样算子可以由一个降维矩阵替代。对高维矢量的压缩过程被重写为：

通常，D称为测量(观测)矩阵。一般来说，高斯矩阵、伯努利矩阵被经常选为测量矩阵。

在本发明中，每个电磁矢量传感器EMVS的六个分量都是在空域同步压缩。同步压缩矩阵的定义为：

根据性质有：

其中，实际上，解释为高斯过程n(t)通过线性系统F的结果。因此，也是一个高斯的过程，其均值为零，协方差矩阵为σ²(FF^H)。设定D行是相互正交的，则存在因此，故：

所以x(t)的协方差为：

与R_y相比，R_x的维度要小得多。因此，对R_x的矩阵分解会大大降低计算复杂度。同样，R_x通过下式进行估计：

对进行特征分解，获得信号子空间E_s，它由K个最大的特征值所对应的特征向量组成。显然，E_s与张成相同的子空间。即存在一个非奇异的矩阵满足：

2D-DOA估计：

大多数现有基于压缩采样的估计器都采用凸或非凸优化工具求解稀疏解，尽管它们通常提供了超分辨率的估计性能，但由于涉及网格搜索，它们的效率往往低下。本发明提出了一种类似ESPRIT的估计器，它提供了参数估计的封闭解。该估计器的2D-DOA估计基于以下旋转不变性原理：

其中，diag{·}表示对角化运算，q＝1,2,…,6，定义：

其中，i_6,q表示为I₆的第q行。再次使用性质得到：

即存在所以，公式(16)被重写为：

将公式(15)的结果代入公式(19)，得到：

其中，上标(·)^-1表示逆。等价地，公式(20)表述为：

其中，上标表示伪逆。可以看出，从的特征分解，可以得到对Φ^(1,2)的估计和T的估计(由表示)。用替换T，并计算公式(21)(第一行除外)的左侧部分，可以得出Φ^(1,3)Φ^(1,4)Φ^(1,5)Φ^(1,6)的估计，从它们的对角线元素可以获得对(由表示)的估计。

根据公式(1)，b_k重新表示为：

易知：

由于||b_k(1)||²是常数，则：

令u_k、v_k、w_k的估计分别为2D-DOA估计通过下式计算得到：

实施例

为了验证本发明框架的有效性，采用蒙特卡罗方法来评估估算性能。在此，设定MIMO基站由一个M行M列的URA组成，每个传感器都是共址式EMVS。假设K＝3个远场源信号，其参数分别为θ＝(10°,50°,30°)，φ＝(-20°,25°,55°)，γ＝(10°,30°,60°)，η＝(20°,40°,55°)。此外，假设已收集L个样本。仿真的每个图的结果都依赖于200次独立的试验。在仿真中，信号与噪声比(SNR)被定义为SNR＝10lg(||Y-N||²/||N||²)。性能评估采用两种方法：均方根误差(RMSE)和平均运行时间。在本发明的框架中，定义压缩因子α为，N₁＝round(αM²)，D是随机生成的随机矩阵。

首先，所提估计器的2D-DOA估计的散点图结果由图2给出，其中M＝10、α＝0.5、L＝500和SNR＝10dB。可以清楚地看到，所有的参数都正确估计和自动配对。结果表明，本框架方案是适用于大规模MIMO系统的。

其次，图3中给出了2D-DOA估计在不同的SNR时的RMSE性能，其中M＝10和L＝500。值得注意的是，当SNR增加时，所有算法都会提供更好的RMSE性能。然而，本发明所提的方法的RMSE比ESPRIT更高，这表明ESPRIT提供更精确的估计性能。这是由于本发明的估计器中的压缩过程会降低阵列DOF，从而导致估计精度降低。此外，图3表明，较小的α导致更高的RMSE，因为它意味着牺牲了更多的DOF。

最后，测试了在不同的基站规模M条件下参数估计的性能比较，其中SNR＝0dB和L＝500。图4显示了平均RMSE比较，图5给出了运行时间(运行时间即纵坐标)的比较。很明显，RMSE在M增加时都有所改进。与之前的观测结果类似，ESPRIT的RMSE性能优于本发明所提的估算器，而较小的α会导致2D-DOA估计精度更加糟糕。另一方面，附图5表明，本发明所提的框架需要的运行时间比ESPRIT少，特别是在M相对较大时，这种优势更为明显。应该强调的是，在大规模MIMO的配置下(例如，M>16)，ESPRIT和本发明所提算法之间的RMSE距离小于一个数量级，并且随着M的增加，差距会缩小。但是，它们之间的运行时间距离一般大于一个数量级(α＝0.8的情况除外)，并且随着M的增加，差距会增加。因此，所提的估计器在大规模MIMO配置方面更具吸引力，因为它提供的估计精度与ESPRIT相当，但其计算负担比后者低得多。实际工程中，α应该被合理地选择以平衡估计的准确性和计算的复杂性。

本发明提供了一种毫米波大规模MIMO系统二维波达角快速估计方法，具体实现该技术方案的方法和途径很多，以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。本实施例中未明确的各组成部分均可用现有技术加以实现。

Claims (4)

1.毫米波大规模多输入多输出系统二维波达角快速估计方法，其特征在于，包括如下步骤：

步骤1，对MIMO基站的接收信号进行空间压缩采样；

步骤2，对样本协方差矩阵进行分解获得信号子空间；

步骤3，进行二维波达角估计；

步骤1还包括：用N表示MIMO基站的共址电磁矢量传感器EMVS接收阵元，N为整数；用K表示与MIMO基站上行链路中的远场用户源信号的个数，K为整数；第n个电磁矢量传感器EMVS的坐标r_n为其中，n＝1,2,…,N，(·)^T表示转置运算；

步骤1还包括：令θ_k，φ_k，γ_k，η_k分别表示需要估计的k个用户源信号的俯仰角、方位角、辅助极化角和极化相位差，含噪的阵列响应y(t)为：

其中，k＝1,2,…,K，t表示快拍索引；和⊙分别表示克罗内克积和KhatriRao积：表示第k个源信号对应的导向矢量，表示以自然指数为底、以-j2πτ_1,k/λ为指数的数，表示复数标识符，g_k为中间变量，具体为λ是载波波长；为第k个源信号对应的极化响应矢量，b_k(q)表示b_k的第q个元素，q＝1,2,…,6；s_k(t)表示第k个源信号，n(t)表示阵列噪声，表示阵列导向矩阵，表示阵列的极化响应矩阵，如果阵列噪声n(t)满足均值为零，方差为σ²的高斯分布，且源信号s(t)与噪声是不相关的；

步骤1还包括：定义空间下采样算子空间下采样算子将N维矢量线性映射为N₁维矢量，N₁<N，即：

其中，为压缩空间的响应矢量；具体地，对高维矢量的压缩过程通过一个降维矩阵实现：

是一个降维矩阵，M<N；定义同步压缩矩阵F：

其中，对MIMO基站的阵列接收数据进行同步压缩后得到：

其中，I₆表示6×6的单位矩阵，为中间变量，

步骤2包括：为高斯过程n(t)通过同步压缩矩阵F的结果，它也是一个高斯的过程，其均值为零，协方差矩阵为σ²(FF^H)，上标(·)^H表示共轭转置运算；设定D行是相互正交的，则存在则得到：

其中，E{·}表示数学期望；令所以x(t)的协方差R_x为：

当存在L个观测样本时，t＝1,2,…,L，R_x通过下式进行估计：

其中，表示R_x的估计值；对进行特征分解，获得信号子空间E_s，信号子空间E_s由K个最大的特征值所对应的特征向量组成；且存在一个非奇异的矩阵满足：

步骤3包括：设计估计器，估计器的二维波达角估计基于以下旋转不变性原理：

其中，以对角矩阵的形式返回B的第q行，即diag{·}表示对角化运算，q＝1,2,…,6，/表示除法。

2.根据权利要求1所述的方法，其特征在于，步骤1中，设定EMVS阵列具有任意的几何形状。

3.根据权利要求2所述的方法，其特征在于，步骤3还包括：定义选择性矩阵J_q如下：

其中，i_6,q表示为I₆的第q行，得到：

即存在公式(10)被重写为：

4.根据权利要求3所述的方法，其特征在于，步骤3还包括：将公式(9)的结果代入公式(13)，得到：

其中，上标(·)^-1表示逆；等价地，公式(14)表述为：

其中，上标表示伪逆；从的特征分解，能够得到对Φ^(1,2)的估计和T的估计，T的估计为用替换T，并计算公式(15)除第一行除外的左侧部分，得出Φ^(1,3)、Φ^(1,4)、Φ^(1,5)、Φ^(1,6)的估计，从Φ^(1,3)、Φ^(1,4)、Φ^(1,5)、Φ^(1,6)的对角线元素获得对的估计

将b_k重新表示为：

得到：

其中，*表示矢量交叉乘，(·)^*表示共轭；||·||表示向量的F-范数；由于||b_k(1)||²是常数，则：

将的估计值带入(18)中，获得u_k、v_k、w_k的估计值，并分别用表示，二维波达角估计通过下式计算得到：

其中，和分别表示θ_k的估计值和φ_k的估计值。