Enhancing the Accuracy and Robustness of a Compressive Sensing Based Device-Free Localization by Exploiting Channel Diversity

As an emerging and promising technique, device-free localization (DFL) estimates target positions by analyzing their shadowing effects. Most existing compressive sensing (CS)-based DFL methods use the changes of received signal strength (RSS) to approximate the shadowing effects. However, in changing environments, RSS readings are vulnerable to environmental dynamics. The deviation between runtime RSS variations and the data in a fixed dictionary can significantly deteriorate the performance of DFL. In this paper, we introduce ComDec, a novel CS-based DFL method using channel state information (CSI) to enhance localization accuracy and robustness. To exploit the channel diversity of CSI measurements, the DFL problem is formulated as a joint sparse recovery problem that recovers multiple sparse vectors with common support. To solve this problem, we develop a joint sparse recovery algorithm under the variational Bayesian inference framework. In this algorithm, dictionaries are parameterized based on the saddle surface model. To adapt to the environmental changes and different channel characteristics, dictionary parameters are modelled as tunable parameters. Simulation results verified the superior performance of ComDec as compared with other state-of-the-art CS-based DFL methods.


Introduction
As a main piece of the context information, location information is essential for the location-based services (LBS) (all the used abbreviations are explained at the end of the article) in many context-aware applications. Nowadays, localization techniques have been an active field in pervasive and ubiquitous computing. Generally, the target localization methods can be classified as: device-free [1][2][3] and device-based [4][5][6]. Device-based methods need objects to attach assistant wireless devices for signal transmitting or receiving. However, in many applications, such as emergency rescue, intruder detection, and smart homes, it is difficult to attach objects with any transceivers. In this case, the device-based method will be infeasible.
Device-free localization (DFL) methods do not need to attach objects with any assistant devices. This approach has become a crucial component in many context-aware applications. As the transceiver-free objects cannot be directly perceived, DFL methods estimate their positions by analyzing their influences on surrounding radio environments. Among current DFL methods, some rely on specialized hardware, such as radar-based methods [7], camera-based methods [8], and infrared-based methods [9]. These methods need an extensive deployment of dedicated devices in the monitoring area and may involve privacy issues. In recent years, DFL methods based on existing infrastructures (e.g., WSNs/WiFi) have attracted a lot of research interests. They do not require dedicated hardware and only make use of the measurement information available from the already deployed wireless devices.
As a common link measurement, the received signal strength (RSS) is nearly ubiquitously available from the standard radio transceivers [10,11]. The RSS measurement information has been widely used in DFL. When targets located at different positions, the RSS readings will be different. To exploit the location dependence of RSS, RSS readings on multiple wireless links are recorded before and after targets entering into the monitoring area. To establish the relationship between RSS variations and target locations, the fingerprinting technique has been introduced [12]. Typically, the fingerprinting-based DFL consists of an online and offline phase. In the offline training phase, by gathering the RSS variations caused by a target at every possible location, a radio map can be built. In the online locating phase, target locations are estimated by matching current RSS variations with the fingerprints in radio map. However, as a major drawback of the fingerprinting-based DFL methods, gathering fingerprints to build a radio map is labor-intensive and time-consuming.
As an alternative to fingerprinting-based methods, model-based methods have been widely used in DFL. They map RSS variations to target positions by theoretical or empirical shadowing models. Based on these models, the dictionary (a.k.a. radio map) can be built without site survey. Unfortunately, in changing environments, RSS is extremely sensitive to environmental dynamics, such as temperature, humidity, electromagnetic characteristics, and pedestrians around [13]. In this case, shadowing effects cannot be well approximated by RSS variations. The spatial and temporal environmental dynamics may result in mismatches between runtime RSS measurements and the data in a fixed dictionary. In [14], a dictionary refinement algorithm is introduced, which alleviates the dictionary mismatches by iteratively refining the model-based dictionary. However, the real-time optimization of dictionary parameters leads to a high computational complexity. To avoid this, in ComDec, the dictionary parameters are modelled as tunable parameters to adapt to the changes of environment, and we do not need to explicitly estimate their values.
Recent years, in target localization, the channel state information (CSI) [15] has been exploited, which is a measurement information from PHY layer. As a fine-grained value, CSI depicts the channel quality on multiple orthogonal subcarriers. In wireless environments, due to the frequency-independent attenuation and frequency-selective fading, each channel will exhibit a unique amplitude and phase. Hence, the dictionaries corresponding to different channels are different. Over the past few years, in many device-based localization methods, CSI measurements has been leveraged [16]. Recently, more efforts have been paid in device-free technique. In [17], the fine-grained subcarrier information is exploited based on the multitask Bayesian compressive sensing (MBCS). However, the DFL method proposed in [17] do not provide a solution to the dictionary mismatch problem caused by environmental dynamics.
Compared with RSS measurements, CSI measurements are more robust and suitable for being utilized in DFL. In this paper, a novel ComDec method is proposed. It can enhance the accuracy and robustness of the CS-based DFL in changing environments. In ComDec, to enrich the measurement information, CSI measurements are collected from multiple frequency bands. Moreover, to reduce wireless sensor nodes, the compressive sensing (CS) theory [18] is applied in ComDec by taking advantage of the spatial sparsity of target localization. The CS-based DFL problem in multi-channel scenario is formulated as a joint sparse recovery problem. Furthermore, to bypass the dictionary training and retraining works, the dictionaries are built based on the saddle surface model [19]. We treat dictionary parameters as tunable parameters to adapt the changes of environment and different channel characteristics. Afterwards, to recover the sparse vectors of multiple channels, we develop a joint sparse recovery algorithm under the variational Bayesian inference framework [20]. The main contributions of this work are as follows: • To enhance the localization accuracy and robustness of CS-based DFL, a novel ComDec method is proposed, which leverages the channel diversity of CSI measurements. In ComDec, the CS-based DFL problem is extended to multi-channel scenario. It is formulated as a joint sparse recovery problem that recovers multiple jointly sparse vectors over two known dictionaries.

