Next Article in Journal
Real Time 3D Facial Movement Tracking Using a Monocular Camera
Next Article in Special Issue
Fuzzy Adaptive Cubature Kalman Filter for Integrated Navigation Systems
Previous Article in Journal
Test-Retest Reliability of an Automated Infrared-Assisted Trunk Accelerometer-Based Gait Analysis System
Previous Article in Special Issue
The Vector Matching Method in Geomagnetic Aiding Navigation
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

State Estimation for a Class of Non-Uniform Sampling Systems with Missing Measurements

School of Electronics Engineering, Heilongjiang University, Harbin 150080, China
*
Author to whom correspondence should be addressed.
Sensors 2016, 16(8), 1155; https://doi.org/10.3390/s16081155
Submission received: 10 June 2016 / Revised: 13 July 2016 / Accepted: 18 July 2016 / Published: 23 July 2016
(This article belongs to the Special Issue Advances in Multi-Sensor Information Fusion: Theory and Applications)

Abstract

:
This paper is concerned with the state estimation problem for a class of non-uniform sampling systems with missing measurements where the state is updated uniformly and the measurements are sampled randomly. A new state model is developed to depict the dynamics at the measurement sampling points within a state update period. A non-augmented state estimator dependent on the missing rate is presented by applying an innovation analysis approach. It can provide the state estimates at the state update points and at the measurement sampling points within a state update period. Compared with the augmented method, the proposed algorithm can reduce the computational burden with the increase of the number of measurement samples within a state update period. It can deal with the optimal estimation problem for single and multi-sensor systems in a unified way. To improve the reliability, a distributed suboptimal fusion estimator at the state update points is also given for multi-sensor systems by using the covariance intersection fusion algorithm. The simulation research verifies the effectiveness of the proposed algorithms.

1. Introduction

Recently, the estimation problems for multi-rate non-uniform sampling systems have attracted much attention due to wide applications in parameter identification [1], industrial detection [2], target tracking and signal processing [3,4,5,6,7,8]. Differently from the single-rate uniform sampling systems, the design on multi-rate non-uniform sampling systems is more complex and challenging.
For multi-rate non-uniform sampling systems, the multi-sensor information fusion filters have been presented based on the data block method in [3,4,5] where the effect of the system noise on modeling is ignored, which brings the errors in modeling. By considering the system noise in modeling, an optimal estimator in the linear minimum variance sense is presented in [6]. The modeling methods in [4,5,6] are, respectively, adopted to obtain the distributed fusion filters by the weighting sums of the local filters in [7,8]. In the literature above, the state models are all from weighted average of states in a data block. When there are multiple measurement samples within a state update period, a state estimator is designed by the measurement augmentation method in [9], which brings expensive computational cost. The discretization of continuous systems is adopted in the multi-rate processing in [10,11].
The distributed fusion filters are designed for multi-rate multi-sensor systems in [12,13]. In [12], an estimator is proposed for systems with uncertain observations. In [13], a “dummy” measurement is used to transform multi-rate into single-rate for each sensor and the packet dropouts are considered when the measurements are transmitted by networks. For the multi-sensor networked systems with packet dropouts, a two-stage distributed fusion estimation algorithm is proposed by using a multi-rate scheme to reduce communication cost in [14], where sensors collect measurements from their neighbors to generate their own local estimates, and then local estimates are collected to form a fused estimate. The multi-rate H filtering problem for the norm-bounded uncertain systems with packet dropouts is investigated by the state augmentation method in [15]. The multi-rate filter in the least mean square sense is designed in [16] and can obtain more accurate estimate than the H filter in [15]. In these studies, the sampling periods of individual sensors are uniform and integer times of the state update period though different sensors have different sampling rates. The two-sensor multi-rate fusion estimation problem for wireless sensor networks is investigated in [17]. Then, the results in [17] are improved and extended in [18]. However, the sampling period of each sensor is still positive integer times of the state update period. A novel method that is named direct estimation of sampling time in solving asynchronous track-to-track fusion problem in [19] is used to predict the pseudo-synchronized state estimates of all the sensors that possess their own sampling rates for the start of the next fusion period. A suboptimal hierarchical fusion estimator is designed in [20] for a clustered sensor network by using covariance intersection fusion algorithm, where local estimators and the fusion center are allowed to be asynchronous. Information fusion estimation problem in multi-sensor networked systems is also explored in [21,22,23,24] where packet dropouts, random delays and missing measurements are considered. However, the multi-rate asynchronous estimation problems are not taken into account.
In practice, not only different sensors can have different sampling rates but also the same sensor can also have asynchronous sampling rates. Moreover, the state estimators need to be obtained in real time for many practical applications. Early work has been done in [25] where the real time estimation can be obtained not only at the state update points but also at the measurement sampling points. However, the single sensor optimal estimation problem is only considered in [25] in the case that the sensor randomly samples the measurement once at most in a state update period. Motivated by the above discussion, we carry out our work in this paper. For systems with the uniformly updated state and randomly sampled measurements, where missing measurements are also considered, a new state space model is developed to depict the dynamics at the measurement sampling points within a state update period by weighting the endpoint states of a state update period. Differently from [25] where sensor randomly samples one measurement at most in a state update period, sensor randomly samples one or multiple measurements or nothing in a state update period in this paper. A non-augmented recursive estimator at the measurement sampling points and at the state update points is designed by applying projection theory. Compared with the augmented method [9], the proposed algorithm can reduce the computational burden and provide the estimates at the measurement sampling points within a state update period. At last, the optimal fusion and distributed suboptimal fusion estimators are given to deal with the fusion estimation problems for multi-sensor systems.

2. Problem Formulation

