A Novel Mobile Localization Method for Distributed Sensor Network with Non-Line-of-Sight Error Mitigation

Abstract

Mobile localization using distributed sensor network has attracted significant research interest in recent years. One of the main challenges for mobile location estimation is the presence of the non-line-of-sight (NLOS) propagation. It results in erroneous measurements containing NLOS error due to the reflection and diffraction. The NLOS error could seriously reduce the localization accuracy. This paper presents a novel mobile localization method for distributed sensor network in NLOS propagation scenarios. The aim of the method is to mitigate the NLOS error and enhance the localization accuracy. The transition between the LOS and NLOS is described by two-state Markov model. The interacting multiple model frame is employed to estimate the position of unknown node. The probability data association algorithm is used to filter the estimated location. The simulation results show that the proposed method outperforms the traditional filter method in NLOS scenarios.

1. Introduction

Mobile location estimation has attracted significant attention in the recent years. Different localization strategies and methods have been proposed. Recent advances in microelectronic, distributed computing, and wireless communication technologies have allowed rapid development of distributed sensor network [1]. The distributed sensor network (DSN) has a wide range of applications such as indoor location service, environment monitoring, and battlefield surveillance. As an important application, location service by the sensor network has received considerable research interest, since the Global Positioning System (GPS) has relatively higher cost and energy consumption.

The distributed sensor network employs the beacon nodes whose locations are usually assumed to be known a priori to estimate the position of unknown node. There are several fundamental approaches for implementing the localization system such as received signal strength indicator (RSSI), time of arrival (TOA), time difference of arrival (TDOA), and angle of arrival (AOA). Several papers have studied the localization method based on the above measures [2, 3]. However, most of the articles assume that the measurements contain only white Gaussian noise (WGN); that is, the signal propagation state is line-of-sight (LOS) path between the beacon nodes and unknown node. An accurate localization result can be obtained by using the filtering techniques [4, 5] in LOS environment. However, the LOS path may not exist for some special circumstances such as the indoor and forest environments in which the direct propagation path may be blocked by obstacles. The obstacle results in erroneous measurements containing non-line-of-sight (NLOS) error due to the reflection and diffraction. The accuracy of mobile location estimation is a very difficult problem for traditional localization methods because the wireless environment is rough. NLOS error is one of the challenges for accurate localization in range based localization methods.

There are two steps to deal with the NLOS propagation problem. The first is to identify the channel conditions. The second is to mitigate the measurement noise, especially the NLOS error. In [6], the authors propose a sequential probability ratio test based NLOS identification algorithm. This algorithm needs less parameters of the measurement model. The binary hypothesis test and generalized likelihood ratio methods [7] are employed to identify the propagation condition. In [8], the proposed method uses the redundant information in the measurements when more than the minimum number of beacon nodes are present. In [9], a low complexity algorithm is proposed to estimate the channel condition. The proposed method combines the measured TOA and RSS to compute the probability of each channel condition using empirical a priori statistical channel model information. Nonparametric machine learning techniques [10] are proposed to perform NLOS identification. The support vector machine classifier which need not the statistical model of the features is used to distinguish between NLOS and LOS conditions. In [11], a min-max algorithm is used to detect the NLOS condition. The min-max method constructs a bounding box for each beacon node and the estimated distance. This method examines if the estimated location is out of the boundary to decide the propagation condition. After identifying the propagation condition, the NLOS mitigation method is the indispensable step. Some methods discard the NLOS measurements and only the LOS measurements are used to locate the position of unknown node. Some methods provide weighting to reduce the effect of NLOS error on the accuracy of position estimates.

In [12], the localization problem has been approximated by a convex optimization problem using the semidefinite programming relaxation technique. According to the prior probabilities and distributions of the NLOS error, three different cases are addressed and three strategies are proposed, respectively. One-step and two-step localization algorithms [13] are proposed. The constrained optimization is utilized to minimize the cost function which takes account of the NLOS error in these methods. In [14], a residual weighting algorithm (Rwgh) is proposed to mitigate the NLOS error. This method uses the sum of squared residuals of a least squares estimation as the indicator of the accuracy of calculated node coordinates. The least squares method is used to estimate the position according to all possible combinations of measurements and then the estimated location is computed as a weighted combination of these intermediate estimates. Recently, the Kalman frame based interacting multiple model (IMM) methods are investigated. A Kalman based interacting multiple model (KIMM) smoother [15] is proposed for mobile location estimation in rough wireless environment. The extended Kalman [16] based and fuzzy logic [17] based IMM smoother are proposed to solve the nonlinear measurement problem. The NLOS mitigation algorithm [18] for tracking in mixed LOS/NLOS environments is proposed. In this method, the pruning techniques are proposed to reduce the number of states and M-estimate based robust Kalman filter is used to reduce the NLOS effect. In [19], the cubature Kalman filter is used for nonlinear Bayesian estimation to improve the accuracy and filter stability. Most of the above methods need the prior information of the NLOS error.