•
To simultaneously recover the jointly sparse vectors, we develop a novel joint sparse recovery algorithm. The joint sparsity of the sparse vectors is induced by a novel two-layer hierarchical prior model. Then, the common support set of the sparse vectors is estimated by inferring the posteriors of the hidden variables that defined in the proposed prior model.

•
To mitigate the influence of environmental dynamics in changing environments, the dictionary parameters with respect to multiple channels are modelled as tunable parameters to adapt the environmental changes and different channel characteristics. In this way, the dictionary mismatch problem can be solved without the need of explicitly estimating the dictionary parameters.

•
To reduce the computational complexity, we introduce four methods in the proposed joint sparse recovery algorithm. Among them, the grid pruning method can improve the convergence speed of the proposed joint sparse recovery algorithm.
The remainder of this paper is organized as follows. Section 2 presents an overview of the related works on multi-target DFL. Section 3 gives the signal model and formulates the CS-based DFL problem as a joint sparse recovery problem. Section 4 proposes a novel joint sparse recovery algorithm to recover the jointly sparse vectors. Section 5 validates the proposed ComDec method with extensive numerical simulations. Finally, Section 6 concludes this paper.

Related Work
Beginning with the initial papers of Youssef et al. [1] and Zhang et al. [2], numerous research works on DFL have been carried out [17,21]. Typically, there are four types of DFL methods: (1) geometry-based methods, (2) fingerprinting-based methods, (3) radio tomographic imaging (RTI)-based methods, and (4) CS-based methods. The geometry-based methods estimate target locations based on the geometry information of wireless links. They have a restriction on target spacing and require to know the deployment information of wireless nodes. Fingerprinting-based methods estimate target locations by matching runtime measurements with fingerprints in a radio map. They can achieve enhanced accuracy, but a site survey is needed for radio map building and retraining. RTI-based methods do not require offline training effort. They treat target locations as the attenuation images of distorted links and can achieve an improved performance. However, they need a dense deployment of wireless links. It may lead to high hardware cost and great energy consumption. CS-based methods formulate the DFL problem as a sparse signal reconstruction problem. Compared with the aforementioned DFL methods, the CS-based DFL methods can ensure a high accuracy with fewer measurements. LCS [22] is a CS-based DFL method, in which the model-based dictionary satisfies the restricted isometry property (RIP) with high probability. E-HIPA [23] is a representative CS-based DFL method. It adopts an adaptive orthogonal matching pursuit (OMP) algorithm to estimate the target number and location vector. DR-DFL [14] introduced a real-time dictionary refinement algorithm for CS-based DFL. It can mitigate the influence of environmental dynamics in changing environments. The dictionary refinement is realized by optimizing the environment-related dictionary parameters. However, explicitly estimating the dictionary parameters may lead to a high computational complexity.
To leverage CSI for DFL, Pilot [24] regards the correlations of CSI as fingerprints and uses the maximum a posteriori probability (MAP) estimator to estimate target location. MonoPHY [25] is a fingerprinting-based DFL method that leverages CSI to locate a target with only one stream. Gao et al. [26] transformed CSI measurements into a radio image and adopted the machine learning method to estimate the position of a person. The aforementioned DFL method leverage CSI by fingerprinting technique. They seek appropriate location-dependent CSI features for DFL and trying to build a robust and precise relationship between CSI measurements and target position. However, to make full use of the channel diversity of CSI measurements is still a challenging issue in CS-based DFL.