Consider the following non-uniform sampling discrete time-invariant linear stochastic system with missing measurements.
x ( k + 1 ) = Φ x ( k ) + Γ w ( k )
y ( k j ) = ξ ( k j ) H x ( k j ) + v ( k j )
where x ( k ) R p is the system state at the moment k T , the Italic T is the state sampling period, and y ( k j ) R q is the measurement of the sensor at the moment k j . x ( k j ) is the state measured by y ( k j ) . w ( k ) R r and v ( k j ) R q are the system noise and measurement noise. Φ , Γ and H are constant matrices with suitable dimensions. k j is the sampling time of the jth measurement. The variable ξ ( k j ) is a Bernoulli distributed stochastic variable that takes values on 1 and 0 with the probability Prob { ξ ( k j ) = 1 } = γ ¯ , γ ¯ [ 0 , 1 ] . ξ ( k j ) = 1 means that the measurement is received, while ξ ( k j ) = 0 means that the measurement is missing. In the whole text, E denotes the mathematical expectation and the superscript Roman T denotes the transpose. We easily obtain the results:
E { ξ ( k j ) } = γ ¯ ,   E { ( ξ ( k j ) γ ¯ ) 2 } = γ ¯ ( 1 γ ¯ ) ,   E { ξ ( k j ) ( 1 ξ ( k j ) ) } = 0 , E { ξ ( k i ) ξ ( l j ) } = γ ¯ 2 ( k i l j )
For brevity of notations, only linear time-invariant systems are considered. However, the results derived later can be easily extended to linear time-varying systems. We adopt the following standard assumptions.
Assumption 1. 
w ( k ) and v ( k j ) are uncorrelated white noises with zero means and covariance matrices Q w and Q v .
Assumption 2. 
The initial state x ( 0 ) is uncorrelated with w ( k ) , v ( k j ) and ξ ( k j ) , and satisfies
E { x ( 0 ) } = μ 0 ,   E { ( x ( 0 ) μ 0 ) ( x ( 0 ) μ 0 ) T } = P 0
Assuming that the sampling time of the sensor is known, we consider a class of non-uniform sampling discrete-time systems with missing measurements where the state updates uniformly and the sensor samples randomly in this paper. An example of sampling time versus sensor map is depicted in Figure 1. It can be seen from Figure 1 that there are two cases when the sensor samples the measurement data. Case 1 is that there are the measurement samples in a state update period ( ( k 1 ) T , k T ] , for example, only one measurement sample at the time instants T, 2T, 3T, 4T, 9T, 10T and within the interval ( 4 T , 5 T ] with the sampling time instants k 1 = 1 , k 2 = 2 , k 3 = 3 , k 4 = 4 , k 11 = 9 , k 12 = 10 , k 5 = 4.85 , two samples in ( 7 T , 8 T ] with the sampling time instants k 9 = 7.6 , k 10 = 8 , and three samples within ( 5 T , 6 T ] with the sampling time instants k 6 = 5.35 , k 7 = 5.65 , k 8 = 6 . Case 2 is that there are no samples within the state update period ( ( k 1 ) T , k T ] , for example, not any measurement sample within ( 6 T , 7 T ] .
Our aim is to find the state filters at state update points and measurement sampling points based on the received measurements.

3. System Modeling

Assume that the sensor starts to sample the measurement at the initial moment. When the sensor samples a measurement at the state update point, it is natural to take the system state at the state update point as the state at the measurement sampling point. When the sensor samples a measurement between the two state update points, we adopt the method in [9] where the weighted value of the two adjacent states is taken as the state at the measurement sampling point.
In order to facilitate the algorithm, let N k is the number of measurement samples within a state update interval ( ( k 1 ) T , k T ] , and S k = l = 1 k N l is the total number of measurement samples before the moment k T . The i ( 1 i N k ) denotes the ith sample within the interval ( ( k 1 ) T , k T ] . k = k S k 1 + i / T where denotes the minimal integer not less than . Let α i ( k ) = k k S k 1 + i / T , and then we have k S k 1 + i = k α i ( k ) . It is easily known that 0 α i ( k ) < 1 . To simplify the expression, the sampling period T will be omitted under no confusion in the subsequent text.
The state at the measurement sampling point within the interval ( k 1 , k ] is given by
x ( k α i ( k ) ) = ( 1 α i ( k ) ) x ( k ) + α i ( k ) x ( k 1 )
Accordingly, the measurement equation is expressed as
y ( k α i ( k ) ) = ξ ( k α i ( k ) ) H x ( k α i ( k ) ) + v ( k α i ( k ) )
Especially, when the measurement is sampled at the state update point, we have α i ( k ) = 0 . Then the state at the measurement sampling point can be reduced as x ( k ) . The measurement equation can be reduced as y ( k ) = ξ ( k ) H x ( k ) + v ( k ) .
Assuming i and i 1 are two adjacent sampling points in ( k 1 , k ] . From Equation (5), we have
x ( k ) = 1 1 α i 1 ( k ) x ( k α i 1 ( k ) ) α i 1 ( k ) 1 α i 1 ( k ) x ( k 1 )
Substituting Equation (7) into Equation (5) yields
x ( k α i ( k ) ) = ( 1 α i ( k ) ) ( 1 α i 1 ( k ) ) x ( k α i 1 ( k ) ) + ( α i ( k ) α i 1 ( k ) ) ( 1 α i 1 ( k ) ) x ( k 1 )
Next, we construct the state space model at the measurement sampling points within a state update period. The following Theorem 1 gives the result.
Theorem 1. 
Under Assumptions 1 and 2, the state space model at the measurement sampling points within the interval ( k 1 , k ] can be set up as follows:
x ( k α i ( k ) ) = Φ i ( k ) x ( k 1 ) + Γ i ( k ) w ( k 1 )
y ( k α i ( k ) ) = ξ ( k α i ( k ) ) H x ( k α i ( k ) ) + v ( k α i ( k ) )
where the coefficient matrices Φ i ( k ) and Γ i ( k ) are defined by
Φ i ( k ) = m = 0 i 1 β i m ( k ) Φ + l = 0 i 1 ( m = 0 l 1 β i m ( k ) ) ( 1 β i l ( k ) ) I ,   Γ i ( k ) = m = 0 i 1 β i m ( k ) Γ
with β i ( k ) = ( 1 α i ( k ) ) / ( 1 α i 1 ( k ) ) , ( i > 1 ) , β 1 ( k ) = 1 α 1 ( k ) and m = 0 1 β i m ( k ) = 1 .
Proof. 
From Equation (8), we have
x ( k α i ( k ) ) = β i ( k ) x ( k α i 1 ( k ) ) + ( 1 β i ( k ) ) x ( k 1 )
From Equation (5), we obtain
x ( k α 1 ( k ) ) = β 1 ( k ) x ( k ) + ( 1 β 1 ( k ) ) x ( k 1 )
By iterating Equation (12) and using Equation (13), we have
x ( k α i ( k ) ) = β i ( k ) x ( k α i 1 ( k ) ) + ( 1 β i ( k ) ) x ( k 1 ) = β i ( k ) β i 1 ( k ) x ( k α i 2 ( k ) ) + β i ( k ) ( 1 β i 1 ( k ) ) x ( k 1 ) + ( 1 β i ( k ) ) x ( k 1 ) = = β i ( k ) β i 1 ( k ) β 1 ( k ) x ( k ) + β i ( k ) β i 1 ( k ) β 2 ( k ) ( 1 β 1 ( k ) ) x ( k 1 ) + + β i ( k ) ( 1 β i 1 ( k ) ) x ( k 1 ) + ( 1 β i ( k ) ) x ( k 1 ) = m = 0 i 1 β i m ( k ) ( Φ x ( k 1 ) + Γ w ( k 1 ) ) + l = 0 i 1 ( m = 0 l 1 β i m ( k ) ) ( 1 β i l ( k ) ) x ( k 1 ) = Φ i ( k ) x ( k 1 ) + Γ i ( k ) w ( k 1 )
where Φ i ( k ) and Γ i ( k ) are defined by Equation (11). The proof is completed.
Remark 1. 
The model proposed in Theorem 1 establishes the relationship between the state at the measurement sampling points within the interval ( k 1 , k ] . It includes the single rate uniform sampling system as the special case, i.e., N k = 1 and α 1 ( k ) = 0 . Thus, it generalizes the standard single rate discrete-time state space model.
To design the filtering algorithm, we first give a lemma to be used in the subsequent sections.
Lemma 1. 
For system Equation (9) under Assumptions 1 and 2, the state second-order moment matrix q x ( k α i ( k ) ) = E { x ( k α i ( k ) ) x T ( k α i ( k ) ) } is computed by
q x ( k α i ( k ) ) = Φ i k q x ( t 1 ) ( Φ i k ) T + Γ i k Q w ( Γ i k ) T
The initial value is q x ( 0 ) = μ 0 μ 0 T + P 0 .
Proof. 
Substituting Equation (9) into q x ( k α i ( k ) ) = E { x ( k α i ( k ) ) x T ( k α i ( k ) ) } , we easily obtain Equation (15).

