大规模多输入多输出(MIMO，multiple-input multiple-output) 技术作为5G移动通信的关键技术之一，逐渐成为国内外研究的热点^[1-2].在大规模MIMO系统中，随着基站端天线数大幅增加，各个信道之间渐近正交^[3]，传统的最小均方误差(MMSE，minimum mean square error) 检测方案可用于大规模MIMO系统中，并获得较好的检测性能.但是，MMSE检测算法引入了复杂的矩阵求逆运算，其计算复杂度随着发射天线数呈立方增长，在大规模MIMO中难以快速有效地实现.

为了降低计算复杂度，首先介绍了一些矩阵求逆简化算法，如基于多项式展开的Neumann级数展开和多级线性接收，以及基于迭代方法的Richardson算法和松弛迭代算法来简化或替代矩阵求逆，比较了它们在硬判决方案下的检测性能和计算复杂度.在计算复杂度为O(K²) 时，松弛迭代算法的检测性能最优，多级线性接收次之，Richardson算法和Neumann级数展开的检测性能相对较差.在此基础上，进一步提出了将其运用于软判决中，在降低复杂度的同时获得相对于硬判决更好的检测性能，最后通过仿真验证了检测性能.

1 系统模型

考虑大规模MIMO系统的上行链路，该系统由一个部署N根天线的基站和K个单天线用户组成.令x=[x₁, x₂, …, x_k]^T为所有用户同时发送的K×1维符号矩阵，H∈ℂ^N×K为瑞利衰落信道矩阵，则基站端接收到的N×1维信号矢量可记为

$ \mathit{\boldsymbol{y = Hx + n}} $

(1)

通过MMSE算法检测发送信号矢量$\mathit{\boldsymbol{\hat x}}$可表示为

$ \mathit{\boldsymbol{\hat x = }}{\mathit{\boldsymbol{W}}^{ - 1}}{\mathit{\boldsymbol{y}}_{{\rm{MF}}}} $

(2)

其中y_MF=H^Hy为匹配滤波器的输出.令G=H^HH，则MMSE的滤波矩阵W^-1表示为

$ {\mathit{\boldsymbol{W}}^{ - 1}} = {\left( {\mathit{\boldsymbol{G}} + {\sigma ^2}{\mathit{\boldsymbol{I}}_K}} \right)^{ - 1}} $

(3)

在大规模MIMO系统中，当基站端天线数远大于单天线用户数(N≥K) 时，矩阵W最后趋近于对角矩阵^[3].基于此重要特征，可仅对W的对角元素求逆，来取代MMSE精确矩阵求逆，以减少随着基站天线数成立方增长的计算复杂度，即

$ {\mathit{\boldsymbol{W}}^{ - 1}} = {\left( {{\mathit{\boldsymbol{D}}^H}\mathit{\boldsymbol{D}} + {\sigma ^2}{\mathit{\boldsymbol{I}}_K}} \right)^{ - 1}} $

(4)

其中D为矩阵W的对角矩阵，但在实际系统实现中，这种粗略的近似方法会导致巨大的性能损失.

2 基于MMSE准则的线性硬检测 2.1 Neumann级数展开

为了平衡计算复杂度和检测性能，Wu等^[4]在单载波频分复用(SC-FDMA，single-carrier frequency-division multiple access) 大规模MIMO系统中，提出了Neumann级数展开，将W矩阵求逆作如下简化为

$ {\mathit{\boldsymbol{W}}^{ - 1}} = \sum\limits_{n = 0}^\infty {{{\left( {{\mathit{\boldsymbol{X}}^{ - 1}}\left( {\mathit{\boldsymbol{X}} - \mathit{\boldsymbol{W}}} \right)} \right)}^n}{\mathit{\boldsymbol{X}}^{ - 1}}} $

(5)

其中当n→∞时，(I_K-X^-1W)=0_K，级数收敛.通过分解正规化格拉姆矩阵W=D+E，其中D为矩阵W的对角矩阵，E为空心矩阵，将式(5) 中的X由D取代，展开前t项为

