Integrated Location Sensing and Communication for Ultra-Massive MIMO With Hybrid-Field Beam-Squint Effect

Hybrid beamforming architecture has been widely considered for massive MIMO and UM-MIMO systems due to their low hardware cost achieved by using a much smaller number of radio frequency (RF) chains than antenna elements. [18]. Nevertheless, this architecture significantly decreases the effective baseband measurement dimension of received signals, bringing great challenges to full-dimensional channel state information (CSI) acquisition. To this end, compressive sensing (CS) technique, as an effective tool to recover the sparse signals via underdetermined measurements, has attracted extensive attention for channel estimation in mmWave and THz bands [18, 19, 20]. Under the assumption of the far-field planar wavefront, the mmWave massive MIMO channel exhibits an obvious sparsity pattern in the angular domain. Therefore, the channel estimation can be formulated as a larger number of antennasal sparse signal recovery problem, and solved resorting to the off-the-shelf CS-based algorithms [18, 21, 22]. However, when extended to UM-MIMO systems, the far-field approximation no longer holds, leading to a non-negligible energy dispersion and deteriorating its inherent sparsity. For this purpose, the authors of [23] introduced a subarray piecewise approximate far-field-based sparse channel estimation scheme, which handled the channels of different subarrays independently with the far-field approximation. However, it overlooked the geometrical correlations of different subarrays, thus increasing the parameters to be estimated, and only performing well at a moderate medium antenna array size. Furthermore, a more complex Fresnel approximation was developed based on second-order Taylor expansion, and [25, 24] proposed higher-dimensional projection matrices to enable sparse channel estimation at the cost of increased computational complexity. This is due to the reason that candidate columns steering vectors of projection matrices were determined by both the angular and the range parameters.

Most of the aforementioned schemes have focused on leveraging the structured common sparsity of the spatial-frequency channel to enhance estimation performance, yet neglecting the beam-squint effect. To address this issue, [26] revealed the bilinear pattern of HFBS, which implies that the sparse support set of channels in both the angle and the range domains can be regarded as a linear function against frequency, and developed a bilinear pattern based structured CS algorithm for channel estimation. Additionally, [9, 27, 28] developed a frequency-dependent projection matrix to correct the beams-squint effect and maintained the spatial-frequency consistency in the corresponding transform domain. Besides, for the problem of dispersed beams caused by the beam-squint effect, a common solution is to employ true-time-delay (TTD) to replace phase shifters to generate frequency-dependent phase shift but inducing higher hardware cost [10]. To retain the original cost-effective hybrid beamforming architecture, [29] proposed analog beam broadening strategies to divide the whole array into several virtual subarrays to generate a wider beam and provide constant beamforming gain as much as possible over the whole operating frequency band, thus mitigating the impact of far-field beam-squint to some extent. Furthermore, [30] studied a spatial frequency modulated continuous wave (FMCW)-based phase shift design for massive phased arrays to overcome the near-field beam misfocus problem, and achieve a uniform gain over a wide bandwidth and a higher rate than the standard design in the near-field propagation region.

In 6G mobile communication systems, the use of higher frequency bands, wider bandwidth, and larger number of antennas will pave the way for advanced sensing capabilities with high accuracy and resolution. As a result, future 6G ISAC systems are likely to support a variety of sensing applications, including ultra-high accuracy localization and tracking, simultaneous imaging, mapping, and localization, as well as enhanced human sensing, gesture, and activity recognition[17]. Localization is a typical application of ISAC. Therefore, the concept of integrated localization and sensing communication (ILSC) is proposed to further specify and refine the broader concept of ISAC. In the context of sensing and localization, research efforts towards SLAM and ISAC techniques are well underway. Standardized cellular network localization methods allocate dedicated reference signals, namely the downlink positioning reference signal (PRS) and uplink sounding reference signal (SRS), to facilitate positioning [31]. These signals are instrumental in deriving location information by measuring parameters such as reference signal received power (RSRP), time difference of arrival (TDoA), angle of departure (AoD), angle of arrival (AoA), and multi-cell round trip time (Multi-RTT). Moreover, these measurements further enable precise UT localization, achieving accuracy ranging from a few meters to a few decimeters depending on the specific deployment scenarios [32]. In [33, 34], the authors proposed a joint localization and location-aware beamforming scheme for a semi-passive RIS-assisted system for the single-user and multiple-user scenarios, where the angular parameters were utilized to position the UTs. However, only the LoS path was considered and the delay information was also dismissed. In [35], the author design a joint users’ activity and channel estimation scheme under near-field scenario and sense the users’ location with the TDoA and AoA of the LoS links between multiple subarrays and users. For the near-field propagation region of UM-MIMO systems, the accurate spherical wavefront presents a new way for sensing and localization. Particularly, [36, 37] pioneered to develope a modified two-dimensional version of the multiple signal classification (MUSIC) algorithm for localizing the signal sources with the spherical wavefront sampled by a passive sensor array. Then, [38] derived the theoretical error bound for target tracking by only harnessing the curvature of arrival in phase-difference measurements in the near-field region, and revealed that the range information tends to decrease with increasing target distance. The authors in [39] further proposed a CS-based reconfigurable intelligent surfaces (RIS) aided localization scheme for near-field region depending on the double anchors. Moreover, [12] presented a near-field wideband sensing scheme, which unveiled that coherent processing of spatial domain and frequency measurements contributes to better sensing accuracy compared with [36, 37, 38, 39].

