Sparse Bayesian Learning-Based 3D Spectrum Environment Map Construction —— Sampling Optimization, Scenario-Dependent Dictionary Construction and Sparse Recovery

As shown in Fig.1, the region of interest (ROI) is discretized into several small cubes. Each cube is colored according to its RSS, where red cubes represent high RSS values and blue cubes represent low RSS values. The ROI constitutes a spectrum tensor $\bm{\chi}\in{\Re^{{N_{x}}\times{N_{y}}\times{N_{z}}}}$ in the 3D space, where ${N_{x}}$, ${N_{y}}$, ${N_{z}}$ indicate the grid number along $x$, $y$, $z$ dimensions, respectively. Technically, 3D SEM construction in this paper aims to recovery RSS values of all $N={N_{x}}\times{N_{y}}\times{N_{z}}$ cubes based on the known RSS values of sampling cubes.

A sparse signal vector $\bm{\omega}={\left[{{\omega_{1}},{\omega_{2}},\ldots,{\omega_{n}},\ldots,{\omega_{N}}}\right]^{\text{T}}}\in{\mathbb{R}^{N\times 1}}$ can be defined as

where $P_{n}^{t}$ is the transmitting power of the $n$th RF transmitter. $\omega$ is a $K$-sparse signal vector with $\left\|\bm{\omega}\right\|_{0}=K$. That is if there are $K$ stationary RF transmitters denoted as $\left\{{{{\mathbf{T}}_{k}}}\right\}_{k=1}^{K}$, where ${{\mathbf{T}}_{k}}=\left({x_{k}^{t},y_{k}^{t},z_{k}^{t}}\right)$ is the location of $k$th RF transmitter, we have $K\left({K\ll N}\right)$ nonzero elements in $\omega$.

Suppose that we select $M$ SLs from all $N$ cubes denoted by $\left\{{{{\mathbf{S}}_{m}}}\right\}_{m=1}^{M}$ and ${{\mathbf{S}}_{m}}=\left({x_{m}^{s},y_{m}^{s},z_{m}^{s}}\right)$. The sampling rate is $r=M/N$. The Euclidean distance from the $k$th transmitter to the $m$th SL can be written as

Under the line-of-sight (LOS) condition, the RSS from the $k$th transmitter to the $m$th SL can be approximated by using the free space propagation model as

where ${G_{t}}$ and ${G_{r}}$ are the antenna gain of transceivers, $P_{k}^{t}$ denotes the transmitting power, $\lambda$ is the wavelength of carrier, and $\eta$ is the path loss exponent. Since the RSS sensed at each SL may include the receiving power of several RF transmitters [15], the total RSS can be approximately expressed as

It should be mentioned that the realistic propagation model is very complex due to reflection, diffraction and other propagation phenomena. The above data recovery method for the unsampled cubes is only suitable for the LOS propagation scenarios.

As shown in Fig.2 (a), the spectrum tensor $\bm{\chi}$ is firstly vectorized into ${\mathbf{X}}\in{\mathbb{R}^{N\times 1}}$. Compared with the total number of discretized cubes in the ROI, the number of stationary RF transmitters is much small ($K\ll N$). Since the spectrum vector ${\mathbf{X}}$ has high spatial correlation, it can be represented by the product of a sparse dictionary $\bm{\varphi}\in{\mathbb{R}^{N\times N}}$ and the sparse signal $\omega$ as

where the element ${\varphi_{i,j}}$ of sparse dictionary is defined as the channel gain or propagation path loss between the $i$th and the $j$th cubes. Let us set the RSS vector of SLs as $\bm{t}\in{\mathbb{R}^{M\times 1}}$. The measurement matrix $\bm{\psi}\in{\mathbb{R}^{M\times N}}$ can be defined as

where each row of $\bm{\psi}$ has an element of 1 denoting the SL’s position in the ROI. Then, we can have

where $\bm{\varepsilon}\in{\mathbb{R}^{M\times 1}}$ is the measurement noise that obeys the zero-mean Gaussian distribution with variance ${\sigma^{2}}$, and ${\mathbf{\Phi}}$ is a sensing matrix.

