Multi-Sensor Consensus Estimation of State, Sensor Biases and Unknown Input

This paper addresses the problem of the joint estimation of system state and generalized sensor bias (GSB) under a common unknown input (UI) in the case of bias evolution in a heterogeneous sensor network. First, the equivalent UI-free GSB dynamic model is derived and the local optimal estimates of system state and sensor bias are obtained in each sensor node; Second, based on the state and bias estimates obtained by each node from its neighbors, the UI is estimated via the least-squares method, and then the state estimates are fused via consensus processing; Finally, the multi-sensor bias estimates are further refined based on the consensus estimate of the UI. A numerical example of distributed multi-sensor target tracking is presented to illustrate the proposed filter.


Introduction
The well-known Kalman filter is optimal given a linear Gaussian model, but its performance will deteriorate in the presence of unknown biases, such as the unknown and time-varying delays in chemical processes [1][2][3][4], faults or failures in fault-tolerant diagnosis and control systems [5,6], registration errors in multi-sensor fusion [7][8][9][10], or inertial drift in navigation [11][12][13][14]. In general, such a bias is represented as an unknown input (UI) to the nominal model. To the best of our knowledge, approaches to UI modeling and the design of corresponding filters typically fall into one of the following categories.
The first type of UI is zero-mean random noise with unknown covariance. In the case of the stationary noise processes in linear dynamic systems, the covariance can be identified via Bayesian or maximum likelihood estimation. A corresponding filter has been applied for orbit determination for near-earth satellites [15,16]. Recently, this type of filter design has been extended to time-varying covariance [17] and jump Markov stochastic systems [18]. Furthermore, an M-robust estimator [19] has been derived for the simultaneous adaptive estimation of unknown states and observation noise statistics. An auto-covariance least-squares method [20] has been presented to achieve lower-variance covariance estimates along with the necessary and sufficient conditions for the uniqueness of the estimated covariances. By treating unknown covariances as missing data, the adaptive estimation problem has been transformed into a problem of the joint optimization of state estimation and parameter identification, allowing it to be solved in an expectation-maximization (EM) iterative processing framework [21].
The second type of UI is unknown but deterministic. Using least-squares estimation and moving-window hypothesis testing, the corresponding filters can cope with cases in which the UI is piecewise constant [22] or a sum of basis functions with piecewise-constant weights [23]. The problem diag{·} indicates a diagonal matrix; col{·} denotes column augmentation; E{·} denotes the operator of mathematical expectation; the superscript i indicates a matrix or variable related to the ith sensor; rank{·} denotes the rank of the specified matrix; the subscripts "i"and "j" indicate the ith and jth subblocks of a matrix, respectively; the superscripts "∧" and "˜" indicate the estimate and residual, respectively, of a random vector; the subscripts "k|k + 1"represents the estimate of a vector at the time k using measurements up to time k + 1 and denotes the Kronecker product.

Problem Formulation
Sensor bias (SB) is widespread in multi-sensor systems and originates from two types of sources: SBs of the first type are local SBs (LSBs), which are independent among sensors, such as calibration error, navigation bias and sensor faults/drift, whereas sensor biases of the second type are global SB (GLSB, also called UI in this work),which are common to all sensors, such as electronic countermeasures (ECMs) in target tracking. A GSB model has been proposed based on dynamic evolution driven by both independent LSBs and a common UI [41]. However, the corresponding filter has two shortcomings: first, all sensors must have an identical measurement matrix, and second, the system state and the UI cannot be estimated simultaneously because the SB estimation is obtained by deriving the equivalent GSB model neglecting system state and UI. Here, we present the two assumptions adopted in our work. Assumption 1. All measurements made by sensors in the network are continuously interrupted by the GSB, which includes LSBs and GLSB.