Overview of Multi-Target Device-Free Localization
The proposed ComDec method is a CS-based DFL method that aims to estimate the number and locations of multiple transceiver-free objects in changing environments. Figure 1 shows an example of the CS-based DFL. As can be seen, multiple transceiver-free objects are randomly distributed in an l × l two-dimensional (2D) monitoring area A. Owing to the inherent spatial sparsity of the targets in A, the position information of multiple targets can be considered as a sparse signal. To achieve this, we discretize A into N equal-sized grids. When a target is located in grid n, we regard the grid center as its position. In this case, the target number and position information can be encoded in a sparse location vector, i.e., where θ ∈ R N×1 is a K-sparse vector. θ n ∈ {0, 1} is its n-th component, which indicates whether a target is located in grid n. If there is a target, θ n = 1; otherwise, θ n = 0. Thus, the location information of K targets can be denoted as L = {(x n , y n ) |θ n = 1, n ∈ {1, ..., N} }, where (x n , y n ) is the coordinate of grid n. Moreover, the target number can be denoted as K = θ 0 . When the target number K is far less than the grid number N (K N), the sparsity of θ can be ensured. As a matter of fact, the grid number N is usually a big number, while there is only a small number of targets in the monitoring area.  As seen from Figure 1, A is surrounded by several wireless nodes. Between these nodes, M bidirectional wireless links are established to cover A. They travel through A and sense the target-induced shadowing effects in the electromagnetic field. When there exists a target, some of the wireless links will be shadowed, and the signal power or other features of radio signals will be affected due to the scattering, reflection, refraction, and absorption of the signals. Fortunately, the changes of the signal features are closely related to target positions. Therefore, by analyzing the target-induced shadowing effects, the location information of multiple targets can be inferred.

CSI Collection and Feature Extraction
RSS is a coarse-grained measurement information from MAC layer. It represents the overall signal power across all subchannels. Due to the multipath fading, RSS is unreliable and varies with time even in a static environment. This may result in limited localization accuracy, especially in changing environments. To enhance the localization performance of CS-based DFL, ComDec adopts CSI to characterize target-induced perturbations. CSI is a fine-grained measurement information from PHY layer. It can provide the phase and strength information of the signals on different subchannels. For link m, the CSI on each subchannel is a complex value which is defined as where H f m denotes the CSI measurement corresponding to the f -th channel of link m. |H f m | and ∠H m denote the amplitude and phase response, respectively. F is the number of channels. The amplitude response |H f m | is the change of amplitude of link m on channel f . By converting it from linear space to logarithmic space, the corresponding power fading can be written asH f m = 20 log |H f m |/10 3 (dBm) [27]. We collect the CSI measurements from F channels, and a set of power fading information with channel diversity can be obtained. They can provide redundancy information to alleviate the location ambiguity that is incurred by environmental dynamics.

Problem Formulation
To characterize shadowing effects, we measure the the change of power fading on each link. For each channel, the change of power fading on link m is The changes of power fading are measured on M links at runtime. For each link, we can obtain F measurements from F channels. Thus, the measurement information can be represented as where Y ∈ R M×F is the measurement matrix.
If any two targets are located sparsely [23] in the monitoring area, S f m can be represented as the summation of attenuations that occur in each cell [28]. Therefore, ∆H f m can be expressed as whereh f m,n is the shadowing loss that caused by a target located in grid n. For M links, the measurement vector y f can be written as where Φ f ∈ R M×N is a dictionary. Its (m, n)-th element ish f m,n . f ∈ R M×1 represents the noise vector of channel f , and its m-th component is f m . In ComDec, the DFL problem can be viewed as a problem of reconstructing the location vector θ from the measurements Y. Theoretically, this problem can be solved by existing joint sparse recovery methods. Nonetheless, as a common shortcoming, they require the true dictionaries {Φ f } F f =1 to be known in advance. However, in changing environments, it is impossible for us to accurately estimate these dictionaries.
Due to the difference in channel characteristics, the dictionaries with respect to different channels are different. Furthermore, in changing environments, the dictionary for an individual channel may differ when observed at different times. In this case, the fingerprinting method will be infeasible because it may cumulate the effects of dictionary mismatches on multiple channels and deteriorate the localization performance significantly. As an alternative, we adopt the saddle surface model to characterize the shadowing effect and establish multiple model-based dictionaries for multiple channels. To adapt to the changes of environment and different channel characteristics, the environmental parameters in these dictionaries are considered as adjustable. For simplicity, we denote φ f m,n =h f m,n , and thus Φ f can be expressed as Figure 2 depicts the spatial impact area of a wireless link. As can be seen, the spatial impact area is an ellipse area, and we set up an U-V coordinate system in the area. Based on the saddle surface model, φ f m,n can be expressed as where (U m,n , V m,n ) denotes the coordinate of grid n. λ 1 and λ 2 denote the semi-major and semi-minor axes of the spatial impact area, respectively. In the saddle surface model, γ f m represents the maximum shadowing loss, and ρ f m ∈ (0, 1] denotes the shadowing ratio, which represents the normalized shadowing loss at the midpoint of LOS path.
Wireless node Spatial impact area Transceiver-free object Let ω f = γ f · θ and υ f = γ f ρ f · θ be the unknown sparse vectors. Based on (6) and (9), y f can be rewritten as where Ψ ∈ R M×N and Ψ ∈ R M×N are known dictionaries. Their (m, n)-th elements are computed as ψ m,n = U 2 m,n λ 2 1 and ψ m,n = 1 − U 2 m,n λ 2 1 − V 2 m,n λ 2 2 , respectively. We can determine their values before the localization stage according to the deployment of wireless links. It is noteworthy that the location vector θ and sparse vectors S = {ω f , υ f } F f =1 share a common sparse support set T ⊆ {1, 2, ..., N}. In other words, their nonzero entries are concentrated at some common locations.
As pointed out earlier, the values of {γ f , ρ f } F f =1 are closely related to the environmental characteristics, which are different for different channels and times. Let γ=[γ 1 , γ 2 , ..., γ F ] and ρ=[ρ 1 , ρ 2 , ..., ρ F ]. We define two matrices Θ, Θ ∈ R N×F , which are constructed as follows: Both Θ and Θ are K jointly sparse matrices. This means there are at most K rows in them that have nonzero elements. The column vectors of Θ and Θ share the common support T, which is the index set of the grids where a target exists. The cardinality of T is K := |suppΘ| = |suppΘ |, which also denotes the target number. Based on (10), (11) and (12), the measurement matrix Y can be expressed as where E ∈ R M×F is the matrix of measurement noises, whose f -th column is f . Θ and Θ are the unknown matrices. By (13), the CS-based DFL problem is recast as a problem that needs to recover Θ and Θ simultaneously over two known dictionaries. It should be pointed out that the problem is different from the conventional joint sparse recovery problem, which only needs to reconstruct a single sparse matrix over a single sparsifying dictionary. However, in this problem, there are two different sparse matrices to be estimated. Therefore, existing joint sparse recovery algorithms cannot be directly applied in ComDec. In this context, the key issue of ComDec is to design a joint sparse recovery algorithm to estimate Θ and Θ simultaneously. It should be noted that, the channel diversity of CSI measurements can be exploited to enhance the accuracy and robustness of the joint sparse recovery. Figure 3 shows an architectural overview of ComDec. The proposed ComDec consists of four main modules: dictionary construction, measurement information collection, joint sparse recovery, and location estimation. The process of target counting and localization is illustrated as follows: First, in the dictionary construction module, we establish two constant dictionaries Ψ and Ψ according to the saddle surface model. Then, at runtime, CSI measurements are collected from M links. Each stream contains the CSI readings of F channels. Afterwards, in the joint sparse recovery module, the posteriors of all hidden variables are inferred by a joint sparse recovery algorithm. Finally, with the knowledge of the posteriors, in the location estimation module, the common support set T can be estimated, and thus the estimated Cartesian coordinates of multiple targets can be obtained.

Target Localization via Variational Bayesian Inference
In this section, we develop a novel joint sparse recovery algorithm. First, to induce the joint sparsity of the jointly sparse vectors S, a two-layer hierarchical prior model is introduced. Then, by using the variational Bayesian inference technique, we infer the posteriors of the hidden variables that defined in the proposed hierarchical prior model. Finally, based on the posteriors of S, target counting and localization are implemented.

Hierarchical Prior Model
The joint sparsity of S is induced by a non-separable sparsity inducing prior model [29]. The graphical model for the joint sparse recovery is shown in Figure 4, which describes the dependencies between random variables. In the first layer, we regard Moreover, the Gaussian-inverse-Gamma prior is imposed on each sparse vector to encourage its sparsity. ω f is treated as a Gaussian random variable, whose prior distribution is can be given as where N · 0, α −1 n denotes the Gaussian distribution with zero mean and variance of where ρ f is regarded as an unknown deterministic parameter, a multivariate Gaussian prior is also imposed on υ f . The variance of υ f n can be given as (ρ f ) 2 α −1 n . We set ρ f = 1 to accommodate the worst case of the variance. In this situation, the prior of υ f can be expressed as The variances of υ   To allow conjugate-exponential analysis [20], an independent identically distributed (i.i.d.) Gamma distribution is imposed on α, which can be expressed as where Gamma (· |c, d n ) denotes the Gamma distribution with parameters c and d n , We set the hyperparameters c and d n to very small values (e.g., 10 −6 ) to provide a non-informative hyperprior over α n . The Gamma distribution is generally chosen as the prior for the inverse variance of a Gaussian distribution, because it is the conjugate prior of the Gaussian distribution. In this case, the associated Bayesian inference can be performed in closed form [4].
When a hierarchical Gaussian prior model imposed on ω f , the true prior distribution of ω f can be computed by marginalizing the parameter α, i.e., In this case, the true distribution of ω f is a Student-t pdf, with mean κ = 0, parameter λ = c d and degrees of freedom v = 2c. According to the property of the Student-t distribution, when v is small, this distribution will exhibit very heavy tails. Thus, it favours sparse solutions, which include only few of the basis functions and prunes the remaining basis functions by setting the corresponding weights to very small values. In this case, a sparse vector ω f can be induced with the hierarchical Gaussian prior model. As the case of ω f , we can also induce the sparsity of υ f by the proposed hierarchical Gaussian prior model. As illustrated earlier, f follows an i.i.d. Gaussian distribution. The prior distribution of f is defined as where β f is the inverse variance of f . We treat β f as a hidden variable and define β = {β f } F f =1 . As f follows a Gaussian distribution, a Gamma prior is also imposed on each β f , i.e., where a f and b f are deterministic parameters of the Gamma distribution. We denote To assume uninformative priors for β, the hyperparameters a and b are also set to very small values (e.g., 10 −6 ).
In the proposed graphical model, the observed variables are y = {y f } F f =1 , and the hidden variables are z {α, ω, υ, β}. To estimate the jointly sparse vectors, we need to infer the posterior distributions of z based on the predefined prior evidence and the measurement data. In addition, the deterministic parameters of the prior model are Ω {c, d, a, b}, which are fixed at small values to allow uninformative hyperpriors for α and β.