This paper recovers the spectrum data by two steps. Firstly, the sparse signal $\omega$ is recovered by the noisy measurement vector $\bm{t}$ and the sensing matrix ${\mathbf{\Phi}}$. Then, the spectrum vector ${\mathbf{X}}$ is constructed according to (5). The sparse signal recovery (SSR) is equivalent to solving the following ${l_{0}}$-minimization problem as

The ${l_{0}}$-minimization problem of (8) is NP-hard, and can be equivalent to a ${l_{1}}$-minimization problem as [28]

The authors in [29] have proved that if ${\mathbf{\Phi}}$ satisfies the condition of restricted isometry property, $\omega$ can be recovered by applying CS algorithms. In this paper, we focus on using the sparse Bayesian theory to solve the recovery problem of the SEM construction.

The proposed 3D SEM construction scheme is illustrated in Fig. 3, which mainly contains the scenario-dependent dictionary construction, the measurement matrix optimization and the 3D SEM construction. Firstly, combined with RT technology and interpolation algorithm, we take the factor of scenario into account and build a scenario-dependent sparse dictionary $\bm{\varphi}$. According to the maximum mutual information criterion, we design the selection scheme of SLs and obtain the measurement matrix $\bm{\psi}$. Then, based on the sparse dictionary and measurement matrix, the sparse signal $\omega$ can be recovered based on the SBL algorithm. Finally, we utilize the sparse dictionary and sparse signal to construct the full SEM.

In order to recover the sparse signal $\omega$ for the SEM construction, we adopt the SBL theory and hierarchical sparse probabilistic model as follows.

1) Noise Model: The sparse regression model (7) is defined in a zero-mean Gaussian noise $\bm{\varepsilon}$ with unknown variance ${\sigma^{2}}$, and then the Gaussian likelihood of $\bm{t}$ can be written as

A Gamma distribution is then posed on $\beta\left({\beta=\left({\sigma^{2}}\right)^{-1}}\right)$, which is a conjugate prior to the Gaussian distribution and can greatly simplifies the analysis [30], as

with

where $c>0$ is the shape parameter, and $d>0$ is the scale parameter, $\Gamma(c)=\int_{0}^{\infty}{{t^{c-1}}{e^{-t}}dt}$.

2) Hierarchical Sparse Prior Model: To induce the sparsity of $\omega$, we deploy a sparseness-promoting prior on it. In the Relevance Vector Machine, $\omega$ is a two-layer hierarchical sparse prior, which shares the same properties with Laplace prior while enabling convenient computation [24]. Specifically, each element of $\omega$ is first posed a zero-mean Gaussian prior

where $\bm{\alpha}=\left[{{\alpha_{1}},{\alpha_{2}},\ldots,{\alpha_{N}}}\right]^{\text{T}}$. Then, to complete the specification of this hierarchical prior, we consider the Gamma hyperpriors over $\bm{\alpha}$ as

The graphical model is shown in Fig.4. The overall prior $p\left(\bm{\omega}\right)$ can be obtained by computing the marginal integral of hyper-parameters in $\bm{\alpha}$ as

Since the integral is computable for Gamma $p\left(\bm{\alpha}\right)$, the true prior $p\left(\bm{\omega}\right)$ is a Student-t distribution which can promote sparsity on $\omega$[31].

3) Sparse Bayesian Inference: Following the Bayesian inference, the posterior distribution over all unknowns is desired as

with the decomposition of the ‘weight posterior’ and the ‘hyper-parameter posterior’. It can be inferred that the weight posterior of $\omega$ is Gaussian

with

where ${\mathbf{\Lambda}}={\text{diag}}\left[{{\alpha_{1}},{\alpha_{2}},\ldots,{\alpha_{N}}}\right]$. The SBL considers the signal recovery from the perspective of statistics. With sparse prior of $\omega$ and compressive samples $\bm{t}$, the posteriori probability density function of sparse vector $\omega$ can be inferred. We further estimate the $\omega$ with the mean $\bm{\mu}$ and evaluate the accuracy of the recovery by the variance ${\mathbf{\Sigma}}$. To calculate $\bm{\mu}$ and ${\mathbf{\Sigma}}$, we estimate probabilistic model hyperparameters $\bm{\alpha}$ and $\beta$, the details of estimation will be discussed in the Section III-D.

