Two-Stage Channel Estimation for Semi-Passive RIS-Assisted Millimeter Wave Systems

Abstract

In a reconfigurable intelligent surface (RIS) assisted millimeter Wave (mmWave) communication system, the channel coefficient increases exponentially with the number of RIS elements which results in expensive pilot overhead. Most previous works have proposed some channel estimation algorithms for the estimation accuracy of cascaded channels, which have improved the estimation accuracy, but the pilot overhead is discouraging in the estimation process. To improve the channel estimation accuracy with reduced pilot overhead, we propose a two-stage channel estimation protocol by exploiting semi-passive elements and the coherent time difference of the channel, where the quasi-static channel between the base stations (BS) and RIS is estimated at the RIS, and the user (UE)-RIS time-varying channel is estimated at the BS. In the first stage, we formulate the BS-RIS channel estimation as a mathematical optimization problem by an iterative weighting method and then propose a gradient descent (GD)-based algorithm to solve it. In the second stage, we first transform the received the UE-RIS signal model into an equivalent parallel factor (PARAFAC) tensor model and estimate the UE-RIS channel by the least-squares (LS) algorithm. The simulation results show that the proposed method has better estimation accuracy than the LS, compression sensing (CS) and minimum mean square error (MMSE) methods with less pilot overhead, and the spectral efficiency is improved by at least 10.5% compared to the other three methods.

1. Introduction

The emerging reconfigurable intelligent surface (RIS) technology has been identified as a promising key technology for future 6G communications. RIS consists of a large number of low-cost passive reflective elements in a planar array, which enhance communication performance by subtly changing the direction of signal transmission [1,2,3,4]. Specifically, by digital reconfigurable/programmable meta-surface technology, each passive element appropriately adjusts the amplitude and phase of the incident signal to flexibly configure the wireless channel between the transmitter and receiver, and reduce fading impairment and interference problems in the wireless channel, thus enhancing communication performance. Existing research indicates that RIS can significantly improve communication rates [5,6], cancel interference, and expand network coverage, etc. However, the performance gain achieved by RIS relies on accurate channel state information (CSI). Consequently, accurate CSI estimation is critical for RIS-assisted communication systems.

However, the following two problems make acquiring CSI extremely difficult. Firstly, due to a large number of reflective elements of the RIS [7], the channel coefficients that need to be estimated subsequently increase, which leads to an expensive pilot overhead. Secondly, conventional RIS lacks active transmitters and receivers with signal processing capability, and it is difficult to estimate the exact the base station (BS)-RIS and RIS-user equipment (UE) channels separately, further complicating channel estimation. As a result, channel estimation is challenging for RIS-aided wireless communications.

There has been limited research on channel estimation for RIS-assisted communication systems to obtain accurate CSI. Most works are devoted to estimating cascaded channels [8,9,10,11,12,13,14,15]. The works [8,9] estimated cascaded channels by using the conventional least-squares (LS) algorithm, which was computationally simple, but it was sensitive to noise and does not measure with high accuracy in the case of few pilots. In [10], the least squares algorithm was improved by the proposed minimum mean square error (MMSE) channel estimator, which reduced the effect of noise on the channel estimation and improved the estimation accuracy. The work [11] described the channel estimation problem as an estimation error minimization problem and proposed a Lagrange multiplier and a pairwise ascending one to solve this problem, which improved the accuracy to some extent but increased the computational complexity. The work [12] proposed a three-stage channel estimation framework to improve the estimation accuracy by exploiting the channel consistency among UEs, but the pilot overhead still grows proportionally to the RIS elements. The work [13,14] proposed a sparse matrix decomposition and complementary channel estimation method based on sparse matrix decomposition to reduce the pilot overhead by exploiting channel sparsity. In addition, the work [15] proposed compressed sensing(CS) techniques to estimate channels for high-frequency millimeter wave(mmWave), exploiting the sparsity of mmWave to reduce the pilot overhead, but the computational complexity grows exponentially with the number of RIS elements.We summarize previous work in

Table 1. Summary of previous work.

RIS Configuration	Previous Work
Passive RIS	Channel estimation based on conventional least-squares (LS) algorithm for the RIS-assisted MIMO system [8,9].
	The MMSE algorithm is used to estimate the cascade channel for the RIS-assisted MIMO system [10].
	Channel estimation based on the methods of Lagrange multipliers and a dual ascent-based algorithm for the RIS-assisted MIMO system [11].
	Reference user-based channel estimation by using the common BS- RIS channel of the RIS-assisted MISO system [12].
	Channel estimation based on the methods of a sparse matrix decomposition and complementary channel estimation method for the RIS-assisted MIMO system [13,14].
	Compression-sensing-based channel estimation based on sparse re-presentation of cascaded channel [15].
Semi-passive RIS	Channel estimation based on compressed sensing for the RIS-assist-ed SISO system [16,17,18].

.

Despite the high channel estimation performance achieved in the above works, it is difficult to estimate the BS-RIS channel and the UE-RIS channel separately due to the passive nature of the RIS. In other words, if each RIS element can sense and reflect the signal, CSI can be easily estimated separately at the RIS. However, the method requires a large number of semi-passive elements with radio frequency (RF) chains, which results in higher hardware cost and power consumption. To tackle this problem, a small number of semi-passive components were proposed in [16,17,18] for receiving and processing the pilot signals in the RIS, and the BS-RIS channel and the UE-RIS channel were estimated separately based on the received pilot signals by the CS algorithm. However, the CS algorithm usually involves complex non-convex optimization problems, and it is difficult to obtain optimal solutions.

Based on the above discussion, we can observe that the problem of obtaining accurate CSI with low pilot overhead requires more research. The paper investigates the channel estimation problem for a semi-passive RIS-assisted mmWave multiple-input multiple-output (MIMO) system. We propose a two-stage channel estimation scheme to achieve high accuracy channel estimation with low pilot overhead. The main contributions of this paper are as follows:

• For the problem of expensive pilot overhead due to many channel coefficients, we introduce semi-passive elements to assist channel estimation and then propose a channel estimation protocol based on the coherence time difference between BS-RIS quasi-static channels and time-varying UE-RIS channels, in which the BS-RIS channels are estimated at long time scales using semi-passive elements and the BS estimates the UE-RIS channels at short time scales.The proposed channel estimation protocol is able to reduce the average pilot overhead at long time scales.
• To estimate BS-RIS channel, we propose an iterative re-weighting-based super-resolution algorithm to estimate the BS-RIS channel. We transform the BS-RIS channel estimation problem to the optimization problem of a new objective function by an iterative weighting method, which is the weighted summation of the sparsity and the data fitting error, Then we propose a gradient descent method to solve the objective function problem and update the weight parameters at each iteration to balance the sparsity with the data fitting error. During the iterative process, the estimated parameters move gradually to the neighborhood of the true value.Compared to traditional algorithms, the proposed algorithm is able to converge the estimates to near the true value and achieve accurate estimates.
• To estimate the time-varying channel of the UE-RIS, we propose a LS algorithm based on parallel factor(PARAFAC) decomposition to estimate the time-varying channel of the UE-RIS. We transform the received signal model into an equivalent PARAFAC tensor model, then obtain the the UE-RIS channel by LS algorithm.The proposed algorithm has higher robustness compared to the traditional algorithm by using PARAFAC decomposition, which avoids the problem of non-existence of matrix inverse.

