Millimeter-wave Channel Estimation with Intelligent Reflecting Surface Assisted Based on Vector Approximate Message Passing

1.   引言

顺应5G无线通信的新时代，对通信系统提出了更高速率、更大容量、更低时延的新要求。随着Sub-6 GHz频段的商用，毫米波也逐渐引起关注并投入研究，为在广域级和局域级提供无缝连接，实现毫米波频谱^[1]的利用至关重要。针对5G毫米波传播损耗高、绕射和衍射能力弱、覆盖相对受限等挑战，智能反射面(Intelligent Reflecting Surface, IRS)技术作为一种新技术，通过低功耗甚至无功耗反射电磁波，降低信号损耗，大大地提高了毫米波通信的信道质量，以此更加适用于非视距传播(non-Line-of-Sight, nLoS)环境中。由于IRS低功耗的特点，使得其在很大程度上优于中继设备。IRS不受接收噪声的影响^[2]、具有全频带响应、部署灵活性和兼容性^[3]等众多优点。将IRS与无人机相结合辅助通信将大大地提高通信质量和速率。针对更多基站的部署将增大电磁辐射水平的问题，IRS的智能性使得部署IRS设备将降低电磁辐射水平。同时IRS也可通过调制来实现低复杂度的大规模多输入多输出(Multiple Input Multiple Output, MIMO)技术。结合IRS的诸多优点，使得该技术在5G毫米波以及未来的6G时代获得越来越多的关注和研究。

对于IRS辅助通信的信道估计问题，文献[4]采用行列块稀疏度来联合估计级联信道，首先利用公共列块稀疏性来估计由出射角阵列控制向量跨越的公共子空间，再利用块稀疏性来制定基于多个测量向量的多用户联合稀疏矩阵恢复问题。但通常优化变量在非凸且难以求解的模拟多用户联合稀疏矩阵恢复问题中是耦合的，此时可通过一种基于替代优化^[4]和迭代加权算法^[4]的方法来有效地解决问题。文献[5]提出了一种基于双线性广义近似消息传递(Bilinear Generalized Approximate Message Passing, BiG-AMP)的稀疏矩阵分解算法来解决信道的稀疏问题，并且利用黎曼梯度算法实现稀疏矩阵补全，此算法的性能会随着在稀疏矩阵分解阶段估计的随机变量的增加而变差。对比于BiG-AMP算法，有学者提出双线性自适应矢量近似消息传递算法^[6](Bilinear Adaptive Vector Approximate Message Passing, BAdVAMP)来解决IRS信道估计，此算法的性能优于前者。在对稀疏信道进行信道估计时，一般可运用压缩感知算法^[7,8]。文献[9]通过压缩感知和深度学习对连接到控制器的IRS单元进行控制，从而估计所有经过IRS单元上的信道。文献[10]中利用压缩感知算法的同时，结合AMP算法对IRS辅助的级联信道进行信道估计。在IRS辅助通信系统中，优化线性预编码算法^[11]和迭代算法^[12]同样也是解决信道估计的方法。在毫米波通信系统中，有研究者提出在IRS辅助毫米波网络中使用波束搜索^[13]方法找到所需用户的最佳波束，以增强毫米波的信道质量。但对于以上所有算法都没有针对IRS辅助毫米波通信系统的信道估计问题进行探索，本文将对此系统环境的信道估计问题进行研究。

针对IRS系统中级联信道的复杂性以及IRS辅助毫米波通信的信道稀疏性，本文结合Khatri-Rao积^[14]、克罗内克积和矢量近似消息传递算法(Khatri-Rao and Vector Approximate Message Passing, KR-VAMP)提出了KR-VAMP算法来解决系统的信道估计问题。由于加入IRS使得系统信道为级联信道，增加了信道的复杂性，本文通过Khatri-Rao积和克罗内克积^[15]将级联信道${\boldsymbol{H}}$转换为稀疏信号恢复问题。VAMP结合了近似消息传递算法^[16](Approximate Message Passing, AMP)的高效的迭代性质以及因子传递的状态演化方式，提供了准确的信道估计值，并且降低了训练迭代次数，为解决IRS辅助毫米波系统的信道估计问题提供了良好的解决方案。