The traditional sparse dictionary usually adopts the free-space propagation model without considering the scene information, which is not suitable for urban environments with a lot of buildings. The propagation occurs direct, reflection, and diffraction phenomena in the actual environment. Accordingly, a scenario-dependent dictionary is constructed here by analyzing the characteristics of propagation channel. The RT technique is based on Geometrical Optics and the Uniform Theory of Diffraction. It has been used to predict all the possible propagation path parameters in a given geographic map [32].

Assume that the position of $m$th transmitting cube and $n$th receiving cube are $\bm{m}=[{m_{x}},{m_{y}},{m_{z}}]$ and $\bm{n}=[{n_{x}},{n_{y}},{n_{z}}]$, respectively. The propagation distance of direct path is obtained as

For the indirect path, we define the intersection coordinates of scatterers as $h=({H_{x}},{H_{y}},{H_{z}})$. Thus, the Euclidean distance of $s$th path ray between $m$ and the scatterer, and the Euclidean distance of $s$th ray between the scatterer and $m$ can be respectively calculated as

The proposed RT-based dictionary construction method includes three steps, i.e., decomposition of ray source, tracking rays, and superposition of the filed strength. Firstly, the ray source is decomposed with direct ray, reflection ray and diffraction ray. If the ray arrives at the receiving field position $n$ in LOS propagation, the field intensity of $m$ arriving at $n$ is

where $l$ is the wave number. $E_{1{\text{m}}}$ is the electric field intensity of $1{\text{ m}}$ away from $m$, and $d_{mn}$ is the propagation distance of the direct path in (19). For the reflected ray, the electric field intensity can be expressed as

where $R$ is the reflection coefficient. The electric field intensity of the diffraction path can be expressed as

where $D$ is the diffraction coefficient.

Then, according to (22)-(24), the final received field strength of $n$ can be obtained by vector superposition of all the field strengths

where ${N_{p}}$ is the total number of effective rays. ${{\mathbf{E}}_{s}}$ is the electric field intensity of direct path, reflection path or diffraction path. It should be mentioned that the direct ray disappears when the link between $m$ and $n$ is blocked. Then the parameter calculation of direct ray can be ignored [32]. Therefore, the total average RSS of $n$ is

where $\lambda$ is the wavelength. ${G^{m}}$ and ${G^{n}}$ are the antenna gains of the transmitter at $m$ and receiver at $n$ respectively. Thus, the channel gain in dB between $m$ and $n$ can be calculated, i.e., the element ${\varphi_{mn}}$ in the $m$th row and the $n$th column of matrix $\bm{\varphi}$, as

where ${p^{m}}$ is transmitting power of the $m$th transmitting cube.

Finally, we can predict the channel gain between any two cubes in the ROI and construct a scenario-dependent dictionary matrix. In this paper, we only calculate a small part of dictionary matrix and achieve the whole matrix by interpolation. The inverse distance weighted interpolation is adopted [33]. It assumes that each known value has a local influence with respect to the distance, which is consistent with the radio propagation principle.

Let us set the data obtained by RT as ${\varphi_{g}},g=1,2,\ldots,{N_{0}}$, corresponding to the element in the ${g_{x}}$ th row and the ${g_{y}}$ th column of matrix $\bm{\varphi}$. The unknown value ${\varphi_{ij}}$ in the $i$th row and the $j$th column can be obtained by

where $d_{ij}^{g}=\sqrt{\left({i-{g_{x}}}\right)^{2}+\left({j-{g_{y}}}\right)^{2}}$ is the distance between interpolation point and the known point, and $p$ is the distance exponent.