However, the aforementioned schemes necessitated the allocation of dedicated radio resources for positioning signals. To minimize resource overhead and achieve highly integrated functionality, [13] proposed to utilize the beam sweeping procedure in the broadcast synchronization signal block (SSB), where the angular parameters can be extracted for reconstructing images of the propagation environment and providing reliable UT localization service simultaneously. Yet, this method only peroformed well for indoor scenario with smooth concrete wall and required the accurate prior information of the initial positions of UTs and base station (BS), significantly limiting its application scope. Additionally, [9] introduced an integrated uplink channel estimation and localization scheme for the RIS-assisted UM-MIMO system, which synthetically exploited the TDoA and AoA information of hybrid-field spherical wave for accurate UT localization based on the estimated channels, and the localization results were further used to improve the channel estimation performance in turn. However, this method relied on double anchors including the BS and RIS, and also overlooked the design of data transmission for achieving ISAC architecture. Meanwhile, the authors in [14] optimized the ISAC signal to maximize the near-field sensing performance subject to a minimum communication rate while requires the priori target location information and only considered the ideal line-of-sight (LoS) propagation condition. In [15], the UAV-enabled ISAC systems are shown to effectively address the requirements of periodic sensing scenarios. In addition, the authors in [15, 16] have investigated the problem of resource allocation between communication and sensing in ISAC. Additionally, it may be beneficial to explore network-level cooperative ISAC, which leverages multi-cell cooperation to enhance both sensing and communication coverage. As discussed in [40], network-level ISAC can provide additional degrees of freedom for achieving better integration between sensing and communication, thereby improving overall system performance. By coordinating across multiple ISAC base stations, such a network can optimize resource allocation and interference management, offering promising opportunities for overcoming limitations in traditional ISAC setups.

Based on the review above and Table I, there still exists a gap among previous works in investigating ILSC performance with only a single UM-MIMO anchor and considering both HF and beam squint effects simultaneously.

In this work, we focus on achieving ILSC function for time-division duplexing (TDD) mmWave UM-MIMO systems. To be specific, the main contributions of this paper can be summarized as follows:

• •

  We conceive a novel frame structure and signal processing flow designed for ILSC. The proposed architecture consists of a training stage and a data transmission stage. During the training stage, UTs transmit the uplink reference signals to acquire CSI. Moreover, without the need for additional radio resources, a BS can simultaneously sense the positions of scatterers based on the estimated channels and further provide accurate UT localization service by exploiting the hybrid-field propagation characteristic. In the following downlink data transmission stage, the BS further employs location information to design beamforming schemes for more efficient payload data transmission.

• •

  We utilize a projection-based Bayesian CS channel estimation scheme to combat the HFBS effect. By utilizing an elaborate frequency-dependent hybrid-field projection matrix, we ensure that the channel vectors (from one UT’s antenna to all BS’s antennas) represented in the projected domain have the identical sparsity patterns across different sub-carriers and different antennas on the UT side. On this basis, we further consider the generalized multiple measurement vector approximate message passing (GMMV-AMP) algorithm to solve this CS problem, where the expectation-maximization (EM) technique is combined to adaptively learn the unknown hyper-parameters across all subcarriers.

• •

  We propose to simultaneously locate the scatterers and UT based on the acquired CSI, by leveraging the hybrid-field propagation characteristic of UM-MIMO systems. Exploiting acquired CSI, we simultaneously determine locations of scatterers and a UT without requiring additional radio resources. Specifically, we extract the targets’ angle and distance information from the spherical electromagnetic wavefront based on the estimated channels to sense the scatterers’ locations. Then, by treating these sensed scatterers as virtual anchors (VAs), we can coarsely estimate the UT’s location according to the propagation geometry. To further enhance the sensing performance, we utilize super-resolution relative delay estimation so that the UT’s and scatterers’ location can be refined through an iterative procedure.

