Practical Performance Analysis for Multiple Information Fusion Based Scalable Localization System Using Wireless Sensor Networks

In practical localization system design, researchers need to consider several aspects to make the positioning efficiently and effectively, e.g., the available auxiliary information, sensing devices, equipment deployment and the environment. Then, these practical concerns turn out to be the technical problems, e.g., the sequential position state propagation, the target-anchor geometry effect, the Non-line-of-sight (NLOS) identification and the related prior information. It is necessary to construct an efficient framework that can exploit multiple available information and guide the system design. In this paper, we propose a scalable method to analyze system performance based on the Cramér–Rao lower bound (CRLB), which can fuse all of the information adaptively. Firstly, we use an abstract function to represent all of the wireless localization system model. Then, the unknown vector of the CRLB consists of two parts: the first part is the estimated vector, and the second part is the auxiliary vector, which helps improve the estimation accuracy. Accordingly, the Fisher information matrix is divided into two parts: the state matrix and the auxiliary matrix. Unlike the theoretical analysis, our CRLB can be a practical fundamental limit to denote the system that fuses multiple information in the complicated environment, e.g., recursive Bayesian estimation based on the hidden Markov model, the map matching method and the NLOS identification and mitigation methods. Thus, the theoretical results are approaching the real case more. In addition, our method is more adaptable than other CRLBs when considering more unknown important factors. We use the proposed method to analyze the wireless sensor network-based indoor localization system. The influence of the hybrid LOS/NLOS channels, the building layout information and the relative height differences between the target and anchors are analyzed. It is demonstrated that our method exploits all of the available information for the indoor localization systems and serves as an indicator for practical system evaluation.


Introduction
Locating a target using a wireless sensor network (WSN) is an efficient way to support multiple Internet of Things (IoT) applications, and many measurement and sensing techniques are proposed [1]. The techniques of measurement, such as angle-of-arrival (AOA), time-of-arrival In this paper, we propose a general analysis method to evaluate the practical WSN localization system, which can fuse multiple pieces of information. The first major contribution is that we construct a scalable framework to model the multiple information fusion in the localization problem generally. We derive the Fisher information matrix (FIM) based on the proposed abstract function of all of the wireless localization system model instead of using a specific wireless propagation model. Furthermore, we employ an extendable estimated random vector θ, which contains all of the unknown parameters, e.g., the previous estimated state, the current state and the unknown indeterministic parameters that may influence the estimation accuracy. Then, we divide θ into two parts: the estimated state vector, which indicates the final position estimation, and the auxiliary vector, which helps improve the estimation accuracy. When θ only contains the state vectors or together with NLOS indicators, the our method is equivalent to the above-mentioned CRLBs. However, when a large number of parameters is involved in the estimation, e.g., the map information, the previous state, the relative height difference or the prior information of the unknown random parameters, our method can still provide practical optimal performance of the localization systems. The major advantage is that it is suitable for complicated and dynamic environment and fully considers the prior information, hybrid unknown factors and the recursive feature of the tracking algorithms.
The second contribution is that we employ the proposed CRLB framework to analyze the TOA range-based indoor localization system as a case study. The simulation environment considers all of the possible factors, e.g., the target-anchor geometry effect, the building layout, the relative height differences between the target and anchors, the NLOS transmission channel, the related prior information and the recursive feature of the tracking algorithm. The impact of each factor for the estimation accuracy is illustrated in the simulation. We use the spatial position error distribution (SPED) as a metric to evaluate the performance. The SPED indicates not only the target-anchor geometry effect, but also illustrates the impacts of height differences between the target and anchors, NLOS transmission and building layout information for the indoor localization. In addition, the impacts of the related prior information are also evaluated via SPED. To indicate which factor is important to the location estimation, we numerically evaluate the NLOS ranging, the height difference and the prior information in multiple scenarios. For the dynamic moving trajectories, the recursive form of the CRLB is applied. Finally, the estimation improvement using multiple anchors is also presented in the simulation. The results indicate that the NLOS ranging measurement mainly influences the estimation accuracy, and the prior information of the NLOS channel and the target position play important roles for improving the estimation. The relative height differences also degrade the estimation accuracy if we ignore them. However, with reliable prior information, we can make the estimation accuracy approach the location performance without any relative height differences. In general, our proposed method is suitable to exploit all of the available information to analyze the performance of the WSN localization system.
The rest parts of this paper is organized as follows: Section 2 provides the WSN localization system model; Section 3 introduces the scalable CRLB framework; then, we use the CRLB to derive the formulation of a practical TOA system in Section 4; the simulation evaluations and analysis are illustrated in Section 5; and Section 6 concludes the whole paper. To make the content more clear, we list all of the fundamental notations in Table 1 for the mathematical formulation.

System Model
In the WSN localization system, the mobile device with an unknown position is called the target, such as a mobile sensor node or robot. The position state of the target is denoted by where p X t and p Y t are the coordinates in the two-dimensional positioning system, and T is the transpose operator. The wireless sensor devices with known positions, which measure the ranges (or distances) to the target, are called anchors. For each anchor, the position is denoted by a j = [a X j a Y j ] T , where a X j and a Y j are the coordinates. In this paper, we assume the WSN localization system uses the time-of-arrival (TOA) method to measure the distance between the target and anchor. Thus, we will have the relative factors that may influence the TOA and the final location estimation.

Time-Of-Arrival Ranging
In the TOA measurement method, the distance between the target and anchor is calculated according to the wireless propagation time. Consider a synchronous wireless communication, where clocks at the target and anchors are strictly synchronized. Anchors send the ranging measurements periodically via a time-division multiplexing (TDM) method based on the related WSN protocol. The target receives a radio signal transmitted from one anchor via a single propagation path. Let τ j t be the time delay of the received signal from anchor j at time t: where c = 3 × 10 8 m/s is the propagation speed of the signal, a j = [a X j a Y j ] T is the anchor position and || · || denotes the distance between two positions; l j t ≥ 0 is the range drift, which is caused by the NLOS effect. The indicator l j t = 0 for the LOS propagation, whereas l j t > 0 for the NLOS propagation. For many indoor systems, the TOA ranging measurement is obtained through the packet transmission time based on the network protocol; thus, the TOA is denoted as the time observation instead of the received waveform, as depicted in [17,18]: where z j t is the range measurement for anchor j, the measurement function h j t ((d(x t , k), l j t )) = cτ j t and v j t is the measurement noise for anchor j. The measurement noise v j t follows the zero-mean Gaussian distribution v j t ∼ N (0, R j ), where R j is the variance of the ranging noise.

TOA Noise
The localization performance mainly relies on the measurement noise. Thus, it is necessary to list the factors that are related to the TOA noise. TOA relies on the quality of detecting the direct path (DP) signal. Thus, the TOA ranging not only tolerates the attenuation of the signal strength, but also depends on the peak of all received multipath signals related to the bandwidth. Thus, the radio interference, the multipath effect and the scenario variability affect the ranging error. In addition, the TOA ranging noise also depends on the system noise. The system noise comes from the unsynchronized signals and the disordered received signals, which are transmitted from different sensors simultaneously. Apart from the NLOS effect, the comprehensive impacts on the ranging noise can be generally modeled as the normal distribution [23].

Relative Height
Another important factor that should be considered in the real-world system is that the localization happens in the 3D world. In many theoretical research works, it is assumed that the anchors and targets are on the same plane in many real localization applications. The goal is to calculate 2D positions, the X − Y coordinates of the target, which is already implemented in our daily apps. In this case, the height difference between anchor and target is ignored. However, the existing height difference in the real world actually affects the accuracy in the 2D position estimation. Here, we define the height difference between the anchor and target as relative height, which is a positive variable in the Z axis. Then, we use k t ≥ 0 to denote the coordinate of the target in the Z axis and a Z j to denote the coordinate of the anchor. Thus, the relative height is k t − a z j at time t. If the relative height is zero or assumed to be zero in the simulation, we define the range measurement as 2D-ranging. The distance of 2D-ranging for each anchor is formulated as: where d j (x t , k t ) denotes the distance function from anchor j to the target. If the relative height between the anchor and target is not zero, which is always applicable in the real case, the measurement depends on 3D coordinates, then we define the range measurement as 3D-ranging. The 3D-ranging for each anchor is formulated as: Figure 1 illustrates the difference between 2D-ranging and 3D-ranging. Suppose a person or a robot is carrying a mobile device and walking in the building. Anchor 1 and Anchor 2 are deployed on the roof, and the distance to the target depends on the target position on the ground and the relative height. Even if the target is just below the anchor, the measurement is still not zero. Anchor 3 is deployed on the same plane of the target. In this case, the range measurement is 2D-ranging. However, 2D-ranging is an ideal case, and the anchors cannot be always on the same plane of the targets. No matter whether 2D-ranging or 3D-ranging, the position estimation is still for 2D in the playing field.

Scalable Framework of FIM
We attempt to construct the scalable analytical framework according to the Bayesian estimation process. The Bayesian estimation is an efficient way to fuse multiple pieces of information. Firstly, we build the predict-update model of the Bayesian estimation framework. Then, we derive the scalable Fisher information matrix.