4. Optimal State Estimator

In this section, we will present our estimation algorithm at the state update points and measurement sampling points based on the developed model in Theorem 1.

4.1. Estimator Design

Theorem 2. 
For system Equations (9) and (10) under Assumptions 1 and 2, the optimal estimators at the state update point and measurement sampling points in the interval ( k 1 , k ] ( N k 0 ) are computed by
x ^ ( k α i ( k ) | k α i r ( k ) ) = Φ i ( k ) x ^ ( k 1 | k α i r ( k ) ) + Γ i ( k ) w ^ ( k 1 | k α i r ( k ) ) , r = 0 , 1
where it is a filter for r = 0 or a predictor for r = 1 . The fixed-point smoothers for the state and white noise are respectively computed by
x ^ ( k 1 | k α i ( k ) ) = x ^ ( k 1 | k α i 1 ( k ) ) + K x ( k 1 | k α i ( k ) ) ε ( k α i ( k ) )
w ^ ( k 1 | k α i ( k ) ) = w ^ ( k 1 | k α i 1 ( k ) ) + K w ( k 1 | k α i ( k ) ) ε ( k α i ( k ) )
where the innovation sequence ε ( k α i ( k ) ) and its covariance Q ε ( k α i ( k ) ) are computed by
ε ( k α i ( k ) ) = y ( k α i ( k ) ) γ ¯ H x ^ ( k α i ( k ) | k α i 1 ( k ) )
Q ε ( k α i ( k ) ) = γ ¯ 2 H P x ( k α i ( k ) | k α i 1 ( k ) ) H T + γ ¯ ( 1 γ ¯ ) H q x ( k α i ( k ) ) H T + Q v
The smoothing gain matrices K x ( k 1 | k α i ( k ) ) for the state and K w ( k 1 | k α i ( k ) ) for the white noise are computed by
K x ( k 1 | k α i ( k ) ) = γ ¯ [ P x ( k 1 | k α i 1 ( k ) ) Φ i ( k ) T + P x w ( k 1 | k α i 1 ( k ) ) Γ i ( k ) T ] H T Q ε 1 ( k α i ( k ) )
K w ( k 1 | k α i ( k ) ) = γ ¯ [ P w x ( k 1 | k α i 1 ( k ) ) Φ i ( k ) T + P w ( k 1 | k α i 1 ( k ) ) Γ i ( k ) T ] H T Q ε 1 ( k α i ( k ) )
The smoothing error covariance matrices P x ( k 1 | k α i ( k ) ) and P w ( k 1 | k α i ( k ) ) , and the cross-covariance matrix P x w ( k 1 | k α i ( k ) ) for the state and white noise are computed by
P x ( k 1 | k α i ( k ) ) = P x ( k 1 | k α i 1 ( k ) ) K x ( k 1 | k α i ( k ) ) Q ε ( k α i ( k ) ) K x T ( k 1 | k α i ( k ) )
P w ( k 1 | k α i ( k ) ) = P w ( k 1 | k α i 1 ( k ) ) K w ( k 1 | k α i ( k ) ) Q ε ( k α i ( k ) ) K w T ( k 1 | k α i ( k ) )
P x w ( k 1 | k α i ( k ) ) = P x w ( k 1 | k α i 1 ( k ) ) K x ( k 1 | k α i ( k ) ) Q ε ( k α i ( k ) ) K w T ( k 1 | k α i ( k ) )
The filtering ( r = 0 ) and prediction ( r = 1 ) error covariance matrices P x ( k α i ( k ) | k α i r ( k ) ) for the state are computed by
P x ( k α i ( k ) | k α i r ( k ) ) = Φ i ( k ) P x ( k 1 | k α i r ( k ) ) Φ i ( k ) T + Γ i ( k ) P w ( k 1 | k α i r ( k ) ) Γ i ( k ) T + Φ i ( k ) P x w ( k 1 | k α i r ( k ) ) Γ i ( k ) T + Γ i ( k ) P w x ( k 1 | k α i r ( k ) ) Φ i ( k ) T ,   r = 0 , 1
The initial values are P x ( k 1 | k α 0 ( k ) ) = P x ( k 1 | k 1 ) , P w ( k 1 | k α 0 ( k ) ) = Q w , P x w ( k 1 | k α 0 ( k ) ) = 0 , x ^ ( k 1 | k α 0 ( k ) ) = x ^ ( k 1 | k 1 ) and w ^ ( k 1 | k α 0 ( k ) ) = 0 .
It is clear that we have the estimator at the state update point as x ^ ( k | k ) = x ^ ( k α N k ( k ) | k α N k ( k ) ) when α N k ( k ) = 0 or x ^ ( k | k ) = x ^ ( k α N k + 1 ( k ) | k α N k ( k ) ) with α N k + 1 ( k ) = 0 , and β N k + 1 ( k ) = 1 / ( 1 α N k ( k ) ) when α N k ( k ) 0 .
Proof. 
Remark 2. 
In Theorem 2, if N k = 1 and α 1 ( k ) = 0 , we have the filtering algorithm at the state update points for single-rate systems. If N k = 0 , i.e., no samples in the interval ( k 1 , k ] , the predictor is used to estimate the state at the state update point k based on the estimator at the state update point k 1 .
Remark 3. 
It can be easily known that the proposed non-augmented recursive estimator has the computational order of magnitude O ( N k p 3 ) in the interval ( k 1 , k ] . Compared with the augmented estimator [9] with the computational order of magnitude O ( N k 3 q 3 ) , our algorithm can obviously reduce the computational cost with the increase of the number N k of measurement samples in the interval ( k 1 , k ] for a deterministic system with the fixed p and q. Meanwhile, it is more important that our filter can give the estimates not only at the state update points but also at the measurement sampling points within a state update period in real time. However, the estimator of [9] can only give the estimates at the state update points.

4.2. Computational Procedures of the Estimator

Based on Theorem 1 and Theorem 2, the computational procedures of our estimator at the state update points and at the measurement sampling points are addressed as follows:
Step 1. k = 1 , set the initial values x ^ ( 0 | 0 ) = μ 0 , P x ( 0 | 0 ) = P 0 , P w ( 0 | 0 ) = Q w , P x w ( 0 | 0 ) = 0 and q x ( 0 ) = μ 0 μ 0 T + P 0 .
Step 2. Construct the model within a state update period ( k 1 , k ] by Theorem 1.
Step 3. Compute the state estimators at the state update points and at the measurement sampling points:
(a)
If N k 0 (i.e., there are samples within a state update period), obtain the state estimates x ^ ( k α i ( k ) | k α i ( k ) ) and the covariance matrices P x ( k α i ( k ) | k α i ( k ) ) by Theorem 2.
(b)
If N k = 0 (i.e., no sample within a state update period), obtain the estimate by prediction method in Remark 2, i.e., x ^ ( k | k ) = Φ x ^ ( k 1 | k 1 ) , P x ( k | k ) = Φ P x ( k 1 | k 1 ) Φ T + Γ Q w Γ T .
Step 4. k = k + 1 , set P x ( k 1 | k α 0 ( k ) ) = P x ( k 1 | k 1 ) , P w ( k 1 | k α 0 ( k ) ) = Q w , P x w ( k 1 | k α 0 ( k ) ) = 0 , x ^ ( k 1 | k α 0 ( k ) ) = x ^ ( k 1 | k 1 ) and w ^ ( k 1 | k α 0 ( k ) ) = 0 .
Go to Step 2.

5. Multi-Sensor Case

