A Physics-based and Data-driven Approach for Localized Statistical Channel Modeling

The technique of beamforming is exploited in the massive MIMO of 5G systems, yielding high-resolution multi-path channels in the angular domain. Thus, we consider the channel modeling of the massive MIMO downlink communication system with beamforming, where the signal is transmitted from a base station to the user equipment through several propagation paths, as shown in Fig. 1.

Similar to the 3GPP’s technical report [11], we focus on the tilt angle of departure (AoD), azimuth AoD and channel gain from the base station to the user equipment with one receive antenna, ignoring the arrival angle and delay. Suppose the uniform rectangular array of the base station contains $N_{T}=N_{x}\times N_{y}$ antennas. The CIR of antenna $(x,y)$ from the base station to the user equipment is given by

The definition of notation in (II-A) is summarized in Table I. The tilt AoD and azimuth AoD over the free space are discretized into $N_{V}$ angles and $N_{H}$ angles, respectively. If there does not exist a path in the angles $\left(\theta_{i},\varphi_{j}\right)$, the corresponding channel gain $\alpha_{i,j}(t)$ is zero, otherwise $\alpha_{i,j}(t)>0$ for the path $\left(\theta_{i},\varphi_{j}\right)$ in the channel. It is readily known that the channel gain ${\alpha_{i,j}}(t)$ is sparse, and contains multiple propagation paths in the angular domain.

The channel gain ${\alpha_{i,j}}(t)$ consists of path loss and the shadowing effect. Path loss is determined by the physical environment (distance, carrier frequency, buildings) which is assumed to be relatively static, while shadowing is caused by the obstacles and is usually modeled as the log-normal distribution. Thus, we assume ${\alpha_{i,j}}(t)$ follows the log-normal distribution, with its mean representing path loss and its covariance representing the shadowing effect. $\omega_{i,j}(t)$ is the random phase error between different angles caused by the reflection, diffraction and scattering effect of electromagnetic waves, while $\omega_{x,y}(t)$ is the random phase error of different antennas caused by the imperfect hardware of the antenna array. Inspired by [11], we assume $\omega_{i,j}(t)$ follows the uniform distribution between $-\pi$ and $\pi$, and $\omega_{x,y}(t)$ follows the Gaussian distribution with zero mean and covariance $\sigma^{2}$.

The channel measurement of our proposed LSCM is the RSRP measured from multiple beams. In the downlink of 5G cellular systems, the RSRP can be measured from synchronization signals block (SSB) beams or channel state information-reference signal (CSI-RS) beams. Denote the precoding matrix of the $m_{\rm th}$ beam as $\bm{W}^{(m)}\in\mathbb{C}^{N_{x}\times N_{y}}$, whose size is the same as that of the antenna array. The entry of $\bm{W}^{(m)}$ is $\left(w_{x,y}^{(m)}\right)=\left(e^{j\phi_{x,y}^{(m)}}\right)$. Denote the CIR matrix from the base station to the user equipment as $\bm{H}\in\mathbb{C}^{N_{x}\times N_{y}}$, where the entry of $\bm{H}$ is $h_{x,y}(t)$ in (II-A). The RSRP of the $m_{\rm th}$ beam at time $t$ is defined as

where $P$ denotes the transmit power. The quality of the channel can be represented through ${rsrp}_{m}(t)$, and a larger value of ${rsrp}_{m}(t)$ reflects better channel quality.

Figure 2 shows an example of channel measurement of RSRP measured from five beams. These five beam patterns are generated according to different beamforming weights, and the directions of main lobe of these five beams are different. The receiver at the target grid is highlighted as a blue point. The target grid can receive five RSRP measurements, corresponding to the received power of five beams. To be specific, the strongest value of RSRP measurements is -59.61 dB in Fig. 2(d), because its main lobe’s direction is $\left[\bar{\theta},\bar{\varphi}\right]=[0,30]$, which is close to the angle of departure of the propagation path $\left[\widehat{\theta},\widehat{\varphi}\right]=[-2.84,23.82]$.

As for ${rsrp}_{m}(t)$, it is a random variable inferred from (2), since $h_{x,y}(t)$ contains three independent random variables, i.e., $\omega_{x,y}(t)$, $\omega_{i,j}(t)$ and ${\alpha_{i,j}}(t)$. Theorem 1 shows the relationship between the first-order statistics of ${rsrp}_{m}(t)$ and ${\alpha_{i,j}}(t)$, i.e., their expectation $\mathbb{E}\left({rsrp}_{m}(t)\right)$ and $\mathbb{E}\left({\alpha_{i,j}}(t)\right)$, implying the statistical relationship of the RSRP and angular information.

The statistical relationship between the expectation of RSRP measurements ${\rm{RSRP}}_{m}$ and the first-order statistics of channel gain $X_{i,j}$ is connected by the coefficient ${\rm A}_{i,j}^{(m)}$, as revealed by Theorem 1. According to this finding, it is possible to extract the angular information from RSRP measurements if we can calculate $X_{i,j}$ from ${\rm{RSRP}}_{m}$. The APS can be expressed by $X_{i,j}$, where the subscripts $i$ and $j$ indicate the tilt angle $\theta_{i}$ and azimuth angle $\varphi_{j}$, respectively.