Bayesian Model
According to the Bayesian estimation framework, the relationship between the estimated state x t and the measurement z t follows: x where (5) is the prediction function. In (5), the target's movement is based on the transition function f t (), and q t is the prediction noise, which follows the normal distribution N (0, Q t ), where Q t indicates the covariance matrix. In order to generally denote the wireless localization model, we propose an abstract measurement function of the unknown vectors: where (6) is the abstract measurement function, which is a general expression based on several range measurement techniques. In (6), z t = [z 1 t . . . z j t . . . z N t ] T is the measurement vector, and N denotes the number of anchors. Note that z t can represent the RSS vector or TOA vector, and (6) can be rewritten according to the different measurement techniques. Then, T is the nonlinear observation function, which relates to the actual received signal waveforms at the target from the anchors T is the ranging noise, which is assumed as independent noise; d() = [d 1 () . . . d j () . . . d N ()] T represents the distance vector between the target and anchors. The auxiliary parameters are defined as the optional unknown factors that may influence the estimation. In the position estimation, the auxiliary parameters are not necessarily calculated. However, the information of such parameters can improve the estimation accuracy. We define two kinds of auxiliary parameters that may influence z t in (6). The first one is the nonlinear auxiliary vector k = [k 1 . . . k j . . . k N k ] T , which affects the actual distance together with the position state x t , and N k denotes the number of elements of k. The second one is the linear auxiliary vector l = [l 1 . . . l j . . . l N l ] T , which affects the observed measurement z t and is independent of x t , where N l ≤ N indicates the number of the parameters. Since (6) is a general expression for several range measurement techniques, it should be specified and rewritten in a particular system. In this paper, (6) is simplified as the TOA formulation, which is discussed in the following sections. According to the Bayesian theorem, the posterior probability of x t is expressed as p(x t |z t , x t−1 ) = p(x t |x t−1 )p(z t |x t ), where t − 1 indicates the previous time interval, and p(x t |x t−1 ) is the prior probability [24].

Recursive FIM
Our analysis fully considers all of the possible unknown random factors that may influence the position estimation; hence, the parameter vector includes: the current state x t , the previous state x t−1 and auxiliary parameter vectors k and l. Thus, θ is expressed as: The CRLB, which is given by the inverse of the Fisher information matrix (FIM), sets the lower limit for the variance (or covariance matrix) of any unbiased estimators of an unknown parameter (or unknown parameters) [25]. If p(θ, z t ) denotes the joint probability density function (PDF) of observations z t and the state θ, then the score function is defined as the gradient of its log-likelihood: is the operator of first order partial derivatives. The FIM, J(θ), is the covariance of the score function: where E {·} indicates the expectation operator. Additionally, the CRLB is just the inverse of FIM, and the estimation covariance cannot be lower than it: where "A B" should be interpreted as matrix A − B is non-negative definite.
Since p(θ, z t ) = p(z t |θ)p(θ) is based on the Bayesian theorem, it is easily seen that J(θ) can be decomposed into two parts: where J D (θ) represents the information obtained from measurement data and J P (θ) represents the prior information. Firstly, we use the notations h = h t (d(x t , k), l), h j = h j t (d(x t , k), l) and decompose J D using the chain rule as: where H = [∇ θ h] and J h is the FIM conditioned on h: The matrix H is further decomposed into four components: where Then, H is formulated as: For J h , we can use diagonal matrices of order N to represent it: where the diagonal term λ j depends on h j t () and will be derived based on the typical system later. Then, J D , which is a (4 + N k + N l ) × (4 + N k + N l ) matrix, is written as: where: Without prior information, J D cannot indicate the CRLB directly, because det(J D ) = 0 and J D is not reversible. However, the CRLB can also be attained through a calculation rule. The calculation rule is illustrated in the next subsection. Next, we derive J P . The prior probability for θ is extended as p(θ) = p(x t |x t−1 )p(k)p(l), then the prior information is written as: where p(k) and p(l) are the independent prior information of x t and x t−1 . If we decompose θ into two sub-vectors: the state vector [x t x t−1 ] T and the auxiliary vector [k l] T . Then, J P can be formulated as: where J P 11 is the recursive form of x t and x t−1 , which is formulated by Tichavsky et al. [15]: where: where J(x t−1 ) is the previous FIM of x t−1 . Additionally, J P 12 are 0 matrices, since p(k) and p(l) are independent of x t and x t−1 : The prior distribution p(x t |x t−1 ) is also independent of l and k; thus, J P 21 = J T P 12 = 0. Finally, the last element J P 22 is expressed as: where J K and J L are the FIMs conditioned on k and l, respectively: The formulations of J K and J L are based on the features of k and l, which are discussed in the later sections. Then, substitute (17) and (20) into (11); we obtain: Using the form of the Schur complement of the submatrix [26], we decompose the FIM into two parts: the state matrix J S and the auxiliary matrix J A : where: Additionally, the formulation of each element can be found in (18), (22) and (25).

Calculation Rule
Equation (27) only holds when all of the elements in θ are to be estimated and the prior information for the whole θ is available. If the CRLB is used to analyze a specific localization system, not all of the elements are necessary for θ, and some vectors are absent sometimes. For instance, for the non-recursive scenario, the system does not consider x t−1 . In addition, when the system has a deterministic value of the auxiliary vectors, k and l are not estimated and useless for J(x t ). Thus, the calculation principle for our CRLB is that: when any vector in θ is absent, the related matrix in (27) turns into 0, and we will treat such a 0 matrix as the empty matrix, then we mitigate the empty matrix for further calculation.
Here, we consider two examples. The first example is that the estimated vector θ does not include x t−1 and k, then M 12 , M 22 , J(x t−1 ), J K , D 13 , D 33 and D 34 are 0 and mitigated from (27): (29) which has the same form as the formulation of E-CRLB in the wideband NLOS environment [18] and L-CRLB for the linear formulation [19].
The second example is that the estimated vector θ does not include k and l, and J A is empty, then only J S remains, which is the form for the recursive nonlinear filter [15]:

Application to a TOA WSN Localization System
In this section, we apply the proposed CRLB to analyze the time-of-arrival (TOA)-based WSN localization systems. Note that this paper provides a scalable framework for theoretical analysis rather than a particular analysis. Thus, our goal is not to find the major impacts or factors via mathematical derivation, but through the numerical simulations. In addition, according to the formulation of (27), many factors are involved in the calculation. It is not easy to analyze them via a single formulation. Thus, after enumerating the important factors in the TOA system, we will evaluate the performance via numerical analytical simulations. The impacts of the practical conditions, e.g., the real 3D anchor deployment and the NLOS transmission based on the building layout, are formulated. The related parts in (26) are derived for both cases with and without prior knowledge.

Relative Height
As mentioned before, no matter whether 2D-ranging or 3D-ranging, the position estimation is still for 2D in the playing field, which means we just want to obtain x t = [p X t p Y t ] T and not k t . However, for CRLB analysis, we still consider if k t is involved in the calculation and check whether the estimation accuracy is influenced. Then, H t is obtained: where: 1×N is expressed as: Since x t and k t are within the same function d(x t , k t ), H t and K are both derived from d(x t , k). In order to analyze the impacts of x t and k t , we decompose the matrix and obtain (31) and (33). For a 3D positioning case, we can combine H t and K together to indicate the 3D position state. However, it is more convenient to separate H t and K to analyze the estimation performance of x t in a 2D positioning case.
If any prior information is unknown to the system, substitute (31) and (33) into (27), then we obtain: For the prior information of k t , we assume that the target is always above the ground, which is k t ≥ 0. Thus, we apply the Gamma distribution to indicate the potential distribution of k t , , α k is the shape parameter and β k is the rate parameter. For the Gamma distribution, J k is complicated. To obtain an analytical expression, we assume α k > 2 for simplicity. Then, the Gamma function is is derived as: Use the property Γ(α k ) = α k Γ(α k − 1), and substitute it into (36); we obtain:

NLOS Impact
We use the vector l to represent the NLOS indicator. In the TOA system, the NLOS delays the wireless packet propagation linearly. Thus, l is the linear auxiliary vector. We assume there are N l ≤ N NLOS measurements, and the drift for each measurement is independent of the others, then L = [∇ l h] N l ×N is formulated as: where I N l is the identity matrix of order N l , and the rest is a N l × (N − N l ) zero matrix due to the independent condition for the LOS measurement. As mentioned before, the zero matrix is mitigated during the calculation, then L = I N l . Note that L indicates the sub-matrix of the FIM, which means that the NLOS condition is unknown to the system. Then, the system should also estimate the NLOS parameters together with the position state in the real environment. In this case, the estimated error is propagated and degrades the position estimation. Since the range drift for the NLOS is also nonnegative, we still use the Gamma distribution as the prior information l m ∼ G(a m , b m )(l m ) = (b m ) am Γ(a m ) l a m −1 m exp(−b m l m ) according to [17], where a m ≥ 2 is the shape parameter, b m is the rate parameter and m is the m-th NLOS measurement. Similar to (37), we obtain J L : where J L is the prior matrix for L, which illustrates that the system can use the prior information to determine whether a measurement is LOS or NLOS if this information is unknown. In this case, the estimation error can be reduced to some extent.

Simulation
We set up several WSN localization simulations to evaluate the analytical performance. In each simulation, we consider several different factors, e.g., the recursive process during the target tracking, estimations with and without considering l and k t . To make the results clear, we mark the CRLBs for different situations by adding superscripts and subscripts, which can be depicted as CRLB ...
... . The subscripts indicate the considered vectors, including the state vector and the auxiliary vector. This means that the system contains other factors in the ranging measurement, e.g., relative height or NLOS measurement or both. The superscripts indicate the available prior information of the related vectors.
For instance, if we want to simulate the estimation with the NLOS range drift, the results of the CRLB are marked by CRLB x t ,l . Additionally, if the prior information of x t is attained in the simulation, the results are marked by CRLB x t x t ,l . For the recursive estimation, we use the notation CRLB x t−1 ,l x t ,l to indicate the results, which considers the prior state x t−1 . Note that, the superscript of x t in CRLB x t ,... to analyze the dynamic continuous target tracking scenarios. Some CRLB notations that appear in the following sections are listed below.
• CRLB x t : Classical CRLB without any prior information, relative height or NLOS measurement • CRLB x t ,k t : CRLB when the relative height exists • CRLB x t ,l : CRLB with the NLOS measurement • CRLB x t ,k t ,l : CRLB both contains the relative height and NLOS measurement To approach the real environment, we initially set the related parameters according to the test building in [27]. Since the localization system can be affected by many factors, we tune the parameters to provide a comprehensive analysis in the following simulations, e.g., k t , l and the number of anchors. To better understand the performances in multiple environments, we uniformly locate the anchors in the SPED evaluation. Additionally, we evaluate the performances considering the relative height, NLOS and recursive estimation based on randomly-deployed anchors and Monte Carlo simulation to draw a general performance. In the real experiment, the anchors are randomly deployed throughout the whole building.

Spatial Position Error Distribution
In the first simulation, a 100 × 100 m 2 playing field with four anchors is constructed. All of the anchors are deployed on the roof being 2.5 m high, and the target is about 0.5∼1 m. To approach the real applications, we use several statistical results according to [27]. We set the relative height as the constant value 1.5 m. The range error for each anchor follows zero-mean Gaussian distribution v j t ∼ N (0, R j t ), where R j t is the variance of v j t and is set to 5 2 . The range drift for the NLOS measurements is set to 2 m. For the prior information, k t ∼ G(2.5, 2)(k). The prior information of the NLOS range drift l m is l m ∼ G(3.5, 1.8)(l m ). For the position state prior information, we assume that the prediction function is a linear static identity matrix with the zero-mean Gaussian prediction noise q t ∼ N (0, Q t ), where Q t = diag(σ 2 x , σ 2 y ) is the covariance of q t . We assume σ x = σ y = 2 m. We apply the CRLB to indicate the optimal squared error, which is tr(J −1 (x t )). To illustrate the target-anchor geometry effect for the 2D localization system in the playing field, we employ tr(J −1 (x t )) to depict the spatial position error distribution (SPED) [28]. The SPED is defined as the distribution of the position error for every possible target position, which estimates the positions point by point in the whole playing field and draws statistical results. Thus, every position in the playing field is assumed as the target position during the simulation. The SPED is derived according to all of the statistical results of the estimation error of all of the positions in the playing field. It illustrates how the performance changes according to the target-anchor (RX-TX) geometry effect. When the anchor positions and the error model change, the SPED changes accordingly, which helps to understand the relationship between the anchor deployment and the algorithm. The SPED results are represented by the contours in Figure 2, where the four anchor positions are marked by triangles.
In Figure 2a, the SPED is drawn based on the classical CRLB (CRLB x t ) in which the unknown vector is only x t , and no relative height or NLOS propagation is introduced. When both k t and l are introduced in the simulation, the dimension of the estimated vector is increased again. Then, a special pattern of the contours is drawn in Figure 2b, which demonstrates that the estimation error is propagated and increased if more unknown factors are introduced. Each anchor seems to be a center of an independent contour area, and the estimation error is extremely larger than Figure 2a, which indicates the strong uncertainty and a special geometric relationship between the target position and the anchor positions. Since more factors cannot be ignored in the complicated environment, the localization problem turns out to be a high dimensional estimation problem. It is quite possible that the squared error of the high dimensional estimation is larger than the low dimensional estimation, as the probability of the wrong estimation increases if more unknown parameters appear. Additionally, it also imposes new error on the original x t estimation error. However, such a high dimensional localization problem cannot happen in the real world. On the one hand, the calculation complexity is increased dramatically in the real system, and only x t is useful, which is a waste of computational resources. On the other hand, the algorithm designers ignore some unimportant factors, simplify the calculation and increase the estimation accuracy based on the prior knowledge.    Then, the prior information of x t , k t and l is used in the CRLB to improve the estimation accuracy. The SPEDs with prior information are depicted in Figure 2c,d. When the prior state information of x t is introduced in Figure 2c, the estimation error is reduced effectively from more than 5 m in Figure 2a to 1.4 m below. In addition, the position error is almost the same everywhere, although the geometric pattern is similar to Figure 2a, which indicates that the geometric impact is reduced by using the prior information of x t . In Figure 2d, we assume that all of the prior information is available, then the estimation error is reduced to below 1.4 m, just as Figure 2c. Although Figure 2d still has a geometric pattern, which is similar to Figure 2c, the errors in different positions are almost the same. Thus, the geometric effect is actually reduced. In addition, since the estimation error is not further reduced, the results in Figure 2d can be the limits based on the prior information. It is also demonstrated that the prior information of x t significantly reduces the estimation error and is the most important prior information especially in the recursive estimation. For the other two parts of the prior information, they are used for the improvement of the impacts of the relative height and the NLOS effects.

