ABSTRACT

There are large numbers of coherent signals in massive MIMO system. Simple decoherence processing leads to the lower accuracy of the 2-D DOA (Two-Dimensional Direction of Arrival). This paper will propose a modulus equation constraints based dimension reduction MUSIC algorithm which can greatly improve the performance of two-dimensional DOA estimation of coherent signals. The algorithm transforms the two-dimensional DOA estimation problem into an optimization problem, and uses the modulus equation constraints to define the additional conditions and impose strong constraints on the direction vector so that the optimization equation is solved more close to the optimal solution. The results of theoretical analysis and simulation experiments show that the proposed DOA algorithm has high reliability and precision. Such algorithm is able to meet the requirements of DOA estimation performance in massive MIMO system, and can also provide high feasibility and practicability.

Keywords:Massive MIMO, Two-Dimensional Direction of Arrival Estimation, Dimensionality Reduction, MUSCI Algorithm, Modular Equality Constraints

大规模MIMO系统中基于模等式约束的 降维去相干DOA估计

孙鹏帅，余小游^*，林培英，杜青松，马和峰，田丽佳，蒋娅林

湖南大学信息科学与工程学院，湖南 长沙

收稿日期：2018年6月4日；录用日期：2018年6月21日；发布日期：2018年6月28日

摘 要

大规模MIMO系统中存在大量相干信号，简单的去相干处理导致二维波达方向(DOA, Direction of Arrival)估计的精度较低。为此，本文提出一种能够大幅度提高相干信号二维DOA估计性能的模等式约束降维MUSIC算法。该算法将二维DOA估计问题转化为优化问题，并采用模等式约束法定义附加条件，对方向矢量施加较强的约束，使优化方程求解更接近最优解。理论分析和仿真实验结果表明，本文提出的去相干DOA算法的可靠性与精度高，满足大规模MIMO系统的DOA估计性能需求，具有较强的可行性和实用性。

关键词 :大规模MIMO，二维波达方向估计，降维，MUSCI算法，模等式约束

Copyright © 2018 by authors and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY).

http://creativecommons.org/licenses/by/4.0/

1. 引言

大规模MIMO系统中存在大量相干信号，如何得到高精度的二维DOA估计，已经成为大规模MIMO系统研究领域的一个热点问题 [1] 。针对信号源的二维DOA估计，目前已经有了很多研究成果，Tayem and Kwon [2] 提出基于Ｌ型阵列，采用传播方法的波达方向估计，虽然克服了相位模糊问题，但参数匹配困难。王伟提出的降维估计算法 [3] ，能自动匹配参数；但此算法没有将回波信号数据矩阵的维数降至最低。Liang et al. [4] 提出一种自动配对的联合估计算法，但由于需要多次特征值分解和复杂的峰谱搜素过程，导致算法运算量很大。Wang et al. [5] 将SUMWE算法应用于L型阵列，提出了一种有效计算2-D DOA估计方法，该算法首先利用L型阵列接收数据的互相关矩阵的线性组合获得零空间，然后利用一维DOA技术独立估计出俯仰角和方位角。Nie and Li [6] 提出一种计算有效性的方法，利用阵列的共轭对称性增加阵列孔径，提高了估计精度，该方法只用了互相关矩阵的第一列、第一行和对角元素，降低了计算的复杂度。上述后两种算法都是基于信号源是不相关的假设，不适用于相关的信号源。Wei and Guo [7] 通过最小化由两个信号源协方差矩阵构造的代价函数实现多个相干信号的方位角和俯仰角配对，但是该方法构造的信号源协方差矩阵引入了噪声，估计的性能还需要进一步改善。杨艳飞等 [8] 提出了一种基于L型阵列的到达方向估计改进算法，用水平线阵和垂直线阵构建的L型阵列取代两个水平线阵构建的L型阵列，并将L型阵分成四个子阵，计算四个互相关矩阵，构造成特殊的大矩阵，对该矩阵进行特征值分解获得信号子空间估计，再用ESPRIT算法进行方位角和仰角的估计。梁浩等 [9] 提出针对L型阵列的二维角估计问题，在ESPRIT算法基础上设计降维矩阵及回波数据的降维变换，将高维回波数据降至低维信号空间，来获取空间二维角的估计。张艳萍 [10] 等针对二维信号在空间中传播遇到的相干问题进行了降维以及解相干处理，在一定程度了降低了信号谱峰搜索的复杂度以及提升了信号的估计精确度。

为了获得二维DOA估计的最佳性能，降维MUSIC算法是值得深入探讨的算法。基于传统MUSIC算法的二维DOA估计的分辨率较高，但精度较低，其根本原因是MUSIC算法中的二次优化方程中的方向矢量没有受到约束，导致谱峰搜索和参数匹配的效率降低。为此，本文在深入分析降维MUSIC算法中二次优化方程的优化条件的基础上，通过建立基于模等式约束的优化条件，提出基于模等式约束的降维去相干MUSIC算法，该算法可以在降低运算复杂度的同时，提高二维DOS的估计精度。仿真实验结果表明，本文提出的模等式约束去相干降维MUSIC算法的精度和可靠性均有较大幅度的改善，能够满足大规模MIMO系统的实际需要。

2. 信号模型与相关性分析

两个信号 $s_1\left(t\right)$ 、 $s_2\left(t\right)$ 的相关程度由相关系数 衡量，定义为 [11] ：

$r=\frac{\mathrm{cov}\left(s_1,s_2\right)}{\sqrt{D\left(s_1\right)D\left(s_2\right)}}=\frac{E\left\{\left[s_1-E\left\{s_1\right\}\right]\left[s_2-E\left\{s_2\right\}\right]\right\}}{\sqrt{E\left\{{\left[s_1-E\left\{s_1\right\}\right]}^2\right\}E\left\{{\left[s_2-E\left\{s_2\right\}\right]}^2\right\}}}$ (1)

当 $r=0$ 时，信号不相关；当 $0<r<1$ 时，信号相关；当 $r=1$ 时，信号相干。

由上述定义可以得到如下推论：1) 频率不相同的信号是不相关的；2) 同频但初相位不同的信号是相关的。不包括特殊情况：初相位之差为±π/2时，信号不相关；初相相差±π时，信号相干；3) 频率和初相位都相同的信号是相干的。

假设L型均匀等距线阵的X轴和Y轴上均有M个阵元，则共有2M-1个阵元。阵元间距为d，其天线坐标系模型如图1所示。假设空间中有K个远场窄带相干信号入射到此L型阵列上，其二维DOA为 $\left({\theta }_k,{\phi }_k\right)$ ， $k=1,2,\cdots ,K$ 。其中 ${\theta }_k,{\phi }_k$ 分别表示第K个信号的到达俯仰角和方位角 $\left({\theta }_k\in \left[0,\text{π}\right],{\phi }_k\in \left[0,\text{π}/2\right]\right)$ 。以坐标原点O处的阵元为参考阵元，则X轴上的阵列接收信号可表示为：

$X\left(t\right)=A_{x}S\left(t\right)+N_x\left(t\right)$ (2)

式中， $S\left(t\right)={\left[s_1\left(t\right),s_1\left(t\right),\cdots ,s_K\left(t\right)\right]}^{\text{T}}$ 和 $N_x\left(t\right)$ 分别表示空间中存在的相应的相干信号源矢量和接收噪声矩阵。

同样，Y轴上的阵列接收信号表示为：

$Y\left(t\right)=A_{y}S\left(t\right)+N_y\left(t\right)$ (3)

式中， $S\left(t\right)$ 和 $N_y$ 分别表示空间中存在的相应的相干信号源矢量和接收噪声矩阵。

导向矢量矩阵 $A_x\in C^{M\times K}$ 和 $A_y\in C^{M\times K}$ 相等，可表示为 [12] ：

图1. L型阵列天线坐标示意图

$\begin{array}{l}{a}_x\left({\theta }_k\right)={\left[a_{x,1}\left({\theta }_k\right),a_{x,2}\left({\theta }_k\right),\cdots ,a_{x,M}\left({\theta }_k\right)\right]}^{\text{T}}={\left[1,{\text{e}}^{-j2\text{π}d\cos {\phi }_1\sin {\theta }_1/\lambda },\cdots ,{\text{e}}^{-j2\text{π}\left(M-1\right)d\cos {\phi }_k\sin {\theta }_k/\lambda }\right]}^{\text{T}}\\ a_y\left({\phi }_k\right)={\left[a_{y,1}\left({\phi }_k\right),a_{y,2}\left({\phi }_k\right),\dots a_{y,M}\left({\phi }_k\right)\right]}^{\text{T}}={\left[1,{\text{e}}^{-j2\text{π}d\sin {\phi }_1\sin {\theta }_1/\lambda },\cdots ,{\text{e}}^{-j2\text{π}\left(M-1\right)d\sin {\phi }_k\sin {\theta }_k/\lambda }\right]}^{\text{T}}\end{array}$ (4)

假设第K个入射信号的入射方向与x轴夹角为 ${\alpha }_k$ ，与y轴夹角为 ${\beta }_k$ ，则有：

$\{\begin{array}{l}\cos {\alpha }_k=\cos {\phi }_k\sin {\theta }_k\\ \cos {\beta }_k=\sin {\phi }_k\sin {\theta }_k\end{array}$ (5)

3. 去相干算法

由于阵列接收到的信号为窄带相干信号，传统的DOA估计方法无法进行估计，需要利用前后向空间平滑技术对阵列接收信号X，Y进行解相干。先对X轴上的阵列信号进行前后向空间平滑处理，然后将X轴上的均匀线阵(M个阵元)分成相互交错的P个子阵，则每个子阵的阵元数为 $m=M-P+2$ ，取第一个子阵作为参考，如图2所示。

则对于K个子阵有数据模型：

$x_k\left(t\right)=\left[x_k{x}_{k+1}\cdots x_{k+m-1}\right]=A_{xk}{D}^{\left(k-1\right)}s\left(t\right)+n_k\left(t\right)$ (6)

其中 $D=\text{diag}\left({\eta }_1,{\eta }_1,\cdots ,{\eta }_K\right)$ ， ${\eta }_i=\exp \left(\left(j2\text{π}d/\lambda \right)\sin {\alpha }_i\right)$ ，于是该子阵数据协方差矩阵为：

$R_{xk}=A_{xk}{\eta }^{\left(k-1\right)}{R}_s{\left(D^{\left(k-1\right)}\right)}^H{A}_{xk}^H+{\sigma }^{2}I$ (7)

前向空间平滑方法对满秩协方差矩阵的恢复是通过求各子阵协方差矩阵的均值来实现的，即取前向平滑修正的协方差矩阵为：

$R_x^f=\frac{1}{p}{\displaystyle {\sum }_{i=1}^p{R}_{ip}}=A_{xk}{R}_{xs}^f{A}_{xk}^H+{\sigma }^{2}I$ (8)

同理，后向平滑修正的协方差矩阵为：

$R_x^b=\frac{1}{p}{\displaystyle {\sum }_{i=1}^p{R}_{xp}^b}=A_{xk}{R}_{xs}^b{A}_{xk}^H+{\sigma }^{2}I$ (9)

由文献 [1] 可知 $R_x^b$ 是 $R_x^f$ 的共轭倒序阵，它们之间满足共轭倒序不变性。因此可得X轴阵列接收的信号经前后向空间平滑处理后得到的数据协方差矩阵为：

$R_x^{fb}=\frac{1}{2}\left(R_x^f+R_x^b\right)$ (10)

图2. 空间平滑算法模型

同理，对Y轴阵列接收到的信号经前后空间平滑处理后得到的数据协方差矩阵为:

$R_y^{fb}=\frac{1}{2}\left(R_y^f+R_y^b\right)$ (11)

于是，L型阵列得到的信号协方差矩阵为：

$Z={\left[R_x^{fb}{R}_y^{fb}\right]}^{\text{T}}$ (12)

4. 降维MUSIC算法

根据一维MUSIC求解空间谱算法，对L型阵列信号协方差矩阵Z进行特征值分解，可得 $Z=U_s{\Lambda }_s{U}_s^H+U_n{\Lambda }_n{U}_n^H$ ，其中 ${\text{Λ}}_s$ 和 ${\text{Λ}}_n$ 为特征值矩阵，分别对应K个大特征值，和( $2M-1-K$ )个小特征值； $U_s$ 为K个大特征值所对应的特征矢量矩阵，其张成的空间为信号子空间， $U_n$ 为( $2M-1-K$ )个小特征值所对应的特征矢量矩阵，其张成的空间为噪声子空间。容易证明信号子空间的方向矢量 $a_x\left({\theta }_k\right)\otimes a_y\left({\phi }_k\right)$ 和噪声子空间 $U_n$ 是相互正交的，经典子空间类DOA估计方法正是利用这一性质对2维DOA进行估计求解的。

传统2维MUSIC算法的空间谱估计公式 [13] 为：

$P_{2dmusic}=\frac{1}{{\left[a_x\left({\theta }_k\right)\otimes a_y\left({\phi }_k\right)\right]}^H{U}_n{U}_n^H\left[a_x\left({\theta }_k\right)\otimes a_y\left({\phi }_k\right)\right]}$ (13)

但是该算法需要进行2维角度搜索，运算量大。为了将2维角度估计问题降维求解，首先将其转化为二次优化问题，定义检测量为：

$Q\left({\theta }_k,{\phi }_k\right)={\left[a_x\left({\theta }_k\right)\otimes a_y\left({\phi }_k\right)\right]}^H{U}_n{U}_n^H\left[a_x\left({\theta }_k\right)\otimes a_y\left({\phi }_k\right)\right]$ (14)

将式(13)作如下变形：

$Q\left({\theta }_k,{\phi }_k\right)=a_y^H\left({\phi }_k\right){\left[a_x\left({\theta }_k\right)\otimes I_M\right]}^H{U}_n{U}_n^H\left[a_x\left({\theta }_k\right)\otimes I_M\right]a_y\left({\phi }_k\right)$ (15)

令

$G\left({\theta }_k\right)={\left[a_x\left({\theta }_k\right)\otimes I_M\right]}^H{U}_n{U}_n^H\left[a_x\left({\theta }_k\right)\otimes I_M\right]$ (16)

则2维DOA估计问题可等价为求解：

$\min a_y^H\left({\phi }_k\right)G\left({\theta }_k\right)a_y\left({\phi }_k\right)=\min Q\left({\theta }_k,{\phi }_k\right)$ (17)

为使得式(17)能够通过1维搜索即可得到 ${\theta }_k$ 和 ${\phi }_k$ ，文献 [14] 对式(17)增加了附加条件： $e_1^H{a}_y\left({\phi }_k\right)=1$ ，其中： $e_1={\left[1,0,\cdots ,0\right]}^{\text{T}}$ ，进而将角度估计问题转化为角度优化问题来求解。但该条件相当于仅对方向向量 $a_y$ 的第一个元素 $a_y\left(1\right)=1$ 进行了限制。经分析，方向向量还应该满足以下条件：

$a_y\left(y\right)=\exp \left[-j2\text{π}d\left(m-1\right)\cos {\beta }_k/\lambda \right]={\left[\exp \left(-j2\text{π}d\cos {\beta }_k/\lambda \right)\right]}^{m-1}=a_y^{y-1}=\{\begin{array}{l}1,m=1\\ a_y\left(2\right),m=2\\ a_y\left(2\right)a_y\left(y-1\right),m\ge 3\end{array}$ (18)

其中， $a_y\left(2\right)=\exp \left(-j2\text{π}d\cos {\beta }_k/\lambda \right)$ 。文献 [10] 的运算模型中方向向量 $a_y$ 是待求解向量，但对 $a_y$ 的限定条件中，仅限定其首元素为1，而对该向量中其它元素应该满足的条件未作限制，这在一定程度上降低了估计的精度甚至会造成错误的估计。

5. 改进的降维MUSIC (MRD-MUSIC)算法

为了加强对待求解向量中元素的约束，将模约束的稳健Capon算法 [15] 引入到降维MUSIC算法中，因此约束条件可以表示为：

$\{\begin{array}{l}{\min }_w{a}_y^{H}G\left({\theta }_k\right)a_y\\ \text{s}\text{.t}\text{.}{\Vert a_y\left({\phi }_k\right)\Vert }^2=M\end{array}$ (19)

为了求解该优化问题，构造代价函数：

$L\left(a_y,{\theta }_k\right)=a_y^{H}G\left({\theta }_k\right)a_y-\lambda \left(a_y^H{a}_y-M\right)$ (20)

其中λ为一常数。利用拉格朗日乘数法求 $L\left(a_y\left({\phi }_k\right),{\theta }_k\right)$ 的极值，对式(20)求导，并令结果为零，可得：

$\frac{\partial }{\partial a_y}L\left(a_y,{\theta }_k\right)=2G\left({\theta }_k\right)\cdot a_y-2\lambda \cdot a_y=0$ (21)

即

$G\left({\theta }_k\right)\cdot a_y=\lambda \cdot a_y$ (22)

显然， $a_y$ 为 $G\left({\theta }_k\right)$ 的特征值 所对应的特征向量，由于

$a_y^H\cdot G\left({\theta }_k\right)\cdot a_y=\lambda \cdot a_y^H\cdot a_y=\lambda \cdot M$ (23)

因此，求目标函数 $a_y^{H}G\left({\theta }_k\right)a_y$ 的最小值，相当于求解λ最小。所求的 $a_y$ 正好为 $G\left({\theta }_k\right)$ 的最小特征值所对应的特征向量，

即

$a_y=P_{\min }\left[G\left({\theta }_k\right)\right]$ (24)

其中， $P_{\min }\left(G\right)$ 表示求解G的最小特征值所对应的特征向量。

将上述结果代入到原优化问题的目标函数，可得：

${\widehat{\theta }}_k=\arg \min P_{\min }^H\left[G\left({\theta }_k\right)\right]\cdot G\left({\theta }_k\right)P_{\min }\cdot \left[G\left({\theta }_k\right)\right]=\arg \max \frac{1}{P_{\min }^H\left[G\left({\theta }_k\right)\right]\cdot G\left({\theta }_k\right)P_{\min }\cdot \left[G\left({\theta }_k\right)\right]}$ (25)

根据公式(22)可知，在 $\left({\theta }^{\circ },{180}^{\circ }\right]$ 范围内对 ${\theta }_k$ 进行搜索，其中较大的K个峰值点所对应的角度值就是所求的各个入射信号的 ${\theta }_k$ 。

通过以上求解可以得到所有的θ角，然后通过θ角进一步估计所有的φ值。因此首先求解入射信号与y轴之间夹角的方向向量 $a_y\left({\phi }_k\right)$ 。通过以上的分析可知， $G\left({\widehat{\theta }}_k\right)$ 特征值分解对应的最小特征向量 $P_{\min }\left[G\left({\widehat{\theta }}_k\right)\right]$ 即是 $a_y\left({\phi }_k\right)$ 。之后通过最小二乘法求解即可。由于θ和φ两角存在一一对应的关系，因此在求解的过程中并不会存在额外的角度配对问题。

6. 实验仿真与结果分析

仿真过程中，阵列为L型均匀阵列，共128个阵元，阵元间距为半波长。均方根误差是描述估计性能的重要指标，设蒙特卡罗实验次数为Num，第Num次试验中2维DOA的实际值为 $\left(\theta ,\phi \right)$ ，估计值为 $\left(\widehat{\theta },\widehat{\phi }\right)$ ，则角度的均方根误差定义为：

$\text{RMSE}=\frac{1}{K}{\displaystyle {\sum }_{k=1}^K\sqrt{\frac{1}{Num}{\displaystyle {\sum }_{num=1}^{Num}{\left({\widehat{\theta }}_{k,num}-{\theta }_k\right)}^2}+{\left({\widehat{\phi }}_{k,num}-{\phi }_k\right)}^2}}$ (26)

式中，RMSE表示方位角和俯仰角的估计均方误差。

利用K个窄带入射信号进行实验，进行θ角度搜索时，取其中较大的K个峰值点所对应的角度值作为角度测量值。如果所有角度测量与实际值的误差均在1˚范围内，则定义为实验成功，否则为失败。经过大量试验后，成功概率定义为试验成功次数所占总实验的比值。

二维角度估计：假设有3个远场窄带相干信号，水平角和俯仰角分别为(10˚, 15˚)、(20˚, 40˚)、(30˚, 85˚) (图3)，x轴方向阵元数为8，y轴方向阵元数为8，快拍数L为100，信噪比取10dB，蒙特卡罗实验次数为500次。

可靠性：针对2D-MUSIC去相干算法，文献 [10] 去相干算法和本文算法成功率进行对比实验，信噪比SNR范围为-6dB到6dB。在每个信噪比条件下对3个窄带信号进行500次蒙特卡罗实验，得到文献 [10] 和本文MRD-MUSIC算法的成功率对比。

从图4中可以看出，在低信噪情况下，文献 [10] 算法具有较低的成功率，而本文算法的成功率接近2D-MUSIC算法。且在−2 dB时候成功率达到了98%，之后一直保持在较高的成功率。

DOA精度：令信噪比SNR范围从−10 dB到10 dB，快拍数为1000，进行500次蒙特卡罗仿真实验。将文献 [10] 去相干算法和本文算法的均方根误差曲线进行对比。对比如图5所示。

由图5可知，在低信噪比时，本文算法仍具有较低的RMSE，同时随着信噪比的增加文献 [10] 的RMSE在逐渐降低，但是估计效果相对于本文算法依然较差。

同时改变X轴和Y轴的阵元数M，使它们从3变化到20，对比文献 [10] 算法和本文算法的阵元数、均方根误差曲线。

由图6可知，随着阵元数的增加，文献 [10] 和本文所提出算法的RMSE都逐渐降低，但是本文算法在相同阵元数的情况下具有较低的RMSE。

图3. 俯仰角和方位角估计

图4. 成功率对比

图5. RMSE随信噪比变化曲线

图6. RMSE随阵元数变化曲线

7. 结论

本文提出了一种改进的去相干降维MUSIC算法，该算法首先对接收到的信号进行解相干，然后利用二次优化算法将2维DOA估计降为1维DOA估计，并且入射角φ的求解是通过相应θ带入公式求得的，因此不存在角度配对问题。同时该算法加强了方向向量的约束，使得估计值更加精确。仿真实验验证了本文针对相干信号估计的有效性和正确性，因此是一种高效的去相干降维DOA算法。

文章引用

孙鹏帅,余小游,林培英,杜青松,马和峰,田丽佳,蒋娅林. 大规模MIMO系统中基于模等式约束的降维去相干DOA估计
Modulus Equation Constraints Based Decoherent Dimension Reduction 2-D DOA Estimation for Massive MIMO Systems[J]. 无线通信, 2018, 08(03): 123-132. https://doi.org/10.12677/HJWC.2018.83014

参考文献