Assumption 2.
The network topology is fixed, and the target can be monitored by all sensors throughout the entire process. There is no fusion centre or lead sensor in this sensor network. In other words, all sensors are equal in status, and they reach a network consensus by exchanging information with their neighbours.
Consider a scenario of collaborative target tracking within a sensor network, as shown in Figure 1. A target is continuously moving within a surveillance zone. Each sensor node obtains raw measurements of the target corrupted by the GSB. At each step, each sensor obtains local estimates of the target state and GSB, exchanges its local system-state estimate with its neighbors, and obtains fused estimates of the target state and UI via network consensus. In other words, this sensor fusion process can be regarded as the joint distributed estimation of the target state, GSB and UI without the limitation that all sensor measurement matrices must be identical. Consider a network of s sensor nodes. The network topology is represented by an undirected graph G = (V, E , A), where V = {1, 2, ..., s} denotes the set of nodes; E ⊂ V × V is the set of permissible communication links; and A = a ij is the graph topology matrix, the elements of which are defined as the set of neighbors of the ith sensor, and let J i = N i ∪ {i} denote the inclusive neighbor set of sensor i. Here, the presented network is assumed to be a connected graph. Otherwise, at least two separate sub-networks would exist, and hence, information could only be fused within each sub-network instead of throughout the entire network because information could not be exchanged between the two.
Here, we formulate the following problem of distributed sensor fusion in the presence of time-varying sensors subject to a common UI: where x k ∈ n is the state vector; y i k ∈ m i is the measurement vector of the ith sensor; .., s}; d k ∈ q is the common UI; and ζ k , w i k and v i k are all independent, zero-mean, Gaussian, white-noise components with covariances of Q k , S i k and R i k , respectively. The initial target state x 0 and ith GSB state b i 0 are Gaussian distributed with respective known means ofx 0 andb i 0 and associated covariances of and G i k are known and have appropriate dimensions. Under the assumption that m i > p i ≥ q i , the G i k are of full column rank.
Our aim is to design an unbiased minimum-variance filter to jointly estimate the target state x k , the GSBs b i k and the UI d k given the distributed measurements recorded up to time k.

Remark 1.
The GSB b i k in Equation (3) represents a generalized bias which includes multiple bias types as the special case: (1). If F i k = I, G i k = 0, S i k = 0, then b i k+1 = b i k represents a constant-value bias. (2). If F i k = 0, G i k = 0, then b i k+1 = F i k b i k + w i k represents a time-varying bias with known dynamic model. (3). If both F i k = 0, G i k = 0, then b i k+1 = w i k represents a zero-mean random error. (4). If G i k = 0, then b i k represents the dynamically evolving bias driven by arbitrary UI. The second type of the UI mentioned in the introduction is equivalent to the case (1). The third type of the UI mentioned in the introduction is equivalent to the case (4). The first and fourth types of the UI mentioned in the introduction can be contributed to the case (4) when these two types of UI can be decoupled from Equation (3).
The above four types of UI can be accommodated within the single-model system proposed in this work. However, for the fifth type of UI mentioned in the introduction, it is usually regarded as a randomly switching parameter obeying a known Markov chain which is classified to the multiple-models system, so this type of UI cannot be accommodated within the presented model in this work. Nevertheless, for this type of UI, if GSB model can be established for the each mode of Markov chain. According to Equation (3), GSB model can be written as where r represents the number of states of the Markov chain. Meanwhile, according to Equation (10), the decoupling condition of the corresponding GSB model can be also confirmed as These decoupling conditions for a UI obeying a Markov chain are more stringent than those considered in this work, as Equation (5) must hold for every r. Then the fifth type of UI can be accommodated within the proposed framework.

Remark 2.
There are two distinct strategies of networked sensor fusion. The first is the centralized strategy [41,42], i.e., all sensors transmit their measurements to a fusion centre, which is responsible for integrated data processing for state estimation and/or parameter identification. The second is the distributed strategy [25], i.e., each sensor processes its own measurements and then shares its processing results with its neighbors in an iterative manner to achieve consistent fusion. Obviously, the distributed fusion strategy allows the computational burden to be shared among sensors and consequently is more desirable for large-scale networks. However, distributed fusion is much more complex, involving both local processing and global fusion. Currently, distributed fusion in the presence of time-varying sensor biases driven by a common UI remains an open problem.

Remark 3.
Concerning our proposed problem, there exist two possible solution approaches. One is iterative optimization between state estimation and UI identification, as in the well-known expectation-maximization (EM) scheme [25]. As shown later via simulation, the EM scheme is not desirable here because it must treat the UI as a constant or slowly varying parameter in each iterative window. The other approach is UI decoupling plus state estimation, which serves as the basis for our proposed method, as derived later. The main technical difficulty of this approach is to determine how to adaptively decouple the UI and then fuse the local state and UI estimates.

Main Results
Because of the presence of the UI, the target state and GSB cannot be estimated directly. In this section, an equivalent UI-free dynamic model of the bias is derived and the necessary and sufficient conditions for the UI-free model are explored. Afterwards, optimal estimators, in the sense of minimum variance, are proposed for the target state x k and the local GSBs b i k .

Local LMMSE Filter for System State and Local SB
Lemma 1. For matrices Γ, Ψ and I with proper dimensions, if the matrix Ψ is of full column rank, the equality will hold if and only if rank{ΓΨ} = rank{Ψ} Proof. See Appendix A.
For convenience of derivation, the model defined in Equations (1) and (2) can be rewritten as follows: where (8) and (9). If and only if

Theorem 1. (UI decoupling). Consider the system model defined in Equations
then the UI will be decoupled from z i k+1 , i.e., there will exist the following UI-free dynamic model of the sensor bias:

Remark 4.
If the matrices Γ and Ψ satisfy Equation (7), then the row rank of the matrix Γ should be no less than the column rank of the matrix Ψ in Lemma 1. As shown in Equation (10), this means that the maximum number of UIs that can be decoupled from the system must be fewer than the number of independent measurement components of all sensors.
As shown in Equation (11), the measurements y i k are introduced into the equivalent bias dynamic model to decouple Here, M i k+1 = 0 means that the process noisew i k is correlated with the measurement noise v i k . Thus, designing an unbiased minimum-variance linear filter for the system defined in Equations (1) and (2) is equivalent to finding the MMSE estimation for the system defined in Equations (9) and (11) with correlated process noise and measurement noise. (9) and (11) have the following recursion relations:

Theorem 2. The MMSE estimators for the UI-free system defined in Equations
where the estimator gain is Proof. The termC i k y i k+1 is an additive known input. Equations (9) and (11) are linear, with additional Gaussian noise. Hence, the optimal filter with correlated process noise and measurement noise [43] can be directly utilized.
The framework and pseudo-code for the proposed local MMSE estimation are presented in Figure 2 and Table 1, respectively.

Equivalent system state evolution model
Equivalent UI-free state evolution model

Joint estimation of state and biases
) Target state evolution model Sensor bias evolution model Figure 2. The scheme of joint MMSE estimation of the state and the SB. Table 1. Pseudo-code for the local MMSE estimation.

Check the Applicability of the Method
Step 1: Check the existence condition given by Equation (10). If the condition is satisfied, then perform the following computation. Otherwise, the proposed method is not applicable.

Offline Equivalence Transformation
Step 2: Construct the equivalent system-state evolution model by definingF i k , A i k ,C i k ,B i k , andw i k using Equations (12)-(16).
Step 3. Determine the parameters Π i k and M i k using Equations (12) and (18).
Step 8: Set k = k + 1 and go to Step 1.

Consensus Estimation of the Target State and the UI
Given the local estimateẑ i k+1|k+1 of the ith sensor derived in the previous section, each sensor exchanges its own estimates with its neighbours to reach consistent consensus estimates of x k and d k . The corresponding consensus estimates consist of two parts: an average consensus filter for x k and a distributed least-squares (DLS) estimate for d k .