• •

  We introduce a hybrid beamforming design for downlink data transmission that leverages the location sensing results and overcomes the HFBS effect. Specifically, we first develop an analog beamforming codebook. Within this codebook, each beamforming vector is associated with a pair of estimated distance and angle parameters, and thus can focus the signal energy on the sensed position under the hybrid-field effect. To further mitigate the beam-squint effect, each beam is designed to be broadened along both the distance and angle domains simultaneously to cover the beam-squint region, thus providing constant beamforming gain as much as possible over the whole bandwidth and enhance spectral efficiency (SE). Furthermore, we employ the simultaneous orthogonal matching pursuit (SOMP) algorithm to select optimal analog beamforming vectors and design the corresponding digital precoder to enable multi-stream transmission.

Throughout this paper, scalar variables are represented by normal-face letters, while column vectors and matrices are indicated by boldface lower and upper-case letters, respectively. The transpose, Hermitian transpose, and inversion of matrices are denoted by $(\cdot)^{\rm T}$, $(\cdot)^{\rm H}$, and $(\cdot)^{-1}$, respectively. In addition, $|\cdot|$ represents the modulus, $\|\cdot\|_{p}$ denotes the $\ell_{p}$-norm, and $\|\cdot\|_{\rm F}$ signifies the Frobenius-norm. $[\cdot]_{i}$ denotes the $i$-th element of a vector, and $[\cdot]_{i,j}$ represents the $i$-th row and $j$-th column element of a matrix. Furthermore, $|\cdot|_{c}$ indicates the cardinality of a set. $\mathcal{CN}(\cdot,\mu,\Sigma)$ represents the Gaussian distribution with mean $\mu$ and variance $\Sigma$, while $\mathcal{U}[a,b]$ stands for a uniform distribution between $a$ and $b$. $\mathbb{E}(\cdot)$ and $\mathbb{V}(\cdot)$ calculate the mean and variance of a variable, respectively. $\partial(\cdot)$ denotes the first-order partial derivative operation, and ${\rm diag}(\mathbf{a})$ constructs a diagonal matrix with elements from vector $\mathbf{a}$ placed along its diagonal. Finally, $\mathbf{I}_{n}$ denotes the identity matrix of size $n\times n$, and $\mathbf{0}_{n\times m}$ represents the $n\times m$ zero matrix.

At this stage, multiple UTs transmit mutually orthogonal reference signals in the $Q$ consecutive time slots in the uplink to acquire CSI. For a specific UT, in order to estimate the $m$-th subcarrier’s CSI, the received baseband signal $\tilde{\mathbf{y}}_{\rm UL}[m,q]\in\mathbb{C}^{N_{s}^{p}}$ at the BS in the $q$-th time slot can be formulated as

where $1\leq q\leq Q$ and $1\leq m\leq M$. In (II-A), $\mathbf{W}_{\rm UL}[m,q]$ represents the BS’s receive combining matrix $\mathbf{W}_{\rm UL}[m,q]=\mathbf{W}^{\rm RF}_{\rm UL}[q]\mathbf{W}^{\rm BB}_{\rm UL}[m,q]\in\mathbb{C}^{N_{\rm BS}\times N_{s}^{p}}$, in which $\mathbf{W}^{\rm RF}_{\rm UL}[m,q]\in\mathbb{C}^{N_{\rm BS}\times N_{\rm BS}^{\rm RF}}$ and $\mathbf{W}^{\rm BB}_{\rm UL}[m,q]\in\mathbb{C}^{N_{\rm BS}^{\rm RF}\times N_{s}^{p}}$ are the analog and digital receive combining matrices respectively, $N_{s}^{p}$ is the number of independent received signal streams on each subcarrier. By activating only one RF chain in one time slot for UT, the RF reference signal vector transmitted by the UT can be denoted as $\mathbf{f}_{\rm UL}^{\rm RF}[q]\in\mathbb{C}^{N_{\rm UT}\times 1}$. $P_{t}^{\rm UL}$ is the total uplink transmit power, while $s[m,q]$ is the baseband reference signal with $s[m,q]s^{*}[m,q]=1$. In addition, $\mathbf{H}_{\rm UL}[m]\in\mathbb{C}^{N_{\rm BS}\times N_{\rm UT}}$ stands for the corresponding uplink frequency-domain channel matrix, whose modeling approach will be detailed in Sec. II-B, and $\tilde{\mathbf{n}}[m,q]\in\mathbb{C}^{N_{\rm BS}\times 1}$ is the complex additive white Gaussian noise (AWGN) with the covariance matrix $\sigma^{2}\mathbf{I}_{N_{\rm BS}}$, i.e., ${\mathbf{n}}[m,q]\sim\mathcal{CN}(\mathbf{0}_{N_{\rm BS}},\sigma^{2}\mathbf{I}_{N_{\rm BS}})$. Due to the constant modulus constraint of the adopted fully-connected RF phase shift network at the BS and UT, the uplink transmit and combining matrices can be expressed as $\left[\mathbf{f}_{\rm UL}^{\rm RF}[q]\right]_{i_{1}}=\frac{1}{\sqrt{N_{\rm UT}}}e^{j\vartheta^{1}_{i_{1}}}$ and $\left[\mathbf{W}_{\rm UL}^{\rm RF}[q]\right]_{i_{1},i_{2}}=\frac{1}{\sqrt{N_{\rm BS}}}e^{j\vartheta^{2}_{i_{1},i_{2}}}$ respectively, where $\vartheta^{1}_{i_{1}}$ and $\vartheta^{2}_{i_{1},i_{2}}$ denote the phase value connecting the $i_{1}$-th antenna and the $i_{2}$-th RF chain, respectively.

Considering the channel reciprocity in TDD systems, we focus on the formulation of the uplink channel matrix next. In an mmWave UM-MIMO system with a hyper-scale antenna aperture and ultra-wide bandwidth, the channel exhibits an intrinsic HFBS effect. Our study specifically targets scenarios that are representative of indoor and short-range localization and communication environments, rather than urban macro-cell (UMa) environments. In these indoor settings, both near-field and far-field effects are present, along with significant beam squint effects due to the large antenna aperture and wide bandwidth, which are distinct from the purely far-field propagation considered in UMa scenarios. Under the spherical wavefront assumption, the spatial-frequency domain hybrid near- and far-field channel matrix $\mathbf{H}_{\rm UL}[m]$ can be accurately modeled as (2) shown at the top of next page, which is a sum of the contributions of $(L+1)$ multipath components (MPCs) including one LoS MPC and several non-LoS (NLoS) MPCs reflected by $L$ scattering clusters [42]. Particularly, for the LoS component $\mathbf{H}_{\rm UL}^{\rm LoS}[m]$, its $(n_{1},n_{2})$-th element can be further written as

where $\beta_{\rm LoS}$ is the large-scale fading gain of the LoS link, $\lambda_{m}$ represents the wavelength of the $m$-th subcarrier, and $r_{n_{1},n_{2}}$ denotes the physical distance between the $n_{1}$-th antenna element of the BS and the $n_{2}$-th antenna element of the UT.

For the $l$-th NLoS component, it can be decomposed into the product of frequency-dependent array manifold of BS and UT sides under the spherical wave assumption, which can be further concretely presented as [25]

Herein, $r_{l}^{\rm BS}$ and $\theta_{l}^{\rm BS}$ denote the distance and AoA between the $l$-th scatterer and the reference point of BS’s antenna array, respectively. Without loss of generality, the reference point could be set to the center of the antenna array. More generally, $r_{l,n_{1}}^{\rm BS}$ represents the distance between the $l$-th scatterer and the $n_{1}$-th element of BS’s antenna array. Under the spherical wave assumption, $r_{l,n_{1}}^{\rm BS}$ is a function of both $r_{l}^{\rm BS}$ and $\theta_{l}^{\rm BS}$, i.e., $r_{l,n_{1}}^{\rm BS}=\sqrt{\left({r_{l}^{\rm BS}}\right)^{2}+\delta^{2}_{n_{1}}d^{2}-2r_{l}^{\rm BS}\delta_{n_{1}}d\sin(\theta_{l}^{\rm BS})}$, in which $\delta_{n_{1}}=(2n_{1}-N_{\rm BS}-1)/2$ with $n_{1}=1,2,\dots,N_{\rm BS}$, and $d$ denotes adjacent antenna spacing commonly with half-wavelength. Similarly, $\mathbf{a}_{\rm UT}(\theta_{l}^{\rm UT},r_{l}^{\rm UT},m)$ can be expressed as

where $r_{l}^{\rm UT}$ and $\theta_{l}^{\rm UT}$ denote the distance and AoD between the reference point of UT’s antenna array and the $l$-th scatterer, respectively. Additionally, in (2), $\beta_{\rm NLoS}^{l}$ and $\alpha_{l}\sim\mathcal{CN}(0,\sigma_{\alpha}^{2})$ are the large-scale fading gain of the NLoS link and small-scale channel gain corresponding to the $l$-th NLoS component, respectively.

Note that in practical environments, multi-hop paths experience significant path loss, as shown in [43], where each hop results in approximately 6 dB reflection loss and 12 dB transmission loss in a building scatterer scenario. Consequently, these multi-hop signals are treated as noise, having a minimal effect on both communication and sensing capabilities. Therefore, in the modeling presented in this paper, only single-hop NLoS components are considered.

