Environment Sensing Considering the Occlusion Effect: A Multi-View Approach

The emergence of innovative wireless communication technologies, such as ultra-massive multiple input multiple output (MIMO) technology, intelligent reflective surface (IRS), and wireless artificial intelligence (AI) [1, 2, 3], etc., provides more possibilities for the development of future wireless communication systems design. In the foreseeable future wireless communication application scenarios, autonomous driving [4], intelligent robot localization [5], and unmanned aerial vehicle (UAV) control [6] require not only wireless broadband connections, but also accurate environmental information, including the location, shape, and electromagnetic (EM) characteristics, etc. of objects in the environment. Therefore, integrated sensing and communication (ISAC), as a research hotspot for the next generation wireless communication system, aims to use the wireless communication equipment and infrastructure to achieve environment sensing.

Such kind of environment sensing has long been accomplished by traditional radar technology and its later evolution of joint radar and communication [7]. The latter usually aims to achieve the sharing of software and hardware equipment and time-frequency resources between radar equipment and communication equipment. Compared with traditional radar or early joint radar and communication systems, efficient use of communication signals in the existing communication systems for environment sensing, and efficient exploitation of the environment sensing results to enhance communication, are two major design goals of the ISAC systems. For the former goal, one straightforward realization is to fully exploit the environmental information embedded in the received communication data, due to the simple fact that the distribution of environmental scatterers dramatically affects the wireless multipath channels. For the latter, direct channel prediction (reconstruction) can be conducted based on the sensed distribution and characteristics of scatterers within the environment thus to enhance the performance of communication [8, 9, 10]. In this paper, we mainly focus on the first part and study how to perform efficient and accurate environment sensing based on the uplink data received from multiple users, in a way compatible with existing communication systems thus to achieve a smooth integration of communication and sensing.

One major challenge to do this is the large number of unknown variables embedded in the generally complicated environment, including the unknown location and EM characteristics of the environmental scatterers themselves, and the very common occlusion effect within and between scatterers, which blocks the propagation of EM waves and makes portions of some scatterers or even all of them invisible to the receivers. In such a case, it is hard for a single user or base station (BS) to fully sense the entire environment and the different views of the environment from different users and/or BS(s) should be exploited. However, a challenge still remains since in each separate sensing (transmitting and receiving), we generally cannot know which scatterers or which part of them participates in the final generation of the received signal in an environment yet to know, and it differs from user to user and BS to BS. This motivates us to develop a practical and efficient multi-view sensing algorithm to deal with the above problem.

So far, the design of ISAC system has raised extensive research effort. In [11], the author comprehensively discussed the novel applications, key performance requirements, challenges, and future research directions of ISAC system designs. In [12], the author compared the difference between the radar sensing scheme and the ISAC scheme. The author showcased that a common signal wave should be jointly designed to realize the integration of sensing and communication, and the advantages and disadvantages of various design methodologies are also discussed.

Rapid progresses have been witnessed in the practical ISAC system and algorithm design recently. [13] designed an ISAC system based on the wireless cellular networks. The author jointly designed signal waveforms suitable for sensing and communication based on 5G NR (New Radio) and its frame structure. Also under the 5G framework, [14] proposed a multi-beam based sensing and communication scheme. In addition to the beams used for communication, some other ones are used to sense the environment around the BS. Based on the 4G and 5G frameworks, [15] studied the processing principles, implementation challenges, and performance of orthogonal frequency division multiplexing (OFDM)-based radars, and achieved excellent sensing performance by using high channel bandwidth and configurable sub-carriers. In a Wi-Fi communication context, based on the Doppler frequency shift embedded in the channel state information (CSI), [16] achieved the sensing of human body posture. In the regime of internet of vehicles, [17] designed a fast beam alignment and tracking algorithm based on a hardware test platform to achieve vehicle tracking while performing millimeter wave communication with low latency and high data rate. As the research interest on ISAC keeps on increasing, more and more potential technologies are also introduced into this field, such as channel modeling [18], joint beam optimization [19] and machine learning [20], etc.

As for the environment sensing system design, some works have made full use of the sparsity of environmental information and achieved better and lower-complexity solutions based on the compressed sensing (CS) [21, 22] framework. In the field of microwave computational imaging, the utilization of the intrinsic sparsity of scatterers within the environment is key to their effective detection. This type of algorithm usually treats the scatterer as a sparsely distributed point cloud, thereby convert the imaging problem into a CS reconstruction problem, which is solved by orthogonal matching pursuit (OMP) [23], approximate message passing (AMP) [24], generalized approximate message passing (GAMP) [25] and other widely used methods.

Recently, in [26] and [27], the authors innovatively used BS to sense actively, and based on the block sparse Bayesian learning (BSBL) algorithm, proposed an IRS-assisted microwave computing imaging method. Although the above-mentioned methods effectively exploit the sparsity of the environment but ignore the occlusion effect caused by the scatterers. They ideally believe that EM signals can be freely transmitted to any location in the environment, and there still lacks sufficient consideration of the occlusion effect in the ISAC system design.

In this paper, based on the existing wireless communication framework, we design a multi-view environment sensing scheme taking into full consideration the occlusion effect in an outdoor scenario.

Different from the above-mentioned application scenarios, we have achieved ISAC by noninvasively making use of the multi-user uplink communication data. At the same time, we have designed an effective sensing algorithm called GAMP-MVSVR which refers to the generalized-approximate-message-passing-based multi-view sparse vector reconstruction. In the uplink communication scenario, we consider an environment sensing scenario with single or multiple BSs, where users and BSs perform multi-view and collaborative sensing, and a centralized processor retrieving the environmental information from all the channel information obtained at the BSs. As a special case of the multi-BS scheme, a single BS itself can also complete all the sensing tasks, and we compare the performance differences between the two schemes.

Our design is depicted as follows. In one or more time slots, the BS estimates the multi-user uplink channel after detecting user pilot symbols (sounding signals). We propose an occlusion detection model based on the geometric location of the cloud points and convert the environment sensing problem into a CS reconstruction problem with the occlusion effect. Then we propose a probabilistic reasoning model based on the factor graph by analyzing the relationship among multipath channels, environmental information, and occlusion effects. Based on the well-known AMP method, the proposed GAMP-MVSVR algorithm achieves environment sensing by iteratively estimating the scattering coefficients of the cloud points and their occlusion relationship. In each iteration, the occlusion relationship between the cloud points of the sparse environment is recalculated according to the proposed occlusion detection rule, and in turn, used to estimate the scattering coefficients of the cloud points. Our proposed algorithm achieves improved sensing performance with multi-BS collaboration in addition to the multi-views from the users. To the best of our knowledge, there has been little work on the design of such an ISAC system in the literature.

The main contributions of this paper are summarized as follows:

• •

  We propose an ISAC scheme, which makes full use of the different views of the environment from different users and/or BS(s) to achieve improved sensing performance with multi-BS collaboration in addition to the multi-views from the UEs.

• •

  We propose an effective environment sensing algorithm called GAMP-MVSVR that considers the occlusion effect, where a multi-layer factor graph is constructed to iteratively estimate the scattering coefficients of the cloud points and derive their occlusion relationship based on a simple occlusion detection rule.

• •

  We analyze the impact of the multi-view schemes on environment sensing range and system performance. Extensive simulation results verify the convergence and effectiveness of the proposed algorithm.

The rest of this paper is organized as follows. Section II presents the environment setting and system model in the uplink communication scenario. Section III proposes the environment sensing problem formulation. Section IV proposes the GAMP-MVSVR algorithm under the occlusion effect. In section V, we analyze the impact of the multi-view schemes on environment sensing range and system performance. Finally, section VII presents the numerical results, and section VIII concludes the paper.

Notations: Fonts $a$, $\bm{a}$ and $\mathbf{A}$ represent scalars, vectors and matrices, respectively. Notations $\mathbf{A}^{\rm{T}}$ and $\|\mathbf{A}\|_{F}$ denote transpose and Frobenius norm of $\mathbf{A}$, respectively. $\mathbf{A}_{i,:}$ and $\mathbf{A}_{:,j}$ represent the $i$-th row and $j$-th column of matrix $\mathbf{A}$, respectively. $\odot$ represents the Hadamard product between two matrices. Finally, notation ${\rm diag}(\bm{a})$ represents a diagonal matrix with the entries of $\bm{a}$ on its main diagonal, and $\delta(\cdot)$ is the Dirac delta function.

As shown in Fig. 1, we consider that in the outdoor wireless communication scenario, multiple BSs are deployed and there are multiple active user equipments (UEs). There are some buildings (target objects, serve as scatterers) in the outdoor environment that affect the wireless communication channel. We assume w.l.o.g. that the channel is a quasi-static channel where one block of transmission time $T_{\rm b}$ is much shorter than the channel coherence time $T_{\rm c}$, $T_{\rm b}<<T_{\rm c}$, and the UE is stationary within the coherence time of the multipath channel estimation by UE pilots. In the uplink communication scenario, multiple single-antenna UEs send uplink communication data to multi-antenna BSs. The transmitted signal is scattered by the scatterer, and received by BSs through multipath channel. Therefore, the received signal of the BS contains environmental information.