The rest of this paper is presented as follows: Section 2 introduces the MIMO system model with the assistance of semi-passive RIS and the channel transmission protocol. In Section 3, we propose a channel estimation algorithm for the BS-RIS channel. In Section 4, we propose a channel estimation algorithm for the UE-RIS channel. The simulation results and the conclusion remarks are provided in Section 5 and  Section 6, respectively.

2. Channel Model and Channel Estimation Protocol

In this section, we present the channel model of the RIS-assisted mmWave communication system and the proposed channel estimation protocol.

2.1. Channel Model

As shown in Figure 1, we consider a semi-passive RIS-assisted mmWave MIMO communication system, which consists of an M antenna base station, a semi-passive RIS, and K single-antenna users. The RIS consists of $N_0$ semi-passive elements and $N_1$ passive elements, while the semi-passive elements have an RF chain capable of receiving and processing signals when it is turned on, among $N_0+N_1=N$, $N_0<<N$. Due to the presence of semi-passive elements, the RIS has two modes: receive mode and reflect mode. In the receive mode, the semi-passive elements are used for channel estimation. In the reflection mode, the semi-passive elements and passive elements reflect the signal to the receiver. In this paper, we assume that the channel between BS and UE is blocked. Then there are two channels in this path, which are the BS-RIS channel $G\in {\mathbb{C}}^{N\times M}$ and the UE-RIS channel $h_{r,k}\in {\mathbb{C}}^{N\times K}$. In addition, since the adjacent semi-passive RIS elements have stronger channel correlation with each other, we assume that the semi-passive elements are continuously distributed over the RIS to better estimate the CSI.

We consider that both the BS and the RIS are uniform planar arrays (UPA) [19]; then the channel between the BS-RIS can be expressed as

$$G=\sum _{l_1=1}^L{a}_{l_1}{a}_{RIS}\left({\phi }_{R,l_1},{\theta }_{R,l_1}\right)a_{BS}^H\left({\phi }_{T,l_1},{\theta }_{T,l_1}\right),$$

(1)

where L is is the number of paths between the BS-RIS, $a_{l_1}$ is the complex gain of the $l_1$th path, ${\phi }_{l_1}$ and ${\phi }_{l_2}$ is the azimuth angle, ${\theta }_{l_1}$ and ${\theta }_{l_2}$ is the elevation angle, $a_{RIS}\left({\phi }_{R,l_1},{\theta }_{R,l_1}\right)$, and $a_{BS}^H\left({\phi }_{T,l_1},{\theta }_{T,l_1}\right)$ is the antenna array response of the base station and RIS. The specific expressions can be expressed as

$$\begin{array}{c}\hfill a_{RIS}\left({\phi }_{R,l_1},{\theta }_{R,l_1}\right)=\left[1,e^{j\frac{2\pi }{\lambda }d sin\left({\phi }_{R,l_1}\right)cos\left({\theta }_{R,l_1}\right)},\dots ,e^{j\frac{2\pi }{\lambda }d{N}_{1}sin\left({\phi }_{R,l_1}\right)cos\left({\theta }_{R,l_1}\right)}\right]\\ \hfill \otimes \left[1,e^{j\frac{2\pi }{\lambda }d sin\left({\phi }_{\left.R,l_1\right)}\right)},\dots ,e^{j\frac{2\pi }{\lambda }d{N}_2 sin\left({\phi }_{\left.R,l_1\right)}\right)}\right],\end{array}$$

(2)

$$\begin{array}{c}\hfill a_{BS}\left({\phi }_{T,l_1},{\theta }_{T,l_1}\right)=\left[1,e^{j\frac{2\pi }{\lambda }d sin\left({\phi }_{T,l_1}\right)cos\left({\theta }_{T,l_1}\right)},\dots ,e^{j\frac{2\pi }{\lambda }d{M}_{1}sin\left({\phi }_{T,l_1}\right)cos\left({\theta }_{T,l_1}\right)}\right]\\ \hfill \otimes \left[1,e^{j\frac{2\pi }{\lambda }d sin\left({\phi }_{T,l_1}\right)},\dots ,e^{j\frac{2\pi }{\lambda }d{M}_2 sin\left({\phi }_{T,l_1}\right)}\right],\end{array}$$

(3)

where $\lambda$ is the wavelength, $d=\frac{\lambda }{2}$ is the spacing between the antennas, $M=M_1\times M_2$, $N=N_1\times N_2$, and ⊗ is the Kronecker product. Similarly, we consider a uniform linear array (ULA) at the UE; the channel between the UE-RIS can be expressed as

$$h_{r,k}=\sum _{l_{2,k}=1}^{L_{r,k}}{\alpha }_{l_{2,k}}{a}_{RIS}\left({\phi }_{R,l_{2,k}},{\theta }_{R,l_{2,k}}\right){\alpha }_t^H\left({\varphi }_{T,l_{2,k}}\right),$$

(4)

where $L_{r,k}$ is is the number of paths between the UE-RIS, ${\alpha }_{l_2}$ is the complex gain of the $l_2$th path between the UE-RIS, ${\phi }_{l_1}$ and ${\phi }_{l_2}$ is the azimuth angle, ${\theta }_{l_1}$ and ${\theta }_{l_2}$ is the elevation angle, $a_{BS}\left({\phi }_{R,l_1},{\theta }_{R,l_1}\right)$, and $a_{RIS}^H\left({\phi }_{T,l_1},{\theta }_{T,l_1}\right)$ is the antenna array response of the base station and RIS. Similarily, the specific expressions can be expressed as

$$a_{UE}\left({\varphi }_{T,l_{2,k}}\right)=\left[1,e^{j\frac{2\pi }{\lambda }d sin\left({\varphi }_{T,L_{2,k}}\right)}\dots e^{j\frac{2\pi }{\lambda }d(K-1) sin\left({\varphi }_{T,l_{2,k}}\right)}\right],$$

(5)

$$\begin{array}{c}\hfill a_{RIS}({\phi }_{R,l_{2,k}},{\theta }_{R,l_{2,k}})=[1,e^{j\frac{2\pi }{\lambda }d{N}_1 sin\left({\phi }_{R,l_{2,k}}\right)cos\left({\theta }_{R,l_{2,k}}\right)},\dots ,e^{j\frac{2\pi }{\lambda }d{N}_1 sin\left({\phi }_{R,l_{2,k}}\right)cos\left({\theta }_{R,l_{2,k}}\right)}]\\ \hfill \otimes [1,e^{j\frac{2\pi }{\lambda }d sin\left({\phi }_{R,l_{2,k}}\right)},\dots ,e^{j\frac{2\pi }{\lambda }d{N}_2 sin\left({\phi }_{R,l_{2,k}}\right)}],\end{array}$$

(6)

where K is the number of users.

2.2. Channel Estimation Protocol

By exploiting the channel coherence time differences and semi-passive elements, we propose a two-stage channel estimation protocol in a semi-passive RIS architecture. Firstly, because the RIS and BS locations are relatively fixed and the channel change rate is relatively slow, the BS-RIS channel is quasi-static.Then the estimated BS-RIS channel can be estimated once in a long time. In addition, due to the mobility of the UE, the UE-RIS channel is time-varying. Then the UE-RIS needs to be estimated frequently in a short time. Secondly, because of the semi-passive elements, both the BS-RIS channel and the UE-RIS channel can be estimated by RIS, but the UE-RIS has a faster channel change speed compared to the BS-RIS. If the UE-RIS time-varying channel is estimated on RIS, then the semi-passive elements need to be turned on frequently to achieve the optimal channel estimation performance, but this will increase the burden of RIS. Therefore, we only estimate the quasi-static channel of the BS-RIS on RIS in the first stage, and the time-varying channel of the UE-RIS is estimated by the BS in the second stage. In general, estimating the frequency of the BS-RIS is much lower than estimating the frequency of the UE-RIS, which reduces the average pilot overhead from a long time scale. Compared to the cascaded channel estimation method that estimates $MNK$ coefficients, the proposed channel estimation protocol can recover the BS-RIS and the UE-RIS channels by estimating only $N_{0}M+NK$ coefficients. Therefore, the required pilot overhead is further reduced.

As shown in Figure 2. There are T sub-frames during the signal transmission, where the first $T_1\ll T$ sub-frames estimates the BS-RIS channel and the second $T_2$ sub-frames are used to estimate the time-varying channel of the UE-RIS, and the channel is kept constant between each small time block. Specifically, in the $T_1$ sub-frames of the first stage, we consider the downlink communication, the BS transmits the pilot signal to the semi-passive RIS, the passive element turns off, the semi-passive element turns on and receives the pilot signal, and feeds the BS-RIS channel estimation result to the BS. In the $T_2$ sub-frames of the second stage, the BS estimates the UE-RIS channel by the uplink pilot signal sent by the UE. After the channel estimation is completed, the semi-passive element turns off the sensing mode and is used to reflect the signal together with the passive element during the data transmission. Based on this transmission protocol, we will discuss how to effectively estimate the BS-RIS channel versus the UE-RIS channel in the next two sections.

3. Stage 1: Estimation of the BS-RIS Channel

When the RIS is equipped with a small number of semi-passive elements, the CS algorithm can recover the complete channel of the BS-RIS using channel sparsity. However, this method usually requires solving a complex non-convex problem, it is difficult to obtain an optimal solution. In light of this, we propose an iterative re-weight-based channel estimation algorithm, which can estimate the BS-RIS channel with high accuracy by converting the non-convex problem into a new objective function by the iterative re-weight-based method and optimally solving it using the gradient descent method.

3.1. Downlink Pilot Transmission and Optimization Formulation

3.1.1. Downlink Pilot Transmission

In the first stage, all the $N_0$ semi-passive elements operate in the sensing mode, while other passive components are off and do not participate in signal reflection. The BS sends pilot signal $x\in {\mathbb{C}}^{M\times 1}$ to the RIS, then the signal is received at the semi-passive through a random combining matrix W. Thus, the received signal at the semi-passive after Q time slots is summarized as

$$Y_{RIS}=W^H\overline{G}X+N,$$

(7)

$$\overline{G}=\mathsf{\Gamma }G,$$

(8)

where $Y_{RIS}=\left[y_1,y_2,\dots y_Q\right]\in {\mathbb{C}}^{N_0\times Q}$, $W\in {\mathbb{C}}^{N_0\times N_0}$ is the random composite matrix, $X=\left[x_1,x_2,\dots x_Q\right]\in {\mathbb{C}}^{M\times Q}$, $\mathsf{\Gamma }\in {\mathbb{C}}^{N_0\times N}$ the selection matrix for selecting the semi-passive elements on the RIS, $\overline{G}\in {\mathbb{C}}^{N_0\times M}$ is the channel between the $\mathrm{BS}$ and the RIS semi-passive elements, and $N\in {\mathbb{C}}^{N_0\times Q}$ is the additive white Gaussian noise (AWGN) with $\mathcal{C}\mathcal{N}\left(0,{\sigma }_n^2\right)$.

As shown in Equation (1), to accurately estimate the BS-RIS channel, we need to estimate the complex gain of all paths ${\alpha }_{l_1}$, the departure angle ${\alpha }_t({\phi }_{T,l_1}$, ${\theta }_{T,l_1})$, and the arrival angle $a_r({\phi }_{\mathrm{R},l_1}$, ${\theta }_{R,l_1})$. In the framework of mmWave, the channel G of the BS-RIS can be transformed into

$$G=A_R\left(\mathsf{\Phi }\right)diag\left(a\right)A_T^H\left(\mathsf{\Psi }\right),$$

(9)

where $\mathsf{\Phi }={\left[{\phi }_{R,1},{\theta }_{R,1},{\phi }_{R,2},{\theta }_{R,2},\dots ,{\phi }_{R,L},{\theta }_{R,L}\right]}^T$, $\mathsf{\Psi }={\left[{\phi }_{T,1},{\theta }_{T,1},{\phi }_{T,2},{\theta }_{T,2},\dots ,{\phi }_{T,L},{\theta }_{T,L}\right]}^T$, $a={[a_1,a_2,\dots ,a_L]}^T$, $A_R\left(\mathsf{\Phi }\right)=\left[a_{RIS}\left({\phi }_{R,1},{\theta }_{R,1}\right),a_{RIS}\left({\phi }_{R,2},{\theta }_{R,2}\right),\dots ,a_{RIS}\left({\phi }_{R,L},{\theta }_{R,L}\right)\right]$, $A_T\left(\mathsf{\Psi }\right)=\left[a_{BS}\left({\phi }_{T,1},{\theta }_{T,1}\right),a_{BS}\left({\phi }_{T,2},{\theta }_{T,2}\right),\dots ,a_{BS}\left({\phi }_{T,L},{\theta }_{T,L}\right)\right]$.

3.1.2. Optimization Formulation

The above channel $\overline{G}$ estimation problem can be expressed as

$$\underset{a.\mathsf{\Phi }.\mathsf{\Psi }}{min}{\parallel a\parallel }_0\mathrm{s}.\mathrm{t}{∥Y_{RIS}-W_t^H\mathsf{\Gamma }{A}_{R}a{A}_T^{H}X∥}_F\le \epsilon ,$$

(10)

where $\epsilon$ is a custom error threshold. Obviously, this is a sparse signal recovery problem for which the CS algorithm has good results, but its estimation accuracy is greatly affected by off-grid effects and the base mismatch problem, as well as the need to solve non-convex optimization problems with high computational complexity. To over-fit the data during the channel estimation process and cause inaccurate estimation, we add a regularization parameter and reformulate the above-obtained signal as

$$\begin{array}{cc}\hfill \underset{a.\mathsf{\Phi }.\mathsf{\Psi }}{min}\mathcal{L}\left(a,\mathsf{\Phi },\mathsf{\Psi }\right)& ={\parallel a\parallel }_0\hfill \\ \hfill & +\xi {∥{\widehat{Y}}_{RIS}-W_t^H\mathsf{\Gamma }{A}_{R}a{A}_T^{H}X∥}_F^2,\hfill \end{array}$$

(11)

where $\xi >0$ is used to weigh the sparse constraint term against the error constraint term. Since the above problem is a non-convex optimization problem, it is difficult to find an optimal solution. Based on ref. [20], we use the logarithm and function instead of $L_0$ particularization; then we can obtain

$$\begin{array}{cc}\hfill \underset{a.\mathsf{\Phi }.\mathsf{\Psi }}{min}\mathcal{L}\left(a,\mathsf{\Phi },\mathsf{\Psi }\right)& =\sum _{l=1}^{L}log\left({\left|a_l\right|}^2+\delta \right)\hfill \\ \hfill & +\xi {∥{\widehat{Y}}_{RIS}-W_t^H\mathsf{\Gamma }{A}_{R}a{A}_T^{H}X∥}_F^2,\hfill \end{array}$$

