{'fraction_non_alphanumeric': 0.03703187112071971, 'fraction_numerical': 0.011762791454541227, 'mean_word_length': 4.471635716102773, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 13, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Mobile Localization in Non-Line-of-Sight Using Constrained Square-Root Unscented Kalman Filter

Localization and tracking of a mobile node (MN) in non-line-of-sight (NLOS) scenarios, based on time of arrival (TOA) measurements, is considered in this work. To this end, we develop a constrained form of square root unscented Kalman filter (SRUKF), where the sigma points of the unscented transformation are projected onto the feasible region by solving constrained optimization problems. The feasible region is the intersection of several discs formed by the NLOS measurements. We show how we can reduce the size of the optimization problem and formulate it as a convex quadratically constrained quadratic program (QCQP), which depends on the Cholesky factor of the \textit{a posteriori} error covariance matrix of SRUKF. As a result of these modifications, the proposed constrained SRUKF (CSRUKF) is more efficient and has better numerical stability compared to the constrained UKF. Through simulations, we also show that the CSRUKF achieves a smaller localization error compared to other techniques and that its performance is robust under different NLOS conditions.


I. INTRODUCTION
N ETWORK-based radio localization has received great attention in recent years due to limitations of the global positioning system (GPS) in indoor places and dense urban areas, and finds numerous applications in surveillance, security, etc. [1]. In this technology, radio signals exchanged between a mobile node (MN) and fixed reference nodes (RN) with known positions 1 , are exploited to determine the unknown location of the MN. There are several different measurements that can be used for localization, e.g., time of arrival (TOA), time difference of arrival (TDOA), received signal strength (RSS), angle of arrival (AOA), and a hybrid of these data. Amongst different localization techniques, TOA-based methods in which the MN and RNs are synchronized are usually preferred, especially in the context of IEEE 802. 15.4a which exploits ultra band-with (UWB) technology [2]. Indeed, measurement of TOA can be done accurately with UWB signaling due to its fine timing resolution and robustness against multipath and fading.
One of the main challenges in radio localization is the non-line of sight (NLOS) problem, which occurs due to the blockage of the direct sight between the MN and RNs. In a NLOS situation, either due to reflection of the radio waves through scatterers or penetration through blocking objects, the travel time of the received signals increases [3], [4], [5]. Consequently, the NLOS error of each measured TOA needs to be modelled as a random variable with a positive bias, which can be quite large [6]. The first step in dealing with the NLOS problem is to detect the NLOS measurements and, if necessary, to discard them. To this end, some techniques estimate the variance of each measurement, and if it is above a given threshold, the corresponding link is identified as NLOS [7], [8], [9]. NLOS identification techniques using signal features have also been proposed for UWB applications [4], [10], [5]. If, however, the NLOS measurements cannot be discarded due to an insufficient number of LOS measurements for unambiguous localization, the next step is to mitigate their effect through further processing.
There are numerous works focusing on NLOS mitigation for the localization of stationary nodes, which are mostly based on (memoryless) constrained optimization techniques. In these approaches, the position of the MN is constrained to be within the convex hull formed by the intersection of multiple discs, each disc being centered at one of the NLOS RNs and with a radius equal to the corresponding measured range. By restricting the MN position in this way and by employing the LOS measurements in the cost function, the unknown location can be found through solving a constrained optimization problem. For a survey on TOA-based memoryless localization in NLOS scenarios, see [6] and the references therein.
For an MN with available dynamic model, filtering techniques are preferred compared to memoryless methods. This is especially the case when data from inertial measurements units (IMU) are used in parallel with range information for tracking purposes [11], [12]. Some methods apply Kalman filter preprocessing on measured TOAs to smooth out the effect of the variances of the NLOS biases, while scaling the covariance matrix in an extended Kalman filter (EKF) to further mitigate their means [7], [8], [9]. However, these approaches can only achieve a moderate performance for large NLOS biases. In [13], [14], it is assumed that the mean and variance of the NLOS biases are known; in practive, however, this information is not available accurately beforehand unless prior field measurements are obtained.
Some other approaches regard the NLOS bias as a nuisance parameter and try to estimate its distribution using Kernel density estimation (KDE). In [15], a robust semi-parametric EKF is proposed for NLOS mitigation of an MN. The performance of this technique is improved by the interacting multiple model (IMM) algorithm in [16]. Although considered for TOA measurements, these techniques are also suitable when AOA, arXiv:1405.0212v1 [stat.AP] 1 May 2014 RSS or a hybrid of these are employed. However, in addition to high computational cost, the performance of KDE still depends on how well it can model the distribution of the NLOS biases. It is claimed that for cellular applications, the performance is only satisfactory when the ratio of NLOS to LOS measurements is less than a half and a higher ratio might result in divergence of KDE algorithms [16].
In some other techniques, the random NLOS biases are considered as parameters in the state vector, to be jointly estimated with other state parameters [17], [18], [19], [20], while the NLOS bias variation over time is modelled as a random walk. The technique in [17] uses EKF, while [18] and [19] use particle filters (PFs) that generally have a high computational cost. In [20], an improved EKF is used where bound constraints on the NLOS biases are enforced for improving the localization accuracy. Although the above techniques can mitigate the effect of NLOS biases to some extent, their performances might not be good due to the mismatch between the random walk model and the physical reality, which is unavoidable considering the unpredictable nature of the biases.
In this work, we propose an efficient square root unscented Kalman filter (SRUKF) with convex inequality constraints for localization of an MN in NLOS situation. The proposed constrained SRUKF (CSRUKF) is mainly based on a combination of the SRUKF in [21] for unconstrained problems and the constrained UKF in [22]. In our proposed algorithm, similar to some memoryless approaches, the NLOS measurements are removed from the observation vector and are employed instead to form a closed convex constraint region [6]. At each time step, we use a SRUKF to estimate the state vector and compute the Cholesky factor of the error covariance matrix. To impose the constraints onto the estimated quantities, as proposed in [22], the sigma points of the unscented transformation may need to be projected onto the feasible region by solving a convex quadratically constrained quadratic program (QCQP). However, we show that the projection can be done in a more efficient and numerically stable way by solving a QCQP with reduced size, in which the cost function depends on the Cholesky factor of the a posteriori error covariance matrix, readily obtained from the SRUKF. Through simulations, our proposed algorithm is shown to achieve a good localization performance under different NLOS scenarios. In particular, in severe NLOS conditions, i.e., when there is only one or no LOS measurement, our method achieves a superior performance compared to other benchmark approaches. Another positive aspect is its robustness to false alarm (FA) errors in NLOS identification, which makes it suitable for practical applications where such errors may be inevitable.
The organization of the paper is as follows: In Section II, the system model is described and the problem formulation is presented. The proposed constrained SRUKF algorithm is developed in Section III, along with a discussion of computational complexity. The simulation results and comparisons with different algorithms are given in Section IV. Finally, Section V concludes the paper.
Notation: Small and capital bold letters represent vectors and matrices respectively. The vector 2-norm operation is denoted by · , and (·) T and (·) −1 stand for matrix trans-pose and inverse operations, respectively. A diagonal matrix with entries x 1 , . . . , x M on the main diagonal is denoted by diag(x 1 , . . . , x M ). For i ≤ j, q(i : j) denotes a vector of size j − i + 1 obtained by extracting the i-th to j-th entries of vector q, inclusively. The symbol I denotes an identity matrix of appropriate dimension. For a positive semi-definite Hermitian matrix R, R 1/2 denotes its unique positive semidefinite definite square root, i.e. such that R 1/2 R 1/2 = R [23].

II. SYSTEM DESCRIPTION AND PROBLEM STATEMENT A. System Model
Consider a network of M fixed RNs and one MN, distributed on a 2-dimensional (2D) plane and exchanging timing signals via wireless links. With reference to a Cartesian coordinate system in this plane, let a i ∈ R 2 denote the known position vector of the i-th RN, where i ∈ {1, . . . , M }, while x k ∈ R 2 and v k ∈ R 2 denote the unknown position and velocity vectors of the MN at discrete time instant k, respectively. Let the state vector be s k = [x T k , v T k ] T ∈ R 4 , which includes the position and velocity components of the MN. The motion model is assumed to be where the matrices F and G are and δt is the time step duration. The vector w k−1 ∈ R 2 in (1) is a zero-mean white Gaussian noise process with diagonal covariance matrix Q = σ 2 w I. In this work, we consider TOA-bsed localization, in which the range between the MN and each RN is obtained by multiplying the time of flight of the radio wave by the speed of light. If the MN and RNs are accurately synchronized, then a one-way ranging scheme can be used; otherwise, a two-way ranging protocol may be employed where the relative clock offsets are removed from the TOA measurements [24]. Let L k and N k denote the index sets of the RNs that are identified as LOS and NLOS nodes at time instant k, respectively. The range measurements can thus be represented as where h i (s k ) = x k − a i , n i k is the measurement noise and b i k is a positive random NLOS bias, which is usually considered independent from n i k . The noise terms n i k , for i = 1, . . . , M , are modelled as independent white Gaussian processes, with zero-mean and known variance σ 2 n . The distribution of each NLOS biases b i k is time-varying due to the movement of the MN and other objects in the area. In the literature, different distributions have been considered for these biases, for instance: exponential [25], [26], shifted Gaussian [14], and uniform [4] are widely employed. However, having a priori knowledge about the distribution of the NLOS biases requires preliminary field measurements, which may not be possible in practical applications. Therefore, in this work, we do not make any specific assumption about the distribution of the NLOS biases, although we suppose that the NLOS links are identified accurately at every time instant.

B. Problem Formulation
The state vector s k and the NLOS biases b i k for i ∈ N k are the unknown parameters in the above model. Representing the NLOS biases by a simple dynamic model such as a random walk, as considered in [18], [19], may not be an accurate approximation, and thus this approach is not used in this work. To simplify the problem and reduce the number of unknowns, we eliminate the NLOS measurements from the observation vector, and instead use the information carried out by the biases to restrict the position of the MN within a certain range. For instance, in many applications, it can be assumed that the TOA measurement noise n i k is small compared to b i k (especially in UWB ranging), which implies that b i k + n i k ≥ 0 [6]. In light of (3), this assumption is equivalent to which is obviously a convex constraint as in [27]. If the small noise assumption cannot be made, e.g., in narrowband systems where TOA-based ranging measurement errors are relatively large, the constraints in (4) may not be satisfied. To avoid this limitation, we can generalize the latter as where ≥ 0 is a small number to ensure that the MN is located inside a disc with radius r i k + σ n . Note that even if the bias is zero, i.e., in LOS situation, it is more likely that the MN satisfies the constraint in (5) as compared to (4). Therefore, we propose to use the constraint in (5) throughout this work due to its robustness against measurement noise and false alarm (FA) error in NLOS identification. In the sequel, the feasible region refers to the convex feasible set formed by the intersection of the discs in (5), denoted by D k , hence At every time instant k, let us remove the NLOS measurements from the observations in (3) and only keep the LOS measurements, i.e., r i k for all i ∈ L k . The remaining LOS range measurements can be represented by the vector z k ∈ R |L k | . Note that in the worst case, where all the measurements are identified as NLOS, the vector z k is empty. The state space model and constraints can thus be expressed as where h(s k ) and n k are vectors whose entries are h i (s k ) and n i k for every i ∈ L k , respectively. Under our previous assumptions on the measurement noise n i k in (3), the covariance matrix of n k is positive-definite diagonal, i.e.
The constraints in (6c) are only on the first two elements of the state vector, i.e., x k , as we have a 2D positioning scenario herein. Note that if the constraints in (6c) are removed from the state model, then an ordinary nonlinear filtering technique such as EKF can be used. This approach is also known as EKF with outlier rejection since the NLOS measurements are regarded as outliers and therefore discarded.
In minimum mean square error (MMSE) estimation, e.g., Kalman-type filters, one tries to find the conditional mean and covariance matrix of the state vector s k given the measurements up to current time instant k, with the conditional probability density function (PDF) f (s k |z 1 , . . . , z k ). However, when extra information about the state vector is available in the form of inequality constraints, the probability that the MN is outside the feasible region should be zero. Hence a truncated or constrained conditional PDF, f c (.|.), can be defined as where β x k ∈D k f (s k |z 1 , . . . , z k )ds k is a normalization constant. Therefore, one can estimate the state vector by finding the conditional mean of s k with truncated PDF aŝ and the covariance matrix of the constrained state estimate can be found througĥ This idea is known as PDF truncation, where the distribution of the state vector given the measurements is forced to be zero outside the feasible region [28]. For a linear dynamic model with zero-mean Gaussian measurement and process noises, where the state vector is subject to linear inequality constraints, closed form expressions forŝ k andΣ k in (8)- (9) have been obtained using PDF truncation along with the Gaussian assumption [29]. For nonlinear inequality constraints, it is proposed to do a Taylor series linearization of the constraints around the current state estimate and then apply the aforementioned method; however, this approach may not be accurate [30]. In general cases with nonlinear inequality constraints, PDF truncation requires multidimensional Monte Carlo integration which becomes computationally expensive as the size of the state vector grows. Therefore, PDF truncation may not be a good approach to solve our problem. In the following section, we show how we can efficiently approximateŝ k andΣ k using an alternative approach that combines the SKURF [21] for unconstrained problems with the projection-based constrained UKF in [22].

III. CONSTRAINED NONLINEAR FILTERING WITH SIGMA POINT PROJECTION
Another family of methods for imposing inequality constraints on the state vector are the projection-based techniques, in which the unconstrained state estimate, obtained through a Kalman-type filter, is projected onto the feasible region by solving an optimization problem [28]. However, by this approach, one cannot estimate the constrained error covariance matrix of the state, i.e.,Σ k , accurately. Therefore, in addition to the unconstrained state estimate, some representative sample points of the conditional PDF f (s k |z 1 , . . . , z k ) need to be projected onto the feasible region. For instance, the sigma points of the unscented transformation (UT) can give good statistical information about the mean and the error covariance matrix of the state estimate [31]. Based on this idea, in [22], a constrained UKF technique has been proposed in which the sigma points of the UKF violating the constraints are projected onto the feasible region. However, due to the dependence of the projection function on the inverse of the a posteriori error covariance matrix, the method in [22] may become numerically unstable [30]. In the following subsections, we first describe a variation of the SRUKF that is better suited to our specific problem; then, to overcome the above mentioned numerical issue, we design a more efficient and numerically reliable method for projecting the sigma points generated from the a posteriori estimates, onto the feasible region; finally, we summarize our algorithm and comment on its numerical complexity.

A. Unconstrained SRUKF Algorithm
The proposed algorithm in this part is based on the SRUKF presented in [21] with slight modification such that the algorithm is more efficient and numerically reliable. Let s k−1|k−1 be the estimated state and Σ k−1|k−1 be the estimated error covariance matrix of the state, based on the available measurements up to current time instant k − 1. Let U k−1|k−1 be the upper triangular Cholesky factor of Σ k−1|k−1 , i.e., Then, for the next time instant, the a priori estimate of the state vector and the corresponding error covariance matrix, denoted as s k|k−1 and Σ k|k−1 , respectively, can be obtained through prediction as Alternatively, the computation of (11) can be avoided as only the Cholesky factor of the a priori covariance matrix, denoted by U k|k−1 is required [21]. To this aim, let us rewrite (11) as If we compute the QR factorization of the second matrix on the right hand side of (12), we obtain U k|k−1 : where by definition, the function qr{.} returns the upper triangular factor of the QR factorization of its matrix argument.
With the help of U k|k−1 , the sigma points of the SRUKF are generated as proposed in [21], i.e.: where N is the dimension of the state vector (in this work, N = 4), (U T k|k−1 ) j denotes the j-th column of matrix U T k|k−1 , and η α is a tuning parameter which controls the spread of the sigma points. To better understand the geometric meaning of parameter η α , we can assume that s k|k−1 and Σ k|k−1 obtained through the proposed filter are approximately equal to the mean and covariance matrix of the conditional PDF , which is the weighted squared distance between s k and s k|k−1 . Suppose that the parameter η α in (14) is chosen such that Pr(η k ≤ η α ) = α, where 0 < α ≤ 1 represents a desired confidence level. Then, the region of R N defined by η k ≤ η α represents a confidence ellipsoid, on the boundary of which the sigma points in (14) fall. For example, if α = 0.9, the probability for s k to lie inside the ellipsoid delimited by the sigma points with the corresponding η α is 90%. If we assume that f (s k |z 1 , . . . , z k−1 ) is approximately Gaussian, then the random variable η has a Chi-square distribution with N degrees of freedom and it becomes easy to find a value for η α corresponding to a certain ellipsoid with confidence level α 2 .
The generated sigma points are transformed through the nonlinear measurement function as Then, the mean, cross-covariance matrix, and error covariance matrix of the transformed sigma points can be estimated by means of weighted sums as in [32]: where (.) means that the same argument contained in the previous parenthesis is used, R is the covariance matrix of the measurement noise n k in (6a) and the weights w (j) appearing in these expressions are defined in a similar way as in [33]: and therefore satisfy 2N j=0 w (j) = 1. If the weight w (0) in (19) is negative, it is possible that the covariance matrix obtained through (18) becomes indefinite (i.e., with negative eigenvalues). However, by choosing a sufficiently large value of α, we can guarantee that η α ≥ N ; in turn, this implies that w (0) ≥ 0 and the covariance matrix (17) then becomes positive definite. In this work, we are interested in projecting the sigma points that are far away from the mean to cover a large confidence ellipsoid, and it is therefore legitimate to consider ellipsoids with larger confidence levels, so that the above issue can be naturally avoided 3 . In our dynamic model, with N = 4 and based on the Chi-square assumption for η k , it follows that if α > 0.6, then η α > N and the positive definiteness of (17) is guaranteed.
For numerical stability, instead of forming P z k|k−1 explicitly, its Cholesky factor is calculated [21]. Specifically, if we let 2N, (20) then the upper triangular Cholesky factor of P z k|k−1 , denoted by U z k is obtained through It is proposed in [21] to first compute the Kalman gain and then, the a posteriori state estimate and the Cholesky factor of the error covariance matrix can be updated through where cholupdate{U k|k−1 , K k U T z k , −1} is the consecutive downdates of the Cholesky factor of U T k|k−1 U k|k−1 using the Herein, however, we propose a more efficient and numerically reliable way to compute s k|k and U k|k . Instead of the Kalman gain K k , we compute which can be obtained by solving multiple triangular linear systems T k U z,k = Σ s,z k|k−1 . Then, it follows from (22) that K k = T k U −T z k . Substituting this expression into (23) we obtain where the vector y k U −T z k (z k −ẑ k|k−1 ) can be obtained by solving the triangular linear system 3 In [31], a scaled version of unscented transformation has been proposed to capture higher moments of the nonlinear function where the generated sigma points are located in the vicinity of each other. This method also guarantees positive definiteness of the covariance matrix. However, our problem is not highly nonlinear and we are interested to generate sigma points that might be far away from one another, therefore, our parameter selection is different from [31] and [21]. 4 In Matlab, the built-in function Cholupdate can be employed to do rank-1 Cholesky update or downdate, indicated by the third argument of the function.
From (24) and (25), it follows that the covariance matrix can be updated as hence the Cholesky factor of Σ k|k can be computed as Compared to the algorithm in [21], this modified algorithm for the estimation of s k|k and U k|k saves about 2N |N k | 2 flops at each time step k. It is also more numerically reliable as it avoids solving some linear systems, which could be illconditioned, and computing some matrix-matrix multiplications.
Note that if all the measurements at time instant k are in NLOS, then the measurement vector z k is empty, hence we will use the predicted state in (10) and the Cholesky factor of the predicted covariance matrix in (13) to replace the a posteriori state vector in (26) and Cholesky factor of the error covariance matrix in (28), respectively.

B. Imposing the Constraints on the Estimates
Up to this point, the a posteriori state estimate and the Cholesky factor of the a posteriori error covariance matrix have been obtained using a SRUKF without taking the constraints (6c) into account. To impose the constraints on the estimated state and error covariance matrix, similar to [22], a new set of sigma points are generated according to The generated sigma points form an uncertainty ellipsoid with s k|k at its centre as illustrated in Fig. 1 for N = 2. After the generation of sigma points s (j) k|k with desired confidence ellipsoid, those which violate the constraints are projected onto the convex feasible region through where W k is a symmetric positive definite (SPD) weighting matrix [30], [32]. One reasonable choice is W k = Σ −1 k|k , which gives the smallest estimation error covariance matrix when a linear KF is applied to a system with linear dynamic equations and with zero-mean Gaussian observation and excitation noises [34]. The optimization problem in (30) is a quadratically constrained quadratic program (QCQP), which is convex since W k is SPD and the constraints are convex [35]. As the constraints are only on the first two elements of the state vector, it is possible to reduce the size of the QCQP problem. A conventional way to do so is as follows. Suppose that q(1 : 2) is fixed. Then we can find the optimal q(3 : N ), which is a function of q(1 : 2). By substituting the optimal q(3 : N ) into the cost function, we obtain a QCQP, which only involves the unknown q(1 : 2).
However, in the above approach, we first need to find the matrix W k through an inverse operation which is both unnecessarily costly and numerically unstable if the covariance matrix Σ k|k is ill-conditioned. To avoid these shortcomings, we propose to use an idea from [36] in order to reformulate and reduce the size of the convex QCQP problem in (30) such that it can be solved in a more numerically reliable way. Recalling that Σ k|k = U T k|k U k|k , the objective function in (30) can be expressed as (q − s To avoid the inverse operation, we define Then we have q = s (j) We partition the lower triangular matrix U T k|k as follows: where L 11 ∈ R 2×2 and L 22 ∈ R (N −2)×(N −2) are lower triangular. Then from (32) we have Using (31) and (34), we can reformulate the QCQP problem (30) as This 2D convex QCQP problem can now be solved efficiently using iterative techniques [35].
After finding the optimal u(1 : 2), we can compute the optimal q using (32) and the fact that the optimal u(3 : N ) = 0 as follows: The above approach for reducing the size of the QCQP problem (30) not only avoids a matrix inverse computation, which may cause numerical instability (see [36]), but it is also computationally efficient. This approach is even more suitable when a SRUKF is employed since the Cholesky factor U k|k of Σ k|k is readily provided in (28).
After finding the projected sigma points through (37), the mean and covariance matrix may be estimated through weighted averaging As before, instead of (39) we compute the Cholesky factor U P k|k of Σ P k|k : As observed in Fig. 1, the projected sigma points have a different mean and covariance matrix. The weighted average of the sigma points achieved through this technique lies inside the feasible region as the average of selected points in a convex feasible region must lie in it [37]. Furthermore, the covariance matrix of the error is generally reduced as the sigma points have moved closer to each other.
Finally, in the next iteration of the unconstrained SRUKF, the constrained a posteriori state estimate s P k|k and the corresponding error covariance matrix Σ P k|k replace s k|k and Σ k|k , respectively as

C. Algorithm Summary and Computational Analysis
The proposed algorithm is summarized in Algorithm 1. Prediction of s k|k−1 using (10), and U k|k−1 using (13).

12:
Estimate the a posteriori mean s k|k using (26) and Cholesky factor of a posteriori covariance matrix U k|k using (28). 13: end if 14: Generate the sigma points using (29). 15: For every sigma point whose first two elements fall outside D k solve (36) and find the projected point (37). 16: Estimate s P k|k using (38) and U P k|k using (40). 17: Replace s P k|k and U P k|k as the a posteriori estimates, i.e., (41) and (42). 18: end for Our algorithm consists of two main parts; SRUKF and projection of sigma points, which will be discussed in more details below.
The SRUKF requires less computations compared to an ordinary UKF, thus it is more efficient. The computational analysis of SRUKF is discussed in [21]. The main computationally demanding task in the SRUKF algorithm is finding the Cholesky factor of a matrix of size R N ×N , where O(N 3 ) computations are needed using the QR factorization approach.
The QCQP in (36) is a convex optimization problem, which is not NP hard [35, p.153], and can be solved in polynomial time using an extended optimization package in Matlab such as Sedumi [38]. Since u(1 : 2) ∈ R 2 , the optimization problem can be solved with low cost. The optimization in (36) has to be done for 2N + 1 sigma points at most, however, these calculations can be performed in parallel and independently of each other; hence our technique is suitable for parallel processing. The computational cost of the algorithm depends on the number of sigma points in (29) that fall inside the feasible region, as the projection operation needs not to be applied on them. By tuning the parameter α we achieve a trade-off between accuracy and computational cost. If α is chosen to be small, then it is more likely that many sigma points will fall inside the feasible region, so the optimization problem does not need to be solved for them, resulting in a lower computational cost. However, selecting a small α may degrade the localization performance as the unconstrained estimated quantities remain unchanged after applying the constraints. On the other hand, selecting a large α increases the computational cost but at the same time may result in sampling many of the non-local points, and thus the linearisation of h(s k ) might be inaccurate [31]. In our simulations, it is observed that selecting 0.65 ≤ α ≤ 0.85 can offer a reasonable trade-off in terms of accuracy and computational cost. To model the range measurement, the true distance between each RN and MN is perturbed with a zero-mean Gaussian noise. We consider two different measurement noise scenarios: large noise σ n = 100m and small noise σ n = 10m. The large noise assumption can model general applications like narrowband cellular mobile positioning, while the small noise assumption is suitable for localization applications with accurate ranging, e.g., IEEE 802.15.4.a. Note that the accuracy of UWB ranging can be improved by increasing the bandwidth of the system [39]. We also perturb some of the measurements by NLOS biases which are modelled as exponential random variables with parameter γ = 500m. We consider three different scenarios in which, out of the total four measurements, the number of LOS ones |L k | is set to be 2, 1, or 0.
For the proposed CSRUKF we consider = 3 (for the feasible region in (5)) and α = 70%, which corresponds to η α = 4.8784 by assuming a Gaussian posterior PDF. Note that for CSRUKF, all the sigma points violating the constraints are projected onto the feasible region. For solving the QCQP problem, we use the optimization toolbox Yalmip [40] and Sedumi solver [38].
In order to see if projecting all the sigma points is necessary to achieve a good result in NLOS scenarios, we first consider the common projection technique where only the a posteriori state estimate of a KF is projected onto the feasible region [34]. Therefore, s k|k obtained through the SRUKF is projected onto the feasible region, thus the new a posteriori state estimate satisfies the constraints, however, the a posteriori estimate of the covariance matrix is not changed as compared to the unconstrained case. This approach has in general a lower computational cost compared to the proposed CSRUKF algorithm since at most one projection operation needs to be done at each iteration. We denote this approach by projection Kalman filter (PKF) and for solving the optimization problem we follow the similar procedure as done for CSRUKF.
For comparison purposes, we consider the conventional techniques proposed in [8], [41], [42], in which the range measurements are processed using a KF and then the smoothed range measurements are used in an EKF where the diagonal elements of the covariance matrix corresponding to the NLOS measurements are scaled for further mitigation of NLOS bias. While these approaches differ slightly in terms of preprocessing and variance calculation, we consider the simple one in [41] denoted by smooth EKF (SEKF) with scaling factor 1.5 and assume that the NLOS identification and variance calculation are done without error.
The Cramer-Rao lower bound (CRLB) analysis in NLOS shows that if no prior statistics about the distribution of the NLOS bias is available then the optimal strategy is to discard the NLOS measurements and only use LOS ones [43]. If prior statistics are available then the NLOS measurements should also be used to achieve a lower MSE. However, this bound can only be practical if there are enough LOS measurements for unambiguous localization, hence, for low number of LOS, i.e., less than two it can not be useful. Even though the posterior Cramer-Rao bound (PCRB) on positioning RMSE has been derived approximately in [44], [45], these derivations are based on the assumption that the NLOS bias has a Gaussian distribution with known mean and variance. Evaluating this PCRB for other NLOS distributions such as exponential is even more challenging. Since in this paper, there is no information about the distribution of the NLOS biases, except that they are positive, the mentioned lower bound is still loose and can not accurately show the lowest possible error in estimating the state vector. Due to these limitations in finding a lower bound on the positioning RMSE, we consider an ideal situation where the mean and variance of the NLOS biases are known. Then we apply a conventional EKF to the complete measurement vector r k = [r 1 k , r 2 k , . . . , r M k ] T . Therefore, the mean of the bias is subtracted from the NLOS measurements, and the covariance matrix R r of the measurement error r k is scaled according to the variance of the NLOS bias. For instance if i ∈ N k then R r (i, i) = σ 2 n + σ 2 b , where σ 2 b is the variance of the NLOS bias. Although this approach, which is denoted by bias-aware EKF (BEKF), is not optimal when the mean and variance of the NLOS bias are known (due to the non-Gaussian errors), it can be regarded as a benchmark for comparison with our method.
To evaluate the performance of the algorithms in different scenarios, we perform T = 500 Monte Carlo (MC) trials for each scenario and consider different trajectories at each trial. Let x t k and x t k|k denote the true state vector and its estimate at the k-th time step of the trajectory over the t-th Monte Carlo trial, respectively. The performance metrics are the cumulative distribution function (CDF) of positioning error e k , expressed as and the root mean square error (RMSE) of position estimate at time step k, defined as where the quantities in (43) and (44) are evaluated approximately using MC trials.
In the following, we compare the effect of measurement noise and NLOS error on the performance of different techniques in each considered scenario. We assume that the initial estimate s 0|0 is normally distributed with mean equal to the true state s 0 and covariance matrix Σ 0|0 = diag([10 4 , 10 4 , 10 2 , 10 2 ]).

A. Large Measurement Noise
In the first scenario, we consider the case of a narrowband ranging application where the noise variance is relatively high, i.e., σ n = 100m is considered. The RMSE versus time step is illustrated in Fig. 2 for the scenarios |L k | = 2,1, and 0. The corresponding CDF of the positioning error is also plotted for each scenario.
As observed, when at least two LOS measurements are available, the RMSE of CSRUKF, PKF and that of the BEKF follow each other closely. The RMSE of SEKF is higher than the RMSE of the other methods because the mean values of the NLOS biases are not mitigated efficiently. The RMSE of our proposed technique is almost the same as that of an EKF with outlier rejection, however, the latter is not shown for the sake of figure visibility. Note that for two or more LOS measurements, the EKF with outlier rejection achieves nearly optimal RMSE because the non-linearity of the range function is not very high, thus the EKF with Taylor series approximation can achieve a MSE close to the PCRB. Therefore, our algorithm has almost the same performance as the optimal approaches in this scenario, which is a satisfactory result in terms of RMSE. The reason is that it is quite unlikely that the sigma points in (29) fall outside the feasible region, thus the constraints do not change the location of the sigma points significantly. For more number of LOS measurements, the same behaviour is observed.
If less than two LOS measurements are available, then the EKF with outlier rejection diverges due to the lack of measurements in the observation vector (see in [18]). The performance of the PKF is worse than the SEKF, which shows that by only projecting the state estimate, the performance may be poor in the absence of enough LOS measurements. This is because the size of the feasible region is large due to the large measurement noise variance, and the first two elements of the sigma points, i.e., the ones related to the position, lie onto the boundary of the feasible region. Note that due to the large NLOS biases, it is generally expected that the unknown position lies inside the feasible region formed by NLOS measurements rather than on its boundary. The proposed CSRUKF, however, obtains the best performance among these approaches, and its RMSE is close to BEKF. This shows that it is necessary to project all the sigma points in (29) to achieve a good performance.

B. Small Measurement Noise
In the second scenario, we consider a case where the noise variance is relatively small σ n = 10m, which can model the UWB signalling applications. The RMSE is illustrated in Fig.  (   3 for the scenarios |L k | = 2, 1, and 0, respectively. The CDF is also illustrated in Fig. 3 for the considered scenarios.
In the presence of two LOS measurements, the RMSE of CSRUKF, PKF, and BEKF are close to each other as observed in Fig. 3. These RMSEs are close to the RMSE of an EKF with outlier rejection, although, similar to before, we have not shown the latter in this figure. The RMSE of SEKF is slightly higher which is due to the same reason mentioned earlier.
When |L k | = 1 and |L k | = 0, the EKF with outlier rejection diverges while the SEKF converges to a slightly lower RMSE than the initial RMSE. The PKF has a decent performance as compared to the large noise scenario considered earlier. This is because the size of the feasible region is smaller in this case, hence, after doing the projection of the state estimate onto the boundary, there is generally a lower uncertainty in positioning. The performance of CSRUKF is better than PKF and SEKF with a noticeable margin. Surprisingly, the performance of the CSRUKF is even better than BEKF, although our algorithm does not use the mean and variance of the NLOS biases, because we assume they are not available. One of the reasons that our method performs well under small measurement noise scenario is that the size of the feasible region is small as compared to the large measurement noise, therefore, the MN position is restricted to be within a smaller area. Hence, through processing consecutive range measurements, the sigma points can move closer and closer to each other such that their weighted average is very close to the true position. This shows the superior performance of our algorithm in scenarios where the NLOS bias is much larger than the measurement noise, e.g., UWB ranging applications.

C. Robustness to Errors in NLOS Identification
In this part, we want to analyse the performance of our technique in the presence of NLOS identification errors, i.e., FA and missed-detection (MD), which are inevitable in some applications.
To see the effect of FA in NLOS identification, we assume that we have one LOS and three NLOS RNs. However, due to the FA, the LOS link is also wrongly detected as being NLOS. Therefore, CSRUKF, and PKF wrongly remove the LOS measurements from the measurement vector, however, they employ the wrongly detected measurement to impose a constraint on the state vector. Since the parameter can increase the chance that a LOS measurement also satisfies the constraint in (5), it is expected that FA does not severely degrade the performance of our proposed technique. The simulation results are shown in Fig. 4, where it is observed that our proposed algorithm CSRUKF is robust against FA error in NLOS identification and outperforms the SEKF. Note that the BEKF algorithm suffers significantly from the FA NLOS identification error.
If the NLOS links are regarded as LOS ones, i.e., in the presence of NLOS MD error, all the Kalman-type filters have to use a biased measurement in their observation vector, and thus it is not surprising that their performances are degraded. Therefore, for our algorithm to perform well in most of the times, the threshold used for NLOS identification should change such that the probability of MD becomes very small.

V. CONCLUSION
A constrained square-root unscented Kalman filter (CSRUKF) with projection technique was considered in this paper for the aim of TOA-based localization of an MN in NLOS scenarios. The NLOS measurements were removed from the measurement vector, instead they were employed to impose quadratic constraints onto the position coordinates of the mobile terminal. The sigma points of the UKF which violated the constraints were projected on the feasible region by solving a convex quadratically constrained quadratic program (QCQP). As compared to other constrained UKF techniques, we considered a square root filter and avoided computing the inverse of state covariance matrix both in the Kalman filter and in the optimization steps, thus our approach has better numerical stability and lower computational cost. Through simulations, it was shown that, our algorithm performed better than other approaches in different NLOS scenarios. In particular, the performance was excellent when a small measurement noise variance was considered, thus, our technique is even more suitable for high resolution TOA-based UWB localization. Another advantage of our technique is its robustness to false alarm error in NLOS identification. The proposed filter can be extended to the case that the information of an IMU is fused with the range measurement for more accurate mobile localization.