Variational Bayesian Inference
In the joint sparse recovery module, the key task is to infer the posteriors of z. Afterwards, based on these posteriors, the target number and locations can be estimated. For this objective, the variational Bayesian inference technique is adopted, which is applied due to it can deal with complicated Bayesian models [20]. Based on (10) and (19), the likelihood function of channel f can be written as According to the chain rule of probability, the joint probability density function (PDF) of y and z can be written as For an arbitrary density function q (z), the evidence p (y) = p (y, z)dz can be decomposed as where F (q; Ω) is a lower bound of ln p (y), and KL (q p ) represents the Kullback-Leibler divergence (KLD) between the approximated posterior q (z) and the exact posterior p (z |y; Ω ). The proposed joint sparse recovery algorithm maximizes ln p (y) iteratively. At each iteration, we set KL (q p ) = 0 to minimize the KLD and update q (z) accordingly. This will lead to the lower bound F (q, Ω) increasing to ln p (y). Meanwhile, the updating of q(z) may enlarge ln p (y) and lead to a new non-negative KLD. We will minimize the new KLD and update q(z) in the next iteration. By doing so, the log-likelihood ln p (y) will be maximized, and the approximated posterior q (z) can be optimized iteratively.
However, p (z |y; Ω ) cannot be computed analytically. Thereby, directly updating q (z) is intractable. To bypass this difficulty, we resort to the variational approximation method. It assumes that the posteriors of z are independent, i.e., By applying the above assumption, in each iteration, the log-posterior of z i ∈ z can be approximated as the expectation of the joint PDF with respect to other hidden variables (z − z i ). More specifically, the log-posteriors of z are approximated as ln q(α) = ln p (y, z) q(ω)q(υ)q(β) + ξ, where ξ denotes a normalizing constant. It makes the corresponding q (·) a true PDF. The update rule for each hidden variable is derived below. In (27), the terms independent of ω f can be treated as a constant value. Keeping only the terms that are related to ω f , ln q(ω f ) can be given as Note that q(ω f ) and q(υ f ) follow the multivariate Gaussian distribution. We assume they have the following forms after the updating where µ f ω and µ f υ are the mean vectors, Σ f ω and Σ f υ are the covariance matrices. The update of the posterior distribution is equivalent to seeking appropriate values for the parameters in the approximated posterior distribution. Our goal is to learn the values of the mean vectors and covariance matrices based on the prior distributions and likelihood function. Substituting (14) and (21) into (31), after some rearrangement, yields where In the same manner, ln q(υ f ) can be given as Substituting (15) and (21) into (37), after some rearrangement, yields where Keeping only the terms of (29) that are related to β f , ln q(β f ) can be given as We assume that the posterior of β f follows a Gamma distribution, i.e., whereã f andb f denote the deterministic parameters of the updated posterior distribution. To infer them, we substitute (20) and (21) into (41). After some rearrangement, the posterior can be given as whereã Similarly, we only keep the terms that are related to α in (30), and thereby ln q (α) can be given as The approximated posterior of α is assumed to be a multivariate Gamma distribution, i.e., wherec and {d n } N n=1 are deterministic parameters. Substituting (14), (15) and (16) into (46), after some rearrangement, yields The notation [·] n denotes the n-th component of the input vector, and [·] n,n denotes the (n, n)-th entry of the input matrix. Based on the results of posterior inference, the required expectations can be calculated as where α n =c d n , n = 1, 2, ..., N.

Joint Sparse Reconstruction
According to the above update rules, we can successively update the posteriors of hidden variables z. The jointly sparse vectors S can be reconstructed according to these posteriors. The algorithm of reconstructing S is summarized as follows: 1.

