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

TL;DR

This paper introduces a 3D spectrum environment map construction method using sparse Bayesian learning, optimizing sensor placement and dictionary construction for improved accuracy with fewer sensors in cognitive radio networks.

Contribution

The paper proposes a scenario-dependent dictionary, an MMI-based sensor layout optimization, and a MMD clustering-based SBL algorithm for efficient 3D SEM construction.

Findings

Reduces the number of required spectrum sensors.

Achieves higher accuracy at low sampling rates.

Demonstrates effectiveness through campus scenario simulations.

Abstract

The spectrum environment map (SEM), which can visualize the information of invisible electromagnetic spectrum, is vital for monitoring, management, and security of spectrum resources in cognitive radio (CR) networks. In view of a limited number of spectrum sensors and constrained sampling time, this paper presents a new three-dimensional (3D) SEM construction scheme based on sparse Bayesian learning (SBL). Firstly, we construct a scenario-dependent channel dictionary matrix by considering the propagation characteristic of the interested scenario. To improve sampling efficiency, a maximum mutual information (MMI)-based optimization algorithm is developed for the layout of sampling sensors. Then, a maximum and minimum distance (MMD) clustering-based SBL algorithm is proposed to recover the spectrum data at the unsampled positions and construct the whole 3D SEM. We finally use the simulation data of the campus scenario to construct the 3D SEMs and compare the proposed method with the state-of-the-art. The recovery performance and the impact of different sparsity on the constructed SEMs are also analyzed. Numerical results show that the proposed scheme can reduce the required spectrum sensor number and has higher accuracy under the low sampling rate.

Tables

Table 1. TABLE I: The main simulation parameters

Parameters	Value
Region of interest	$100\text{m}\times 100\text{m}\times 50\text{m}$
SEM tensor size	$10\times 10\times 10$
Granularity of SEM tensor	$10\text{m}\times 10\text{m}\times 5\text{m}$
Number of RF transmitters ($K$)	4, 8, 12, 16
Transmitting frequency ($f$)	1 GHz
Transmitting power ($P^t$)	2 W
Positions of RF transmitters	Random generated in the ROI
Sampling rate and number of RF transmitters in SEM construction performance	$r=0.1$, $K=4$
Algorithms for comparison	Random-SBL, Random-CSBL, Random-MSBL, MMI-SBL, MMI-CMSBL, Random-SWOMP

Equations