In this paper, we consider mobile localization based on range measurements in a clutter environment. The measurements contain the WGN and NLOS error. Most of the IMM based methods assume that the NLOS error obeys Gaussian distribution. In this paper, a modified IMM method which employs the probability data association method is proposed. The proposed method focuses on the NLOS mitigate problem without the requirement of accurately distinguishing between LOS and NLOS conditions. Therefore, the proposed method has much applicability.

The rest of the paper is organized as follows. Section 2 introduces some preliminaries such as measurement model, state model, and probabilistic data association method. Section 3 derived the proposed method for NLOS environment. In Section 4, simulation results are presented. The conclusions are given in Section 5.

2. Background

In this section, we introduce some preliminaries for the proposed method. The basic setting of the proposed method is as follows. N beacon nodes are randomly deployed in the field. The beacon nodes transmit the wireless signal. The unknown node receives it and estimates the distance between the beacon node and unknown node.

2.1. Measurement Model

At time k, the distance between ith $(i=1,\dots ,N)$ beacon node and unknown node can be modeled as

$$\begin{array}{cc}& {\widehat{d}}_i\left(k\right)=d_i\left(k\right)+v_i\left(k\right),\hspace{1em}k=\mathrm{1},\dots ,K\\ d_i\left(k\right)=\sqrt{{\left(x\left(k\right)-x_i\left(k\right)\right)}^{\mathrm{2}}+{\left(y\left(k\right)-y_i\left(k\right)\right)}^{\mathrm{2}}},\end{array}$$

(1)

where $d_i(k)$ is the true distance between ith beacon node and unknown node, $[x(k),y(k)]$ is the true position of unknown node at time k, $[x_i(k),y_i(k)]$ is the position of ith beacon node at time k, and $v_i(k)$ is the noise. For LOS condition, $v_i(k)$ is the white Gaussian noise; that is, $v_i~N(\mathrm{0},{\sigma }_{g_i}^{\mathrm{2}})$ . For NLOS condition, $v_i(k)$ is the convolution of the WGN and NLOS error [20].

The probability density function (PDF) of $v_i(k)$ is given by

$$\begin{array}{cc}& f\left(v_i\right)=\{\begin{array}{cc}g\left(v_i\right),\hfill & \text{LOS}\hfill \\ {{\displaystyle \int }}_{-\infty }^{\infty }h\left(x\right)g\left(v_i-x\right)dx,\hfill & \text{NLOS},\hfill \end{array}\end{array}$$

(2)

where $g(·)$ is the PDF of measurement noise; we assume that the measurement noise obeys zero means Gaussian distribution with variance ${\sigma }_{g_i}^{\mathrm{2}}$ ; that is, $N(\mathrm{0},{\sigma }_{g_i}^{\mathrm{2}})$ . $h(·)$ is the PDF of NLOS error.

The NLOS error obeys different distributions such as Gaussian, Uniform, or Exponential distribution for different environments. If the NLOS error obeys Gaussian distribution, that is, $N({\mu }_{N_i},{\sigma }_{N_i}^{\mathrm{2}})$ , the probability density function of $v_i(k)$ is

$$\begin{array}{cc}& f\left(v_i\right)=\frac{\mathrm{1}}{\sqrt{\mathrm{2}\pi \left({\sigma }_{g_i}^{\mathrm{2}}+{\sigma }_{N_i}^{\mathrm{2}}\right)}}\exp \left(-\frac{{\left(v_i-{\mu }_{N_i}\right)}^{\mathrm{2}}}{\mathrm{2}\left({\sigma }_{g_i}^{\mathrm{2}}+{\sigma }_{N_i}^{\mathrm{2}}\right)}\right).\end{array}$$

(3)

If the NLOS error obeys Exponential distribution, that is, $E(u,b)$ , the probability density function of $v_i(k)$ is

$$\begin{array}{cc}& f\left(v_i\right)=u\exp \left(u\left(b-v_i\right)+\frac{u^{\mathrm{2}}{\sigma }_{g_i}^{\mathrm{2}}}{\mathrm{2}}\right)Q\left(\frac{b-v_i+u{\sigma }_{g_i}^{\mathrm{2}}}{{\sigma }_{g_i}}\right),\end{array}$$

(4)

where $Q(·)$ is the Q-function and it is defined as $Q(x)=(\mathrm{1}/\mathrm{2}\pi ){\int }_x^{\infty }\exp (-u^{\mathrm{2}}/\mathrm{2})du$ .

If the NLOS error obeys Uniform distribution, that is, $U(a,b)$ , the probability density function of $v_i(k)$ [20] is

$$\begin{array}{cc}& f\left(v_i\right)=\frac{\mathrm{1}}{b-a}\left[Q\left(\frac{v_i-b}{{\sigma }_{g_i}}\right)-Q\left(\frac{v_i-a}{{\sigma }_{g_i}}\right)\right].\end{array}$$

(5)

The cumulative distribution function (CDF) for $h(v_i)~N(\mathrm{2},{\mathrm{4}}^{\mathrm{2}})$ , $h(v_i)~E(\mathrm{2,1})$ , and $h(v_i)~U(\mathrm{0,3})$ is shown in Figure 1.

Figure 1

The CDF for different NLOS error.

2.2. State Model

The state vector of mobile unknown node at time k can be expressed as

$$\begin{array}{cc}& X\left(k\right)={\left[x\left(k\right),\dot{x}\left(k\right),y\left(k\right),\dot{y}\left(k\right)\right]}^T,\end{array}$$

(6)

where $(x(k),y(k))$ and $(\dot{x}(k),\dot{y}(k))$ denote the position and velocity of the unknown node.

The mobile unknown node is moving on a 2D-plane and the state model is defined as

$$\begin{array}{cc}& X\left(k+\mathrm{1}\right)=\varphi X\left(k\right)+G\gamma \left(k\right),\end{array}$$

(7)

where the state transition matrix ϕ is defined as $\varphi =\left[\begin{array}{cccc}\mathrm{1}& \mathrm{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{1}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right]$ , $G=\left[\begin{array}{cc}{t}^{\mathrm{2}}/\mathrm{2}& \mathrm{0}\\ \begin{array}{c}t\\ \begin{array}{c}\mathrm{0}\\ \mathrm{0}\end{array}\end{array}& \begin{array}{c}\mathrm{0}\\ \begin{array}{c}{t}^{\mathrm{2}}/\mathrm{2}\\ t\end{array}\end{array}\end{array}\right]$ , t is the sample period, and $\gamma (k)$ is the random process noise modeled as zero mean with variance $Q(k)$ . It obeys zero mean and variance ${\sigma }_{\gamma }^{\mathrm{2}}$ .

The switching of propagation condition (LOS or NLOS) is treated as a switching-mode system. In this paper, we employ the two-state Markov process to describe the switching. The two-state Markov is modeled as

$$\begin{array}{cc}& \Pi =\left[\begin{array}{cc}{\pi }_{\mathrm{11}}& {\pi }_{\mathrm{12}}\\ {\pi }_{\mathrm{21}}& {\pi }_{\mathrm{22}}\end{array}\right],\end{array}$$

(8)

where ${\pi }_{ij}$ denotes the Markov transition probability from state i to state j.

The measurement equation is defined as

$$\begin{array}{cc}& Z\left(k\right)=H·X\left(k\right)+\xi \left(k\right),\end{array}$$

(9)

where $H=[\mathrm{1,0}]$ and $\xi (k)$ is the noise.

2.3. Probabilistic Data Association (PDA) Method

The probabilistic data association method is a statistical approach to solve the uncertainty of the measurement [21, 22]. Rather than possibly erring by choosing the data closet to what is expected in order to update the state, PDA hedges its bets by weighting the influence of the various candidate measurements [23]. Since the NLOS error results in the larger accidental error, we employ the PDA method to reduce the effect of NLOS error.

The updated state can be obtained via the combination with different weights:

$$\begin{array}{cc}& \widehat{X}\left(k+\mathrm{1}\hspace{0.17em}|\hspace{0.17em}k+\mathrm{1}\right)=\widehat{X}\left(k+\mathrm{1}\hspace{0.17em}|\hspace{0.17em}k\right)+K\left(k+\mathrm{1}\right){\displaystyle \sum }_{i=\mathrm{1}}^{m_k}{\beta }_i{v}_i\left(k+\mathrm{1}\right),\end{array}$$

(10)

where $\widehat{X}(k+\mathrm{1}\hspace{0.17em}|\hspace{0.17em}k)$ is the predicted state, $v_i(k+\mathrm{1})$ is the residual of ith measurement, and ${\beta }_i$ is the probability that ith measurement is obtained from LOS.

The updated covariance is given by

$$\begin{array}{cc}& P\left(k+\mathrm{1}\hspace{0.17em}|\hspace{0.17em}k+\mathrm{1}\right)={\beta }_{\mathrm{0}}P\left(k+\mathrm{1}\hspace{0.17em}|\hspace{0.17em}k\right)+\left(\mathrm{1}-{\beta }_{\mathrm{0}}\right)P_c+\tilde{P}\left(k+\mathrm{1}\right),\end{array}$$

(11)

where $P_c=P(k+\mathrm{1}\hspace{0.17em}|\hspace{0.17em}k)-K(k+\mathrm{1})HP(k+\mathrm{1}\hspace{0.17em}|\hspace{0.17em}k)$ ,

$$\begin{array}{l}\tilde{P}\left(k+\mathrm{1}\right)=K\left(k+\mathrm{1}\right)[{\displaystyle \sum }_{i=\mathrm{1}}^{m_k}{\beta }_i{v}_i\left(k+\mathrm{1}\right){\beta }_i^T\hfill \\ \phantom{aaaaaaaaa}-v\left(k+\mathrm{1}\right)v_{}^T\left(k+\mathrm{1}\right)\phantom{{\displaystyle \sum }_{i=\mathrm{1}}^{m_k}}]K_{}^T\left(k+\mathrm{1}\right),\hfill \end{array}$$

(12)

and ${\beta }_{\mathrm{0}}$ denotes the probability that none of the measurements is obtained from LOS.

The weights are given by

$$\begin{array}{cc}& {\beta }_i\left(k+\mathrm{1}\right)=\text{Pr}\left({\theta }_i\left(k+\mathrm{1}\right)\hspace{0.17em}|\hspace{0.17em}Z\left(k+\mathrm{1}\right)\right),\hspace{1em}i=\mathrm{0,1},\dots ,R,\end{array}$$

(13)

where ${\theta }_i(k+\mathrm{1})$ denotes that $\widehat{Z}(k+\mathrm{1})$ is obtained from LOS measurement with smallest error. ${\theta }_{\mathrm{0}}(k+\mathrm{1})$ denotes that $\widehat{Z}(k+\mathrm{1})$ is obtained from all NLOS measurements.

The weights can be rewritten as

$$\begin{array}{l}{\beta }_i\left(k+\mathrm{1}\right)\hfill \\ \hspace{1em}=\frac{\mathrm{1}}{c\left(k+\mathrm{1}\right)}p\left(Z\left(k+\mathrm{1}\right)\hspace{0.17em}|\hspace{0.17em}{\theta }_i\left(k+\mathrm{1}\right),m\left(k+\mathrm{1}\right),Z\left(k\right)\right)\hfill \\ \hspace{1em}\hspace{1em}\times \text{Pr}\left({\theta }_i\left(k+\mathrm{1}\right)\hspace{0.17em}|\hspace{0.17em}m\left(k+\mathrm{1}\right),Z\left(k\right)\right),\hfill \end{array}$$

(14)

where $c(k+\mathrm{1})$ is the normalization factor. $p(·)$ is the probability density of the accepted or rejected condition $[{\theta }_i(k+\mathrm{1}),m(k+\mathrm{1}),Z(k)]$ . $\text{Pr}(·)$ is the prior probability.

3. Proposed NLOS Localization Method

In this section, we introduce the details of the proposed method. As shown in Figure 2, the mode probability is firstly computed. Then the mixed state and its covariance are calculated. The state is simultaneously estimated by two PDA filters based on the corresponding modes (LOS and NLOS). We construct different subgroups of range measurements. Each subgroup is used to estimate the position by least squares method. The mode probabilities of the present condition are updated. Finally estimation results of the two parallel PDA filters are combined with their corresponding mode probability.

Figure 2

The data flow of our proposed method.

The processing of each step at $k+1$ time is as follows.

3.1. Computing the Mode Probability

The mode probability is defined as

$$\begin{array}{cc}& {\alpha }_{ij}\left(k\hspace{0.17em}|\hspace{0.17em}k\right)=\frac{{\pi }_{ij}{\alpha }_i\left(k\right)}{{\stackrel{-}{\alpha }}_j},\end{array}$$

(15)

where ${\pi }_{ij}$ is the state transition probability. ${\stackrel{-}{\alpha }}_j$ is the normalized mode probability and it is expressed as

$$\begin{array}{cc}& {\stackrel{-}{\alpha }}_j={\displaystyle \sum }_{i=\mathrm{1}}^{\mathrm{2}}{\pi }_{ij}{\alpha }_i\left(k\right).\end{array}$$

(16)

3.2. Interaction

This step attempts to mix a priori state estimation for two states at time k. The mixed state estimation denotes

$$\begin{array}{cc}& {\widehat{X}}_{\mathrm{0}j}\left(k\hspace{0.17em}|\hspace{0.17em}k\right)={\displaystyle \sum }_{i=\mathrm{1}}^{\mathrm{2}}{\widehat{X}}_i\left(k\hspace{0.17em}|\hspace{0.17em}k\right){\alpha }_{ij}\left(k\hspace{0.17em}|\hspace{0.17em}k\right).\end{array}$$

(17)

The mixed covariance is given by

$$\begin{array}{l}{P}_{\mathrm{0}j}\left(k\hspace{0.17em}|\hspace{0.17em}k\right)={\displaystyle \sum }_{i=\mathrm{1}}^{\mathrm{2}}{\alpha }_{ij}\left(k\hspace{0.17em}|\hspace{0.17em}k\right)\{\phantom{{\left[{\widehat{X}}_i\left(k\hspace{0.17em}|\hspace{0.17em}k\right)-{\widehat{X}}_{\mathrm{0}j}\left(k\hspace{0.17em}|\hspace{0.17em}k\right)\right]}^T}{P}_i\left(k\hspace{0.17em}|\hspace{0.17em}k\right)\hfill \\ \hspace{0.5em}+\left[{\widehat{X}}_i\left(k\hspace{0.17em}|\hspace{0.17em}k\right)-{\widehat{X}}_{\mathrm{0}j}\left(k|k\right)\right]\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.5em}\times {\left[{\widehat{X}}_i\left(k\hspace{0.17em}|\hspace{0.17em}k\right)-{\widehat{X}}_{\mathrm{0}j}\left(k\hspace{0.17em}|\hspace{0.17em}k\right)\right]}^T\}.\hfill \end{array}$$

(18)

3.3. PDA Filter

( 1) Filter Prediction. The prediction step includes the following operations:

$$\begin{array}{cc}& {\widehat{X}}_j\left(k+\mathrm{1}\hspace{0.17em}|\hspace{0.17em}k\right)=\varphi {\widehat{X}}_{\mathrm{0}j}\left(k\hspace{0.17em}|\hspace{0.17em}k\right)\\ P_j\left(k+\mathrm{1}\hspace{0.17em}|\hspace{0.17em}k\right)=\varphi P_{\mathrm{0}j}\left(k\hspace{0.17em}|\hspace{0.17em}k\right){\varphi }^T+GQ\left(k+\mathrm{1}\right)G^T.\end{array}$$

(19)

( 2) Probability Data Association. We firstly construct $R=C_N^{N-\mathrm{1}}$ different subgroup of the measurements. For mode 1, each subgroup estimates the position of unknown node using the following equation:

$$\begin{array}{cc}& Z_{\mathrm{1}}^r\left(k+\mathrm{1}\right)={\left({\left(A^r\right)}^T{A}^r\right)}^{-\mathrm{1}}{\left(A^r\right)}^T{B}^r.\end{array}$$

(20)

For example, for subset $[\mathrm{1,2},3,\dots ,N-1]$ , the matrices $A^r$ and $B^r$ are expressed as

$$\begin{array}{l}{A}^r=\mathrm{2}\left[\begin{array}{ll}\left(x_{\mathrm{2}}\left(k+\mathrm{1}\right)-x_{\mathrm{1}}\left(k+\mathrm{1}\right)\right)\hfill & \left(y_{\mathrm{2}}\left(k+\mathrm{1}\right)-y_{\mathrm{1}}\left(k+\mathrm{1}\right)\right)\hfill \\ \left(x_{\mathrm{3}}\left(k+\mathrm{1}\right)-x_{\mathrm{2}}\left(k+\mathrm{1}\right)\right)\hfill & \left(y_{\mathrm{3}}\left(k+\mathrm{1}\right)-y_{\mathrm{2}}\left(k+\mathrm{1}\right)\right)\hfill \\ ⋮\hfill & ⋮\hfill \\ \left(x_{N-\mathrm{1}}\left(k+\mathrm{1}\right)-x_{N-\mathrm{2}}\left(k+\mathrm{1}\right)\right)\hfill & \left(y_{N-\mathrm{1}}\left(k+\mathrm{1}\right)-y_{N-\mathrm{2}}\left(k+\mathrm{1}\right)\right)\hfill \end{array}\right],\hfill \\ B^r=\left[\begin{array}{l}{\widehat{d}}_{\mathrm{2}}^{\mathrm{2}}-{\widehat{d}}_{\mathrm{1}}^{\mathrm{2}}-\left(x_{\mathrm{2}}^{\mathrm{2}}+y_{\mathrm{2}}^{\mathrm{2}}\right)+\left(x_{\mathrm{1}}^{\mathrm{2}}+y_{\mathrm{1}}^{\mathrm{2}}\right)\hfill \\ {\widehat{d}}_{\mathrm{3}}^{\mathrm{2}}-{\widehat{d}}_{\mathrm{2}}^{\mathrm{2}}-\left(x_{\mathrm{3}}^{\mathrm{2}}+y_{\mathrm{3}}^{\mathrm{2}}\right)+\left(x_{\mathrm{2}}^{\mathrm{2}}+y_{\mathrm{2}}^{\mathrm{2}}\right)\hfill \\ ⋮\hfill \\ {\widehat{d}}_{N-\mathrm{1}}^{\mathrm{2}}-{\widehat{d}}_{N-\mathrm{2}}^{\mathrm{2}}-\left(x_{N-\mathrm{1}}^{\mathrm{2}}+y_{N-\mathrm{1}}^{\mathrm{2}}\right)+\left(x_{N-\mathrm{2}}^{\mathrm{2}}+y_{N-\mathrm{2}}^{\mathrm{2}}\right)\hfill \end{array}\right].\hfill \end{array}$$

(21)

For mode 2, each subgroup estimates the position of unknown node using the residual weighting algorithm to obtain $Z_{\mathrm{2}}^r(k+\mathrm{1})$ , $r=\mathrm{1},\dots ,R$ .

The innovation (or residual) covariance is given by

$$\begin{array}{cc}& S_j\left(k+\mathrm{1}\right)=H{P}_j\left(k+\mathrm{1}\hspace{0.17em}|\hspace{0.17em}k\right)H^T+{\omega }_j\left(k+\mathrm{1}\right).\end{array}$$

(22)

The measurement residual is

$$\begin{array}{cc}& {\epsilon }_j^r\left(k+\mathrm{1}\right)=Z_j^r\left(k+\mathrm{1}\right)-H{\widehat{X}}_j\left(k+\mathrm{1}\hspace{0.17em}|\hspace{0.17em}k\right),\hspace{1em}r=\mathrm{1},\dots ,R.\end{array}$$

(23)

For mode j, the data association of the measurement residual is calculated as

$$\begin{array}{cc}& {\epsilon }_j\left(k+\mathrm{1}\right)={\displaystyle \sum }_{r=\mathrm{1}}^R{\beta }_j^r\left(k+\mathrm{1}\right){\epsilon }_j^r\left(k+\mathrm{1}\right),\end{array}$$

(24)

where ${\beta }_j^r(k+\mathrm{1})=e_j^r/(q+{\sum }_{r=\mathrm{1}}^R{e}_j^r)$ , ${\beta }_j^{\mathrm{0}}(k+\mathrm{1})=q/(q+{\sum }_{r=\mathrm{1}}^R{e}_j^r)$ ,

$$\begin{array}{cc}& e_j^r=\exp \left(-\frac{{\left[{\epsilon }_j^r\right]}^T{S}_j^{-\mathrm{1}}\left(K+\mathrm{1}\right){\epsilon }_j^r}{\mathrm{2}}\right),\\ q=\lambda {\left|\mathrm{2}\pi S_j\left(k+\mathrm{1}\right)\right|}^{\mathrm{1}/\mathrm{2}}\frac{\mathrm{1}-P_G{P}_D}{P_D},\end{array}$$

(25)

where $P_D$ is the target detection probability, $P_G$ is the probability that the correct measurement falls in the validation date, and $\lambda =R/V$ .

( 3) Filter Update. The update step includes the following operations:

$$\begin{array}{l}{\widehat{X}}_j\left(k+\mathrm{1}\hspace{0.17em}|\hspace{0.17em}k+\mathrm{1}\right)={\widehat{X}}_j\left(k+\mathrm{1}\hspace{0.17em}|\hspace{0.17em}k\right)+K_j\left(k+\mathrm{1}\right){\epsilon }_j\left(k+\mathrm{1}\right)\hfill \\ P_j\left(k+\mathrm{1}\hspace{0.17em}|\hspace{0.17em}k+\mathrm{1}\right)\hfill \\ \hspace{1em}={\beta }_{\mathrm{0}}{P}_j\left(k+\mathrm{1}\hspace{0.17em}|\hspace{0.17em}k\right)+\left(\mathrm{1}-{\beta }_{\mathrm{0}}\right)\hfill \\ \hspace{1em}\hspace{1em}\times \left[P_j\left(k+\mathrm{1}\hspace{0.17em}|\hspace{0.17em}k\right)-K_j\left(k+\mathrm{1}\right)H{P}_j\left(k+\mathrm{1}\hspace{0.17em}|\hspace{0.17em}k\right)\right]\hfill \\ \hspace{1em}\hspace{1em}+K_j\left(k+\mathrm{1}\right)\hfill \\ \hspace{1em}\hspace{1em}\times [{\displaystyle \sum }_{r=\mathrm{1}}^R{\beta }_j^r\left(k+\mathrm{1}\right){\epsilon }_j^r\left(k+\mathrm{1}\right){\left({\beta }_j^r\left(k+\mathrm{1}\right)\right)}^T\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}-{\epsilon }_j^{}\left(k+\mathrm{1}\right){\epsilon }_j^T\left(k+\mathrm{1}\right)\phantom{{\displaystyle \sum }_{r=\mathrm{1}}^R{\left({\beta }_j^r\left(k+\mathrm{1}\right)\right)}^T}]K_j^T\left(k+\mathrm{1}\right).\hfill \end{array}$$

(26)

3.4. Computing the Updated Mode Probability

The updated mode probability is given by

$$\begin{array}{cc}& {\alpha }_j\left(k+\mathrm{1}\right)=\frac{{\Lambda }_j\left(k+\mathrm{1}\right){\stackrel{-}{\alpha }}_j}{\alpha },\end{array}$$

(27)

where $\alpha ={\sum }_{j=\mathrm{1}}^{\mathrm{2}}{\Lambda }_j(k+\mathrm{1}){\stackrel{-}{\alpha }}_j$ , ${\stackrel{-}{\alpha }}_j={\sum }_{i=\mathrm{1}}^{\mathrm{2}}{\pi }_{ij}{\alpha }_i(k+\mathrm{1})$ , and

$$\begin{array}{l}{\Lambda }_j\left(k+\mathrm{1}\right)=\mathbb{N}\left({\epsilon }_j\left(k+\mathrm{1}\right);\mathrm{0},S_j\left(k+\mathrm{1}\right)\right)\hfill \\ =\frac{\mathrm{1}}{\sqrt{\mathrm{2}\pi \left|S_j\left(k+\mathrm{1}\right)\right|}}\exp \left(-\frac{{\epsilon }_j^T\left(k+\mathrm{1}\right){\epsilon }_j\left(k+\mathrm{1}\right)}{\mathrm{2}{S}_j\left(k+\mathrm{1}\right)}\right).\hfill \end{array}$$

(28)

3.5. Combination Output

The combination of the two modes can be expressed as

$$\begin{array}{l}{\widehat{X}}_I\left(k+\mathrm{1}\hspace{0.17em}|\hspace{0.17em}k+\mathrm{1}\right)={\displaystyle \sum }_{j=\mathrm{1}}^{\mathrm{2}}{\widehat{X}}_j\left(k+\mathrm{1}\hspace{0.17em}|\hspace{0.17em}k+\mathrm{1}\right)·{\alpha }_j\left(k+\mathrm{1}\right)\hfill \\ P_I\left(k+\mathrm{1}\hspace{0.17em}|\hspace{0.17em}k+\mathrm{1}\right)\hfill \\ ={\displaystyle \sum }_{j=\mathrm{1}}^{\mathrm{2}}{\alpha }_j\left(k+\mathrm{1}\right)\{\phantom{{\left[{\widehat{X}}_j-{\widehat{X}}_I\right]}^T}{P}_j\left(k+\mathrm{1}\hspace{0.17em}|\hspace{0.17em}k+\mathrm{1}\right)\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.17em}\hspace{0.17em}+\left[{\widehat{X}}_j\left(k+\mathrm{1}\hspace{0.17em}|\hspace{0.17em}k+\mathrm{1}\right)-{\widehat{X}}_I\left(k+\mathrm{1}\hspace{0.17em}|\hspace{0.17em}k+\mathrm{1}\right)\right]\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.17em}\hspace{0.17em}\times [{\widehat{X}}_j\left(k+\mathrm{1}\hspace{0.17em}|\hspace{0.17em}k+\mathrm{1}\right)\hfill \\ {\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.17em}\hspace{0.17em}\hspace{1em}\hspace{0.5em}-{\widehat{X}}_I\left(k+\mathrm{1}\hspace{0.17em}|\hspace{0.17em}k+\mathrm{1}\right)]}^T\}.\hfill \end{array}$$

(29)

The output of the proposed method is

$$\begin{array}{cc}& {\left[\widehat{x}\left(k+\mathrm{1}\right),\widehat{y}\left(k+\mathrm{1}\right)\right]}^T=\Omega ·{\widehat{X}}_I\left(k+\mathrm{1}\hspace{0.17em}|\hspace{0.17em}k+\mathrm{1}\right),\end{array}$$

(30)

where $\Omega =\left[\begin{array}{cccc}\begin{array}{c}\mathrm{1}\\ \mathrm{0}\end{array}& \begin{array}{c}\mathrm{0}\\ \mathrm{0}\end{array}& \begin{array}{c}\mathrm{0}\\ \mathrm{1}\end{array}& \begin{array}{c}\mathrm{0}\\ \mathrm{0}\end{array}\end{array}\right]$ .

4. Simulation Results

In this section, we evaluate the performance of the proposed method and compare it with maximum likelihood (ML) method [24], residual weighting (Rwgh) algorithm, ML based Kalman filter (ML-KF) [25], Rwgh based Kalman filter (Rwgh-KF), and interacting multiple model based Kalman filter (IMM-KF) method via the computer simulations. We consider a 2-dimensional square region of size 100 m by 100 m, where 6 beacon nodes and obstacles are randomly deployed. The configuration of this experiment is shown in Figure 3. The communication ranges of the sensor nodes are assumed to be identical and equal to 150 m. The measurement noise $g(v_i)$ is modeled as zero mean white Gaussian noise with variance ${\sigma }_{g_i}^{\mathrm{2}}=\mathrm{1}$ . That is, $g(v_i)$ ~ $N(\mathrm{0},{\mathrm{1}}^{\mathrm{2}})$ . We evaluate the performance of the methods when NLOS error obeys Exponential, Uniform, and Gaussian distributions, respectively. The detection probability is $P_D=0.99$ . The random process noise is modeled as zero mean with variance $Q(k)=\mathrm{diag}\{\mathrm{0}.{\mathrm{5}}^{\mathrm{2}},\mathrm{0}.{\mathrm{5}}^{\mathrm{2}}\}$ . The two-state Markov transition probability matrix is

$$\begin{array}{cc}& \Pi =\left[\begin{array}{cc}\mathrm{0.8}& \mathrm{0.2}\\ \mathrm{0.2}& \mathrm{0.8}\end{array}\right].\end{array}$$

(31)

Figure 3

The configuration of this experiment.

The propagation state with respect to all the beacon nodes is shown in Figure 4 for each time step. The value “1” denotes that the propagation state is NLOS. The value “0” denotes that the propagation state is LOS. We can see that the propagation state varies with time.

Figure 4

The propagation state for each time step.

The following results are obtained from 500 Monte Carlo runs. We employ the root mean square error (RMSE) to evaluate the performance of these methods:

$$\begin{array}{cc}& \text{RMSE}=\sqrt{\frac{\mathrm{1}}{K}{\displaystyle \sum }_{k=\mathrm{1}}^K{\left(x\left(k\right)-\widehat{x}\left(k\right)\right)}^{\mathrm{2}}+{\left(y\left(k\right)-\widehat{y}\left(k\right)\right)}^{\mathrm{2}}},\end{array}$$

(32)

where $[x(k),y(k)]$ is the true position of the unknown node at time k. $[\widehat{x}(k),\widehat{y}(k)]$ is the estimated position of the unknown node at time k.

Figure 5 shows localization result of our proposed method, when the NLOS error obeys Exponential distribution $(E(1/\mathrm{2,1}))$ and the measurement noise is WGN $(N({\mathrm{0,1}}^2))$ . It can be seen that the proposed algorithm is able to estimate the position of unknown node accurately. The relationship between the RMSE and time step for the six methods is shown in Figure 6. The simulation results show that the proposed method owns the best performance. Most of the RMSE is lower than 3 m. The ML method has the worst performance. Rwgh method achieves higher localization accuracy than ML method.

Figure 5

The estimated location of the proposed method.

Figure 6

RMSE versus time step for the six methods.

In Figures 7 and 8, we assume that the NLOS error obeys Exponential distribution, that is, $E(u,1)$ .

Figure 7

The average RMSE versus the parameter u.

Figure 8

The average RMSE versus the standard variance of measurement noise.

Figure 7 shows the impact of parameter u of NLOS error on average RMSE, when the NLOS error obeys Exponential distribution $(E(u,1))$ and the measurement noise is WGN $(N({\mathrm{0,1}}^2))$ . The average RMSE is reduced according to increasing the parameter u. Rwgh and Rwgh-KF have similar performance. The average RMSE of IMM-KF method is lower than ML, ML-KF, Rwgh, and Rwgh-KF methods. Our proposed method has the best performance in comparison with other methods.

Figure 8 shows the relationship between RMSE and the standard variance of measurement noise when NLOS error obeys Exponential distribution, that is, $E(\mathrm{2,1})$ . It is obvious that the proposed method achieves better results than other methods. So the proposed method is robust to the measurement noise.

In Figure 9, we assume that the NLOS error obeys the Uniform distribution, that is, $U(\mathrm{0,1})$ . It shows the cumulative distribution function (CDF) of the localization error for the six methods. We can observe that around 90% localization error of the proposed method achieved an accuracy of less than 10 m. In contrast, 90% localization error of IMM-KF, Rwgh, Rwgh-KF, ML, and ML-KF was achieved at 14 m, 22.5 m, 22 m, 34 m, and 32 m, respectively.

Figure 9

The CDF versus the localization error.

In Figures 10 and 11, we assume that the NLOS error obeys the Gaussian distribution, that is, $N({\mu }_{N_i},{\sigma }_{N_i}^{\mathrm{2}})$ .

Figure 10

The average RMSE versus the mean of NLOS error.

Figure 11

The average RMSE versus the standard variance of NLOS error.

To further study the performance of the proposed method with the NLOS error, the average RMSE with different mean of NLOS error ${\mu }_{N_i}$ is shown in Figure 10, where measurement noise is $N({\mathrm{0,1}}^2)$ and the standard variance of NLOS error ${\sigma }_{N_i}^{}=\mathrm{3}$ . It can be observed that all the methods yield larger localization error as ${\mu }_{N_i}$ increases. The proposed method significantly outperforms other methods and it achieves the highest localization accuracy.

The impact of the standard variance of the NLOS error on the average RMSE is studied in Figure 11, where measurement noise is $N({\mathrm{0,1}}^2)$ and the mean of the NLOS error ${\mu }_{N_i}=\mathrm{2}$ . It indicates that the larger ${\sigma }_{N_i}^{}$ results in the higher localization error. The ML method outperforms the Rwgh method, when ${\sigma }_{N_i}^{}\le \mathrm{3}$ . But the average RMSE of Rwgh is lower than ML, when ${\sigma }_{N_i}^{}>\mathrm{3}$ . And IMM-KF method has better performance than ML, ML-KF, Rwgh, and Rwgh-KF methods. On the whole, the proposed method has the best performance.

5. Conclusion

In this paper, we present a novel approach to deal with non-line-of-sight error. This method does not require the statistical models for NLOS error. The two-state Markov model is firstly used to describe the transition between the LOS and NLOS. The IMM frame based probability data association filter is proposed to mitigate the NLOS error. Simulation results illustrate that the performance of the proposed method outperforms other methods in mixed LOS/NLOS environments.