Time-of-arrival source localization based on weighted least squares estimator in line-of-sight/non-line-of-sight mixture environments

Abstract

In this article, we propose a line-of-sight/non-line-of-sight time-of-arrival source localization algorithm that utilizes the weighted least squares. The proposed estimator combines multiple sorted measurements using the spatial sign concept, Mahalanobis distance, and Stahel–Donoho estimator, that is, assigning less weight to the samples as they are far from the center of inlier distribution. Also, the eigendecomposition Kendall’s $\tau$ covariance matrix is utilized as the scatter measure instead of the conventional median absolute deviation. Thus, the adverse effects by outliers can be attenuated effectively. To validate the superiority of the proposed methods, the root mean square error performances are compared with that of the existing algorithms via extensive simulation.

Introduction

The aim of the source localization system is to find a geometrical point of intersection using the measurements from each receiver, such as the time difference of arrival (TDOA), time of arrival (TOA), or received signal strength (RSS). Localizing a point source often requires passive and stationary sensors commonly utilized in the areas of radar, sonar, global positioning system, video conferencing, and telecommunication. Even though location estimation problems have been investigated extensively in the existing literature,^1–9 there are still some unresolved problems. One of the key challenges of the localization problem is to estimate the position of the source in dense cluttered non-line-of-sight (NLOS) environments.^10,11 Therefore, we concentrate on localization using robust statistics.

The motivation of this article is as follows. The weighted least squares (WLS) estimator utilizes the weight matrix which is determined by the inverse of the covariance of the measurements. It is well known that the estimation performance of the WLS estimator is closer to the efficiency than the least squares (LS) estimators. ¹² However, the algorithm which utilizes the covariance information of the measurement in the line-of-sight (LOS)/NLOS mixture state has not yet been reported. That is, the WLS estimator in the LOS/NLOS mixture situation, in which case the weight matrix is determined using the covariance, has not yet been developed. Thus, we employ the WLS algorithm for the LOS/NLOS mixture state, where the error distribution of LOS/NLOS state is not known. Although the WLS estimator in the LOS/NLOS mixture situation was already developed in Park and Chang, ¹³ the notable difference between the existing WLS and proposed method is that the proposed method does not require the noise variance information unlike the existing WLS algorithm. Also, the proposed method does not need the statistical testing to discern the outliers. The weight is determined as the inverse of the noise variance in the conventional weighted average method. However, the error distribution of LOS/NLOS state is generally difficult to be estimated. Thus, weights are determined through the spatial sign concept and Mahalanobis distance^14,15 in the proposed weighted average method by attenuating the effects of outliers. Also, the Stahel–Donoho (SD) estimator has been widely used in the robust statistics.^16–18 It was the first estimator with a breakdown point (i.e. the maximum fraction of outliers that the estimator can withstand) close to 50% for any dimension. It has good robustness, as shown in Van Aelst et al., ¹⁶ Wilcox, ¹⁷ and Maronna and colleagues,^18,19 which enables the estimator to be useful for multivariate robust estimation. However, the SD estimator has not yet been adopted in the LOS/NLOS mixture localization context. Thus, we apply the SD estimator to the LOS/NLOS mixture localization problem in this work. The weight in the weighted average using the SD estimator is determined inversely proportional to the outlyingness (the difference between the measurement and model). However,the computation of the outlyingness of the SD estimator requires intensive computational load.^20,21 Hence, we utilize the Mahalanobis distance as the outlyingness measure instead of the outlyingness measure of the existing SD estimator to reduce the computational complexity. Also, the eigendecomposition Kendall’s $\tau$ covariance matrix is adopted to estimate the covariance matrix for the distribution of the inliers. ¹⁴ The proposed methods show the superior root mean square error (RMSE) performances compared to that of the existing methods. The proposed methods outperformed the M-estimator, ²² JMAP-ML,^23,24 and approximated maximum likelihood (ML) ²⁵ estimators in entire NLOS noise regimes. Furthermore, the proposed localization method which uses the eigendecomposition Kendall’s $\tau$ covariance matrix was more superior to that using the median absolute deviation (MAD) as the number of LOS/NLOS sensors increases. The organization of this article is listed as follows: section “Background” investigates important concepts used in this article. Section “Problem formulation” explains the LOS/NLOS mixture source localization problem to be solved in this article. In section “Review of the existing algorithms,” the details of the Kendall’s $\tau$ covariance matrix and SD estimator are dealt with. The proposed localization methods using the spatial sign concept, Mahalanobis distance, and SD estimator are addressed in section “Proposed TOA source localization in LOS/NLOS mixture environments.” The mean square error (MSE) performance of the proposed method is analyzed in section “MSE performance analysis.” The estimation performances of the proposed methods are evaluated via simulation results in section “Simulation results,” comparing them with those of the existing algorithms. Finally, the conclusion is presented in section “Conclusion.”

Background

TOA

The ranging makes use of synchronized transmission and travel time of the signal, whose velocity is known, between the transmitter and receiver, hence time interval between the signal transmission and reception time is utilized to estimate the TOA. In order to obtain the signal reception time, the direct (first) signal path should be determined. The two-step TOA algorithm first estimates the block where the direct signal path exists. ²⁶ Then, it determines the chip position where the direct path is present. The time interval from transmission to reception is converted to distance by multiplying the signal speed and time interval ( ${\tau }_{\mathrm{TOA}}$ ). The performance of TOA estimation method using the ultra-wideband (UWB) signal can be degraded when the synthetic data are utilized. Also, TOA estimation algorithm has disadvantages that the synchronization of the clock and estimation of the light speed are required and the positioning accuracy can be considerably deteriorated in the presence of clock imperfections.

LOS and NLOS

LOS propagation is a characteristic of electromagnetic radiation or acoustic wave propagation. LOS is a type of propagation that can transmit and receive data without any sort of an obstacle between transmitter and receiver. NLOS propagation is the radio transmission across a path that is partially obstructed, usually by a physical object. Obstacles that commonly cause NLOS conditions include buildings, trees, hills, and mountains. Some of these obstructions reflect certain radio frequencies, while some simply absorb or garble the signals. The signal of NLOS path makes the TOA estimation and localization performances much degraded.

Signal attenuation model

Let us consider the UWB signal. ²⁷ In a UWB channel, there are frequency-dependent path loss (PL) defined as follows

$$\mathrm{PL}(f,d)=E\{\underset{f-\Delta f/2}{\overset{f+\Delta f/2}{\int }}|\mathrm{H}(\tilde{f},d){|}^{2}d\tilde{f}\}$$ (1)

where $\mathrm{H}(f,d)$ is the transfer function from antenna connector to antenna connector, d is the distance between the transmitter and receiver, $\Delta f$ is chosen small enough so that diffraction coefficients, dielectric constants can be considered constant within that bandwidth, and the total path loss is obtained by integrating over the whole bandwidth of interest. To simplify computations, we assume that the path loss as a function of the distance and frequency can be written as a product of the terms

$$\mathrm{PL}\begin{array}{c}(f,d)=\mathrm{PL}(f)\mathrm{PL}(d)\hfill \end{array}$$ (2)

It is found in Kunisch and Pamp ²⁸ that

$$\sqrt{\mathrm{PL}(f)}\propto f^{-m}$$ (3)

with m varying between 0.8 and 1.4. The path loss also depends on the distance in which the distance dependence is usually modeled as a power decay law

$$\mathrm{PL}(d)=\mathrm{P}{\mathrm{L}}_0+10n\mathrm{lo}{\mathrm{g}}_{10}(\frac{d}{d_0})$$ (4)

where n is the path loss exponent, $\mathrm{P}{\mathrm{L}}_0$ is the path loss at the reference distance, and $d_0$ is set to 1 m.

Problem formulation

The main idea behind the TOA-based source localization method is to find the position of a source accurately using multiple circles whose centers are the locations of sensors. In the LOS/NLOS mixture source localization context, the measurement equation is represented as

$$r_{i,j}=d_i+n_{i,j}=\sqrt{{(x-x_i)}^2+{(y-y_i)}^2}+n_{i,j}$$ (5)

where $n_{i,j}~(1-ϵ)N(0,{\sigma }_1^2)+ϵN({\mu }_2,{\sigma }_2^2),i=1,2,\dots ,M,j=1,2,\dots ,P$ with M and P denoting the number of sensors and number of samples in the ith sensor.^29–31 Details for equation (5) can be obtained in Park and Chang. ¹³ Squaring equation (5) and rearranging yield the following equation

$$\begin{array}{c}{x}_{i}x+y_{i}y-0.5R+m_{i,j}=0.5({x_i}^2+{y_i}^2-r_{i,j}^2),\hfill \\ i=1,2,\dots ,M,j=1,2,\dots ,P\hfill \end{array}$$ (6)

where $R=x^2+y^2,m_{i,j}=-d_i{n}_{i,j}-\frac{1}{2}{n}_{i,j}^2$ . For convenience, equation (6) can be simply represented in a matrix form as

$$\mathbf{Ax}+{\mathbf{q}}_j={\mathbf{b}}_j,j=1,\dots ,P$$ (7)

where ${\mathbf{q}}_j=[m_{1,j},\dots ,m_{M,j}{]}^T,\mathbf{x}=[xyR{]}^T,$

$$\begin{array}{c}\mathbf{A}=\left(\begin{array}{ccc}{x}_1\hfill & y_1\hfill & -0.5\hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill & ⋮\hfill \\ \multicolumn{1}{c}{x_M}& y_M\hfill & -0.5\hfill \end{array}\right),\mathrm{and}\hfill \\ {\mathbf{b}}_j=\frac{1}{2}\left(\begin{array}{c}{x}_1^2+y_1^2-r_{1,j}^2\hfill \\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{x_M^2+y_M^2-r_{M,j}^2}\end{array}\right).\hfill \end{array}$$

Then, the WLS location estimate is obtained as given by

$${\hat{\mathbf{x}}}_1=({\mathbf{A}}^T\mathbf{WA}{)}^{-1}{\mathbf{A}}^T\mathbf{W}\{\mathbf{f}({\mathbf{b}}_{1:P}^s)\}.$$ (8)

The purpose of this article is to find the multivariate function (f) and weight matrix (W) using ${\mathbf{b}}_{1:P}^s$ , for which the RMSE of the position estimate is minimized, where ${\mathbf{b}}_{1:P}^s=[(b_{1,1:P}^s{)}^T\cdots (b_{M,1:P}^s{)}^T{]}^T$ , $b_{i,1:P}^s$ denotes the sorted samples of $[b_{i,1}\cdots b_{i,P}]$ in the ascending order and $b_{i,j}$ is the jth transformed sample of the ith sensor.

Review of the existing algorithms

We first briefly review the concept of the spatial sign, eigendecomposition Kendall’s $\tau$ covariance matrix, and SD estimator.

Eigendecomposition Kendall’s $\tau$ covariance matrix

The so-called spatial sign is the function defined as follows

$$\mathbf{S}(\mathbf{x})=\{\begin{array}{cc}\frac{\mathbf{x}}{\Vert \mathbf{x}\Vert },\hfill & \mathbf{x}\ne 0\hfill \\ \multicolumn{1}{c}{0,}& \mathbf{x}=0.\hfill \end{array}$$ (9)

The Kendall’s $\tau$ covariance matrix is found as follows ¹³

$$\mathrm{TCM}=\frac{1}{P(P-1)}\sum _{i=1}^P\sum _{j=1}^P\{\mathbf{S}({\mathbf{b}}_i^s-{\mathbf{b}}_j^s){\mathbf{S}}^T({\mathbf{b}}_i^s-{\mathbf{b}}_j^s)\}$$ (10)

where ${\mathbf{b}}_j^s$ denotes the jth column of ${\mathbf{b}}_{1:P}^s$ . As can be seen from equations (9) and (10), as the denominator ( $\Vert {\mathbf{b}}_i^s-{\mathbf{b}}_j^s\Vert$ ) has larger values, the impact from the corresponding samples is smaller. The Kendall’s $\tau$ covariance matrix is further improved by the following estimation procedure: ¹⁴

1. Find the eigenvector estimates of the TCM, that is, U.

2. Estimate the eigenvalues of ${\mathbf{U}}^T{\mathbf{b}}_1^s,\dots ,{\mathbf{U}}^T{\mathbf{b}}_P^s$ using the MAD. ³² Write $\Lambda =\mathrm{diag}({\lambda }_1,\dots ,{\lambda }_M)$ for the estimates.

3. The covariance matrix estimate for the distribution of the inlier measurements is ${\mathbf{C}}_{\tau }=\mathbf{U}\Lambda {\mathbf{U}}^T.$

Hereafter, we name the above covariance matrix as the eigendecomposition Kendall’s $\tau$ covariance matrix.

SD estimator

The SD estimator, proposed independently by Stahel and Donoho,^19–21 is a widely used robust estimator of multivariate location and scatter. It is the first robust estimator with breakdown point close to 50% for any dimension. Let ${\mathbf{B}}^s$ be sorted data matrix consisting of $({\mathbf{b}}_{1:P}^s{)}^T=[{\mathbf{b}}_1^s\cdots {\mathbf{b}}_P^s{]}^T$ . Then, for any $\mathbf{y}\in R^M$ , the SD outlyingness is defined as

$$\mathrm{re}(\mathbf{y},{\mathbf{B}}^s)=\mathrm{su}{\mathrm{p}}_{\mathbf{a}\in S_M}\frac{|{\mathbf{y}}^T\mathbf{a}-\mu ({\mathbf{B}}^s\mathbf{a})|}{\sigma ({\mathbf{B}}^s\mathbf{a})}$$ (11)

where $S_M=\{\mathbf{a}\in R^M:\Vert \mathbf{a}\Vert =1\}$ , $\mu$ is the mean, $\sigma$ is the standard deviation, and sup denotes the supremum. ¹⁵ The SD location estimator is defined as the weighted average as follows

$$T_{\mathrm{SD}}=\frac{\sum _{i=1}^P{w}_i{\mathbf{b}}_i^s}{\sum _{i=1}^P{w}_i}$$ (12)

where $w_i=w(\mathrm{r}{\mathrm{e}}_i)$ , $\mathrm{r}{\mathrm{e}}_i=\mathrm{re}({\mathbf{b}}_i^s,{\mathbf{B}}^s)$ and $w(\mathrm{r}{\mathrm{e}}_i)=I_{(\mathrm{r}{\mathrm{e}}_i\le c)}+(c/\mathrm{r}{\mathrm{e}}_i{)}^2{I}_{(\mathrm{r}{\mathrm{e}}_i>c)}$ for some threshold c. ¹⁶ Also, $I_{(A)}$ is the indicator function where it is 1 when A is true otherwise 0. The SD estimator weights the measurements depending on a measure of outlyingness (the difference between the measurement and model), that is, the weight is given such that it is inversely proportional to the outlyingness in the SD estimator.

Proposed TOA source localization in LOS/NLOS mixture environments

The proposed methods utilize the weighted average in which the weight is determined using spatial sign concept, Mahalanobis distance, and SD estimator to fuse the multiple sorted measurements from each sensor. Then, the covariance matrix is found for these outlier-removed samples.

Fusion of multiple measurements

The weight determination algorithm in the fusion of multiple sorted measurements is based on (1) the spatial sign concept, (2) Mahalanobis distance, and (3) SD estimator. As in equation (8), as the accuracy of $\mathbf{f}({\mathbf{b}}_{1:P}^s)$ is higher, the localization can be performed more accurately. In this work, three methods are utilized to estimate $\mathbf{f}({\mathbf{b}}_{1:P}^s)$ more effectively.

Weighted average based on the spatial sign concept

The fusion of multiple sorted measurements which uses the weighted average based on the spatial sign concept is obtained as follows

$$\mathbf{f}({\mathbf{b}}_{1:P}^s)={\mathbf{b}}_c=\frac{\sum _{j=1}^P{w}_{c,j}{\mathbf{b}}_j^s}{\sum _{j=1}^P{w}_{c,j}}$$ (13)

where $w_{c,j}=1/\Vert {\mathbf{b}}_j^s-\mathrm{med}\{{\mathbf{b}}_{1:P}\}\Vert$ , $\mathrm{med}({\mathbf{b}}_{1:P})$ is $[\mathrm{med}(b_{1,1:P}),\dots ,\mathrm{med}(b_{M,1:P}){]}^T$ , and med is the abbreviation for the median. Let us assume the ascending ordered P components of $\{a_{(1)},\dots ,a_{(P)}\}$ . The median is defined as $a_{((P+1)/2)}$ if P is odd, or $\{a_{(P/2)}+a_{(P/2+1)}\}/2$ if P is even. Note that the weight in the conventional weighted average method is determined as 1/ ${\sigma }_j^2$ , where ${\sigma }_j^2$ is the variance of ${\mathbf{b}}_j^s$ . However, the noise variance of NLOS noise is in general not known, thus the conventional weighted average cannot be applied to the LOS/NLOS mixture localization context. The proposed weighted average using the spatial sign assigns the less weight for the sample which is far from the center of inliers (median). It is well known that the median can be used as the center of the inliers when outliers are contaminated if the contamination ratio does not exceed 50%. ³² Thus, the fused estimate ( ${\mathbf{b}}_c$ ) can attenuate the adverse effects of outliers. The reason why that the sorting is performed for raw measurements is that the effects of outliers, which are ordered in the both ends, can be attenuated more effectively than the non-sorted case because the weights given to the sorted outliers are smaller than that of the non-sorted outliers.

Weighted average based on the Mahalanobis distance

The fusion of multiple sorted measurements which uses the weighted average based on the Mahalanobis distance is obtained as follows

$$\mathbf{f}({\mathbf{b}}_{1:P}^s)={\mathbf{b}}_m=\frac{\sum _{j=1}^P{w}_{m,j}{\mathbf{b}}_j^s}{\sum _{j=1}^P{w}_{m,j}}$$ (14)

where $w_{m,j}=\frac{1}{{({\mathbf{b}}_j^s-\mathrm{med}\{{\mathbf{b}}_{1:P}\})}^T{\mathbf{C}}_{\tau }^{-1}({\mathbf{b}}_j^s-\mathrm{med}\{{\mathbf{b}}_{1:P}\})}$ .

We can attenuate the adverse effects caused by outliers using the Mahalanobis distance that is represented as $\sqrt{{({\mathbf{b}}_j^s-\mathrm{med}\{{\mathbf{b}}_{1:P}\})}^T{\mathbf{C}}_{\tau }^{-1}({\mathbf{b}}_j^s-\mathrm{med}\{{\mathbf{b}}_{1:P}\})}$ because the weight given to the sample is smaller as it is far from the center of inliers (median).

Weighted average based on the SD estimator

The fusion of multiple sorted measurements which uses the weighted average based on the SD estimator is obtained as follows

$$\mathbf{f}({\mathbf{b}}_{1:P}^s)={\mathbf{b}}_{\mathrm{sd}}=\frac{\sum _{j=1}^P{w}_{\mathrm{sd},j}{\mathbf{b}}_j^s}{\sum _{j=1}^P{w}_{\mathrm{sd},j}}$$ (15)

where $w_{\mathrm{sd},j}=I_{(\mathrm{r}{\mathrm{e}}_j\le c)}+(c/\mathrm{r}{\mathrm{e}}_j{)}^2{I}_{(\mathrm{r}{\mathrm{e}}_j>c)}$ , $\mathrm{r}{\mathrm{e}}_j=\sqrt{{({\mathbf{b}}_j^s-\mathrm{med}\{{\mathbf{b}}_{1:P}\})}^T{\mathbf{C}}_{\tau }^{-1}({\mathbf{b}}_j^s-\mathrm{med}\{{\mathbf{b}}_{1:P}\})}$ , and c is determined by tuning procedure. That is, the weight is set to one when the outlyingness (Mahalanobis distance) is smaller than the threshold (c) and determined such that it is inversely proportional to the outlyingness when the outlyingness is larger than the threshold. The difference between the existing SD estimator and the proposed estimator is that the parameter a in equation (11) does not need to be determined because the Mahalanobis distance is utilized instead of equation (11) as the outlyingness measure. The SD estimator is known that the computational load is intensive.^19,21 As can be seen from equation (15), as the dimension of measurement is larger, the determination of a requires the higher computational load because the number of all cases should be considered to find the supremum of the outlyingness. Although the approximation algorithm where the outlyingness is determined by searching over a large number of directions perpendicular to hyperplanes that pass through observations is used, the computational complexity is still much high.^19,21 Thus, we utilize the Mahalanobis distance as the outlyingness measure instead of equation (11) to reduce the computational load of existing SD estimator.

Computation of the covariance matrix for the distribution of inlier measurements

In the previous fusion method, the effects of outliers were removed. In this section, the covariance matrix for these outlier-removed (inlier) measurements is obtained. By the definition of the WLS method, although the covariance matrix for the actual fused estimate (equations (13)–(15)) should be utilized, it is much difficult to determine the corresponding covariance matrix directly. Therefore, we approximate the covariance matrix for the fused estimate as the covariance matrix for the distribution of inlier measurements ( ${\mathbf{C}}_{\tau }$ ) because both covariance matrices are that from the distribution for the inlier measurements. The covariance matrix for inlier measurements is obtained using the procedure of section “Eigendecomposition Kendall’s $\tau$ covariance matrix.”

The WLS method in the LOS/NLOS mixture environments

The location estimator using the WLS estimator in the LOS/NLOS mixture environments (LOS/NLOS-WLS estimator) can be divided into three types as follows

$$\begin{array}{c}\mathrm{LOS}/\mathrm{NLOS}-\mathrm{WLS}-\mathrm{C}:({\mathbf{A}}^T{\mathbf{C}}_{\tau }^{-1}\mathbf{A}{)}^{-1}{\mathbf{A}}^T{\mathbf{C}}_{\tau }^{-1}{\mathbf{b}}_c\hfill \\ \mathrm{LOS}/\mathrm{NLOS}-\mathrm{WLS}-\mathrm{M}:({\mathbf{A}}^T{\mathbf{C}}_{\tau }^{-1}\mathbf{A}{)}^{-1}{\mathbf{A}}^T{\mathbf{C}}_{\tau }^{-1}{\mathbf{b}}_m\hfill \\ \mathrm{LOS}/\mathrm{NLOS}-\mathrm{WLS}-\mathrm{SD}:({\mathbf{A}}^T{\mathbf{C}}_{\tau }^{-1}\mathbf{A}{)}^{-1}{\mathbf{A}}^T{\mathbf{C}}_{\tau }^{-1}{\mathbf{b}}_{\mathrm{sd}}.\hfill \end{array}$$ (16)

The $\mathrm{LOS}/\mathrm{NLOS}-\mathrm{WLS}-\mathrm{C}$ , $\mathrm{LOS}/\mathrm{NLOS}-\mathrm{WLS}-\mathrm{M}$ , and $\mathrm{LOS}/\mathrm{NLOS}-\mathrm{WLS}-\mathrm{SD}$ denote the LOS/NLOS-WLS estimator based on the spatial sign concept, Mahalanobis distance, and SD estimator, respectively. The first-step WLS estimate, ${\hat{\mathbf{x}}}_1$ (LOS/NLOS-WLS-C, LOS/NLOS-WLS-M, and LOS/NLOS-WLS-SD), can be further improved using the two-step WLS estimator ³ which is represented as follows

$${\hat{\mathbf{x}}}_2=({\mathbf{H}}^T{\mathbf{C}}_{\hat{\mathbf{h}}}^{-1}\mathbf{H}{)}^{-1}{\mathbf{H}}^T{\mathbf{C}}_{\hat{\mathbf{h}}}^{-1}\hat{\mathbf{h}}$$ (17)

where

$$\begin{array}{c}\hat{\mathbf{h}}={[[{\hat{\mathbf{x}}}_1{]}_1^2[{\hat{\mathbf{x}}}_1{]}_2^2{[{\hat{\mathbf{x}}}_1]}_3]}^T,{\mathbf{C}}_{\hat{\mathbf{h}}}=\mathbf{D}({\mathbf{A}}^T{\mathbf{C}}_{\tau }^{-1}\mathbf{A}{)}^{-1}\mathbf{D},\mathrm{and}\hfill \\ \mathbf{H}=\left(\begin{array}{cc}1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 1\hfill \\ \multicolumn{1}{c}{1}& 1\hfill \end{array}\right)\hfill \end{array}$$

$[·{]}_k$ means the kth element of [·], x, y are substituted as $[{\mathbf{x}}_1{]}_1,[{\mathbf{x}}_1{]}_2$ in the computation of ${\mathbf{C}}_{\hat{\mathbf{h}}}$ and $\mathbf{D}=\mathrm{diag}[2x2y1]$ . The final closed-form two-step WLS source location estimate is found as follows

$${\hat{\mathbf{x}}}_f={[\mathrm{sgn}({[{\hat{\mathbf{x}}}_1]}_1)\sqrt{{[{\hat{\mathbf{x}}}_2]}_1}\mathrm{sgn}({[{\hat{\mathbf{x}}}_1]}_2)\sqrt{{[{\hat{\mathbf{x}}}_2]}_2}]}^T$$ (18)

where $\mathrm{sgn}(·)$ denotes the sign function.

MSE performance analysis

In this section, the MSE of the proposed localization algorithm is analyzed. Let us consider the WLS method using the spatial sign concept because the MSEs of other algorithms can be derived in the similar manner. The variance of the second-step solution of the proposed WLS method using the spatial sign concept can be represented in the following

$$\mathrm{Var}({\hat{\mathbf{x}}}_2)=({\mathbf{H}}^T{\mathbf{C}}_q^{-1}\mathbf{H}{)}^{-1}$$ (19)

where ${\mathbf{C}}_q=\mathbf{D}({\mathbf{A}}^T{\mathbf{C}}_{\tau }^{-1}\mathbf{A}{)}^{-1}\mathbf{D}$ . When assuming ${\mathbf{C}}_{\tau }$ is similar to the variance of ${\mathbf{b}}_c$ in equation (19)

$${\mathbf{C}}_q\approx \mathbf{D}({\mathbf{A}}^T\{\mathrm{Var}({\mathbf{b}}_c){\}}^{-1}\mathbf{A}{)}^{-1}\mathbf{D}$$ (20)

Because the bias is small, the sum of MSE is approximated as $\mathrm{tr}\{\mathrm{Var}({\hat{\mathbf{x}}}_2)\}$ and $\mathrm{tr}[·]$ is the trace of matrix. Then, the MSE of the WLS estimator using the spatial sign concept is obtained as

$$\mathrm{MSE}({\hat{\mathbf{x}}}_2)\approx \mathrm{tr}[({\mathbf{H}}^T{\mathbf{C}}_q^{-1}\mathbf{H}{)}^{-1}]$$ (21)

The $\mathrm{MSE}({\hat{\mathbf{x}}}_2)$ can be rewritten using the pseudo-inverse and property of trace as follows

$$\begin{array}{c}\mathrm{MSE}({\hat{\mathbf{x}}}_2)\approx \mathrm{tr}[({\mathbf{H}}^{†}{\mathbf{C}}_q\{{\mathbf{H}}^{†}{\}}^T)]\hfill \\ =\mathrm{tr}\{\mathrm{Var}({\mathbf{b}}_c)\{{\mathbf{A}}^{†}{\}}^T\mathbf{D}({\mathbf{H}}^{†}{)}^T{\mathbf{H}}^{†}\mathbf{D}{\mathbf{A}}^{†}\}\hfill \end{array}$$ (22)

where ${\mathbf{H}}^{†}$ is the pseudo-inverse of H. $\mathrm{Var}({\mathbf{b}}_c)$ is the diagonal matrix because the measurement of respective sensor is uncorrelated and obtained as follows

$$\mathrm{Var}({\mathbf{b}}_c)=\mathrm{diag}[{\sigma }_{b,1}^2\cdots {\sigma }_{b,M}^2]$$ (23)

where ${\sigma }_{b,i}^2$ = $\frac{1}{w^2}\sum _{j=1}^P[w_j^2{d}_i^2{\sigma }_1^2·I(j\in \mathrm{LOS})+w_j^2(\frac{1}{2}{\sigma }_2^4+d_i^2{\sigma }_2^2+{\mu }_2^2{\sigma }_2^2)·I(j\in \mathrm{NLOS})]$ , $w=\sum _{j=1}^P\frac{1}{\Vert {\mathbf{b}}_j^s-\mathrm{med}\{{\mathbf{b}}_{1:P}\}\Vert }$ , and $w_j=\frac{1}{\Vert {\mathbf{b}}_j^s-\mathrm{med}\{{\mathbf{b}}_{1:P}\}\Vert }$ . The weight w and $w_j$ were approximated as the constant in the derivation of equation (23) although they are actually random variables for the analytical tractability. The MSE( ${\hat{\mathbf{x}}}_2$ ) is proportional to the magnitude of ${\sigma }_{b,i}^2(i=1,\dots ,M)$ as can be seen from equations (22) and (23) because $({\mathbf{A}}^{†}{)}^T\mathbf{D}({\mathbf{H}}^{†}{)}^T{\mathbf{H}}^{†}\mathbf{D}{\mathbf{A}}^{†}$ is assumed to be positive definite matrix. The weight for the NLOS state is much small as can be seen from Figure 8, thus the RMSE as a function of NLOS noise is almost constant. Also, the RMSE increases as the contamination ratio ( $ϵ$ ) is larger because the sample number which belongs to the NLOS set increases.

Simulation results

In this section, the RMSE performances of the proposed LOS/NLOS mixture source localization methods were compared with those of the M-estimator, ²² JMAP-ML method,^23,24 and approximated ML estimator. ²⁵ In this simulation settings, the source was assumed to be located within a $4\times {10}^6(2000\times 2000)$ ${\mathrm{m}}^2$ region to determine the performance over the entire area. Note that the number of sensors used in this experiment was seven. Next, 10 different source locations were generated with a uniform distribution. The sensors were located as shown in Figure 1. In all, 200 Monte-Carlo simulations were performed for each given standard deviation of the NLOS noise. The standard deviation of the LOS noise of all of the sensors was assumed to be identical. In addition, the single and omni-directional source was assumed to be in the stationary state. The RMSE average was calculated as follows

$$\begin{array}{c}\mathrm{RMSE}\mathrm{average}=\hfill \\ \sqrt{\frac{\sum _{i=1}^{10}\sum _{k=1}^{200}[{({\hat{x}}^k(i)-x(i))}^2+{({\hat{y}}^k(i)-y(i))}^2]}{10\times 200}}\hfill \end{array}$$ (24)

where ${\hat{x}}^k(i),{\hat{y}}^k(i)$ is the estimated position of the source in the ith position set and kth iteration and x(i) and y(i) indicate the ith true position of the source.

Figure 1.

Deployment of sensors.

The localization accuracy as a function of the standard deviation of the NLOS noise is shown in Figure 2 when the radius of the sensor network was 1000 m. In Figure 2(a), the contamination ratio ( $ϵ$ ) was 20%, the standard deviation of the LOS noise ( ${\sigma }_1$ ) was 0.5 m, the bias of the NLOS noise ( ${\mu }_2$ ) was 300 m, sensors 3 and 7 were the LOS/NLOS sensors, and the remaining sensors were LOS sensors at which the number of measurements in each sensor was 20. The threshold c used in the determination of the weight based on the SD estimator was set to 0.5 in the simulation. The RMSE averages of the proposed methods, that is, the LOS/NLOS-WLS-C, LOS/NLOS-WLS-M, and LOS/NLOS-WLS-SD algorithms were superior to those of the approximated ML, M-estimator, and JMAP-ML algorithm. However, the localization performance of the proposed methods using the eigendecomposition Kendall’s $\tau$ covariance matrix was similar to that using the MAD ³² (LOS/NLOS-WLS-C (MAD)). The localization performance of the LOS/NLOS-WLS-M and LOS/NLOS-WLS-SD algorithms are the same because the LOS/NLOS-WLS-SD algorithm is similar to the LOS/NLOS-WLS-M method except for the case where the outlyingness is smaller than the threshold and the threshold utilized in the simulation is small. In Figure 2(b), the contamination ratio was 30%, and the remaining conditions are the same with those of Figure 2(a). Again, the RMSE averages of the proposed methods outperform that of the approximated ML, M-estimator, and JMAP-ML method. Figures 3 and 4 assume the same condition as Figure 2 with the exception that sensors 3, 6, and 7 are the LOS/NLOS sensors in Figure 3 and sensors 3, 5, 6, and 7 are the LOS/NLOS sensors in Figure 4. The proposed LOS/NLOS-WLS methods outperformed the approximated ML, M-estimator, and JMAP-ML method. Note that the proposed methods using the eigendecomposition Kendall’s $\tau$ covariance matrix are more superior to that using the MAD as the number of LOS/NLOS sensors increases as can be seen from Figures 2–4.

Figure 2.

Comparison of RMSE averages of the proposed estimators with that of the existing methods when the sensors 3 and 7 are LOS/NLOS sensors and the remaining sensors are LOS sensors: (a) $ϵ$ : 20%, ${\sigma }_1$ : 0.5 m, ${\mu }_2$ : 300 m and (b) $ϵ$ : 30%, ${\sigma }_1$ : 0.5 m, ${\mu }_2$ : 300 m.

Figure 3.

Comparison of RMSE averages of the proposed estimators with that of the existing methods when the sensors 3, 6, and 7 are LOS/NLOS sensors and the remaining sensors are LOS sensors: (a) $ϵ$ : 20%, ${\sigma }_1$ : 0.5 m, ${\mu }_2$ : 300 m and (b) $ϵ$ : 30%, ${\sigma }_1$ : 0.5 m, ${\mu }_2$ : 300 m.

Figure 4.

Comparison of RMSE averages of the proposed estimators with that of the existing methods when the sensors 3, 5, 6, and 7 are LOS/NLOS sensors and the remaining sensors are LOS sensors: (a) $ϵ$ : 20%, ${\sigma }_1$ : 0.5 m, ${\mu }_2$ : 300 m and (b) $ϵ$ : 30%, ${\sigma }_1$ : 0.5 m, ${\mu }_2$ : 300 m.

Figure 5 shows the RMSE averages of the proposed algorithms as a function of the number of sensors. In this case, the number of sensors increases from 5 to 9 and the number of LOS/NLOS sensors is fixed to three and the number of LOS sensors is increased. The standard deviation of the LOS noise was 0.5 m, that of the NLOS noise was 100 m, the bias was 300 m, and the contamination ratio was 20%. We can see that the RMSE averages of the localization decreased as the number of LOS sensors increased. Meanwhile, Figure 6 shows the RMSE averages of the proposed methods as a function of the number of sensors when the number of the LOS/NLOS sensors increases. The number of the LOS/NLOS sensors is one when the number of sensors is 5 and then increases in parallel with the number of sensors. The RMSE averages of the proposed methods decrease as the number of LOS/NLOS sensors increases. Figure 7 shows the RMSE averages as a function of the radius of the sensor network and the RMSE averages decrease as the radius of the sensor network increases. The results of Figures 5–7 were consistent with that of Godrich et al., ³³ that is, as the number of sensors is larger, the geometric dilution of precision (GDOP) is lower, on the other hand, as the number of sensors is smaller, the GDOP gets higher. The GDOP metric has been used as the indicator of the localization accuracy for given deployment of GPS and the localization accuracy is higher as the GDOP is lower. Also, the GDOP is low when the sensors are far apart, on the contrary, the GDOP is high when the sensors are close together. Figure 8 is the comparison of the weights for the sorted measurements and non-sorted measurements. The weights for sorted measurements are smaller for outliers located on both sides than those for non-sorted measurements. Thus, the outliers in the sorted case can be removed more effectively than the non-sorted case. Table 1 shows the comparison of the computational time of M-estimator and proposed SD estimator using Mahalanobis distance when the contamination ratio is 20% and sensors 3, 6, and 7 are the LOS/NLOS sensors. The M-estimator is computationally efficient. ³⁴ It can be seen that the computational time of the SD estimator using Mahalanobis distance is slightly larger than that of M-estimator.

Figure 5.

Comparison of RMSE averages of the proposed estimators as a function of the number of sensors (when the number of LOS sensors increases).

Figure 6.

Comparison of RMSE averages of the proposed estimators as a function of the number of sensors (when the number of LOS/NLOS sensors increases).

Figure 7.

Comparison of RMSE averages of the proposed estimators as a function of the radius of sensor network.

Figure 8.

Comparison of weights for sorted and non-sorted measurements.

Table 1.

Comparison of computational time.

Algorithm	Computational time (s)
M-estimator	3.3522×10−4
SD estimator using Mahalanobis distance	1.2×10−3

SD: Stahel–Donoho.

Conclusion

Robust TOA source localization methods for multiple measurements using the LOS/NLOS-WLS approach were proposed. The proposed methods utilized the weighted average based on the spatial sign concept, Mahalanobis distance, and SD estimator to fuse multiple sorted measurements from each sensor. Also, we used the eigendecomposition Kendall’s $\tau$ covariance matrix as the scale measure for the distribution of the inlier samples. With the above fused estimate and eigendecomposition Kendall’s $\tau$ covariance matrix, the LOS/NLOS-WLS estimators were obtained. The proposed closed-form LOS/NLOS-WLS methods outperformed the approximated ML, M-estimator, and JMAP-ML method in the entire NLOS noise regimes.