The designed frame structure for ILSC is depicted in Fig. 2. Its time-frequency radio resources can be divided into multiple resource elements to convey both reference signals and payload data. Specifically, the first $Q=N_{\rm UT}Q_{\rm BS}$ symbols are allocated for channel training, while the remaining resources are reserved solely for payload data transmission. During the full frame, the discrete Fourier transform (DFT) length of an OFDM symbol is set to $M_{p}=M_{\text{CP}}$. Consequently, the symbol duration of the reference signal becomes $(M_{p}+M_{\text{CP}})/B_{w}$, where $M_{\text{CP}}$ represents the CP length, and $B_{w}$ is the system bandwidth. Furthermore, the reference signals in the training stage consists of $Q_{\text{BS}}$ subframes. Within each subframe, the BS analog receive combiner $\mathbf{W}_{\text{UL}}^{\text{RF}}[q]$ remains constant, and its phase shift of $(i_{1},i_{2})$-th elements $\vartheta^{2}_{i_{1},i_{2}}$ for $1\leq i_{1}\leq N_{\rm BS}$ and $1\leq i_{2}\leq N_{\rm BS}^{\rm RF}$ are uniformly distributed in $\mathcal{U}[0,2\pi)$. Meanwhile, the digital receive combiner at the BS adopts an identity matrix $\mathbf{W}^{\text{BB}}_{\text{UL}}[m,q]=\mathbf{I}_{N_{\text{BS}}^{\text{RF}}}$, where $N_{s}^{p}=N_{\text{BS}}^{\text{RF}}$. In this setting, UTs transmit $N_{\text{UT}}$ orthogonal reference signals in each subframe, which could be designed as a DFT matrix as $\mathbf{F}_{\text{UL}}=\left[\mathbf{f}_{\text{UL}}^{\text{RF}}[1],\mathbf{f}_{\text{UL}}^{\text{RF}}[2],\dots,\mathbf{f}_{\text{UL}}^{\text{RF}}[N_{\text{UT}]}\right]\in\mathbb{C}^{N_{\text{UT}}\times N_{\text{UT}}}$ to meet the prerequisites of constant modulus and orthogonality.

Based on the above assumptions, by collecting $\tilde{\mathbf{y}}_{\rm UL}[m,(q_{1}-1)N_{\rm UT}+q_{2}]$ for all $Q_{\rm BS}$ subframes together (where $q_{1}$ is the subframe index and $q_{2}$ is the symbol index in one subframe), the aggregate received signals $\tilde{\mathbf{Y}}_{\rm UL}[m]\in\mathbb{C}^{P\times N_{\rm UT}}$ can be written as

where $\tilde{\mathbf{Y}}_{\rm UL}[m]$ successively selects every symbol of the same subframes $1\leq q_{2}\leq N_{\rm UT}$ as its columns, and stacks different subframes $1\leq q_{1}\leq Q_{\rm BS}$ on top of each another, and $P=Q_{\rm BS}N_{\rm BS}^{\rm RF}$ denotes the effective measurement length. Besides, $\mathbf{W}_{\rm UL}=\left[\mathbf{W}_{\rm UL}^{\rm RF}[1],\mathbf{W}_{\rm UL}^{\rm RF}[2],\dots,\mathbf{W}_{\rm UL}^{\rm RF}[Q_{\rm BS}]\right]\in\mathbb{C}^{N_{\rm BS}\times P}$, $\mathbf{S}[m]=\sqrt{P_{t}^{\rm UL}/M}{\rm diag}\left\{s[m,1],s[m,2],\dots,s[m,N_{\rm UT}]\right\}$ is the stacked baseband reference signal, and $\tilde{\mathbf{N}}[m]\in\mathbb{C}^{N_{\rm BS}\times N_{\rm UT}}$ stands for the corresponding AWGN for all these $Q_{\rm BS}$ subframes. Finally, the received baseband signal in (6) is post-processed by multiplying $\mathbf{S}^{\rm H}[m]\mathbf{F}_{\rm UL}^{\rm H}$ right-hand, i.e.,

in which ${\mathbf{N}}[m]=\mathbf{W}^{\rm H}_{\rm UL}\tilde{\mathbf{N}}[m]\mathbf{S}^{\rm H}[m]\mathbf{F}_{\rm UL}^{\rm H}$ is the effective noise term.

In the case of hybrid far- and near-field propagation effect, the planar wave approximation based projection matrix fails to exactly capture the inherent characteristics of the channel, leading to a degradation in the sparsity of the channels represented in such a projection matrix, as illustrated in Fig. 3(a). To address this issue, a more accurate Fresnel approximation [6] is developed for the distance $r_{l,n_{1}}^{\rm BS}$ based on a second-order Taylor expansion as

Based on this approximation, the BS array manifold of each MPC depends on not only its AoD/AoA but also on the distance between the BS and scatterers. To this end, [25, 9] proposed a non-uniform two-dimensional sampling lattice that considers both angles and distances as

where $\rho N_{\rm BS}$ and $S$ are the number of grids along the angles and distances, respectively, and $R_{r}=2D_{\rm BS}^{2}/\lambda_{c}$ is the Rayleigh distance corresponding to the antenna array aperture $D_{\rm BS}$. Furthermore, $\eta$ is a scaling factor utilized to control the sampling density.

In this way, the projection matrix $\boldsymbol{\Phi}[m]\in\mathbb{C}^{N_{\rm BS}\times\rho N_{\rm BS}S}$ can be designed as

in which the frequency-selective array manifold $\mathbf{a}_{\rm BS}(\theta,r,m)$ sampled from the defined 2D lattice forms the $f(\theta,r)$-th candidate columns of the projection matrix $\boldsymbol{\Phi}[m]$. It is essential to note that when $s=0$, $r_{n,s}\rightarrow\infty$, the defined sampling grids match the traditional virtual-angular lattices in the far-field. This alignment indicates that the projection matrix $\boldsymbol{\Phi}[m]$ serves as a generalized representation of channels in both far-field and near-field. Furthermore, in contrast to the frequency-flat projection matrix leading to the offset of sparse support sets across different subcarriers [25], the proposed matrix considers the frequency-dependent behavior of the array manifold to ensure an identical support set across the entire frequency spectrum, as illustrated in Fig. 3(b).

As a result, the NLoS components of the HFBS channel can be readily represented in a sparse manner using the defined projection matrix. Furthermore, while the LoS component in (3) cannot be straightforwardly decomposed into the product of two array manifold vectors, it can be approximated as follows using the piece-wise subarray planar wave for UT’s array as depicted in Fig. 3(c)

where $G$ indicates the number of subarrays divided for UT’s array, ${\mathbf{a}_{\rm UT}^{(g)}}(m)=\mathbf{a}_{\rm UT}(\theta_{g}^{\rm UT},\infty,m)$ and ${\mathbf{a}_{\rm BS}^{(g)}}(m)=(\theta_{g}^{\rm BS},r_{g}^{\rm BS},m)$ denote the UT’s array steering vectors of the $g$-th subarray and the BS’s array steering vectors, respecitvely. Under this approximation, the LoS component can be represented sparsely with $G$ components as

Therefore, the channel estimation in (7) is equivalent to a structured sparse signal recovery problem

where $\mathbf{A}[m]=\mathbf{W}^{\rm H}_{\rm UL}\boldsymbol{\Phi}[m]\in\mathbb{C}^{P\times\rho N_{\rm BS}S}$ is the effective sensing matrix and $\mathbf{H}_{\rm UL}^{P}[m]\in\mathbb{C}^{\rho N_{\rm BS}S\times N_{\rm UT}}$ denotes the high-dimensional projected sparse channel matrix.

It is well-known that AMP is an efficient iterative Bayesian algorithm for CS problems, which employs low-complexity message passing for approximating the traditional minimum mean square error (MMSE) estimation [44]. For notation brevity, the superscript, subscript, and index $m$ have been omitted for ${\mathbf{Y}}_{\rm UL}[m]$ and $\mathbf{H}_{\rm UL}^{P}[m]$ hereinafter, and the different entries of matrix is distinguished by its subscript, e.g., $h_{k,n}=\left[\mathbf{H}\right]_{k,n}$. Starting with the MMSE estimate of $\mathbf{H}$, we can obtain its posterior mean as

where $p({h}_{k,n}|\mathbf{Y})$ denotes the marginal posterior distribution expressed as follows:

(15) requires the high-dimensional integral of joint posterior distribution $p(\mathbf{H}|\mathbf{Y})$ w.r.t. each variable entry of $\mathbf{H}$ except ${h}_{k,n}$. Based on the Bayesian theorem, $p(\mathbf{H}|\mathbf{Y})$ can be represented as

Under the assumption of AWGN, the above likelihood function $p(\mathbf{Y}|\mathbf{H})$ obeys complex Gaussian distribution

By following the classic AMP literature [45], we consider the prior distribution of $\mathbf{H}$ to follow an i.i.d. Bernoulli-Gaussian distribution for effectively capturing its sparsity, i.e.,