$ {\mathit{\boldsymbol{W}}^{ - 1}} = \sum\limits_{n = 0}^{t - 1} {{{\left( { - {\mathit{\boldsymbol{D}}^{ - 1}}\mathit{\boldsymbol{E}}} \right)}^n}{\mathit{\boldsymbol{D}}^{ - 1}}} $

(6)

其中当t→∞时，(-D^-1E)ⁿ=0_K，级数收敛.对较小的t值，Neumann级数展开能以较小的复杂度近似计算.比如，当t=1时，得到W₁^-1=D^-1，相较于MMSE精确矩阵求逆时，计算复杂度从O(K³) 减少到O(K)；当t=2时，W₂^-1=D^-1-D^-1ED^-1，其计算复杂度为O(K²).因此，可在较低计算复杂度的情况下，获得二阶的Neumann级数近似；当t=3时，W₃^-1=D^-1-D^-1ED^-1+D^-1ED^-1ED^-1，计算复杂度为O(K³)，与MMSE精确矩阵求逆的计算复杂度相同，低复杂度的优势就不那么明显了.

2.2 多级线性接收

2.1节中提到的Neumann级数展开，当展开级数t=3时，计算复杂度与MMSE精确矩阵求逆的计算复杂度相同，但检测性能尚有很大的提升空间.为了解决这一问题，下面提出了多级线性接收来近似MMSE精确矩阵的求逆.与Neumann级数展开算法一样，多级线性接收也是基于多项式展开^[5]，可通过展开较低的阶数，达到很好的检测性能.利用多级线性接收，求逆矩阵W^-1可表示为

$ {\mathit{\boldsymbol{W}}^{ - 1}} = {\left( {\mathit{\boldsymbol{G}} + {\mathit{\boldsymbol{\sigma }}^2}{\mathit{\boldsymbol{I}}_K}} \right)^{ - 1}}\sum\limits_{n = 0}^t {{\mathit{\boldsymbol{\omega }}_n}{{\left( \mathit{\boldsymbol{G}} \right)}^n}} $

(7)

其中：ω_n为第n级优化权重，t为选择的展开级数.

定义ω=[ω₀, ω₁, …, ω_n, …, ω_t]^T，则优化权矢量可计算为

$ \mathit{\boldsymbol{\omega }} = {\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}^{ - 1}}\mathit{\boldsymbol{c}} $

(8)

其中Φ为(t+1)×(t+1) 维方阵，其中元素为

$ {\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_{ij}} = {\rm{tr}}\left[ {{{\left( \mathit{\boldsymbol{G}} \right)}^{i + j + 2}}} \right] + {\mathit{\boldsymbol{\sigma }}^2}{\rm{tr}}\left[ {{{\left( \mathit{\boldsymbol{G}} \right)}^{i + j + 1}}} \right] $

(9)

c为(t+1)×1的列向量，且元素为

$ {\mathit{\boldsymbol{c}}_i} = {\rm{tr}}\left[ {{{\left( \mathit{\boldsymbol{G}} \right)}^{i + 1}}} \right] $

(10)

多级线性接收的计算复杂度由2部分构成，一部分为优化权重的计算，复杂度为(t+1)².另一部分为G的求和以及匹配滤波输出的计算，复杂度为tK².因此，总的计算复杂度为(t+1)²+tK²，与MMSE比较，复杂度由O(K³) 减少到O(K²).

2.3 Richardson算法

不同于2.1节和2.2节介绍的Neumann级数展开和多级线性接收，是基于多项式展开来近似计算求逆矩阵W^-1.下面将介绍的Richardson算法通过迭代的方法，在接收端估计发送信号矢量$\mathit{\boldsymbol{\hat x}}$时，直接计算W^-1y_MF.式(2) 也可表示为

$ \mathit{\boldsymbol{W\hat x}} = {\mathit{\boldsymbol{y}}_{{\rm{MF}}}} $

(11)

于是把简化矩阵求逆问题转换为求解Rs=b线性方程的问题.其中R为对称正定矩阵，b为N×1的测量向量，s为解向量.根据文献[6]，Richardson算法可描述为