Local LS Estimator of UI
Let us rewrite the sensor bias evolution model defined in Equation (3) as in Section 2. It is easily confirmed that E ξ i k = 0. We can rewrite Equation (24) as the following compact expression: : j ∈ J i Furthermore, based on the local GSB estimates of the ith node and its neighbours, the local DLS estimate of d k and the corresponding covariance of the ith sensor are obtained as follows:

Average Consensus Filter for the Target State and the UI
Now, estimates of the target state and the UI have been obtained, but they are local. In other words, the estimates from different nodes are not identical to each other. Therefore, an average consensus fusion (ACF) is presented here to obtain a uniform result at each node via iterative interaction.
Let W (l) = w ij (l) denote the linear weight matrix in the lth step of iteration. Here, w ij (l) is the measurement weight of node j at node i, which satisfies represent the weight matrices for the target state and the UI, respectively.
Given the consensus covariances of the state estimate and bias estimate of node i after the lth iteration, denoted by P xx,i k|k (l) and P dd,i k (l), respectively, we design the weight matrices as follows: such that local estimates of higher accuracy will be assigned higher weights in the weighted fusion process.  (1)- (3), with topology G = (V, E , A) and weight matrices w x ij (l) and w d ij (l) for the target state and UI, respectively, at the lth step of iteration. The ACFs for x k and d k arex where P xx,i Proof. See Appendix C.
In the ith node, Equations (27)-(30) will remain iterative until tr P xx,i k|k (l + 1) − tr P xx,i k|k (l) tr P xx,i k|k (l) and tr P dd,i k (l + 1) − tr P dd,i k (l) where η x and η d are positive thresholds for system state and UI, respectively.

Refinement of the Local Bias Estimates
After the consensus iteration ofx i k|k andd i k−1 , the sensor biases can be refined according to whered i k−1 is the consensus estimate at the kth sampling instant. The covariance ofb i k|k is as follows: Thus, the estimatesb i k|k , P bb,i k|k ,x i k|k , P xx,i k|k ,d i k−1 and P dd,i k−1 will be output as the final results at the kth sampling instant.
The framework and pseudo-code for the iterative consensus process are presented in Figure 3 and Table 2, respectively.

Formulate the Local Information Before Transmission
Step 1: Calculate the estimatorsb i k using Equation (24).

Exchange Information with All Neighbours of the ith Sensor
Step 2: Exchange G i k ,x i k|k , P xx,i k|k , andb i k with all neighbouring sensors.

Formulate the Local LS Estimators of the UI
Step 3: Define the parameters G k and B k .
Step 4: Use the exchanged information to calculated i k−1 and P dd,i k−1 via Equations (25) and (26).

Refine the Local Sensor Bias Estimators after Consensus
Step 9: Refine the estimates b i k|k and P i,bb k|k using Equations (33) and (34).
Step 10: Set k = k + 1 and go to Step 1.

Remark 5.
In this paper, a fusion algorithm for multiple cooperative sensors subject to interference from GSBs is proposed. However, when this algorithm is applied to continuous target tracking in a wireless sensor networks (WSN), a new problem should be considered. Because of the sensors' limited detection capabilities and low energy storage, a WSN is divided into many clusters. Each cluster forms a sub-WSN and is responsible for tracking the target for a certain period of time. When the target moves out of the surveillance area of the current cluster, one or more new clusters must be generated dynamically; thus, a strategy for generating such clusters must be proposed. This strategy must specify how to choose the leader and followers in each cluster and how to transmit estimate information between clusters. This problem will be an interesting topic for future research regarding the joint optimization of signal processing and distributed cooperative detection. Currently, this topic has been tried in our research group [44].

Remark 6.
For the case of multiple-target tracking with dense clutter, the association relationships between the targets and measurements are unknown. Therefore, the problem of data association must be considered when our proposed algorithm is applied for multiple-target tracking. To this end, the optimal estimation for a single target must be extended to the joint optimization of state estimation and target identification for multiple targets. However, because of the limited detection capabilities of nodes in a WSN, as also mentioned in Remark 5, a WSN is often divided into many clusters, each of which monitors only a relatively small region. Consequently, the probability of multiple targets being present in any single cluster is statistically low. Therefore, our algorithm is focused only on single-target tracking.