where $0\leq\lambda_{k,n}\leq 1$ denotes the non-zero probability of $h_{k,n}$, and $\mu_{k,n}$, $\gamma_{k,n}$ are the prior mean and variance of a Gaussian distribution, respectively.

The joint distribution $p(\mathbf{H},\mathbf{Y})$ can be represented by a bipartite graph $\mathcal{G}=(\mathcal{V},\mathcal{F},\mathcal{E})$, which consists of variable nodes $\mathcal{V}$, factor nodes $\mathcal{F}$, and the corresponding edges $\mathcal{E}$ [46]. This suggests the use of message passing algorithms contributes to rigorous solution of the marginal posterior distribution in (15), and the messages approximation can further reduce the complexity of high-dimensional integral. Therefore, it yields [19]

where $t$ indicates the iteration index, $\Sigma_{k,n}^{t}$ and $R_{k,n}^{t}$ are the message parameters updated at variable nodes $\mathcal{V}$ as

while $V_{p,n}^{t}$ and $Z_{p,n}^{t}$ are the message parameters updated at factor nodes $\mathcal{F}$ as

where $\hat{h}_{k,n}^{t}$ and $\hat{v}_{k,n}^{t}$ are defined in (24). By combining (18) and (19), the marginal posterior distribution can be merged as

where

Then, the posterior mean and variance can be explicitly calculated as

The message switching procedures mentioned in (20)-(24) constitute the fundamental steps of the AMP algorithm. However, several concerns arise when applying the AMP algorithm to practical systems: (a) The hyper-parameters in the prior assumption and the noise variance are often unknown in practice. (b) The column-wise independent solution to the channel estimation problem neglects the common support structure that arises from spatial and frequency correlations. (c) The AMP algorithm is primarily designed for large system limits and typically assumes the use of i.i.d. Gaussian sensing matrices. This assumption is unrealistic for the problem described in (13).

To overcome the first challenge, EM can be exploited to learn the unknown hyper-parameters $\boldsymbol{\theta}=\{\lambda_{k,n},\mu_{k,n},\gamma_{k,n},\sigma^{2},\forall k,n\}$ as [47]

where $\mathbb{E}\left\{\cdot|\mathbf{Y};\boldsymbol{\theta}^{t}\right\}$ denotes the expectation with respect to the posterior distribution $p(\mathbf{H}|\mathbf{Y};\boldsymbol{\theta}^{t})$, and its approximation is concluded in (22). By updating one of the elements in (25b) at a time and holding others constant, the learning rules of the hyper-parameters can be summarized as [19]

To further address the second concern aforementioned, it is reasonable to assign a common non-zero probability $\lambda_{k,n}$ for each different columns and subcarriers of $\mathbf{H}$, which indicates that the learning rule in (26a) and can be refined as

Moreover, the damping operation [48] could be introduced to prevent AMP from divergence due to the nonideal sensing matrix $\mathbf{A}[m]$ and the medium size problem mentioned in the concern (c). To sum up, the overall procedure of the AMP-EM-based HFBS channel estimation algorithm is presented in Alg. 1.

Under the spherical wave assumption, the HFBS channel depends not only on the angle of the target, but also on its distance. Therefore, it is feasible to simultaneously derive the angle and distance features of the target from the channel state information. According to the analysis in Section III-B, under ideal conditions, the non-zero support set of sparse channel $\mathbf{H}_{\rm UL}^{P}[m]$ exactly reflects the approximate LoS components and the NLoS MPCs of the HFBS channel, and their angle and distance parameters are associated with the spatial spatial sampling grid defined in (9). In this way, we can map the location of the center of UT’s subarrays and the scatterers base on the channel estimation results.

One of the premise to achieve location mapping lies on accurate knowleadge about MPCs’ number. To this end, we estimate the effective MPCs’ number $L^{\prime}=G+L$ based on the minimum description length (MDL) principle, which selects the parametric probability model that fits the dataset best based on the maximum likelihood criterion, so as to estimate the number of signal sources without the need of empirical threshold [50].

Firstly, in order to overcome the rank-deficient problem of covariance matrix caused by coherent multipath signals, the covariance matrix of received pilot signals is calculated with the spatial smoothing technique [49] as

where $\mathbf{R}_{\rm yy}[m]\approx{\mathbf{Y}}_{\rm UL}^{\rm H}[m]{\mathbf{Y}}_{\rm UL}[m]$ is the estimated covariance matrix of the original signal, while $\mathbf{I}_{N_{\rm UT}-K,i}=\left[\mathbf{0}_{(N_{\rm UT}-K)\times i}\;\mathbf{I}_{N_{\rm UT}-K}\;\mathbf{0}_{(N_{\rm UT}-K)\times(K-i)}\right]\in\mathbb{C}^{(N_{\rm UT}-K)\times N_{\rm UT}}$. By the spectrum decomposition theorem, if there are $L^{\star}$ of the multipath, $\mathbf{R}_{\rm yy}^{K}$ can be decomposed into