The layout of sampling sensors, i.e., the measurement matrix $\bm{\psi}$, has a significant effect on the SEM construction accuracy. In [34, 35], the authors used the mutual information between the predicted sensors’ observation and the current target location distribution to optimize the layout.

Different from traditional methods using the random selection, we model the sparse signal reconstruction problem as an information theory problem in communication channel, where sparse signal $\omega$ is the input to the channel and the SLs’ RSS samples $t$ is the output. The SLs are tasked to observe in order to increase the information (or to reduce the uncertainty) about the $\omega$’s state. We select SLs based on the MMI criterion.

We define the index set of candidate SLs for selection is ${\mathbf{S}}$. The subset of indices for the determined SLs is expressed as ${{\mathbf{S}}_{k}}$. When the SLs have different observation angles and perceptual uncertainties, the information gain attributable to different SLs can be quite different [36, 37]. Then, the SLs set ${\mathbf{\varsigma}}$ $\left({{\mathbf{\varsigma}}\subset{\mathbf{S}},\left|{\mathbf{\varsigma}}\right|=M}\right)$ is chosen when their observations ${\bm{t}_{\mathbf{\varsigma}}}$ minimizes the expected conditional entropy of the posterior distribution of $\omega$, as given by

which is equivalent to maximize the entropy reduction of $\omega$. We maximize mutual information

The derivation of formula (31) is given in AppendexA. In the initialization stage, there is no prior information related to the sparse signal $\omega$. According to SBL, we assign the same variance $\alpha$, i.e., 1. Then, the formula (31) can be further converted to

with

where ${\raise 3.01385pt\hbox{$\alpha$}\!\mathord{\left/{\vphantom{\alpha\beta}}\right.\kern-1.2pt}\!\lower 3.01385pt\hbox{$\beta$}}$ is a sufficiently small number denoted the ratio of noise variance to sparse signal variance [36]. Therefore, by ignoring the constant terms, the final objective function of MMI-based sampling is asymptotically approaches to

where ${\mathbf{\Phi}}=\bm{\psi}\bm{\varphi}$. $\psi=\left({\mathbf{I}_{N}}\right)_{\mathbf{\varsigma}\bm{.}}$ is consisted of the rows indexed by set ${\mathbf{\varsigma}}$ in ${{\mathbf{I}}_{N}}$. Accordingly, the problem can be further expressed as

We solve the above problem by the greedy algorithm. Let ${{\mathbf{\Phi}}_{t}}$ denote the sensing matrix after the $t$th SL’s selection as

where ${i_{t}}\left({{i_{t}}\in{\mathbf{\varsigma}},{i_{t}}\in{\mathbf{S}}}\right)$ is the index of the $t$th selected SL and ${\bm{\varphi}_{{i_{t}}}}$ is the ${i_{t}}$th row vector of the sparse dictionary matrix $\bm{\varphi}$. The step-by-step maximization process is considered by the greedy method. The matrix can be expanded as

Therefore, the $t$th row of $\bm{\psi}$ is determined according to the following formula,

After $M$ recursive iterations, the optimal $M$ SLs can be finally selected. The measurement matrix $\bm{\psi}$ is obtained.

In Section II-B, we recover $\omega$ with the mean $\mu$ and evaluate the recovery accuracy by the variance ${\mathbf{\Sigma}}$. The hyperparameters are estimated by maximum a posterior (MAP) probability [24] as

which is equivalent to maximize the product of marginal likelihood $p\left({\bm{t}|\bm{\alpha},\beta}\right)$ and the priors over hyperparameters in the logarithmic case. By ignoring the irrelevant terms, we can obtain the objective equation

We exploit expectation-maximization (EM) to solve $\alpha$ and $\beta$. For $\alpha$, the update procedure is equivalent to maximize $E_{\bm{\omega}|{\mathbf{t}},\bm{\alpha},\beta}\left[{\log p\left({\bm{\omega}|\bm{\alpha}}\right)p\left(\bm{\alpha}\right)}\right]$, and then the update rule can be derived through differentiation