Impact of Relative Height
In this simulation, we evaluate the impact of the relative height independently. For the indoor localization, the relative height cannot be too high unless the room is tall enough. Thus, we choose the value of k t between 0.5 m and 10 m. The simulation evaluates CRLB x t ,k t and also CRLB x t as a comparison. To illustrate the general relationship between the relative height and the measurement noise, we also tune the variance of the noise from 0.5 2 to 5 2 and k t from 0.5 m and 10 m. In this simulation, we run 1000 Monte Carlo experiments. In each experiment, we choose a random position for analysis. The averaged results are listed in Figure 3.
It is clearly observed that the estimation error rises according to the increase of k t . However, the increased value is quite small. Take Figure 3c for instance; the increased error is only 0.06 m when k t is tuned from 0.5 m to 10 m, and the estimation error is 3.43 m by that time. Therefore, the value of k t does not affect the estimation accuracy too much. On the other hand, the existence of k t does affect the performance no matter what the value of CRLB x t is. Compare the results in Figure 3c; k t affects the estimation with a typical certain value. When k t is small, the average RMSE of CRLB x t is 1.4390 m, but the average RMSE of CRLB x t ,k t =0.5 is 1.6856 m and CRLB x t ,k t =10 is 1.7137 m. According to the analytical results, when the relative height is introduced into the simulation, there is a gap between 2D localization without the relative height. Additionally, the gap becomes larger with the increased value of the measurement noise. Take v j t ∼ N (0, 0.5 2 ) for instance; when k t = 0, the average RMSE is only 0.84 m. When v j t ∼ N (0, 5 2 ), the average RMSE without the relative height is about 8.3 m, which is almost 10-times that of the low noise. Thus, the relative height really affects the estimation no matter the value.
In addition, we evaluate the impact of the prior knowledge of k t in the multiple noise environment. The simulation results are presented in Figure 4. It is clearly observed that the prior knowledge can significantly improve the RMSE. If v j t ∼ N (0, 5 2 ), the prior information of k t can even reduce 1 m RMSE. Additionally, if the error is small, the improvement is only a little bit. However, with the increased value of the relative height, the RMSE rises accordingly even with the prior information. For the real application, the relative height can be within 10 m, in which the improvement can effectively reduce the RMSE to a reasonable range. Therefore, even if the impact of the relative height is limited, it is still necessary to employ the prior information.

Impact of NLOS
The analysis of NLOS measurements has been mentioned in several literature works [17,18]. Here, we evaluate the impact of the number of NLOS measurements through simulations. The playing field is still 100 × 100 m 2 with 16 randomly-deployed anchors. We set up 1000 Monte Carlo runs, and a random set of anchors is chosen to have NLOS measurements. We tune the number of NLOS measurements from 1 to 16 and evaluate the performance of CRLB x t and CRLB l x t ,l . The averaged results are depicted in Figure 5. It is observed that the RMSE of CRLB x t ,l increases according to the rise of the NLOS measurement number. It is also possible to have useless estimation when all of the measurements are NLOS, just as mentioned in [17]. When there is only one NLOS measurement, the estimation error reaches 3 m. Additionally, the error is increased to 5 m if the number of NLOS measurements is changed to 16. In real applications, the number of NLOS measurements is unknown, and the system has to calculate all of the parameters from all of the measurements. Thus, the estimation error is quite high without the help of prior information. In addition, with the prior information of l, the estimation error is almost a constant value. Thus, the prior information of the NLOS measurement can successfully reduce the NLOS effect, which is applicable for real applications. Such prior information can be a combination of the building layout information and the wireless propagation model, e.g., map-matching-based algorithms [13]. Additionally, several literature works have proposed many algorithms that either use pre-assumptions or estimate the parameters to obtain the NLOS information.
The NLOS-related methods are beyond the scope of this paper; please refer to other research works in [20,23,29].

Map Assist Localization System
In the following simulation, we evaluate a practical scenario where the layout of the building is involved. The playing field is still 100 × 100 m 2 . There are four big rooms located at four corners of the playing field. The area for each room is 40 × 40 m 2 . The rest of the playing field is the hallways. Four anchors are deployed inside each room and placed at the four corners. Thus, 16 anchors are uniformly deployed in the playing field, which is depicted in Figure 6a. The triangles mark the positions of anchors. The LOS measurements can only be obtained in one room with four associate anchors, and others measurements are NLOS. For the positions in the hallways, all of the measurements are NLOS rangings.
All of the anchors are deployed on the roof with relative height 1∼2 m to the target. The distributions of the range error and the related prior information are the same in Figure 2. The range drift for the NLOS measurements depends on the signal transmission from the anchor to the target. According to our previous research, the average positive NLOS drift of the signal through one wall is set to 2 m and to 5m for the signal through two or three walls if the signal can be detected [30]. The prior information of the NLOS range drift l m depends on the signal transmission. Here, we use the conclusions based on the experimental models in [27]: the distribution for the signal transmission through one wall is l m ∼ G(3.5, 1.8)(l m ); the distribution for the signal transmission through more than two walls is l m ∼ G(4.5, 2.2)(l m ). In this simulation, we evaluate the SPED of CRLB x t ,k t ,l , CRLB k t ,l x t ,k t ,l and CRLB x t ,k t ,l x t ,k t ,l , which means that every position is calculated, and the error contours are drawn to indicate the target-anchor geometry effect. The results are depicted in Figure 6.  (d) Figure 6. The simulation of a building layout. (a) Anchor deployment of the playing field; For numerical comparison, the RMSE in Figure 6b is higher than in Figure 6c,d, which is more than 3.39 m in the central area. Additionally, the error becomes higher and higher when the position is approaching the corner, which is more than 8 m. Due to the lack of prior information, the geometric shape does not have special characteristics, which are related to rooms or corridors. The contours are almost like rectangles located in the center of the playing field. When the prior information of k t and l is introduced, the accuracy is significantly improved, which is reduced to 2.55 m on average. The geometric shapes of the contours are different in the rooms and hallways. This indicates that the localization algorithms using the prior knowledge of NLOS conditions based on the building layout information and the NLOS identification and mitigation methods can reasonably improve the estimation performance. Thus, the layout information in the building map is an important information source for localizations. When the prior information of x t is introduced in the estimation as indicated in Figure 6d, the RMSE is further reduced, which is 1.235 m in almost all of the playing field where the target-anchor geometry effect is reduced effectively.

Bayesian-Based Target Tracking Estimation
In this simulation, we evaluate the performance of the recursive Bayesian estimation for target tracking. Since the impacts of k t and l are extensively analyzed above, this section mainly focuses on analyzing the impact of x t−1 . The playing field is still 100 × 100 m 2 , and 16 anchors are deployed randomly in the field. We evaluate the CRLB in three scenarios: The first one considers NLOS and relative height measurement; the second scenario considers only the relative height; the third scenario only has NLOS measurement. In each scenario, both CRLB without prior information and the recursive CRLB with prior information based on (27) are analyzed. Therefore, we evaluate CRLB x t ,k t ,l and CRLB x t−1 ,k t ,l x t ,k t ,l in the first scenario, CRLB x t ,k t and CRLB x t−1 ,k t x t ,k t in the second scenario and CRLB x t ,l and CRLB x t−1 ,l x t ,l in the third scenario. We assume that the prior information is unknown initially in the recursive estimation. This information can only be obtained recursively after the initial estimation is obtained. Then, the recursive estimation is applied.
We run 1000 Monte Carlo simulations, and the target moves in a separate random path in each simulation. In addition, the target can also be static. Since x t−1 can also be estimated in the static scenario and be used for recursive estimation, the analysis results are the same as the dynamic target tracking scenarios. The estimation results are averaged and represented by the RMSE in Figure 7. There are three solid parallel straight lines, which indicate the estimations without prior information: CRLB x t ,k t ,l , CRLB x t ,k t and CRLB x t ,l . The three other dashed curves illustrate the recursive estimations according to time steps, which are CRLB The Bayesian recursive estimation method with related prior information effectively reduced the estimation as indicated in Figure 7. The RMSEs of the three curves gradually converge to low values according to time steps. The impacts of the relative heights and the NLOS drifts still degrade the estimation performance. Even with the recursive estimation, the estimation error cannot be further reduced, where the CRLB x t−1 ,k t ,l x t ,k t ,l is 0.5 m larger than CRLB x t−1 ,k t x t ,k t when t = 20.