(12)

where $\delta >0$ is to ensure that the weight is not 0. Furthermore, solving the above problem remains a non-convex problem, and we use an iterative surrogate function instead of a logarithmic function. Then the minimization problem of $\mathcal{L}\left(z,{\theta }_R,{\theta }_T\right)$ is equivalent to the iterative optimization problem for the surrogate function.

$$\begin{array}{cc}\hfill \underset{a.\mathsf{\Phi }.\mathsf{\Psi }}{min}\mathcal{F}\left(a,\mathsf{\Phi },\mathsf{\Psi }\right)& ={\xi }^{-1}{a}^H{\mathsf{\Lambda }}^{\left(i\right)}a\hfill \\ \hfill & +{∥Y_{RIS}-W_t^H\mathsf{\Gamma }{A}_{R}a{A}_T^{H}X∥}_F^2,\hfill \end{array}$$

(13)

where

$${\mathsf{\Lambda }}^{\left(i\right)}=diag\left(\frac{1}{\left|{\widehat{a}}_1\left(i\right)\right|+\delta }\frac{1}{\left|{\widehat{a}}_2\left(i\right)\right|+\delta }\cdots \frac{1}{\left|{\widehat{a}}_L\left(i\right)\right|+\delta }\right),$$

(14)

where $\widehat{a}\left(i\right)$ is the value of the value of the ith iteration. To obtain the optimal $\mathcal{F}\left(a,{\theta }_R,{\theta }_T\right)$, Equation (16) can be rewritten as

$$\begin{array}{cc}\hfill \mathcal{F}\left(a,\mathsf{\Phi },\mathsf{\Psi }\right)& ={\xi }^{-1}{a}^H{\mathsf{\Lambda }}^{\left(i\right)}a+\sum _{p=1}^Q{∥y_p-W_t^H\mathsf{\Gamma }{A}_{R}a{A}_T^{H}X∥}_2^2\hfill \\ \hfill & ={\xi }^{-1}{z}^H{\mathsf{\Lambda }}^{\left(i\right)}a+\sum _{p=1}^Q{\left(y_p-{\Xi }_{p}a\right)}^H\left(y_p-{\Xi }_{p}a\right)\hfill \\ \hfill & =a^H\left({\xi }^{-1}{\mathsf{\Lambda }}^{\left(i\right)}+\sum _{p=1}^Q{\Xi }_p^H{\Xi }_p\right)a\hfill \\ \hfill & -a^H\left(\sum _{p=1}^Q{\Xi }_p^H{y}_p\right)-\left(\sum _{p=1}^Q{y}_p^H{\Xi }_p\right)a\hfill \\ \hfill & +\sum _{p=1}^Q{y}_p^H{y}_p,\hfill \end{array}$$

(15)

where ${\Xi }_p=W_t^H\mathsf{\Gamma }{A}_{R}diag\left(A_T^H{x}_p\right)$. Then taking the partial derivative of this gives

$$\begin{array}{cc}\hfill \frac{\partial \mathcal{F}\left(a,\mathsf{\Phi },\mathsf{\Psi }\right)}{\partial a}=& a^H\left({\xi }^{-1}{\mathsf{\Lambda }}^{\left(i\right)}+\sum _{p=1}^Q{\Xi }_p^H{\Xi }_p\right)\hfill \\ \hfill & -\left(\sum _{p=1}^Q{y}_p^H{\Xi }_p\right),\hfill \end{array}$$

(16)

By setting the derivative to zero, the minimum point a and the corresponding minimum value of $a_{\mathrm{opt}}\left(\mathsf{\Phi },\mathsf{\Psi }\right)$ can be obtained as

$$a_{opt}^{\left(i\right)}\left(\mathsf{\Phi },\mathsf{\Psi }\right)=\left({\xi }^{-1}{\mathsf{\Lambda }}^{\left(i\right)}+\sum _{p=1}^Q{\Xi }_p^H{\Xi }_p\right)\left(\sum _{p=1}^Q{\Xi }_p^H{y}_p\right),$$

(17)

Substituting this into the original Equation (15) yields ${\mathcal{F}}_{opt}\left({\theta }_R,{\theta }_T\right)$ the optimal value of

$$\begin{array}{cc}\hfill {\mathcal{F}}_{\mathrm{opt}}^{\left(i\right)}\left(\mathsf{\Phi },\mathsf{\Psi }\right)=& -{\left(\sum _{p=1}^Q{\Xi }_p^H{y}_p\right)}^H{\left({\xi }^{-1}{\mathsf{\Lambda }}^{\left(i\right)}+\sum _{p=1}^Q{\Xi }_p^H{\Xi }_p\right)}^{-1}\hfill \\ \hfill & \left(\sum _{p=1}^Q{\Xi }_p^H{y}_p\right)+\sum _{p=1}^Q{y}_p^H{y}_p,\hfill \end{array}$$

(18)

In this way, we obtain an estimation of the channel gain a. Then the channel estimation problem is further simplified to the estimation of parameters $\mathsf{\Phi }$ and $\mathsf{\Psi }$. The specific solving process will be described in the next section.

3.2. Propose Super-Resolution Channel Estimation Scheme

In the previous subsection, we transformed the BS-RIS channel estimation problem into a new optimization problem. Based on the optimization problem derived above, we propose an iterative re-weight-based channel estimation algorithm, as shown in steps 7 to 16 of Algorithm 1. In Equation (12), the former term constrains the sparsity of the estimation results, and the latter constrains the fitting error to avoid over-fitting. Specifically, a larger $\xi$ can lead to more accurate solutions, while a smaller $\xi$ will lead to over-sparse solutions and poor estimation results. Therefore, to have a good estimate for the iteration closer to the true value and obtain a smaller error, we update the regularization parameter in the iteration $\xi$. The updated equation is shown below

$$\xi =min\left(\frac{R}{r_i^2},{\xi }_{max}\right),$$

(19)

where R is a custom scale factor, which is determined by factors such as signal-to-noise ratio and frequency interval. Here, ${\xi }_{max}$ is used to limit the range of $\xi$ to ensure that the algorithm is feasible, and $r_i$ is the residual from the previous step, i.e.,

$$r_i={∥Y-W^H\mathsf{\Gamma }{A}_R\left({\widehat{\mathsf{\Phi }}}^i\right)diag\left({\widehat{a}}_i\right)A_R\left({\mathsf{\Psi }}^i\right)X∥}_F,$$

(20)

Moreover, to reduce the impact of quantization error on the channel estimation, we use the gradient descent method in machine learning to traverse the imaginary angle domain to obtain ${\mathsf{\Phi }}_R^{(i+1)}$ that ${\mathsf{\Psi }}_R^{(i+1)}$ which achieves super-resolution [21], with the following updated equation

$${\mathsf{\Phi }}^{(i+1)}={\mathsf{\Phi }}^{\left(i\right)}-\eta {\nabla }_{{\mathsf{\Phi }}_R}{\mathcal{F}}_{opt}^{\left(i\right)}\left({\mathsf{\Phi }}^{\left(i\right)},{\mathsf{\Psi }}^{\left(i\right)}\right),$$

(21)

$${\mathsf{\Psi }}^{(i+1)}={\mathsf{\Psi }}^{\left(i\right)}-\eta {\nabla }_{{\mathsf{\Psi }}_R}{\mathcal{F}}_{opt}^{\left(i\right)}\left({\mathsf{\Phi }}^{\left(i\right)},{\mathsf{\Psi }}^{\left(i\right)}\right),$$

(22)

where $\eta$ is the step size to make sure ${\mathcal{F}}_{\mathrm{opt}}^{\left(i\right)}\left({\mathsf{\Phi }}^{(i+1)},{\mathsf{\Psi }}^{(i+1)}\right)\le {\mathcal{F}}_{\mathrm{opt}}^{\left(i\right)}\left({\mathsf{\Phi }}^{\left(i\right)},{\mathsf{\Psi }}^{\left(i\right)}\right)$, and ∇ is gradient calculation operation. In addition, the number of reflected paths of mmWave are generally unknown in practical applications. During the iterative process of the proposed scheme, paths with too little gain will be considered as noise and pruned out, which makes the result more sparse, and the true path number is finally obtained by iterative pruning.

In the iterative update process, if randomly selected initial values cause significant computational complexity, then choosing the initial values for the iterations can reduce the number of iterations and thus achieve a reduction in complexity, as discussed in the next subsection.

3.3. SVD Algorithm of Preconditioning

To make the calculation easier, we adopt the singular value decomposition (SVD) algorithm to pre-process the received signal [21], as shown in steps 1 to 6 in Algorithm 1, to find the angular domain grid of arrival/departure(AoA/AoD) that is closest to the true value. Compared with setting the initial value arbitrarily, the pre-process can significantly reduce the computational complexity. The SVD decomposition of the received signal yields

$$Y_{RIS}=U\mathsf{\Sigma }{V}^H,$$

(23)

where $Y_{RIS}$ is the signal received by the semi-passive element, $\mathsf{\Sigma }=diag\left({\xi }_1,{\xi }_2,\dots ,{\xi }_{min(N_0,Q)}\right)$$\in {\mathbb{R}}^{N_0\times Q},\mathrm{U}$ is the left singular matrix, $\mathrm{V}$ is the right singular matrix, and $U^{H}U=I_{N_0\times N_0}$, and $V^{H}V=I_{Q\times Q}$. According to Equation (10), it is obtained that

$$Y_{RIS}=\left(W_t^H\mathsf{\Gamma }{A}_R\right)a{\left(X{A}_T\right)}^H+N,$$

(24)

Algorithm 1 Two-stage channel estimation algorithm

Stage 1: Estimation of the BS-RIS Channel.

Input: Receive signal $Y_{RIS}=[y_1,y_2,\dots ,y_Q]$, combination matrix W, pilot signal $X=\left[X_1,X_2,\dots X_Q\right]$, BS-RIS channel G, selection matrix $\mathsf{\Gamma }$, trimming threshold ${\Sigma }_t$, termination thresholds $\epsilon$ and the number of paths $L_{init}$ to detect. Output: Estimated $A\left(\mathsf{\Phi }\right)$ with $A\left(\mathsf{\Psi }\right)$ and the path gain for each path.

1. $Y_{RIS}=U\mathsf{\Sigma }{V}^H$

2. Take the first $L_{init}$ columns of U, V and $L_{init}$ largest singular values.

3. for $i=1,2,3,\dots ,N_{it}$ do

4. Calculated from Equations (26) and  (27).

5. end for

6. Output ${\widehat{\phi }}_R^{\left(0\right)}$, ${\widehat{\theta }}_R^{\left(0\right)}$ and ${\widehat{\phi }}_T^{\left(0\right)}$, ${\widehat{\theta }}_T^{\left(0\right)}$.

7. Initialize $a_{\mathrm{opt}}$ according to Equation (17).

8. Repeat:

9. Update according to Equation (19) $\xi$.

10. Calculated from Equation (18) ${\mathcal{F}}_{\mathrm{opt}}\left(\mathsf{\Phi },\mathsf{\Psi }\right)$.

11. Iterate according to Equations (21) and  (22) to find the optimal ${\widehat{\phi }}_R^{\left(i\right)}$, ${\widehat{\theta }}_R^{\left(i\right)}$ and ${\widehat{\phi }}_T^{\left(i\right)}$, ${\widehat{\theta }}_T^{\left(i\right)}$.

12. Calculate the path gain according to Equation (17).

13. Trim path number l if $a_l^{(i+1)}<{\mathsf{\Sigma }}_t$.

14. until $L^i=L^{(i+1)}$ and ${∥a^{(i+1)}-a^i∥}_2<\epsilon$

15. $\mathsf{\Phi }={\widehat{\mathsf{\Phi }}}^{\left(\mathrm{last}\right)}$, $\mathsf{\Psi }={\widehat{\mathsf{\Psi }}}^{\left(\mathrm{last}\right)}$, $a=a^{\mathrm{last}}$.

16. Reconstructed from Equation (9) ${\mathrm{G}}_{es}$.

Stage 2: Estimation of the UE-RIS Channel.

Input: Receive signal $Y_k=\left[Y_1,Y_2,\dots Y_T\right]$, combination matrix W, pilot signal $X=\left[X_1,X_2,\dots X_T\right]$, the BS-RIS estimated channel $G_{es}$, and the UE-RIS channel $h_{r,k}$.

Output: $h_{r_{e}s}$.

17. The PARAFAC decomposition channel problem is obtained according to Equation (34).

18.Estimate $h_{r,es}$ according to Equation (35).

In SVD decomposition, the size of the singular value is proportional to the amount of information it bears, i.e., the larger the singular value, the closer the association with the true information. Furthermore, when the noise is small, the maximum L singular values are approximated by the number of paths, such that when i = 1, 2, $\cdots L$, we have

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{u}_i\approx W^H\mathsf{\Gamma }{a}_R\left({\phi }_{R,l_i},{\theta }_{R,l_i}\right)/{∥W^H\mathsf{\Gamma }{a}_R\left({\phi }_{R,l_i},{\theta }_{R,l_i}\right)∥}_2\hfill \\ {\xi }_i\approx \mid a_{l_i}\mid {∥W^H\mathsf{\Gamma }{a}_R\left({\phi }_{R,l_i},{\theta }_{R,l_i}\right)∥}_2{∥X^H{a}_T\left({\phi }_{R,l_i},{\theta }_{R,l_i}\right)∥}_2\hfill \\ v_i\approx X^H{a}_R\left({\phi }_{R,l_i},{\theta }_{R,l_i}\right)/{∥X^H\mathsf{\Gamma }{a}_R\left({\phi }_{R,l_i},{\theta }_{R,l_i}\right)∥}_2,\hfill \end{array}\right.\end{array}$$

(25)

where $u_i$ and $v_i$ are the i-th column of the left singular matrix U and the right singular matrix V, respectively. At this point, a crude estimate of the AOD and AOA can be obtained for

$$A_R({\widehat{\phi }}_{R,l_i}^{\left(0\right)},{\widehat{\theta }}_{R,l_i}^{\left(0\right)})=\underset{({\phi }_R,{\theta }_R)\in {\mathsf{\Omega }}_1}{\arg \max }{u}_i{W}^H\mathsf{\Gamma }{A}_R({\phi }_{R,l_i},{\theta }_{R,l_i}),$$

(26)

$$A_T({\widehat{\phi }}_{T,l_i}^{\left(0\right)},{\widehat{\theta }}_{T,l_i}^{\left(0\right)})=\underset{({\theta }_T,{\theta }_T)\in {\mathsf{\Omega }}_2}{\arg \max }{v}_i{X}^H{A}_T({\phi }_{T,l_i},{\theta }_{T,l_i}),$$