We consider two schemes: a single-BS scheme and a multi-BS joint scheme to achieve multi-view environment sensing. The specific analysis is as follows.

• •

  As shown in Fig. 1, the single-BS environment sensing scheme is treated as a special case of the multi-BS joint scheme. When multiple UEs send uplink data to BS 1, the uplink communication channel is affected by the target scatterer. In the scenario, due to the occlusion effect, the multipath channel of UE 1 is only affected by the scattering of target scatterers 4 and 5, and the multipath channel of UE 2 is only affected by the scattering of target scatterers 5. Therefore, not all scatterers in the environment will affect the multipath channel of the same user, and the occlusion effect makes it hard for a single user to fully sensed the entire environment. The different views of the environment from different users should be exploited to deal with the occlusion effect.

• •

  Fig. 1 shows the scenario where multiple BSs are used for joint environment sensing, where the multi-BS scheme is performed. Sensing by multiple BSs will effectively expand the sensing range. In addition to the same user multi-view method as in the single-BS scheme, the different views of the environment from different BSs should be exploited. For example, in Fig. 1, BS 2 cannot receive the scattered signal from the target scatterer 3, and BS 1 and 3 also have scatterers that cannot be observed. Therefore, a multi-BS joint scheme should be designed to achieve the sensing of the entire environment.

For the receiving antennas deployed on the BS, the single-BS scheme and the multi-BS joint scheme have the same system model. Let the number of users in the environment be $N_{\rm{u}}$, and the total number of receiving antennas of all BSs is $N_{\rm{R}}$. In the uplink communication scenario, the channel from the user to the BS is mainly composed of two parts: the line-of-sight channel ${\bf{H}}^{{\rm{LOS}}}$ from the user to the BS, and the multipath channel ${\bf{H}}^{{\rm{S}}}$ caused by the scatterer.

We discretize environmental information and treat the environmental information in the entire environment as point clouds. Each point in the point cloud represents the environmental information of the small cube with length $l_{\rm{s}}$, width $w_{\rm{s}}$, and height $h_{\rm{s}}$ around it, and these small cubes are called pixels. The inside of each pixel may be empty, or there may be scatterers. The length, width and height of the environment to be sensed are $L_{\rm{s}}$, $W_{\rm{s}}$ and $H_{\rm{s}}$ respectively, then the number of point clouds is ${N_{\rm{s}}}=\left({{{L_{\rm{s}}}}\mathord{\left/{\vphantom{{{L_{\rm{s}}}}{{l_{\rm{s}}}}}}\right.\kern-1.2pt}{{l_{\rm{s}}}}}\right)\times\left({{{W_{\rm{s}}}}\mathord{\left/{\vphantom{{{W_{\rm{s}}}}{{w_{\rm{s}}}}}}\right.\kern-1.2pt}{{w_{\rm{s}}}}}\right)\times\left({{{H_{\rm{s}}}}\mathord{\left/{\vphantom{{{H_{\rm{s}}}}{{h_{\rm{s}}}}}}\right.\kern-1.2pt}{{h_{\rm{s}}}}}\right)$. Fig. 2 shows the method of discretizing environmental information into point clouds. We use the scattering coefficient ${x_{{n_{\rm{s}}}}}$ to represent the scattering coefficient of the small cube where the $n_{\rm{s}}$-th cloud point is located. If the small cube is empty, then ${x_{{n_{\rm{s}}}}}=0$. Otherwise, ${x_{{n_{\rm{s}}}}}>0$. Therefore, the environmental information to be sensed is expressed as ${\bm{x}}={\left[{{x_{1}},{x_{2}},\cdots,{x_{{N_{\rm{s}}}}}}\right]^{\rm{T}}}.$

Multiple users in the space share time-frequency resources. On a single subcarrier, the antennas of all BSs are analyzed independently. The signal received at the $n_{\rm{R}}$-th receiving antenna can be expressed as

$$\small{\bf{Y}}\left({{n_{\rm{R}}}}\right)={\bf{S}}\left({{{\bf{H}}^{\rm{S}}}\left({{n_{\rm{R}}}}\right)+{{\bf{H}}^{{\rm{LOS}}}}\left({{n_{\rm{R}}}}\right)}\right)+{\bf{W}},$$

(1)

where ${\bf{S}}\in\mathbb{C}^{N_{\rm{T}}\times N_{\rm{u}}}$ represents codewords with the length of $N_{\rm{T}}$ sent by $N_{\rm{u}}$ users, ${\bf{Y}}\left({{n_{\rm{R}}}}\right)\in\mathbb{C}^{N_{\rm{T}}\times 1}$ represents the received signal of the $n_{\rm{R}}$-th BS antenna, and $\bf{W}$ represents noise. The term ${{\bf{H}}^{{\rm{LOS}}}}\left({{n_{\rm{R}}}}\right)\in\mathbb{C}^{N_{\rm{u}}\times 1}$ represents the channel coefficient of the line-of-sight path from the user to the receiving antenna, where the free space channel is decomposed into amplitude attenuation and phase shift, $h^{{\rm{LOS}}}_{n_{\rm{u}}}\left({{n_{\rm{R}}}}\right)=\alpha^{\rm{LOS}}\cdot e^{j\varphi^{\rm{LOS}}}$, $\alpha^{\rm{LOS}}$ is the line-of-sight channel amplitude attenuation, and $\varphi^{\rm{LOS}}$ is the line-of-sight channel phase shift. The line-of-sight channel occlusion can be detected through beam scanning, received power, etc. The term ${{\bf{H}}^{{\rm{S}}}}\left({{n_{\rm{R}}}}\right)\in\mathbb{C}^{N_{\rm{u}}\times 1}$ represents the channel coefficient of the multipath channel from the user to the receiving antenna.

As mentioned in Section I-B, we consider that in ISAC system design, the utilization of the intrinsic sparsity of objects or scatterers within an environment is key to the effective detection. We separate the free space channel (microwave background field) and scatterers in the multipath channel model to obtain sparse environmental information. In this way, the occlusion relationship between the scatterer and the antenna can be modeled by the existence of the free space channel between them. The multipath channel from the user to the receiving antenna ${{\bf{H}}^{{\rm{S}}}}\left({{n_{\rm{R}}}}\right)$ is expressed as

		$\displaystyle{{\bf{H}}^{{\rm{S}}}}\left({{n_{\rm{R}}}}\right)={\bf{\tilde{H}}}\left({{n_{\rm{R}}}}\right){\bm{x}}$		(2a)
	$\displaystyle=$	$\displaystyle\left[{{\bf{H}}\left({{n_{\rm{R}}}}\right)\odot{\bf{V}}\left({{n_{\rm{R}}}}\right)}\right]{\bm{x}}$		(2b)
	$\displaystyle=$	$\displaystyle\left({{{\bf{H}}^{{\rm{U\rightarrow s}}}}\odot{{\bf{V}}^{{\rm{U\rightarrow s}}}}}\right){\rm{diag}}\left({{{\bf{H}}^{{\rm{s\rightarrow B}}}}\left({{n_{\rm{R}}}}\right)\odot{{\bf{V}}^{{\rm{s\rightarrow B}}}}\left({{n_{\rm{R}}}}\right)}\right){\bm{x}},$		(2c)

where the free space channel coefficient ${\bf{\tilde{H}}}\left({{n_{\rm{R}}}}\right)$ with occlusion is expressed as the Hadamard product of the free space channel coefficient ${\bf{H}}\left({{n_{\rm{R}}}}\right)\in\mathbb{C}^{N_{\rm{u}}\times N_{\rm{s}}}$ and the occlusion matrix ${\bf{V}}\left({{n_{\rm{R}}}}\right)\in\left\{0,1\right\}^{N_{\rm{u}}\times N_{\rm{s}}}$. The element 0 in the occlusion matrix ${\bf{V}}\left({{n_{\rm{R}}}}\right)$ indicates that the path from $n_{\rm{u}}$-th user to $n_{\rm{R}}$-th receiving antenna is blocked, and the element 1 indicates that the path from $n_{\rm{u}}$-th user to $n_{\rm{R}}$-th receiving antenna is not blocked. ${{\bf{H}}^{{\rm{U\rightarrow s}}}}\in\mathbb{C}^{N_{\rm{u}}\times N_{\rm{s}}}$ represents the channel coefficient from the user to the point cloud position, where $h^{{\rm{U\rightarrow s}}}_{n_{\rm{u}},n_{\rm{s}}}=\alpha^{\rm{U\rightarrow s}}\cdot e^{j\varphi^{\rm{U\rightarrow s}}}$, $\alpha^{\rm{U\rightarrow s}}$ represents the channel amplitude attenuation from the user to the point cloud position, $\varphi^{\rm{U\rightarrow s}}$ represents the channel phase shift from the user to the point cloud position. ${{\bf{H}}^{{\rm{s\rightarrow B}}}}\left({{n_{\rm{R}}}}\right)\in\mathbb{C}^{N_{\rm{s}}\times 1}$ represents the channel coefficient from the point cloud position to the BS receiving antenna, where $h^{{\rm{s\rightarrow B}}}_{n_{\rm{s}}}\left({{n_{\rm{R}}}}\right)=\alpha^{\rm{s\rightarrow B}}\cdot e^{j\varphi^{\rm{s\rightarrow B}}}$, $\alpha^{\rm{s\rightarrow B}}$ represents the channel amplitude attenuation from the point cloud position to the BS receiving antenna, $\varphi^{\rm{s\rightarrow B}}$ represents the channel phase shift from the point cloud position to the BS receiving antenna, ${\bf{H}}\left({{n_{\rm{R}}}}\right)={{\bf{H}}^{{\rm{U\rightarrow s}}}}{\rm{diag}}\left({{{\bf{H}}^{{\rm{s\rightarrow B}}}}\left({{n_{\rm{R}}}}\right)}\right)$. ${{\bf{V}}^{{\rm{U\rightarrow s}}}}\in\left\{0,1\right\}^{N_{\rm{u}}\times N_{\rm{s}}}$ represents the occlusion matrix from the user to the point cloud position, and its distribution reflects the different views from different users. When ${{\bf{V}}^{{\rm{U\rightarrow s}}}}$ contains an all-zero column, it means that all users cannot sense the corresponding pixel, and the corresponding pixel is out of the sensing range. ${{\bf{V}}^{{\rm{s\rightarrow B}}}}\left({{n_{\rm{R}}}}\right)\in\left\{0,1\right\}^{N_{\rm{s}}\times 1}$ represents the occlusion matrix from the point cloud position to the receiving antenna, and its distribution reflects the different views from different BSs. When ${{\bf{V}}^{{\rm{s\rightarrow B}}}}\left({{n_{\rm{R}}}}\right)$ contains a zero element, it means that the receiving antenna cannot sense the corresponding pixel, and the corresponding pixel is also out of the sensing range, ${\bf{V}}\left({{n_{\rm{R}}}}\right)={{\bf{V}}^{{\rm{U\rightarrow s}}}}{\rm{diag}}\left({{{\bf{V}}^{{\rm{s\rightarrow B}}}}\left({{n_{\rm{R}}}}\right)}\right)$.

In (2), the model is derived as follows. First, in (2a), we decompose the multipath channel ${{\bf{H}}^{\rm{S}}}\left({{n_{\rm{R}}}}\right)$ in (1) into free space channel with occlusion ${\bf\tilde{H}}\left({{n_{\rm{R}}}}\right)$ and environmental information $\bm{x}$. Then, in (2b), we decompose free space channel with occlusion ${\bf\tilde{H}}\left({{n_{\rm{R}}}}\right)$ into free space channel ${\bf H}\left({{n_{\rm{R}}}}\right)$ and occlusion matrix ${\bf V}\left({{n_{\rm{R}}}}\right)$. Finally, in (2c), we further decompose the free space channel with occlusion ${\bf\tilde{H}}\left({{n_{\rm{R}}}}\right)$ into the channels from the user, the base station to the point cloud position ${{\bf{H}}^{{\rm{U\rightarrow s}}}}$, ${{\bf{H}}^{{\rm{s\rightarrow B}}}}\left({{n_{\rm{R}}}}\right)$ and their corresponding occlusion matrices ${{\bf{V}}^{{\rm{U\rightarrow s}}}}$, ${{\bf{V}}^{{\rm{s\rightarrow B}}}}\left({{n_{\rm{R}}}}\right)$.

We detect the occlusion effect between scatterers based on the geometric positions between pixels. As shown in Fig. 3, there is a pixel C with a scattering coefficient greater than the threshold $\eta$ between pixels A and B. We believe that pixel C will cause an occlusion effect. Let the position of pixel A be the origin of coordinates, then the position vectors of pixel B and C are expressed as $\bm{b}$ and $\bm{c}$ respectively. According to the following conditions, it is considered that pixel C blocks the line-of-sight path between pixels A and B.

• •

  The distance $d$ between pixel C and the line-of-sight path between pixels A and B is less than the threshold $l$,

  $$\small d=\frac{{\left\lvert{{\bm{c}}\times{\bm{b}}}\right\rvert}}{{\left\lvert{\bm{b}}\right\rvert}}<l,$$

  (3)

  where the detection threshold $l$ should be selected according to the side length or diagonal length of the small cube, that is, $l\in\left[\min\left\{l_{\rm{s}},w_{\rm{s}},h_{\rm{s}}\right\},\left(l_{\rm{s}}^{2}+w_{\rm{s}}^{2}+h_{\rm{s}}^{2}\right)^{1/2}\right]$.

• •

  The angle between vector $\bm{b}$ and vector $\bm{c}$ is an acute angle, which means that pixel C is located in the same direction of the line-of-sight path from pixel A to pixel B, but not pixel A is located between the pixel B and the pixel C. It is expressed as

  $$\small{\bm{b}}\cdot{\bm{c}}>0.$$

  (4)

• •

  Pixel C is located between pixel A and pixel B, which is expressed as

  $$\small\left\lvert{{\bm{c}}\cdot{\bm{b}}}\right\rvert<\left\lvert{\bm{b}}\right\rvert^{2}.$$

  (5)

According to the above occlusion detection method, we design the environment sensing algorithm. However, the proposed occlusion detection method is not the only method, and the proposed environment sensing algorithm does not depend on the uniqueness of the occlusion detection method.

As mentioned in Section I, since the distribution of scatterers in the environment is sparse, according to the CS theory, the optimization problem of solving environmental information is expressed as

$$\small{\bm{\hat{x}}}=\arg\mathop{\min}\limits_{\bm{x}}{\left\|{\bm{x}}\right\|_{1}}\quad{\rm{s}}{\rm{.t}}{\rm{.}}\quad{\left\|{{\bf{Y}}-{\bf{S}}\left({{{\bf{H}}^{\rm{S}}}+{{\bf{H}}^{{\rm{LOS}}}}}\right)}\right\|_{2}}\leq\varepsilon,$$

(6)

where $\varepsilon$ is a slack variable. In the system model of Section II-B, We can calculate the free space channel coefficients ${{\bf{H}}^{{\rm{LOS}}}}\left({{n_{\rm{R}}}}\right)$ and ${{\bf{H}}}\left({{n_{\rm{R}}}}\right)$ through the basic free space path propagation model. In addition, some empirical models, simplified models, or interpolation methods are also suitable, and the appropriate model needs to be selected according to the practical environment. The BS estimates the multipath channel ${{\bf{H}}^{{\rm{S}}}}\left({{n_{\rm{R}}}}\right)$ by detecting the UE pilot and calculating the relationship with the received signal ${\bf{Y}}\left({{n_{\rm{R}}}}\right)$. Therefore, in general, in order to achieve accurate channel estimation (e.g. least-squares channel estimation), the pilot length should be greater than the number of users, that is, $N_{\rm T}>N_{\rm u}$. The proposed algorithm does not rely much on the method of channel estimation. Note that the channel estimation is only the first step of the overall integrated sensing and communication process. The subsequent environment sensing problem studied here is always an underdetermined CS problem with given channel matrix. Due to the large propagation range of wireless signals, according to different precision requirements, the number of pixels $N_{\rm s}$ will be greater than $N_{\rm u}N_{\rm R}$. At the same time, based on the sparsity of the environmental information $\bm{x}$, the optimization problem (6) is expressed as the CS reconstruction problem equation

	$\displaystyle{\left[{\begin{array}[]{*{20}{c}}{{{\bf{H}}^{\rm{S}}}\left(1\right)}\\ \vdots\\ {{{\bf{H}}^{\rm{S}}}\left({{n_{\rm{R}}}}\right)}\\ \vdots\\ {{{\bf{H}}^{\rm{S}}}\left({{N_{\rm{R}}}}\right)}\end{array}}\right]_{{N_{\rm{u}}}{N_{\rm{R}}}\times 1}}$	$\displaystyle={\left[{\begin{array}[]{*{20}{c}}{{\bf{\tilde{H}}}\left(1\right)}\\ \vdots\\ {{\bf{\tilde{H}}}\left({{n_{\rm{R}}}}\right)}\\ \vdots\\ {{\bf{\tilde{H}}}\left({{N_{\rm{R}}}}\right)}\end{array}}\right]_{{N_{\rm{u}}}{N_{\rm{R}}}\times{N_{\rm{s}}}}}{\left[{\bm{x}}\right]_{{N_{\rm{s}}}\times 1}},$		(17)
	$\displaystyle\Rightarrow{{\bf{H}}^{\rm{S}}}$	$\displaystyle={\bf{\tilde{H}}}{\bm{x}},$		(18)