Numerical Example
Consider examples of distributed networks tracking one constant-velocity target. The target state is [x k ,ẋ k , y k ,ẏ k ] T , where [x k ,ẋ k ] T and [y k ,ẏ k ] T represent the target's Cartesian position and velocity on the x and y axes, respectively. The sampling interval is T = 1 s, and the number of the samples is 60. Each sensor node synchronously measures the position of the target at each sampling instant. The measurements of each sensor node are corrupted by a local dynamic system bias, and all local system biases originate from a common UI. This situation can be modelled using Equations (1) and (3) with the following parameters:

A Distributed Network of Twelve Sensors
Consider an example of a distributed network of twelve sensors, and the matrices N i k for the 12 sensor nodes are The dimensions of the surveillance area are 16 km × 14 km. The topological map is depicted in Figure 4a, and the sensors' true locations and the target trajectory are depicted in Figure 4b.
The UI, plotted in Figure 5, is The histograms shown in Figure 6 depict the numbers of steps of iteration for different thresholds. Clearly, the number of iterations required for consensus increases with a decreasing threshold. Moreover, when the threshold is greater than 0.1, the number of iteration steps varies only slowly. However, when the threshold is less than 0.1, the number of iteration steps increases sharply, and the computational cost is dramatically increased. Considering both calculation and accuracy, the consensus threshold values for both target state and UI estimation are chosen to be η x = η d =0.1.   In accordance with Theorem 2, the rank of the filter matrix is rank H i k+1Ḡ i k =rank Ḡ i k = 1, which satisfies Equation (10).
The proposed algorithm is compared with the EM algorithm [25] in this simulation. The iteration threshold for the EM algorithm is δ = 10 −4 , and the processing window length is l = 5. We also apply the standard Kalman filter, neglecting the presence of a UI. The estimation of the Kalman filter shows that the fusion performance will deteriorate considerably when bias exists in the sensor measurements.
Because of the page limitation and the similarity in the data variation trends among the simulated sensors, only the estimation results for the fourth sensor are presented here. Figures 7 and 8 show the target-state estimation errors for the Kalman filter (called 'KF'), the EM method, and the proposed method without consensus (called 'DP no consensus') and with consensus (called 'DP consensus'). It can be seen that in the DP methods, the UI is decoupled from the system. Furthermore, the estimation error of the DP consensus method is less than that of the DP no consensus method, and the estimation accuracy is also improved for the DP consensus method. It is also observed that the estimation results based on the KF and EM methods are not as good as those achieved using the DP methods. Two reasons for these findings can be identified. One is that the UI is ignored in the KF method, causing the KF algorithm to become invalid. The other is that the UI is assumed to take a constant value in the iterative process of the EM method. Similar conclusions can also be drawn from the results of the sensor bias estimation, as shown in Figure 9.  The estimation errors for the UI based on the DP methods with and without consensus are plotted in Figure 10. It can be seen that the UI estimation error is greatly reduced by the consensus processing.   Tables 3 and 4, respectively. It can be seen from Table 3 that the mean estimation error values of the DP consensus method are all smaller than those of the KF, EM and DP no consensus methods. As shown in Table 4, the peak estimation error values are also smallest for the DP consensus must. Thus, the proposed method is proven to be effective.
The run times are also compared as listed in Table 5. It can be seen that the run time for the DP no consensus method is merely 0.33 s because of its non-iteration operation. After iteration, the run time of the DP consensus method is 4.89 s, which is nevertheless shorter than the 6.57 s required for the EM method.