$ {\mathit{\boldsymbol{s}}^{\left( {t + 1} \right)}} = {\mathit{\boldsymbol{s}}^{\left( t \right)}} + \mathit{\boldsymbol{\omega }}\left( {\mathit{\boldsymbol{b}} - \mathit{\boldsymbol{R}}{\mathit{\boldsymbol{s}}^{\left( t \right)}}} \right) $

(12)

其中：s^(t)为通过Richardson算法t次迭代得到的线性方程的解向量，ω为松弛因子，控制Richardson算法收敛的快慢.

由于在大规模MIMO系统中，当基站端的天线数远大于单天线用户数时，W=G+σ²I_K为对称正定的^[6].利用Richardson公式估计发送信号矢量$\mathit{\boldsymbol{\hat x}}$为

$ {{\mathit{\boldsymbol{\hat x}}}^{\left( {t + 1} \right)}} = {{\mathit{\boldsymbol{\hat x}}}^{\left( t \right)}} + \mathit{\boldsymbol{\omega }}\left( {{\mathit{\boldsymbol{y}}_{{\rm{MF}}}} - \mathit{\boldsymbol{W}}{{\mathit{\boldsymbol{\hat x}}}^{\left( t \right)}}} \right) $

(13)

其中$\mathit{\boldsymbol{\hat x}}$ ^(t)为通过Richardson算法t次迭代估计出的发送信号矢量.松弛因子ω一般选择为0 < ω < 2/λ₁，λ₁为W的最大特征值.

Richardson算法的计算复杂度由2部分构成，一部分为K×K的矩阵W和K×1的矢量$\mathit{\boldsymbol{\hat x}}$ ^(t)相乘，另一部分为松弛参数ω和K×1的矢量$\mathit{\boldsymbol{\hat y-W}}{{\mathit{\boldsymbol{\hat x}}}^{\left( t \right)}}$相乘.因此，每次迭代总的计算复杂度为K²+K，与MMSE比较，复杂度由O(K³) 降为了O(K²).

2.4 松弛迭代算法

下面介绍的松弛迭代算法^[7](RI，relaxation iterative) 和2.3节介绍的Richardson算法都是基于迭代方法来取代复杂的矩阵求逆.但松弛迭代算法能以更少的迭代次数达到接近MMSE的检测性能.

如2.3节式(13) 所示，把矩阵求逆的问题转化成了求解Rs=b线性方程的问题.将K×K维对称正定矩阵分解为R=D-L-U，其中D、-L和-U分别为R的对角矩阵、严格下三角矩阵和严格上三角矩阵.松弛迭代算法可描述为

$ {\left( {\mathit{\boldsymbol{D - }}\frac{1}{\mathit{\boldsymbol{\omega }}}{\mathit{\boldsymbol{L}}^{\rm{H}}}} \right)^{ - 1}}\mathit{\boldsymbol{b}} $

(14)

其中s^(t)为通过松弛迭代算法t次迭代得到的线性方程的解向量.

利用上述的松弛迭代算法求解$\mathit{\boldsymbol{W\hat x}} = {\mathit{\boldsymbol{y}}_{{\rm{MF}}}}$方程，来估计发送信号矢量$\mathit{\boldsymbol{\hat x}}$.对称正定矩阵W可分解为W=D-L-L^H，其中D、-L和-L^H分别为W的对角矩阵、严格下三角矩阵和严格上三角矩阵.因此，通过松弛迭代算法t次迭代估计得到发送信号矢量为

$ {{\mathit{\boldsymbol{\hat x}}}^{\left( {t + 1} \right)}} = \mathit{\boldsymbol{B}}{{\mathit{\boldsymbol{\hat x}}}^{\left( t \right)}} + \mathit{\boldsymbol{C}} $

(15)

其中：B=(ωD-L^H)^-1[(ω-1)L^H+ωL]，C=(D-1/ωL^H)^-1y_MF，松弛参数一般选为ω≥1/2.