Based on (3), we reshape $X_{i,j}$ with the subscripts $i=1,\cdots,N_{V},j=1,\cdots,N_{H}$ as a vector, i.e.,

where the length of ${\bf x}$ is $N=N_{V}\times N_{H}$. If the path $\left(\theta_{i},\varphi_{j}\right)$ does not exist, then $\mathbb{E}\left({\alpha_{i,j}}(t)\right)=0$. As there are only several propagation paths from the base station to the user equipment, the channel vector ${\bf x}$ is sparse.

Suppose there are $M$ measured beams in total, and define ${\bf y}\in\mathbb{R}^{M\times 1}$ as the vector of the expectation of RSRP measurements from different beams, i.e.,

Define the vector ${\bf a}^{(m)}\in\mathbb{R}^{N\times 1}$ as

and the coefficient matrix ${\bf A}\in\mathbb{R}^{M\times N}$ is

Then we can recast (3) as

To construct the localized statistical channel model, we want to recover the sparse channel vector $\bf x$ from the expectation of RSRP measurements $\bf y$ and the coefficient matrix $\bf A$. According to (10), as the expectation of channel gain in terms of different angles ${x}_{n}\geq 0,\ n=1,\cdots,N$, is non-negative, we can formulate the optimization problem as

where $K$ is the maximum number of non-zero entries, representing the maximum number of propagation paths.

The OMP is an iterative algorithm, which selects one column into the set of indexes of non-zero entries $\mathcal{S}$ at each iteration. The selection criteria is to select the column that is best correlated with the residual of current iteration. Problem (11) can be solved by the non-negative OMP (NNOMP), as shown in Algorithm 1. The biggest difference between NNOMP and OMP lies in Step 5 of Algorithm 1, where the subproblem of non-negative least squares (NNLS) is solved instead of the unconstrained least squares (ULS). After solving the NNLS subproblem, the entries of $\bf x$ in the complement set $\mathcal{S}^{c}$ are forced zero to ensure the sparsity. A fast implementation of NNOMP is proposed in [12] based on the active-set algorithm. The role of $\max\left({\bf A}^{T}{\bf r}_{k}\right)<0$ is also explained in [12] by using the Karush-Kuhn-Tucker (KKT) condition. The condition ${\bf a}_{n}^{T}{\bf r}_{k}>0$ means that the $n_{\rm th}$ column is the descending column, which can be added into the set $\mathcal{S}$ and reduce the residual.

However, in our case, because our coefficient matrix $\bf A$ is not column-normalized, this sparse recovery problem is intractable. The magnitude of each column varies significantly, which influences the accuracy of finding the correct column that is best correlated with the residual. The NNOMP algorithm prefers the large-magnitude columns of coefficient matrix $\bf A$, resulting in a low probability to select the correct index of non-zero entries. For example, although there exists one column that is parallel to the residual, the NNOMP algorithm will still select the wrong column with large magnitude.

To deal with this problem, we design the WNOMP algorithm, which is shown in Algorithm 2. The WNOMP algorithm can strike a good balance between the effect of column magnitude and correlation with the residual by setting the dynamic weight for each selection. The weight of magnitude should be smaller than that of the correlation. If we define the column-wise normalized coefficient matrix as

where ${\bf a}_{n}\in\mathbb{R}^{M\times 1}$ is the $n_{\rm th}$ column vector of coefficient matrix $\bf A$. The dynamic weight $\lambda_{k}$ can be calculated as

The dynamic weight $\lambda_{k}$ plays an important role in the selection of Step 4 of Algorithm 2. An intuitive idea to reduce the influence of column is to use the column-normalized matrix $\widehat{\bf A}$, i.e., the selection in Step 4 of Algorithm 2 is modified as $i=\arg\max_{n}\left(\left({{\bf a}_{n}/\|{\bf a}_{n}\|_{2}}\right)^{\rm T}{\bf r}_{k}\right)$. However, this modification is not good, because the propagation paths are more likely to exist in the angles with large antenna gain of the transmit antenna array in the base station, which means the large-magnitude column corresponds to the non-zero entries with a large probability. As for the selection $i=\arg\max_{n}\left({{\bf a}_{n}}^{\rm T}{\bf r}_{k}\right)$, it can be viewed as the multiplication of the modification case and $\|{\bf a}_{n}\|_{2}$, where the correlation with residual contributes little to the selection because the value of $\|{\bf a}_{n}\|_{2}$ is much larger than $\left({{\bf a}_{n}/\|{\bf a}_{n}\|_{2}}\right)^{\rm T}{\bf r}_{k}$. As a result, the WNOMP considers these two terms separately, and use $\lambda_{k}$ to balance them. The Step 7 of Algorithm 2 compresses the set $\mathcal{S}$ by updating it with the support of $\bf x$, because some components of $\mathcal{S}$ may vanish due to the activation of the non-negative constraints.