$${\omega }_n=\{\begin{array}{ll}{P}_n^t,& \text{if there is a RF transmitter in the }n\text{th cube,}\\ 0,& \text{otherwise.}\end{array}$$

$$d_{m,k}={\Vert {\mathbf{T}}_k-{\mathbf{S}}_m\Vert }_2.$$

$$P_{m,k}^r=\frac{G_t{G}_r{\lambda }^2{P}_k^t}{{\left(4\pi \right)}^2{d}_{m,k}^{\eta }},$$

$$t_m=\sum _{k=1}^K{P}_{m,k}^r.$$

$$\mathbf{X}=\mathbf{\phi }\mathbf{\omega },$$

$${\psi }_{i,j}=\{\begin{array}{ll}1,& \text{if the }i\text{th SL is at the }j\text{th cube,}\\ 0,& \text{otherwise.}\end{array}$$

$$\mathbf{t}=\mathbf{\psi }\mathbf{X}+\mathbf{\epsilon }=\mathbf{\psi }\mathbf{\phi }\mathbf{\omega }+\mathbf{\epsilon }=\mathbf{\Phi }\mathbf{\omega }+\mathbf{\epsilon },$$

$$\begin{gathered}\bm{\hat{\omega}}=\arg\min{\left\|\bm{\omega}\right\|_{0}},\hfill\\ \begin{array}[]{*{20}{c}}{{\text{s}}{\text{.t}}{\text{.}}}&{\bm{t}={\mathbf{\Phi}}\bm{\omega}+\bm{\varepsilon}}\end{array}.\hfill\\ \end{gathered}$$

$$\begin{gathered}\bm{\hat{\omega}}=\arg\min{\left\|\bm{\omega}\right\|_{1}},\hfill\\ \begin{array}[]{*{20}{c}}{{\text{s}}{\text{.t}}{\text{.}}}&{\bm{t}={\mathbf{\Phi}}\bm{\omega}+\bm{\varepsilon}}\end{array}.\hfill\\ \end{gathered}$$

$$p\left(\mathbf{t}|\mathbf{\omega },{\sigma }^2\right)={\left(2\pi {\sigma }^2\right)}^{-M/2}\exp \left\{-\frac{{\Vert \mathbf{t}-\mathbf{\Phi }\mathbf{\omega }\Vert }^2}{2{\sigma }^2}\right\}.$$

$$p\left(\beta ;c,d\right)=\text{Gamma}\left(\beta |c,d\right),$$

$$\text{Gamma}\left(\beta |c,d\right)=\Gamma (c{)}^{-1}{d}^c{\beta }^{c-1}{e}^{-d\beta },$$

$$p\left(\mathbf{\omega }|\mathbf{\alpha }\right)=\prod _{i=1}^N\mathcal{N}\left({\omega }_i|0,{\alpha }_i^{-1}\right),$$

$$p\left(\mathbf{\alpha };a,b\right)=\prod _{i=1}^N\text{Gamma}\left({\alpha }_i|a,b\right).$$

$$p\left(\mathbf{\omega }\right)=\int p\left(\mathbf{\omega }|\mathbf{\alpha }\right)p\left(\mathbf{\alpha }\right)d\mathbf{\alpha }.$$

$$p\left(\mathbf{\omega },\mathbf{\alpha },\beta |\mathbf{t}\right)\equiv p\left(\mathbf{\omega }|\mathbf{t},\mathbf{\alpha },\beta \right)p\left(\mathbf{\alpha },\beta |\mathbf{t}\right),$$

$$\begin{array}{c}p\left(\mathbf{\omega }|\mathbf{t},\mathbf{\alpha },\beta \right)=\frac{p\left(\mathbf{t}|\mathbf{\omega },\mathbf{\alpha },\beta \right)p\left(\mathbf{\omega }|\mathbf{\alpha }\right)}{p\left(\mathbf{t}|\mathbf{\alpha },\beta \right)}\\ =\mathcal{N}\left(\mathbf{\omega }|\mathbf{\mu },\mathbf{\Sigma }\right),\end{array}$$

$$\begin{array}{c}\mathbf{\mu }=\beta \mathbf{\Sigma }{\mathbf{\Phi }}^T\mathbf{t},\text{hfill}\\ \mathbf{\Sigma }={\left(\beta {\mathbf{\Phi }}^T\mathbf{\Phi }+\mathbf{\Lambda }\right)}^{-1},\text{hfill}\end{array}$$

$$d_{mn}=\sqrt{{(m_x-n_x)}^2+{(m_y-n_y)}^2+{(m_z-n_z)}^2}.$$

$$d_{h,m}^{\text{s}}=\sqrt{{(H_x-m_x)}^2+{(H_z-m_z)}^2+{(H_z-m_z)}^2},$$

$$d_{n,h}^{\text{s}}=\sqrt{{(H_x-n_x)}^2+{(H_z-n_z)}^2+{(H_z-n_z)}^2}.$$

$${\mathbf{E}}_{\text{LOS}}=E_{1\text{m}}\frac{e^{-jl{d}_{mn}}}{d_{mn}},$$

$${\mathbf{E}}_s^{\text{R}}={\mathbf{E}}_{\text{LOS}}R\frac{{\text{e}}^{-lj(d_{h,m}^{\text{s}}+d_{n,h}^{\text{s}})}}{d_{h,m}^{\text{s}}+d_{n,h}^{\text{s}}},$$

$${\mathbf{E}}_s^{\text{D}}=\frac{{\mathbf{E}}_{\text{LOS}}}{d_{n,h}^{\text{s}}}D\sqrt{\frac{d_{n,h}^{\text{s}}}{d_{h,m}^{\text{s}}\cdot (d_{h,m}^{\text{s}}+d_{n,h}^{\text{s}})}}\cdot {\text{e}}^{-jl(d_{h,m}^{\text{s}}+d_{n,h}^{\text{s}})},$$

$$\begin{array}[]{*{20}{c}}{{{\mathbf{E}}_{n}}=\sum\limits_{s=1}^{{N_{P}}}{{{\mathbf{E}}_{s}}},}&{{{\mathbf{E}}_{s}}\in\left\{{{{\mathbf{{\mathbf{E}}}}_{{\text{LOS}}}},{\mathbf{E}}_{s}^{\text{R}},{\mathbf{E}}_{s}^{\text{D}}}\right\},}\end{array}$$

$$p^n=G^n{G}^m{\left(\frac{\lambda }{4\pi }\right)}^2\left|\frac{{\mathbf{E}}_n}{{\mathbf{E}}_{1\text{m}}}\right|\text{,}$$

$${\phi }_{mn}=-10{\log }_{10}(p^n/p^m\text{)}\text{.}$$

$${\phi }_{ij}=\frac{\sum _{g=1}^{N_0}{\left(d_{ij}^g\right)}^{-p}{\phi }_g}{\sum _{g=1}^{N_0}{\left(d_{ij}^g\right)}^{-p}},$$

$$\hat{\varsigma }=\arg \underset{\varsigma \subset \mathbf{S}}{\min }H\left(\mathbf{\omega }|{\mathbf{t}}_{\varsigma }\right),$$

$$\hat{\varsigma }$$

Full text

††publicationid: pubid: 0000–0000/00$00.00 © 2021 IEEE

I Introduction

With the increase of various electronic devices, i.e., radio, radar, navigation and so on, the electromagnetic environment becomes significantly complex [1, 2, 3]. The Spectrum Environment Map (SEM) visualizes the spectrum related information, including the time, frequency and the received signal strength (RSS) of signals, and the locations of the sensor devices, on a geographical map [4, 5]. It is useful for the abnormal spectral activity detection, radiation source localization, radio frequency (RF) resource management, and so on. However, constructing an accurate SEM for the scenario with lots of buildings is still difficult. This is because the sensor device number and sampling time are limited in practice. Besides, with the development of space-air-ground integrated communication networks, electronic devices are distributed in the three-dimensional (3D) space.

Many SEM reconstruction methods have been proposed recently [6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16]. They can be divided into two major categories, i.e., the direct construction methods driven by massive measured data and the indirect construction methods driven by channel propagation characteristics.

The direct construction methods typically employ interpolation algorithms to recovery the missing data by mining the correlation between the given measured data. In [17], the authors utilized the Kriging interpolation to construct the two-dimensional (2D) SEM at different frequencies. A tensor completion method was used in [9] to recovery the missing data in both the spatial domain and temporal domain. Machine learning techniques have also been adopted for the data-driven REM construction. For example, from the perspective of image processing, the authors in [11] developed a generative adversarial network (GAN) to obtain the SEM based on the sampling data. Deep neural networks were used in [10] to ”learn” the intricate underlying structure from the given data and constructed the SEM. Nevertheless, the aforementioned methods mainly focus on the 2D SEM re-construction and can only achieve satisfactory performance by using a large amount of sampling data.

The indirect construction methods can greatly reduce the number of sampling data by using the rule of wireless signal propagation [18]. In [12, 13], the limited sampling data was used to estimate the transmitters’ information, and then the missing data was recovered by using the ideal propagation model. Considering the realistic propagation model, the authors in [19] modeled the spectrum data as the tomographic accumulation of spatial loss field.

However, in many practical cases, the number of sampling sensors is very limited and the spatial sampling rate may be much lower than the Nyquist rate. To tackle this issue, the compressed sensing (CS) technique is applied, which decomposes the measurement into the linear superposition of the sensing matrix and the sparse signal to realize data recovery under sparse sampling. For example, the authors in [20] constructed the SEM based on least absolute shrinkage and selection operator (LASSO) with random sparse sensing data. In [15], the authors proposed an improved orthogonal matching pursuit algorithm (OMP) to build compressed SEM with right-triangular (QR) pivoting based sampling locations optimization. However, the traditional CS approaches, such as LASSO [21], OMP [22] and linear programming [23], are point estimation for sparse signals. Besides, the SEM sensing matrix usually has high spatial correlation, which would greatly deteriorate the recovery performance from the noisy measurements. The sparse Bayesian learning (SBL) [24, 25, 26] can recover the exact sparse signal under the high correlated sensing matrix. A SBL-based SEM construction algorithm was proposed in [27]. The authors used a Laplacian function to describe the signal propagation model for sparse dictionary construction, which does not consider the realistic propagation scenario, especially the impact of buildings. Furthermore, the authors arranged the sampling sensors randomly without considering the sampling location (SL) optimization for different scenarios.

To fill these gaps, we propose a novel 3D SBL-based SEM construction scheme with optimized sampling positions of sensors and a scenario-dependent dictionary with full consideration of scenario characteristics. The main novelties and contributions of this paper are summarized as follows:

• •

A channel propagation dictionary design method for SEM construction under sparse sampling is proposed. Combined with the ray tracing (RT) simulation technology and the spatial interpolation algorithm, the scenario-dependent SEM construction dictionary is obtained.

• •

A maximum mutual information (MMI)-based measurement matrix optimization architecture is proposed. By modeling the construction problem as a communication channel model, the objective function is derived from the SBL framework and solved by the greedy algorithm, which improving the efficiency of spectrum data acquisition.

• •

An improved SBL algorithm based on cluster analysis is developed for 3D compressed spectrum recovery. The maximum and minimum distance (MMD) clustering and dynamic threshold pruning are combined with the SBL, which can recover sparse signals and achieve accurate SEM construction.

The rest of this paper is organized as follows. Section II gives the 3D SEM construction model and the sparse sampling model. In Section III, the details of proposed 3D SEM construction scheme is given and demonstrated. Then, Section IV presents the simulation and comparison results and Section V gives some conclusions.

II System Model

II-A 3D SEM Model

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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.

II-B Sparse Sampling Model

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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.

III 3D SEM Construction based on SBL

III-A An Overview of the Proposed 3D SEM Construction Scheme

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

[TABLE]

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

[TABLE]

with

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

with

[TABLE]

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.

III-B Scenario-dependent Sparse Dictionary Construction

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

[TABLE]

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

[TABLE]

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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.

III-C MMI-based Sampling Optimization

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

[TABLE]

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

[TABLE]

[TABLE]

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

[TABLE]

with

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

III-D SEM Construction with Improved SBL

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

IV Simulation Results And Validations

IV-A Experiment Setup

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.

IV-B Sparse Signal Recovery Performance

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

[TABLE]

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.

IV-C SEM Construction Performance

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

[TABLE]

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.

IV-D Impact of Sparsity

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.

V Conclusion

In this paper, we have investigated the issue of 3D SEM construction based on SBL, which offers high value for efficient application of CR. We have first formulated 3D SEM construction as a sparse sampling problem by exploiting the underlying sparse nature of 3D spectrum situation. Then, we have proposed a 3D scenario-dependent SEM construction scheme, which is composed of three components: sparse dictionary construction, sampling optimization and spectrum situation recovery. Firstly, considering the complexity of electromagnetic propagation environment in the actual scenario, we have designed a scenario-dependent sparse dictionary for SEM construction based on the channel model. Secondly, we have developed MMI-based sampling architecture to obtain optimized measurement matrix. Based on SBL framework, we have derived the optimization function of MMI sampling and have solved it by greedy algorithm. Finally, due to the ineffectiveness of traditional SBL recovery algorithm in 3D SEM construction, a tailored MMD-clustering based SBL algorithm has been proposed. The sparse signal recovered can achieve high precision by dynamic threshold pruning. We have also compared the sparse signal recovery performance and SEM construction performance among six methods, i.e., Random-SBL, Random-CSBL, Random-MSBL, MMI-SBL, MMI-CMSBL and Random-SWOMP. The impact of sparsity on situation recovery precision has been studied. Simulations have demonstrated the superiority of the proposed 3D SEM construction scheme.

Appendix A Proof of the equation (31)

In this part, we present the derivation of (31) in the text. Firstly, we introduce the Mahalanobis transformation lemma,

Lemma 1. An arbitrary $N$ dimension Gaussian distribution $\vartheta\sim\mathcal{N}\left({\mu,\Sigma}\right)$, we call $\rho={\Sigma^{-\frac{1}{2}}}\left[{\vartheta-\mu}\right]$ the Mahalanobis transformation, where

[TABLE]

which means that ${\rho_{i}}$ is the standard Gaussian distribution $\mathcal{N}\left({0,1}\right)$.

Proof: We have

[TABLE]

[TABLE]

we compute the Jacobian determinant

[TABLE]

According to Lemma of substitution of variables

[TABLE]

where $\vartheta\sim f\left(\vartheta\right)$ and inverse function $\vartheta=h\left(\rho\right)$. $\left|{\mathbf{J}}\right|$ is the absolute value of the Jacobian determinant. Then the distribution of $\rho$ is

[TABLE]

The lemma has been proven. Then, for continuous random variables $\omega$, we have its conditional entropy

[TABLE]

The posterior distribution of $\omega$ obeys the multidimensional Gaussian distribution $\mathcal{N}(\bm{\mu},{\mathbf{\Sigma}})$ with mean and covariance given by (18). Therefore, we use a posterior distribution to approximate the distribution of $\omega$, as SBL use a posterior mean to estimate $\omega$, so we have

[TABLE]

By using Lemma 1, we then have

[TABLE]

Accordingly, we can obtain the conditional entropy term in the formula (31), as given by

[TABLE]

The derivation of $H\left(\bm{\omega}\right)$ can be obtained by analogy.

Appendix B Proofs of the equations (41) and (42)

To obtain Eq. (41), we first define

[TABLE]

with $con$ being the item constant to $\bm{\alpha}$. Then, we let $\partial\mathcal{L}\left(\bm{\alpha}\right)/\partial\alpha_{i}=0$ to find the stationary point of $\alpha_{i}$, as

[TABLE]

And we obtain

[TABLE]

To obtain Eq. (42), we define

[TABLE]

Then, we let $\partial\mathcal{L}\left(\beta\right)/\partial\beta=0$ to find the stationary point of $\beta$, and we obtain

[TABLE]

where

[TABLE]