where ${\Sigma_{i,i}}$ represents the $i$ th diagonal element of ${\mathbf{\Sigma}}$. We maximize $E_{\bm{\omega}|\bm{t},{\mathbf{\alpha}},\beta}\left[{\log p\left({\bm{t}|\bm{\omega},\beta}\right)p\left(\beta\right)}\right]$ to update $\beta$ as

The derivation of (41) and (42) are given in AppendixB. Continue the above iterations between (18), (41) and (42) until the convergence condition is satisfied, and then we can obtain the solution of $\omega$ by $\mu$.

Furthermore, in the traditional SBL algorithm, a fix threshold $thre\_\alpha^{-1}$ is set to prune the small values of the recovered sparse signal in each iteration. The ${\mu_{i}}$ equals to zero when $\alpha_{i}^{-1}<thre\_\alpha^{-1}$, and thus the corresponding $i$th column in ${\mathbf{\Phi}}$ can be pruned out. The recovery algorithm based on the traditional SBL is unable to accurately restore the locations of all sources in the 3D SEM. Besides, it is difficult to determine the value of $thre\_\alpha^{-1}$ under different scenarios. However, the positions of the recovered signal sources are usually adjacent or close to the real signal sources in Cartesian coordinate. Therefore, we develop the clustering-based SBL (CSBL) algorithm and use adaptive threshold to improve the recovery accuracy.

Firstly, we define an adaptive dynamic threshold truncation in the pruning step of algorithm iterations, which is

Through the adaptive dynamic threshold truncation, the $\alpha_{i}^{-1}$ near the source points are selected to the most extent and the cube whose power is below threshold can be abandoned. The pruning rule is

Then, we propose to combine clustering algorithm with the SBL. Considering the number of transmitters is unknown, we employ the MMD clustering algorithm. It can adaptively determine the cluster seeds when the number of clusters is unknown and improve the efficiency of partitioning the dataset. The MMD clustering is a trial-based clustering algorithm in pattern recognition. It takes the furthest object as the clustering center based on Euclidean distance, and can avoid the situation that the cluster seeds may be too close when the initial value is selected by k-means method.

At the end of the SBL iteration process, the solution $\omega$ may only have a small number of significant coefficients. The rest of small coefficients, which contribute very little to the sparse signal, is negligible. Therefore, we first approximate $\omega$ by neglecting the small coefficients as

where $\delta$ is a negative sparsity threshold. Accordingly, the $q$ non-zero item set $\Im=\left\{{{\zeta_{1}},{\zeta_{2}},\ldots,{\zeta_{q}}}\right\}$ of $\omega$ estimated preliminarily is obtained, which represents the cubes in 3D space. Due to the characteristics of clusters they are distributed in space, the candidate set $\Im$ will be divided into $K$ clusters denoted by ${\aleph_{1}},{\aleph_{2}},\ldots,{\aleph_{K}}$ adaptively by MMD clustering. Then, by weighting the items in $\Im$ with the averaging rule, the cluster centers are obtained as

where $local_{t}^{x}$, $local_{t}^{y}$ and $local_{t}^{z}$ are the $x$, $y$, $z$ coordinates of the $t$ th cube of set ${\aleph_{k}}$ in 3D space. $local_{k}^{{x_{0}}}$, $local_{k}^{{y_{0}}}$ and $local_{k}^{{z_{0}}}$ are the $x$, $y$, $z$ coordinates of the cluster center of set ${\aleph_{k}}$, which is also the re-estimated non-zero position in the updated sparse signal ${\bm{\omega}^{*}}$. Simultaneously, we update the sparse coefficient as

Eventually, the updated sparse signal ${\bm{\omega}^{*}}$ is obtained. According to the recovered sparse signal, the SEM construction model is

where ${\mathbf{X}}$ is the spectrum vector we constructed.