4.
If a convergence criterion has been met, terminate and choose the posterior means of ω and υ as the estimation of S. Otherwise, go to step 1).
Based on the above joint sparse recovery algorithm, we can obtain a reconstructed sparse The computational cost of the algorithm is dominated by the matrix inversion operations in (35) and (39), whose computational complexities are O N 3 F . Moreover, the matrix-vector multiplications in (36) and (40), as well as the matrix multiplication operations in (45) can also incur heavy computational burden. Their computational complexities are O N 2 MF . Thus, when applying ComDec in large-scale areas (where N is large), the joint sparse recovery algorithm will be computationally expensive. In this paper, we adopt the following four methods to alleviate the computational burden of the proposed algorithm.
In step 1, the covariance matrices Σ f ω and Σ f υ are computed by (35) and (39), which contain matrix inversion operations. Using the matrix inversion lemma [30], the covariance matrices can be evaluated as where With these matrix transformations, we only need to compute the inversions of Ξ f ω and Ξ f υ , whose computational complexities are O M 3 (M N). Note that, Λ is an N × N diagonal matrix. We can easily obtain its inversion.
As mentioned before, the expectation evaluations of µ f ω and µ f υ can also lead to a high computational cost. To reduce their computational complexities, we reformulate (36) and (40) to cast them as a problem of solving the following linear systems of equations: The equations can be solved by the conjugate gradient (CG) algorithm. Theoretically, we can reach the exact solution after N iterations. However, this may incur a considerable computational burden. Fortunately, in practice, a few iterations is sufficient to obtain a good solution. If W (W N) iterations are required, the computational complexity of the CG algorithm will be O (W N log N).
In step 2, the computational cost is mainly attributed to the matrix multiplications in the last term of (45). To mitigate this, (45) can be rewritten as Here, matrix multiplications are replaced by simple element-wise operations. As a consequence, the computational complexity of step 2 can be reduced to O (NMF). By using the CG method and the above reformulations, the computational cost per iteration of the proposed joint sparse recovery algorithm can be reduced to O N M 2 F . To further reduce the computational load and speed up the convergence, we conduct real-time grid pruning according to the posterior of α. In each iteration, when α n is sufficiently large, {ω will be driven to zero. This imply that the contribution of grid n to the signal power fading is negligible. In this case, we can remove grid n from the grid set Π, which is defined as the set of grids where a target possibly exists. With the reduction of grids in Π, the computational load of the next iteration will decrease.
We denote the initial grid set by Π (0) {1, 2, ..., N}. Then, the grid set in each iteration can be updated as where τ denotes the iteration number and α th is the threshold of α n . For the reconstruction of S, the selection of α th provides a trade-off between the localization accuracy and the convergence speed. After grid pruning, the dimension of θ can be reduced to N = |Π (τ) |. In addition, α, S, Ψ, and Ψ also should be pruned accordingly. By virtue of the local convergence property of variational Bayesian inference [20], the proposed joint sparse recovery algorithm is guaranteed to be convergent. The stop criterion is set as the residual error Res smaller than a pre-determined threshold r th , where the Res is defined as (62) Furthermore, to prevent the computational load from being excessively high, we set a maximum iteration number τ max . In this situation, the iterative algorithm will stop when Res becomes smaller than r th or τ reaches τ max .

Target Counting and Localization
Based onŜ, the target number and locations can be estimated. However, these vectors are not strictly sparse. They may contain many negligible but nonzero coefficients. Thus, we estimate their common support by the following formulae: The common support of the jointly sparse vectors is estimated based on µf ω , in whichf is the channel that has the minimum residual error. In (63), µ th is the sparsity threshold. We use it to filter out the small coefficients in µf ω . After that, we can estimate the target number and locations asK = |T| and L = {(x n , y n )|n ∈T}, respectively.

Simulation Setup
The numerical simulations are carried out in MATLAB R2015b 64 bit version running on a PC with i7-8550U and 16 GB memory. In our simulations, the channel state information is assumed to be collected from the wireless devices that implement an orthogonal frequency division multiplexing (OFDM) system. The wireless devices are assumed to mode at 2.4 GHz with 20 MHz bandwidth. Table 1 summarizes the default values of simulation parameters. The numerical simulations are conducted in a complex radio transmission environment. In this transmission environment, as the impact of environmental dynamics, the environmental parameters are changing with time. Moreover, the measurement of each channel is corrupted by addictive white Gaussian noise. The signal-to-noise ratio (SNR) of each channel is defined as where σ f is the standard deviation of the measurement noise vector. In our simulations, all results are averaged over T = 10 3 Monte Carlo trials. For each trial, the localization error E t is computed as where L k andL k represent the true and estimated Cartesian coordinates of target k. We use the "average localization error" (Avg.Error) and "root-mean-square error" (Rms.Error) as the metrics to measure the localization accuracy. The Avg.Error is defined as The Rms.Error is defined as To evaluate the target counting performance, we introduce another metric Prob.CoC. It is the probability of correct counting (i.e.,K = K). We compare the performances of ComDec with other three CS-based DFL approaches: (i) LCS with the GMP algorithm (LCS-GMP) as the CS recovery algorithm [22], (ii) E-HIPA with the adaptive OMP algorithm (E-HIPA-OMP) as the CS recovery algorithm [23], and (iii) DR-DFL with the VEM algorithm for dictionary refinement and sparse recovery (DR-DFL-VEM) [14]. Table 2 reports the computational complexities and localization accuracies of multiple DFL methods. As can be seen from this table, the ComDec method can achieve the lowest average localization error (Avg.Error). It can enhance the accuracy and robustness of CS-based DFL by exploiting channel diversity. The computational complexity of the ComDec method is higher than the LCS and E-HIPA methods. The inferiority is mainly attributed to the estimation of the posterior distributions of sparse vectors. It noteworthy that the E-HIPA and LCS methods reconstruct the unknown sparse vector by using the OMP and GMP algorithms, respectively. They are greedy algorithms for sparse recovery, and can only provide a point estimation of the unknown sparse signal. In contrast, the proposed sparse recovery algorithm and VEM algorithm reconstruct sparse signal from a Bayesian perspective. They can provide a posterior belief (distribution function) for the values of the sparse signal, and therefore can achieve an improved accuracy. In real deployment, the sparse recovery algorithm is usually conducted on a server, where the computational resources are abundant. The server collects CSI measurements from all wireless nodes and estimates the target location accordingly. The simulation flowchart is shown in Figure 6. In the upper part of the flowchart, the measurement data is generated. In the lower part of the flowchart, the target locations can be estimated.

Impact of the Number of Channels
In this subsection, the impact of channel number on the performance of ComDec is evaluated. Theoretically, the accuracy of joint sparse recovery is closely related to the channel number. If we increase F, more useful information will be provided. Consequently, the target counting and localization performance will be improved.
In Figure 7, F varies from 1 to 25, and the Avg.Error of ComDec decreases as F increases. It is noteworthy that, the Avg.Error decreases dramatically when F < 10, while decreases slowly when F > 10. As in the case of Avg.Error, the Prob.CoC also increases as F increases, and its growth rate drops gradually. This is because the increase of F leads to a more complex hierarchical prior model, which contains (3F + 1) hidden variables. When reconstructing S, the posterior distributions of these hidden variables should be inferred. Thus, a large F may deteriorate the accuracy of the posterior inference. From this perspective, increasing the channel number has a negative effect on sparse recovery. In fact, when F > 20, the negative effect of increasing F will offset the positive effect contributed by channel diversity. Hence, we set F = 15 in our simulation to achieve a trade-off between model simplicity and the performance of target counting and localization.

Effectiveness of ComDec in Changing Environments
To mitigate the dictionary mismatches introduced by environmental dynamics, the proposed ComDec method combines the environmental parameters and the location vector to form a set of jointly sparse vectors. We treat these vectors as random variables and learn their values by variational Bayesian inference. This makes the dictionary parameters can adapt to the environmental changes and different channel characteristics. In this subsection, we test the effectiveness and robustness of ComDec in changing environments.
Firstly, we set up a simulated changing environment. To simulate the environmental dynamics in changing environments, we add Gaussian noises to the environmental parameters (γ f and ρ f ). Based on these noisy environmental parameters and the saddle surface model, the dictionaries {Φ f } F f =1 can be built. After that, we can obtain the measurement vectors {y f } F f =1 according to (6). In our simulation, the values of γ f and ρ f are calculated by where i γ and i ρ are additive white Gaussian noises whose variances are 0.6 and 0.06, respectively. The setting of the noise variances is according to the real data in changing environments. t represents the times of environmental changes and also denotes the number of additive noises added to the original environmental parameters (γ  Figure 8 depicts the variations of two dictionary atoms in the simulated changing environment. As expected, the accumulation of additive noises can lead to a deviation of dictionary atoms, which is a simulation of the effect of environmental dynamics in real-world environments. Time t=10 Figure 8. The values of the dictionary atoms: (a) φ. 5 10 ; (b) φ 12 42 , when t = 1, 5, and 10.
Secondly, in the simulated changing environment, the localization performances of multiple CS-based multi-target DFL methods are evaluated. Figure 9 compares the localization performances when t increases from 1 to 10. The proposed ComDec method that using the variational Bayesian inference technique for joint sparse recovery (ComDec-VBI) can achieve the minimum Avg.Error. It should be noted that, its localization performance is robust to the environmental changes. In contrast, with the increase of t, the localization accuracies of other three DFL methods deteriorate gradually. This reveals that these methods are sensitive to environmental dynamics. The simulation results demonstrate the impact of environmental dynamics on Avg.Error and verify the accuracy and robustness of the ComDec method. Finally, we investigate the sparse reconstruction performances of the CS recovery algorithms that adopted in the CS-based DFL methods. We run E-HIPA, LCS, DR-DFL, and the ComDec method with F = 5 and F = 15 to reconstruct the target location vector. Figure 10 shows an example of the recovered sparse vectors of different DFL approaches. As can be seen, the recovered sparse vectors corresponding to E-HIPA and DR-DFL have many negligible but nonzero coefficients. Although the reconstructed sparse vector in Figure 10c does not have small coefficients, the indices of the nonzero components are different from those of the original location vector. Fortunately, the recovered sparse vector corresponding to our proposed ComDec method has a few small coefficients, and the indices of their significant coefficients are the same as those in the original location vector. Based on them, the target number and target locations can be estimated correctly.  Figure 11 shows the target positions estimated by multiple CS-based DFL methods. As can be seen, the ComDec method can accurately estimate target positions, whereas other DFL methods have 1 or 2 incorrect estimated positions. This is because the proposed ComDec method can eliminate the position ambiguities introduced by environmental dynamics by taking advantage of the channel diversity.

Localization Error vs. SNR
In this subsection, the localization performances of different DFL methods under different SNRs are evaluated. In cellular communication (e.g., LTE), the typical SNR value is below 0 dB. However, in short range communication protocols such as 802.11 b/g, the typical SNR value is from 10-15 dB (low signal) to 40 dB (excellent signal). In fact, the scenario that is chosen for our simulations is the short-range communication scenario. Thus, the SNR in this simulation is varied from 0 dB to 40 dB. Figure 12 shows the Avg.Error of the ComDec method (F = 10, 15, and 20) and other DFL methods. As can be seen, the Avg.Error of all DFL methods decrease with the increase of SNR. Generally, the localization error is mainly caused by the measurement noise and the environmental dynamics in changing environments. As the ComDec and DR-DFL have countermeasures to mitigate the dictionary mismatches introduced by environmental dynamics, their Avg.Error decrease dramatically when SNR increases from 0 dB to 25 dB. Furthermore, as can be seen from Figure 12, by increasing the channel number, the accuracy improvement of ComDec in low SNR case (SNR < 20 dB) is much larger than the improvement in high SNR case (SNR > 20 dB). This demonstrates that the channel diversity can help to reduce the location uncertainty that introduced by measurement noises especially in the low SNR case. Moreover, the DR-DFL performs better than LCS in the high SNR case, while performs worse in low SNR cases. This is because, in the low SNR case, the measurement noise can significantly affect the dictionary refinement process. It may lead to a wrong estimation of the dictionary parameters, and the dictionary mismatches will deteriorate the localization performance seriously. Fortunately, in the high SNR case, the DR-DFL can correctly estimate the dictionary parameters and performs better than LCS.

Localization Error vs. Number of Targets
In this subsection, we examine the performances of multiple CS-based DFL methods when K varies from 2 to 20 at a step of 2. In this simulation, we set the environmental change time t = 1 to simulate the environmental dynamics in changing environments. With the increase of K, the joint sparsity level of the jointly sparse vectors S changes accordingly. This will make the localization accuracies of all CS-based DFL decrease. Figure 13 plots the Avg.Error as a function of K. As we can see, the Avg.Error of all DFL methods increase as K increases, and the ComDec is more accurate than other DFL methods. Furthermore, we can also observe that the DR-DFL and ComDec perform better than the E-HIPA and LCS. This verifies that the Bayesian recovery algorithm can provide a more accurate sparse reconstruction than the greedy algorithm for CS recovery.

Conclusions and Future Work
In this paper, a novel multi-channel CS-based DFL method, ComDec is proposed. It can solve the dictionary mismatch problem caused by environmental dynamics in changing environments. The key novelty of ComDec is the making use of the channel diversity of CSI measurements in CS-based DFL. Moreover, the dictionary parameters of multiple channels are assumed to be adjustable. They can adapt to the changes of environment and different channel characteristics. In this manner, we can avoid the site survey process that is typically adopted in fingerprinting-based DFL and improve the robustness of DFL against environmental dynamics. We formulate the CS-based DFL problem as a joint sparse recovery problem, and develop a novel joint sparse recovery algorithm under the variational Batesian inference framework. Furthermore, four methods are presented to reduce its computational complexity. Finally, numerical simulation results validate the effectiveness of the ComDec method.
As future work, we will further investigate the dictionary mismatch problem in CS-based multi-target DFL caused by more realistic and complex environmental dynamics. Additionally, we intend to design and implement an efficient and accurate DFL framework that can utilize the phase information of CSI measurements.
Author Contributions: All authors contributed extensively to the work in this paper. D.Y., Y.G., and N.L. conceived the proposed scheme, derived inference procedure and conducted numerical simulation. D.Y. analyzed the data and wrote the manuscript. Y.G., N.L. and X.Y. provided valuable comments and contributed to the revision of the manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: