1.   引 言

大规模MIMO和毫米波通信技术作为5G移动通信系统的关键技术，由于其在基站端配置大量天线，导致系统较高的功率消耗；毫米波由于其传输损耗高，并且易受障碍物的阻挡，使得通信的覆盖范围低^[1]。为了降低通信系统的功率消耗，增加通信的覆盖范围，可重构智能表面(Reconfigurable Intelligent Surface, RIS)具有低功耗、低成本且易于部署等优点，可通过调节其反射元件的相位达到对无线环境进行智能调控，从而改变无线通信环境，已成为6G移动通信系统的潜在关键技术之一^[2]。RIS的研究方面主要有：RIS辅助无人机通信^[3]，RIS辅助通信系统能效^[4]，全息多输入多输出表面(Holographic MIMO Surface, HMIMOS)在通信系统中的应用^[5]，RIS辅助物理层安全^[6]，RIS辅助通信感知一体化^[7]等。

基于RIS辅助的毫米波通信系统研究已受到学术界和产业界的广泛关注，其带来通信系统的覆盖率增强，抗衰落和抗干扰等性能提高，然而这些性能提高取决于精确的信道状态信息的获取。但是由于RIS的无源反射特性，难以对信号进行处理。因此，对基于RIS辅助的通信系统的信道估计产生了很大的挑战。

目前，RIS辅助的通信系统的信道估计主要有以下两个方面：一方面为级联信道估计，即联合估计基站(Base Station, BS)到RIS的信道，RIS到用户(User Equipment, UE)之间的信道；另一方面为分离信道估计，即分开估计BS-RIS, RIS-UE信道。在级联信道估计中，文献[8]采用并行因子求解信道，由于该算法未考虑毫米波信道的角域稀疏性，系统的信道估计开销增加。毫米波信道通常被认为是稀疏的，可用压缩感知理论来求解信道估计。文献[9]利用级联信道的稀疏结构将RIS辅助的信道估计问题转化为压缩感知问题，并采用正交匹配追踪算法进行求解。压缩感知类算法是利用信号的角域稀疏性来进行信道角度信息估计，对连续分布的角度进行离散化会产生误差，从而导致信道估计的精度降低。由于RIS在实际的部署中，BS-RIS的相对安装位置往往固定不变，可假设BS-RIS之间的信道是固定不变的。而用户随时可以移动，即 RIS-UE之间的信道是时变的。因此，对BS-RIS和RIS-UE的信道进行分离信道估计更符合实际需求。在分离信道估计方中，文献 [10]采用双时间尺度进行分离信道估计，降低了总体导频开销。为了提高信道估计精度，文献[11]采用基于原子范数最小化(Atomic Norm Minimization, ANM)方法作用于信道角度参数的估计。基于原子范数最小化方法由于对信道角度进行连续的稀疏恢复，因此估计精度有所提升。但该文献考虑的是RIS采用均匀线性阵列(Uniform Linear Array, ULA)结构，不符合实际情况。

为了提高RIS辅助下信道估计的性能，一些研究通过为RIS的少量元素安装射频(Radio Frequency, RF)链来辅助信道估计^[12-14]。文献[12]采用2维多信号分类算法(Two Dimensional MUlti-SIgnal Classification algorithm, 2D-MUSIC)对信道的离开角和到达角度同时进行估计。文献[13]采用求根多重信号分类算法(ROOT MUlti-SIgnal Classification algorithm, ROOT-MUSIC)估计信道的角度参数。以上文献的方法是基于MUSIC算法，其在信噪比较高的情况下估计性能较好。文献[14]采用2维快速傅里叶算法对信道的角度进行估计，并通过信号补零以获得更加精确的信道估计值。采用2维快速傅里叶算法(Two Dimensional Fast Fourier Transform, 2D-FFT)求解信道的角度参数，算法的复杂度较低，但是该算法信道估计的精度与信号补零的个数有关，补零个数的增加虽然可以提高信道估计的精度，但是也会提高算法的复杂度。

为提高RIS辅助的毫米波系统的信道估计的性能，并降低算法的复杂度，考虑给RIS部分无源器件安装RF链，以分离估计BS-RIS, RIS-UE信道。采用基于原子范数最小化的稀疏重构算法，利用信号协方差矩阵的Toeplitz结构特性，有效提高了信道的估计精度。为了降低问题求解的复杂度，采用交替方向乘子法(Alternating Direction Method of Multipliers, ADMM)对信道的参数进行求解，通过联合优化信道参数和迭代步长以获得更加精确的信道估计值，并通过所获得的角度参数和信道矩阵参数获得路径增益，从而最终得到信道估计值。仿真结果表明，本文算法在保持较低计算复杂度的同时，提高了系统信道估计性能。

2.   系统模型

考虑RIS辅助的毫米波上行通信系统的信道估计，系统模型如图1所示。BS和UE分别配备N根和K根天线^[8]，RIS由${N_{{{\rm{RIS}}}}}$个反射元件组成，${N_{{{\rm{RIS}}}}} = {N_1} \times {N_2}$, ${N_1}$, ${N_2}$表示RIS的水平和垂直侧反射元件数目。

图 1  RIS辅助的毫米波上行通信系统模型

下载: 全尺寸图片 幻灯片

由于RIS只反射信号，难以对信号进行处理，因此，考虑在RIS的部分无源器件上安装L型的RF链路。假设UE-BS的直接链路被障碍物所阻挡，UE-BS的通信采用RIS提供的反射路径来辅助通信。系统采用时分双工模式，根据用户发送导频符号进行上行信道估计，下行信道的CSI通过信道的互易性获得。假设RIS辅助的毫米波通信信道的训练阶段由T个连续的时间帧组成，每个时间帧包含Q个时隙^[13]。因此，在第t个时间帧中，导频符号经RIS反射到BS的接收信号${{{\boldsymbol{Y}}}_t}$表示为

其中，${{{\boldsymbol{H}}}_{{\rm{1}}}} \in {{{\boldsymbol{C}}}^{{N_{{{\rm{RIS}}}}} \times {{{K}}}}}$表示UE-RIS之间的传输链路，${{{\boldsymbol{H}}}_2} \in {{{\boldsymbol{C}}}^{N \times {N_{{{\rm{RIS}}}}}}}$表示RIS-BS之间的传输链路，${{\boldsymbol{\varPhi }}} \in {{{\boldsymbol{C}}}^{{N_{{{\rm{RIS}}}}} \times {N_{{{\rm{RIS}}}}}}}$表示RIS的相移矩阵，${{\boldsymbol{\varPhi}}} = {{\rm{diag}}}({a_1}{{{\rm{e}}}^{{{\rm{j}}}{w_1}}}, {a_2}{{{\rm{e}}}^{{{\rm{j}}}{w_2}}}, \cdots, {a_{{N_{{{\rm{RIS}}}}}}}{{{\rm{e}}}^{{{\rm{j}}}{w_{{N_{{{\rm{RIS}}}}}}}}})$，${a_n} \in [0, 1]$表示RIS元件的反射幅度，${w_n} \in \left[ {0, 2\pi } \right)$表示RIS元件的反射相位。为简化研究，本文只考虑信号RIS的相位变化。用户发送的正交导频信息${{{\boldsymbol{X}}}_t} \in {{{\boldsymbol{C}}}^{K \times Q}}$，$Q$为导频符号的长度。${{{\boldsymbol{N}}}_t} \in {{{\boldsymbol{C}}}^{K \times Q}}$表示加性高斯白噪声，${{{\boldsymbol{N}}}_t} \sim {{\rm{CN}}}(0, {{{\boldsymbol{\sigma}}} ^2})$。

3.   信道模型

由于毫米波信道的稀疏性，本文采用S-V信道模型。考虑BS和UE端配备ULA, RIS端配备均匀平面阵列(Uniform Plane Array, UPA)。由于RIS上安装L型RF链，将信道的估计问题分为两部分：第1部分为UE到RIS的信道估计，第2部分为RIS到BS的信道估计。UE-RIS和RIS-BS的信道分别表示为

其中，${L_1}$, ${L_2}$分别表示UE-RIS, RIS-BS的路径个数，${{\boldsymbol{\alpha}} _{{l_1}}}$, ${{\boldsymbol{\beta}} _{{l_2}}}$分别表示UE-RIS的第${l_1}$条路径的复增益，RIS-BS的第${l_2}$条路径的复增益。$({{\boldsymbol{\theta}} _{{l_1}}}, {{\boldsymbol{\phi}} _{{l_1}}})$表示RIS处的第${l_1}$条路径的到达角(Angle Of Arrival, AOA)，$({{\boldsymbol{\varphi}} _{{l_1}}})$表示UE处的第${l_1}$条路径的离开角(Angle Of Departure, AOD)。$({{\boldsymbol{\theta}} _{{l_2}}}, {{\boldsymbol{\phi}} _{{l_2}}})$表示RIS处的第${l_2}$条路径的AOD，$({{\boldsymbol{\varphi}} _{{l_2}}})$表示BS处的第${l_2}$条路径的到达角AOA。其中信道${{{\boldsymbol{H}}}_1}$的天线阵列响应表示为

其中，${{{\boldsymbol{\alpha}}} _1}({{\boldsymbol{\theta}} _{{l_1}}}, {{{\boldsymbol{\phi}}} _{{l_1}}}) = {1 \mathord{\left/{\vphantom {1 {\sqrt {{N_1}} }}} \right.} {\sqrt {{N_1}} }}[1, {{{\rm{e}}}^{{{\rm{j}}}\frac{{2\pi }}{\lambda }d\sin ({{\boldsymbol{\theta}} _{{l_1}}})\sin ({{{\boldsymbol{\phi}}} _{{l_1}}})}}, \cdots , {{{\rm{e}}}^{{{\rm{j}}}\frac{{2\pi }}{\lambda }({N_1} - 1)d\sin ({{\boldsymbol{\theta}} _{{l_1}}})\sin ({{{\boldsymbol{\phi}}} _{{l_1}}})}}]^{{\rm{T}}}$, ${{{\boldsymbol{\alpha}}} _2}({{{\boldsymbol{\phi}}} _{{l_1}}}) = {1 \mathord{\left/{\vphantom {1 {\sqrt {{N_2}} }}} \right.} {\sqrt {{N_2}} }}[1, {{{\rm{e}}}^{{{\rm{j}}}\frac{{2\pi }}{\lambda }d\cos ({{{\boldsymbol{\phi}}}_{{l_1}}})}}, \cdots , {{{\rm{e}}}^{{{\rm{j}}}\frac{{2\pi }}{\lambda }({N_2} - 1)d\cos ({{{\boldsymbol{\phi}}} _{{l_1}}})}}]^{{\rm{T}}}$。$\lambda$表示毫米波波长，$d$表示天线阵列的间距。${{{\boldsymbol{H}}}_2}$中的 ${{{\boldsymbol{b}}}_{{\rm{BS}}}}({{{\boldsymbol{\varphi }}}_{{l_2}}})$，${{\boldsymbol{b}}}_{{\rm{RIS}}}^{{\rm{H}}} ({{{\boldsymbol{\boldsymbol{\theta}} }}_{{l_2}}}, {{{\boldsymbol{\phi}}} _{{l_2}}})$和${{{\boldsymbol{H}}}_1}$中的${{\boldsymbol{\alpha }}}_{{\rm{UE}}}^{{\rm{H}}}({{{\boldsymbol{\varphi }}} _{{l_1}}})$，${{{\boldsymbol{\alpha }}}_{{\rm{RIS}}}}({{\boldsymbol{\theta}} _{{l_1}}}, {{{\boldsymbol{\phi}}} _{{l_1}}})$表示形式相同，只需将${{{\boldsymbol{\varphi }}} _{{l_2}}}$用${{{\boldsymbol{\varphi }}} _{{l_1}}}$代替，$({{\boldsymbol{\theta}} _{{l_2}}}, {{{\boldsymbol{\phi}}} _{{l_2}}})$ 用$({{\boldsymbol{\theta}} _{{l_1}}}, {{{\boldsymbol{\phi}}} _{{l_1}}})$代替即可。

4.   信道估计

本文采用分离信道估计方法，根据式(2)可知，只需要估计信道的角度参数以及路径增益，即可恢复信道矩阵。由于安装RF链，RIS有两种工作模式，反射模式和接收模式。在估计UE-RIS信道时，RIS工作在接收模式，可利用RF链对RIS接收的信号进行处理，从而在RIS端进行估计该信道；在估计RIS-BS的信道时，RIS工作在反射模式下，BS接收RIS反射的信号，从而在BS端对该信道进行估计^[14]。

4.1   ${{{\boldsymbol{H}}}_1}$信道矩阵参数估计

首先对UE-RIS的信道${{{\boldsymbol{H}}}_1}$中的角度参数以及路径增益进行估计，即角度参数需要估计UE的AOD和RIS的AOA。由于RIS采用UPA结构，估计AOA需要估计RIS的方位角和仰角两种角度信息。传统的基于原子范数最小化的方法不能直接对3维角度进行估计。因此，利用RIS上安装的L型RF链的水平方向阵列估计RIS的方位角和垂直方向阵列估计RIS的仰角，在各个方向上只需要估计两种角度信息。因此，本文的信道矩阵参数估计方法步骤如下：

首先证明RIS上RF链接收信号的协方差矩阵满足双重Toeplitz结构，采用基于解耦原子范数最小化算法进行求解；其次，为了降低该 半正定规划(Semi-Definite Programming, SDP)问题求解的复杂度，采用改进的快速ADMM算法进行求解，该改进算法采用动量梯度下降法进行信道参数变量更新，并通过联合优化迭代步长和信道矩阵参数，提高信道估计的精度，通过该改进算法即得到UE处的AOD和RIS的方位角的2维角度信息和信道矩阵参数；最后，通过所获得的角度参数和信道矩阵参数估计路径增益。根据以上步骤，可获得UE处的AOD和RIS处的AOA，即可恢复信道矩阵。${{{\boldsymbol{H}}}_1}$信道矩阵参数的估计过程如图2所示。

图 2  ${{{\boldsymbol{H}}}_1}$信道矩阵参数的估计过程

下载: 全尺寸图片 幻灯片

4.1.1   ${{{\boldsymbol{H}}}_1}$信道矩阵角度参数估计

RIS端部分无源器件上配置“L”型RF链路，该RF链在水平和竖直方向上的导向矢量分别表示为

因此，RIS上射频链在水平方向接收的信号${{\boldsymbol{y}}}_t^x$表示为

其中，${{\boldsymbol{H}}}_1^x$为UE到RIS水平RF链的信道矩阵，${{\boldsymbol{H}}}_1^x = \sqrt {{{ K{N_{{\rm{RIS}}}}}/ {{L_1}}}} \sum\nolimits_{{l_1} = 1}^{{L_1}} {{a_{{l_1}}} {{\boldsymbol{\alpha }}}_{{\rm{RIS}}}^x( {{\boldsymbol{\theta}} _{{l_1}}} ){{\boldsymbol{\alpha }}}_{{\rm{UE}}}^{{\rm{H}}} ( {{{\boldsymbol{\varphi}}} _{{l_1}}} ) = } {{\boldsymbol{A}}}_{{\rm{RIS}}}^x{{\rm{diag}}}( {{\boldsymbol{A}}} ) {{\boldsymbol{A}}}_{{\rm{UE}}}^{{\rm{H}}}$ , ${{\boldsymbol{A}}} = \sqrt {{{K{N_{{\rm{RIS}}}}} \mathord{\left/{\vphantom {{K{N_{{\rm{RIS}}}}} {{L_1}}}} \right.} {{L_1}}}} [{a_1}, {a_2}, \cdots , {a_{{L_1}}}]$, ${{\boldsymbol{A}}}_{{\rm{RIS}}}^x = [{{\boldsymbol{\alpha }}}_{{\rm{RIS}}}^x({{\boldsymbol{\theta}} _1}), {{\boldsymbol{\alpha }}}_{{\rm{RIS}}}^x({{\boldsymbol{\theta}} _2}), \cdots , {{\boldsymbol{\alpha }}}_{{\rm{RIS}}}^x({{\boldsymbol{\theta}} _{{L_1}}})]$, ${{\boldsymbol{A}}}_{{\rm{UE}}} = [{{\boldsymbol{\alpha }}}_{{\rm{UE}}}^{}({{{\boldsymbol{\varphi}}} _1}), {{\boldsymbol{\alpha }}}_{{\rm{UE}}}^{}({{{\boldsymbol{\varphi}}} _2}), \cdots , {{\boldsymbol{\alpha }}}_{{\rm{UE}}}({{{\boldsymbol{\varphi}}} _{{L_1}}})]$，${{\boldsymbol{N}}}_t^x$表示高斯白噪声。对式(5)进行向量化得

其中，$\otimes$, $( \circ )$, $(*)$, ${{\rm{vec}}}(·)$分别表示矩阵的Kronecher积、Khatri-Rao积、共轭、向量化。${{{\boldsymbol{I}}}_{{{{nx}}}}}$表示单位矩阵，采用T个时间帧进行训练时，式(6)可以写为

其中，${{{\boldsymbol{Y}}}^{{\boldsymbol{x}}}} \in {\mathbb{C}^{TQ{n_x} \times 1}}$, ${{\boldsymbol{W}}} \in {\mathbb{C}^{TQ{n_x} \times K{n_x}}}$。RIS水平方向上接收信号的协方差矩阵为

其中，${{\boldsymbol{T}}} \in {\mathbb{C}^{K{n_x} \times K{n_x}}}$为半正定具有双重Toeplitz结构的矩阵。由于${{\boldsymbol{T}}}$表达式中含有UE-RIS水平链路信道矩阵的离开角和到达角参数。因此通过求解${{\boldsymbol{T}}}$，可得到信道矩阵的角度参数。由于${{\boldsymbol{T}}}$满足双重Toeplitz结构，因此，将求解UE-RIS水平链路的信道估计问题转化为基于解耦原子范数最小化的SDP问题

其中，${{\boldsymbol{T}}}({{{\boldsymbol{u}}}_x}) \in {\mathbb{C}^{{n_x} \times {n_x}}}$, ${{\boldsymbol{T}}}({{{\boldsymbol{u}}}_y}) \in {\mathbb{C}^{K \times K}}$分别代表以${{{\boldsymbol{u}}}_x}$, ${{{\boldsymbol{u}}}_y}$作为第1行得到的Toeplitz矩阵，这两个矩阵为半正定的Hermitian-Toeplitz矩阵，${\lambda _1}$为正则化参数，用于平衡信号稀疏性和信道估计的准确性^[15]，$\left\| {{{{\boldsymbol{Y}}}^x} - {{\boldsymbol{WH}}}_1^x} \right\|_{{\rm{F}}}^2$ 为噪声干扰项，${\Vert \cdot\Vert }_{{{\rm{F}}}}$表示矩阵的Frobenius范数。通过式(9)，可以将角度信息${{\boldsymbol{\varphi}}}$和${\boldsymbol{\theta}}$进行解耦。根据原子范数最小化算法要求，其任意两角度之间应满足以下的角度分辨率，即${\varDelta _{\min , {\boldsymbol{\theta}} }} = \mathop {\min }\limits_{{l_1} \ne {l_2}} \left| {\sin ({{\boldsymbol{\theta}} _{{l_1}}}) \;-\; \sin ({{\boldsymbol{\theta}} _{{l_2}}})} \right| \;\ge\; {4 \mathord{\left/{\vphantom {4 {{n_x}}}} \right.} {{n_x}}}$, ${\varDelta _{\min , {{\boldsymbol{\varphi}}} }} \;=\; \mathop {\min }\limits_{{l_1} \ne {l_2}} \left| {\sin ({{{\boldsymbol{\varphi}}} _{{l_1}}}) - \sin ({{{\boldsymbol{\varphi}}} _{{l_2}}})} \right| \ge {4 \mathord{\left/{\vphantom {4 K}} \right.} K}$。

式(9)可通过 MATLAB中凸优化工具箱(ConVeX, CVX)进行求解，由于通过CVX求解的算法复杂度为求解的半正定矩阵维度的3.5次方^[16]，因此，本文采用改进ADMM算法以降低算法的复杂度。ADMM算法的核心在于：通过对问题进行分布式优化，即对一个多变量优化问题进行解耦，通过对每个单变量进行迭代求解来简化问题本身的计算量^[17]。因此，通过引入一个新的变量${{\boldsymbol{U}}}$，将原始的优化变量的约束去掉，即将式(9)改写为

式(10)对应的增广拉格朗日函数为

其中，${{\boldsymbol{F}}}$为式(10)的目标函数。$\left\langle {{{\boldsymbol{a}}}, {{\boldsymbol{b}}}} \right\rangle$ 表示矩阵${{\boldsymbol{a}}}$与${{\boldsymbol{b}}}$的内积，即$\left\langle {{{\boldsymbol{a}}}, {{\boldsymbol{b}}}} \right\rangle = {{\rm{tr}}}({{{\boldsymbol{b}}}^{{\rm{H}}}}{{\boldsymbol{a}}})$，${{\boldsymbol{\varLambda }}}$为拉格朗日乘子，$p$为惩罚项，$p = 1$。通过式(11)，将式(10)转化为无约束优化问题(先忽略${{\boldsymbol{U}}} \succcurlyeq 0$约束，在接下来求解${{\boldsymbol{U}}}$时考虑这一约束)。因此，根据ADMM算法，优化变量的更新步骤为

由于在式(11)中$({{\boldsymbol{H}}}_1^x, {{{\boldsymbol{u}}}_x}, {{{\boldsymbol{u}}}_y})$的优化变量分别为${{\boldsymbol{L}}}$的凸函数。因此，可对这些优化变量进行求导，定义${{{\boldsymbol{U}}}^t} = \left[ \begin{array}{l}\;{{\boldsymbol{U}}}_0^t\;\;\;\;{{\boldsymbol{U}}}_1^t \\{\left( {{{\boldsymbol{U}}}_1^t} \right)^{{\rm{H}}}}{{\boldsymbol{U}}}_2^t \\ \end{array} \right]$, ${{{\boldsymbol{\varLambda }}}^t} = \left[ \begin{array}{l}\;{{\boldsymbol{\varLambda }}}_0^t\;\;\;\;{{\boldsymbol{\varLambda }}}_1^t \\{\left( {{{\boldsymbol{\varLambda }}}_1^t} \right)^{{\rm{H}}}}{{\boldsymbol{\varLambda }}}_2^t \\ \end{array} \right]$，接下来推导式(12a)的具体表达式，即式(11)对${{\boldsymbol{H}}}_1^x$求导，得

若令式(13)求导结果为0，则可直接计算得到${{\boldsymbol{H}}}_1^x$的更新表达式，但是会涉及矩阵的求逆运算，增加运算复杂度。因此，为了避免矩阵的求逆运算，本文采用动量梯度下降法来更新信道参数。由于动量梯度下降法的迭代步长影响信道估计参数的精确度，当选取的迭代步长较小时，导致迭代次数较高，算法复杂度升高；当迭代步长较大时，又会降低信道估计的精确性。本文采用对迭代步长和信道参数进行联合优化，迭代步长的选择采用Armijo准则，令步长因子${\eta _{{{t}}}} = {\beta ^{{{{\rm{m}}}_{{{t}}}}}}$，其中${{{\rm{m}}}_{{{t}}}}$的选取满足

其中，$\beta \in (0, 1)$, $\xi \in (0, 0.5)$，${{{\boldsymbol{d}}}_t}$为式(11)函数的下降方向。采用动量梯度下降法对信道矩阵参数进行更新，其更新公式为

其中，$\delta$为动量项，且初始化为0，$\gamma$为关于动量项的超参数，$\gamma = 0.9$。因此，本文联合优化迭代步长和信道参数的计算步骤如算法1所示。

算法1　联合优化迭代步长和信道参数算法
输入：初始化${{\boldsymbol{H}}}_1^x$, $\beta = 0.5$, $\xi = 0.2$，容许误差$\varepsilon$, $t = 1$；
(1)计算${ {\boldsymbol{g} } }_t^{} = {{\text{∇}} _{ { {({ {\boldsymbol{H} } }_{ {\rm{1} } }^x{ {\rm{)} } } }^t} } }{ {\boldsymbol{L} } }$，如果$\left\| {{{{\boldsymbol{g}}}_t}} \right\| \le \varepsilon$，输出${({{\boldsymbol{H}}}_1^x)^t}$作为最优解；
(2)令${{{\boldsymbol{d}}}_t}{{\boldsymbol{ = }}} - {{{\boldsymbol{g}}}_t}$；
(3)根据式(14)确定迭代步长${\beta ^{{m_t}}}$；
(4)计算式(15)，$t = t + 1$，转步骤(1)。
输出：${({{\boldsymbol{H}}}_1^x)^t}$, ${\beta ^{{m_t}}}$。

下载: 导出CSV 

| 显示表格

根据算法1获得${{\boldsymbol{H}}}_1^x$, ${\beta ^{{m_t}}}$的表达式。接下来对$({{{\boldsymbol{u}}}_x}, {{{\boldsymbol{u}}}_y})$进行更新，令式(11)对$({{{\boldsymbol{u}}}_x}, {{{\boldsymbol{u}}}_y})$的导数为0

其中，${{{\boldsymbol{T}}}^*}$是${{\boldsymbol{T}}}$的伴随矩阵，${{{\boldsymbol{e}}}_1} = {\left[ {1, 0, \cdots , 0} \right]^{{\rm{T}}}}$。接下来计算${{\boldsymbol{U}}}$的更新表达式，由于式(11)中只有最后一项与${{\boldsymbol{U}}}$有关，并考虑${{\boldsymbol{U}}} \succcurlyeq 0$约束，即寻找使式(11)最后一项最小化的半正定矩阵，令${{\boldsymbol{\varXi }}} = \left[ \begin{array}{l}{{\boldsymbol{T}}}({{{\boldsymbol{u}}}_y})\;\;\;{({{\boldsymbol{H}}}_1^x)^{{\rm{H}}}} \\\;\;{{\boldsymbol{H}}}_1^x\;\;\;\;\;{{\boldsymbol{T}}}({{{\boldsymbol{u}}}_x}) \\ \end{array} \right] - {p^{ - 1}}{{\boldsymbol{\varLambda }}}$，对${{\boldsymbol{\varXi }}}$进行特征值分解，即${{\boldsymbol{\varXi }}} = {{\boldsymbol{E\varSigma }}}{{{\boldsymbol{E}}}^{{\rm{H}}}}$，因此，${{\boldsymbol{U}}}$的更新表达式为

其中，${{{\boldsymbol{\varSigma }}}_ + }$表示将${{\boldsymbol{\varSigma }}}$中所有负的特征值置为0后，得到的新的特征值矩阵。因此，根据式(15)—式(17)迭代更新$({{\boldsymbol{H}}}_1^x, {{{\boldsymbol{u}}}_x}, {{{\boldsymbol{u}}}_y}, {{\boldsymbol{U}}})$的值。

与采用CVX求解器求解的算法复杂度相对比，本节采用改进的ADMM算法，本算法利用动量梯度下降法对信道参数进行更新，避免矩阵的求逆过程，并通过对迭代步长和信道矩阵参数的联合优化，提高信道估计的精度。因此，本文所提算法在保持信道估计的精确性的同时，算法复杂度降低，运行速度更快。

根据估计得到的${{\boldsymbol{\tilde T}}}({{{\boldsymbol{u}}}_x})$, ${{\boldsymbol{\tilde T}}}({{{\boldsymbol{u}}}_y})$，对${{\boldsymbol{\tilde T}}}({{{\boldsymbol{u}}}_x})$, ${{\boldsymbol{\tilde T}}}({{{\boldsymbol{u}}}_y})$进行范德蒙德分解，可求解得到角度信息的表达式

由于式(18)中范德蒙德分解的形式与MUSIC算法中接收信号协方差矩阵的形式一样。因此，采用Root-MUSIC算法，转化为多项式的求根问题，即可得到$\tilde {{\boldsymbol{\varphi}}}$和${\boldsymbol{\tilde \theta}}$的角度信息，进而根据matrix pencil算法进行角度的自动配对^[18]。根据上述的基于解耦原子范数最小化算法和改进的ADMM算法可以得到$\tilde {{\boldsymbol{\varphi}}}$和${\boldsymbol{\tilde \theta}}$的角度信息，同理可以用相同的方法，得到RIS上的仰角${{\boldsymbol{\tilde \phi}}}$的估计值。

4.1.2   ${{{\boldsymbol{H}}}_1}$路径增益估计

根据4.1.1节改进的ADMM算法可求得信道的估计值${{\boldsymbol{\tilde H}}}_1^x$和角度的估计值$({{\boldsymbol{\tilde \varphi}}} , {\boldsymbol{\tilde \theta}} )$，因此，根据式(5)${{\boldsymbol{H}}}_1^x$的定义可知，信道的路径增益${{{\tilde {\boldsymbol{A}}}}}$可通过信道角度参数和信道矩阵参数估计，即${{{\tilde {\boldsymbol{A}}}}}$表达式为

其中，${(\cdot)}^{†}$表示矩阵的Moore-Penrose伪逆，其中${({{\boldsymbol{A}}}_{{\rm{RIS}}}^x)^\dagger } \;=\;{({({{\boldsymbol{A}}}_{{\rm{RIS}}}^x)^{{\rm{H}}}}\;({{\boldsymbol{A}}}_{{\rm{RIS}}}^x))^{ - 1}}\;{({{\boldsymbol{A}}}_{{\rm{RIS}}}^x)^{{\rm{H}}}}$, ${({{\boldsymbol{A}}}_{{\rm{UE}}}^{{\rm{H}}})^\dagger } = {({{\boldsymbol{A}}}_{{\rm{UE}}}^{}{{\boldsymbol{A}}}_{{\rm{UE}}}^{{\rm{H}}})^{ - 1}}{{\boldsymbol{A}}}_{{\rm{UE}}}^{}$。

4.1.3   信道矩阵重建

根据4.1.1节的解耦原子范数最小化算法和改进的ADMM算法可获得UE处的AOD和RIS处的AOA；根据4.1.2节方法可获得${{{\boldsymbol{H}}}_1}$的路径增益。因此，根据式(2)，可获得重建${{{\boldsymbol{\tilde H}}}_1}$矩阵为

4.2   ${{{\boldsymbol{H}}}_2}$信道估计

根据上文中的算法可以得到UE-RIS信道矩阵的估计值${{{\boldsymbol{\tilde H}}}_1}$，设此时的RIS工作在反射模式下，BS接收到RIS的水平和垂直RF链反射的信号，从而在BS端对该信道进行估计。BS端接收的信号表达式为

其中，${{{\boldsymbol{\tilde Y}}}_t} \in {\mathbb{C}^{N \times Q}}$表示BS端接收的信号，${{\boldsymbol{Y}}}_t^{{\rm{RIS}}}$表示RIS上配备的RF链的接收信号。要估计RIS-BS的信道矩阵${{{\boldsymbol{H}}}_2}$，可根据4.1节所提出的D-ANM+改进的ADMM算法进行求解，即得到${{{\boldsymbol{H}}}_2}$的估计值

因此，本文所提基于RIS辅助的毫米波信道估计的D-ANM+改进ADMM算法流程图如算法2所示。

算法2 D-ANM+改进ADMM信道估计算法
输入：BS, RIS, UE处的天线数目、导频信号、RIS射频链接收 　的信号、噪声的功率，最大迭代次数${s_{\max }}$；
计算UE-RIS的信道矩阵${{\boldsymbol{H}}}_1^{}$：
(1)根据式(11)求解得到RIS水平射频链路接收信号的协方差矩阵；
(2)根据解耦原子范数最小化方法，并采用改进的ADMM算法求 　解式(17)，其信道参数的迭代过程如下：
While $s \lt {s_{\max }}$
(a)根据算法1更新${({{\boldsymbol{H}}}_1^x)^t}$, ${\beta ^{{m_t}}}$；
(d)根据式(16)更新$({{{\boldsymbol{u}}}_x}, {{{\boldsymbol{u}}}_y})$；
(e)根据式(17)更新${{\boldsymbol{U}}}$；
(f)根据式(12c)更新${{\boldsymbol{\varLambda }}}$；
(g)$s = s + 1$
End while
(3)根据步骤(2)求解得到信道的估计值${{\boldsymbol{\tilde H}}}_1^x$, ${{\boldsymbol{\tilde T}}}({{{\boldsymbol{u}}}_x})$, ${{\boldsymbol{\tilde T}}}({{{\boldsymbol{u}}}_y})$；
(4)采用Root-MUSIC算法求解角度和matrix pencil进行角度匹配 　求解信道矩阵${{\boldsymbol{H}}}_1^{}$的UE端的离开角估计值$\tilde {{\boldsymbol{\varphi}}}$和RIS端的水平射频 　链的AOA的方位角估计值${\boldsymbol{\tilde \theta}}$；并以相同的思路估计得到RIS端垂 　直射频链的AOA的仰角估计值${{\boldsymbol{\tilde \phi}}}$；
(5)根据式(19)，求解信道矩阵${{\boldsymbol{H}}}_1^{}$的路径增益${{\boldsymbol{\tilde \alpha}}}$；
(6)根据以上得到的信道参数值$({\boldsymbol{\tilde \theta}} , {{\boldsymbol{\tilde \phi}}} , {{\boldsymbol{\tilde \varphi}}} , {{\boldsymbol{\tilde \alpha}}} )$，得到估计的信道矩 　阵${{\boldsymbol{\tilde H}}}_1^{}$；
RIS-BS的信道矩阵${{{\boldsymbol{H}}}_2}$的估计参照步骤(1)～(6)的思路。
输出：UE-RIS估计的信道矩阵${{\boldsymbol{\tilde H}}}_1^{}$，RIS-BS估计的信道矩阵 　${{{\boldsymbol{\tilde H}}}_2}$。

下载: 导出CSV 

| 显示表格

5.   算法复杂度分析

本文所提D-ANM+改进ADMM信道估计算法的复杂度主要由信道矩阵${{{\boldsymbol{H}}}_1}$和${{{\boldsymbol{H}}}_2}$中的角度参数估计和路径增益估计两部分组成。由于求解信道矩阵${{{\boldsymbol{H}}}_2}$的算法与信道矩阵${{{\boldsymbol{H}}}_1}$的算法相同，因此只分析求解${{{\boldsymbol{H}}}_1}$的算法复杂度。由于采用动量梯度下降法对信道矩阵参数进行更新，在更新时，需要求解${{\boldsymbol{H}}}_1^x$的梯度，该过程仅涉及矩阵的乘法运算，其复杂度为$O({n_x}K{L_1} + {n_x}K)$，对$({{{\boldsymbol{u}}}_x}, {{{\boldsymbol{u}}}_y}, {{\boldsymbol{\varLambda }}})$更新时，只涉及每个矩阵的加减运算，计算复杂度可忽略，对${{\boldsymbol{U}}}$进行更新时，由于需要进行特征值分解，计算复杂度为$O({n_x}{({n_x} + K)^2})$。因此所提出的算法的复杂度为$O({n_x}K{L_1} + {n_x}K + {n_x}{({n_x} + K)^2})$；其次，采用ROOT-MUSIC算法和matrix pencil算法对角度进行求解和配对，其求解的计算复杂度表示为$O((n_x^2{K^2}){L_1})$；最后求解信道的路径增益，其计算的复杂度为$O({({L_1}{n_x})^3})$。由于${L_1} \ll \left\{ {N, K, {N_{{\rm{RIS}}}}} \right\}$, ${n_x} \ll {N_{{\rm{RIS}}}}$。因此，求解${{{\boldsymbol{H}}}_1}$的算法复杂度为$O({n_x}{({n_x} + K)^2})$。本文算法和文献[8, 12-14]的算法复杂度对比分析如表1所示。

表 1  本文算法与相关文献算法的复杂度对比

算法	信道估计类别	算法复杂度
D-ANM+改进ADMM	分离信道估计	$O({n_x}{({n_x} + K)^2})$
ALS[8]	级联信道估计	$O(2{N^3} + 4{N^2}P(M + K) - NP(M + K))$
2D-MUSIC[12]	分离信道估计	$O(n_x^2(T + {D_1}{D_2}) + {N^3})$
ROOT-MUSIC[13]	分离信道估计	$O(n_x^2T + {N^3})$
2D-FFT[14]	分离信道估计	$O(n_x^2(M + K))$

下载: 导出CSV 

| 显示表格

6.   仿真分析

本节中通过仿真实验来评估本文所提的D-ANM+改进ADMM信道估计算法与现有文献中的ALS算法^[8]、2D-MUSIC算法^[12]、ROOT-MUSIC算法^[13]、2D-FFT算法^[14]系统性能的对比。在仿真中，设置BS端天线的数目$N = 32$，其中RIS的水平/垂直射频链的数目${n_x} = {n_y}$，导频符号数$Q = 20$, $T = 50$。UE-RIS和RIS-BS信道的路径${L_1} = {L_2} = 4$。根据原子范数算法的角度分辨率要求，UE-RIS和RIS-BS信道的不同路径之间的正弦差$\sin ({{\boldsymbol{\theta}} _1}) - \sin ({{\boldsymbol{\theta}} _2}) \ge {4 \mathord{\left/{\vphantom {4 N}} \right.} N}, {4 \mathord{\left/{\vphantom {4 K}} \right.} K}, {4 \mathord{\left/{\vphantom {4 {{N_{{\rm{RIS}}}}}}} \right.} {{N_{{\rm{RIS}}}}}}$。

采用均方根误差(Root Mean Squared Error, RMSE)评价信道矩阵的角度参数的性能指标，表达式为

其中，$R$表示蒙特卡罗运行的次数，$R = 500$，$({\boldsymbol{{\theta}}} , {\boldsymbol{\varphi}} )$表示RIS的方位角和仰角真实的角度值，$({\tilde{\boldsymbol{\theta}}} , {{\tilde {\boldsymbol{\varphi}}}} )$表示RIS的方位角和仰角估计的角度值。

采用两种系统性能评估指标对本文算法与相关文献算法的信道估计性能进行对比，第1种性能指标为：归一化均方误差(Normalized Mean Square Errror, NMSE)，其表达式为

其中，${{\boldsymbol{H}}}$表示真实的信道值，${{\boldsymbol{\tilde H}}}$表示估计的信道值。

第2种性能指标为频谱效率(Spectral Efficiency, SE)，由于本文不涉及发送端的基带预编码器和接收端的合并器设计，因此，对信道矩阵${{\boldsymbol{H}}}$进行奇异值分解，假设预编码器${{\boldsymbol{W}}}$和合并器${{\boldsymbol{F}}}$的取值分别为所估计的信道矩阵${{\boldsymbol{H}}}$的${L_1}$个左右奇异向量^[11]。考虑RIS的离散相移矩阵${{\boldsymbol{\varPhi }}}$，从$\left\{ - {1 \mathord{\left/{\vphantom {1 {\sqrt {{N_{{\rm{RIS}}}}} }}} \right.} {\sqrt {{N_{{\rm{RIS}}}}} }}, {1 \mathord{\left/{\vphantom {1 {\sqrt {{N_{{\rm{RIS}}}}} }}} \right.} {\sqrt {{N_{{\rm{RIS}}}}} }} \right\}$中选择RIS的相移矩阵${{\boldsymbol{\varPhi }}}$的每个元素。因此，根据文献[11]，频谱效率的表达式为

其中，信道矩阵${{\boldsymbol{H}}} = {{{\boldsymbol{\tilde H}}}_2}{{\boldsymbol{\varPhi }}}{{{\boldsymbol{\tilde H}}}_1}$, ${{\boldsymbol{\hat H}}} = {{\boldsymbol{H}}} - {{\boldsymbol{\tilde H}}}$表示信道估计的误差，${\rm{var}}(\cdot)$表示求解矩阵的协方差。

图3给出了克拉美罗下界(Cramer-Rao Lower Bound, CRLB)^[19]、本文所提D-ANM+改进ADMM算法、2D-MUSIC算法和2D-FFT算法估计得到的RIS的角度的均方根误差性能对比。由图可知，本文所提算法的角度估计的均方根误差和CRLB相接近，这是由于本文采用解耦原子范数最小化算法，该算法是利用信号的协方差矩阵的Toeplitz特性，对连续的角度进行稀疏化，避免了网格失配问题，因此提高了角度估计的精度。随着信噪比的增加，本文算法角度估计的均方根误差和CRLB逐渐接近。而2D-MUSIC算法是将连续的角度离散化，角度的估计精度与离散化的角度之间的间隔有关，本文采用${0.1^ \circ }$的角度间隔来搜索。2D-FFT算法根据信号空间谱进行谱峰搜索以获取信道的角度估计值，由于采用信号补零方法增加信号的频谱点，因此，信道估计的性能和信号补零的个数有关，本文假设对信号补零后的维度为256。因此，本文算法角度估计值达到了良好的性能。

图 3  不同算法求解RIS的角度均方根误差性能对比

下载: 全尺寸图片 幻灯片

图4给出了UE-RIS的信道${{{\boldsymbol{H}}}_1}$, RIS-BS的信道${{{\boldsymbol{H}}}_2}$估计的归一化均方误差系统的性能对比。由图4可知，随着RIS上配置的射频链个数的增加，信道估计的性能也逐渐增加。这是由于射频链数目的增加，信道估计的准确性增加。随着信噪比的增加，信道${{{\boldsymbol{H}}}_1}$和信道${{{\boldsymbol{H}}}_2}$的估计性能逐渐增加，这是由于信噪比的增加，信道环境变好，因此信道估计性能增加。该结果验证了本文所提的信道估计算法在UE-RIS和RIS-BS分离信道估计时均能保持良好的系统性能，由此验证了本文算法的可行性。

图 4  信道${{{\boldsymbol{H}}}_1}$与信道${{{\boldsymbol{H}}}_2}$估计性能对比

下载: 全尺寸图片 幻灯片

图5当信噪比为20 dB时，本文提出的改进ADMM算法与文献[11]采用CVX求解器算法的运行时间对比。由图可知，随着用户天线数的增加，改进的ADMM算法的运行时间曲线的斜率远小于采用CVX求解的斜率，即改进ADMM算法的运行速度比CVX求解器快得多。采用CVX求解的SDP问题的算法是采用内点算法求解的，其算法复杂度为$O((n_x^{3.5}{K^{3.5}}) {\log _2}({1 \mathord{\left/{\vphantom {1 \varepsilon }} \right.} \varepsilon }))$^[16]，$\varepsilon$为目标的精度。采用ADMM求解SDP问题的算法复杂度为$O({n_x}{({n_x} + K)^2})$，由于ADMM在每次迭代中的参数更新都是以封闭形式进行的，所以运行速度较快。因此，本文所提出的算法运行时间较快，并且比较适合大型阵列。

图 5  本文算法与CVX算法的运行时间对比

下载: 全尺寸图片 幻灯片

图6给出了当RIS上配置的射频链数目分别为4, 10, 14时，系统的频谱效率性能的变化情况。由图可知，随着射频链数目的增加，频谱效率也逐渐增加。当射频链数目增加时，能接收到更多的信道信息，因此系统的信道估计性能提高，频谱效率增加。

图 6  RIS射频链数目对频谱效率性能的影响

下载: 全尺寸图片 幻灯片

图7给出了当RIS上元件数目为16, 32, 64时，对信道估计的归一化均方误差性能的影响。由图7可知，当RIS上元件数增加时，信道估计的性能在逐渐下降。这是由于当RIS上元件数增加时，其信道估计所需的参数也在增加，从而导致系统性能下降。而且当信噪比逐渐增加时，信道估计的性能也逐渐上升。这是由于当信噪比增加时，通信信道的状态良好，信道的估计性能提高。

图 7  RIS元件数对信道${{\boldsymbol{H}}}$归一化均方误差性能的影响

下载: 全尺寸图片 幻灯片

图8和图9分别给出了当RIS上无源元件配置的射频链数目为8时，本文算法与现有ALS算法、ROOT-MUSIC算法、2D-FFT算法、2D-MUSIC算法的归一化均方误差和频谱效率系统性能的对比。由图可知，2D-MUSIC算法与本文算法的系统性能相近，而2D-MUSIC算法采用的是信号子空间与噪声子空间正交的性质得到信号的谱峰，需要更多的观测数据，更多的观测数据会导致导频开销增加。本文采用文献[12]中的2D-MUSIC算法仅使用了2D-MUSIC算法估计信道的角度参数，未采用其使用深度学习方法估计信道增益的方案。因为深度学习方法会带来更大的计算开销。2D-FFT算法由于信号的空间谱存在栅栏效应，采用信号补零方式以提高信道估计的精度，当信号补零后的维度增加时，信道估计的精确性进一步提高，但信号的维度增加，会增加信号的处理时间。为避免MUSIC算法的2维谱峰搜索，采用ROOT-MUSIC算法以多项式求根方式求解信道的角度参数，其角度估计精度较低。因此，本文所提算法可以提高信道估计的精确性，并且达到了系统性能和复杂度的平衡。

图 8  本文算法与其他文献算法的归一化均方误差系统性能对比

下载: 全尺寸图片 幻灯片

图 9  本文算法与其他算法频谱效率性能对比

下载: 全尺寸图片 幻灯片

7.   结论

本文针对RIS辅助的毫米波通信系统提出一种低复杂度的信道估计算法。通过在RIS部分的无源器件上配备射频链，以分离估计BS-RIS, RIS-UE的信道。所提算法首先通过解耦原子范数最小化将信道的离开角和到达角2维角度估计问题解耦为两个1维的角度估计问题，其次采用动量梯度下降法避免ADMM算法中信道矩阵参数更新时的矩阵求逆运算，并通过联合优化迭代步长和信道参数，提高信道参数的估计精度，最后根据角度参数和信道矩阵参数获得信道路径增益。通过仿真实验验证本文所提算法在保持较低复杂度的同时提高了系统的信道估计性能。