(27)

where ${\mathsf{\Omega }}_1=\left\{(\frac{i}{N_3},\frac{j}{N_4})\mid i=1,\dots ,N_3-1;j=1,\dots ,N_4-1\right\}$,$N_3\times N_4=N_0$. Similarly, ${\mathsf{\Omega }}_2=\left\{(\frac{i}{M_1},\frac{j}{M_2})\mid i=1,\dots ,M_1-1;j=1,\dots ,M_2-1\right\}$, $M_1\times M_2=M$.

By the SVD method, a coarse estimate around the exact value is obtained, and using it as the initial value of the iteration can move the estimate to the vicinity of the true value faster, thus reducing the computational complexity.

3.4. Pilot Overhead and Computational Complexity

3.4.1. Pilot Overhead

In the first stage of estimating the BS-RIS channel, we use $N_0$ semi-passive elements to recover the $MN$ coefficients of the BS-RIS, then RIS obtains at $N_{0}Q$ measurements for each sub-frame. We need at least $\left[\frac{MN}{N_{0}Q}\right]$ sub-frames. The minimum pilot overhead for the first stage is ${\tau }_1=T_1=Q\left[\frac{MN}{N_{0}Q}\right]$.

3.4.2. Computational Complexity

In the first stage of the BS-RIS channel estimation, the algorithmic complexity of the super-resolution algorithm for solving the optimization problem is $\mathcal{O}\left(MT\left(M+N_0\right)N_0^2{M}^2\right)$, and the complexity is further reduced by the SVD algorithm to $\mathcal{O}\left(MT\left(M+N_0\right)N_0^2{L}^2\right)$.

4. Stage 2: Estimation of the UE-RIS Channel

In this section, based on the BS-RIS channels estimated in the previous stage $G_{es}$, we propose a channel estimation algorithm based on PARAFAC decomposition for estimating the UE-RIS channel ${\widehat{h}}_{r,es}$.

4.1. Uplink Pilot Transmission and Problem Formulation

4.1.1. Uplink Pilot Transmission

In the second stage, we consider an uplink transmission where the user transmits signals, and the semi-passive elements turn off the received signal mode and engage with the passive elements to reflect the signals, reflecting the user-transmitted pilot signals to the BS. To reduce the impact of inter-signal interference on the estimation performance, we use orthogonal pilot $x_k\in {\mathbb{C}}^{K\times 1}$ to perform channel estimation, where the orthogonal pilot satisfies