In the preceding section, we have presented a non-augmented optimal recursive estimator for single sensor system. In this section, we will discuss how to use the proposed algorithm to solve the multi-sensor case.
For a multi-sensor system, the state equation is the same as Equation (1) and the measurement equations are given as follows:
y ( l ) ( k j ( l ) ) = ξ ( l ) ( k j ( l ) ) H ( l ) x ( k j ( l ) ) + v ( l ) ( k j ( l ) )
where the superscript l denotes the lth sensor, L is the number of sensors, y ( l ) ( k j ( l ) ) is the measurement at the sampling time k j ( l ) for the lth sensor, k j ( l ) is the sampling time of the jth measurement, v ( l ) ( k j ( l ) ) is the measurement noise with zero mean and covariance Q v ( l ) , and H ( l ) is the measurement matrix. The variable ξ ( l ) ( k j ( l ) ) is a Bernoulli distributed stochastic variable that takes values on 1 and 0 with the probability Prob { ξ ( l ) ( k j ( l ) ) = 1 } = γ ¯ ( l ) , γ ¯ ( l ) [ 0 , 1 ] . We assume that ξ ( l ) ( k j ( l ) ) is independent of w ( k ) and v ( l ) ( k j ( l ) ) and x ( 0 ) . Moreover ξ ( l ) ( k j ( l ) ) satisfies that
E { ξ ( l ) ( k j ( l ) } = γ ¯ ( l ) ,   E { ( ξ ( l ) ( k j ( l ) γ ¯ ( l ) ) 2 } = γ ¯ ( l ) ( 1 γ ¯ ( l ) ) ,   E { ξ ( l ) ( k j l ) γ ¯ ( l ) } = 0 , E { ξ ( l ) ( k j ( l ) ) ( 1 ξ ( l ) ( k j ( l ) ) ) } = 0 ,   E { ξ ( l ) ( k i ( l ) ) ξ ( s ) ( m j ( s ) ) } = γ ¯ ( l ) γ ¯ ( s ) ( l s  or  k i ( l ) m j ( s ) )
We will adopt two types of methods to deal with the estimation problem of multi-sensor case: optimal fusion estimator and suboptimal fusion estimator.

5.1. Optimal Fusion Estimator

In fact, we can reorder the measurements y ( l ) ( k j ( l ) ) of all sensors in each state update period ( k 1 , k ] according to the order of sampling time k j ( l ) . If the sampling time of measurements from some sensors is same, i.e., sampling at the same time point, they will be combined to an augmented measurement at this sampling time. These reorder measurements can be considered from a certain single sensor. Then, our estimation algorithm in Section 4 can be applied to the optimal fusion estimation problem of multiple sensors.
When all sensors work healthily, the proposed optimal fusion algorithm can obtain the best accuracy. However, it has bad reliability, which means that a faulty sensor can lead to the optimal fusion estimator to diverge [26]. To improve the reliability, we give the following distributed suboptimal fusion algorithm.

5.2. Suboptimal Fusion Estimator

Firstly, apply our estimation algorithm in this paper to obtain the local estimators at the state update points based on the measurements from individual sensors, and then apply the covariance intersection (CI) fusion algorithm [27,28] to obtain the distributed suboptimal fusion estimator at the state update points. The reason that we adopt CI algorithm is that it avoids the computation of cross-covariance matrices between any two local estimators [26]. The CI fusion estimator can be given as follows [27,28]:
x ^ CI ( k | k ) = l = 1 L ω ( l ) ( k ) P CI ( k | k ) [ P x ( l ) ( k | k ) ] 1 x ^ ( l ) ( k | k )
P CI ( k | k ) = [ l = 1 L ω ( l ) ( k ) [ P x ( l ) ( k | k ) ] 1 ] 1
where P CI ( k | k ) is the upper bound of variance of the CI fusion estimator. The optimal weighted coefficients ω ( l ) ( k ) satisfy 0 ω ( l ) ( k ) 1 and l = 1 L ω ( l ) ( k ) = 1 , and can be obtained by solving the following optimization problem:
min l = 1 L ω ( l ) ( k ) = 1 , 0 ω ( l ) ( k ) 1 tr ( P C I ( k | k ) )
The nonlinear optimization problem Equation (29) does not generally have the analytical solutions. It can be solved by numerical methods, which can be carried out by the function “fmincon” in Matlab optimization toolbox. It only involves the scalar optimization problem.
The proposed distributed CI fusion estimator has the advantage of good robustness and flexibility, i.e., good reliability since it has the distributed parallel structure that is convenient for detection and isolation of a faulty sensor. Certainly, when all sensors work healthily, it has worse accuracy than the optimal fusion estimator. However, it has better accuracy than local estimators [27,28].
Remark 4. 
The proposed optimal fusion algorithm is a non-augmented method. Certainly, we can also use the centralized fusion way [9] that combines the measurements from all sensors within a state update period into an augmented measurement and then use standard Kalman filter to solve. However, it is well known that this way has the expensive computational burden with the increase of the number of sensors [26].

6. Simulation Research

In this section, a spring-mass system that has been widely adopted in many studies [29,30,31] will be used to verify the effectiveness of our algorithms:
x ˙ ( t ) = [ 0 0 1 0 0 0 0 1 k 1 + k 2 M 1 k 2 M 1 μ M 1 0 k 2 M 2 k 2 M 2 0 μ M 2 ] x ( t ) + [ 0 0 1 M 1 1 M 2 ] w ( t )
where x ( t ) = [ x 1 ( t ) x 2 ( t ) x ˙ 1 ( t ) x ˙ 2 ( t ) ] T , x 1 and x 2 are the positions of mass M 1 and mass M 2 , respectively; k 1 and k 1 are the spring constants of spring 1 and spring 2, respectively; and μ denotes the viscous friction coefficient between the mass blocks and the horizontal surface. Moreover, the measurement equations are given as Equation (27) with L = 3. v ( i ) ( k j ( i ) ) , i = 1 , 2 , 3 are independent Gaussian noises with zero means and variances Q v ( i ) and uncorrelated with the process noise w ( k ) with zero-mean and variance Q w . We set Q w = 2 , Q v ( 1 ) = 2 , Q v ( 2 ) = 1 , Q v ( 3 ) = 3 , H ( 1 ) = [ 1 1 1 0 ] , H ( 2 ) = [ 0 1 0 0 ] , H ( 3 ) = [ 0 0 1 1 ] , γ ¯ ( 1 ) = 0.7 , γ ¯ ( 2 ) = 0.9 , γ ¯ ( 3 ) = 0.8 , M 1 = 1 , M 2 = 0.5 , k 1 = 1 , k 2 = 1 , μ = 0.5 , and the sampling period T = 0.1 s , we obtain the following system parameter matrices by discretization:
Φ = [ 0.9902 0.0049 0.0972 0.0002 0.0096 0.9903 0.0003 0.0948 0.1941 0.0969 0.9416 0.0047 0.1891 0.1894 0.0095 0.8955 ] ,   Γ = [ 0.0049 0.0097 0.0975 0.1900 ]
We set the initials as x ( 0 ) = 0 and P 0 = 0.1 I 4 . We take 100 sampling data.
Figure 2 gives the sketch map of the sampling data of the three sensors within the time interval 0–20. Figure 3 gives the tracking performance of the optimal fusion estimator, where the solid curves denote the true values and the dashed ones denote the estimators. Figure 4 gives the comparisons curves of the traces of variances for local estimators and the fusion estimators. From Figure 4, we can see that our optimal fusion estimator has best accuracy and CI fusion estimator has better accuracy than local estimators.
Then Monte Carlo simulation is carried out to compare the algorithms in our paper and [9]. The mean square errors (MSEs): MSE i , s ( k ) = 1 N j = 1 N ( x i j ( k ) x ^ i , s j ( k | k ) ) 2 , i=1, 2, 3, 4; s=our paper, [9]) will be used to evaluate the estimation performance of the two estimation algorithms, where the subscript i denotes the ith components and s denotes our filter or that in [9], respectively, and N is the number of Monte Carlo runs.
In Figure 5, the sensor 1 is randomly employed to compare the algorithms in our paper and [9]. The comparison of MSEs simulated by 100 Monte Carlo runs is shown in Figure 5. Because Yan, L. et al. [9] does not consider the effect of missing measurements, from Figure 5, we can see that the proposed algorithm has better accuracy than the one in [9] when the missing measurements occur in sensors.
Next, we will make the comparison with [9] when there are no missing measurements. The sensor 1 is employed with γ ¯ ( 1 ) = 1 in simulation. Figure 6 gives the comparison of MSEs simulated by 100 Monte Carlo runs for the estimators in this paper and [9]. In Figure 6, the estimates at the state update points and the measurement sampling points within a state update period are all shown by using the our proposed algorithm. Moreover, from Figure 6, we can see that the proposed algorithm has the same accuracy as the augmented one in [9] at the state update points when there are no missing measurements. It is also significant from Remark 3 that our estimator has smaller computational burden than [9]. To compare the computation time between the two algorithms, we use the “unifrnd” function in Matlab to sample 2, 3, 4, and 5 measurements within a state update period, respectively, and give the average of 50 run time of the two algorithms in Table 1. The used computer is Lenovo with the CPU speed of 3.20 GHz (Intel (R) Core (TM) i5-3470 CPU) and the used Matlab version is Matlab R2012a. From Table 1, we can see that our proposed algorithm can save the run time compared with that in [9] with the increase of the number of sampling points within a state update period.

7. Conclusions

A non-augmented recursive estimator has been designed for a class of non-uniform sampling systems with missing measurements, where the system state is updated uniformly and the measurements are sampled randomly. A state space model is constructed to depict the dynamics at the measurement sampling points within a state update period. The proposed estimator dependent on the missing rate can provide the state estimates not only at the state update points but also at the measurement sampling points within a state update period. Compared with [9], our algorithm has better accuracy when there are missing measurements and the same accuracy at the state update points when there are no missing measurements. For multi-sensor systems with missing measurements, our algorithm can also deal with the optimal fusion estimation by reordering the measurements from all sensors. A distributed suboptimal fusion estimator is also given using the covariance intersection (CI) fusion algorithm.

Acknowledgments

This work was supported by the National Natural Science Foundation of China under Grant No. 61174139, 61573132, the Heilongjiang Province Outstanding Youth fund (No. JC201412), Chang Jiang Scholar Candidates Program for Provincial Universities in Heilongjiang (No. 2013CJHB005), and the Program for Graduate Innovation Scientific Research of Heilongjiang University (YJSCX2015-001HLJU).

Author Contributions

Shuli Sun proposed the algorithm and wrote the paper. Honglei Lin derived the algorithm and performed the simulation work.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A. The Proof of Theorem 2

By using the projection theory [32], we readily have Equations (16)–(19), where the filtering gain matrices for the state and input white noise are, respectively, defined as
K x ( k 1 | k α i ( k ) ) = E { x ( k 1 ) ε T ( k α i ( k ) ) } Q ε 1 ( k α i ( k ) )
K w ( k 1 | k α i ( k ) ) = E { w ( k 1 ) ε T ( k α i ( k ) ) } Q ε 1 ( k α i ( k ) )
Substituting Equation (10) into Equation (19), the innovation ε ( k α i ( k ) ) can be rewritten as
ε ( k α i ( k ) ) = ( ξ ( k α i ( k ) ) γ ¯ ) H x ( k α i ( k ) ) + γ ¯ H x ˜ ( k α i ( k ) | k α i 1 ( k ) ) + v ( k α i ( k ) )
where the prediction error x ˜ ( k α i ( k ) | k α i 1 ( k ) ) = x ( k α i ( k ) ) x ^ ( k α i ( k ) | k α i 1 ( k ) ) . Substituting Equation (A3) into the innovation sequence covariance Q ε ( k α i ( k ) ) = E { ε ( k α i ( k ) ) ε T ( k α i ( k ) ) } leads to Equation (20).
Subtracting Equation (16) with r = 1 from Equation (9) yields the prediction error equation
x ˜ ( k α i ( k ) | k α i 1 ( k ) ) = Φ i ( k ) x ˜ ( k 1 | k α i 1 ( k ) ) + Γ i ( k ) w ˜ ( k 1 | k α i 1 ( k ) )
where the smoothing errors x ˜ ( k 1 | k α i 1 ( k ) ) = x ( k 1 ) x ^ ( k 1 | k α i 1 ( k ) ) and w ˜ ( k 1 | k α i 1 ( k ) ) = w ( k 1 ) w ^ ( k 1 | k α i 1 ( k ) ) . Using Equations (A3) and (A4), we have
E { x ( k 1 ) ε T ( k α i ( k ) ) } = γ ¯ E { x ( k 1 ) ( Φ i ( k ) x ˜ ( k 1 | k α i 1 ( k ) ) + Γ i ( k ) w ˜ ( k 1 | k α i 1 ( k ) ) ) T } H T = γ ¯ [ P x ( k 1 | k α i 1 ( k ) ) Φ i ( k ) T + P x w ( k 1 | k α i 1 ( k ) ) Γ i ( k ) T ] H T
where x ( k 1 ) = x ^ ( k 1 | k α i 1 ( k ) ) + x ˜ ( k 1 | k α i 1 ( k ) ) , x ^ ( k 1 | k α i 1 ( k ) ) x ˜ ( k 1 | k α i 1 ( k ) ) and x ^ ( k 1 | k α i 1 ( k ) ) w ˜ ( k 1 | k α i 1 ( k ) ) have been used, where the symbol denotes orthogonality. Substituting Equation (A5) into Equation (A1) yields Equation (21). Similarly, we prove Equation (22) holds.
Subtracting Equation (17) from x ( k 1 ) yields the smoothing error equation
x ˜ ( k 1 | k α i ( k ) ) = x ˜ ( k 1 | k α i 1 ( k ) ) K x ( k 1 | k α i ( k ) ) ε ( k α i ( k ) )
The smoothing error covariance matrix is derived as follows:
P x ( k 1 | k α i ( k ) ) = E { x ˜ ( k 1 | k α i ( k ) ) x ˜ T ( k 1 | k α i ( k ) ) } = P x ( k 1 | k α i 1 ( k ) ) + K x ( k 1 | k α i ( k ) ) Q ε ( k α i ( k ) ) K x T ( k 1 | k α i ( k ) ) E { x ˜ ( k 1 | k α i 1 ( k ) ) ε T ( k α i ( k ) ) } K x T ( k 1 | k α i ( k ) ) K x ( k 1 | k α i ( k ) ) E { ε ( k α i ( k ) ) x ˜ T ( k 1 | k α i 1 ( k ) ) }
Using Equation (A1) and note that
E { x ˜ ( k 1 | k α i 1 ( k ) ) ε T ( k α i ( k ) ) } = E { x ( k 1 ) ε T ( k α i ( k ) ) } = K x ( k 1 | k α i ( k ) ) Q ε ( k α i ( k ) )
Substituting Equation (A8) into Equation (A7) leads to Equation (23). Similarly, we can prove Equations (24) and (25) hold.
Subtracting Equation (16) with r = 0 from Equation (9) yields the filtering error equation
x ˜ ( k α i ( k ) | k α i ( k ) ) = Φ i ( k ) x ˜ ( k 1 | k α i ( k ) ) + Γ i ( k ) w ˜ ( k 1 | k α i ( k ) )
Substituting Equations (A4) and (A9) into P x ( k α i ( k ) | k α i r ( k ) ) = E { x ˜ ( k α i ( k ) | k α i r ( k ) ) x ˜ T ( k α i ( k ) | k α i r ( k ) ) } , r = 0 , 1 leads to Equation (26). This proof is completed.

References

  1. Ding, F.; Liu, G.; Liu, X. Partially coupled stochastic gradient identification methods for non-uniformly sampled systems. IEEE Trans. Autom. Control 2010, 55, 1976–1981. [Google Scholar] [CrossRef]
  2. Li, W.; Shah, S.; Xiao, D. Kalman filters in non-uniformly sampled multirate systems: For FDI and beyond. Automatica 2008, 44, 199–208. [Google Scholar] [CrossRef]
  3. Yan, L.; Liu, B.; Zhou, D. Asynchronous multirate multisensor information fusion algorithm. IEEE Trans. Aerosp. Electron. Syst. 2007, 43, 1135–1146. [Google Scholar] [CrossRef]
  4. Yan, L.; Zhou, D.; Fu, M.; Xia, Y. State estimation for asynchronous multirate multisensor dynamic systems with missing measurements. IET Signal Process. 2010, 4, 728–739. [Google Scholar] [CrossRef]
  5. Yan, L.; Li, X.; Xia, Y.; Fu, Y. Modeling and estimation of asynchronous multirate multisensor system with unreliable measurements. IEEE Trans. Aerosp. Electron. Syst. 2015, 51, 2012–2026. [Google Scholar] [CrossRef]
  6. Lin, H.; Ma, J.; Sun, S. Optimal state filters for a class of non-uniform sampling systems. J. Syst. Sci. Math. Sci. 2012, 32, 768–779. [Google Scholar]
  7. Mahmoud, M.S.; Emzir, M.F. State estimation with asynchronous multirate multismart sensors. Inform. Sci. 2012, 196, 15–27. [Google Scholar] [CrossRef]
  8. Ma, J.; Lin, H.; Sun, S. Distributed fusion filter for asynchronous multi-rate multi-sensor non-uniform sampling systems. In Proceedings of the 15th International Conference on Information Fusion, Singapore, Singapore, 9–12 July 2012; pp. 1645–1652.
  9. Yan, L.; Xiao, B.; Xia, Y.; Fu, M. State estimation for a kind of non-uniform sampling dynamic system. Int. J. Syst. Sci. 2013, 44, 1913–1924. [Google Scholar] [CrossRef]
  10. Hu, Y.; Duan, Z.; Zhou, D. Estimation fusion with general asynchronous multirate sensors. IEEE Trans. Aerosp. Electron. Syst. 2010, 46, 2090–2101. [Google Scholar] [CrossRef]
  11. Xu, S.; Lin, X.; Fu, M.; Zhao, D. Multi-rate sensor adaptive data fusion algorithm based on incomplete measurement. In Proceedings of the IEEE International Conference on Automation and Logistics, Chongqing, China, 15–16 August 2011; pp. 448–491.
  12. Peng, F.; Sun, S. Distributed fusion estimation for multisensor multirate systems with stochastic observation multiplicative noises. Math. Probl. Eng. 2014, 2014. [Google Scholar] [CrossRef]
  13. Ma, J.; Sun, S. Distributed Fusion Filter for Multi-rate Multi-sensor Systems with Packet Dropouts. In Proceedings of the 10th World Congress on Intelligent Control and Automation, Beijing, China, 6–8 July 2012; pp. 4502–4506.
  14. Zhang, W.; Feng, G.; Yu, L. Multi-rate distributed fusion estimation for sensor networks with packet losses. Automatica 2012, 48, 2016–2028. [Google Scholar] [CrossRef]
  15. Liang, Y.; Chen, T.; Pan, Q. Multi-rate stochastic H filtering for networked multisensor fusion. Automatica 2010, 46, 437–444. [Google Scholar] [CrossRef]
  16. Geng, H.; Liang, Y.; Zhang, X. Linear-minimum-mean-square-error observer for multi-rate sensor fusion with missing measurements. IET Control Theory Appl. 2014, 8, 1375–1383. [Google Scholar] [CrossRef]
  17. Zhang, W.; Liu, S.; Chen, M.; Yu, L. Fusion estimation for two sensors with nonuniform estimation rates. In Proceedings of the 51st IEEE Conference on Decision Control, Maui, HI, USA, 10–13 December 2012; pp. 4083–4088.
  18. Zhang, W.; Liu, S.; Yu, L. Fusion estimation for sensor networks with nonuniform estimation rates. IEEE Trans. Circuits Syst. I Reg. Pap. 2014, 61, 1485–1498. [Google Scholar] [CrossRef]
  19. Talebi, H.; Hemmatyar, A.M.A. Asynchronous track-to-track fusion by direct estimation of time of sample in sensor networks. IEEE Sens. J. 2014, 14, 210–217. [Google Scholar] [CrossRef]
  20. Song, H.; Zhang, W.; Yu, L. Hierarchical fusion in clustered sensor networks with asynchronous local estimates. IEEE Signal Process. Lett. 2014, 21, 1506–1510. [Google Scholar] [CrossRef]
  21. Chen, B.; Zhang, W.; Yu, L. Distributed fusion estimation with missing measurements, random transmission delays and packet Dropouts. IEEE Trans. Autom. Control 2014, 59, 1961–1967. [Google Scholar] [CrossRef]
  22. Ma, J.; Sun, S. Centralized fusion estimators for multisensor systems with random sensor delays, multiple packet dropouts and uncertain observations. IEEE Sens. J. 2013, 13, 1228–1235. [Google Scholar] [CrossRef]
  23. Ma, J.; Sun, S. Optimal linear estimators for systems with random sensor delays, multiple packet dropouts and uncertain observations. IEEE Trans. Signal Process. 2011, 59, 5158–5192. [Google Scholar]
  24. Moayedi, M.; Foo, Y.K.; Soh, Y.C. Adaptive Kalman filtering in networked systems with random sensor delays, multiple packet dropouts and missing measurements. IEEE Trans. Signal Process. 2010, 58, 1577–1588. [Google Scholar] [CrossRef]
  25. Lin, H.; Sun, S. State estimation for a class of non-uniform sampling systems. In Proceedings of the 34th Control Conference, Hangzhou, China, 28–30 July 2015; pp. 2024–2027.
  26. Sun, S.; Deng, Z. Multi-sensor optimal information fusion Kalman filter. Automatica 2004, 40, 1017–1023. [Google Scholar] [CrossRef]
  27. Julier, S.J.; Uhlmann, J.K. General decentralized data fusion with covariance intersection. In Handbook of Multisensor Data Fusion; CRC Press: Boca Raton, FL, USA, 2001; pp. 319–342. [Google Scholar]
  28. Deng, Z. Information Fusion Estimation Theory with Application; Science Press: Beijing, China, 2012. [Google Scholar]
  29. Song, H.; Yu, L.; Zhang, W. Networked H filtering for linear discrete-time systems. Inf. Sci. 2011, 181, 686–696. [Google Scholar] [CrossRef]
  30. Gao, H.; Chen, T. H estimation for uncertain systems with limited communication capacity. IEEE Trans. Autom. Control 2007, 52, 2070–2084. [Google Scholar] [CrossRef]
  31. Ma, J.; Sun, S. Optimal linear estimators for multi-sensor stochastic uncertain systems with packet losses of both sides. Digit. Signal Process. 2015, 37, 24–34. [Google Scholar] [CrossRef]
  32. Anderson, B.D.O.; Moore, J.B. Optimal Filtering; Prentice-Hall: Englewood Cliffs, NJ, USA, 1979. [Google Scholar]
Figure 1. Illustration of non-uniform sampling.
Figure 1. Illustration of non-uniform sampling.
Sensors 16 01155 g001
Figure 2. Samples of sensors 1–3 in [1,20]: (a) samples of the first sensor; (b) samples of the second sensor; and (c) samples of the third sensor.
Figure 2. Samples of sensors 1–3 in [1,20]: (a) samples of the first sensor; (b) samples of the second sensor; and (c) samples of the third sensor.
Sensors 16 01155 g002
Figure 3. Optimal fusion estimator: (a) tracking for the position of M1; (b) tracking for the position of M2; (c) tracking for the velocity of M1; and (d) tracking for the velocity of M2.
Figure 3. Optimal fusion estimator: (a) tracking for the position of M1; (b) tracking for the position of M2; (c) tracking for the velocity of M1; and (d) tracking for the velocity of M2.
Sensors 16 01155 g003
Figure 4. Comparison of the traces of variances of local estimators and fusion estimators.
Figure 4. Comparison of the traces of variances of local estimators and fusion estimators.
Sensors 16 01155 g004
Figure 5. Comparison of MSEs of the estimators in this paper and [9] for sensor 1: (a) MSEs of the position estimators of M1; (b) MSEs of the position estimators of M2; (c) MSEs of the velocity estimators of M1; and (d) MSEs of the velocity estimators of M2.
Figure 5. Comparison of MSEs of the estimators in this paper and [9] for sensor 1: (a) MSEs of the position estimators of M1; (b) MSEs of the position estimators of M2; (c) MSEs of the velocity estimators of M1; and (d) MSEs of the velocity estimators of M2.
Sensors 16 01155 g005
Figure 6. Comparison of MSEs of the estimators in this paper and [9] for sensor 1 without missing measurements: (a) MSEs of the position estimators of M1; (b) MSEs of the position estimators of M2; (c) MSEs of the velocity estimators of M1; and (d) MSEs of the velocity estimators of M2.
Figure 6. Comparison of MSEs of the estimators in this paper and [9] for sensor 1 without missing measurements: (a) MSEs of the position estimators of M1; (b) MSEs of the position estimators of M2; (c) MSEs of the velocity estimators of M1; and (d) MSEs of the velocity estimators of M2.
Sensors 16 01155 g006
Table 1. Average time of 50 runs.
Table 1. Average time of 50 runs.
Number of SamplesOur Algorithm[9]
20.062064 s0.066064 s
30.078872 s0.086958 s
40.092339 s0.112633 s
50.103354 s0.146654 s

Share and Cite

MDPI and ACS Style

Lin, H.; Sun, S. State Estimation for a Class of Non-Uniform Sampling Systems with Missing Measurements. Sensors 2016, 16, 1155. https://doi.org/10.3390/s16081155

AMA Style

Lin H, Sun S. State Estimation for a Class of Non-Uniform Sampling Systems with Missing Measurements. Sensors. 2016; 16(8):1155. https://doi.org/10.3390/s16081155

Chicago/Turabian Style

Lin, Honglei, and Shuli Sun. 2016. "State Estimation for a Class of Non-Uniform Sampling Systems with Missing Measurements" Sensors 16, no. 8: 1155. https://doi.org/10.3390/s16081155

APA Style

Lin, H., & Sun, S. (2016). State Estimation for a Class of Non-Uniform Sampling Systems with Missing Measurements. Sensors, 16(8), 1155. https://doi.org/10.3390/s16081155

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop