LEO Satellite-Enabled Grant-Free Random Access with MIMO-OTFS

We consider a LEO satellite-enabled grant-free random access system with MIMO-OTFS, where $U$ single-antenna devices adopting OTFS modulation intend to communicate with a LEO satellite. The satellite is equipped with $N_{a}=N_{y}\times N_{z}$ uniform planar array (UPA) antennas and works with the regenerative payload, which allows it for on-board processing of baseband signals. In a given time interval, only $U_{a}$ devices are active and transmit signal to satellite utilizing the same time-frequency resources. We denote $\lambda_{u}$ as the activity indicator of the $u$-th device, with $\lambda_{u}=1$ if active and $\lambda_{u}=0$ otherwise. We follow the recommendations of 3GPP[9] and assume that the satellite precompensates Doppler shift for the center of the beam and broadcasts a common delay to all devices. The residual delay and Doppler shift (called differential delay and Doppler shift in 3GPP terminology) are handled by satellite in the uplink transmissions.

To exploit the sparsity in the angular domain, the two-dimensional discrete Fourier transform (2D-DFT) is applied for the channel of $u$-th device along the space dimension, and then the delay-Doppler-angle domain channel is given by

where $a_{y}=0,\dots,N_{y}-1$ and $a_{z}=0,\dots,N_{z}-1$ are indexes along angular domain and $\Pi_{N}(x)\triangleq\frac{1}{N}\sum_{i=0}^{N-1}e^{-\bar{\jmath}2\pi\frac{x}{N}i}$. From (II-C), $H_{u,a_{y},a_{z}}^{\mathrm{DDA}}[k,l_{u,i},l]^{\mathrm{DDA}}$ has dominant elements only if $k\approx NT_{\text{sym}}(\nu_{u,i})_{1/T_{\text{sym}}}$, $l_{u,i}\approx M\Delta f(\tau_{u,i})_{T}$, $a_{y}\approx N_{y}\Omega_{y_{u}}/2$, and $a_{z}\approx N_{z}\Omega_{z_{u}}/2$. Therefore, the channel in the delay-Doppler-angle domain has the 3D-structured sparsity[5].

The received symbols from all the devices in the delay-Doppler-angle domain is given by

where $Z^{\mathrm{DDA}}$ is the noise in the delay-Doppler-angle domain. To facilitate the following analysis, we rewrite (II-C) into the matrix form for each $l$ as

where the elements of $\mathbf{Z}_{l}$ are independent zero-mean Gaussian noises with variance $\sigma^{2}$; the $(a_{z}+1+N_{z}a_{y})$-th column of $\mathbf{Y}_{l}$ is $\mathbf{Y}_{a_{y},a_{z}}^{\mathrm{DDA}}[:,l]$; $\mathbf{X}_{l}=[\mathbf{X}_{0,l},\dots,\mathbf{X}_{U-1,l}]\in\mathrm{C}^{N\times UMN}$, where $\mathbf{X}_{u,l}\in\mathrm{C}^{N\times MN}$ is the pilot matrix of the $u$-th device in delay dimension $l$. Note that $\mathbf{X}_{u,l}$ is a block circulant matrix due to the 2D circular convolution in (II-C), and its sub-matrix $\mathbf{X}_{u,l,l^{\prime}}$, $l^{\prime}=0,\dots,M-1$ is circulant matrix, formed by $\left[X_{u}^{\mathrm{DD}}[0,(l-l^{\prime})_{M}],\dots,X_{u}^{\mathrm{DD}}[-1,(l-l^{\prime})_{M}]\right]$. $\Lambda=\mathrm{diag}([\lambda_{0}\mathbf{1}_{MN},\cdots,\lambda_{U-1}\mathbf{1}_{MN}])$ is a diagonal matrix composed of the activity indicator of all the devices and $\mathbf{1}_{MN}\in R^{MN}$ is a all-one vector; Finally, the $(a_{z}+1+N_{z}a_{y})$-th column of $\tilde{\mathbf{H}}$ is given by

where $\mathrm{vec}(\cdot)$ denotes the vectorization of a matrix. We denote the channel matrix $\Lambda\tilde{\mathbf{H}_{l}}$ as $\mathbf{H}_{l}$. Then, (10) can be written as

Therefore, the joint device activity detection and channel estimation in this scenario turns into a sparse signal reconstruction problem, and can be separately and parallelly handled along delay dimension for each $\mathbf{H}_{l}$. In addition, the 2D burst block sparsity[7] appears in $\mathbf{H}_{l}$, which is different from those in the conventional grant-free random access systems[3, 4] and can be utilized to improve the estimation performance further.

In this work, we estimate $\mathbf{H}_{l}$ adopting the minimum mean square error (MMSE) rule, represented as

Based on the estimation, the active device can be detected by comparing the energy of each device’s channel with a pre-defined threshold, i.e.,

where $\hat{H}_{l,i,j}$ is the $(i,j)$-th element of $\hat{\mathbf{H}}_{l}$. Note that since each $\mathbf{H}_{l}$ can be estimated separately, we omit the index $l$ in the following content.

Since channel matrix $\mathbf{H}$ has the 2D burst block sparsity, i.e., its non-zero elements are statistically dependent along the row and column, we consider the cross-correlation-based Gaussian prior for the $(i,j)$-th element of $\mathbf{H}$, given by

where $\mathrm{\mathbf{A}}$ is the precision matrix whose $(i,j)$-th element $\alpha_{i,j}$ is the local precision hyperparameter of $h_{i,j}$; $\mathrm{\mathbf{B}}\in[0,1]^{D_{x}\times D_{y}}$ is the kernel and its $(p,q)$-th element $\beta_{p,q}$ represents the relevence between $h_{i,j}$ and $h_{i+p,j+q}$. We denote $\theta_{i,j}$ as the prior variance of $h_{i,j}$. From (18), $\theta_{i,j}$ not only depends on local hyperparameter $\alpha_{i,j}$ of $h_{i,j}$, but also on the hyperparameters of its neighborhoods along the row and column of $\mathbf{H}$. Therefore, such coupled prior has the potential to capture the unknown 2D sparsity structure in $\mathbf{H}$. Then, following the conventional SBL, we adopt the Gamma distribution as the hierarchical prior for the local precision hyperparameter $\mathrm{\mathbf{A}}$, i.e.,

where $\Gamma(c)\triangleq\int_{0}^{\infty}t^{c-1}e^{-t}dt$. Note that this hierarchical Bayesian modeling encourages the sparseness in the estimation since the overall prior of $\mathbf{H}$ (i.e., integrate out $\alpha_{i,j}$) is a Student-$t$ distribution that is sharply peaked at zero. Based on the above prior distribution, it is possible to find the analytical solution of (16). However, the matrix inversion will inevitably be involved, which leads to the high computational complexity. Fortunately, the prior and the likelihood function of $\mathbf{H}$ can be fully factorized, leading to the following posterior distribution, i.e.,

where $\mathbf{y}_{j}$ and $\mathbf{h}_{j}$ are the $j$-th column of $\mathbf{Y}$ and $\mathbf{H}$, respectively, and $\mathbf{g}_{j}=\mathbf{X}\mathbf{h}_{j}$. Hence, we can resort to the GAMP algorithm to estimate the posterior mean and variance of the channel matrix with lower complexity.

The messages passed within the GAMP algorithm are approximated by Gaussian distributions, and then only the means and variances of the messages need to be calculated and preserved, which reduces the complexity significantly. We next outline the MMSE estimation of channel matrix $\mathbf{H}$ in the $t$-th iteration by following the GAMP algorithms[10]. The detailed derivations are omitted for brevity. Here, the hyperparameters {$\mathbf{A}$, $\mathbf{B}$, $\sigma$} are assumed to be known and fixed, and we will update them later.

In this subsection, we utilize the EM algorithm together with deep learning to learn the hyperparameters {$\mathbf{A}$, $\mathbf{B}$, $\sigma$}, where $\mathbf{A}$ and $\sigma$ are updated by the EM-based rule and $\mathbf{B}$ is learned by a convolutional neural network (CNN), as shown in Fig. 1. We then explain it in detail. Following the conventional SBL, we place a Gamma hyperprior over $\sigma^{-2}$, i.e.,

Based on the estimation in the $t$-th iteration, the precision matrix and the noise variance can be updated by

where the kernel $\mathrm{\mathbf{B}}$ is assumed to be known in this step. The optimization in (28) is difficult to get the analytic solution since the hyperparameters are entangled with each other due to the logarithm term. Although the gradient descent methods can be adopted to get the optimal solution, it involves higher computational complexity as compared with an analytical update rule. We derive a closed-form sub-optimal solution (see [7] for more detail) by examining the first derivative of the objective function with respect to $\alpha_{i,j}$ and by proper scaling, given as

where $\omega_{i,j}^{t}=\sum_{p=-D_{x}}^{D_{x}}\sum_{q=-D_{y}}^{D_{y}}\beta_{p,q}\mathrm{E}[\left|h_{i-p,j-q}\right|^{2}]$ and posterior second moment $\mathrm{E}[\left|h_{pq}\right|^{2}]=\left|\hat{h}^{t}_{p,q}\right|^{2}+\left|\hat{\gamma}^{t}_{p,q}\right|^{2}$, which are the outputs of the GAMP algorithm in the $t$-th iteration. Then, computing the first derivative of the objective function in (29) with respect to $\sigma^{2}$ and setting it equal to zero, we get

where

By examining (18) and (30), the update of precision matrix $\mathbf{A}$ is related to the 2D convolution between kernel $\mathbf{B}$ and the posterior second moment, and the update of prior variance $\theta_{i,j}$ is related to the cross-correlation between the kernel and precision matrix, i.e.,

The only difference between the 2D convolution and the cross-correlation is whether the kernel is flipped within the computations, which motivates us to utilize the CNN to learn kernel $\mathbf{B}$ and then update precision matrix $\mathbf{A}$. Here, the trainable filter of the convolutional layer can be seen as flipped kernel $\mathbf{B}$. Combining this idea with the EM update rule and the GAMP algorithm leads to the deep learning-based GAMP algorithm (DL-GAMP), and we summarize the computation process as follows.

As shown in Fig. 1, in the $t$-th iteration, the posterior second moment computed by the GAMP algorithm in the $(t-1)$-th iteration is firstly inputted to a convolutional layer and then the rectified linear unit (ReLU)[11] will ensure a positive output $\omega_{i,j}^{t}$. Next, precision $\alpha_{i,j}^{t+1}$ is updated according to (30), and then is inputted to another convolutional layer followed by the ReLU to get the update of prior variance $\theta_{i,j}^{t+1}$ (similar to (33)). Simultaneously, noise standard variance $\sigma^{t+1}$ is given by (31). Finally, the output of GAMP in this iteration {$\hat{h}_{i,j}^{t+1}$, $\hat{\gamma}_{i,j}^{t+1}$} and the update of {$\theta_{i,j}^{t+1}$, $\sigma^{t+1}$} computed by the CNN and the EM approach will be inputted to the next iteration of the DL-GAMP. The algorithm will run until the maximum iteration $T_{\text{max}}$ is reached, and then the active device is detected according to (17). Note that the filters in the two convolutional layers are not restricted to be identical since they will be learned from the training samples. Besides, the matrix multiplication in the GAMP algorithm dominates other computations, like the 2D convolution in CNN, leading to complexity $\mathcal{O}(UMN^{2}N_{y}N_{z})$ in each iteration when each channel matrix $\mathbf{H}_{l}$ is parallelly estimated by DL-GAMP. If $\mathbf{H}_{l}$ is estimated in a serial manner, the complexity will be $\mathcal{O}(UM^{2}N^{2}N_{y}N_{z})$ in each iteration.