where ${\bf{\tilde{H}}}$ is formed by splicing the multipath channel matrix ${{\bf{\tilde{H}}}\left({{n_{\rm{R}}}}\right)}$ between $N_{\rm{R}}$ receiving antennas and $N_{\rm{u}}$ users. In the multi-BS joint environment sensing scheme, $N_{\rm{R}}$ receiving antennas include the number of all antennas on multiple BSs.

Different from the traditional CS reconstruction problem, solving the environment sensing problem with occlusion effect requires not only solving the environment information ${\bm{x}}$ but also needs to consider the unknown occlusion matrix ${\bf{V}}$. Traditional CS reconstruction algorithms, such as GAMP algorithm [25], can solve environmental information in an ideal (ignoring occlusion, ${\bf{\tilde{H}}}={\bf{H}}$) situation. In this way, since the accurate measurement matrix ${\bf{\tilde{H}}}$ cannot be obtained, the environmental information ${\bm{x}}$ cannot be accurately solved. Compared with the GAMP algorithm, the Bilinear GAMP algorithm [28] can solve the bilinear problem, that is, simultaneously solve the occlusion matrix ${\bf{V}}$ and the environmental information ${\bm{x}}$. However, it does not consider the relationship between the occlusion effect and the environmental information, the Bilinear GAMP algorithm cannot constrain the uniqueness and accuracy of the occlusion matrix ${\bf{V}}$ and the environmental information ${\bm{x}}$.

Fig. 4 shows a sketch of the proposed GAMP-MVSVR algorithm. In Section IV-B, we propose a factor graph based on the decomposition result of the joint posterior probability. Based on the sum-product algorithm (SPA) algorithm, we calculate the message passed in the factor graph, and define the prior probability and likelihood function. In Section IV-C, we approximate the message from factor node to variable node. Through the approximate process, we use the result in the previous iteration to calculate its corresponding multipath channel. We compare it with the channel ${{\bf{\hat{H}}}^{{\rm{S}}}}$ estimated by the pilot and calculate the residual to guide this iteration. In Section IV-D, we approximate the message from variable node to factor node. Through the approximate process, according to the occlusion effect, we estimate ${\bf{\tilde{H}}}$ and ${\bm{x}}$ based on the residual. In Section IV-E, We update the intermediate variables in Section IV-C and IV-D, and achieve the iterative execution of the estimation methods in Section IV-C and IV-D.

As mentioned in Section III-B, let ${{\bf{\hat{H}}}^{\rm{S}}}$ be the multipath channel estimated by the pilot, then ${{\bf{\hat{H}}}^{\rm{S}}}={{\bf{H}}^{\rm{S}}}+{\bf{W}_{\rm e}}$, where ${\bf{W}_{\rm e}}$ represents the estimation error. The estimation error ${\bf{W}_{\rm e}}$ is caused by two main reasons, including the multipath channel estimation error caused by pilot estimation and the free space channel estimation error caused by modeling estimation. In order to simplify the analysis, we assume that the noise ${\bf W}_{\rm e}$ obeys a Gaussian distribution. The joint posterior distribution of the occluded multipath channel ${\bf{\tilde{H}}}$ and the environmental information ${\bm{x}}$ is expressed as

$$\small\begin{split}&{p_{{\bm{x}},{\bf{\tilde{H}}}|{{\bf{\hat{H}}}^{\rm{S}}}}}\left({{\bm{x}},{\bf{\tilde{H}}}|{{\bf{\hat{H}}}^{\rm{S}}}}\right)\\ =&p_{{{\bf{\hat{H}}}^{\rm{S}}}|{\bm{x}},{\bf{\tilde{H}}}}\left({{{\bf{\hat{H}}}^{\rm{S}}}|{\bm{x}},{\bf{\tilde{H}}}}\right){p_{\bm{x}}}\left({\bm{x}}\right){p_{{\bf{\tilde{H}}}|{\bm{x}}}}\left({{\bf{\tilde{H}}}|{\bm{x}}}\right)/{{p_{{\bf{\hat{H}}}^{\rm{S}}}}\left({{\bf{\hat{H}}}^{\rm{S}}}\right)}\\ \propto&p_{{{\bf{\hat{H}}}^{\rm{S}}}|{\bf{\tilde{H}}}{\bm{x}}}\left({{{\bf{\hat{H}}}^{\rm{S}}}|{\bf{\tilde{H}}}{\bm{x}}}\right){p_{\bm{x}}}\left({\bm{x}}\right){p_{{\bf{\tilde{H}}}|{\bm{x}}}}\left({{\bf{\tilde{H}}}|{\bm{x}}}\right)\\ =&p_{{{\bf{\hat{H}}}^{\rm{S}}}|{{\bf{H}}^{\rm{S}}}}\left({{{\bf{\hat{H}}}^{\rm{S}}}|{{\bf{H}}^{\rm{S}}}}\right){p_{\bm{x}}}\left({\bm{x}}\right){p_{{\bf{\tilde{H}}}|{\bm{x}}}}\left({{\bf{\tilde{H}}}|{\bm{x}}}\right),\end{split}$$

(20)

where ${{p_{{\bf{H}}^{\rm{S}}}}\left({{\bf{H}}^{\rm{S}}}\right)}$ is a constant. To simplify the expression, let $N_{\rm{H}}={N_{\rm{u}}}{N_{\rm{R}}}$ denote the length of the multipath channel vector ${{\bf{\hat{H}}}^{\rm{S}}}$. As shown in Fig. 5, according to the (20), the posterior distribution is expressed in the form of a factor graph.

As shown in Fig. 5, the circles represent the variable nodes $\bf{\tilde{H}}$ and $\bm{x}$, the squares represent the factor node $p_{{{\bf{\hat{H}}}^{\rm{S}}}|{{\bf{H}}^{\rm{S}}}}({{{\bf{\hat{H}}}^{\rm{S}}}|{{\bf{H}}^{\rm{S}}}})$, ${p_{\bm{x}}}\left({\bm{x}}\right)$ and ${p_{{\bf{\tilde{H}}}|{\bm{x}}}}({{\bf{\tilde{H}}}|{\bm{x}}})$. The lines between the nodes represent that the variables are included in the function of the factor nodes. Specifically, the environmental information prior probability function ${p_{\bm{x}}}\left({\bm{x}}\right)$ contains the environmental information variable $\bm{x}$. The relationship between the environmental information variable $\bm{x}$ and the channel with occlusion $\bf{\tilde{H}}$ conforms to the occlusion detection model function ${p_{{\bf{\tilde{H}}}|{\bm{x}}}}({{\bf{\tilde{H}}}|{\bm{x}}})$. The multipath channel observation function $p_{{{\bf{\hat{H}}}^{\rm{S}}}|{{\bf{H}}^{\rm{S}}}}({{{\bf{\hat{H}}}^{\rm{S}}}|{{\bf{H}}^{\rm{S}}}})$ contains the environmental information variable $\bm{x}$ and the free space channel $\bf{\tilde{H}}$. The specific function is expressed as follows: ${{\bf{\hat{H}}}^{\rm{S}}}$ can be considered as the observation of ${{\bf{H}}^{\rm{S}}}$ under the interference of additive white Gaussian noise, i.e.,

$$\small p_{{\hat{h}^{\rm{S}}_{n_{\rm{H}}}}|{h_{n_{\rm{H}}}^{\rm{S}}}}\left({{\hat{h}_{n_{\rm{H}}}^{\rm{S}}}|{h_{n_{\rm{H}}}^{\rm{S}}}}\right)={\mathcal{N}}\left({\hat{h}_{n_{\rm{H}}}^{\rm{S}};{h_{n_{\rm{H}}}^{\rm{S}}},{\sigma^{\rm{w}}}}\right),$$

(21)

where ${\sigma^{\rm{w}}}$ is the variance of Gaussian white noise. Due to the sparseness of environmental information $\bm{x}$, we assume that the $\bm{x}$ obeys a Bernoulli-Gaussian distribution,

$$\small{p_{\bm{x}}}\left({\bm{x}}\right)=\left({1-\lambda}\right)\delta\left(x\right)+\lambda{\cal N}\left({x;{\theta^{\rm{x}}},{\sigma^{\rm{x}}}}\right),$$

(22)

where $\lambda$ is the sparsity of the environment, $\theta^{\rm{x}}$ and $\sigma^{\rm{x}}$ represent the mean and variance of the scattering coefficient of the environmental scatterer. ${p_{{\bf{\tilde{H}}}|{\bm{x}}}}\left({{\bf{\tilde{H}}}|{\bm{x}}}\right)$ represents the occlusion effect of environmental information $\bm{x}$ on the channel, i.e.,