全文结构安排如下：第2节介绍IRS辅助毫米波系统的信道模型，以及AMP算法模型。第3节对KR-VAMP算法进行分析和推导。第4节通过KR-VAMP算法和其他算法进行对比，通过仿真验证KR-VAMP算法的有效性。第5节对全文进行总结并给出下一步工作方向。

符号说明：符号${\left( \cdot \right)^{\text{T}}}$, ${\left( \cdot \right)^ * }$, ${\left( \cdot \right)^{\text{H}}}$以及${\text{vec}}\left( \cdot \right)$分别表示转置、共轭、共轭转置和向量化运算符。${{\hat {\boldsymbol{A}}}}$, ${{\boldsymbol{A}}'}$和${{\boldsymbol{A}}^ + }$分别代表矩阵${\boldsymbol{A}}$的估计值、求导以及左逆表达式，${{\boldsymbol{A}}^ + }$可由${\left( {{{\boldsymbol{A}}^{\text{T}}}{\boldsymbol{A}}} \right)^{ - 1}}{{\boldsymbol{A}}^{\text{T}}}$计算得出。$\otimes$和$\odot$分别为克罗内克积和Khatri-Rao积。

2.   IRS辅助通信的系统模型

为增强毫米波通信的信道质量，引入IRS辅助毫米波通信。如图1所示的IRS辅助通信的信道示意图，当基站与用户之间的视距(Line-of-Sight, LoS)传播受阻时，通过IRS辅助毫米波信号的nLoS信道传播。基站通过微控制器^[17]对IRS系统的每个单元进行控制调节，本文主要研究IRS辅助通信的下行链路系统。在此通信系统中，有$M$条发送天线，IRS由$N$个单元组成，接收端为$K$个用户。

图  1  IRS辅助通信示意图

下载: 全尺寸图片 幻灯片

为简化模型，假设系统为单用户系统。此系统以毫米波信号为例，基站与IRS之间的信道${{\boldsymbol{h}}_{{\text{IB}}}} \in {\mathbb{C}^{N \times M}}$，IRS与用户之间的信道为${{\boldsymbol{h}}_{{\text{UI}}}} \in {\mathbb{C}^{1 \times N}}$，且基站与用户之间受干扰物阻挡，无法进行LoS通信。在本文中，利用酉矩阵^[18]构建IRS，其中IRS系统的反射系数为

其中，${\beta _n} \in \left[ {0,1} \right]$为反射幅度，${\theta _n} \in \left[ {0,}\;\; {2\pi } \right)$为反射相位。根据毫米波的特性以及基本信道模型，假设基站与IRS之间有${L_1}$条路径，IRS与用户之间有${L_2}$条路径，路径的衰落系数由$\partial$表示，假设环境中的电磁波传播路径较为稳定，$\partial$不随时间变化，$\partial$根据路径距离的不同而不同，发射端的方向矢量和接收端的方向矢量分别用${{\boldsymbol{a}}_{\text{t}}}$和${{\boldsymbol{a}}_{\text{r}}}$表示，则可将信道${{\boldsymbol{h}}_{{\text{IB}}}}$和信道${{\boldsymbol{h}}_{{\text{UI}}}}$分别表示为

基站与IRS之间的位置关系图如图2(a)所示，假设${{\boldsymbol{f}}_{dN}}\left( x \right) \triangleq \dfrac{1}{N}{\left[ {1, \cdots ,{{\text{e}}^{{\text{j}}\frac{{2\pi }}{\lambda }d\left( {nx} \right)}}, \cdots ,{{\text{e}}^{{\text{j}}\frac{{2\pi }}{\lambda }d\left( {\left( {N - 1} \right)x} \right)}}} \right]^{\text{T}}}$，即在式(3)中${{\boldsymbol{a}}_{\text{t}}}\left( \phi \right) = {{\boldsymbol{f}}_{dM}}\left( \phi \right)$，其中，基站处的出发角$\phi = \cos \theta \cos \varphi$。${{\boldsymbol{a}}_{\text{r}}}\left( {{\vartheta _l},{\omega _l}} \right) = {{\boldsymbol{a}}_{{y}}}\left( v \right) \otimes {{\boldsymbol{a}}_{{z}}}\left( w \right)$，此时$\vartheta$表示IRS处到达角，$\omega$为仰角，其中$v = \cos \theta \cos \varphi$, $w = \sin \theta$, ${{\boldsymbol{a}}_{{y}}}\left( v \right)$和${{\boldsymbol{a}}_{{z}}}\left( w \right)$可分别表示为${{\boldsymbol{a}}_{{y}}}\left( v \right){\text{ = }} {{\boldsymbol{f}}_{p{N_1}}}\left( v \right)$, ${{\boldsymbol{a}}_{{z}}}\left( w \right){\text{ = }}{{\boldsymbol{f}}_{p{N_2}}}\left( w \right)$，其中${N_1}$和${N_2}$分别为IRS横纵元素数量，且满足${N_1} \times {N_2} = N$，则方向矢量${{\boldsymbol{a}}_{\text{r}}}\left( {{\vartheta _l},{\omega _l}} \right)$可表示为

图  2  基站、IRS与用户之间的位置关系图

下载: 全尺寸图片 幻灯片

IRS与用户之间的位置关系图如图2(b)所示，即在式(4)中，${{\boldsymbol{a}}_{\text{t}}}\left( {{\vartheta _l},{\omega _l}} \right) = {{\boldsymbol{a}}_{{y}}}\left( v \right) \otimes {{\boldsymbol{a}}_{{z}}}\left( w \right)$，此时$\vartheta$表示IRS处出发角，其中$v = \sin \theta \cos \varphi$, $w = \cos \theta$。则方向矢量${{\boldsymbol{a}}_{\text{t}}}\left( {{\vartheta _l},{\omega _l}} \right)$同理于式(5)。

由于毫米波信道的稀疏散射特性，路径数量远小于${{\boldsymbol{h}}_{{\text{IB}}}}$或${{\boldsymbol{h}}_{{\text{UI}}}}$的维度，即${{\boldsymbol{h}}_{{\text{IB}}}}$和${{\boldsymbol{h}}_{{\text{UI}}}}$可构建完备矩阵，分别表示为

其中，${{\boldsymbol{A}}_{\text{I}}}$和${{\boldsymbol{A}}_{\text{B}}}$可分别表示为

式(6)和式(7)所示的稀疏信道的信道系数${\boldsymbol{\varPsi}} \in {\mathbb{C}^{{N_G} \times {M_G}}}$, ${\boldsymbol{\varGamma}} \in {\mathbb{C}^{{N_G} \times 1}}$，其中${N_{{\text{G}},1}} \times {N_{{\text{G}},2}} = {N_{\text{G}}}$，并且满足${N_{\text{G}}} \gg N$, ${M_{\text{G}}} \gg M$。

在如图1所示的IRS辅助毫米波通信的级联信道${\boldsymbol{H}}$可表示为

其中，${\text{diag}}\left( {{{\boldsymbol{h}}_{{\text{UI}}}}} \right){{\boldsymbol{h}}_{{\text{IB}}}}$用${\boldsymbol{G}}$表示，结合Khatri-Rao积^[14]以及克罗内克积^[15]的特性可将${\boldsymbol{G}}$表示为

其中，${\boldsymbol{D}}$表示${\left( {{\boldsymbol{A}}_{\text{I}}^{\text{H}} \odot {\boldsymbol{A}}_{\text{I}}^{\text{T}}} \right)^{\text{T}}}$，即级联信道${\boldsymbol{H}}$可用${{\boldsymbol{V}}^{\text{T}}} \cdot {\boldsymbol{D}} \cdot \left( {{{\boldsymbol{\varGamma}} ^ * } \otimes {\boldsymbol{\varPsi}} } \right) \cdot {\boldsymbol{A}}_{\text{B}}^{\text{H}}$表示。假定基站的发射信号为${\boldsymbol{X}}$，系统信道存在均值为0， 方差为$\xi _\varSigma ^{ - 1}$的噪声，噪声信号用${\boldsymbol{\varSigma}}$表示。结合式(11)可将用户接收信号表示为

3.   AMP算法分析

在非凸实例中，AMP不仅可以使用可测试的贝叶斯最优条件精确地描述方法的性能，而且承接了迭代阈值算法的优势，也融入了近似消息传递算法的精髓。本文将结合AMP算法解决信道${\boldsymbol{H}}$的估计问题。假设${{\boldsymbol{H}}^0}$为无噪声条件下的信道冲激响应，即在有噪声环境中，接收信号为${\boldsymbol{Y}} = {{\boldsymbol{H}}^0}{\boldsymbol{X}} + {\boldsymbol{\varSigma}}$，其中${\boldsymbol{\varSigma}}$是噪声信号，且服从均值为0，方差为$\xi _\varSigma ^{ - 1}$的高斯分布。根据AMP算法的特性和迭代阈值原理可将残差函数表示为

其中，$k$为函数迭代次数，伪数据${{\boldsymbol{r}}_k} = {{\boldsymbol{\hat H}}_k} + {{\boldsymbol{X}}^{\text{T}}}{{\boldsymbol{R}}_k}$为信道估计值与响应残差的和，${\gamma _k}$为阈值函数中的阈值精度，本文设定其迭代规律如式(16)所示，${{\boldsymbol{\hat H}}_k} = g\left( {{{\boldsymbol{r}}_{k - 1}},{\gamma _{k - 1}}} \right)$为信道的估计值，${\alpha _k} = \left\langle {g'}\left( {{\boldsymbol{r}}_{k - 1}}, {\gamma _{k - 1}} \right) \right\rangle$可表示为

其中，${M'} = \left\| {\boldsymbol{H}} \right\|_{\text{F}}^{\text{2}}$，阈值函数$g$可由最大后验估计(Maximum A Posteriori, MAP)确定，假设${\boldsymbol{H}}$的先验概率为$p\left( {\boldsymbol{H}} \right) = \displaystyle\prod\nolimits_{m = 1}^M {p\left( {{H_m}} \right)}$，利用MAP可将阈值函数$g$表示为

由于信道估计值由真实信道冲激响应与信道噪声组成，可假设${{\boldsymbol{r}}_k}$表示为真实信道冲激响应和干扰的和，即${{\boldsymbol{r}}_k} = {{\boldsymbol{H}}^0} + N\left( {0,{\tau _k}{\boldsymbol{I}}} \right)$，${\tau _k}$为干扰信号的平均功率。即${{\boldsymbol{\hat H}}_k}$的均方误差可表示为

设定初始条件$\varepsilon \left( {{\tau _{ - 1}},{\gamma _{ - 1}}} \right) = \mu$，随着初始值$\mu$的变化，${{\boldsymbol{\hat H}}_k}$的MMSE也将产生规律性的变化，第4节将给出具体仿真验证结果。式(16)中${\tau _k}$和${\gamma _k}$满足的迭代条件为

结合传统AMP算法分析和描述，AMP算法的仿真流程如下。

步骤1　设置迭代总次数$K$以及初始$g\left( {{{\boldsymbol{r}}_0},{\gamma _0}} \right)$，${{\boldsymbol{R}}_0} = 0$，迭代次数初始值设定为$k = 1$。

步骤2　${{\boldsymbol{\hat H}}_k} = g\left( {{{\boldsymbol{r}}_{k - 1}},{\gamma _{k - 1}}} \right)$，函数$g$由MAP确定。

步骤3　${\alpha _k} = \left\langle {{g'}\left( {{{\boldsymbol{r}}_{k - 1}},{\gamma _{k - 1}}} \right)} \right\rangle$。

步骤4　${{\boldsymbol{R}}_k} = {\boldsymbol{Y}} - {{\boldsymbol{\hat H}}_k}{\boldsymbol{X}} + \dfrac{1}{\delta }{{\boldsymbol{R}}_{k - 1}}\left\langle {g'}\left( {{\boldsymbol{r}}_{k - 1}}, {\gamma _{k - 1}} \right) \right\rangle$。

步骤5　${{\boldsymbol{r}}_k} = {{\boldsymbol{\hat H}}_k} + {{\boldsymbol{X}}^{\text{T}}}{{\boldsymbol{R}}_k}$。

步骤6　确定${\gamma _k}$值。

步骤7　将迭代次数$k + 1$，直至迭代次数等于$K + 1$，得到${{\boldsymbol{\hat H}}_K}$，否则循环步骤2—步骤7。

4.   KR-VAMP算法分析

本节将对KR-VAMP 算法的因子转换方式进行分析。首先在式(12)中，定义${\boldsymbol{\varPhi}} \triangleq {{\boldsymbol{\varGamma}} ^ * } \otimes {\boldsymbol{\varPsi}}$, ${\boldsymbol{\varOmega}} \triangleq \left( {{{\boldsymbol{X}}^{\text{T}}} \otimes {{\boldsymbol{V}}^{\text{T}}}} \right) \cdot \left( {{\boldsymbol{A}}_{\text{B}}^ * \otimes {\boldsymbol{D}}} \right)$，假设在加性高斯白噪声(Additive White Gaussian Noise, AWGN)的环境中，区别于AMP, VAMP算法使用带有缩放值节点的循环因子图^[19]，并且使用带有矢量节点的非循环图，因此命名为“矢量AMP”，即VAMP。为实现VAMP，定义起始因子为

因子转化图如图3所示，黑色方框表示因子节点，白色圆形表示变量节点，通过图3因子图可知$p\left( {{\boldsymbol{\varPhi}} |{\boldsymbol{\varGamma}} ,{\boldsymbol{\varPsi}} } \right) = \delta \left( {{\boldsymbol{\varPhi}} - {\boldsymbol{\varGamma }} \otimes {\boldsymbol{\varPsi}} } \right)$，其中$\delta$为狄拉克三角分布函数，可将等效因子图表达式为

图  3  因子转换图

下载: 全尺寸图片 幻灯片

在因子转换图中，设定变量节点${\boldsymbol{\varGamma}}$处的近似置信${b_{{\text{app}}}}\left( {\boldsymbol{\varGamma}} \right) = N\left( {{\boldsymbol{\varGamma}} ;{{\hat {{\varGamma}} }_k},\eta _{{{\varGamma}} k}^{ - 1}{\boldsymbol{I}}} \right)$，其中${\hat \varGamma _k} = {\text{E}} \left[ {\varGamma |{b_{{\text{sp}}}}\left( \varGamma \right)} \right]$, $\eta _{\varGamma k}^{ - 1} = \left\langle {{\text{diag}}\left( {{{\rm{cov}}} \left[ {{\boldsymbol{\varGamma}} |{b_{{\text{sp}}}}\left( {\boldsymbol{\varGamma}} \right)} \right]} \right)} \right\rangle$，积和置信度${b_{{\text{sp}}}}\left( {\boldsymbol{\varGamma}} \right)$正比于$p\left( {\boldsymbol{\varGamma}} \right)N\left( {{\boldsymbol{\varGamma}} ;{{\boldsymbol{r}}_{\varGamma k}},\xi _{\varGamma k}^{ - 1}{\boldsymbol{I}}} \right)$。初始${u_{\delta \to {\boldsymbol{\varGamma}} }}\left( {\boldsymbol{\varGamma}} \right) = N \left( {{\boldsymbol{\varGamma}} ;{{\boldsymbol{r}}_0},\xi _0^{ - 1}{\boldsymbol{I}}} \right)$，即从变量节点${\boldsymbol{\varGamma}}$处到因子节点$\delta$处的转换关系可表示为

同理设定变量节点${\boldsymbol{\varPsi}}$处的近似置信${b_{{\text{app}}}}\left( {\boldsymbol{\varPsi}} \right) = N\left( {{\boldsymbol{\varPsi}} ;{{\hat {\boldsymbol{\varPsi}} }_k},\eta _{\varPsi k}^{ - 1}{\boldsymbol{I}}} \right)$，其中${\hat {\boldsymbol{\varPsi}} _k} = {\text{E}}\left[ {{\boldsymbol{\varPsi}} |{b_{{\text{sp}}}}\left( {\boldsymbol{\varPsi}} \right)} \right]$，$\eta _{\varPsi k}^{ - 1} = \left\langle {{\text{diag}}\left( {{{\rm{cov}}} \left[ {{{\boldsymbol{\varPsi}} }|{b_{{\text{sp}}}}\left( {{\boldsymbol{\varPsi}} } \right)} \right]} \right)} \right\rangle$，积和置信度${b_{{\text{sp}}}}\left( {\boldsymbol{\varPsi}} \right)$同理正比于$p\left( {\boldsymbol{\varPsi}} \right)N\left( {{{\boldsymbol{\varPsi}} };{{\boldsymbol{r}}_{{\varPsi }k}},\xi _{{\varPsi }k}^{ - {\text{1}}}{\boldsymbol{I}}} \right)$，并且此因子图中设定初始${u_{\delta \to {{\varPsi} }}}\left( {{\boldsymbol{\varPsi}} } \right) = N\left( {{{\boldsymbol{\varPsi}} };{{\boldsymbol{r}}_0},\xi _0^{ - 1}{\boldsymbol{I}}} \right)$，即从变量节点$\varPsi$到因子节点$\delta$的转换关系为

在图3所示的因子图中，可将${u_{\delta \to {\boldsymbol{\varPhi}} }}\left( {\boldsymbol{\varPhi}} \right)$表示为$N\left( {{{\boldsymbol{\varPhi }}};{{\boldsymbol{r}}_{{{\boldsymbol{\varPhi}} }k}},\xi _{{{\boldsymbol{\varPhi}} }k}^{ - 1}{\boldsymbol{I}}} \right)$，结合式(21)和式(22)，${{\boldsymbol{r}}_{{\boldsymbol{\varPhi}} k}}$和$\xi _{{\boldsymbol{\varPhi}} k}^{ - 1}$的表达式分别为

变量节点$\varPhi$处的近似置信${b_{{\text{app}}}}\left( {{\boldsymbol{\varPhi}} } \right) = N \left( {{\boldsymbol{\varPhi}} }; {{{\hat {\boldsymbol{\varPhi}} }}_k},\eta _{{{\boldsymbol{\varPhi}} }k}^{ - 1}{\boldsymbol{I}} \right)$，其中${{\hat {\boldsymbol{\varPhi}} }_k} = {\text{E}}\left[ {{{\boldsymbol{\varPhi}} }|{b_{{\text{sp}}}}\left( {{\boldsymbol{\varPhi}} } \right)} \right]$, $\eta _{{{\boldsymbol{\varPhi}} }k}^{ - 1} = \left\langle {\text{diag}}\left( {{\rm{cov}}} \left[ {{{\boldsymbol{\varPhi}} }|{b_{{\text{sp}}}}\left( {{\boldsymbol{\varPhi }}} \right)} \right] \right) \right\rangle$, ${\hat {\boldsymbol{\varPhi}} _k}$表达式中的${b_{{\text{sp}}}}\left( {\boldsymbol{\varPhi}} \right)$正比于$N\left( {{\boldsymbol{Y}};{{\boldsymbol{\varOmega }}}{\text{vec}}\left( {{\boldsymbol{\varPhi}} } \right),\xi _{{\boldsymbol{\varSigma}} }^{ - 1}{\boldsymbol{I}}} \right)N\left( {{{\boldsymbol{\varPhi}} };{{\boldsymbol{r}}_{{{\boldsymbol{\varPhi}} }k}},\xi _{{{\boldsymbol{\varPhi}} }k}^{ - 1}{\boldsymbol{I}}} \right)$，通过基本变化可将${b_{{\text{sp}}}}\left( {\boldsymbol{\varPhi}} \right)$表示为

结合${b_{{\text{sp}}}}\left( {\boldsymbol{\varPhi}} \right)$表达式以及因子图，可得${\text{vec}}\left( {\hat {\boldsymbol{\varPhi}} } \right)$的表达式为${\text{vec}}\left( {{\hat {\boldsymbol{\varPhi}} }} \right) = {\left( {{{{\boldsymbol{\varOmega}} }^{\text{T}}}{\xi _{{\boldsymbol{\varSigma}} }}{{\boldsymbol{\varOmega}} } + {\xi _{{{\boldsymbol{\varPhi}} }k}}} \right)^{ - 1}}\Bigr( {{{\boldsymbol{\varOmega}} }^{\text{T}}}{\xi _{{\boldsymbol{\varSigma}} }}{{\boldsymbol{\varOmega}} }{{{\boldsymbol{\varOmega}} }^{\text{ + }}}{\boldsymbol{Y}} + {\xi _{{{\boldsymbol{\varPhi}} }k}}{{\boldsymbol{r}}_{{{\boldsymbol{\varPhi}} }k}} \Bigr)$，$\hat {\boldsymbol{\varPhi}}$可表示为${\text{reshape}}\left( {{\text{vec}}\left( {\hat {\boldsymbol{\varPhi }}} \right),N_{\text{G}}^{\text{2}},{M_{\text{G}}}} \right)$，${\text{reshape}}$函数为MATLAB中的矩阵变形函数。信道的估计值${\boldsymbol{\hat H}}$为${{\boldsymbol{V}}^{\text{T}}}{\boldsymbol{D}}\hat {\boldsymbol{\varPhi}} {\boldsymbol{A}}_{\text{B}}^{\text{H}}$，则在式(16)中阈值函数$g$的表达式为

结合式(16)和式(26)，可将${{\boldsymbol{\hat H}}_k}$的均方误差计算出，结合VAMP算法分析和描述，KR-VAMP算法的仿真流程如下。

步骤1　利用Khatri-Rao积以及克罗内克积的特性将信道${\boldsymbol{G}}$简化，得到简化后的${\boldsymbol{H}}$和${\boldsymbol{Y}}$。

步骤2　设置初始$\mu$, $\xi _{\boldsymbol{\varSigma}} ^{ - 1}$以及迭代总次数$K$，迭代次数初始值设定为$k = 0$。

步骤3　${\tau _k} = \xi _{\boldsymbol{\varSigma}} ^{ - 1} + \dfrac{{{M_{\text{G}}}}}{M}\varepsilon \left( {{\tau _{k - 1}},{\gamma _{k - 1}}} \right)$。

步骤4　${{\boldsymbol{r}}_k} = {{\boldsymbol{H}}^0} + N\left( {0,{\tau _k}{\boldsymbol{I}}} \right)$。

步骤5　${\gamma _k} = {\left( {{\tau _k}} \right)^{ - 1}}$。

步骤6　${{\boldsymbol{\hat H}}_{k + 1}} = g\left( {{{\boldsymbol{r}}_k},{\gamma _k}} \right)$，阈值函数$g$如式(26)所示。

步骤7　$\varepsilon \left( {{\tau _{k + 1}},{\gamma _{k + 1}}} \right) = \dfrac{1}{{{M_{\text{G}}}}}{\rm{E}}\left[ \left\| g\left( {{\boldsymbol{H}}^0} + N\left( 0, {\tau _k}{\boldsymbol{I}} \right), {\gamma _k} \right) - {{\boldsymbol{H}}^0} \right\|^2 \right]$，得到${{\boldsymbol{\hat H}}_{k + 1}}$的均方误差$\varepsilon \left( {\tau _{k + 1}}, {\gamma _{k + 1}} \right)$。

步骤8　将迭代次数$k + 1$，直至其等于$K$，得到${{\boldsymbol{\hat H}}_K}$以及${{\boldsymbol{\hat H}}_k}$的均方误差$\varepsilon \left( {{\tau _k},{\gamma _k}} \right)$，否则循环步骤2—步骤7。

由上述AMP算法的仿真流程和KR-VAMP算法的仿真流程可知，阈值函数$g$的确实对整个算法是至关重要的。

5.   仿真结果与分析

本文主要验证IRS辅助毫米波通信系统的有效性，以及KR-VAMP算法与其他AMP算法在IRS辅助毫米波的系统环境中信道误差数值。结果主要对比了AMP算法、BiG-AMP算法^[5]和KR-VAMP算法的信道误差。仿真数据设置IRS的单元数为16，各IRS单元之间的距离为15 mm。基站处的天线数为2，天线之间的距离设置为15 mm。完备矩阵维度${N_{\text{G}}}$为40，${M_{\text{G}}} = {M^2}$，即设定为4。毫米波信号的功率为50 dBm。

在图4所示的3维坐标中，基站、IRS和用户的坐标分别为$\left( {0,0,0} \right)$, $\left( {5,50,10} \right)$和$\left( {20,y,0} \right)$，假设用户在如图4所示的水平位置移动，其中$y$为用户的移动变量。

图  4  基站、IRS和用户的位置示意图

下载: 全尺寸图片 幻灯片

在上述的前提条件下，当噪声平均功率$\xi _{\boldsymbol{\varSigma}} ^{ - 1}$的值为30 dBm。图5比较了不同迭代初始值$\mu$的条件下，各种算法的信道误差值。由图5可知随着$\mu$值的增加，信道均方误差也逐渐增大，且各算法具有相似的变化趋势。相比于AMP算法，BiG-AMP算法，本文所提KR-VAMP算法最优，在$\mu$值为–30 dB (10^–3)时，MMSE的值可达到–33.8 dB。

图  5  不同迭代初始值条件下的均方误差

下载: 全尺寸图片 幻灯片

当迭代初始值$\mu$的值为10 dBm。图6比较了在不同噪声平均功率条件下，各算法的信道均方误差，由图6可知随着噪声平均功率的增大，信道误差逐渐增大，各算法之间的MMSE差值也增大，在此情况下，KR-VAMP算法性能最佳，当$\xi _{\boldsymbol{\varSigma}} ^{ - 1}$值为–10 dB (10^–1)时，MMSE的值可达到–29 dB。随着$\xi _{\boldsymbol{\varSigma}} ^{ - 1}$增大到20 dB (10²)时，KR-VAMP算法相较于AMP算法和BiG-AMP算法，MMSE差值分别为6.6 dB和3.8 dB。

图  6  不同噪声功率条件下的均方误差

下载: 全尺寸图片 幻灯片

综合上述两个变量，当噪声平均功率$\xi _{\boldsymbol{\varSigma}} ^{ - 1}$的值为30 dBm，迭代初始值$\mu$的值为37 dBm的条件下，图7仿真了MMSE随着迭代次数的变化，由此可见，随着迭代次数的增加，信道误差降低，并且迭代次数增加到10，MMSE也趋于平缓，达到–36 dB。同时也可以观察到，KR-VAMP算法使得此系统的信道均方误差最低，迭代次数较于其他算法也较少。

图  7  算法不同迭代次数下的均方误差

下载: 全尺寸图片 幻灯片

6.   结束语

本文利用IRS辅助毫米波通信，以此降低毫米波在大气中的损耗，提高毫米波nLoS传播路径的性能。在此系统环境中，本文基于VAMP算法提出了KR-VAMP算法，通过Khatri-Rao积和克罗内克儿积将级联信道转换为稀疏信号恢复问题，以及VAMP算法良好的迭代阈值优势，减少了训练迭代次数，提高了系统的性能，降低了整个系统的信道误差。最后仿真验证，与传统的AMP相比，使用KR-VAMP算法可以更好地降低系统的信道误差，对毫米波信道具有一定的性能提升。在本文信道估计过程所使用的VAMP算法还存在一定的复杂度。于是，后续工作将重点围绕降低算法复杂度问题展开研究。