$$x_{k_1}{x}_{k2}^H=\left\{\begin{array}{cc}{P}_{MS},\hfill & \mathrm{if}{k}_1=k_2\hfill \\ 0,\hfill & \mathrm{otherwise},\hfill \end{array}\right.$$

(28)

where $P_{MS}$ is the transmit power of the UE; ${(.)}^H$ is the Matrix conjugate transpose. The uplink pilot transmission frame consists of t sub-frames, and each sub-frame lasts for K time slots. In the t-th sub-frame, the reflection coefficient vector at RIS is $\mathsf{\Theta }\in {\mathbb{C}}^{N\times 1}$. We keep the RIS phase shift constant for K time slots in each sub-frame. In one sub-frame, the signal received by the BS is modeled as

$$y_k=W_1^H{G}_{es}diag\left(\mathsf{\Theta }\left(1\right)\right)h_{r,k}{x}_k+n_k,k=1,2\dots ,K,$$

(29)

where $W\in {\mathbb{C}}^{M\times M}$ is the composite matrix at the BS, $h_{r,k}\in {\mathbb{C}}^{N\times K}$ is the first k UE-RIS equivalent channel of the UE-RIS, $G_{es}\in {\mathbb{C}}^{N\times M}$ is the BS-RIS channel estimated in the first stage, $n_k$ is AWGN with $\mathcal{C}\mathcal{N}\left(0,{\sigma }_n^2\right)$, and $\mathsf{\Theta }\left(1\right)={\left[{\beta }_1{e}^{j{\theta }_1\left(1\right)},\dots ,{\beta }_n{e}^{j{\theta }_n\left(1\right)},\dots ,{\beta }_N{e}^{j{\theta }_N\left(1\right)}\right]}^T$$\in {\mathbb{C}}^{N\times 1}$ is a diagonal matrix of the reflection coefficients of the RIS with phase shift, the ${\beta }_n\in [0,1](n=1,2,\dots ,N)$ is the amplitude reflection coefficient, and ${\theta }_n\in [0,2\pi ](n=1,2,\dots ,N)$ is the modulation phase. Then the signal received at t sub-frames is

$${\overline{Y}}_t=W_t^H{G}_{es}diag\left(\mathsf{\Theta }\left(t\right)\right){\left(X{h}_{r,k}{}^T\right)}^T+N_{t,k},$$

(30)

where ${\overline{Y}}_t={[y_{k,1},y_{k,2},\dots y_{k,t}]}^T$,$X={[x_{k,1},x_{k,2},\dots x_{k,t}]}^T$, and $N_k={[n_{k,1},n_{k,2},\dots n_{k,t}]}^T$. We consider $\mathsf{\Theta }$ for a DFT matrix satisfying $\mathsf{\Theta }*{\mathsf{\Theta }}^H=tI$ in this paper.

4.1.2. Problem Formulation

The above problem can be translated into an optimization problem

$${\widehat{h}}_{r,es}=\underset{h_{r,k}}{argmin}{∥{\overline{Y}}_t-W_t^H{G}_{es}\mathsf{\Theta }{Z}^T∥}_F^2,$$

(31)

where $Z=X{h}_{r,k}{}^T$. Since the UE-RIS is a low-dimensional channel and the LS algorithm has good results in dealing with low-dimensional channels [22], then the solution of Equation (31) obtained by least squares is ${\widehat{h}}_{r,\theta s}={\left(A^{H}A\right)}^{-1}{A}^H{\overline{Y}}_k$, where $\mathrm{A}=W_t^H{G}_{es}$. However, the channel $G_{es}$ is a sparse channel in mmWave systems, i.e., $Rank\left(G_{es}\right)=L<min(M,N)$ and $Rank\left(A^{H}A\right)<min(M,N)$, then the matrix inverse does non-existent when using the LS algorithm for channel estimation, which leads to inaccurate estimation results. To solve this problem, we propose a least squares method based on PARAFAC decomposition to estimate the UE-RIS channel in the next subsection, as shown in steps 17 to 18 of Algorithm 1.

4.2. The LS Algorithm Based on PARAFAC Decomposition

Based on the PARAFAC decomposition, the matrix ${\overline{Y}}_k$ can be considered as a three-way tensor $Y\in T\times M\times K$ of the k the first positive matrix slice [23], then the noise-free received signal tensor can be obtained ${\overline{Y}}_k$ of the $(\ell ,t,k)$ term

$${\left[{\overline{Y}}_t\right]}_{\ell ,t,k}=\sum _{n=1}^N{g}_{\ell ,n}{z}_{t,n}{s}_{k,n},$$

(32)

where $g_{l,n}={\left[W_t^H{G}_{es}\right]}_{l,n}={\left[G\right]}_{l,n}$, $z_{t,n}={\left[Z\right]}_{t,n}$, $s_{k,n}={\left[S\right]}_{k,n}$. Exploiting the trilinearity of the PARAFAC decomposition, Equation (28) can be transformed as [24]

$${\overline{Y}}_t=X{h}_{r,k}{}^T{\left(\mathsf{\Theta }\diamond \left(W_t^H{G}_{es}\right)\right)}^T+B_t,$$

(33)

where ⋄ is the Khatri–Rao product, $B_t$ is the tensor unfolding of the noise. Since the $G_{\mathrm{es}}$ was estimated in the previous stage, solving the above problem can be translated into the following function

$${\widehat{h}}_{r,es}=\underset{h_{r,k}}{argmin}{∥{\overline{Y}}_t-X{h}_{r,k}{}^T{\left(\mathsf{\Theta }\diamond \left(W_t^H{G}_{es}\right)\right)}^T∥}_F^2,$$

(34)

where ${\widehat{h}}_{r,es}$ is the estimated value of the UE-RIS channel. The solution is given by

$${\widehat{h}}_{r,es}=X^{†}{\overline{Y}}_t{\left[{\left(\mathsf{\Theta }\diamond \left(W_t^H{G}_{es}\right)\right)}^T\right]}^{†},$$

(35)

where † is the pseudo-inverse of the matrix. To ensure the accuracy of the estimated solution, it is necessary to satisfy $rank\left(\mathsf{\Theta }\right)+rank\left(W_t^H{G}_{es}\right)\ge N+1$. In other words, $t\ge N-L_{r,k}+1$. The detailed proof process is in the literature [24] and will not be overly described here. The proposed two-stage channel estimation algorithm is summarized in Figure 3 and Figure 4.

4.3. Pilot Overhead and Computational Complexity

4.3.1. Pilot Overhead

In the second stage of estimating the time-varying channel UE-RIS, the RIS-UE channel has $NK$ coefficients to be estimated. Since the $\mathrm{BS}$ is estimated in a sub-frame of K time slot to obtain $MK$ coefficients, we need at least $N-L_{r,k}+1$ sub-frames. Therefore, the channel estimation for the second stage has a pilot overhead of ${\tau }_2=T_2=(N-L_{r,k}+1)K$.

4.3.2. Computational Complexity

In the second stage of estimating the UE-RIS channel estimation, the algorithm complexity is determined by the LS channel estimation algorithm based on the PARAFAC decomposition, which is $\mathcal{O}\left(MT(2M+1)+N^2(M+1)+\right.MN)$.

5. Simulation Results

In this section, we present simulation results to verify the effectiveness of the proposed channel estimation algorithm. In the simulations, we use normalized mean square error (NMSE) and average spectrum efficiency (ASE) to evaluate the channel estimation algorithm performance, whose expression is

$$NMSE\left(G_{es}\right)=10{log}_{10}\left(\frac{1}{R}\sum _{r=1}^R\frac{\parallel G_{es}{-G\parallel }_F^2}{{\parallel G\parallel }_F^2}\right),$$

(36)

$$NMSE\left({\widehat{h}}_{r,k}\right)=10{log}_{10}\left(\frac{1}{R}\sum _{r=1}^R\sum _{k=1}^K\frac{{∥{\widehat{h}}_{r,k}-h_{r,k}∥}_F^2}{{∥h_{r,k}∥}_F^2}\right),$$

(37)

$$\begin{array}{c}\hfill NMSE\left(\overline{H}\right)=10{log}_{10}\left(\frac{1}{R}\sum _{r=1}^R\frac{\parallel \overline{H}{-H\parallel }_F^2}{{\parallel H\parallel }_F^2}\right),\end{array}$$

(38)

$$ASE=E\left\{{log}_2\left(det\left(1+{\left(N_R{N}_T({\sigma }^2+NMSE)\right)}^{-1}H{H}^H\right)\right)\right\},$$

(39)

where R is the number of Monte Carlo simulations and $\widehat{G}$ is the result of the r-th the BS-RIS channel estimation, ${\widehat{h}}_{r,k}$ is the result of the r-th UE-RIS channel estimation, $\overline{H}=\widehat{G}diag\left({\widehat{h}}_{r,k}\right)$ is the r-th estimated cascaded channel, and $H=Gdiag\left(h_{r,k}\right)$ is the true channel.

5.1. Parameter Setting and Simulation Analysis

The simulation parameters are set as follows: $M=64$, $N=256$, $K=10$, $L=3$, $L_{r,k}=6$, $Q=64$, the AoA and AoD parameters are uniformly generated from $[-\pi /2,\pi /2]$, and the number of Monte Carlo simulations is 2000. The signal-to-noise ratio (SNR) is defined by $\mathrm{SNR}=\frac{1}{{\sigma }_n^2}$, where ${\sigma }_n^2$ is the noise variance. MATLAB R2019a is used for simulation. To verify the superiority of the proposed algorithm, we compare the following algorithms:

LS: The LS estimator from the literature [8] to estimate the channel.

MMSE: Based on the LS channel estimation algorithm, using the MMSE estimator from the literature Equation (10) to estimate the channel.

CS: Solving the non-convex optimization problem of Equation (12) using the CS algorithm which can be found in the literature [15].

Oracle LS: By assuming that the exact AOA and AOD are known to the RIS or BS and that the received signal is estimated to be a known common LS, the subspace is upper bounded with respect to the performance of the solution, which cannot realistically be achieved.

Figure 5 shows the number of semi-passive elements $N_{semi}=\sqrt{N_0}$ versus NMSE. We choose the UPAs with the sizes of $16\times 16$ and $24\times 24$ for comparison, and it can be seen that as the number of RIS elements increases, the estimation performance improves slightly, but the estimation performance decreases slightly as the number of RIS elements increases, which is due to the decrease in the proportion of semi-passive elements. It suggests that the number of semi-passive parts can be raised properly in actual applications to improve estimation performance.When the semi-passive elements are 64 and $SNR=0,10,20$ dB, the estimation accuracy of the proposed algorithm can reach ${10}^{-1.5}$, ${10}^{-3}$, and ${10}^{-4}$ orders of magnitude, which shows the superiority of the proposed algorithm in the channel estimation accuracy. Considering the trade-off between complexity and the performance affected by the number of semi-passive elements, we take the number of semi-passive elements to be 64 in the later simulations.

Figure 6 shows the pilot overhead versus NMSE for the estimated BS-RIS channel with semi-passive elements of 16 and 64. We see that the scheme proposed in this paper performs better than the CS-based scheme for the same pilot overhead conditions. Moreover, the difference between the two grows as the training time increases. The CS-based method does not improve with increasing pilot overhead, mainly due to the limitation of the per-discrete AOD/AOA resolution. When the pilot overhead is 100, the estimation accuracy of the algorithm proposed in this paper can reach ${10}^{-3.5}$, ${10}^{-4.5}$ orders of magnitude; however, the estimation accuracy of CS is ${10}^0$, ${10}^0$ orders of magnitude, which shows the superiority of the algorithm proposed in this paper, which can still achieve a more desirable channel estimation with less pilot. Furthermore, the error of the scheme proposed in this paper gradually decreases as the number of semi-passive elements increases from 16 to 64, which is due to the increase in the proportion of semi-passive elements.

Figure 7 and Figure 8 show the relationship between NMSE and signal-to-noise ratio for the BS-RIS and the UE-RIS channel estimation. We compare the NMSE of the proposed channel estimation with the LS, CS, and Oracle LS algorithms. Among them, the LS and the MMSE algorithms are carried out under the condition that all elements of the RIS are semi-passive elements. As can be seen from the figure, except for the CS algorithm, the rest of the algorithms decrease with the increase of SNR because the compressed sensing algorithm adopts a robust algorithm and has low sensitivity to noise. Because the real channel information is not necessarily in the preset dictionary, quantization errors are caused, resulting in poor estimation performance. The LS algorithm is sensitive to noise. Although the estimation decreases with the increase of SNR, the overall estimation performance is poor. When $SNR=-15$ dB, the proposed algorithm can reach the estimation accuracy of ${10}^{-0.5}$ and ${10}^{-1}$ orders of magnitude for the BS-RIS channel and the RIS-UE channel, respectively, which has obvious advantages over the LS, MMSE, and CS algorithms. When $SNR=25$ dB, the proposed algorithm can reach ${10}^{-5}$ orders of magnitude, respectively, and the estimation accuracy has obvious advantages over the other three algorithms. In general, the proposed algorithm has lower NMSE than the OMP, LS, and MMSE algorithms. However, the Oracle LS algorithm is more accurate for estimation because RIS uses the LS algorithm to estimate the perfect AOD and AOA. The algorithm is an ideal situation and has superior performance, but it cannot be applied in practice. Figure 9 shows the relationship between NMSE and SNR for cascaded channel estimation. It can be seen that for the cascaded channel, the estimation accuracy of the proposed division algorithm can reach ${10}^{-0.5}$, ${10}^{-5}$ orders of magnitude for $SNR=-15$ dB and $SNR=25$ dB, respectively, which are higher than those of the CS, LS, and MMSE algorithms for both high and low SNR cases, and is consistent with the results in Figure 6 and Figure 7. Therefore, we can see the superiority of the performance of the proposed algorithm.

Figure 10 shows the ASE of several methods for various SNRs, with the ideal CSI case serving as an upper bound. The proposed algorithm achieves a better spectral efficiency (SE) than other algorithms and is close to the ideal value. This is because the suggested channel estimation algorithm is better than other algorithms in terms of estimate accuracy, as shown in Figure 9. The CS algorithm grows slowly at SNRs of 5 to 10, which is caused by the stabilization of NMSE performance. When the SNR is 5, the performance of the proposed algorithm outperforms the CS, LS, and MMSE algorithms by 90.2%, 86.3%, and 10.5%, respectively. In addition, we present the SE performance without RIS, which verifies that the presence of RIS does improve the SE performance of mmWave communication because of the additional beamforming gain due to RIS. As a result, it can be demonstrated that the proposed scheme may successfully improve the SE. In addition, the quantified data when SNR = 10 dB are given in

Table 2. The average spectrum efficiency of different algorithms when SNR = 5 dB.

Algorithm	Average Spectrum Efficiency (bps/Hz)
Perfect Channel	13.089
Proposed	12.911
LS	5.725
CS	8.053
MMSE	11.695

.

5.2. Pilot Overhead and Computational Complexity Analysis

5.2.1. Pilot Overhead

To calculate the minimum pilot overhead, i.e., the average pilot overhead in the $T_L$ time period, the average pilot overhead is $\tau =\frac{T_1}{T_L}{\tau }_1=\frac{T_1}{T_L}M\left[\frac{N}{N_0}\right]$, due to $T_L>>T_1$, then $\frac{T_1}{T_L}{\tau }_1\approx 0$, the minimum pilot overhead of the proposed algorithm is $\tau +{\tau }_2=(N-L_{r,k}+1)K$. In

Table 3. Pilot overhead comparison of different channel estimation schemes.

Algorithm	Minimum Pilot Overhead
Proposed	$(N-L_{r,k}+1)K$
LS	$KN$
CS	$K\left[\frac{8L{L}_{r,k}-2}{M}\right]$
MMSE	$KN$

, we compare the minimum pilot overhead of the proposed two-stage channel estimation with the LS, MMSE, and CS algorithms. As shown in Figure 11, Figure 12 and Figure 13, we plot the relationship between the pilot overhead and the number of BS antennas, the number of RIS elements, and the number of UEs. From the figures, we can see that the pilot overhead is lower than the LS and MMSE algorithms. Specifically, when the number of BS antennas is more than 30, the pilot overhead of the proposed algorithm is 472, which is lower than the LS algorithm and the MMSE algorithm. When the number of RIS elements as well as the number of UEs increases, the pilot overhead of the proposed algorithm, the LS algorithm, and the MMSE algorithm also increases, but the proposed pilot overhead is lower than these two algorithms. In addition, the CS algorithm uses the sparsity of the channel, and the pilot overhead is only related to the pre-defined grid matrix. Although the proposed algorithm has a pilot overhead higher than the CS algorithm, the proposed algorithm is more accurate than the CS algorithms in estimating. This also validates the superiority of the proposed algorithm.

5.2.2. Computational Complexity

In the first stage of the BS-RIS channel estimation, the algorithmic complexity of the super-resolution algorithm for solving the optimization problem is $\mathcal{O}\left(MT\left(M+N_0\right)N_0^2{M}^2\right)$, and the complexity is further reduced by the SVD algorithm to $\mathcal{O}\left(MT\left(M+N_0\right)N_0^2{L}^2\right)$. In the second stage of estimating the UE-RIS channel estimation, the algorithmic complexity is determined by the LS channel estimation algorithm based on the PARAFAC decomposition with complexity $\mathcal{O}\left(MT(2M+1)+N^2(M+1)+MN\right)$. In addition, we give the computational complexity of the compression-aware CS algorithm $\mathcal{O}\left(2T{N}_{0}M\right)$; the computational complexity of the LS algorithm is $\mathcal{O}\left(T{\left(NM\right)}^2\right)$.

6. Conclusions and Future Work

6.1. Conclusions

Aiming at the channel estimation problem of RIS-assisted MIMO systems, this paper proposes a two-stage channel estimation scheme. Firstly, we propose a two-stage channel estimation protocol to reduce the pilot overhead. The core idea of the protocol is to introduce semi-passive components to assist channel estimation and to use the coherent time difference between the BS-RIS channel and the RIS-UE channel to reduce the pilot overhead on a long time scale. Then, to estimate the BS-RIS channel accurately, we propose a super-resolution channel estimation algorithm. The algorithm moves the estimated AOA/AOD value to close to the real values through gradient descent and then obtains the accurate channel gain values based on the closed solution of the optimization problem. The algorithm avoids off-grid effects and power leakage and realizes the super-resolution estimation of the BS-RIS channels. Finally, to estimate the RIS-UE channel accurately, we propose an LS algorithm based on PARAFAC decomposition. The algorithm effectively utilizes the Khatri–Rao structure of the combined channel matrix to improve the estimation accuracy. Simulation results show that the proposed scheme can achieve higher accuracy than LS, CS, and MMSE schemes with less lead time overhead, and the spectral efficiency is improved by at least 10% compared to the other three algorithms. It also proves the superiority of the proposed scheme.

6.2. Future Work

The scheme proposed in this research can be utilized to greatly improve the performance of single RIS-assisted mmWave systems. We will look into the channel estimation challenge for multi-RIS-assisted mmWave systems in the future. Furthermore, in the entire application system, beamforming and RIS reflection coefficient design are critical. As a result, we will optimize the design of the beamforming vector and the RIS reflection coefficient matrix in tandem to improve system performance in the future study.