由式(15) 可知，松弛迭代算法每次迭代需要总的计算复杂度为K²+2K.与MMSE比较，复杂度由O(K³) 降为了O(K²).尽管相对于Richardson算法，松弛迭代算法每次迭代的计算复杂度更高，但是它能以较少的迭代次数获得更优的检测性能.各简化算法的检测性能和计算复杂度对比分析如表 1所示.

3 基于MMSE准则的线性软检测

前面介绍的各种矩阵求逆算法都是以硬判决为核心的检测技术.为了充分利用信道编译码的软信息，减少硬判决的信息损失，笔者提出了将各求逆简化算法应用于软判决中，进一步提高检测性能.

下面首先介绍基于多项式展开的矩阵求逆简化算法如何通过信号检测器产生软输出值.令A_i为A=W^-1H^H的第i行向量，其中W^-1为式(6) 和式(7) 中Neumann级数展开和多级线性接收的近似逆矩阵，则2种算法下对应的A_i分别可表示为

$ {\mathit{\boldsymbol{A}}_i} = {\left[ {\sum\limits_{n = 0}^{t - 1} {{{\left( { - {\mathit{\boldsymbol{D}}^{ - 1}}\mathit{\boldsymbol{E}}} \right)}^n}{\mathit{\boldsymbol{D}}^{ - 1}}} {\mathit{\boldsymbol{H}}^{\rm{H}}}} \right]_i} $

(16)

$ {\mathit{\boldsymbol{A}}_i} = {\left[ {\sum\limits_{n = 0}^t {{\mathit{\boldsymbol{\omega }}_n}{{\left( \mathit{\boldsymbol{G}} \right)}^n}} {\mathit{\boldsymbol{H}}^{\rm{H}}}} \right]_i} $

(17)

针对{x_i}_i=1^K，简化后的MMSE检测器的输出为

$ \begin{array}{*{20}{c}} {{{\mathit{\hat x}}_i} = {\mathit{\boldsymbol{A}}_i}\mathit{\boldsymbol{y = }}{\mathit{\boldsymbol{A}}_i}{\mathit{\boldsymbol{h}}_1}{x_1} + \cdots + {\mathit{\boldsymbol{A}}_i}{\mathit{\boldsymbol{h}}_K}{x_K} + {\mathit{\boldsymbol{A}}_i}\mathit{\boldsymbol{z = }}}\\ {{\rho _i}{x_i} + {\mathit{\boldsymbol{I}}_j} + \tilde z,j \ne i} \end{array} $

(18)

其中：MMSE简化算法均衡后得到的纯信号项为ρ_ix_i=A_ih_ix_i，MMSE简化算法均衡后得到的来自其他天线的残留干扰项为${\mathit{\boldsymbol{I}}_j} = \sum\limits_{j = 1, j \ne i}^K {{\mathit{\boldsymbol{A}}_i}{\mathit{\boldsymbol{h}}_j}{x_j}} $，MMSE简化算法均衡后得到的纯噪声项为$\tilde z = {A_i}\mathit{\boldsymbol{z}}$.根据式(18)，可获得第i层上的信号与干扰和噪声功率之间的信号与干扰加噪声比(SINR，signal to interference plus noise ratio) 为

$ \begin{array}{*{20}{c}} {{\delta _i} = \frac{{E{{\left\{ {\left| {{\rho _i}{x_i}} \right|} \right\}}^2}}}{{E\left\{ {{{\left| {{I_j}} \right|}^2}} \right\} + E{{\left\{ {\left| {\tilde z} \right|} \right\}}^2}}} = }\\ {\frac{{{{\left| {{\rho _i}} \right|}^2}{E_x}}}{{{{\left| {{\mathit{\boldsymbol{A}}_i}{\mathit{\boldsymbol{h}}_j}} \right|}^2}{E_x} + {{\left\| {{\mathit{\boldsymbol{A}}_i}} \right\|}^2}{\sigma ^2}}},j \ne i} \end{array} $

(19)

假设式(19) 中的干扰和噪声项(NPI，noise-plus-interference) 分量是独立和高斯分布的，则可用1个均值为0、方差为σ_i²=|A_ih_j|²E_x+‖A_i‖²σ²的高斯随机变量近似表示为${\mathit{\boldsymbol{I}}_j} + \tilde z$.

由于基于迭代方法的Richardson算法和松弛迭代算法直接估计发送矢量$\mathit{\boldsymbol{\hat x}}$，避免了对MMSE精确的矩阵求逆操作，无法直接得到W^-1，而W^-1矩阵对于软检测中的NPI评估不可或缺.针对以上问题，Dai等^[8]提出了一种近似估计求逆矩阵的方法.分别组合式(3)、式(13) 和式(15)，且令y_MF=I_K，利用Richardson算法和松弛迭代算法分别估计MMSE求逆矩阵${{\mathit{\boldsymbol{\hat W}}}_{{\rm{inv}}}}$为

$ {\left( {{{\mathit{\boldsymbol{\hat W}}}_{{\rm{inv}}}}} \right)^{\left( {t + 1} \right)}} = {\left( {{{\mathit{\boldsymbol{\hat W}}}_{{\rm{inv}}}}} \right)^{\left( t \right)}} + \mathit{\boldsymbol{\omega }}\left[ {{\mathit{\boldsymbol{I}}_K} - \mathit{\boldsymbol{W}}{{\left( {{{\mathit{\boldsymbol{\hat W}}}_{{\rm{inv}}}}} \right)}^{\left( t \right)}}} \right] $

(20)

$ {\left( {{{\mathit{\boldsymbol{\hat W}}}_{{\rm{inv}}}}} \right)^{\left( {t + 1} \right)}} = \mathit{\boldsymbol{B}}{\left( {{{\mathit{\boldsymbol{\hat W}}}_{{\rm{inv}}}}} \right)^{\left( t \right)}} + \mathit{\boldsymbol{C}}{\mathit{\boldsymbol{I}}_K} $

(21)

其中：B=(ωD-L^H)^-1[(ω-1)L^H+ωL]，C=(D-1/ωL^H)^-1I_K.若采用对角近似作为估计求逆矩阵的初始值，${\left( {{{\mathit{\boldsymbol{\hat W}}}_{{\rm{inv}}}}} \right)^{\left( 0 \right)}} = {\mathit{\boldsymbol{D}}^{-1}}{\mathit{\boldsymbol{I}}_K}$，则近似NPI的方差可表示为

$ \sigma _i^2 = \sum\limits_{i \ne j}^K {{{\left| {{\mathit{\boldsymbol{V}}_{ij}}} \right|}^2} + {\mathit{\boldsymbol{U}}_{ij}}{\sigma ^2}} $

(22)

其中：V=D^-1G，U=D^-1GD^-1.且有ρ_ix_i=V_iix_i.因此，可得下面的条件概率密度函数(PDF，probability density function) 为

$ {f_{\hat x}}\left( {{{\hat x}_i}\left| {{x_i}} \right.} \right) = \frac{1}{{\sqrt {2\pi \sigma _i^2} }}\exp \left( { - \frac{{{{\left| {{{\hat x}_i} - {\rho _i}{x_i}} \right|}^2}}}{{2\sigma _i^2}}} \right) $

(23)

令S_{l, i}⁺和S_{l, i}^-分别为第i个符号的第l位取值为1和0时对应的调制星座集，x_i的第l位的对数似然比(LLR，log likelihood ratio) 表示为