Multiple Anchors
The relationship between the CRLB and the number of anchors is illustrated in Figure 8. We randomly deploy multiple anchors in the playing field. The number of anchors is adapted from four to 30. Both k t and l are considered in this simulation. We assume that all of the measurements are NLOS. In this simulation, we evaluate CRLB x t ,k t ,l , CRLB k t ,l x t ,k t ,l , CRLB x t ,k t ,l and CRLB x t−1 ,k t ,l x t ,k t ,l . With a few anchors, CRLB x t ,k t ,l contains large error, which achieves more than 15 m. When the number of anchors is four, the prior information of k t and l can reduce much of the error as indicated by the curve of CRLB k t ,l x t ,k t ,l . Furthermore, the recursive estimation based on CRLB x t ,k t ,l can reduce much of the error, which is even smaller than CRLB k t ,l x t ,k t ,l . The error of the recursive estimation based on CRLB x t ,k t ,l is only a little lower than CRLB x t ,k t ,l . It is demonstrated that compared to the prior information of k t and l, x t−1 is the dominant factor for reducing the estimation error. However, when the number of anchors is increased, all of the curves converge gradually to a low value. In this case, the prior information does not improve much for the estimation. Thus, even if the prior information is not available, using more observations still improves the estimation accuracy effectively.

Practical Evaluation
For practical usage, we employ the CRLB to evaluate a reference system. In this system, we deployed 17 wireless sensor nodes either along the corridor or in the offices of the research building. A robot carrying a sensor node as a target moved along the corridor of the building with constant speed while recording its own positions [31]. All sensors are integrated with the nanoPAN 5375 RF (Nanotron Tech. GmbH, Berlin, Germany) module with a 2.4-GHz transceiver and a 1-Mbps data rate for range measurement, the LPC 2387 as the micro-controller (Nanotron Tech. GmbH, Berlin, Germany) and the CC1101 900-MHz transceiver (Nanotron Tech. GmbH, Berlin, Germany) as the radio transceiver for communication. The data collected from sensor nodes are also range measurement values, which are based on TOA. Figure 9 depicts the map of our experimental building. The triangles, which are randomly deployed, mark the sensor node positions.

Discussion
In this section, we use our proposed method to evaluate the optimal performance in multiple environments. It can be observed that our method can adapt the estimated parameters to fit the practical environment. In addition, the scalable architecture can effectively fuse multiple pieces of information, which outperforms other methods, which only consider parts of the information. In Table 2, we compare our method to three mainly used CRLBs, which are generalized CRLB [17], conditional CRLB [16] and equivalent CRLB [18]. The first row in Table 2 indicates the considered information, and the word "prior" means the related prior information. As illustrated in the table, the generalized CRLB can only analyze the NLOS effect. The conditional CRLB provides the formulation of the state vector instead of other factors. The equivalent CRLB exploits more information. However, it is a closed form formulation, which is not scalable and cannot fuse the relative height information.

Conclusions
In this paper, we propose a scalable analyzing method for a WSN localization system, which can fuse multiple heterogeneous information to indicate the optimal performance. Theoretically, we divide the estimated vector θ into three parts: the estimated state vector and two auxiliary vectors. The recursive formulation of FIM is provided, which fully considers all of the possible factors that may influence the estimation accuracy, and it exploits all of the available information to derive the fundamental limits. It is a suitable tool to indicate the optimal estimation bound of practical systems.
We employ our theoretical contribution to analyze the TOA range-based WSN localization system. The impacts of the height k t and NLOS range drift l are considered as the auxiliary vectors. The target-anchor geometry, prior information and the recursive form are extensively analyzed in the simulation. In addition, we employ a real test-bed and practical data to evaluate the overall performance. In the simulation, we find that both the relative height and NLOS impact can heavily degrade the estimation performance. However, many pieces of available information can improve it, e.g., the prior state distribution, the prior knowledge of the building and multiple anchors.
According to the simulation and experimental demonstration, the proposed CRLB is a general framework for analyzing the WSN localization systems, and it is not restricted to any specific technique. Future work will use the scalable CRLB to exploit other localization techniques and in other complicated environments to find potential factors that may influence the performance.