Elliptic and hyperbolic localizations using minimum measurement solutions

Signal Processing 167 (2020) 107273

Contents lists available at ScienceDirect

Signal Processing journal homepage: www.elsevier.com/locate/sigpro

Elliptic and hyperbolic localizations using minimum measurement solutions Sanaa S.A. Al-Samahi a,b,1, Yang Zhang a,1, K.C. Ho a,∗ a b

Electrical Engineering and Computer Science Department, University of Missouri, Columbia, MO 65211, United States Alkhwarizmi College, University of Baghdad, Baghdad, Iraq

a r t i c l e

i n f o

Article history: Received 12 January 2019 Revised 27 August 2019 Accepted 29 August 2019 Available online 2 September 2019

a b s t r a c t Localization of an object using a number of sensors is often challenged by outlier observations and solution finding. This paper derives a new algebraic positioning solution using a minimum number of measurements, and from which to develop an outlier detector and an object location estimator. Two measurements are sufficient in 2-D and three in 3-D to yield a solution if they are consistent. The derived minimum measurement solution is exact and reduces the computation to the roots of a quadratic equation. The solution derivation leads to simple criteria to ascertain if the line of positions from two measurements intersect. The intersection condition enables us to establish an outlier detector based on graph processing. By partitioning the overdetermined set of measurements first to obtain the individual minimum measurement solutions, we propose a best linear unbiased estimator to form the final location estimate. Analysis supports the proposed estimator in reaching the CRLB accuracy under Gaussian noise. A measurement partitioning scheme is developed to improve performance when the noise level becomes large. We mainly use elliptic time delay measurements for presentation, and the derived results are applicable to the hyperbolic time difference measurements as well. Both the 2-D and 3-D scenarios are considered. © 2019 Elsevier B.V. All rights reserved.

1. Introduction Determining the location of an object using the measurements from a number of spatially distributed sensors has been a subject for research over the years [1–11]. Commonly used measurements are time delay for active scenario [12] and time differences for passive environment [13]. The research interest is in part due to the fundamental nature of the problem that has a wide range of applications, as well as the non-linear nature of the problem that makes it interesting to tackle. In a practical environment, often not all the positioning measurements are good and some maybe outliers. Outliers can come from non-line of sight (NLOS) propagation [14–16], the environment or the malfunctioning of some sensors [17,18]. Outliers could cause tremendous damage to positioning accuracy. Various techniques have been proposed to mitigate the resulting detrimental effect [19–21]. The work [21] used a residual weighting scheme to de-emphasize the measurements that may come from NLOS.

∗

Corresponding author. E-mail addresses: [email protected] (S.S.A. Al-Samahi), yzk48@mail. missouri.edu (Y. Zhang), [email protected] (K.C. Ho). 1 The first and second authors have equal contributions. https://doi.org/10.1016/j.sigpro.2019.107273 0165-1684/© 2019 Elsevier B.V. All rights reserved.

In the paper [22], the authors applied the residual test to identify and discard the NLOS measurements but the performance is highly sensitive to the detection threshold that was determined experimentally. The paper [23] formulated several hypotheses based on the Maximum Likelihood approach to detect the NLOS measurements and the performance is dependent on the availability of noise power and the NLOS error statistics. Another study [24] developed an outlier rejection method for device-free localization using RSS but it only functions properly under small noise. The paper [25] applies the Random Sample Consensus (RANSAC) technique for angle of arrival (AOA) or time difference of arrival (TDOA) localization in the presence of outliers. RANSAC is a competitive approach to handle the presence of outlier measurements for robust estimation. Nevertheless it is computationally demanding and requires the knowledge of the maximum level of non-outlier measurement error. Although the measurement equations are nonlinear, exact algebraic positioning solution exists in the specific case of critically determined scenario in which the number of measurements is equal to the number of unknowns. In such a case, the pioneer work [26] illustrated the location fix is a focus of a conic formed by the measurements. Fang [27] later showed that navigation fix of an object on earth by TDOA can be reduced to the solution of a

2

S.S.A. Al-Samahi, Y. Zhang and K.C. Ho / Signal Processing 167 (2020) 107273

quadratic equation by using the station baseline plane as reference, but location fix in the 3-D space requires solving a quartic equation. Using an intermediate variable, the study in [28] developed a solution based on Spherical-Intersection that requires the root of a quadratic equation only. Nevertheless, the method is not robust and fails to produce a reasonable solution for some sensor arrangements, due to the need to invert a matrix formed directly by the sensor positions. More recently, the papers [30,29] examine the location fix and provide a solution from a statistical point of view. In practice more measurements than unknowns are used to mitigate the noise effect and increase the positioning accuracy. Some exact solutions for the overdetermined case can be found in [31–33], which formulate the optimization problem as a polynomial system that is solved by numerical computer algebra methods or polynomial continuation techniques. Nevertheless, this paper provides an exact and simple algebraic solution which is more computationally efficient and does not rely on numerical technique. This paper first derives an algebraic positioning fix for the critically determined situation using elliptic measurements [34]. Comparing to the solutions in [27,30], we reduce complexity by applying the roots of a quadratic polynomial only rather than quartic, do not introduce extraneous solutions, and provide a direct and simpler derivation with geometric interpretation. In addition, the solution is exact, more general and without the use of station baseline plane as reference. Unlike the previous work [28], the proposed solution works for arbitrary sensor arrangements without having robustness issue (except linear arrangement in the 3D scenario in which position fix is impossible). Most important, the new proposed solution enables us to deduce the conditions for the intersection of line of positions (LOPs) defined by the measurements in 2-D and 3-D for both elliptic and hyperbolic localizations, which have not appeared before in the literature. The paper [29] did a thorough analysis and provided the compatibility conditions of two TDOAs in 2-D and their results are consistent with ours. An immediate use of the intersection condition is that it enables us to detect which of the measurements are outliers by taking a clustering approach from the graph theory. The clustering technique was used in the papers [35,36], but for time of arrival (TOA) measurements with circle loci only. Here the outlier detection is applied to 2-D and 3-D for elliptic and hyperbolic localization. The resulting outlier detector has lower complexity, is less restrictive to apply and can perform better than comparative methods [21–25]. We next propose a new estimator for the common scenario of having more measurements than unknowns. The proposed estimator partitions the measurements, generates the individual critically determined solutions and combines them together using the Best Linear Unbiased Estimator (BLUE) to form the final. The estimator is algebraic and in closed-form, does not approximate the measurement equation and performs better than the existing closedform solutions. Theoretical analysis supports the proposed method in achieving the CRLB performance under Gaussian noise, before the thresholding effect [37] sets in. A partitioning scheme is also developed to increase the noise tolerance of the final solution, based on the volume of the η-confidence ellipsoid. The proposed estimator shares some similarity with the divide and conquer framework by Abel [38]. The theory from [38] requires all measurements be independent and the total of them be an integer multiple of the minimum number of measurements in order to achieve the optimum performance. The proposed estimator does not have these limitations and reaches the CRLB performance for Gaussian noise as supported by analysis and simulations. Furthermore, we develop a measurement partitioning scheme to improve performance in the high noise region, which is new and has not been considered before in the literature.

The paper uses elliptic time delay measurements [34] for presentation. All the developments and results apply to hyperbolic TDOA measurements as well, with the changes indicated wherever necessary. Both 2-D and 3-D localizations are included. The contributions of the paper include •

•

•

•

•

The minimum measurement solutions for elliptic and hyperbolic localizations that are algebraic, closed-form, robust and require the roots of a quadratic equation only; The intersection conditions for two LOPs in 2-D and in 3-D, and their application for outlier detection using the spectral graph technique; An estimator based on BLUE to combine individual minimum measurement solutions for the more common scenario of having more measurements than unknowns (over-determined situation); Performance analysis to show that the estimator by combining individual minimum measurement solutions is able to reach the CRLB performance under Gaussian noise; A scheme using the η-confidence ellipsoid to select the individual measurement solutions to combine that provides higher noise tolerance.

We first introduce the localization scenario in Section 2, and derive the new critically determined algebraic solutions and deduce the conditions for LOP intersection in Section 3. Section 4 deduces the outlier measurement detection method and evaluates the performance. Section 5 devises the positioning estimator in the overdetermined situation from the minimum measurement solutions, conducts analysis and proposes a measurement partitioning scheme to improve performance at high noise level. Section 6 presents simulations and Section 7 gives the conclusion. We shall use the common notations that bold lower and upper case letters represent column vectors and matrices. a(i) and a(i, j) are the i-th and the (i, j)-th elements of a and A. a(i: j) is a subvector containing the i-th to the j-th element of a, and A(i: j, k: l) is a submatrix constructed from the elements of A in rows i to j and columns k to l. IN is an identity matrix of size N, the subscript may be omitted if the size is clear from the context. 0 represents a matrix of zero that may not necessarily be square. diag{ (•), , (∗ )} is a block diagonal matrix with diagonal blocks (•), , (∗ ). det(•) and trace(•) are the determinant and trace operations and • is the Euclidean norm. (•)† is the pseudo-inverse of the matrix (•). (•)o denotes the true value of (•) and (• ) = (• ) − (• )o. a is the ceiling operation for the scalar a and a the floor operation. a! is the factorial of the positive integer a. 2. Localization scenario The scenario consists of one transmitter at s0 ∈ RN and M receivers at si ∈ RN for locating an object at uo ∈ RN , i = 1, 2, · · · , M, such as in UWB indoor localization [34,39]. N is the dimension of localization which is either 2 or 3. Both the transmitter and receiver positions s0 and si are known exactly. The transmitter sends out a signal which is reflected or re-transmitted by the object [34,39] and arrives at the receiver si , as illustrated by the dotted line of Fig. 1. The resulting signal propagation time (range) measurement at si is modeled as

di = uo − si  + uo − s0  + ni

(1)

for i = 1, 2, · · · , M and ni is the measurement noise. Note that we have multiplied the time with the signal propagation speed that is assumed known in (1). The propagation time can be determined by time stamping the transmitted signal or through cross-correlating the reflected signal and the direct signal from the transmitter to the receiver [34]. Each di defines an elliptic (N=2) curve or an

S.S.A. Al-Samahi, Y. Zhang and K.C. Ho / Signal Processing 167 (2020) 107273

3

Squaring both sides of (4) once more yields a quadratic equation in u, which can be expressed in a compact matrix form as

 

[ uT 1 ] A

u

s2

s0



sM

ellipsoidal surface (N=3) for the object location. The intersection among them yields the object position estimate. For simplicity, we use d = [ d1 , d2 , · · · , dM ]T to represent the collection of the measurements. The associated noise vector n is modeled as a zero mean Gaussian random vector with covariance matrix Q. In addition to elliptic measurement, a slight modification of (1) represents TDOA measurement. In this case, s0 denotes the reference sensor of TDOA and the addition becomes subtraction:

di = u − si  − u − s0  + ni



αi αTi − I ki αTi + sT0

k i αi + s 0 . k2i − s0 2

(6)

( α 1 − α2 ) T

k1 − k2



uT

1

T

= 0.

(7)

The solution is the intersection between the ellipse (6) with i=1 (or 2) and the line (7), giving two points at most as illustrated in Fig. 1. Putting the unknown u as u = [ x y ]T in (7) and expressing x in terms of y lead to

Fig. 1. Elliptic positioning geometry in 2-D view.

o



A=

It is the standard form of an ellipse for 2-D (N = 2) or an ellipsoid for 3-D (N = 3). The 2-D case requires two measurements d1 and d2 . Subtracting (4) with i=1 and that with i=2 yields the straight line

s1

o

u = 0, 1

(2)

for i = 1, 2, · · · , M. di is the TDOA between sensor pair (si , s0 ) and it forms a hyperbola (N=2) or a hyperboloid (N=3) in which the object lies. 3. Minimum measurement solution We shall derive the explicit algebraic solution for the object position, using a minimum number of measurements that can yield a finite number of solutions. The minimum number of measurements needed is equal to the number of unknowns, which is 2 for N=2 and 3 for N=3. The minimum measurement solution has been investigated by many researchers [27,28,30,29], for the TDOA case. We focus on the elliptic case instead. The proposed solution here is more general, can work with any sensor geometry and has low complexity where it requires the root of a quadratic equation. We begin with the solution derivation for the elliptic positioning, extend it for hyperbolic localization next and then deduce the intersection conditions.

 

⎡

 

u y =B 1 1

, B=⎣

α2 (2 )−α1 (2 ) α1 (1 )−α2 (1 )

1 0

k2 −k1

α1 (1 )−α2 (1 )

0 1

⎤

⎦.

(8)

Substituting (8) into (6) with i = 1 yields a quadratic equation in y,



 

y

T

1 H y

The roots are



y=

−h(1,2 )±

1

H = BT AB.

= 0,

√

− 2hh((21,,22)) ,

h(1,2 )2 −h(1,1 )h(2,2 ) h ( 1,1 )

(9)

h ( 1, 1 ) = 0 ;

,

(10)

h ( 1, 1 ) = 0 .

Putting (10) back to (8) completes the solution. If h(1, 1 ) = 0, the line and the ellipse touch each other and we have only one solution. There are two solutions when h(1, 1) = 0 and only one of them is the actual object location estimate. The solution ambiguity can be resolved by knowing the region where the object lies [2]. The 3-D case requires three measurements d1 , d2 and d3 . With d1 and d2 (7) now represents a plane. Using d1 and d3 , we have another plane defined by



( α 1 − α3 ) T

k1 − k3



uT

1

T

= 0.

(11)

Putting the unknown u in its coordinates as u = [x y z]T , we can solve x and y in terms of z from (7) and (11) and hence

 

 

u z =B , 1 1

⎡

(α1 (1 : 2 ) − α2 (1 : 2 ))T ⎣ (α1 (1 : 2 ) − α3 (1 : 2 ))T B=

−1  α2 (3 ) − α1 (3 ) α3 (3 ) − α1 (3 )

⎤

k2 − k1 k3 − k1 ⎦. (12)

I2

3.1. Elliptic positioning We shall show that the solution evaluation reduces to the intersection of an ellipse (2-D) or an ellipsoid (3-D) with a straight line, thereby requiring the roots of a quadratic equation only. Let us first express the measurement equation in a quadratic form representing an ellipse (ellipsoid). The solution obtained from the noisy measurements is represented by u such that it satisfies di = u − si  + u − s0  according to (1). Rewriting the equation as di − u − s0  = u − si  and squaring both sides give

2di ||u − s0 || = di2 + ||s0 ||2 − ||si ||2 + 2(si − s0 )T u.

(3)

Dividing both sides by 2di , which is positive and non-zero, yields

||u − s0 || = αTi u + ki

(4)

where αi and ki are known values equal to

αi

1 = ( si − s0 ) , di

1 ki = (d2 + ||s0 ||2 − ||si ||2 ). 2di i

(12) is a straight line in the 3-D space. The solution is the intersection of the line with an ellipsoid defined by any one of the three measurements. Substituting (12) into (6) (with i=1, 2 or 3) gives a quadratic equation in z,



z

 

1 H z

T

1

= 0,

H = BT AB .

(13)

Putting the roots back to (12) yields two solutions, unless h(1, 1 ) = 0 which corresponds to the plane (7) touches the ellipsoid from d3 . The one corresponding to the object position can be identified based on the expected region of interest. The object location solution is subject to error due to measurement noise. Using the first order analysis where the noise level is not significant such that the bias square is negligible compared to variance, we have from the differential of (1)

di = (ρuo,si + ρuo,s0 )T u

(14)

where ρa,b is a unit vector pointing from b to a, i.e.

(5)

ρa,b =

(a − b ) . ||a − b||

(15)

4

S.S.A. Al-Samahi, Y. Zhang and K.C. Ho / Signal Processing 167 (2020) 107273

Thus the covariance matrix of the location estimate is

cov(u ) = G G=



ρ1

−1

QG

···

−T

ρN

.

(16)

T

(17)

is a square matrix of size N (M = N for minimum measurement solution) and

ρi = ρuo,si + ρuo,s0 .

(18)

The algebraic solution is simple to compute. It is indeed the Maximum Likelihood (ML) estimate when only N measurements are available and the covariance matrix (16) is the CRLB when the measurement noise is Gaussian.

Hence there is intersection of the two LOPs from d1 and d2 if (23) is satisfied. Albeit the condition (20) is obtained for elliptic positioning, it is also applicable for 2-D hyperbolic localization. In 3-D, the intersection of two hyperboloids can be circle, ellipse or hyperbola [42]. As a result, the proposed detection method will not be effective for 3-D hyperbolic positioning. Having derived the minimum measurement solution and the intersection conditions, we shall illustrate the use of the developed results for the detection of outliers and for the localization solution of overdetermined situation. The results in this section can be used individually for either outlier detection or overdetermined solution, and we decouple the presentation of the two applications individually in the following two sections.

3.2. Hyperbolic positioning The measurement equation is (2). Going through the same derivation process, the solution expressions for hyperbolic positioning are the same as for elliptic positioning, the only difference is the sign change of αi and ki that are defined in (5), which will not affect the solution expressions. The covariance matrix of the location estimate remains to be given by (16), with the vector ρi for the matrix G now becomes

ρi = ρuo,si − ρuo,s0 .

(19)

3.3. Intersection condition The curves/surfaces defined by measurements must intersect to ensure a solution. Intersection may not occur when the noise level is large or some measurement is invalid. Revealing whether intersection occurs from two measurements would be useful to determine if solution exists or if a measurement is corrupted to some extent that should not be used. The solution derivations above can give us the condition for intersection that has not appeared in the literature. 2-D case (N=2) The y-coordinate solution (10) must be real. Hence the condition for intersection is that the quantity inside the square-root is positive, i.e.

det(H ) ≤ 0 .

(20)

H is the 2 × 2 matrix defined in (9), and B and A are given by (8) and (6). 3-D case (N=3) Let us consider the two measurements d1 and d2 . The intersection between two ellipsoids can be reduced to that of the plane (7) and one of them. Solving x in terms of y and z (7) gives

 



y u =C z , 1 1



C=

4. Outlier measurement detection The measurements acquired by a set of receivers often contain outliers; due to sensor malfunctioning, NLOS propagation, or other possible reasons. The outlier measurements can reduce considerably the localization accuracy if included in a positioning estimator. It is of practical interest to detect if a measurement is an outlier and remove it from estimation. The outlier detection problem has been considered by many researchers and a number of techniques have been proposed especially for applications that encounter NLOS propagation [21–25]. We shall propose an outlier detection method based on the intersection conditions derived in Section III and the scheme from [35] that was used for 2-D TOA measurements only. Let us be given M measurements, among them an unknown amount of Mo are outliers. The number of non-outlier measurements M − Mo is more than 2 for the 2-D case and 3 for the 3-D case, to ensure localization is possible. We are interested in determining which of the M measurements are outliers. 4.1. Approach We shall cluster the measurements into two groups, one for outliers and the other for non-outliers, such that the non-outlier group has the largest number of intersections among them. For better illustration, Fig. 2 shows an example where the number of measurements is eight and they are represented by circles. A connection between two measurements indicates that their measurement equations intersect. Each of the first six measurements has at least three intersections among them whilst the last two has only one. The clustering approach groups the measurements based on their number of intersections. 4.2. Method

(α1 (1 ) − α2 (1 ))

 −1

(α2 (2 : 3 ) − α1 (2 : 3 ))T I3

k2 − k1

 .

(21)

We perform outlier detection through the construction of an adjacency matrix with the steps shown below [35].

Using (21) in (6), the intersection of the two ellipsoids is defined by



y

z

 

1 H y

z

T

1

= 0,

H = CT AC .

(22)

It has been proven that the intersection of a plane and an ellipsoid forms an ellipse [40]. Thus ensuring intersection for two elliptical measurements requires (22) to form a real ellipse in the y-z plane. The expression (22) is a conic that defines a real ellipse if the following two conditions are met [41]:

det(H(1 : 2, 1 : 2 ) ) > 0 ; trace(H(1 : 2, 1 : 2 )) × det(H ) ≤ 0 .

(23)

Fig. 2. Illustration of the proposed outlier detection method through clustering of the measurements.

S.S.A. Al-Samahi, Y. Zhang and K.C. Ho / Signal Processing 167 (2020) 107273 •

Construct the M × M adjacency matrix  whose elements are either 0 or 1, such that

ψ (i, j ) =

1, 0,

if measurement equations of di and d j intersect otherwise

(24)

•

for i, j = 1, 2, . . . , M, where (20) for 2-D or (23) for 3-D is used to decide if intersection occurs. When there is no outlier, all elements except those in the diagonal will be 1. Compute the vector γ ∗ that maximizes the Rayleigh quotient of ,

γ ∗ = arg max γ

•

γT  γ . γT γ

(25)

The result is the eigenvector of  having the largest eigenvalue. Generate the indication vector v by comparing each element of γ ∗ with the threshold β :

v (i ) =

1, 0,

if γ ∗ (i )  β otherwise

(26)

for i = 1, 2, · · · , M. Essentially, di is an outlier measurement if v(i ) = 0. A smaller value of β makes the decision of outlier measurement more conservative. The outlier detection scheme above is based on the spectral graph theory [43]. Having the adjacency matrix constructed by (24) with elements 1 and 0 only in Step-1, the Rayleigh quotient corresponds to the average number of intersections among the measurements indicated by the indices of the non-zero elements in γ . Thus maximizing the Rayleigh quotient in Step-2 provides the subset of the measurements among them having the largest number of intersections. The resulting γ ∗ is a vector of real numbers. Step-3 converts it to a non-outlier measurement indication vector of 1 and 0 only by using a threshold β . In practice, choosing a suitable β value may be challenging. Clearly v is a function of β . A reasonable choice is to select the threshold through optimization by [35]

β = arg max β

v ( β )T γ ∗ . v ( β )T v ( β )

5

from the actual values much so that their LOP intersections remain with other measurements. The proposed outlier detection method is applicable to both elliptic and hyperbolic positionings. Nevertheless, it is not effective for the hyperbolic case in 3-D due to the large possibilities of the geometry from the intersections of two hyperboloids. 4.3. Multiple objects scenario It may be possible to apply the proposed outlier detection method to the localization scenario of having multiple objects in the presence of outliers. Among the total M measurements, Mi will come from object i and the total number of inlier measurements is K i=1 Mi , where K is the number of objects. The number of outliers  is Mo so that M = Ki=1 Mi + Mo. We shall interpret a little differently when making use of the proposed algorithm in Section IV.B for this scenario to detect the outliers. Let us consider without loss of generality that object 1 has the largest number of inlier measurements, followed by object 2, etc., until object K. With respect to object 1 the rest of the M − M1 measurements are outliers. Using the proposed algorithm will be able to identify the M1 inlier measurements of object 1 and they can be taken out from the total. Applying the proposed algorithm again on the remaining will separate the M2 inlier measurements for object 2. The process continues and we will be able to identify the inlier measurements for all the objects, leaving the outliers at the end. In the situation when the number of objects K is not known, the process can stop if the maximum eigenvalue of  from the Rayleigh quotient evaluation is below certain value. Another possibility is to apply the proposed algorithm is to identify the outliers and the inlier measurements for all K objects at once based on the K largest eigenvalues of  . Outlier detection together with measurement association in the multiple objects scenario is an interesting subject for our future study. Detailed investigation will be performed on the effectiveness of the proposed algorithm to handle this scenario. We shall present the results in the near future when the investigation is complete. 4.4. Simulations

(27)

Setting β from (27) will give an indication vector that has the largest correlation coefficient with γ ∗ . Once the indication vector is obtained, only the non-outlier measurements will be passed to an estimator for obtaining the object location. The existing outlier detection or mitigation techniques from the literature require certain approximations, several iterations or residual error evaluations. We use the approach of measurement equation intersections and adjacency matrix formation. It is computationally efficient and applicable under various environments with little knowledge about the outliers required. The proposed method does not require an outlier not intersecting with any inlier for it to be detected, and it does not require an inlier intersecting with all other inliers (non-outliers). This is illustrated in the example shown in Fig. 2 where the outlier represented by node 7 has intersection with the inlier node 5 and inlier node 5 has no intersection with inlier node 2 and inlier node 1. As long as an outlier has less number of intersections than those among the inliers, the proposed method is able to detect it. This is because the proposed method considers all pairwise intersections in forming the adjacency matrix and maximizes the Rayleigh quotient to collect the measurements together with maximum number of connections among them to form the inlier group. We would like to clarify that the proposed method may not be able to detect NLOS measurements, when they are not deviating

We follow [14,16] to model the outlier scenario. Each of the M measurements has a probability of being an outlier. An outlier will have significant bias in the mean and higher in the noise power. In other words, the probability density of measurement noise is [14,16]

p( n i ) = ( 1 − ε ) N ( 0 , σ 2 ) + ε N ( μ o ,

σo2

σ 2.

σo2 ) , i = 1, 2, · · · , M.

(28)

where |μo | 0 and > The parameter settings for the simulation are as follows: M = 7, ε = 3/7 (3 outliers), μo = (1 + rand ) × true measurement values with rand sampled from U (0, 1 ) for the elliptic and rand sampled from U (−1, 1 ) for the hyperbolic case, σo = 2.18 × σ , and the measurement noise is uncorrelated for elliptic and correlated with a correlation coefficient of 0.5 for hyperbolic cases. The object, transmitter and sensor positions are tabulated in Tables 1 and 2. The sensor positions for the 2-D case is from [22] and those for the 3-D case are random placement. The number of ensemble runs is L = 10 0 0 and in each run which sensor measurements are outliers are assigned randomly. The proposed method in the simulation has two components. The first is to apply the outlier detection method in Section 4, where the detection threshold β is chosen using (27) and it is different for different ensemble runs. The decided non-outlier measurements are next passed to the MLE with Gauss-Newton implementation [37,45] to obtain the location estimate, where it is initialized at the true object position corrupted by uncorrelated zero-

6

S.S.A. Al-Samahi, Y. Zhang and K.C. Ho / Signal Processing 167 (2020) 107273 Table 1 The object, transmitter and receiver positions for the 2-D case. uo

s0

s1

s2

s3

s4

s5

s6

s7

8000 8000

−6000 −6000

4000 6000

−4000 6000

−6000 −1000

−3000 −5000

3000 −6000

6000 0

0 5000

Table 2 The object, transmitter and receiver positions for the 3-D case. uo

s0

s1

s2

s3

s4

s5

s6

s7

1400 1400 1400

0 0 0

641.83 −244.38 577.12

372.39 −752.76 −996.86

981.56 607.20 −886.61

−339.52 −293.63 −363.99

−903.42 −379.67 −538.14

−714.57 891.23 279.46

−392.72 833.71 −202.85

mean Gaussian noise with variance equal to twice the CRLBs of the coordinates. The estimation performance in terms of the root mean-square error (RMSE) is compared with the following methods from the literature: (i) Residual Weighting algorithm (Rwgh) [21], (ii) Residual Test (RT) [22], (iii) approximate Maximum Likelihood detection algorithm (ML) [23], (iv) Intersection of Link approach with Outlier Intersection Rejection (ILOIR) [24], and (v) RANSAC [25]. The RANSAC tailored for localization in [25] requires the number of non-outlier measurements to be at least 5 in 2-D and 6 in 3-D localization. For the RANSAC results presented, the algorithm in [25] is used to detect the outliers and the iterative MLE is followed to obtain the best possible final location estimate. The maximum level of measurement error required for RANSAC is set at 3 times the square-root of the noise power, which covers 99.75% of the possible Gaussian measurement noise values. To gain better understanding of the detection performance of the proposed method, we also provide the percentage of correct detection (PCD ) and the percentage of false alarms (PFD ) defined as

L PCD =

l=1 all outliers detected in run l ×100% total number of ensemble runs L

L PF D =

l=1

(29)

at least one falsely detected outlier in run l ×100%. total number of ensemble runs L (30)

When the noise level is high, it is possible that some non-outlier measurements will be falsely detected as outliers. Nevertheless, the localization performance may still be acceptable. Being conservative of having higher PFD while maintaining 100% PCD by increasing the threshold β may be preferable, since the performance degradation is often much more severe when including outliers instead of removing non-outlier measurements for localization. The localization accuracy in the 2-D scenario with three outliers (ε = 3/7) is shown in Fig. 3. The performance without detecting and removing the outliers is very poor. The proposed method is able to mitigate the outlier measurements and reaches the CRLB. It outperforms the competing methods with ML following closest. At the noise level beyond 103 , the proposed method deviates from the bound while ML has the best result. The lower plot illustrates the detection performance of the proposed method, and the RT and ML methods. The other methods do not handle outliers by detecting and removing them. The correct detection rates correlate well with the localization performance. RT has an inherent large false detection rate even when the noise level is small (see [22] for explanation). This happens to ML as well, due to the unavailability of outlier occurrence probability that can cause the method favoring the hypothesis with less non-outlier measurements [23]. Nevertheless, the false detections of RT and ML do not affect the localization accuracy much, as the falsely detected outlier measurements have little effect on the location estimate.

Table 3 Relative computation times for the proposed, Rwgh, RT, ML, ILOIR and RANSAC methods for outliers detection only. Algorithm

Proposed

Rwgh

RT

ML

ILOIR

RANSAC

Time/2-D Time/3-D

1 1

6.37 7.51

50.09 62.98

8.61 10.94

2.06 2.41

18.74 21.30

Fig. 4 illustrates the performance for 2-D localization with two outliers (ε = 2/7), where the RANSAC results are included. The proposed algorithm remains to perform well until the noise level σ 2 reaches 10,0 0 0, and the observations are consistent with those in Fig. 3. RANSAC provides the CRLB localization accuracy before the noise level approaches 1,0 0 0. Beyond the percentage of correct detection drops below 100%, leading to large localization error. RANSAC has non-zero percentage of false detection even when the noise level is very low, and such a behavior for RANSAC was also indicated in [25]. When the noise level is very high, RANSAC is unable to detect outliers and classifies all measurements as nonoutliners, leading to a decrease in the percentage of false detection. Fig. 5 gives the performance for the 3-D scenario with the challenging situation of three outliers. The proposed method remains to function well. The proposed method is better than ML in the high noise region and deviates from the CRLB later. The localization accuracy relates well with outlier detection performance in the lower plot. Fig. 6 gives the results for 3-D localization with one outlier (ε = 1/7) in which the RANSAC results are illustrated for comparison. The observations are similar to those of the 2-D case in Fig. 4. The proposed method remains to outperform the others and it can tolerate larger noise level in giving 100% correct detection than RANSAC. For completeness, Fig. 7 illustrates the results for hyperbolic positioning in the 2-D scenario with three outlier measurements, where the proposed method stays close to the CRLB until the noise power reaches 10. The proposed method performs as expected and gives consistent behaviors with the elliptic positioning case. Fig. 8 repeated the simulation as in Fig. 5 with the distribution of the noise in (28) changed to Laplacian (non-Gaussian). The proposed method remains to behave properly. Rwgh and ILOIR have lower RMSE values than the other methods at very high noise level due to their considerable amounts of bias. Nevertheless, the error at such high noise level is so large that renders the usefulness of their solutions. We have tested elliptic and hyperbolic localizations with one or two outliers, for Gaussian or Laplacian noise. The observations are consistent with the 3 outlier cases. In addition to performance, the computation advantage of the proposed method is shown in Table 3. It tabulates the relative computation time needed among the methods, obtained by Matlab implementation running on a PC with an i7 processor. The values

S.S.A. Al-Samahi, Y. Zhang and K.C. Ho / Signal Processing 167 (2020) 107273

7

106

4

10

2

RMSE

10

CRLB Without Detection Proposed Method Rwgh RT ML ILOIR

100 -10

0

10

20

30

40

50

60

70

40

50

60

70

2

10log(σ )

Detection Rate (%)

100 80

CD_Proposed Method FD_Proposed Method CD_RT FD_RT CD_ML FD_ML

60 40 20 0 -10

0

10

20

30 2

10log(σ ) Fig. 3. Performance for 2-D elliptic positioning in the presence of three outlier measurements, Gaussian noise.

10 6 CRLB Without Detection 10

Proposed Method

4

RMSE

RANSAC MLE Rwgh RT 10 2

ML ILOIR

10 0 -10

0

10

20

30

40

50

60

70

40

50

60

70

10log(σ 2)

Detection Rate(%)

100 CD_Proposed Method

80

FD_Proposed Method CD_RANSAC

60 40 20

FD_RANSAC CD_RT FD_RT CD_ML FD_ML

0 -10

0

10

20

30 2

10log(σ ) Fig. 4. Performance for 2-D elliptic positioning in the presence of two outlier measurements, Gaussian noise.

8

S.S.A. Al-Samahi, Y. Zhang and K.C. Ho / Signal Processing 167 (2020) 107273

10

6

CRLB Without Detection Proposed Method Rwgh RT ML ILOIR

RMSE

104

102

10

0

-10

0

10

20

30

40

50

60

70

40

50

60

70

2

10log(σ )

Detection Rate (%)

100 80

CD_Proposed Method FD_Proposed Method CD_RT FD_RT CD_ML FD_ML

60 40 20 0 -10

0

10

20

30 2

10log(σ ) Fig. 5. Performance for 3-D elliptic positioning in the presence of three outlier measurements, Gaussian noise.

10 6 CRLB Without Detection Proposed Method

10 4

RMSE

RANSAC MLE Rwgh RT

10

ML

2

ILOIR

10 0 -10

0

10

20

30

40

50

60

70

40

50

60

70

10log(σ 2)

Detection Rate(%)

100 CD_Proposed Method

80

FD_Proposed Method CD_RANSAC

60 40

FD_RANSAC FD_RT CD_RT

20

CD_ML FD_ML

0 -10

0

10

20

30

10log(σ 2) Fig. 6. Performance for 3-D elliptic positioning in the presence of one outlier measurement, Gaussian noise.

S.S.A. Al-Samahi, Y. Zhang and K.C. Ho / Signal Processing 167 (2020) 107273

RMSE

10

9

10

105

CRLB Without Detection Proposed Method Rwgh RT ML ILOIR

100

-30

-20

-10

0

10

20

30

10

20

30

2

10log(σ )

Detection Rate (%)

100 80 60 40 20

CD_Proposed Method FD_Proposed Method CD_RT FD_RT CD_ML FD_ML

0 -30

-20

-10

0 2

10log(σ ) Fig. 7. Performance for 2-D hyperbolic positioning in the presence of three outlier measurements, Gaussian noise.

for Rwgh, RT, ML and ILOIR are from the 3 outlier case (Figs. 3 and 5) and those for RANSAC are from the 2 outlier case in 2-D (Fig. 4) and the 1 outlier case in 3-D (Fig. 6). The proposed method is most efficient, which is about 50 times, 9 times and 2 times better than RT, ML and ILOIR. It is more than 18 times faster than RANSAC. The proposed detector favors not a large ratio of the number of outliers to the number of total measurements. In the simulations above, the detector functions well when the ratio is not larger than 3/7=0.43. 5. Overdetermined solution We shall consider in this section the derivation of the overdetermined solution from the minimum measurement solution. To simplify the exposition, we shall assume the outliers have been detected by the algorithm in the previous section if outliers exist, so that the measurements in obtaining the overdetermined solution presented here are free of outliers. Such a setting also enables us to validate the proposed overdetermined solution is able to reach the CRLB accuracy. In practice, a localization system often has sensor measurements more than the number of unknowns to mitigate the additive noise and increase the positioning accuracy. Having more than the minimum number of measurements also eliminates the solution ambiguity. Most of the solutions proposed over the years in the literature are for such a situation. Some of them are it-

erative such as the iterative implementation of the MLE that requires good initial guesses close to the actual solution [37,45]. Others are explicit closed-form solutions that are applicable when the noise level is not significant [2,11,39], or are convex optimization solutions with semi-definite relaxation that may increase solution bias [6,46]. We shall develop a new solution by combining the individual minimum measurement algebraic solutions developed in Section 3. The proposed solution does not require iteration and can tolerate a larger amount of noise than the existing closed-form solutions from the literature. Under Gaussian noise, we are able to show that the proposed estimator can achieve the CRLB performance. In the following, we shall first identify the correct minimum measurement solution by eliminating the solution ambiguity. The proposed solution is next developed through the BLUE [37] approach. Analysis is performed to show the proposed method is able to give the CRLB performance. Finally, a method to select the individual solutions for BLUE is proposed that can yield better performance when the noise level is large. 5.1. Eliminating minimum solution ambiguity There are two solutions obtained from a minimum number of measurements as shown in Section III. Both are exact but only one corresponds to the object location. In the absence of prior knowledge of the region where the object lies, we shall use a statistical

10

S.S.A. Al-Samahi, Y. Zhang and K.C. Ho / Signal Processing 167 (2020) 107273

RMSE

10

6

104

10

Without Detection Proposed Method Rwgh RT ML ILOIR

2

-10

0

10

20

30

40

50

60

70

40

50

60

70

2

10log(σ )

Detection Rate (%)

100 CD_Proposed Method FD_Proposed Method CD_RT FD_RT CD_ML FD_ML

80 60 40 20 0 -10

0

10

20

30 2

10log(σ ) Fig. 8. Performance for 3-D elliptic positioning in the presence of three outlier measurements, Laplacian noise.

approach based on the residual error obtained from the rest of the measurements to decide which of the two relates to the object location. Let us be given M measurements, denoted by the vector d. Among the M measurements, we select N and use them to produce the two minimum measurement solutions p1 and p2 . We take those N measurements used in finding p1 and p2 out from d and ˜ , which has a length of M − N. Also, denote the rest by the vector d let  Q be the noise covariance matrix of  d, which is essentially Q with the rows and columns corresponding to the N measurements removed. For each of the two solutions, we generate the residual square error

ξ j = ( d − d(p j ))T  Q−1 ( d − d(p j )) , j = 1, 2

(31)

where  d(p j ) is the reconstructed measurement values using pj . If pj is the solution corresponding to the object, ignoring the estimation error, ξ j will follow the central χ 2 distribution with M − N 2 [37]. Given a certain tail degrees of freedom [44], i.e. ξ j ∼ χM−N probability for incorrect decision, say 0.01, we can determine the corresponding threshold T . The proper solution is p j∗ if ξ j∗ ≤ T , j ∗ = 1, 2.

ity, using measurements d(((i−1 )N+1 )modM ) , , d((iN)modM) , where d0 is the same as dM . The modulo operation is to account for the situation that the number of measurements M is not an integer multiple of the localization dimension N where the first few measurements are used twice. For instance, if M = 7 and N = 3, u(2) is the 3-D positioning solution obtained by using the measurements d4 , d5 and d6 , and u(3) is the solution from d7 , d1 and d2 . According to (14) and (16), we can model u(i) as

u(i ) = uo + (G(i ) )−1 n(i )

where is the matrix defined in (17), with the sensor indexes {1, 2, , N} replaced by {((i − 1 )N + 1 )modM, · · · , (iN )modM} and n(i) is the N × 1 noise vector of the measurements used to determine u(i) . We shall obtain the final estimate from (32) using the BLUE [37]. Since n(i) is Gaussian distributed, the BLUE will correspond to the MLE under the linear model (32) [37]. The proposed method consists of the following three steps: •

•

5.2. Proposed final solution Let us first introduce the notation u(i) to denote the minimum measurement solution, also called individual solution for simplic-

(32)

G(i)

•

divide M measurements into L = M/N  sets each has N elements; generate the minimum measurement solution from each set to obtain u(i) , i = 1, 2, · · · , L, where the solution ambiguity is resolved using the procedure described in the previous subsection; combine these individual solutions by BLUE to produce the final estimate.

S.S.A. Al-Samahi, Y. Zhang and K.C. Ho / Signal Processing 167 (2020) 107273

When the noise is not significant, it is irrelevant of how the measurements are divided in the first step. A smart way to partition the measurements, however, can improve the accuracy in the large error region, which will be developed later in this Section. The second step follows the method in Section III. We shall elaborate on the third step in more details. After obtaining the individual solutions u(1) , , u(L) , the matrix equation constructed based on (32) is

h = Huo + B,

(33)



h = u ( 1 )T , u ( 2 )T , · · · , u

 ( L )T T

and B is the LN × LN matrix



B = diag G(1) , G(2) , · · · , G(L )

H = [ IN , IN , · · · , IN ]T ,

,

−1

.

(34)

(35)

The noise vector is

 T  = n(1)T · · · , n(L)T = Cn

(36)

and C is an LN × M matrix given by



C=



IM

ILN−M 0

.

(37)

Applying the Gauss-Markov Theorem [37], the BLUE for uo from (33) is

u = (HT WH )−1 HT Wh,

(38)

W = E[BT BT ]† = B−T C†T Q−1 C† B−1 ,

(39)

where (36) has been used. We use pseudo inverse since E[T ] is singular when M < LN. Appendix A shows that C† is an M × LN matrix equal to



0.5ILN−M C = 0 †

0 I2M−LN



0.5ILN−M . 0

(40)

The diagonal blocks of B contain the true object location that is  not known. We shall replace uo in G(i) by the average Li=1 u(i ) /L.

Using (45) in (43), comparing with (41) validates immediately that the proposed solution is able to achieve the CRLB accuracy:

cov(u ) = CRLB(uo ) .

(46)

We would like to emphasize that (46) is obtained under the condition that the noise level is not exceedingly large so that re placing the true value uo by the average Li=1 u(i ) /L in G(i) for the proposed method has negligible error. 5.4. Individual solution selection The performance analysis shows that as long as the proposed method combines a minimum set of individual solutions that covers all the measurements, the accuracy can reach the CRLB under Gaussian noise. This is the case when the measurement noise level is not significant. Otherwise, the performance will move away from the bound. This is a general phenomenon for a non-linear estimation problem, which is apparent in Figs. 3 to 7. Different groupings of the measurements give different individual solutions for BLUE in the proposed method, which can lead to different behaviors in the large error region. We shall propose a grouping such that the final estimate has higher noise tolerance that the performance will deviate from the CRLB later as the noise power increases. In most cases, we observe that the final solution with the proposed grouping has comparable noise threshold with MLE having performance deviating from the CRLB. The measurement grouping problem with least amount of measurement replication is equivalent to selecting L = M/N  individual solutions to combine, where the union of the measurements that generate them contains all M measurements. The motivation for individual solution selection comes from the fact that the quality of the individual solutions is different, depending on the relative geometry between the object and the sensors for an individual solution. When selecting the better quality individual solution to combine using BLUE, the final solution can have a higher level of noise tolerance in reaching the performance bound. In essence, given M measurements for localization in the N dimensional space, there is a total of CNM individual solutions



5.3. Performance

11



S = ui : i = 1, 2, · · · , CNM .

(47)

When the error is negligible in approximating G(i) by using the  average Li=1 u(i ) /L for uo in forming B, the analysis below shows that the proposed solution attains the CRLB performance under Gaussian noise. The CRLB for the object location with elliptic measurements in Gaussian noise is [11]

We are interested in selecting L = M/N  elements from S for BLUE to generate the final solution. The set of all possible combinations having the L individual solutions produced by all the measurements is denoted by

CRLB(uo ) = [ρ1 ,

Appendix B shows that J is equal to



ρ2 , · · · , ρM ]Q−1 [ρ1 , ρ2 , · · · , ρM ]T

−1

(41)

where ρi is defined in (18). To obtain the covariance matrix of the proposed solution, we subtract both sides of (38) by uo , multiply with its transpose and take expectation, giving



cov(u ) = HT WH

−1

.

(42)

It becomes after substituting the second line of (39)



cov(u ) = HT B−T C†T Q−1 C† B−1 H

−1

.

(43)

From direct evaluation using (34) and (35),



HT B−T = G(1)T , G(2)T , · · · , G(L )T





I 0 = [ρ1 , ρ2 , · · · , ρM ] LN−M 0 I2M−LN



ILN−M . 0

(44)

ρ2 , · · · , ρM ] .



J=

( M − L )! C M C LN−N . (N! )L−1 (M − (L − 1 )N )! L LN−M

(48)

(49)

For instance assuming M = 5. For 2-D positioning (N = 2), the number of individual solutions is CNM = 10, the minimum number of individual solutions to cover all measurements is L = M/N = 3 and the number of possible choices of individual solutions for BLUE is J = 20. In the case of 3-D, the values become 10, 2 and 15. We would like to choose the combination in C such that the final solution has high noise tolerance. The quality of an individual solution can be characterized by the η-confidence ellipsoid [47]. It is the minimum volume ellipsoid centered at the estimate u that contains uo with a confidence level of η, which is equal to

  ξη = v | vT  −1 v ≤ Fχ−12 (η )

(50)

N

Substituting (40) yields

HT B−T C†T = [ ρ1 ,



C = c j = ( j1 , j2 , · · · , jL ) : j = 1, 2, · · · , J .

(45)

where Fχ 2 is the cumulative distribution function of a chi-squared N

random variable with N degrees of freedom.  is the covariance

12

S.S.A. Al-Samahi, Y. Zhang and K.C. Ho / Signal Processing 167 (2020) 107273

40

15 CRLB Proposed Solution Fang Method [27] SX Solution [28]

10 5

CRLB Proposed Solution Proposed Solution (NoSelection) MLE Solution from [11]

30 20

10 log(MSE)

10 log(MSE)

0 -5 -10 -15

10 0 -10

-20

-20

-25 -30 -40

-35

-30

-25

-20

-15

-10

-5

-30 -20

0

-10

0

10 log(σ ) Fig. 9. Performance of the proposed minimum measurement solution in 3-D hyperbolic localization.

40 CRLB Proposed Solution Proposed Solution (NoSelection) MLE Solution from [11]

30

CRLB Proposed Solution Proposed Solution (NoSelection) MLE Solution from [11]

30

10 log(MSE)

10 log(MSE)

30

20

Fig. 11. Performance for 2-D elliptic positioning with 5 measurements.

50 40

10

10 log(σ2 )

2

20 10 0

20 10 0 -10

-10

-20

-20

-30 -40

-30 -20

-10

0

10

20

10 log(σ2 )

matrix of the estimate u. A scalar measure for the solution quality is the volume of the η-confidence ellipsoid given by [47]

vol(ξη ) =

Fχ−1 2 ( η )π

N/2

N

(N/3 )



det  1/2



•

(51)

For each individual solution ui ∈ S, compute the determinant of its covariance matrix defined in (16) and call it Ri ; For each combination choice c j ∈ C, obtain the quality factor

Qj =

L  l=1

R jl ;

-10

0

10

20

10 log(σ2 )

•

The combination choice for the individual solutions to apply the proposed estimator is

c∗ = arg min Q j . c j ∈C

where  (•) is the Gamma function. The determinant of a square matrix is equal to the product of its eigenvalues. Thus the volume is proportional to the determinant of the covariance matrix of the location estimate. A logical approach to improve the performance at higher noise level is to select the individual solutions with smaller ellipsoid volumes (smaller determinants of the covariances) for BLUE. The procedure to select the individual solutions for BLUE is as follows: •

-20

Fig. 12. Performance for 2-D elliptic positioning, averaged over 20 randomly generated localization geometries.

Fig. 10. Performance for 2-D elliptic positioning with 3 measurements.



-30

30

(52)

(53)

In the first step, the uo needed in (18) for the covariance expression will be replaced by the individual solution. In the second step, we use product to emphasize that the ellipsoid volumes of all the individual solutions for BLUE are considered together. Using addition instead of multiplication should also work. The proposed individual solution selection scheme is more suitable when the number of measurements M is not large. Otherwise, the complexity may not be easy to manage as the number of combinations J can increase considerably. The individual solution selection scheme presented in this section is needed only when the measurement noise level is high for the purpose of extending the noise threshold where the estimation accuracy starts deviating from the CRLB. When the noise level is low, it is not necessary to apply this scheme. The minimum solution selection algorithm is summarized below. The solution and analysis presented in this section assume the measurements do not have outliers. When outliers are present, we use the outlier detection method presented in Section IV to remove them before applying the overdetermined solution derived in this

S.S.A. Al-Samahi, Y. Zhang and K.C. Ho / Signal Processing 167 (2020) 107273

60

60 CRLB Proposed Solution Proposed Solution (NoSelection) MLE Solution from [11]

50

50 40

30

10 log(MSE)

10 log(MSE)

40

20 10

30

10 0

-10

-10

-20

-20 -30

-20

-10

0

10

-30 -40

20

CRLB Proposed Solution Proposed Solution (NoSelection) MLE Chan-Ho Solution [2] SCWLS Solution [48]

20

0

-30 -40

-30

-20

-10

0

10

10 log(σ2 )

10 log(σ2 ) Fig. 13. Performance for 3-D elliptic positioning, averaged over 20 randomly generated localization geometries.

Fig. 15. Performance for 2-D hyperbolic positioning, averaged over 20 randomly generated localization geometries.

100

50 CRLB Proposed Solution Proposed Solution (NoSelection) MLE Chan-Ho Solution [2] SCWLS Solution [48]

30 20

80 60

10 log(MSE)

40

10 log(MSE)

13

10 0

CRLB Proposed Solution Proposed Solution (NoSelection) MLE Chan-Ho Solution [2] SCWLS Solution [48]

40 20

-10 -20

0

-30

-20

-40 -30

-25

-20

-15

-10

-5

0

5

10

-40

-30

Fig. 14. Performance for 2-D hyperbolic positioning with the object near the center of the sensors.

section. The proposed overdetermined solution will be degraded if there is some miss in the detection of outliers. The degradation is not coming from the overdetermined solution but rather from the accuracy of the detection method. The degradation will follow the Gauss-Newton MLE localization performance as indicated in Figs. 3–8 when outlier detection fails at high noise level. 6. Simulations We shall present localization performance of the proposed solution for elliptic as well as hyperbolic positionings. The results will be for a single configuration and for the average of 20 configurations. In the latter case, the geometries are generated by placing the sensors and the object at random locations. The number of ensemble runs in each geometry is 20 0 0. For the proposed estimator, no prior knowledge about the region where the object lies is assumed. We use the procedure described at the end of Section V.A to decide which of the two individual solutions corresponds to the location of interest. The threshold T is set to 6.635. This value corresponds to a tail (false detection) probability of 0.01 from the central χ 2 distribution, when the number of measurements is M = N + 1, where N is the dimension of localization. The results from MLE are provided for comparison. It is implemented by the Gauss-Newton iteration with at

-20

-10

0

10

10 log(σ2 )

10 log(σ2 )

Fig. 16. Performance for 3-D hyperbolic positioning, averaged over 20 randomly generated localization geometries. Table 4 Computational complexity comparison for minimum measurement solutions.

Relative Computation Time / 3-D

Proposed

FANG[27]

SX Solution [28]

1

1.48

2.48

least 30 random initializations, where the final MLE solution is the one having the largest likelihood. It is quite computationally intensive albeit expected to yield the CRLB performance. Also included are the closed-form solution from [11] for elliptic, and [2] and [48] for hyperbolic in the literatures. The noise covariance matrix is σ 2 IM for elliptic and σ 2 (IM + 1M 1TM )/2 for hyperbolic positionings [2]. The performance is shown in terms of the Mean-Square Error (MSE) of the uo estimate. 6.1. Minimum measurement solution To validate the proposed minimum measurement solution, Fig. 9 performs comparison with the Fang’s solution [27] and the SX method [28]. 3-D hyperbolic localization configuration is used as the Fang and SX methods were developed for hyperbolic positioning. The object is located at uo = [15, 10, 6]T and sensor positions are s0 = [0, 0, 0]T , s1 = [10, 10, −10]T , s2 = [10, −20, 30]T and s3 = [20, −10, 20]T . In each ensemble run, s1 to s3 were perturbed slightly by adding uniformly distributed random values

14

S.S.A. Al-Samahi, Y. Zhang and K.C. Ho / Signal Processing 167 (2020) 107273 Table 5 Computational complexity comparison for the simulation results in elliptic positioning.

Relative Computation Time / 2-D Relative Computation Time / 3-D

Proposed

Proposed(NoSelection)

MLE

Solution from [11]

1 1

0.88 0.98

26.04 47.13

0.48 0.80

Table 6 Computational complexity comparison for the simulation results in hyperbolic positioning.

Relative Computation Time / 2-D Relative Computation Time / 3-D

Proposed

Proposed(NoSelection)

MLE

Chan-Ho [2]

SCWLS [48]

1 1

0.98 0.95

16.62 20.10

0.12 0.17

1.50 0.67

U[−0.5, 0.5] × 10−5 , otherwise the SX method could not produce a solution. The proposed method and the Fang method give essentially identical result, although the former only needs to solve a quadratic equation while the latter quartic. SX method is not robust and does not provide good estimate for this sensor configuration especially when the noise level is not significant.

6.2. Elliptic positioning We next consider the 2-D scenario where the object is at uo = [−15, 10]T , transmitter at s0 = [20, −30]T , and the receiving sensors at s1 = [35, 15]T , s2 = [−40, −30]T , s3 = [10, 10]T , s4 = [40, −20]T , s5 = [0, −50]T . Fig. 10 illustrates the performance when only the first three receivers are used. There are three individual solutions and only two are selected to combine together to form the final. The proposed solution is able to reach the CRLB accuracy as expected by the analysis. When the proposed scheme in Section V.D is applied to select the better individual solutions for BLUE, rather than just using the two from (d1 , d2 ) and (d3 , d1 ) (indicated by NoSelection in the legend), we are able to extend the reaching of the CRLB from 10 log(σ 2 ) = −15 to 5 and provide comparable performance with the MLE. The proposed solution also behaves better than the closed-form solution from [11] with a gain of about 3dB at 10 log(σ 2 ) = 5. Fig. 11 gives the results when we have all five receivers. The number of individual solutions is 10 and only two among those are selected by the proposed scheme to generate the final. In this particular simulation, the proposed solutions with and without individual solution selection give comparable performance. The proposed estimator has some improvement over the closed-form solution from [11] and similar performance with MLE. Fig. 12 illustrates the averaged results over 20 randomly generated geometries where the transmitter and the receivers are placed randomly in the space [−100, 100]2 from uniform distribution. The number of measurements is 3. To limit the possibility of degenerated geometry, the distance between any two of them is not less than 20. The object is randomly placed in the difference of the area between [−100, 100]2 and [−50, 50]2 . Using individual solution selection does not appear to improve results compared to using individual solutions from (d1 , d2 ) and (d3 , d1 ) to combine as shown in the Figure. However, the difference appears if we use the individual solutions from (d1 , d2 ) and (d2 , d3 ) for NoSelection instead. The results are consistent with the single geometry cases. The proposed algorithm is able to reach the CRLB performance, yields better accuracy than the closed-form solution from [11] by about 5dB at 10 log(σ 2 ) = −5, and has comparable performance with MLE although the complexity is much lower. Fig. 13 repeats the simulation in Fig. 12 for the 3-D scenario, where the transmitter, receivers, object are placed in a similar manner now over the region [−100, 100]3 . The number of measurements is 4. The observations are similar to those of Fig. 12.

The performance improvement over the solution from [11] is very significant. 6.3. Hyperbolic positioning The classical Chan-Ho method from [2] provides a closed-form solution to hyperbolic positioning. Its performance may not be sufficient when the object is near the center of the sensors arranged in a circle or sphere. The SCWLS solution [48] is an alternative, at the expense of higher complexity. We repeat the same simulation experiment in [48], where the object is at [4.9, 5.1]T and the sensors at [0, 0]T , [0, 10]T , [10, 0]T and [10, 10]T . The results are shown in Fig. 14. The Chan-Ho Solution fails. The proposed method does not have difficulty for such a configuration and provides the CRLB performance. It seems to deviate a little earlier from the bound than SCWLS as the noise power increases. To better assess performance, Fig. 15 illustrates the accuracy in 2-D positioning with the average of 20 geometries. There are 4 sensors, giving a total of 3 measurements. Their positions are generated randomly in each geometry over the area of [−50, 50]2 , with the distance between any two larger than 20. The object is placed randomly in the same area. The proposed solution works better than the Chan-Ho method and the SCWLS solution as the noise level increases. It yields comparable performance with MLE having a similar noise threshold. Fig. 16 is the results for the 3-D scenario averaged over 20 geometries. The sensors and the object are placed in a similar manner as in Fig. 15 except the region now is expanded to [−50, 50]3 . The number of sensors is 5 giving 4 measurements for positioning. In this simulation, both Chan-Ho and SCWLS do not work well and the proposed method yields much better performance even without individual solution selection. The use of the individual solution selection scheme extends the noise tolerance by about 5 dB compared to the MLE in this simulation. Algorithm 1 Minimum measurement solution selection. Input: Individual solution ui ∈ S, where S is defined in (47). Step 1. For each ui , obtain Ri = det(cov(ui )) from (16) with uo replaced by ui . Step 2. For each choice c j ∈ C where C is defined in (48), evaluate the quality factor (52). Step 3. The individual solution set for applying BLUE is the choice c∗ obtained from (53). To complete the simulation studies, we compare the computation complexity of the proposed method with other solutions in Table 4 (for Fig. 9), Table 5 (for Figs. 11 and 13), and Table 6 (for Figs. 15, 16). The complexity is obtained from Matlab implementation running on a PC with i7 processor and 8GB of memory. The values shown are computation times relative to the proposed method. It is clear from Table 6 that the proposed minimum measurement solution is more efficient than the other two that take

S.S.A. Al-Samahi, Y. Zhang and K.C. Ho / Signal Processing 167 (2020) 107273

about 1.5 and 2.5 times longer. For elliptic positioning, without solution selection in the proposed method saves a little the computation time as illustrated in Table 5. It runs much faster than the iterative MLE, although it is slower than the comparison algorithm [11]. We have similar observations for hyperbolic positioning in Table 6. The Chan-Ho method is most efficient, but it does not work well as seen in Figs. 14, 15, 16. The proposed solution is more efficient than the SCWLS method for the 2-D scenario.

We have provided a simple and direct derivation of the proposed minimum measurement solution from a geometric perspective for elliptic time delay measurement and apply it for outlier detection and object position estimation. Compared to the previous minimum measurement solutions [27,28], the proposed solution is more computationally attractive that requires the roots of a quadratic equation instead of quartic, and it is more general and robust that can work with arbitrary sensor arrangements. The solution derivation provides the condition for the intersection of two LOPs, and they are exploited to form an adjacency matrix for outlier detection using the spectral graph processing technique. Using a set of the individual minimum measurement solutions that covers all the measurements, a linear estimator based on BLUE is proposed to integrate them together to produce the final location estimate. Such an estimator by forming and combining individual solutions is algebraic and closed-form. Most important it is able to achieve the CRLB performance under Gaussian noise as validated by analysis and simulations. We also proposed an individual solution selection scheme to improve the final estimate by extending the noise level at which performance deviates from the CRLB. The results presented are applicable to elliptic time delay as well as hyperbolic time difference measurements, in both 2-D and 3-D. Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgment The work of Sanaa S. A. Al-Samahi was supported by the Ministry of Higher Education and Scientific Research, Iraq. Appendix A The pseudo-inverse of C is defined as [49]



C† = CT C

−1

CT .

(54)

From the definition of C in (37),



CT C

−1



=

0.5ILN−M 0

out of (L − 1 )N that have been used before. Thus the total number of ordered individual solutions is



L −2



CNM−iN





LN−N CLN−M .

0 I2M−LN

 .

(55)

Thus we obtain (40). Appendix B A set of L = M/N individual solutions that cover all the measurements can be obtained as follows. First, we select N measurements from M to produce the first individual solution. Second, we select another N measurements from the remaining M − N to generate the second individual solution. The process continues until the number of remaining measurements is M − (L − 1 )N. To use these remaining measurements, we choose LN − M measurements

(56)

i=0

Simplifying the first product term gives

M!

(N! )L−1 (M − (L − 1 )N )!

7. Conclusion

15





LN−N CLN−M .

(57)

Taking into consideration that the ordering of the individual solutions is irrelevant, the number of the minimum collection of individual solutions that covers all the measurements is

 LN−N  1 M! C , L −1 L! (N! ) (M − (L − 1 )N )! LN−M

(58)

which can be rewritten as (49). References [1] H. Chen, W.-P. Zhu, W. Liu, M.N.S. Swamy, Y. Li, Q. Wang, Z. Peng, RARE-based localization for mixed near-field and far-field rectilinear sources, Digital Signal Process., Elsevier 85 (2019) 54–61. [2] Y.T. Chan, K.C. Ho, A simple and efficient estimator for hyperbolic location, IEEE Trans. Signal Process. 42 (1994) 1905–1915. [3] X. Guo, L. Chu, N. Ansari, Joint localization of multiple sources from incomplete noisy euclidean distance matrix in wireless networks, Computer Commun., Elsevier 122 (2018) 20–29. [4] C.H. Park, J.H. Chang, Closed-form localization for distributed MIMO radar systems using time delay measurements, IEEE Trans. Wireless Commun. 15 (2) (2016) 1480–1490. [5] S.W. Chen, C.K. Seow, S.Y. Tan, Elliptical lagrange-based NLOS tracking localization scheme, IEEE Trans. Wireless Commun. 15 (5) (2016) 3212–3225. [6] G. Wang, A.M.C. So, Y. Li, Robust convex approximation methods for TDOA-based localization under NLOS conditions, IEEE Trans. Signal Process. 64 (13) (2016) 3281–3296. [7] X. Shi, G. Mao, B.D.O. Anderson, Z. Yang, J. Chen, Robust localization using range measurements with unknown and bounded errors, IEEE Trans. Wireless Commun. 16 (6) (2017) 4065–4078. [8] S. Asgari, J. Stafsudd, R. Hudson, K. Yao, E. Taciroglu, Moving source localization using seismic signal processing, J. Sound and Vibration, Elsevier 335 (2015) 384–396. [9] Y. Wang, K.C. Ho, Unified near-field and far-field localization for AOA and hybrid AOA-TDOA positionings, IEEE Trans. Wireless Commun. 17 (2) (2018) 1242–1254. [10] A. Guerra, F. Guidi, D. Dardari, Single-anchor localization and orientation performance limits using massive arrays: MIMO vs. beamforming, IEEE Trans. Wireless Commun. 17 (8) (2018) 5241–5255. [11] L. Rui, K.C. Ho, Efficient closed-form estimators for multistatic sonar localization, IEEE Trans. Aerosp. Electron. Syst. 51 (1) (2015) 600–614. [12] N.H. Nguyen, K. Dog˘ ançay, Optimal geometry analysis for multistatic TOA localization, IEEE Trans. Signal Process. 64 (16) (2016) 4180–4193. 15 [13] K. Dog˘ ançay, N.H. Nguyen, Low-complexity weighted pseudolinear estimator for TDOA localization with systematic error correction, in: Proc. 2016 24th European Signal Processing Conf. (EUSIPCO), Budapest, 2016, pp. 2086–2090. [14] F. Gustafsson, F. Gunnarsson, Mobile positioning using wireless networks: possibilities and fundamental limitations based on available wireless network measurements, IEEE Signal Process. Mag. 22 (4) (2005) 41–53. [15] J. Liang, D. Wang, L. Su, B. Chen, H. Chen, H.C. So, Robust MIMO radar target localization via nonconvex optimization, Signal Process., Elsevier 122 (2016) 33–38. [16] U. Hammes, A.M. Zoubir, Robust mobile terminal tracking in NLOS environments based on data association, IEEE Trans. Signal Process. 58 (11) (2010) 5872–5882. [17] Z. Yang, et al., Detecting outlier measurements based on graph rigidity for WSN localization, IEEE Trans. Veh. Technol. 62 (1) (2013) 374–383. [18] X. Chenga, G. Mishnea, S. Steinerberger, The geometry of nodal sets and outlier detection, Number Theory, Elsevier 185 (2017) 48–64. [19] M. Compagnoni, A geometrical–statistical approach to outlier removal for TDOA measurements, IEEE Trans. Signal Process. 65 (15) (2017) 3960–3975. [20] X. Shi, G. Mao, Z. Yang, J. Chen, MLE-based localization and performance analysis in probabilistic LOS/NLOS environment, Neurocomputing, Elsevier 270 (2017) 101–109. [21] P.C. Chen, A Non-line-of-sight Error Mitigation Algorithm in Location Estimation, in: Proc. IEEE Wireless Commun. Networking Conf. (WCNC), New Orleans, volume 1, 1999, pp. 316–320. [22] Y.T. Chan, W.Y. Tsui, H.C. So, P.C. Ching, Time-of-arrival based localization under NLOS conditions, IEEE Trans. Veh. Technol. 55 (1) (2006) 17–24. [23] J. Riba, A. Urruela, A Non-line-of-sight Mitigation Technique Based on ML-Detection, in: Proc. IEEE Int. Conf. Acoust. Speech, Signal Process. (ICASSP), Montreal, Canada, volume 2, 2004, pp. 153–156.

16

S.S.A. Al-Samahi, Y. Zhang and K.C. Ho / Signal Processing 167 (2020) 107273

[24] B. Song, W. Xiao, Device-free Localization with Outlier Intersection Rejection, in: Proc. IEEE Int. Conf. Commun. Probl. (ICCP), Beijing, China, 2014, pp. 690–693. [25] C. Luo, J.H. McClellan, Robust Geolocation Estimation Using Adaptive RANSAC Algorithm, in: Proc. IEEE Int. Conf. Acoustics, Speech, Signal Process. (ICASSP), Dallas, 2010, pp. 3862–3865. [26] R.O. Schmidt, A new approach to geometry of range difference location, IEEE Trans. Aerosp. Electron. Syst. AES-8 (5) (1972) 821–835. [27] B.T. Fang, Simple solutions for hyperbolic and related position fixes, IEEE Trans. Aerosp. Electron. Syst. 26 (5) (1990) 748–753. [28] M. Malanowski, K. Kulpa, Two methods for target localization in multistatic passive radar, IEEE Trans. Aerosp. Electron. Syst. 48 (1) (2012) 572–580. [29] M. Compagnoni, R. Notari, F. Antonacci, A. Sarti, A comprehensive analysis of the geometry of TDOA maps in localization problems, Inverse Problem 30 (3) (2014) 1–49. [30] M. Compagnoni, R. Notari, F. Antonacci, A. Sarti, On the statistical model of source localization based on range difference measurements, J. Franklin Institute, Elsevier 354 (2017) 7183–7214. [31] K.C. Ho, Y.T. Chan, Geolocation of a known altitude object from TDOA and FDOA measurements, IEEE Trans. Aerosp. Electron. Syst. 33 (3) (1997) 770–783. [32] S. Shuster, A.J. Sinclair, T.A. Lovell, Initial Relative-orbit Determination Using Heterogeneous TDOA, in: Proc. IEEE Aero. Conf., Big Sky, MT, 2017, pp. 1–7. [33] K.J. Cameron, D.J. Bates, Geolocation with FDOA Measurements via Polynomial Systems and RANSAC, in: IEEE Radar Conf. (RadarCon), Oklahoma City, OK, 2018, pp. 0676–0681. [34] L. Rui, K.C. Ho, Elliptic localization: performance study and optimum receiver placement, IEEE Trans. Signal Process. 62 (18) (2014) 4673–4688. [35] E. Olson, J.J. Leonard, S. Teller, Robust range-only beacon localization, IEEE J. Ocean. Eng. 31 (4) (2006) 949–958. [36] S. Al-Samahi, K.C. Ho, N. Islam, Improving TOA Localization through Outlier Detection Using Intersection of Lines of Position, in: Proc. 23rd Int. Conf. Digital Signal Process., Shanghai, China, 2018.

[37] S.M. Kay, Fundamentals of statistical signal processing, estimation theory, Upper Saddle River, NJ:Prentice Hall, 1993. [38] J.S. Abel, A divide and conquer approach to least-squares estimation, IEEE Trans. Aerosp. Electron. Syst. 26 (2) (1990) 423–427. [39] Y. Zhou, C.L. Law, Y.L. Guan, F. Chin, Indoor elliptical localization based on asynchronous UWB range measurement, IEEE Trans. Instrum. Meas. 60 (1) (2011) 248–257. [40] P.P. Klein, On the ellipsoid and plane intersection equation, Appl. Math. (Irvine) 3 (2012) 1634–1640. [41] B. Spain, Analytical conics, Pergamon Press, New York, 1957. [42] P.P. Klein, On the intersection equation of a hyperboloid and a plane, Appl. Math. (Irvine) 4 (2013) 40–49. [43] D.I. Shuman, et al., The emerging field of signal processing on graphs: extending high-dimensional data analysis to networks and other irregular domains, IEEE Signal Process. Mag. 30 (3) (2013) 83–98. [44] Y. Bar-Shalom, X.R. Li, T. Kirubarajan, Estimation with applications to tracking and navigation: Theory, algorithms and software, 1st ed, John Wiley and Sons, 1985. [45] Y. Wang, K.C. Ho, TDOA positioning irrespective of source range, IEEE Trans. Signal Process. 65 (6) (2017) 1447–1460. [46] C. Wang, F. Qi, G. Shi, J. Ren, A linear combination-based weighted least square approach for target localization with noisy range measurements, Signal Process., Elsevier 94 (2014) 202–211. [47] S. Joshi, S. Boyd, Sensor selection via convex optimization, IEEE Trans. Signal Process. 57 (2) (2009). 4513–462 [48] L. Lin, H.C. So, F.K.W. Chan, Y.T. Chan, K.C. Ho, A new constrained weighted least squares algorithm for TDOA-based localization, Signal Process., Elsevier 93 (2013) 2872–2878. [49] L. Scharf, Statistical signal processing, Addison Wesley, MA:Reading, 1991.