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

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.


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 8 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 8 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.

Problem Formulation
Consider the following non-uniform sampling discrete time-invariant linear stochastic system with missing measurements.
xpk`1q " Φxpkq`Γwpkq (1) where xpkq P R p is the system state at the moment kT, the Italic T is the state sampling period, and ypk j q P R q is the measurement of the sensor at the moment k j . xpk j q is the state measured by ypk j q. wpkq P R r and vpk j q P 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 ξpk j q is a Bernoulli distributed stochastic variable that takes values on 1 and 0 with the probability Prob ξpk j q " 1 ( " γ, γ P r0, 1s. ξpk j q " 1 means that the measurement is received, while ξpk j q " 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: Etξpk j qu " γ, Etpξpk j q´γq 2 u " γp1´γq, Etξpk j qp1´ξpk j qqu " 0, Etξpk i qξpl j qu " γ 2 pk i ‰ l j q 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. wpkq and vpk j q are uncorrelated white noises with zero means and covariance matrices Q w and Q v . Assumption 2. The initial state xp0q is uncorrelated with wpkq, vpk j q and ξpk j q, and satisfies 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 ppk´1qT, kTs, for example, only one measurement sample at the time instants T, 2T, 3T, 4T, 9T, 10T and within the interval p4T, 5Ts 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 p7T, 8Ts with the sampling time instantsk 9 " 7.6, k 10 " 8, and three samples within p5T, 6Ts 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 ppk´1q T, kTs, for example, not any measurement sample within p6T, 7Ts.
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.
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  Our aim is to find the state filters at state update points and measurement sampling points based on the received measurements.

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 k N is the number of measurement samples within a state update interval ( )  Our aim is to find the state filters at state update points and measurement sampling points based on the received measurements.

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 ppk´1q T, kTs, and S k " ř k l"1 N l is the total number of measurement samples before the moment kT. The i (1 ď i ď N k ) denotes the ith sample within the interval ppk´1q T, kTs.
where r˚s denotes the minimal integer not less than˚. Let α then we have k S k´1`i " k´α pkq i . It is easily known that 0 ď α pkq i ă 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 pk´1, ks is given by Accordingly, the measurement equation is expressed as Especially, when the measurement is sampled at the state update point, we have α pkq i " 0. Then the state at the measurement sampling point can be reduced as xpkq. The measurement equation can be reduced as ypkq " ξpkqHxpkq`vpkq.
Assuming i and i´1 are two adjacent sampling points in pk´1, ks. From Equation (5), we have Substituting Equation (7) into Equation (5) yields 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 pk´1, ks can be set up as follows: where the coefficient matrices Φ pkq i and Γ pkq i are defined by From Equation (5), we obtain By iterating Equation (12) and using Equation (13), we have where Φ pkq i and Γ pkq i 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 pk´1, ks. It includes the single rate uniform sampling system as the special case, i.e., N k " 1 and α pkq 1 " 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. (9) under Assumptions 1 and 2, the state second-order moment matrix

Lemma 1. For system Equation
The initial value is q x p0q " µ 0 µ T 0`P 0 .

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.

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 pk´1, ks(N k ‰ 0) are computed bŷ 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 bŷ where the innovation sequence εpk´α The smoothing gain matrices K x pk´1|k´α pkq i q for the state and K w pk´1|k´α pkq i q for the white noise are computed by The smoothing error covariance matrices P x pk´1|k´α pkq i q and P w pk´1|k´α pkq i q, and the cross-covariance matrix P xw pk´1|k´α pkq i q for the state and white noise are computed by The filtering (r " 0) and prediction (r " 1) error covariance matrices P x pk´α pkq i |k´α pkq i´r q for the state are computed by The initial values are P x pk´1|k´α pkq 0 q " P x pk´1|k´1q, P w pk´1|k´α pkq 0 q " Q w , P xw pk´1|kά Proof. See Appendix A.

Remark 2.
In Theorem 2, if N k " 1 and α pkq 1 " 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 pk´1, ks, 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 OpN k p 3 q in the interval pk´1, ks. Compared with the augmented estimator [9] with the computational order of magnitude OpN 3 k q 3 q, our algorithm can obviously reduce the computational cost with the increase of the number N k of measurement samples in the interval pk´1, ks 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.

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 valuesxp0|0q " µ 0 , P x p0|0q " P 0 , P w p0|0q " Q w , P xw p0|0q " 0 and q x p0q " µ 0 µ T 0`P 0 .
Step 2. Construct the model within a state update period pk´1, ks 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 estimateŝ xpk´α pkq i |k´α pkq i q and the covariance matrices P x pk´α pkq i |k´α pkq i q 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.,xpk|kq " Φxpk´1|k´1q, P x pk|kq " ΦP x pk´1|k´1qΦ T`Γ Q w Γ T .

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: where the superscript l denotes the lth sensor, L is the number of sensors, y plq pk plq j q is the measurement at the sampling time k plq j for the lth sensor, k plq j is the sampling time of the jth measurement, v plq pk plq j q is the measurement noise with zero mean and covarianceQ v plq , and H plq is the measurement matrix. The variable ξ plq pk plq j q is a Bernoulli distributed stochastic variable that takes values on 1 and 0 with the probability Prob ! ξ plq pk plq j q " 1 ) " γ plq , γ plq P r0, 1s. We assume that ξ plq pk plq j q is independent of wpkq and v plq pk plq j q and xp0q. Moreover ξ plq pk plq j q satisfies that Etξ plq pk plq j u " γ plq , Etpξ plq pk plq j´γ plq q 2 u " γ plq p1´γ plq q, Etξ plq pk l j q´γ plq u " 0, Etξ plq pk plq j qp1´ξ plq pk plq j qqu " 0, Etξ plq pk plq i qξ psq pm psq j qu " γ plq γ psq pl ‰ s or k plq i ‰ m psq j q We will adopt two types of methods to deal with the estimation problem of multi-sensor case: optimal fusion estimator and suboptimal fusion estimator.

Optimal Fusion Estimator
In fact, we can reorder the measurements y plq pk plq j q of all sensors in each state update period pk´1, ks according to the order of sampling time k plq j . 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.

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] where P CI pk|kq is the upper bound of variance of the CI fusion estimator. The optimal weighted coefficients ω plq pkq satisfy 0 ď ω plq pkq ď 1 and ř L l"1 ω plq pkq " 1, and can be obtained by solving the following optimization problem: 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].

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: where xptq " x 2 ptq ı 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 piq pk piq j q, i " 1, 2, 3 are independent Gaussian noises with zero means and variances Q v piq and uncorrelated with the process noise wpkq with zero-mean and variance Q w . We set Q w " 2, Q v p1q " 2, Q v p2q " 1, Q v p3q " 3, H p1q " " 0 0 1 1 ı , γ p1q " 0.7, γ p2q " 0.9, γ p3q " 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: Φ " We set the initials as xp0q " 0 and P 0 " 0.1I 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]. , 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.   Then Monte Carlo simulation is carried out to compare the algorithms in our paper and [9]. The mean square errors (MSEs): x k x k k N = = −  , 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 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.
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 γ p1q " 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.
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.

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

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.