$ \begin{array}{*{20}{c}} {\eta \left( {{b_{l,i}}} \right) = \ln \frac{{\sum\limits_{x \in S_{l,i}^ + } {{f_{\hat x}}\left( {{{\hat x}_i}\left| x \right.} \right)} }}{{\sum\limits_{x \in S_{l,i}^ - } {{f_{\hat x}}\left( {{{\hat x}_i}\left| x \right.} \right)} }} \approx \ln \frac{{\mathop {\max }\limits_{x \in S_{l,i}^ + } {f_{\hat x}}\left( {{{\hat x}_i}\left| x \right.} \right)}}{{\mathop {\max }\limits_{x \in S_{l,i}^ - } {f_{\hat x}}\left( {{{\hat x}_i}\left| x \right.} \right)}} = }\\ {\frac{1}{{2\sigma _i^2}}\left( {{{\left| {{{\hat x}_i} - x_{i,j,{\rm{opt}}}^ - } \right|}^2} - {{\left| {{{\hat x}_i} - x_{i,j,{\rm{opt}}}^ + } \right|}^2}} \right)} \end{array} $

(24)

其中：x_{i, l, opt}⁺和x_{i, l, opt}^-定义为

$ \begin{array}{*{20}{c}} {x_{i,j,{\rm{opt}}}^ + = \mathop {\arg \min }\limits_{x \in S_l^ + } {{\left| {{{\hat x}_i} - x} \right|}^2}}\\ {x_{i,j,{\rm{opt}}}^ - = \mathop {\arg \min }\limits_{x \in S_l^ - } {{\left| {{{\hat x}_i} - x} \right|}^2}} \end{array} $

(25)

4 仿真结果

为了验证各种矩阵求逆简化算法的检测性能，下面对比分析了Neumann级数展开、多级线性接收、Richardson算法和松弛迭代算法的误比特率(BER，bit error rate)，并以仅对对角矩阵求逆和MMSE算法的检测性能作为参照.设置仿真时的传输信道为准静态瑞利衰落信道，基带信号调制方式为16-正交振幅调制(QAM，quadrature amplitude modulation)，天线规模N×K为128×16.在仿真中，t为多项式展开的项数和迭代的次数.

图 1显示了各简化算法在硬判决中的BER检测性能对比.随着展开级数或迭代次数的不断增加，各简化算法的检测性能也随之提高.在相同的迭代次数下，松弛迭代算法检测性能最好，多级线性接收次之，Richardson算法和Neumann级数展开的检测性能较差.比如在迭代次数t=2，BER为10^-4时，松弛迭代算法需要的信噪比为10 dB，多级线性接收所需为11 dB，而Richardson算法和Neumann级数展开需要超过16 dB.同时，松弛迭代算法和多级线性接收能通过更少的迭代次数，达到接近MMSE算法的检测性能.

图 2显示了各简化算法在软判决中的BER检测性能对比.很明显，在软判决中，各种简化算法的检测性能都有明显地提升.在信噪比为2 dB时，松弛迭代算法和多级线性接收展开级数t=3时的误码率就接近10^-5.而在硬判决中，BER为10^-5时，需要11 dB的信噪比.容易看出，在软判决中，松弛迭代算法和多级线性接收在检测性能上依然优于Richardson算法和Neumann级数展开.但是，在硬判决中性能相对一般的Neumann级数展开在软判决中能通过较少的展开级数即能达到接近MMSE最优的检测性能，而性能较好的Richardson算法在多次迭代后性能依然不理想.且在软判决中，多级线性接收优于松弛迭代算法.

5 结束语

MMSE算法已广泛运用于MIMO通信系统中，但当其涉及大规模矩阵求逆时，计算复杂度太高.在矩阵求逆简化算法中，基于多项式展开的Neumann级数展开和多级线性接收，以及基于迭代方法的Richardson算法和松弛迭代算法，代表了2类典型的以低复杂度为实现目标的解决方案.对比分析2类算法在硬判决方案中的检测性能和计算复杂度之间的关系，在计算复杂度为O(K²) 时，松弛迭代算法的检测性能最优，多级线性接收次之，Richardson算法和Neumann级数展开相对较差.且若将以上简化算法应用于软判决方案中，能进一步使得2类简化算法的检测性能都有明显的提升.与硬判决方案不同的是，在软判决中，多级线性接收略优于松弛迭代算法，Neumann级数展开明显优于Richardson算法.基于多项式展开和迭代方法的这2类算法由于较低的实现复杂度，可作为大规模MIMO系统信号检测实现中的现实解决方案.