ANLoC: An Anomaly-Aware Node Localization Algorithm for WSNs in Complex Environments

Abstract

Accurate and sufficient node location information is crucial for Wireless Sensor Networks (WSNs) applications. However, the existing range-based localization methods often suffer from incomplete and detorted range measurements. To address this issue, some methods based on low-rank matrix recovery have been proposed, which usually assume noises follow single Gaussian distribution or/and single Laplacian distribution, and thus cannot handle the case with wider noise distributions beyond Gaussian and Laplacian ones. In this paper, a novel Anomaly-aware Node Localization (ANLoC) method is proposed to simultaneously impute missing range measurements and detect node anomaly in complex environments. Specifically, by utilizing inherent low-rank property of Euclidean Distance Matrix (EDM), we formulate range measurements imputation problem as a Robust ${\ell }_{2,1}$-norm Regularized Matrix Decomposition (RRMD) model, where complex noise is fitted by Mixture of Gaussian (MoG) distribution, and node anomaly is sifted by ${\ell }_{2,1}$-norm regularization. Meanwhile, an efficient optimization algorithm is designed to solve proposed RRMD model based on Expectation Maximization (EM) method. Furthermore, with the imputed EDM, all unknown nodes can be easily positioned by using Multi-Dimensional Scaling (MDS) method. Finally, some experiments are designed to evaluate performance of the proposed method, and experimental results demonstrate that our method outperforms three state-of-the-art node localization methods.

1. Introduction

Wireless Sensor Networks (WSNs) have made great progress and have been widely used in various fields, such as environmental monitoring, intelligent transportation, target tracking, etc. [1,2,3]. The premise that these applications work well is accurate location information acquisition [4,5]. Up to now, many localization methods have been proposed, which could be divided into two categories [6]. One is called range-based localization method that could achieve more accurate positioning, but the computation and communication overhead is large while some hardware support is also required. The other is named range-free localization method that is generally suitable for low power and cost applications, yet their positioning accuracy is low.

In this paper, we mainly focus on the range-based localization methods, which could be described as follows: in WSNs applications, some sensor nodes are randomly disposed, and a few of them called anchors could get actual location information by GPS device or other equipment. Then, with the pair-wise range measurements between nodes and actual location of anchors, the location of all unknown nodes could be easily estimated. In general, range-based localization methods often depend on a large amount of accurate inter-node distance information. However, in practical situations, limited by the energy of sensors or the distribution of nodes in application scenario, only a small amount of inter-node distance information could be obtained [7]. Additionally, due to complex environments and uncertain hardware abnormality, the range measurements inevitably suffer from some errors, which lead to lower positioning accuracy. Usually, these errors can be regarded as a mixture of complex noise and anomaly. In reference [6], the complex noise is considered to be caused by environmental interference, malicious attacks, or other unpredictable factors, and is assumed to be a mixture of Gaussian noise with single distribution and outlier noise with single distribution. Correspondingly, Xiao et al. [6] proposed a method based on low-rank matrix completion to sift the complex noise by adopting both Frobenius-norm regularization and ${\ell }_1$-norm regularization. However, presetting the complex noise into two known noise types was too arbitrary. Actually, the noise distribution in practical applications was often unknown and showed a wider noise distribution beyond Gaussian and Laplacian ones. On the other hand, reference [8] handled one more general application scenario with the co-existence of complex noise and anomaly nodes, where anomaly nodes were defined as ones with abnormal transmission module or unpredictable hardware defects. In reference [8], a ${\ell }_{2,1}$-norm regularization term was employed to detect the anomaly from the corrupted range measurements. However, similar to reference [6], the so-called complex noise is still limited to be modeled as two known noise types, i.e., the mixture of outlier noise and Gaussian noise.

To address this limitation, we propose a novel Anomaly-aware Node Localization (ANLoC) method to simultaneously position the unknown nodes and probe the abnormal nodes in complex environments. Specifically, a Robust ${\ell }_{2,1}$-norm Regularized Matrix Decomposition (RRMD) model is constructed by introducing Mixture of Gaussian distribution and ${\ell }_{2,1}$-norm into the conventional low-rank Matrix Decomposition (MD) model, which can not only well fit the intrinsic low-rank property of Euclidean distance matrix, but also is robust against the node anomaly and a wider range of complex noise distributions beyond Gaussian and Laplacian noises. Our basic idea is to encode the node anomaly and the complex noise as a structural row/column sparsity matrix and a noise matrix, respectively, with entries that satisfy Mixture of Gaussian (MoG) distribution. Here, we prefer to employ MoG distribution as the general noise model due to its universal approximation property to any continuous distribution [9]. Such idea is inspired by some recent noise modeling works including Low-Rank Matrix Factorization (LRMF) [10] and Low-Rank Representation (LRR) [11], and has been verified to be effective in the complex noise scenarios. Note that although anomaly can also be seen as a continuous distribution that is approximated by MoG in theory, we still explicitly model it by using ${\ell }_{2,1}$-norm. The reason is to explicitly detect the location of anomaly node, and thus provide a basis for troubleshooting. Furthermore, based on the popular Expectation Maximization (EM) method [12], an efficient optimization algorithm for solving the proposed RRMD model is designed to obtain the true underlying Euclidean Distance Matrix (EDM). Finally, the actual coordinates of unknown nodes could be easily estimated by employing Multi-Dimensional Scaling (MDS) method [13], and the abnormal nodes can also be detected.

The primary contributions of our work can be summarized as follows.

• A Robust ${\ell }_{2,1}$-norm Regularized Matrix Decomposition (RRMD) model is proposed to jointly estimate the missing range measurements and detect the node anomaly, which takes advantage of the potential relationship between two tasks which could help each other to achieve more accurate performance.
• The MoG distribution is employed to fit the unknown complex noise, which allows the proposed RRMD model to adaptively handle a wider range of noise beyond the existing methods. Meanwhile, an efficient optimization algorithm is designed to solve the proposed RRMD model by adopting the popular EM method.
• A novel Anomaly-aware Node Localization (ANLoC) method is proposed based on the RRMD model, and extensive experiments verify the superior positioning performance of the ANLoC method in the coexistence of node anomaly and complex noise.

The rest of this paper is organized as follows. In Section 2, we introduce the current research advances about the range-based node localization methods and low-rank matrix decomposition methods, respectively. Section 3 describes the notations used in this paper and some related mathematical foundations. In Section 4, the RRMD model is constructed and an optimization algorithm based on the EM method is designed to solve this model. This section also presents the Anomaly-aware Node Localization (ANLoC) method based on the proposed RRMD and the classic MDS. In Section 5, a series of simulation experiments are conducted to evaluate the performance of our proposed method. Finally, the conclusions are drawn in Section 6.

2. Related Work

2.1. Range-Based Node Localization

At present, the typical range-based localization methods include two steps: (1) using certain ranging methods to measure the distance between nodes, such as Time of Arrival (ToA), Time Difference of Arrival (TDoA) and Received Signal Strength Indicator (RSSI); (2) using the range measurements combined with the location of anchor nodes to calculate the position information of unknown nodes. The popular localization method called Maximum Likelihood (ML) is asymptotically efficient [14] with enough data records. Tomic et al. [15] built a new convex estimator that approximated the ML by applying efficient convex relaxations, which reduced the estimation errors. In references [16,17], WSNs localization problem is treated as a variant of EDM recovery problem or graph implementation problem. By using the range measurements between nodes and introducing the slack variable to convert non-convex quadratic distance constraints into linear constraints, the authors formulated the WSNs localization problem as a Semi-Definite Programming (SDP) problem, and designed an efficient optimization method to solve the proposed problem. References [18,19] employed the classical MDS method to map the distance relationship between wireless sensor nodes to low-dimensional space, and generated a relative coordinate map which fitted well the distance relationship between nodes. Then, a few anchor nodes were used to convert the relative position to the global position. However, the above methods did not work well in real application scenarios. In general, these methods required complete and accurate EDM. Unfortunately, due to the requirement of energy-saving and the effects of potential node anomaly and complex noise, the existing range-based node localization methods often suffer from incomplete and detorted range measurements. In addition, the ML method and SDP method could be only used to solve a small-scale problem because of its high computational complexity.

To address the above issues, some matrix completion based node localization methods have been proposed in recent years. Specifically, by taking advantage of the low-rank characteristics of EDM, Feng et al. [20] firstly formulated the range measurement imputation problem as a low-rank EDM Matrix Completion (MC) model, and designed an efficient optimization method to solve the proposed MC model. However, their work is limited to assume the errors contained in the range measurements are Gaussian noise, ignoring the existence of some complex errors caused by node hardware failure, multipath transmission, etc. In reference [6], another kind of errors called outliers was assumed to obey Laplacian distribution. Then, ${\ell }_1$-norm regularization was introduced to deal with it, which effectively improved the positioning accuracy. However, it was also too simple to preset the complex noise to these two known types, and the actual noises should be more complicated in practical applications. Recently, Liu et al. [21] proposed a Linear Bregman Iteration based matrix completion method to localize node position for WSNs. However, this method did not consider the actual scenario under the co-existence of complex noise and anomaly. More importantly, all these matrix completion based node localization methods inevitably involved Singular Value Decomposition (SVD) operation with heavy computation cost. To address these issues, an ANLoC method is proposed in this paper, which formulate the range measurements imputation problem as low-rank matrix decomposition instead of low-rank matrix completion model. Therefore, in the next section, we will introduce the related work on low-rank matrix decomposition.

2.2. Low-Rank Matrix Decomposition

Low-Rank Matrix Decomposition (LRMD) is an important technique in data science, which can uncover the latent manifold structures of data, and thus obtain a low dimensional compression representation. Recently, by decomposing the target matrix into the product of two low-rank matrices, LRMD has been widely used in various fields such as dimensionality reduction, clustering, and matrix recovery.

The original LRMD [22] model can be formulated as

$$\underset{\mathit{U},\mathit{V}}{min} ||\mathit{M}-\mathit{U}{\mathit{V}}^T||{}_{{\ell }_p}$$

(1)

where $\mathit{M}\in {ℝ}^{m\times n}$ is the target matrix to be approximated, $\mathit{U}\in {ℝ}^{m\times r}$ and $\mathit{V}\in {ℝ}^{n\times r}$ are two low-dimensional matrix variables ($r<min\left(m,n\right)$), ${||·||}_{{\ell }_p}$ denotes the ${\ell }_p$-norm, and ${\ell }_1$-norm and ${\ell }_2$-norm are commonly used to make this model robust to Gaussian noise and outlier noise, respectively. When some elements of $\mathit{M}$ are missing, the original LRMD could be changed into a matrix recovery model by adding an orthogonal projection operator [23], which could be formulated as:

$$\underset{\mathit{U},\mathit{V}}{min} ||P_{\mathsf{\Omega }}{\left(\mathit{M}-\mathit{U}{\mathit{V}}^T\right)||}_{{\ell }_p}$$

(2)

where $P_{\mathsf{\Omega }}(·)$ is an orthogonal projection operator defined as

$${\left[P_{\mathsf{\Omega }}\left(\mathit{M}\right)\right]}_{ij}=\{\begin{array}{cc}{M}_{ij},\hfill & \hfill \left(i,j\right)\in \mathsf{\Omega },\\ 0,\hfill & \hfill otherwise,\end{array}$$

(3)

where $\mathsf{\Omega }\subseteq \left[m\right]\times \left[n\right]\left(\left[m\right]=\left\{1,2,\dots ,m\right\},\left[n\right]=\left\{1,2,\dots ,n\right\}\right)$ represents the index set of sampled elements. This model could be applied to various matrix recovery problems such as recommender system [24] and image representation [25], and has lower computational complexity than low-rank matrix completion model with nuclear norm minimization. However, ${\ell }_1$-norm and ${\ell }_2$-norm regularizations are only optimal when the noise follows a Gaussian or Laplacian distribution, which is not in line with the actual situation. To address this limitation, Meng et al. [10] proposed a noise-tolerant LRMD model based on Mixture of Gaussian distribution, which could be described as:

$$\underset{\mathit{U},\mathit{V}}{min} ||P_{\mathsf{\Omega }}{\left(\mathit{M}-\mathit{U}{\mathit{V}}^T\right)||}_{MoG}$$

(4)

where $||·|{|}_{MoG}$ represents that each element of this matrix is modeled by Mixture of Gaussian distribution. Reference [10] has demonstrated the robustness of the MoG model to unknown noise. However, since this model does not explicitly consider the row/column-wise structural anomaly that the sampled matrix may suffer from, it cannot be directly applied to the node localization of WSNs in this paper.

3. Preliminaries

3.1. Notations

We first introduce the important notations used in this paper. All italic letters denote variables and non-italic letters denote constants. Bold uppercase letters denote matrices and bold lowercase letters denote vectors. Specifically, $X_{ij}$ denotes the scalar in the $i$-th row and $j$-th column of $\mathit{X}$. ${\mathit{x}}^i$ and ${\mathit{x}}_j$ represent the $i$-th row and $j$-th column of the matrix $\mathit{X}$ respectively. Additionally, we denote the ${\ell }_p$-norm, Frobenius-norm, and ${\ell }_{2,1}$-norm of matrix $\mathit{X}\in {ℝ}^{m\times n}$ as ${\left({\displaystyle \sum }_{i=1}^m{\displaystyle \sum }_{j=1}^n{|X_{ij}|}^p\right)}^{1/p},{||\mathit{X}||}_F= {\left({\displaystyle \sum }_{i=1}^m{\displaystyle \sum }_{j=1}^n{X}_{ij}{}^2\right)}^{1/2}$, and ${||\mathit{X}||}_{2,1}= {\displaystyle \sum }_{i=1}^m{\left({\displaystyle \sum }_{j=1}^n{X}_{ij}{}^2\right)}^{1/2}$, respectively. ${\mathit{X}}^T$ denotes the transpose of matrix $\mathit{X}$. $\mathcal{N}\left(0,{\sigma }^2\right)$ denotes a Gaussian distribution with mean 0 and variance ${\sigma }^2$. Finally, $\odot$ represents the Hadamard product of two matrices.

3.2. Mathematical Foundation

Definition 1 

[26]: (Proximal Operator) Let $\mathrm{F}\left(\mathit{X}\right)$ be a real convex function defined on $\mathit{X}\in {ℝ}^{m\times n}$, for any $\tau , \mu >0$ and constant matrix $\mathit{M}\in {ℝ}^{m\times n}$, the proximal operator could be defined as:

$$pro{x}_{\tau \mathrm{F}\left(\mathit{X}\right)}\left(\mathit{M}\right)=\arg \underset{\mathit{X}}{min}\tau \mathrm{F}\left(\mathit{X}\right)+\frac{\mu }{2}||\mathit{X}-{\mathit{M}||}_F^2$$

(5)

Definition 2 

[27]: (Structural Thresholding Operator) For any $\tau >0,\mathit{M}\in {ℝ}^{m\times n}$, the proximal operator of ${\ell }_{2,1}$-norm, i.e., structural thresholding operator, is defined as:

$$pro{x}_{\tau {||\mathit{X}||}_{2,1}}\left(\mathit{M}\right)={\mathcal{J}}_{\tau /\mu }\left(\mathit{M}\right)$$

(6)

where

$${\left({\mathcal{J}}_{\tau /\mu }\left(\mathit{M}\right)\right)}^i=max\left\{||{\mathit{M}}^i||{}_2-\tau /\mu ,0\right\}·{\mathit{M}}^i/||{\mathit{M}}^i||{}_2,i=1,2,\dots ,m$$

(7)

Definition 3 

[28]: (Lipschitz Continuous Gradient) A differentiable convex function $\mathrm{F}\left(\mathit{X}\right)$ defined on ${ℝ}^{m\times n}$ is said to have a Lipschitz continuous gradient $\xi$, i.e., for $\forall {\mathit{X}}_1,{\mathit{X}}_2\in {ℝ}^{m\times n}$, there exists a constant $\xi >0$, such that

$${\left({\mathcal{J}}_{\tau /\mu }\left(\mathit{M}\right)\right)}^i=max\left\{||{\mathit{M}}^i||{}_2-\tau /\mu ,0\right\}·{\mathit{M}}^i/||{\mathit{M}}^i||{}_2,i=1,2,\dots ,m$$

(8)

Theorem 1 

[29]: Let ${\mathrm{F}}_1$ and ${\mathrm{F}}_2$ be two lower semi-continuous convex functions on $\mathit{X}\in {ℝ}^{m\times n}$, and ${\mathrm{F}}_2$ is also a differentiable function with Lipschitz continuous gradient $\xi$. Then, for a convex optimization problem defined as

$$\underset{\mathit{X}}{min}{\mathrm{F}}_1\left(\mathit{X}\right)+{\mathrm{F}}_2\left(\mathit{X}\right)$$

(9)

if ${\mathrm{F}}_1+{\mathrm{F}}_2$ is mandatory and strictly convex, then for any initial value ${\mathit{X}}_0$ and $0<\delta <1/\xi$, the iterative sequence ${\mathit{X}}_{k+1}$ generated by the following Equation (10) is the unique solution of problem (9)

$${\mathit{X}}_{k+1}=arg\underset{\mathit{X}}{min}\delta {\mathrm{F}}_1\left(\mathit{X}\right)+\frac{1}{2}||\mathit{X}-{\left({\mathit{X}}_k-\delta \nabla {\mathrm{F}}_2\left({\mathit{X}}_k\right)\right)||}_F^2$$

(10)

4. Anomaly-Aware Node Localization for WSNs

In this section, we first establish a RRMD model and employ EM method to optimize it, and then a novel ANLoC method is proposed based on the RRMD model.

4.1. Euclidean Distance Matrix Completion

4.1.1. Problem Description and RRMD Model Construction

In a typical application scenario, $n$ sensor nodes are randomly deployed in a specific area, whose coordinates can be formulated as $\mathit{X}=\left[{\mathit{x}}_1,{\mathit{x}}_2,\dots ,{\mathit{x}}_n\right]\in {ℝ}^{d\times n}$ ($d$ usually be 2 or 3, indicating the dimension of the coordinate space). Then, the corresponding EDM matrix $\mathit{D}\in {ℝ}^{n\times n}$ could be calculated by

$$D_{ij}=||{\mathit{x}}_i-{\mathit{x}}_j||{}^2={\mathit{x}}_i{}^T{\mathit{x}}_i+{\mathit{x}}_j{}^T{\mathit{x}}_j-2{\mathit{x}}_i{}^T{\mathit{x}}_j,i,j\in \left\{1,2,\dots ,n\right\}$$

(11)

where $D_{ij}$ represents the squared Euclidean distance between the i-th sensor node and the j-th sensor node. Reference [8] has proofed that the rank of matrix $\mathit{D}$ is at most $d+2$, which indicates the EDM matrix has an inherently strict low-rank property. Usually, in practical applications, only a few sensor nodes could obtain their accurate coordinates by loading GPS, and the position information of other unknown nodes need be indirectly calculated by employing some localization methods. Due to irregular node distribution and energy consumption limitation, only the partial range measurements could be collected. Based on the collected incomplete distance measurements, we can only establish a sampled EDM matrix with missing elements. Our goal is to estimate the coordinates of unknown nodes based on the sampled range measurements and the known coordinates of anchor nodes. Figure 1 illustrates the pipeline of node localization for WSNs in complex environments. As shown in Figure 1, the node localization process involves two main steps: EDM recovery and node positioning. At the first step of EDM recovery, by using our proposed RRMD algorithm, we can impute the missing range measurements and obtain the estimated complete EDM matrix, which is also called the true underlying EDM matrix, and then at the second step of node positioning, based on the estimated true underlying EDM matrix, we can further employ MDS method to position each unknown node. Compared to the second step of node positioning, the first step of EDM recovery is more critical and challenging. Specifically, in practical applications, limited by the complex environmental interference, malicious attacks, hardware malfunction, etc., the existence of errors in range measurements is unavoidable, which leads to the destruction of sampled EDM and reduction in localization accuracy. In this paper, we assume that the sampled EDM matrix $\mathit{M}\in {ℝ}^{n\times n}$ consists of the true underlying EDM component $\widehat{\mathit{D}}\in {ℝ}^{n\times n}$, the row structural anomaly matrix component $\mathit{R}\in {ℝ}^{n\times n}$, the column structural anomaly matrix component $\mathit{C}\in {ℝ}^{n\times n},$and the noise matrix component $\mathit{N}\in {ℝ}^{n\times n}$, which is illustrated in Figure 2. The noise is considered to be caused by environmental interference, malicious attacks, etc., which usually follows continuous complex distributions. Moreover, the existence of the anomaly is caused by abnormal transmission module or other unpredictable hardware defects. Specifically, when the receiving module of the nodes fails, it will cause row structural anomaly in the corresponding row in EDM, and when the sending module of nodes fails, the corresponding column in the EDM will have the column structural anomaly [8]. To address the above problems, we need to design an efficient method to obtain a true underlying EDM. Finally, based on this EDM, the actual coordinates of unknown nodes could be easily estimated by employing the classic MDS method.

Therefore, our primary goal is to design an effective EDM reconstruction method to impute the missing range measurements and de-noise the inaccurate ones. Intuitively, based on the low-rank characteristics of EDM, we can classify the EDM reconstruction problem into a standard matrix decomposition model, which is formulated as

$$\underset{\mathit{U},\mathit{V}}{min} ||P_{\mathsf{\Omega }}{\left(\mathit{M}-\mathit{U}{\mathit{V}}^T\right)||}_F.$$

(12)

However, this model can only be applied to the case of single Gaussian noise, and its reconstruction performance will decline dramatically when the sampled EDM is under co-existence of complex noise and anomaly. To this end, we employ MoG distribution and ${\ell }_{2,1}$-norm to improve the model, and then establish a Robust ${\ell }_{2,1}$-norm Regularized Matrix Decomposition (RRMD) model as follows:

$$\begin{array}{c}\underset{\widehat{\mathit{D}},\mathit{R},\mathit{C},\mathit{N},\mathit{U},\mathit{V},\Pi ,\Sigma }{min}{\lambda }_1{||\mathit{R}||}_{2,1}+{\lambda }_2||{\mathit{C}}^T||{}_{2,1}-log\mathcal{L}\left(P_{\mathsf{\Omega }}\left(\mathit{N}\right)|\Pi ,\Sigma \right),\\ s. t. P_{\mathsf{\Omega }}\left(\mathit{M}\right)=P_{\mathsf{\Omega }}(\widehat{\mathit{D}}+\mathit{N}+\mathit{R}+\mathit{C}),\widehat{\mathit{D}}=\mathit{U}{\mathit{V}}^T,{\displaystyle \sum }_{k=1}^{n_c}{\pi }_k=1,{\pi }_k\ge 0, k=1,2,\cdots ,n_c,\end{array}$$

(13)

where $\mathit{M}\in {ℝ}^{n\times n}$ denotes the sampled matrix, $\mathit{R}\in {ℝ}^{n\times n}$ and $\mathit{C}\in {ℝ}^{n\times n}$ represent row structural anomaly and column structural anomaly, respectively, $\mathit{U},\mathit{V}\in {ℝ}^{n\times \left(d+2\right)}$ are two low-rank components of the true underlying EDM matrix $\widehat{\mathit{D}}$, $log\mathcal{L}\left(P_{\mathsf{\Omega }}\left(\mathit{N}\right)|\Pi ,\Sigma \right)$ is the log-likelihood function, $\Pi =\left\{{\pi }_1,{\pi }_{2,\dots ,}{\pi }_{n_c}\right\}$,$\Sigma =\left\{{\sigma }_1^2,{\sigma }_2^2,\dots ,{\sigma }_{n_c}^2\right\}$ and $n_c$ are the parameters of MoG, and ${\lambda }_1$ and ${\lambda }_2$ are the tunable parameters. The key of this model is that we use MoG distribution to smooth any unknown noises and ${\ell }_{2,1}$-norm to detect node anomaly, respectively. The main motivations for this model are as follows. (1) First of all, in response to the complex noise, we introduce a Mixture of Gaussians distribution, which is treated as a universal approximator to any continuous density function. Therefore, the RRMD model is robust against a wider range of complex noise distributions beyond Gaussian and Laplacian noises. (2) Then, as the optimal measure of row sparsity matrix, ${\ell }_{2,1}$-norm could be used to smooth row-wise anomaly and the transpose of column-wise anomaly, thus the abnormal nodes can be detected. (3) Finally, in order to deal with the challenge of data missing, we choose the classic low-rank MD method, which can impute the missing elements with lower computational costs than the low-rank MC method.

In this paragraph, we introduce how to apply MoG distribution to the RRMD model. Without loss of generality, each element $N_{ij}\left(i,j=1,2,\dots ,n\right)$ in the noise matrix $\mathit{N}$ is considered to be from a Mixture of Gaussians distribution, which is defined as

$$\mathrm{𝕡}\left(N_{ij}\right)~{\displaystyle \sum }_{k=1}^{n_c}{\pi }_k\mathcal{N}\left(N_{ij}|0,{\sigma }_k^2\right)$$

(14)

where $\mathcal{N}\left(N_{ij}|0,{\sigma }_k^2\right)$ denotes the Gaussian distribution with mean 0 and variance ${\sigma }^2$. $n_c$ is the number of Gaussian components and ${\pi }_k\ge 0$ represents the mixing proportion where ${\displaystyle \sum }_{k=1}^{n_c}{\pi }_k=1$. Then, the likelihood function can be written as

$$\mathcal{L}\left(P_{\mathsf{\Omega }}\left(\mathit{N}\right)|\Pi ,\Sigma \right)={{\displaystyle \prod }}_{i,j\in \mathsf{\Omega }}{{\displaystyle \sum }}_{k=1}^{n_c}{\pi }_k\mathcal{N}\left(N_{ij},0,{\sigma }_k^2\right)$$

(15)

where $\Pi =\left\{{\pi }_1,{\pi }_{2,\dots ,}{\pi }_{n_c}\right\}$, $\Sigma =\left\{{\sigma }_1^2,{\sigma }_2^2,\dots ,{\sigma }_{n_c}^2\right\}.$ Usually, we use the log-likelihood function instead of this likelihood function for convenient calculation. Then, our aim is to maximize the log-likelihood function to obtain the MoG parameters. Obviously, we can construct the final objective function expressed as Equation (13) by combining the negative log-likelihood function with the ${\ell }_{2,1}$-norm, and then minimize it to get all unknown variables in RRMD model.

4.1.2. Optimizing RRMD via Expectation Maximization Method

In this section, we employ the popular EM method to solve the proposed RRMD model. Firstly, we explain the motivation for using this method. Obviously, in order to gain the estimated values of parameters $\Pi ,\Sigma$, the log-likelihood function $log\mathcal{L}\left(P_{\mathsf{\Omega }}\left(N\right)|\Pi ,\Sigma \right)$ should be maximized. However, for many specific problems, parameters $\Pi ,\Sigma$ cannot be directly calculated due to complex expressions of log-likelihood functions. Fortunately, the EM method can be employed to solve such problems, which has been proven to be effective by many researches. Specifically, let $Z_{ijk}\in \left\{0,1\right\}$ be a set of hidden variables with ${{\displaystyle \sum }}_{k=1}^{n_c}{Z}_{ijk}=1\left(i,j=1,2,\dots ,n\right)$, and $Z_{ijk}=1$ implies the noise $N_{ij}$ comes from the $k$-th Gaussian component. Therefore, the log-likelihood function can be re-written in a form that is easy to optimize as

$$\begin{gathered} log\mathcal{L}\left(P_{\mathsf{\Omega }}\left(\mathit{N}\right)|\Pi ,\Sigma \right) \\ ={{\displaystyle \sum }}_{i,j\in \mathsf{\Omega }}{{\displaystyle \sum }}_{k=1}^{n_c}{Z}_{ijk}\left(log{\pi }_k-log\sqrt{2\pi }{\sigma }_k-\frac{1}{2{\sigma }_k^2}{\left({\mathit{M}}_{ij}-u_i{}^T{v}_j-R_{ij}-C_{ij}\right)}^2\right) \end{gathered}$$

(16)

Next, we will introduce the specific optimization process of estimating the parameters ($\widehat{\mathit{D}},\mathit{R},\mathit{C},\mathit{N},\mathit{U},\mathit{V},\Pi ,\Sigma$) in Equation (13) by alternately conducting the E and M steps.

• E step:

In the E step, we calculate the conditional expectation ${{\rm Y}}_{ijk}$ of $Z_{ijk}$. The specific calculation formula can be given as follows.

$${{\rm Y}}_{ijk}=E\left(Z_{ijk}\right)= \frac{{\pi }_k\mathcal{N}\left(M_{ij}-{\mathit{u}}_i{}^T{\mathit{v}}_j-R_{ij}-C_{ij},0,{\sigma }_k^2\right)}{{{\displaystyle \sum }}_{k=1}^{n_c}{\pi }_{k}N\left(M_{ij}-{\mathit{u}}_i{}^T{\mathit{v}}_j-R_{ij}-C_{ij},0,{\sigma }_k^2\right)}$$

(17)

• M step:

In the M step, we need to optimize the following optimization model as

$$\begin{gathered} \underset{\mathit{U},\mathit{V},\mathit{R},\mathit{C},\Pi ,\Sigma }{\min }{\lambda }_1{||\mathit{R}||}_{2,1}+{\lambda }_2||{\mathit{C}}^T||{}_{2,1} \\ -{{\displaystyle \sum }}_{i,j\in \mathsf{\Omega }}{{\displaystyle \sum }}_{k=1}^{n_c}{{\rm Y}}_{ijk}\left(log{\pi }_k-log\sqrt{2\pi }{\sigma }_k-\frac{1}{2{\sigma }_k^2}{\left(M_{ij}-{\mathit{u}}_i{}^T{\mathit{v}}_j-R_{ij}-C_{ij}\right)}^2\right) \\ s. t. {{\displaystyle \sum }}_{k=1}^{n_c}{\pi }_k=1,{\pi }_k\ge 0 \end{gathered}$$

(18)

A natural idea is to use alternate iteration methods to solve MoG parameters $\Pi ,\Sigma$ and matrix variables $\mathit{U},\mathit{V},\mathit{R},\mathit{C}$, i.e.,

Update$\mathit{\Pi },\mathit{\Sigma }$:

$${\pi }_k=B_k/{{\displaystyle \sum }}_{k=1}^{n_c}{B}_k,\mathrm{where}{B}_k={{\displaystyle \sum }}_{i,j\in \mathsf{\Omega }}{{\rm Y}}_{ijk}$$

(19)

$${\sigma }_k^2=\frac{1}{B_k}{{\displaystyle \sum }}_{i,j\in \mathsf{\Omega }}{{\rm Y}}_{ijk}{\left(M_{ij}-{\mathit{u}}_i{}^T{\mathit{v}}_j-R_{ij}-C_{ij}\right)}^2$$

(20)

Update$\mathit{U},\mathit{V},\mathit{R},\mathit{C}$:

In order to update $\mathit{U},\mathit{V},\mathit{R},\mathit{C}$, we need to solve the following sub-problem:

$$\underset{\mathit{U},\mathit{V},\mathit{R},\mathit{C}}{\min }{\lambda }_1{||\mathit{R}||}_{2,1}+{\lambda }_2||{\mathit{C}}^T||{}_{2,1}+{\displaystyle \sum }_{i,j\in \mathsf{\Omega }}{{\displaystyle \sum }}_{k=1}^{n_c}{{\rm Y}}_{ijk}\frac{{\left(M_{ij}-{\mathit{u}}_i{}^T{\mathit{v}}_j-R_{ij}-C_{ij}\right)}^2}{2{\sigma }_k^2}$$

(21)

If we let

$$W_{ij}=\{\begin{array}{c}\sqrt{{{\displaystyle \sum }}_{k=1}^{n_c}{{\rm Y}}_{ijk}/{\sigma }_k^2},\left(i,j\right)\in \mathsf{\Omega },\\ 0 , otherwise.\end{array}$$

(22)

Then Equation (21) can be reformulated as follows:

$$\underset{\mathit{U},\mathit{V},\mathit{R},\mathit{C}}{\min }{\lambda }_1{||\mathit{R}||}_{2,1}+{\lambda }_2||{\mathit{C}}^T||{}_{2,1}+\frac{1}{2}||\mathit{W}\odot {\left(\mathit{M}-\mathit{U}{\mathit{V}}^T-\mathit{R}-\mathit{C}\right)||}_F^2$$

(23)

Furthermore, we can optimize Equation (23) by solving the following three sub-problems.

1) Update $\mathit{U},\mathit{V}$ by solve the following sub-problem:

$$\underset{\mathit{U},\mathit{V}}{\min }\frac{1}{2}||\mathit{W}\odot {\left(\mathit{M}-\mathit{U}{\mathit{V}}^T-\mathit{R}-\mathit{C}\right)||}_F^2$$

(24)

Obviously, as a typical weighted matrix decomposition problem, it could be solved by various methods, such as ALS [30], WLRA [31] and DN [32].

2) Update $R$ by conducting the following iteration:

$${\mathit{R}}_t={\mathcal{J}}_{{\lambda }_1{\tau }_R}\left({\mathit{R}}_{t-1}-{\tau }_R\mathit{W}\odot \mathit{W}\odot \left({\mathit{R}}_{t-1}\mathit{+}C+\mathit{U}{\mathit{V}}^T-\mathit{M}\right)\right)$$

(25)

where $t$ represents the $t$-th iteration and ${\tau }_R$ is the step size of the proximal gradient.

3) Update $\mathit{C}$ by conducting the following iteration:

$${\left({\mathit{C}}^T\right)}_t={\mathcal{J}}_{{\lambda }_2{\tau }_C}\left({\left({\mathit{C}}_{t-1}-{\tau }_c\mathit{W}\odot \mathit{W}\odot \left(\mathit{R}+{\mathit{C}}_{t-1}+\mathit{U}{\mathit{V}}^T-\mathit{M}\right)\right)}^T\right)$$

(26)

where $t$ represents the $t$-th iteration and ${\tau }_C$ is the step size of the proximal gradient.

Based on the aforementioned analysis, we summarize the whole procedure in Algorithm 1.

Algorithm 1. Proposed Robust ${\ell }_{2,1}$-norm Regularized Matrix Decomposition (RRMD) Algorithm

Input: The sampled EDM matrix $\mathit{M}$, the index set $\mathsf{\Omega }$, the parameters ${\lambda }_1$, ${\lambda }_2$ and threshold $\theta$, the initial number of Gaussian components $n_c$.

Output: The true underlying EDM matrix $\widehat{\mathit{D}}$, the row structural anomaly $\mathit{R},$ and the column structural anomaly $\mathit{C}$.

1. Randomly initialize ${\mathit{U}}_0,{\mathit{V}}_0,{\Pi }_0,{\Sigma }_0$; initialize ${\mathit{R}}_0,{\mathit{C}}_0$ to zero matrix;

2. $t=1$;

3. While not convergence do

4.  (E-step): update $\mathit{{\rm Y}}$ according to Equation (17);

5.  (M-step): update $\Pi ,\Sigma$ according to Equation (19) and Equation (20);

6.  (M-step): update $\mathit{U},\mathit{V}$ according to optimize Equation (24);

7.  (M-step): update $\mathit{R}$ according to Equation (25);

8.  (M-step): update $\mathit{C}$ according to Equation (26);

9.  (Tuning $n_c$): Let $g_i$ and $g_i$ represent the number of $i$-th and$j$-th Gaussian component respectively. if $|{\sigma }_i{}^2-{\sigma }_j{}^2|/\left({\sigma }_i{}^2+{\sigma }_j{}^2\right)<\theta$, then let ${\pi }_i={\pi }_i+{\pi }_j$, ${\sigma }_i{}^2=\left(g_i{\sigma }_i{}^2+g_j{\sigma }_j{}^2)/(g_i+g_j\right)$, $K=K-1$. Lastly, remove ${\pi }_j$ and ${\sigma }_j^2$ from $\Pi ,\Sigma$, respectively.

10.  $t=t+1$;

11. End while

12. $\mathit{R},\mathit{C} \mathrm{and} \widehat{\mathit{D}}=\mathit{U}{\mathit{V}}^T$.

4.2. Anomaly-Aware Node Localization

Based on the proposed RRMD model, all pair-wise range measurements between nodes can be easily obtained. However, it is only the first step of this range-based localization method. The second step of this method is to calculate the location information of the unknown nodes based on the complete distance information and the actual coordinates of anchor nodes, and which can be implemented by using the classical MDS method [6,8,21].

MDS is a typical low-dimensional embedding method, which is originally proposed to solve the curse of dimensionality in the field of machine learning. The MDS method ensures that the distance between samples in the original high-dimensional space can be preserved in low-dimensional space. The input of this method can be high-dimensional features of samples (the pair-wise distance information can be easily obtained from high-dimensional features) or the distance information between samples, while the output is the low-dimensional features of samples in the specified $d$-dimensional space. Therefore, in the sensor networks localization application, we can employ the MDS method to calculate the relative coordinates of each node in the $d$-dimensional space. Then, based on the actual/absolute coordinates of at least $d+1$ anchor nodes and the relative coordinates obtained, the coordinate transformation matrix between the actual coordinates and the relative coordinates can be calculated (with the details provided in Theorem 2). Finally, the relative coordinates can be mapped to the actual ones by using the coordinate transformation matrix.

Next, we introduce the specific steps to implement the above MDS method. Firstly, we let the relative coordinates and actual coordinates of $n$ nodes be represented as $\mathit{T}=\left[{\mathit{t}}_1,{\mathit{t}}_2,\dots ,{\mathit{t}}_n\right]\in {ℝ}^{d\times n}$, $\mathit{A}=\left[{\mathit{a}}_1,{\mathit{a}}_2,\dots ,{\mathit{a}}_n\right]\in {ℝ}^{d\times n}$, respectively, and assume that nodes $1,2,\dots n_a\left(n_a\ge d+1\right)$ are anchor nodes. Secondly, the classic MDS method is employed to calculate the relative coordinates of all nodes. Thirdly, with the actual coordinates of anchor nodes, we have coordinate transformation matrix:

$$\mathit{Q}=\frac{\left[{\mathit{a}}_2-{\mathit{a}}_1,{\mathit{a}}_3-{\mathit{a}}_1,\dots ,{\mathit{a}}_{n_a}-{\mathit{a}}_1\right]}{\left[{\mathit{t}}_2-{\mathit{t}}_1,{\mathit{t}}_3-{\mathit{t}}_1,\dots ,{\mathit{t}}_{n_a}-{\mathit{t}}_1\right]}$$

(27)

Finally, the coordinate of all unknown nodes can be obtained by

$$({\mathit{a}}_i=\mathit{Q}·\left({\mathit{t}}_i-{\mathit{t}}_1\right)+{\mathit{a}}_1,i=n_a+1,n_a+2,\dots n$$

(28)

Based on the aforementioned analysis, the proposed Anomaly-aware Node Localization (ANLoC) algorithm can be summarized as in Algorithm 2.

Algorithm 2. Anomaly-aware Node Localization (ANLoC) Algorithm

Input: The sampled EDM matrix $\mathit{M}$, the index set $\mathsf{\Omega }$, the parameters ${\lambda }_1$, ${\lambda }_2$ and threshold $\theta$, the initial number of Gaussian components $n_c$. The coordinates of anchor nodes $\left\{{\mathit{a}}_1,{\mathit{a}}_2,\dots ,{\mathit{a}}_{n_a}\right\}\left(n_a\ge d+1\right),$ where $n_a$ denotes the number of anchor nodes.

Output: Coordinates of all unknown nodes $\left\{{\mathit{a}}_i|i=n_a+1,n_a+2,\dots n\right\}$.

1. Calculate the true underlying EDM matrix $\widehat{\mathit{D}}$ by using Algorithm 1;

2. Double centering the matrix $\widehat{\mathit{D}}$:

$\mathit{S}=-0.5\times \mathit{J}\widehat{\mathit{D}}\mathit{J}$, where $\mathit{J}=\mathit{I}-1·1^T/n$ and $\mathit{I}$ is identity matrix;

3. Perform SVD decomposition on matrix $\mathit{S}$:

$\left[\mathit{H},\mathit{\Lambda },\mathit{{\rm K}}\right]=svd\left(\mathit{S}\right)$;

4. Calculate relative coordinates:

$\mathit{T}=\left[{\mathit{t}}_1,{\mathit{t}}_2,\dots ,{\mathit{t}}_n\right]=\sqrt{{\mathit{\Lambda }}_d}·{\mathit{H}}_d^T$

where ${\mathit{t}}_i\in {ℝ}^{d\times 1}$, ${\mathit{\Lambda }}_d= \mathit{\Lambda }\left(1:d,1:d\right)$, ${\mathit{H}}_d=\mathit{H}\left(:,1:d\right)$;

5. Calculate the coordinate transformation matrix $\mathit{Q}$:

$\mathit{Q}=\left[{\mathit{a}}_2-{\mathit{a}}_1,{\mathit{a}}_3-{\mathit{a}}_1,\dots ,{\mathit{a}}_{n_a}-{\mathit{a}}_1\right]/\left[{\mathit{t}}_2-{\mathit{t}}_1,{\mathit{t}}_3-{\mathit{t}}_1,\dots ,{\mathit{t}}_{n_a}-{\mathit{t}}_1\right]$;

6. Calculate and output coordinates of all unknown nodes:

${\mathit{a}}_i=\mathit{Q}·\left({\mathit{t}}_i-{\mathit{t}}_1\right)+{\mathit{a}}_1,i=n_a+1,n_a+2,\dots n.$

Theorem 2: 

For sensor nodes localization in$d$-dimensional space, given the absolute positions$\left\{{\mathit{a}}_1,{\mathit{a}}_2,\cdots ,{\mathit{a}}_{n_a}\right\}$of the$n_a$anchor nodes and the relative coordinates$\left\{{\mathit{t}}_1,{\mathit{t}}_2,\cdots ,{\mathit{t}}_n\right\}$of all the$n \left(n\gg n_a\right)$sensor nodes, if$n_a\ge d+1$, then the relative coordinates$\left\{{\mathit{t}}_{n_a+1},{\mathit{t}}_{n_a+2},\cdots ,{\mathit{t}}_n\right\}$can be transformed to the corresponding absolute positions$\left\{{\mathit{a}}_{n_a+1},{\mathit{a}}_{n_a+2},\cdots ,{\mathit{a}}_n\right\}$.

Proof: 

According to reference [13], the absolute positions ${\mathit{a}}_i$ ($i\ge n_a+1$) can be computed according to

$${\mathit{a}}_i=\mathit{Q}·\left({\mathit{t}}_i-{\mathit{t}}_1\right)+{\mathit{a}}_1,i=n_a+1,n_a+2,\cdots ,n$$

(29)

where $\mathit{Q}\in {ℝ}^{d\times d}$ is the unknown coordinate-transform matrix, which should be determined by the following matrix equation:

$$\mathit{Q}·\left[{\mathit{t}}_2-{\mathit{t}}_1,{\mathit{t}}_3-{\mathit{t}}_1,\cdots ,{\mathit{t}}_{n_a}-{\mathit{t}}_1\right]=\left[{\mathit{a}}_2-{\mathit{a}}_1,{\mathit{a}}_3-{\mathit{a}}_1,\cdots ,{\mathit{a}}_{n_a}-{\mathit{a}}_1\right]$$

(30)

Without loss of generality, let $\widehat{\mathit{T}}=\left[{\mathit{t}}_2-{\mathit{t}}_1,{\mathit{t}}_3-{\mathit{t}}_1,\cdots ,{\mathit{t}}_{n_a}-{\mathit{t}}_1\right]$ and $\widehat{\mathit{A}}=\left[{\mathit{a}}_2-{\mathit{a}}_1,{\mathit{a}}_3-{\mathit{a}}_1,\cdots ,{\mathit{a}}_{n_a}-{\mathit{a}}_1\right]$, then we can see that $\widehat{\mathit{T}}\in {ℝ}^{d\times \left(n_a-1\right)}$ and $\widehat{\mathit{A}}\in {ℝ}^{d\times \left(n_a-1\right)}$. Therefore, the above matrix Equation can be regarded as the following equivalent equations:

$${{\displaystyle \sum }}_{l=1}^d{Q}_{il}·{\widehat{T}}_{lj}={\widehat{A}}_{ij},i=1,2,\cdots ,d,\mathrm{and}j=1,2,\cdots ,n_a-1$$

(31)

Obviously, if $n_a\le d$, then the number of unknown entries (i.e., $d^2$) is greater than the number of Equations (i.e., $d\times \left(n_a-1\right)$), so we cannot obtain the determined coordinate-transform matrix $\mathit{Q}$. On the contrary, if $n_a\ge d+1$, we only need to arbitrarily select $d+1$ anchor nodes, then we can obtain the unique solution of $\mathit{Q}$. □

5. Performance Evaluation

In this section, we first describe the experimental setting and evaluation metrics, and then report some extensive experimental results under the different scenarios.

5.1. Experimental Setting

In order to investigate the performance of our proposed ANLoC, some simulation experiments were conducted. All these experiments were designed based on MATLAB 2017a and run on PC with Intel i5-8400 CPU and 16G RAM. Firstly, we randomly disposed 100 sensor nodes in a square area with $100\times 100$ unit (where the unit can be determined according to actual communication condition, such as meter, decimeter, foot, and inch, etc.), and 6 of them were anchor nodes while others were unknown nodes. Let $\mathit{X}\in {ℝ}^{2\times 100}$ be the actual coordinate matrix of all nodes and $\mathit{D}\in {ℝ}^{100\times 100}$ denote the ground truth EDM matrix between nodes. Secondly, we artificially added some error matrices into matrix $\mathit{D}$ to simulate the complex noises and the structural anomaly, and thus obtained a corrupted matrix ${\mathit{D}}_{error}$. Specifically, the complex noise was set to the mixture of the following components: (1) Gaussian noise with mean 0 and variance 100; (2) Gaussian noise with mean 0 and variance 50; and (3) sparse noise with a pollution ratio 1% randomly generated within the range of $\left[0,10000\right]$. Moreover, the structural anomalies were set to be randomly generated within the range of $\left[0,500\right]$ and with a row-wise pollution ratio 3% and a column-wise pollution ratio 3%. Thirdly, sampled matrix $\mathit{M}$ could be obtained by randomly sampling some elements from ${\mathit{D}}_{error}$. Finally, the true underlying EDM matrix can be reconstructed from $\mathit{M}$ via Algorithm 1, and then the actual coordinates of all nodes could be calculated using Algorithm 2. Moreover, we summarized all the simulation parameters in

Table 1. Simulation parameters.

Name	Value	Name	Value
Size of experimental scenario	$100\times 100$	Gaussian noise 1	Mean 0, variance 100
Number of sensor nodes	100	Gaussian noise 2	Mean 0, variance 50
Number of anchors	6	Range of sparse noise	$\left[0,10000\right]$
Range of row anomaly	$\left[0,500\right]$	Range of column anomaly	$\left[0,500\right]$

.

The proposed ANLoC method involves three hyper-parameters ${\lambda }_1, {\lambda }_2$ and $n_c$. It is important to adjust these parameters to proper values. In order to determine the optimal value of these parameters, we employed a 10-fold cross validation procedure and conduct the grid search method. In particular, firstly, let $n_c=5$, we searched for the optimal solution of ${\lambda }_1$ and ${\lambda }_2$ within the range of $\left\{0.001, 0.01, 0.1, 1, 10,100,1000\right\}$ by minimizing the matrix recovery error of EDM. Then, based on the selected optimal parameter ${\lambda }_1$ and ${\lambda }_2$, $n_c$ was determined within the range of $\left\{1, 2, 3, 4, 5,6,7,8,9,10\right\}$. Last but not least, we compared our method with three different competing methods, including NLIRM method [6], SVT-based method [20], and OptSpace-based method [33]. All the involved parameters in these competing methods were optimized by using the same nested 10-fold cross-validation procedure as in our ANLoC method.

5.2. Evaluation Metrics

We selected the following four evaluation metrics to evaluate the node localization performance:

• EDM recovery error:

  $$e_r=||\widehat{\mathit{D}}-\mathit{D}|{|}_F/{||\mathit{D}||}_F$$

  (32)

  where $\widehat{\mathit{D}}$ denotes the reconstructed matrix and $\mathit{D}$ is the ground truth matrix.
• Localization error:

  $$e_l=||\widehat{\mathit{X}}-{\mathit{X}||}_F/n$$

  (33)

  where $\widehat{\mathit{X}}$ denotes the coordinates calculated by Algorithm 2.
• Anomaly recognition accuracy:

  $$p_{ave}=\left(p_{row}+p_{col}\right)/2$$

  (34)

  where $p_{row}$ and $p_{col}$ represent the row anomaly recognition accuracy and column anomaly recognition accuracy, respectively. Specifically, the $p_{row}$ is calculated by

  $$p_{row}=2\times \frac{r_{pre}·r_{rec}}{r_{pre}+r_{rec}},r_{pre}=\frac{r_{tru}}{r_{all}},r_{rec}=\frac{r_{tru}}{r_{act}}$$

  (35)

  where $r_{all}$ represents the number of abnormal rows identified by the method proposed in this paper, $r_{tru}$ indicates the number of correctly recognized rows in $r_{all}$ and $r_{act}$ represents the number of actual abnormal rows; the $p_{row}$ is calculated by

  $$p_{col}=2\times \frac{c_{pre}·r_{rec}}{c_{pre}+r_{rec}},c_{pre}=\frac{c_{tru}}{c_{all}},c_{rec}=\frac{c_{tru}}{c_{act}}$$

  (36)

  where $c_{all}$ represents the number of abnormal columns identified by the method proposed in this paper, $c_{tru}$ indicates the number of correctly recognized columns in $c_{all},$ and $c_{act}$ represents the number of actual abnormal columns.
• Cumulative distribution of localization errors [21]:

  $$e_{CD}=P\left(\Delta l_i\le \sigma \right)$$

  (37)

  where $\Delta l_i=\sqrt{{\left({\widehat{\mathit{x}}}_i-{\mathit{x}}_i\right)}^2+{\left({\widehat{y}}_i-y_i\right)}^2}$ defines the localization error of $i$-th node, $\left({\mathit{x}}_i,y_i\right)$ and $\left({\widehat{\mathit{x}}}_i,{\widehat{y}}_i\right)$ denote the actual coordinates and the estimated coordinates, respectively.

5.3. Experimental Results

In order to investigate the performance of proposed method in complex environments, we conducted the following four experiments under the different application scenarios.

• Scenario 1: Localization without the noise and anomaly

This experiment assumed that all observed range measurements were accurate. Figure 3 shows that the variation of the EDM recovery error and average localization error of WSNs nodes when the EDM is sampled at different proportions. The horizontal axis of the two subgraphs represents the sampling ratio, while the vertical axis in Figure 3a represents the EDM recovery error, and the vertical axis in Figure 3b represents the localization error of WSNs nodes. Comparing Figure 3a with Figure 3b, it is obvious to find that the EDM recovery error and localization error are consistent with the trend of the sampling rate. Under this scenario, when sampling rate reaches 0.3, the EDM recovery error and the localization error of these four methods reach a very low level. This shows that the existing methods and the method proposed in this paper can achieve good positioning accuracy under the ideal condition.

• Scenario 2: Localization with only complex noise

This experiment aims to examine the performance of four methods with only complex noise. Therefore, we only add complex noise into EDM. After sampling the EDM in different proportions, the experimental results are shown in Figure 4. Under this scenario, the SVT-based method and the OptSpace-based method have large EDM recovery error and localization error even if the sampling rate is high, which indicates that the two methods do not handle the complex noise well. When the sampling ratio reaches 0.3, the EDM recovery error of NLIRM could be stabilized at around 0.01, and the localization error of nodes could be less than 0.03. However, the EDM recovery error and localization error of ANLoC could reach 0.001 and 0.01, respectively. Comparing with Figure 3, we find that both the SVT-based method and the OptSpace-based method could not deal with complex noise at all. The NLIRM has a certain ability to resist noise, but it is slightly worse than our ANLoC. Therefore, our ANLoC could perform best under complex noise than others.

• Scenario 3: Localization with only anomaly

The purpose of this experiment is to examine the performance of detecting anomalous nodes in the environment without any noise. In this case, only anomaly is added into EDM. As shown in Figure 5, compared with Figure 3, the EDM recovery error and the localization error of the four methods are all increased when the EDM is destroyed by the anomaly. However, compared with other three methods, ANLoC could have better performance at a lower sampling rate (0.2). In addition, the ANLoC could also detect abnormal nodes. Specifically, Figure 5c shows that when the sampling rate reaches 0.2, the recognition accuracy of abnormal nodes can reach 100%. In comparison with Figure 3, we find that the performance of all methods has been degraded in varying degrees, but our ANLoC still has the best performance. Therefore, when some sensor nodes have abnormities, the ANLoC has good performance for precisely locating unknown nodes and detecting abnormal nodes.

• Scenario 4: Localization with both complex noise and anomaly

This experiment is designed to evaluate the performance of the ANLoC in a complex environment where complex noise and anomaly coexist. The complex noise and the anomaly are both added to the EDM to simulate the impact of the complex real application scenarios. Experimental results are shown in Figure 6. Obviously, compared with the other three methods, our ANLoC method proposed in this paper not only has the lowest EDM recovery error and localization error, but also has high recognition accuracy of abnormal nodes when the sampling rate reaches 0.3. Comparing with Figure 3, we can find that the performance of the SVT-based method and the OptSpace-based method are much worser in complex environment. The NLIRM has less performance degradation, but its EDM recovery accuracy and localization accuracy are still not as high as our ANLoC. In short, our proposed ANLoC has good robustness to complex application scenarios and could detect faulty nodes well. In addition, Figure 7a shows the cumulative distribution of the localization error for the four methods when the sampling rate is 0.5 under this scenario. The probability of the localization error of the ANLoC being less than 0.5 is 98% and the probability of being less than 1 is 100% while the other three methods are all below 30%. In addition, Figure 7b shows the intuitionistic localization results of all nodes in this case when the sampling rate is 0.5.

5.4. Effects of the Proposed Strategies

To analyze the effects of our proposed strategies, experiments are designed to compare the proposed method with three partial deleted methods. Actually, our RRMD can be treated as LRMD + MoG + ${\ell }_{2,1}$-norm. Then, the remaining three types of partial deleted methods are as follows. (1) Pure LRMD: do not consider any noise or anomaly; (2) LRMD + MoG: only consider complex noise; (3) LRMD + ${\ell }_{2,1}$-norm: only consider the anomaly. What’s more, these experiments are set under the scenario with both complex noise and anomaly. Finally, all experimental results are shown in Figure 8a–c. Obviously, our ANLoC considers both complex noise and anomaly, which leads to its low EDM recovery error and localization error. Moreover, when the sampling rate reaches 0.3, it could accurately detect abnormal nodes. Due to too ideal assumptions, pure LRMD is the worst performer, which means that it cannot be applied in actual situation. Although LRMD + MoG could handle the error of EDM well, it does not have the ability to detect abnormal nodes. Specially, LRMD + ${\ell }_{2,1}$-norm could hardly detect abnormal nodes in this case. Naturally, we speculate the noise that is not unprocessed may have a negative effect on the performance for detecting abnormal nodes in this method. Therefore, another experiment is designed with LRMD + ${\ell }_{2,1}$-norm when the EDM is only destroyed by anomaly, and the result of anomaly recognition is drawn in Figure 8d. The pleasure is that this supplementary experiment validates our conjecture. In conclusion, these experiments demonstrate the effectiveness of the MoG and ${\ell }_{2,1}$-norm strategies.

5.5. Localization for Large-Scale Scenario

In a practical application scenario, there is a case that when the size of the localization scenario is much larger than ranging length of the sensor nodes and only a few range measurements between inter-nodes could be obtained. As a result of that, the sampled EDM based on range measurements will be very sparse, which leads to poor performance of our method. To solve this problem, a large-scale localization method is proposed in this paper. For the convenience of discussion, we simplify the application scenario as shown in Figure 9. Suppose there is a $180\times 100$ unit area divided into three parts I, II, and III from left to right, and among them, there are 80, 20, 80 nodes, respectively. Of the 180 nodes, there are 6 anchor nodes and 6 relay nodes. Then, regions I and II are combined into a localization sub-region $\mathcal{M}$ while regions II and III are combined into a localization sub-region $\mathcal{N}$. The large-scale process is as follows. (1) Calculate the actual coordinates of all nodes within sub-region $\mathcal{M}$ by ANLoC. (2) Based on step one, the actual coordinates of relay nodes are gained. Taking them as anchor nodes, we employ ANLoC again to get the coordinates of the remaining unknown nodes in region $\mathcal{N}$.

Further, in order to verify the feasibility of the above-mentioned large-scale scenario localization method, we designed an experiment with complex noise and anomaly. Moreover, the EDM sampling rates of region $\mathcal{M}$ and $\mathcal{N}$ are both set 0.5. As Figure 10 shows, the experimental results confirm that our large-scale localization method can be well applied to large-scale scenarios.

6. Conclusions

In this paper, we aim to estimate node position of WSNs in complex environments. Considering that the coexistence of potential node anomaly and complex noise in practical applications, we propose a novel Anomaly-aware Node Localization (ANLoC) method to address this task. Specifically, the proposed ANLoC method is divided into two steps. First, a Robust ${\ell }_{2,1}$-norm Regularized Matrix Decomposition (RRMD) model is designed to simultaneously detect node anomaly and impute the missing range measurements. Second, based on the imputed EDM matrix, all unknown nodes can be easily estimated by using the classic MDS method. The extensive experiments demonstrated our proposed ANLoC method consistently outperforms three state-of-the-art localization methods in terms of EDM recovery error, localization error, and anomaly recognition accuracy. Additionally, our proposed ANLoC method also can extend to a large-scale localization scenario. Future work will focus on extending our ANLoC approach to an incremental version, making it suitable for handling dynamic positioning of mobile sensor nodes.