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
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
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.
where
is the system state at the moment
, the Italic
is the state sampling period, and
is the measurement of the sensor at the moment
.
is the state measured by
.
and
are the system noise and measurement noise.
,
and
are constant matrices with suitable dimensions.
is the sampling time of the jth measurement. The variable
is a Bernoulli distributed stochastic variable that takes values on 1 and 0 with the probability
,
.
means that the measurement is received, while
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:
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. and are uncorrelated white noises with zero means and covariance matrices and .
Assumption 2. The initial state is uncorrelated with , and , 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
, for example, only one measurement sample at the time instants
T, 2
T, 3
T, 4
T, 9
T, 10
T and within the interval
with the sampling time instants
,
,
,
,
,
,
, two samples in
with the sampling time instants
,
, and three samples within
with the sampling time instants
,
,
. Case 2 is that there are no samples within the state update period
, for example, not any measurement sample within
.
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 is the number of measurement samples within a state update interval , and is the total number of measurement samples before the moment . The i () denotes the ith sample within the interval . where denotes the minimal integer not less than . Let , and then we have . It is easily known that . To simplify the expression, the sampling period will be omitted under no confusion in the subsequent text.
The state at the measurement sampling point within the interval
is given by
Accordingly, the measurement equation is expressed as
Especially, when the measurement is sampled at the state update point, we have . Then the state at the measurement sampling point can be reduced as . The measurement equation can be reduced as .
Assuming
and
are two adjacent sampling points in
. 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 can be set up as follows:where the coefficient matrices and are defined bywith , and . Proof. From Equation (8), we have
From Equation (5), we obtain
By iterating Equation (12) and using Equation (13), we have
where
and
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 . It includes the single rate uniform sampling system as the special case, i.e., and . 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 is computed byThe initial value is .
Proof. Substituting Equation (9) into , 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 () are computed bywhere it is a filter for or a predictor for . The fixed-point smoothers for the state and white noise are respectively computed bywhere the innovation sequence and its covariance are computed by The smoothing gain matrices for the state and for the white noise are computed by The smoothing error covariance matrices and , and the cross-covariance matrix for the state and white noise are computed by The filtering () and prediction () error covariance matrices for the state are computed by The initial values are , , , and .
It is clear that we have the estimator at the state update point as when or with , and when .
Remark 2. In Theorem 2, if and , we have the filtering algorithm at the state update points for single-rate systems. If , i.e., no samples in the interval , the predictor is used to estimate the state at the state update point based on the estimator at the state update point .
Remark 3. It can be easily known that the proposed non-augmented recursive estimator has the computational order of magnitude in the interval . Compared with the augmented estimator [9] with the computational order of magnitude , our algorithm can obviously reduce the computational cost with the increase of the number of measurement samples in the interval 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. , set the initial values , , , and .
Step 2. Construct the model within a state update period by Theorem 1.
Step 3. Compute the state estimators at the state update points and at the measurement sampling points:
- (a)
If (i.e., there are samples within a state update period), obtain the state estimates and the covariance matrices by Theorem 2.
- (b)
If (i.e., no sample within a state update period), obtain the estimate by prediction method in Remark 2, i.e., , .
Step 4. , set , , , and .
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:
where the superscript
denotes the lth sensor,
is the number of sensors,
is the measurement at the sampling time
for the lth sensor,
is the sampling time of the
jth measurement,
is the measurement noise with zero mean and covariance
, and
is the measurement matrix. The variable
is a Bernoulli distributed stochastic variable that takes values on 1 and 0 with the probability
,
. We assume that
is independent of
and
and
. Moreover
satisfies that
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
of all sensors in each state update period
according to the order of sampling time
. 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]:
where
is the upper bound of variance of the CI fusion estimator. The optimal weighted coefficients
satisfy
and
, 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]. 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:
where
,
and
are the positions of mass
and mass
, respectively;
and
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.
,
are independent Gaussian noises with zero means and variances
and uncorrelated with the process noise
with zero-mean and variance
. We set
,
,
,
,
,
,
,
,
,
,
,
,
,
,
, and the sampling period
, we obtain the following system parameter matrices by discretization:
We set the initials as and . 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):
,
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
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
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.