In this section, the performance of proposed 3D SEM construction method is analyzed and verified by simulations under the campus scenario, as shown in Fig. 5. The ROI is $100{\text{ m}}\times 100{\text{ m}}\times 50{\text{ m}}$. We firstly discretize the ROI into $10\times 10\times 10$ cubes and each cube is $10{\text{ m}}\times 10{\text{ m}}\times 5{\text{ m}}$. Then we construct a spectrum tensor $\bm{\chi}\in{\Re^{10\times 10\times 10}}$. We denote the proposed algorithm as MMI-CMSBL. Six SEM construction algorithms, i.e., Random-SBL, Random-CSBL, Random-MSBL, MMI-SBL, MMI-CMSBL and random stagewise weak orthogonal matching pursuit (Random-SWOMP) [38] are used to construct the 3D SEM. The main simulation parameters are shown in the Table 1, where the random means that SLs are selected randomly from the ROI.

We compare the proposed MMI-CMSBL with Random-SBL, Random-CSBL, Random-MSBL, MMI-SBL, and Random-SWOMP in terms of the SSR performance. The Mean Squared Error (MSE) of sparse signal recovery is defined as

where ${\bm{\omega}^{{\text{est}}}}$ and ${\bm{\omega}^{{\text{true}}}}$ are the estimated sparse signal and true sparse signal.

It can be seen from the simulation results in the Fig. 6 that the MSE of sparse signal recovery decreases as the sampling rates increases. The proposed MMI-CSBL outperforms other algorithms in terms of the convergence speed. The performance of Random-MSBL algorithm is better than that of the Random-CSBL, MMI-SBL and Random-SBL, which reveals that the propagation model-based SBL algorithm brings more performance improvement than traditional SBL in SSR by considering scenario to construct the sparse dictionary. Besides, the SBL has better performance than other CS algorithms when recovering sparse signals if the sensing matrix has high correlation.

Fig. 7 presents the SEM construction results of Random-SBL, Random-CSBL, Random-MSBL, MMI-SBL, MMI-CMSBL and Random-SWOMP. We can also see that the sparse signal recovery accuracy of MMI-CMSBL is higher than other algorithms, and the minimal components in Fig. 7 (b) are reduced obviously through clustering and dynamic threshold operations.

The Root Mean Square Error (RMSE) is used to evaluate the RSS recovery performance. The RMSE is defined as the difference of each cube in RSS between the estimated REM and the real REM, given by

where ${P{{}_{i}^{r}}}^{{\text{est}}}$ and ${P{{}_{i}^{r}}}^{{\text{true}}}$ are the estimated and the actual RSS values at the th cube respectively. It can be seen from Fig. 8 that MMI-CMSBL achieves the best performance compared with other methods even at a very low sampling rate. SWOMP performs worse than other SBL-based algorithms, due to the highly correlation of the sparse dictionary. Furthermore, based on the sparse dictionary constructed by RT technology, the MMI-CMSBL and Random-MSBL algorithm can rapidly converge and accurately recover the spectrum map at a low sampling rate. Comparing Fig. 6 and Fig. 8, we can see that when the MSE of sparse signal is large, the SEM construction performance is unsatisfactory.

The impact of $K$ on the performance of SEM construction and sparse signal recovery are shown in Figs. 9 - 12. In Fig. 9 and 11 show that as the sparsity $K$ increases, the SEM construction error and the MSE of sparse signal recovery with MMI-CMSBL increase. Meanwhile, as shown in Fig. 9, when sampling rate is small, the SEM construction error first increases and then decreases with the increasing of $K$. Fig. 11 shows that the MSE of sparse signal recovery increases as $K$ increases, and we can successfully recover the signal at least when $M>2K\ln\left({N/K}\right)$. In addition, the convergence rate of sparse support distortion improves with the decrease of sparsity $K$.

In Fig. 10 and 12, it can be seen that with the increase of sparsity $K$ the MMI-CMSBL can still maintain excellent performance on SEM recovery and sparse signal recovery compared with other algorithms at a fixed low sampling rate. MMI-CMSBL and Random-MSBL are less influenced by the sparsity $K$ than Random-SBL. Due to the increase of RF transmitters, their interference with each other also increases. Additionally, CSBL can achieve better performance than the traditional SBL. The visualization of SEM construction with different sparsity is shown in Fig. 13.