where $\lambda_{i}$ and $\mathbf{u}_{i}$ are the eigenvalues and eigenvectors of $\mathbf{R}_{\rm yy}^{K}$, respectively. Defines a set of parameters for the received signal model$\boldsymbol{\eta}=\left\{\lambda_{1},\dots,\lambda_{L^{\star}},\mathbf{u}_{1},\dots,\mathbf{u}_{L^{\star}},\sigma^{2}_{n}\right\}$, the estimate of $L^{\star}$ can be obtained based on the MDL rule as

where $P({\mathbf{Y}}_{\rm UL}|{\boldsymbol{\eta}})$ is the joint probability density distribution of the received pilot signals, while $\gamma(L^{\star})=\frac{1}{2}L^{\star}(2N_{\rm UT}-L^{\star}-K)\log N_{\rm BS}$. Under the assumption of complex Gaussian distribution, it can be further expressed as (31), which is shown at the top of the next page, where $\hat{\lambda}_{i}$ is the estimate of eigenvalue by decomposing the estimated covariance matrix $\mathbf{R}_{\rm yy}^{K}$. By calculating the covariance matrix of the received pilot signal on different subcarriers, $M$ independent estimation results can be obtained in parallel. Finally, the fusion of different estimation results can be completed by the mode operation.

When the corresponding row energy of the estimated sparse channel is greater than the noise threshold, the MPCs can be considered to exist. Extract the angle and distance parameters of them

where $\hat{h}_{\rm UL}^{P}(\theta,r)=\sum_{n=1}^{N_{\rm UT}}\sum_{m=1}^{M}\left|\left\{\widehat{\mathbf{H}}_{\rm UL}^{P}[m]\right\}_{f(\theta,r),n}\right|^{2}$. Unfortunately, due to the finite-resolution sampling lattice, the energy leakage is an inevitable negative effect. As a result, the grids set $(\tilde{\varTheta},\tilde{\mathcal{R}})$ exhibits a cluster distribution, where $\hat{L}^{\prime}$ truthful peaks are surrouded by a few ghost points. In order to divide clusters before peak identification, K-Means clustering algorithm [52] can be exploited to select the spatial grids subset $(\tilde{\varTheta}_{l},\tilde{\mathcal{R}}_{l})$ for the $l\,(1\leq l\leq\hat{L}^{\prime})$-th effective MPC. In this case, the angle and distance of $l$-th MPC corresponding to the summit of the $l$-th cluster grids can be represented as

Considering the location sensing geometry illustrated in Fig. 4, where the BS is located at the origin of coordinates, we can coarsely map the position of VAs with the orientation and distances measured above as $(\hat{x}_{l},\hat{y}_{l})=(\hat{r}_{l}^{\rm BS}\cos\hat{\theta}_{l}^{\rm BS},\hat{r}_{l}^{\rm BS}\sin\hat{\theta}_{l}^{\rm BS})$, where we treat the sensed UT subarray centers and scatterers equally as VAs, and they can contribute to UT location sensing hereinafter. In order to distinguish the sensed UT subarray centers from the scatterers, we adopt the following method: VAs satisfying the array aperture constraint are inferred as UT subarray centers, and their coordinate set $\mathcal{G}=\{l|1\leq l\leq\hat{L^{\prime}},(\hat{x}_{l}-\hat{x}^{\star})^{2}+(\hat{y}_{l}-\hat{y}^{\star})^{2}\leq D_{\rm UT}\}$, where $(\hat{x}^{\star},\hat{y}^{\star})$ is the VAs with the maximum energy and thus can be truthfully inferred as LoS components, while $D_{\rm UT}$ is the antenna array aperture of UT. The remains are classified as scatterers and their coordinate set $\mathcal{S}=\{l|1\leq l\leq\hat{L}^{\prime},l\notin\mathcal{G}\}$.

In order to acquire the UT’s position with the aid of sensed VAs, one of the prerequisites lies in acquiring the orientation between the UT and VAs. For this purpose, motivated by the channel eigenvector decomposition in (2) and LoS approximation in (III-B), and the BS steering vector can be estimated as $\mathbf{a}_{\rm BS}(\theta_{l}^{\rm BS},r_{l}^{\rm BS},m)\approx\mathbf{a}_{\rm BS}(\hat{\theta}_{l}^{\rm BS},\hat{r}_{l}^{\rm BS},m)$, thus the angle and distance can be estimated via matching pursuit as