$$\small\begin{split}&{p_{{{\tilde{h}}_{{n_{\rm{H}}},{n_{\rm{s}}}}}|{\bm{x}}}}\left({{{\tilde{h}}_{{n_{\rm{H}}},{n_{\rm{s}}}}}|{\bm{x}}}\right)=\\ &\left\{\begin{array}[]{l}\delta\left({{{\tilde{h}}_{{n_{\rm{H}}},{n_{\rm{s}}}}}-{h_{{n_{\rm{H}}},{n_{\rm{s}}}}}}\right)\quad{v_{{n_{\rm{H}}},{n_{\rm{s}}}}}=f_{\rm{occ}}\left({\bm{x}}\right)=1\\ \delta\left({{{\tilde{h}}_{{n_{\rm{H}}},{n_{\rm{s}}}}}}\right)\quad\quad\quad\quad\;\;\;{v_{{n_{\rm{H}}},{n_{\rm{s}}}}}=f_{\rm{occ}}\left({\bm{x}}\right)=0\end{array}\right.,\end{split}$$

(23)

where $f_{\rm{occ}}\left(\cdot\right)$ represents the occlusion detection function, the specific method is as described in the Section III-A, and the occlusion detection threshold is set to $\eta={\theta^{\rm{x}}}/4$.

Based on the sparse prior information and occlusion relationship of environmental information, we aim to obtain the minimum mean square error (MMSE) estimate of environmental information

$$\small\hat{x}_{{n_{\rm{s}}}}^{{\rm{MMSE}}}=\arg\mathop{\min}\limits_{{{\hat{x}}_{{n_{\rm{s}}}}}}{\mathbb{E}_{\bm{x}}}\left[{{{\left({{x_{{n_{\rm{s}}}}}-{{\hat{x}}_{{n_{\rm{s}}}}}}\right)}^{2}}}\right].$$

(24)

According to the loopy belief propagation (LBP) theory [29], based on the factor graph model, we calculate the MMSE estimated value of the environment information through the SPA [30]. As explained in Appendix A, the method of SPA message passing is as follows. As shown in Fig. 5, the output message of the variable node along the edge (solid line) is obtained by calculating the product of all its input messages. The output message of the factor node is calculated by multiplying and integrating the function of the factor node and its input message. The marginal distribution of the variable is calculated from the product of all messages input to the node. Since it is difficult to accurately calculate the passed message , we approximate the message in the following sections.

In this section, we approximate the message from the factor node to the variable node according to the central limit theorem (CLT) and Gaussian approximation based on the idea of approximate message passing [24, 25, 28]. Based on CLT, when $N_{\rm{s}}\rightarrow\infty$, $N_{\rm{H}}/N_{\rm{s}}$ remains constant, we drop some high-order infinitesimal terms. In the $t$ iteration, the sparsity $\lambda$ in the prior distribution of environmental information ${p_{\bm{x}}}\left({\bm{x}}\right)$ remains constant. As mentioned in Section III-A and Section V, the number of zero elements in the occlusion matrix $\bf{V}$ in the variable $\bf{\tilde{H}}$ remains fixed, and the position of the zero elements is determined by ${\hat{x}_{{n_{\rm{s}}}}}$. Therefore, we consider that variable ${\hat{x}_{{n_{\rm{s}}}}}$ is independent of variable ${\hat{h}_{{n_{\rm{H}}},{n_{\rm{s}}}}}$ in the $t$-th iteration. As explained in Appendix B, the passed message in (77) is approximated as

	$\displaystyle{\mu_{\hat{h}_{{n_{\rm{H}}}}^{\rm{S}}\to{x_{{n_{\rm{s}}}}}}}\left({t,{x_{{n_{\rm{s}}}}}}\right)=\log\int_{h_{{n_{\rm{H}}}}^{\rm{S}}}{{p_{\hat{h}_{{n_{\rm{H}}}}^{\rm{S}}|h_{{n_{\rm{H}}}}^{\rm{S}}}}\left({\hat{h}_{{n_{\rm{H}}}}^{\rm{S}}\left|{h_{{n_{\rm{H}}}}^{\rm{S}}}\right.}\right)}$
	$\displaystyle\times{\cal N}\left({h_{{n_{\rm{H}}}}^{\rm{S}};{{\bar{h}}_{{n_{\rm{H}}}{n_{\rm{s}}},{n_{\rm{H}}}}}\left(t\right){x_{{n_{\rm{s}}}}}+{{\bar{p}}_{{n_{\rm{s}}},{n_{\rm{H}}}}}\left(t\right),}\right.$
	$\displaystyle\quad\quad\quad\left.{\sigma_{{n_{\rm{H}}}{n_{\rm{s}}},{n_{\rm{H}}}}^{\rm{h}}\left(t\right)x_{{n_{\rm{s}}}}^{2}+\sigma_{{n_{\rm{s}}},{n_{\rm{H}}}}^{\rm{p}}\left(t\right)}\right)+c$
	$\displaystyle={F_{{n_{\rm{H}}}}}\left({{{\bar{h}}_{{n_{\rm{H}}}{n_{\rm{s}}},{n_{\rm{H}}}}}\left(t\right){x_{{n_{\rm{s}}}}}+{{\bar{p}}_{{n_{\rm{s}}},{n_{\rm{H}}}}}\left(t\right),}\right.$
	$\displaystyle\quad\quad\quad\left.{\sigma_{{n_{\rm{H}}}{n_{\rm{s}}},{n_{\rm{H}}}}^{\rm{h}}\left(t\right)x_{{n_{\rm{s}}}}^{2}+\sigma_{{n_{\rm{s}}},{n_{\rm{H}}}}^{\rm{p}}\left(t\right);\hat{h}_{{n_{\rm{H}}}}^{\rm{S}}}\right)+c,$		(25)

where

$$\small{F_{{n_{\rm{H}}}}}\left({{\bar{p}},{\sigma^{\rm{p}}};{{\hat{h}}^{\rm{S}}}}\right)=\log\int_{{h^{\rm{S}}}}{{p_{{{\hat{h}}^{\rm{S}}}|{h^{\rm{S}}}}}\left({{{\hat{h}}^{\rm{S}}}|{h^{\rm{S}}}}\right){\cal N}\left({{{\hat{h}}^{\rm{S}}};{\bar{p}},{\sigma^{\rm{p}}}}\right)}.$$

(26)

Adding the $n_{\rm{s}}$-th term to the summation in ${\bar{p}_{{n_{\rm{s}}},{n_{\rm{H}}}}}\left(t\right)$ in (86), then ${\bar{p}_{{n_{\rm{H}}}}}\left(t\right)$ and $\sigma_{{n_{\rm{H}}}}^{\rm{p}}\left(t\right)$ are considered as an estimate of ${h_{{n_{\rm{H}}}}^{\rm{S}}}$, i.e.,

$$\small{\bar{p}_{{n_{\rm{H}}}}}\left(t\right)=\sum_{k=1}^{N_{\rm{s}}}{{{\bar{h}}_{{n_{\rm{H}}}k,{n_{\rm{H}}}}}\left(t\right){{\bar{x}}_{k,{n_{\rm{H}}}}}\left(t\right)},$$

(27)

$$\small\begin{split}\sigma_{{n_{\rm{H}}}}^{\rm{p}}&\left(t\right)=\sum_{k=1}^{N_{\rm{s}}}{\left({\bar{h}_{{n_{\rm{H}}}k,{n_{\rm{H}}}}^{2}\left(t\right)\sigma_{k,{n_{\rm{H}}}}^{\rm{x}}\left(t\right)}\right.}\\ &{\left.{+\sigma_{{n_{\rm{H}}}k,{n_{\rm{H}}}}^{\rm{h}}\left(t\right)\bar{x}_{k,{n_{\rm{H}}}}^{2}\left(t\right)+\sigma_{{n_{\rm{H}}}k,{n_{\rm{H}}}}^{\rm{h}}\left(t\right)\sigma_{k,{n_{\rm{H}}}}^{\rm{x}}\left(t\right)}\right)}.\end{split}$$

(28)

We define a residual variable whose mean value is,

$$\small{\bar{s}_{n_{\rm{H}}}}\left(t\right)={F^{\prime}_{{n_{\rm{H}}},1}}\left({{{\bar{p}}_{{n_{\rm{H}}}}}\left(t\right),\sigma_{{n_{\rm{H}}}}^{\rm{p}}\left(t\right);\hat{h}_{{n_{\rm{H}}}}^{\rm{S}}}\right),$$

(29)

due to the Gaussian property of ${F_{{n_{\rm{H}}}}}$, the variance of the residual is expressed as,

$$\small\sigma_{n_{\rm{H}}}^{\rm{s}}\left(t\right)=-{F^{\prime\prime}_{{n_{\rm{H}}},1}}\left({{{\bar{p}}_{{n_{\rm{H}}}}}\left(t\right),\sigma_{{n_{\rm{H}}}}^{\rm{p}}\left(t\right);\hat{h}_{{n_{\rm{H}}}}^{\rm{S}}}\right).$$

(30)

where ${F^{\prime}_{{n_{\rm{H}}},1}}$ and ${F^{\prime\prime}_{{n_{\rm{H}}},1}}$ represent the first and second derivatives of the first argument of ${F_{{n_{\rm{H}}}}}$, and ${F^{\prime}_{{n_{\rm{H}}},2}}$ represents the first derivative of the second argument of ${F_{{n_{\rm{H}}}}}$.

As explained in Appendix C, we get the final approximation result of the passed message ${\mu_{\hat{h}_{{n_{\rm{H}}}}^{\rm{S}}\to{x_{{n_{\rm{s}}}}}}}\left({t,{x_{{n_{\rm{s}}}}}}\right)$ in (77)

	$\displaystyle{\mu_{\hat{h}_{{n_{\rm{H}}}}^{\rm{S}}\to{x_{{n_{\rm{s}}}}}}}\left({t,{x_{{n_{\rm{s}}}}}}\right)$
	$\displaystyle\approx\left[{{{\bar{s}}_{{n_{\rm{H}}}}}\left(t\right){{\bar{h}}_{{n_{\rm{H}}}{n_{\rm{s}}},{n_{\rm{H}}}}}\left(t\right)+\sigma_{{n_{\rm{H}}}}^{\rm{s}}\left(t\right)\hat{h}_{{n_{\rm{H}}},{n_{\rm{s}}}}^{2}\left(t\right){{\hat{x}}_{{n_{\rm{s}}}}}\left(t\right)}\right]{x_{{n_{\rm{s}}}}}$
	$\displaystyle-\frac{1}{2}\left[{\sigma_{{n_{\rm{H}}}}^{\rm{s}}\left(t\right)\hat{h}_{{n_{\rm{H}}},{n_{\rm{s}}}}^{2}\left(t\right)-\sigma_{{n_{\rm{H}}},{n_{\rm{s}}}}^{\rm{h}}\left(t\right)}\right.$
	$\displaystyle\times\left.{\left({\bar{s}_{{n_{\rm{H}}}}^{2}\left(t\right)-\sigma_{{n_{\rm{H}}}}^{\rm{s}}\left(t\right)}\right)}\right]x_{{n_{\rm{s}}}}^{2}+c,$		(31)

According to (27), ${\bar{p}_{{n_{\rm{H}}}}}$ is considered to be an estimate of ${{\bf{H}}^{\rm{S}}}$, we can express the posterior probability of ${{\bf{H}}^{\rm{S}}}$ as the product of the prior probability ${\cal N}\left({h_{{n_{\rm{H}}}}^{\rm{S}};{{\bar{p}}_{{n_{\rm{H}}}}}\left(t\right),\sigma_{{n_{\rm{H}}}}^{\rm{p}}\left(t\right)}\right)$ and the likelihood ${p_{\hat{h}_{{n_{\rm{H}}}}^{\rm{S}}|h_{{n_{\rm{H}}}}^{\rm{S}}}}({\hat{h}_{{n_{\rm{H}}}}^{\rm{S}}|h_{{n_{\rm{H}}}}^{\rm{S}}})$ as follows

$$\small\begin{split}&{p_{h_{{n_{\rm{H}}}}^{\rm{S}}|{p_{{n_{\rm{H}}}}}}}\left({h_{{n_{\rm{H}}}}^{\rm{S}}|{{\bar{p}}_{{n_{\rm{H}}}}}\left(t\right);\sigma_{{n_{\rm{H}}}}^{\rm{p}}\left(t\right)}\right)\\ &=\frac{1}{C}{p_{\hat{h}_{{n_{\rm{H}}}}^{\rm{S}}|h_{{n_{\rm{H}}}}^{\rm{S}}}}\left({\hat{h}_{{n_{\rm{H}}}}^{\rm{S}}|h_{{n_{\rm{H}}}}^{\rm{S}}}\right){\cal N}\left({h_{{n_{\rm{H}}}}^{\rm{S}};{{\bar{p}}_{{n_{\rm{H}}}}}\left(t\right),\sigma_{{n_{\rm{H}}}}^{\rm{p}}\left(t\right)}\right),\end{split}$$

(32)

where $C=\int_{h^{\rm{S}}}{p_{h_{{n_{\rm{H}}}}^{\rm{S}}|{p_{{n_{\rm{H}}}}}}}\left({h^{\rm{S}}|{{\bar{p}}}\left(t\right);\sigma^{\rm{p}}\left(t\right)}\right)$ is a normalization constant. Then the mean and variance of (32) are the mean and variance of calculated ${{\bf{H}}^{\rm{S}}}$ as follows

$$\small\bar{h}_{{n_{\rm{H}}}}^{\rm{S}}\left(t\right)={\mathbb{E}}\left\{{h_{{n_{\rm{H}}}}^{\rm{S}}|{{\bar{p}}_{{n_{\rm{H}}}}}\left(t\right);\sigma_{{n_{\rm{H}}}}^{\rm{p}}\left(t\right)}\right\},$$

(33)

$$\small\sigma_{{n_{\rm{H}}}}^{{\rm{Hs}}}\left(t\right)={\mathbb{D}}\left\{{h_{{n_{\rm{H}}}}^{\rm{S}}|{{\bar{p}}_{{n_{\rm{H}}}}}\left(t\right);\sigma_{{n_{\rm{H}}}}^{\rm{p}}\left(t\right)}\right\},$$

(34)

where $\mathbb{E}(X)=\int{xf(x){\rm d}x}$ represents the mean of $x$ and $\mathbb{D}(X)=\int{\left[x-\mathbb{E}(X)\right]^{2}f(x){\rm d}x}$ represents the variance of $x$.

Pluging (32), (33) and (34) into (29) and (30). Based on the likelihood in (21), according to the derivation in [25, 28], the mean and variance of the residual are expressed as

$$\small{\bar{s}_{{n_{\rm{H}}}}}\left(t\right)=\frac{{\hat{h}_{{n_{\rm{H}}}}^{\rm{S}}-{{\bar{p}}_{{n_{\rm{H}}}}\left(t\right)}}}{{\sigma_{{n_{\rm{H}}}}^{\rm{p}}\left(t\right)+{\sigma^{\rm{w}}}}},$$

(35)

$$\small\sigma_{{n_{\rm{H}}}}^{\rm s}\left(t\right)=\frac{1}{{\sigma_{{n_{\rm{H}}}}^{\rm{p}}\left(t\right)+{\sigma^{\rm{w}}}}}.$$

(36)

As mentioned in this section, in the $t$-th iteration, we consider that variable $\bm{x}$ is independent of variable $\bf{\tilde{H}}$, so the same as the derivation of (31), we approximate the message ${\mu_{\hat{h}_{{n_{\rm{H}}}}^{\rm{S}}\to{{\tilde{h}}_{{n_{\rm{H}}},{n_{\rm{s}}}}}}}({t,{{\tilde{h}}_{{n_{\rm{H}}},{n_{\rm{s}}}}}})$ in (79) as

	$\displaystyle{\mu_{\hat{h}_{{n_{\rm{H}}}}^{\rm{S}}\to{{\tilde{h}}_{{n_{\rm{H}}},{n_{\rm{s}}}}}}}\left({t,{{\tilde{h}}_{{n_{\rm{H}}},{n_{\rm{s}}}}}}\right)\approx$
	$\displaystyle\left[{{{\bar{s}}_{{n_{\rm{H}}}}}\left(t\right){{\bar{x}}_{{n_{\rm{s}}},{n_{\rm{H}}}}}\left(t\right)+\sigma_{{n_{\rm{H}}}}^{\rm{s}}\left(t\right)\hat{x}_{{n_{\rm{s}}}}^{2}\left(t\right){{\hat{h}}_{{n_{\rm{H}}},{n_{\rm{s}}}}}\left(t\right)}\right]{{\tilde{h}}_{{n_{\rm{H}}},{n_{\rm{s}}}}}$
	$\displaystyle-\frac{1}{2}\left[{\sigma_{{n_{\rm{H}}}}^{\rm{s}}\left(t\right)\hat{x}_{{n_{\rm{s}}}}^{2}\left(t\right)-\sigma_{{n_{\rm{s}}}}^{\rm{x}}\left(t\right)}\right.$
	$\displaystyle\times\left.{\left({\bar{s}_{{n_{\rm{H}}}}^{2}\left(t\right)-\sigma_{{n_{\rm{H}}}}^{\rm{s}}\left(t\right)}\right)}\right]\tilde{h}_{{n_{\rm{H}}},{n_{\rm{s}}}}^{2}+c.$		(37)

In summary, in this section we approximate the message from factor node to variable node as follows. We calculate ${{{\bar{p}}_{{n_{\rm{H}}}}}\left(t\right)}$ (as shown in (27)) to obtain the calculated channel ${\bar{h}_{{n_{\rm{H}}}}^{\rm{S}}\left(t\right)}$ (as shown in (33)). Then we calculate the residual ${\bar{s}_{{n_{\rm{H}}}}}\left(t\right)$ according to ${\bar{h}_{{n_{\rm{H}}}}^{\rm{S}}\left(t\right)}$ and the channel estimation value using the pilot ${\hat{h}_{{n_{\rm{H}}}}^{\rm{S}}}$ (as shown in (35)). Finally, we use the residual and the approximate mean and variance of the passed message to represent the message from the factor node to the variable node (as shown in (31) and (37)).

In this section, we approximate the message from the factor node to the variable node. We plug (31) in (78), the message ${{\mu_{{x_{r}}\to\hat{h}_{{n_{\rm{H}}}}^{\rm{S}}}}\left({t+1,{x_{r}}}\right)}$ from variable node $\bm{x}$ to factor node $p_{{{\bf{\hat{H}}}^{\rm{S}}}|{{\bf{H}}^{\rm{S}}}}({{{\bf{\hat{H}}}^{\rm{S}}}|{{\bf{H}}^{\rm{S}}}})$ is approximated as

$$\small\begin{split}&{\mu_{{x_{{n_{\rm{s}}}}}\to\hat{h}_{{n_{\rm{H}}}}^{\rm{S}}}}\left({t+1,{x_{{n_{\rm{s}}}}}}\right)\\ &=\log{p_{\bm{x}}}\left({{x_{{n_{\rm{s}}}}}}\right)-\frac{1}{{2\sigma_{{n_{\rm{s}}},{n_{\rm{H}}}}^{\rm{r}}\left(t\right)}}{\left({{x_{{n_{\rm{s}}}}}-{{\bar{r}}_{{n_{\rm{s}}},{n_{\rm{H}}}}}\left(t\right)}\right)^{2}}+c\\ &=\log\left({{p_{\bm{x}}}\left({{x_{{n_{\rm{s}}}}}}\right){\cal N}\left({{x_{{n_{\rm{s}}}}};{{\bar{r}}_{{n_{\rm{s}}},{n_{\rm{H}}}}}\left(t\right),\sigma_{{n_{\rm{s}}},{n_{\rm{H}}}}^{\rm{r}}\left(t\right)}\right)}\right)+c,\end{split}$$

(38)

where

	$\displaystyle{{\bar{r}}_{{n_{\rm{H}}},{n_{\rm{s}}}}}\left(t\right)$
	$\displaystyle={{\hat{x}}_{{n_{\rm{s}}}}}\left(t\right)\left({1+\sigma_{{n_{\rm{s}}},{n_{\rm{H}}}}^{\rm{r}}\left(t\right)\sum\limits_{k\neq{n_{\rm{H}}}}{\sigma_{k,{n_{\rm{s}}}}^{\rm{h}}\left(t\right)\left[{\bar{s}_{k}^{2}\left(t\right)-\sigma_{k}^{\rm{s}}\left(t\right)}\right]}}\right)$
	$\displaystyle+\sigma_{{n_{\rm{s}}},{n_{\rm{H}}}}^{\rm{r}}\left(t\right)\sum\limits_{k\neq n_{\rm{H}}}{{{\hat{h}}_{k,{n_{\rm{s}}}}}\left(t\right){{\bar{s}}_{k}}\left(t\right)},$		(39)

$$\small\begin{split}&\sigma_{{n_{\rm{s}}},{n_{\rm{H}}}}^{\rm{r}}\left(t\right)\\ &={\left({\sum\limits_{k\neq{n_{\rm{H}}}}{{{\hat{h}}_{k,{n_{\rm{s}}}}}\left(t\right)\sigma_{k}^{\rm{s}}\left(t\right)-\sigma_{k,{n_{\rm{s}}}}^{\rm{h}}\left(t\right)\left({\bar{s}_{k}^{2}\left(t\right)-\sigma_{k}^{\rm{s}}\left(t\right)}\right)}}\right)^{-1}}.\end{split}$$

(40)

Then the mean ${\bar{x}_{{n_{\rm{s}}},{n_{\rm{H}}}}}\left({t+1}\right)$ and variance $\sigma_{{n_{\rm{s}}},{n_{\rm{H}}}}^{\rm{x}}\left({t+1}\right)$ of the passed message in (38) are expressed as,

$$\small\begin{split}&{{\bar{x}}_{{n_{\rm{s}}},{n_{\rm{H}}}}}\left({t+1}\right)\\ &=\frac{1}{C}\int_{{x_{{n_{\rm{s}}}}}}{{x_{{n_{\rm{s}}}}}{p_{\bm{x}}}\left({{x_{{n_{\rm{s}}}}}}\right){\cal N}\left({{x_{{n_{\rm{s}}}}};{{\bar{r}}_{{n_{\rm{H}}},{n_{\rm{s}}}}}\left(t\right),\sigma_{{n_{\rm{H}}},{n_{\rm{s}}}}^{\rm{r}}\left(t\right)}\right)}\\ &={G_{\bm{x}}}\left({{{\bar{r}}_{{n_{\rm{H}}},{n_{\rm{s}}}}}\left(t\right),\sigma_{{n_{\rm{H}}},{n_{\rm{s}}}}^{\rm{r}}\left(t\right)}\right),\end{split}$$

(41)

$$\small\sigma_{{n_{\rm{s}}},{n_{\rm{H}}}}^{\rm{x}}\left({t+1}\right)=\sigma_{{n_{\rm{H}}},{n_{\rm{s}}}}^{\rm{r}}\left(t\right){G^{\prime}_{\bm{x}}}\left({{{\bar{r}}_{{n_{\rm{H}}},{n_{\rm{s}}}}}\left(t\right),\sigma_{{n_{\rm{H}}},{n_{\rm{s}}}}^{\rm{r}}\left(t\right)}\right),$$

(42)

where $C=\int_{{x_{{n_{\rm{s}}}}}}{{p_{\bm{x}}}\left({{x_{{n_{\rm{s}}}}}}\right){\cal N}\left({{x_{{n_{\rm{s}}}}};{{\bar{r}}_{{n_{\rm{H}}},{n_{\rm{s}}}}}\left(t\right),\sigma_{{n_{\rm{H}}},{n_{\rm{s}}}}^{\rm{r}}\left(t\right)}\right)}$ is a normalization constant, and ${G^{\prime}_{\bm{x}}}$ represents the first derivative of the first argument of the ${{G}_{\bm{x}}}$. The specific derivation of the relationship between (41) and (42) was shown in [25].

Adding the $n_{\rm{H}}$-th term to the summation in ${{\bar{r}}_{{n_{\rm{H}}},{n_{\rm{s}}}}}\left(t\right)$, then ${{\bar{r}}_{{n_{\rm{s}}}}}\left(t\right)$ and $\sigma_{{n_{\rm{s}}}}^{\rm{r}}\left(t\right)$ are expressed as

	$\displaystyle{{\bar{r}}_{{n_{\rm{s}}}}}\left(t\right)$	$\displaystyle={{\hat{x}}_{{n_{\rm{s}}}}}\left(t\right)\left({1+\sigma_{{n_{\rm{s}}}}^{\rm{r}}\left(t\right)\sum_{k=1}^{N_{\rm{H}}}{\sigma_{k,{n_{\rm{s}}}}^{\rm{h}}\left(t\right)\left[{\bar{s}_{k}^{2}\left(t\right)-\sigma_{k}^{\rm{s}}\left(t\right)}\right]}}\right)$
		$\displaystyle+\sigma_{{n_{\rm{s}}}}^{\rm{r}}\left(t\right)\sum\limits_{k=1}^{N_{\rm{H}}}{{{\hat{h}}_{k,{n_{\rm{s}}}}}\left(t\right){{\bar{s}}_{k}}\left(t\right)},$		(43)

$\displaystyle\sigma_{{n_{\rm{s}}}}^{\rm{r}}\left(t\right)$

$\displaystyle={\left({\sum_{k=1}^{N_{\rm{H}}}{{{\hat{h}}_{k,{n_{\rm{s}}}}}\left(t\right)\sigma_{k}^{\rm{s}}\left(t\right)}}\right.}{\left.{{-\sigma_{k,{n_{\rm{s}}}}^{\rm{h}}\left(t\right)\left({\bar{s}_{k}^{2}\left(t\right)-\sigma_{k}^{\rm{s}}\left(t\right)}\right)}}\right)^{-1}},$

(44)

Based on AMP [24, 25, 28], ${{\bar{r}}_{{n_{\rm{s}}}}}\left(t\right)$ and $\sigma_{{n_{\rm{s}}}}^{\rm{r}}\left(t\right)$ are considered as an estimate of $x_{n_{\rm{s}}}$ as follows

$$\small{\hat{x}}_{n_{\rm{s}}}\left(t+1\right)={G_{\bm{x}}}\left({{{\bar{r}}_{{n_{\rm{s}}}}}\left(t\right),\sigma_{{n_{\rm{s}}}}^{\rm{r}}\left(t\right)}\right),$$

(45)

$$\small\sigma_{{n_{\rm{s}}}}^{\rm{x}}\left(t+1\right)=\sigma_{{n_{\rm{s}}}}^{\rm{r}}\left(t\right){G^{\prime}_{\bm{x}}}\left({{{\bar{r}}_{{n_{\rm{s}}}}}\left(t\right),\sigma_{{n_{\rm{s}}}}^{\rm{r}}\left(t\right)}\right),$$

(46)

where ${G^{\prime}_{\bm{x}}}$ represent the first derivative of the first argument of ${G_{\bm{x}}}$. Finally, we calculate the relationship between the mean of the passed message ${\bar{x}_{{n_{\rm{s}}},{n_{\rm{H}}}}}\left({t+1}\right)$ and the mean of the marginal distribution ${\hat{x}_{{n_{\rm{s}}}}}\left({t+1}\right)$ as follows

$$\small\begin{split}&{{\bar{x}}_{{n_{\rm{s}}},{n_{\rm{H}}}}}\left({t+1}\right)\\ &\approx{G_{\bm{x}}}\left({{{\bar{r}}_{{n_{\rm{s}}}}}\left(t\right)-\sigma_{{n_{\rm{s}}}}^{\rm{r}}\left(t\right){{\hat{h}}_{{n_{\rm{H}}},{n_{\rm{s}}}}}\left(t\right){{\bar{s}}_{{n_{\rm{H}}}}}\left(t\right),\sigma_{{n_{\rm{s}}}}^{\rm{r}}\left(t\right)}\right)\\ &\approx{G_{\bm{x}}}\left({{{\bar{r}}_{{n_{\rm{s}}}}}\left(t\right),\sigma_{{n_{\rm{s}}}}^{\rm{r}}\left(t\right)}\right)\\ &-\sigma_{{n_{\rm{s}}}}^{\rm{r}}\left(t\right){{\hat{h}}_{{n_{\rm{H}}},{n_{\rm{s}}}}}\left(t\right){{\bar{s}}_{{n_{\rm{H}}}}}\left(t\right){{G^{\prime}}_{\bm{x}}}\left({{{\bar{r}}_{{n_{\rm{s}}}}}\left(t\right),\sigma_{{n_{\rm{s}}}}^{\rm{r}}\left(t\right)}\right)\\ &={{\hat{x}}_{{n_{\rm{s}}}}}\left({t+1}\right)-{{\hat{h}}_{{n_{\rm{H}}},{n_{\rm{s}}}}}\left(t\right){{\bar{s}}_{{n_{\rm{H}}}}}\left(t\right)\sigma_{n}^{\rm{x}}\left({t+1}\right),\end{split}$$

(47)

in which we drop infinitesimal terms such as $\sigma_{{n_{\rm{s}}},{n_{\rm{H}}}}^{\rm{r}}\left(t\right)-{\sigma_{{n_{\rm{s}}}}^{\rm{r}}\left(t\right)}$ in the second step of approximation, and we perform Taylor series expansion on the first argument of the function at point ${{\bar{r}}_{{n_{\rm{s}}}}}\left(t\right)$ and drop the infinitesimal term in the second step of approximation.

In the $t$-th iteration, we consider that variable $\bm{x}$ is independent of variable $\bf{\tilde{H}}$ and recalculate the relationship between variables $\bf{\tilde{H}}$ and $\bm{x}$ after each iteration. Therefore, we approximate the passed message from variable node $\bf{\tilde{H}}$ to factor node $p_{{{\bf{\hat{H}}}^{\rm{S}}}|{{\bf{H}}^{\rm{S}}}}({{{\bf{\hat{H}}}^{\rm{S}}}|{{\bf{H}}^{\rm{S}}}})$ the same way as (38), and calculate the mean ${\hat{h}_{{n_{\rm{H}}},{n_{\rm{s}}}}}\left({t+1}\right)$ and the variance $\sigma_{{n_{\rm{H}}},{n_{\rm{s}}}}^{\rm{h}}\left({t+1}\right)$ of the marginal distribution of variable $\bf{\tilde{H}}$ the same way as (45) and (46), i.e.,

$$\small{\hat{h}_{{n_{\rm{H}}},{n_{\rm{s}}}}}\left({t+1}\right)={G_{\bf{H}}}\left({{{\bar{q}}_{{n_{\rm{H}}},{n_{\rm{s}}}}}\left(t\right),\sigma_{{n_{\rm{H}}},{n_{\rm{s}}}}^{\rm{q}}\left(t\right)}\right),$$

(48)

$$\small\sigma_{{n_{\rm{H}}},{n_{\rm{s}}}}^{\rm{h}}\left({t+1}\right)=\sigma_{{n_{\rm{H}}},{n_{\rm{s}}}}^{\rm{q}}\left(t\right){G^{\prime}_{\bf{H}}}\left({{{\bar{q}}_{{n_{\rm{H}}},{n_{\rm{s}}}}}\left(t\right),\sigma_{{n_{\rm{H}}},{n_{\rm{s}}}}^{\rm{q}}\left(t\right)}\right),$$

(49)

$$\small\begin{split}&{G_{\bf{H}}}\left({{{\bar{q}}_{{n_{\rm{H}}},{n_{\rm{s}}}}}\left(t\right),\sigma_{{n_{\rm{H}}},{n_{\rm{s}}}}^{\rm{q}}\left(t\right)}\right)\\ &\quad\quad\quad\quad\quad=\frac{1}{C}\int_{\tilde{h}}{{\tilde{h}}_{{n_{\rm{H}}},{n_{\rm{s}}}}}{p_{{{\tilde{h}}_{{n_{\rm{H}}},{n_{\rm{s}}}}}|{\bm{x}}}}\left({{{\tilde{h}}_{{n_{\rm{H}}},{n_{\rm{s}}}}}|{\bm{x}}}\right)\\ &\quad\quad\quad\quad\quad\times{\cal N}\left({{{\tilde{h}}_{{n_{\rm{H}}},{n_{\rm{s}}}}};{{\bar{q}}_{{n_{\rm{H}}},{n_{\rm{s}}}}}\left(t\right),\sigma_{{n_{\rm{H}}},{n_{\rm{s}}}}^{\rm{q}}\left(t\right)}\right),\end{split}$$

(50)

where $C=\int_{\tilde{h}}{{p_{{{\tilde{h}}}|{\bm{x}}}}({{{\tilde{h}}}|{\bm{x}}}){\cal N}}({{{\tilde{h}}};{{\bar{q}}}\left(t\right),\sigma^{\rm{q}}\left(t\right)})$ is a normalization constant, ${G^{\prime}_{\bf{H}}}$ represent the first derivative of the first argument of ${G_{\bf{H}}}$, ${\bar{q}_{{n_{\rm{H}}},{n_{\rm{s}}}}}\left(t\right)$ and $\sigma_{{n_{\rm{H}}},{n_{\rm{s}}}}^{\rm{q}}\left(t\right)$ are considered as an estimate of ${{\tilde{h}}_{{n_{\rm{H}}},{n_{\rm{s}}}}}$ as follows

	$\displaystyle{\bar{q}_{{n_{\rm{H}}},{n_{\rm{s}}}}}\left(t\right)={\hat{h}_{{n_{\rm{H}}},{n_{\rm{s}}}}}\left(t\right)\left({1+\sigma_{{n_{\rm{H}}},{n_{\rm{s}}}}^{\rm{q}}\left(t\right)\sigma_{{n_{\rm{s}}}}^{\rm{x}}\left(t\right)}\right.$
	$\displaystyle\left.{\left({\bar{s}_{{n_{\rm{H}}}}^{2}\left(t\right)-\sigma_{{n_{\rm{H}}}}^{\rm{s}}\left(t\right)}\right)+\sigma_{{n_{\rm{H}}},{n_{\rm{s}}}}^{\rm{q}}\left(t\right){{\bar{x}}_{{n_{\rm{s}}},{n_{H}}}}\left(t\right){{\bar{s}}_{{n_{H}}}}\left(t\right)}\right),$		(51)

$$\small\sigma_{{n_{\rm{H}}},{n_{\rm{s}}}}^{\rm{q}}\left(t\right)={\left({\hat{x}_{{n_{\rm{s}}}}^{2}\left(t\right)\sigma_{{n_{\rm{H}}}}^{\rm{s}}\left(t\right)-\sigma_{{n_{\rm{s}}}}^{\rm{x}}\left(t\right)\left({\bar{s}_{{n_{\rm{H}}}}^{2}\left(t\right)-\sigma_{{n_{\rm{H}}}}^{\rm{s}}\left(t\right)}\right)}\right)^{-1}}.$$

(52)