CN106291453B - 一种基于波束空间变换的二维波达角测定方法 - Google Patents

Abstract

本发明公开了一种基于波束空间变换的二维波达角测定方法，适用于大规模MIMO系统测定非相干信号二维波达角的测定方法。包括：设置一个具有M根天线的均匀矩形阵列，对接收的数据进行波束空间变换；建立变换矢量的协方差矩阵；通过协方矩阵证获取信号子空间；根据波束空间变换的特性和阵列的结构，得到相应的广义选择信号子空间；然后获得旋转矩阵；最后分别确定信号源的俯仰角与方位角。本发明在大规模MIMO场景下具有较小的测定误差，避免了非线性优化和角度搜索。因此，本发明具有精度高，数据处理复杂度低等优点，应用于大规模MIMO系统可有效提高系统的性能和实用价值。

Description

一种基于波束空间变换的二维波达角测定方法

技术领域

本发明属于无线移动通信技术领域，特别是涉及一种采用大规模阵列天线测定非相干信号二维波达角(Direction of Arrival，简称DOA)的测定方法。

背景技术

随着移动通信技术的发展，特别是近来第四代移动通信系统(4G)的发展，人们对通信的速率和质量要求越来越高，对通信业务的需求量越来越大。在第四代移动通信系统中，频谱利用率得到很大的提到，主要得益于多输入多输出(MIMO)天线技术。但是，随着业务量的增加，传统MIMO阵列已经渐渐无法满足传输速率增加的需求。加之越来越多的科研人员对第五代移动通信系统(5G)的探索，一种被称为大规模MIMO的技术越来越得到重视，该技术具有很高的频谱利用率和链路可靠性。在大规模MIMO系统中，基站配置成百上千根天线阵元，同时服务几十个用户。但是大规模MIMO系统的部署受到基站空间的限制。比如，在经典的2.5GHz LTE载波场景中，在水平方向部署32个间距为半波长的阵元占据近1.9米，当天线阵元比较多时，这种线型布局显然不合适。所以二维均匀矩形阵列(URA)更加受到科研人员的青睐。在大规模MIMO系统所具有的众多优越性中，用于提高链路可靠性的三维波束成形技术受到了较大的关注。要实施三维波束成形技术，首先需要对信号源的二维波达角(DOA)进行精确的估计，包括方位角和俯仰角。另一方面，在时分双工(TDD)大规模MIMO系统中，精确的DOA估计在抑制导频污染方面也具有重要的作用。

在近几十年中，人们已经对对于非相干独立信号的波达角测定问题做了大量的研究并提出了不少的方法。例如，ESPRIT(rotational variance technique)方法，MUSIC(multiple signal classfication)方法，基于波束成形(beam-forming)的方法，基于WSF(weighted subspace fitting)的方法和基于最大似然估计(ML)的方法等。但是现有算法大多需要一维或者二维的谱搜索，有时还设置多维非线性优化，导致数据的处理复杂度比较高，难以实现实时的DOA估计。对于大规模MIMO系统来说，其天线数非常大，传统算法的复杂度是难以承受的。另一方满，现在虽然有基于传播算子的低复杂度算法，但是由于其估计精度较低，一般用于DOA预估或者对DOA估计精度要求不高的场景。

发明内容

本发明的目的是公开一种基于波束空间变换的适用于大规模MIMO系统的非相干信源的二维波达角测定方法，以达到有效地提高非相干信源的二维波达角测定的实时性和精度并降低数据处理的计算复杂度的目的，应用于大规模MIMO系统中，可以有效地提高系统的性能及实际应用的价值等目的。

本发明首先对阵列接收到的数据进行波束空间变换并计算变换后的数据协方差矩阵；然后确定信号子空间；再根据波束空间变换的特性获得广义选择信号子空间之间的旋转矩阵；然后对旋转矩阵进行特征值分解并对特征值进行配对；最后根据旋转矩阵的特征值估计信号源的俯仰角和方位角。本发明的方案不仅不需用到一维或者二维的角度搜索，而且同元素空间算法相比具有更低的计算复杂度，并且提高测定的精度和实效率，从而实现发明的目的。因此，本发明方法包括下列步骤：

步骤1：通过均匀矩形天线阵列获取K个非相干信源的接收数据矢量，其中均匀矩形天线阵列包括M_x×M_y个阵元，M_x对应矩形天线阵列的行，M_y对应矩形天线阵列的列；

对接收数据矢量进行波束域变换：基于预设的M_x×P维的波束空间变换矩阵对接收数据矢量进行波束空间变换，得到变换矢量，其中P为压缩系数，且P＜M_x；

步骤2：计算变换矢量在N次采样下的协方差矩阵；

步骤3：对协方差矩阵进行特征值分解，基于前K个最大特征值所对应的特征向量e₁,e₂,…e_K构建信号子空间E_S＝[e₁,e₂,...,e_K]；

步骤4：基于波束空间变换矩阵和天线阵列结构，根据信号子空间E_S构造广义选择信号子空间E_Sf，f＝1,2,3,4；

步骤5：基于E_S1和E_S2构建第一旋转矩阵Ψ₁，基于E_S3和E_S4构建第二旋转矩阵Ψ₂；

步骤6：分别对第一旋转矩阵、第二旋转矩阵进行特征值分解，得到对角矩阵Λ₁、Λ₂，其中对角矩阵Λ₁的对角元素为ψ₁的特征值，对角矩阵Λ₂的对角元素为ψ₂的特征值；

对Λ₁和Λ₂的对角元素进行配对对齐，得到对齐的对角矩阵Λ₁和

步骤7：根据对齐的对角矩阵Λ₁和计算各个信号源的俯仰角方位角其中k＝1,2,...,K：

定义ξ_1,k和ξ_2,k分别为Λ₁和的第k个对角线元素，则 其中u＝2πd/λ，d为相邻阵元之间的间距，λ为载波波长。

其中，波束空间变换矩阵为：

且 表示纬度为M_x×P的实数矩阵，即C的上标用于标识矩阵的纬度，下同；其中P(P＜M_x)为压缩系数。

波束空间变换定义为：

其中x_b(t)表示变换矢量，是块对角矩阵，其块对角元素均为(·)^H表示矩阵共轭转置，函数blkdiag(·)表示生成指定对角线元素的矩阵；x(t)为均匀矩形天线阵列的各个阵元接收到数据按阵元的顺序排成一列，第二行阵列接受的数据排在第一行接收的数据的后边，以此类推，即将M_x×M_y的均匀矩形天线阵列的每一行的接收数据排成一列构成接收数据矢量x(t)。

变换矢量x_b(t)在N次采样下的协方差矩阵为：其中N为采样快拍数，t为采样序号。

广义选择信号子空间E_Sf具体为：E_s1＝ΩE_s、E_s3＝J₃E_s、E_s4＝J₄E_s；

其中

其中Ω和为块对角矩阵，对角元素分别为QF^H和Q，为P(M_y-1)×P(M_y-1)的单位矩阵，为P(M_y-1)×P的零矩阵。Q和F的表达式如下

Q＝J_QQ_d∈C^P×P

F＝diag(exp(-j2π/M_x),…,exp(-j2πP/M_x))∈C^P×P

其中

Q_d＝diag(exp(j2π(M_x-1)/M_x),…exp(j2πP(M_x-1)/M_x))∈C^P×P

在基于E_S1和E_S2构建第一旋转矩阵Ψ₁，基于E_S3和E_S4构建第二旋转矩阵Ψ₂时，可由总体最小二乘法(Total Least-Squares，TLS)得到。旋转矩阵Ψ₁和Ψ₂具有相同的计算方式，这里仅给出Ψ₁的计算步骤：

(a)建立矩阵E_S12＝[E_S1,E_S2]；

(b)定义矩阵并对其进行特征分解得到F＝GΛ_FG^H；其中Λ_F是一个对角阵，G是特征向量矩阵；将G分成4个K×K的子矩阵,即

(c)建立矩阵

在实现旋转矩阵ψ₁和ψ₂的对角矩阵Λ₁和Λ₂的对角元素的配对对齐时，可以通过如下两种方式实现：

方式一：

(a)对旋转矩阵进行特征值分解 特征值均为降序排列；

(b)定义矩阵Ψ₃＝Ψ₁Ψ₂，

(c)对Ψ₃进行特征值分解并定义矩阵设置p＝0；

(d)更新p＝p+1，计算乘积和商其中然后根据最小二乘法匹配和 和其中即根据找到Λ₁与Λ₂之间的对应关系；令其中是对角矩阵；

(e)重复步骤(d)直至p＝K。所得Λ₁与对角线元素是配对的，即[Λ₁]_k,k与对应于同一个信号源。

方式二：

(a)对旋转矩阵进行特征值分解 令u_i＝[u_i1,u_i2,...,u_iK]^T表示Ψ₁的第i个特征值ξ_1,i对应的特征向量，设置p＝0。

(b)更新p＝p+1，令w_p＝Ψ₂u_p＝[w_p1,w_p2,...,w_pK]，计算从Λ₂中选择与差值最小的元素ξ_2,j作为与ξ_1,p配对的元素，令

(c)重复步骤(b)直至p＝K。所得Λ₁与对角线元素是配对的，即[Λ₁]_k,k与对应于同一个信号源。

综上所述，由于采用了上述技术方案，本发明的有益效果是：本发明在大规模MIMO场景下具有较小的测定误差，避免了非线性优化和角度搜索。因此，本发明具有精度高，数据处理复杂度低等优点，应用于大规模MIMO系统可有效提高系统的性能和实用价值。

附图说明

图1是本发明方法的阵列设置示意图；

图2是本发明实施方式仿真运行所得方位角估计误差随信噪比变化轨迹意图；

图3是本发明实施方式仿真运行所得俯仰角估计误差随信噪比变化轨迹意图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，下面结合实施方式和附图，对本发明作进一步地详细描述。

将本发明用于如图1所示的二维均匀矩形阵列，图1中共有M＝100个阵元，在x方向有M_x＝10行，在y轴方向有M_y＝10列。信号波长λ为0.375m，相邻阵元在x与在y方向的间距为d＝λ/2＝0.1875m。本实施例中，有K＝2个窄带非相干以不同方向入射到此阵列，第一个入射信号的方位角和俯仰角分别为25°和50°，第二个入射信号的方位角和俯仰角分别为62°和45°，各阵元上的噪声为加性高斯白噪声(AWGN)，且噪声与信号不相关，信号源均为二进制相位调制(BPSK)信号，阵列的接收信噪比为10dB。

因此，在t时刻，天线阵列接收到的数据矢量为：其中，s_k(t)为第k个信号源，n(t)为独立高斯白噪声，a(θ_k(t),φ_k(t))是相对于对第k个信号的阵列流形矢量，即[a(θ_k(t),φ_k(t))]_m＝exp(iusin(φ_k(t))[(m_x-1)cos(θ_k(t))+(m_y-1)sin(θ_k(t))])，其中m＝(m_y-1)M_x+m_x,m_x＝1,2...M_x,m_y＝1,2...M_y，θ_k和φ_k分别为第k个信号的方位角和俯仰角。

天线阵列接收到的数据矢量可以写成矩阵形式：x(t)＝AS(t)+n(t)，其中A＝[a(θ₁(t),φ₁(t)),a(θ₂(t),φ₂(t))]为阵列流形矩阵，S(t)＝[s₁(t),s₂(t)]^T为信号源矢量。

基于图1所示的天线阵列的接收数据矢量完成二维波达角测定的具体过程如下：

步骤1：对接收数据矢量进行波束域变换：

令P＝7,则变换矩阵为

所以波束空间变换表达式为：

步骤2：计算变换矢量x_b(t)的协方差矩阵：

本实施例中，取快拍数N＝500，协方差矩阵通过下式获得：其中t为采样的序号。

步骤3：确定信号子空间：

对步骤2得到的协方差矩阵进行特征值分解，并利用最大的2个特征值对应的特征向量e₁,e₂作为列建立信号子空间矩阵E_s＝{e₁,e₂}。

步骤4：确定广义选择信号子空间：

广义选择信号子空间的可以表示为：E_s1＝ΩE_s、E_s3＝J₃E_s、E_s4＝J₄E_s；

其中

Ω＝blkdiag(QF^H,QF^H,…,QF^H)∈C^70×70

J₃＝[I₆₃,0_63×7]∈C^63×70

J₄＝[0_63×7,I₆₃]∈C^63×70

上式中，Ω和为块对角矩阵，对角元素分别为QF^H和Q，I₆₃为63×63的单位矩阵，0_63×7为63×7的零矩阵。Q和F的表达式如下

Q＝J_QQ_d∈C^7×7

F＝diag(exp(-j2π/10),…,exp(-j2π×7/10))∈C^7×7

其中

Q_d＝diag(exp(j2π(10-1)/10),…exp(j2π×7(10-1)/10))∈C^7×7

步骤5：确定旋转矩阵：

利用步骤4得到的子矩阵E_S1和E_S2通过总体最小二乘法确定旋转矩阵Ψ₁，同时用同样的方式确定E_S3和E_S4的旋转矩阵Ψ₂，这里仅给出计算Ψ₁的具体步骤：

(1)建立矩阵E_S12＝[E_S1,E_S2]；

(2)定义矩阵进行特征分解得到F＝GΛ_FG^H；其中Λ_F是一个对角阵，G是特征向量矩阵；将G分成4个2×2的子矩阵,即

(3)建立矩阵

步骤6：特征值配对：

(a)对旋转矩阵进行特征值分解 特征值均为降序排列；

(b)定义矩阵Ψ₃＝Ψ₁Ψ₂，

(c)对Ψ₃进行特征值分解并定义矩阵设置p＝0；

(d)更新p＝p+1，计算乘积和商其中然后根据最小二乘法匹配和 和其中即根据找到Λ₁与Λ₂之间的对应关系；令其中是对角矩阵；

(e)重复步骤(d)直至p＝2。所得Λ₁与对角线元素是配对的，即[Λ₁]_k,k与对应于同一个信号源。

步骤7：确定俯仰角和方位：

定义ξ_1,k和ξ_2,k分别为Λ₁，的第k(k＝1,2)个对角线元素，对应于第k个信号，则第k个信号源的方位角和俯仰角为：其中 的上标“^”用于表示该方位角、俯仰角为计算得到，即估计值；

分别得到两个信号源的估计角度：

为了评估方法的性能，经过100次测定所得的两个信号的方位角和俯仰角的统计平均值分别为：

相应的中心方位角和中心俯仰角的均方根误差分别为：

为了进一步验证算法的性能，在天线数为100的情形下，进行100次独立实验验证俯仰角和方位角估计误差随信噪比变化的轨迹，其结果如图2,3所示。本发明在大规模MIMO场景下具有较小的测定误差，避免了非线性优化和角度搜索。

以上所述，仅为本发明的具体实施方式，本说明书中所公开的任一特征，除非特别叙述，均可被其他等效或具有类似目的的替代特征加以替换；所公开的所有特征、或所有方法或过程中的步骤，除了互相排斥的特征和/或步骤以外，均可以任何方式组合。

Claims (4)

1.一种基于波束空间变换的二维波达角测定方法，其特征在于，包括下列步骤：

步骤1：通过均匀矩形天线阵列获取K个非相干信源的接收数据矢量，其中均匀矩形天线阵列包括M_x×M_y个阵元，M_x对应矩形天线阵列的行，M_y对应矩形天线阵列的列；

对接收数据矢量进行波束域变换：基于预设的M_x×P维的波束空间变换矩阵对接收数据矢量进行波束空间变换，得到变换矢量，其中P为压缩系数，且P＜M_x；

所述M_x×P维的波束空间变换矩阵为：

其中j表示虚数单位，p∈{1,2,…,P}；

波束空间变换定义为：其中x_b(t)表示变换矢量，其中且为PM_y×M_xM_y维的块对角矩阵；将M_x×M_y的均匀矩形天线阵列的每一行的接收数据排成一列构成接收数据矢量x(t)；符号(·)^H表示矩阵共轭转置，函数blkdiag(·)表示生成指定对角线元素的矩阵；

步骤2：计算变换矢量在N次采样下的协方差矩阵；

步骤3：对协方差矩阵进行特征值分解，基于前K个最大特征值所对应的特征向量e₁,e₂,…e_K构建信号子空间E_S＝[e₁,e₂,...,e_K]；

步骤4：基于波束空间变换矩阵和天线阵列结构，根据信号子空间E_S构造广义选择信号子空间E_Sf，f＝1,2,3,4；

步骤5：基于E_S1和E_S2构建第一旋转矩阵Ψ₁，基于E_S3和E_S4构建第二旋转矩阵Ψ₂；

步骤6：分别对第一旋转矩阵、第二旋转矩阵进行特征值分解，得到对角矩阵Λ₁、Λ₂，其中对角矩阵Λ₁的对角元素为ψ₁的特征值，对角矩阵Λ₂的对角元素为ψ₂的特征值；

对Λ₁和Λ₂的对角元素进行配对对齐，得到对齐的对角矩阵Λ₁和

步骤7：根据对齐的对角矩阵Λ₁和计算各个信号源的俯仰角方位角其中k＝1,2,...,K：

定义ξ_1,k和ξ_2,k分别为Λ₁和的第k个对角线元素，则 其中u＝2πd/λ，d为相邻阵元之间的间距，λ为载波波长。

2.如权利要求1所述的方法，其特征在于，广义选择信号子空间E_Sf具体为：E_s1＝ΩE_s、E_s3＝J₃E_s、E_s4＝J₄E_s；

Ω＝blkdiag(QF^H,QF^H,…,QF^H)，且块对角矩阵Ω的维度为PM_y×PM_y；

且块对角矩阵的维度为PM_y×PM_y；

且J₃的维度为P(M_y-1)×PM_y；

且J₄的维度为P(M_y-1)×PM_y；

矩阵Q＝J_QQ_d，其中Q_d＝diag(exp(j2π(M_x-1)/M_x),…exp(j2πP(M_x-1)/M_x))，矩阵F＝diag(exp(-j2π/M_x),…,exp(-j2πP/M_x))；

且矩阵Q、J_Q、Q_d、F的维度均为P×P；

为P(M_y-1)×P(M_y-1)的单位矩阵，为P(M_y-1)×P的零矩阵。

3.如权利要求2所述的方法，其特征在于，通过如下步骤对Λ₁和Λ₂的对角元素进行配对对齐，得到对齐的对角矩阵Λ₁和

(a)对旋转矩阵进行特征值分解其中T₁、T₂分别表示矩阵Ψ₁、Ψ₂的特征向量矩阵；

令u_i表示Ψ₁的第i个特征值ξ_1,i对应的特征向量，则u_i＝[u_i1,u_i2,...,u_iK]^T，其中u_i1,u_i2,...,u_iK表示u_i的K个向量元素；并设置p＝0；

(b)更新p＝p+1，令w_p＝Ψ₂u_p，则w_p＝[w_p1,w_p2,...,w_pK]，其中w_p1,w_p2,...,w_pK表示w_p的K个向量元素，计算从Λ₂中选择与差值最小的元素ξ_2,j作为与ξ_1,p配对的元素，令

(c)重复步骤(b)直至p＝K。

4.如权利要求2所述的方法，其特征在于，通过如下步骤对Λ₁和Λ₂的对角元素进行配对对齐，得到对齐的对角矩阵Λ₁和

(a)对旋转矩阵进行特征值分解特征值均为降序排列，其中T₁、T₂分别表示矩阵Ψ₁、Ψ₂的特征向量矩阵；

(b)定义矩阵Ψ₃＝Ψ₁Ψ₂，

(c)对Ψ₃进行特征值分解并定义矩阵设置p＝0；其中T₃表示矩阵Ψ₃的特征向量矩阵，Λ₃表示矩阵Ψ₃的对角矩阵；

(d)更新p＝p+1；计算乘积和商其中然后根据最小二乘法匹配和 和其中即根据找到Λ₁与Λ₂之间的对应关系(x_p,y_p)；令其中为K×K的对角矩阵；

(e)重复步骤(d)直至p＝K。