A Distributed Network of Sixty Sensors
Now we consider an example of distributed network of sixty sensors, and the matrices N i k of the sensor nodes are The dimensions of the surveillance area are 18 km × 16 km. The topological map is shown in Figure 11a. Sensors' true locations and the target trajectory are depicted in Figure 11b.
In accordance with Theorem 2, the rank of filter matrix rank H i k+1Ḡ i k =rank Ḡ i k = 1 for each sensor in the network, which satisfies Equation (10).
The UI, as plotted in Figure 12, is a stochastic noise obeying Gaussian distribution with mean 0 and covariance 15 2 .  The histograms in Figure 13 depict the mean numbers of consensus iteration with different thresholds. Similarly, the consensus iteration number increases with thresholds decreasing, and when the thresholds are less than 0.1, the iteration step number increase sharply. Therefore, we choose the threshold as η x = η d = 0.1. It is interesting to see from the histograms that the consensus iteration needed for UI is roughly one more time than that needed for sensor states. Due to the page limitation and the similar data variation trend in simulation figures, only the estimation results for the forth sensor is presented here. The iteration threshold for EM algorithm is δ = 10 −4 and the processing window length is l = 5 too. Figures 14-17 show the proposed algorithm estimation results are obviously better than those of KF and EM methods. Indeed, The estimation errors are greatly reduced by the DP consensus operation.
The mean and peak values of the estimation errors of the target state [x,ẋ, y,ẏ] T and the bias [b x , b y ] T based on the KF, EM and DP methods are correspondingly computed and listed in Tables 6  and 7, respectively. It can be seen from the Tables that our proposed method is proven to be effective.
The run times listed in Table 8 show that the run time of DP consensus is slightly shorter than that of EM method. DP no consensus DP consensus Figure 17. UI estimation errors of the DP consensus and DP no consensus methods for sensor 4 in the multi-sensor system of 60 sensors.

Computational Burden Analysis
Computation burden is also a critical issue to determine the practical applicability of the proposed method and worth to be discussed. Generally speaking, the computational burden is strongly related to the total sensors' number and the connectivity of the sensor network. Thus, the mean values of run times of DP consensus method are tested with different sensor numbers and connectivity probabilities. In the simulation, the network are generated randomly with different connective probabilities. And the simulation program runs ten times with every connective probability and sensor number. The mean values of run time for one sensor of processing 60 samples are recorded to measure the computational burden.
According to Figures 6 and 13, the thresholds are chosen as η x = η d =0.1 in this simulation. Further, Figure 18 depicted the relationship between the mean of the run time and the number of sensor at different connectivity probabilities ('CP' in Figure 18). It can be seen from Figure 18 that the run time increases with sensor number increase. It means more iterations are needed to reach a consensus when sensor number increase. Particularly, when the number of the sensor are 12 and 60 with connectivity probability of 1, the run times reach 0.08 s and 0.28 s, respectively. Furthermore, when the sensor number fixed, the run time increases with the connectivity probability decreases. This is because that a looser sensor network topology means fewer neighboring sensors, therefore, more jumpers are needed between each pair of sensors. So more iterations are needed to reach consensus, and the computational cost is corresponding high. It is also worth to mention that when the connectivity probability is larger than 0.2, run time almost have no difference. It implies that our proposed algorithm is somewhat robust when the connectivity probability is larger than 0.2.
As Figure 18 shows, with the sensor number increases, the mean values of run time for one sensor increases linearly. Thus, our proposed method can handle a relatively large network in a reasonable period of time.

Conclusions
The joint estimation of system state and sensor bias based on a generalized system model is proposed. The conditions for UI decoupling derived in this paper are helpful for guiding sensor choice and deployment. Local UI estimates are obtained via the LS approach using all estimates from neighboring sensors. Network consensus processing is then applied to obtain more accurate estimates of the system (target) state and the UI. Future work may address distributive and collaborative processing in sensor networks in the case in which bias exists in only a subset of the sensors. GLSB Global sensor bias ACF Average consensus fusion DLS Distributed least squares WSN Wireless sensor networks KF Kalman Filter DP no consensus The proposed method without consensus processing DP consensus The proposed method with consensus processing CP Connectivity probability