Distributed fusion estimation for multisensor systems with non-Gaussian but heavy-tailed noises

ISA Transactions 101 (2020) 160–169

Contents lists available at ScienceDirect

ISA Transactions journal homepage: www.elsevier.com/locate/isatrans

Research article

Distributed fusion estimation for multisensor systems with non-Gaussian but heavy-tailed noises ∗

Liping Yan a,b , , Chenying Di a , Q.M. Jonathan Wu b , Yuanqing Xia a , Shida Liu c a

Key Laboratory of Intelligent Control and Decision of Complex Systems, School of Automation, Beijing Institute of Technology, Beijing 100081, China Department of Electrical and Computer Engineering, University of Windsor, Windsor N9B3P4, Canada c School of Electrical and Control Engineering, North China University of Technology, Beijing 10093, China b

article

info

Article history: Received 22 June 2019 Received in revised form 2 February 2020 Accepted 2 February 2020 Available online 13 February 2020 Keywords: State estimation Information filter Distributed fusion Non-Gaussian disturbance Heavy-tailed noise Multivariate t distribution

a b s t r a c t Student’s t distribution is a useful tool that can model heavy-tailed noises appearing in many practical systems. Although t distribution based filter has been derived, the information filter form is not presented and the data fusion algorithms for dynamic systems disturbed by heavy-tailed noises are rarely concerned. In this paper, based on multivariate t distribution and variational Bayesian estimation, the information filter, the centralized batch fusion, the distributed fusion, and the suboptimal distributed fusion algorithms are derived, respectively. The centralized fusion is given in two forms, namely, from t distribution based filter and the proposed t distribution based information filter, respectively. The distributed fusion is deduced by the use of the newly derived information filter, and it has been demonstrated to be equivalent to the centralized batch fusion. The suboptimal distributed fusion is obtained by a parameter approximation from the derived distributed fusion to decrease the computation complexity. The presented algorithms are shown to be the generalization of the classical Kalman filter based traditional algorithms. Theoretical analysis and exhaustive experimental analysis by a target tracking example show that the proposed algorithms are feasible and effective. © 2020 ISA. Published by Elsevier Ltd. All rights reserved.

1. Introduction Information fusion is an effective technique, which can cope with detection, correlation, estimation by combining useful information of multiple sensors to obtain high-precision and consistent estimated results, and provide effective information for decision making [1]. The key to data fusion is to make the best use of effective data from multiple data sources to estimate the concerned quantities [2]. It was first studied in military applications and developed in some high-technology fields, such as aeronautics and astronautics, robots and manufacturing industry, target tracking and localization, with its higher estimation accuracy and good reliability [1,2]. It is well known that for the linear dynamic systems driven by Gaussian noises, Kalman filter can obtain the optimal estimates in the sense of least mean square error (LMSE) [3,4]. In the framework of Kalman filter, many fusion algorithms have been generated [2,5–8]. Basically, these algorithms can be divided into two categories, such as the centralized fusion and the distributed fusion [9]. If all observation information of multiple local agents ∗ Corresponding author at: Key Laboratory of Intelligent Control and Decision of Complex Systems, School of Automation, Beijing Institute of Technology, Beijing 100081, China. E-mail address: [email protected] (L. Yan). https://doi.org/10.1016/j.isatra.2020.02.004 0019-0578/© 2020 ISA. Published by Elsevier Ltd. All rights reserved.

can be obtained at the same time, the centralized fusion estimation is usually considered as the optimal in the sense of LMSE. Whereas, considering the factors in real applications such as restricted communicating bandwidth and transmission delay, to improve the system robustness in a hostile environment, it is better for each local agent to generate the local estimates by using its local measurements and transmit it to the process center [10,11]. Then, the fusion center will make use of every local estimation to obtain a globally estimation with higher accuracy, this is called distributed fusion. Briefly, the centralized fusion aims at achieving the globally optimal state estimates by dealing with measurements directly, while the key task for the distributed fusion is processing the measurements with higher efficiency and fusing the local estimates for more accurate or reliable state estimation [12]. In several practical applications, there exist unavoidable outliers produced by unreliable source or in the influence of unknown disturbances, the process and measurement noises are both non-Gaussian. For the state estimation of non-Gaussian systems, there are many adaptive algorithms based on modified Kalman filter [13–15]. A more general framework for filter of non-Gaussian systems is based on the mixture of Gaussians [15]. In [14], the statistical fusion approaches based on mixture of Gaussian noises were divided into four categories and a costeffective fusion algorithm has been proposed under Bayesian

L. Yan, C. Di, Q.M.J. Wu et al. / ISA Transactions 101 (2020) 160–169

161

filter formulas for the system with heavy-tailed noise or formulated by Student’s t distribution are not derived. • The distributed fusion for dynamic systems with heavytailed noises is not derived.

Fig. 1. Gaussian distribution and Student’s t distribution with different dof.

In this paper, the information filter for linear time-variant dynamic system with heavy-tailed noises is proposed. Based on the information filter, the distributed fusion algorithm to deal with the estimation problem for multisensor linear time-variant systems is derived, and proven to be equivalent to the derived centralized fusion in the sense of LMSE. We also give a suboptimal distributed fusion to decrease computation complexity. The rest of the paper is organized as follows. In Section 2, the problem is formulated. Section 3 derives the information filter for the linear time-variant dynamic system with heavytailed noises. Section 4 presents the centralized fusion and the distributed fusion algorithms in sequence. Section 5 shows the numerical simulation and Section 6 draws the conclusion. 2. Problem formulation

framework. However, the methods that approximating the nonGaussian noise with Gaussian mixture model might be limited in the cases where outliers exist [13,16]. Instead of Gaussian mixture, an effective way to model nonGaussian noise is using Student’s t distribution [17–19]. As is demonstrated in Fig. 1, compared with the Gaussian distribution, t distribution takes extreme value with higher possibility, which resulting in a heavy-tailed shape and is therefore more robust to outliers, where ν denotes the degrees of freedom (dof) of the t distribution. Multivariate student t distribution based filtering or smoothing for system with heavy-tailed noises are of interests to many researchers in recent years because of its many advantages, such as, it can accommodate outliers, can be obtained with proper computational complexity, and is free of parameter tuning. [18, 20–25]. Literatures [26] and [27] derived the filter and the smoother for linear systems disturbed by heavy-tailed noise. For filter of nonlinear systems, Piché et al. proposed the recursive robust filter and smoother for nonlinear systems by the use of the Student’s t distribution, which is able to solve the problem caused by outliers [18]. An improved Gaussian approximate filtering method for the nonlinear system with heavy-tailed noise is derived by use of the variational Bayesian algorithm [20]. The robust t distribution based nonlinear filtering (RSTNF) and smoothing (RSTNS) algorithms for dynamic systems disturbed by heavy-tailed noise by using the unscented transform (UT) has been derived in [28]. Similar to the RSTNF, by the use of cubature filter, a robust filtering approach for systems with noises modeled by t distribution is presented in [21]. [24] analyzed the stability of heavy-tailed dynamics. [25] studied the modeling and stochastic linearization of heavy-tailed systems. Recently, an unscented particle filter based on multivariate t distribution is studied [22]. The above algorithms only consider the state estimation for filtering or smoothing of heavy-tailed systems with single sensor. For fusion estimation of multisensor systems with heavy-tailed noises, the results are few. References [29] and [16] studied the centralized fusion of linear time-invariant and nonlinear systems with heavy-tailed noises, respectively. Through the above analysis, we find that there are many open problems for filtering and fusion estimation of dynamic systems disturbed by heavy-tailed noises:

One target observed by N sensors can be formulated as, [2,30, 31] xl+1 = Fl xl + wl ,

l = 0, 1, . . . ,

(1)

yi,l = Hi,l xl + vi,l , i = 1, 2, . . . , N ,

(2)

where i = 1, 2, . . . , N denotes the ith sensor. xl ∈ R is the system state at the lth time instant, Fl ∈ Rq×q is the state transition matrix. yi,l ∈ Rmi is the observation of sensor i at time l. Hi,l is the observation matrix of the ith sensor. The system noise wl and the measurement noise vi,l are heavy-tailed noises, which can be modeled by the multivariate t distribution as, q

p(wl ) = St(wl ; 0, Ql , νw ),

(3)

p(vi,l ) = St(vi,l ; 0, Ri,l , νi ),

(4)

where St(·; x¯ , P , ν ) is the multivariate t distribution whose mean is x¯ , the scale matrix is P and degrees of freedom (dof) is ν . For a random vector x with the probability density function of ν St(·; x¯ , P , ν ), ν− P is the covariance matrix [32]. 2 Similarly, we assume that the system state initial value x0 also to be heavy-tailed that meets the multivariate t distribution whose mean is xˆ 0|0 , the scale matrix is P0|0 and the dof is ν0 , i.e., p(x0 ) = St(x0 ; xˆ 0|0 , P0|0 , ν0 ).

(5)

It is assumed that x0 , wl and vi,l are uncorrelated. The goal for our work is to obtain the optimal estimation of xl by utilizing multisensor observations Yl = {yi,t , t = 1, 2, . . . , l; i = 1, 2, . . . , N }. 3. Filter and information filter for single sensor systems with heavy-tailed noises For system (1)–(2), in case of N = 1, it is the single sensor linear time-variant dynamic system. For simplicity, we denote yl = y1,l , Hl = H1,l , vl = v1,l , and Rl = R1,l . Then, system (1)–(2) can be rewritten as xl+1 = Fl xl + wl ,

(6)

yl = Hl xl + vl ,

(7)

where,

• The information filter for dynamic systems with heavytailed noise is not derived. Although in [16], the fusion estimation is obtained by the use of cubature information filter (CIF), the real information

⎧ ⎨p(wl ) = St(wl ; 0, Ql , νw ), p(v ) = St(vl ; 0, Rl , ν1 ), ⎩ l p(x0 ) = St(x0 ; xˆ 0|0 , P0|0 , ν0 ).

(8)

162

L. Yan, C. Di, Q.M.J. Wu et al. / ISA Transactions 101 (2020) 160–169

In this section, the information filter for single sensor system (6)–(8) is derived. Before that, we will introduce three lemmas first. Lemma 1.

Random vector x obeys the multivariate t distribution

St(x; x¯ , P , ν ) if its probability density function (pdf) is:

p(x) =

Γ ( ν+2 2 )

1

Γ ( ν2 ) (νπ ) 2d

√

(

1

det(P)

∆2 1+ ν

)− ν+2 2

,

(9)

where ∆2 = (x − x¯ )T P −1 (x − x¯ ). x¯ is the mean of x, P is the scale matrix, and ν is the degrees of freedom (dof). It has the following properties:

• The covariance of x is ν−ν 2 P. • The distribution of x tends to Gaussian when ν increases to infinity.

• Let z = Ax + B, then p(z) = St(z ; Ax¯ + B, APAT , ν ), where A and B are matrix and vector of proper dimensions.

• If x1 ∈ Rd1 and x2 ∈ Rd2 are jointly t-distributed with pdf of

([ ] [ ] [ x¯ P x1 ; 1 , 11 p(x1 , x2 ) = St ¯ x2

x2

P21

]

P12 ,ν P22

)

,

(10)

where Pii ∈ Rdi ×di are the scale matrices of xi , Pij ∈ Rdi ×dj are

(6)–(8) is given by [28]

⎧ ⎪ ⎪ ⎪xˆ l|l−1 = Fl−1 xˆ l−1|l−1 ⎪ ⎪ ⎪ ⎪ ⎪ Pl|l−1 = Fl−1 Pl−1|l−1 FlT−1 + Q¯ l−1 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ yˆ l|l−1 = Hl xˆ l|l−1 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪xˆ l|l = xˆ l|l−1 + Kl (yl − yˆ l|l−1 ) ⎪ ⎪ ⎪ ⎪ y˜ y˜ ⎪ Pl|l = al (Pl|l−1 − Kl Pl|l−1 KlT ) ⎪ ⎪ ⎪ ⎪ ⎪ x˜ y˜ y˜ y˜ ⎪ ⎪ Kl = Pl|l−1 (Pl|l−1 )−1 ⎪ ⎪ ⎨ y˜ y˜ Pl|l−1 = Hl Pl|l−1 HlT + R¯ l ⎪ ⎪ ⎪ ⎪P x˜ y˜ = Pl|l−1 H T ⎪ ⎪ l l|l−1 ⎪ ⎪ ⎪ (ν0 −2)(ν0 +∆2l ) ⎪ ⎪ al = ν (ν +m −2) ⎪ ⎪ ⎪ √0 0 1 ⎪ ⎪ ⎪ ⎪∆l = (yl − yˆ l|l−1 )T (Ply˜|ly˜−1 )−1 (yl − yˆ l|l−1 ) ⎪ ⎪ ⎪ ⎪ ⎪ νw (ν0 −2) ⎪ ¯ ⎪ ⎪Ql = (νw −2)ν0 Ql ⎪ ⎪ ⎪ ⎪ ⎪ R¯ l = bRl ⎪ ⎪ ⎪ ⎪ ⎩ (ν −2)ν b = ν 0(ν −2)1 , 0 1

where xˆ l|l−1 and xˆ l|l denote the state prediction and the state estimate, Pl|l−1 and Pl|l are the scale matrices of the prediction error and the estimation error. m1 is the dimension of measurement yl . yˆ l|l−1 y˜ y˜ is the measurement prediction, and Pl|l−1 is the scale matrix of its prediction error. Kl is the innovation matrix.

the joint scale matrices of xi and xj , i = 1, 2, then the marginal

Proof. From Lemma 1 [28],

pdfs of x1 and x2 are

xˆ l|l−1 =

{

(17)

∫

p(x1 ) = St(x1 ; x¯ 1 , P11 , ν ),

= Fl−1 xˆ l−1|l−1 ,

(11)

p(x2 ) = St(x2 ; x¯ 2 , P22 , ν ).

xl p(xl |Yl−1 )dxl (18)

Denote y˜ l|l−1 = yl − yˆ l|l−1 , x˜ l|l−1 = xl − xˆ l|l−1 , then

Lemma 3 (Filter for Linear Dynamic System with Heavy-tailed

ν0 − 2 E [˜xl|l−1 x˜ Tl|l−1 |Yl−1 ] ν0 ∫ ν0 − 2 x˜ l|l−1 x˜ Tl|l−1 p(xl |Yl−1 )dxl = ν0 ∫ ν0 − 2 ν0 − 2 = xl xTl p(xl |Yl−1 )dxl − xˆ l|l−1 xˆ Tl|l−1 ν0 ∫ ν0 ∫ ν0 − 2 = xl xTl { St(xl ; Fl−1 xl−1 , Ql−1 , νw ) ν0 × St(xl−1 ; xˆ l−1|l−1 , Pl−1|l−1 , ν0 )dxl−1 }dxl ν0 − 2 − xˆ l|l−1 xˆ Tl|l−1 (19) ν0 ν0 ν0 − 2 = Fl−1 (xˆ l−1|l−1 xˆ Tl−1|l−1 + Pl−1|l−1 )FlT−1 ν0 ν0 − 2 ν0 − 2 νw (ν0 − 2) − xˆ l|l−1 xˆ Tl|l−1 + Ql−1 ν0 (νw − 2)ν0 νw (ν0 − 2) = Fl−1 Pl−1|l−1 FlT−1 + Ql−1 , (νw − 2)ν0 ∫ = E [yl |Yl−1 ] = yl P(yl |Yl−1 )dyl ∫ ∫ = { yl St(yl ; Hl xl , Rl , ν1 )dyl } × St(xl ; xˆ l|l−1 , Pl|l−1 , ν0 )dxl ∫ = Hl xl St(xl ; xˆ l|l−1 , Pl|l−1 , ν0 )dxl (20)

Noises).

= Hl xˆ l|l−1 ,

The conditional pdf p(x1 |x2 ) is given by

Pl|l−1 =

p(x1 |x2 ) = St(x1 ; x¯ 1|2 , P1|2 , ν1|2 ),

(12)

where

ν1|2 = ν + d2 ,

(13)

x¯ 1|2 = x¯ 1 + P12 P22 (x2 − x¯ 2 ), −1

P1|2

(14)

ν + ∆22 −1 T = (P11 − P12 P22 P12 ), ν + d2

(15)

−1 and where ∆22 = (x2 − x¯ 2 )T P22 (x2 − x¯ 2 ).

Lemma 2 (Matrix Inversion Lemma [3]).

Suppose M1 ∈ Rn×n ,

M2 ∈ Rn×p , M3 ∈ Rp×n , M4 ∈ Rp×p , where M1 , M1 + M2 M4−1 M3 , M4 + M3 M1−1 M2 are full rank, then, yˆ l|l−1 (M1 + M2 M4−1 M3 )−1 = M1−1 − M1−1 M2 (M4 + M3 M1−1 M2 )−1 M3 M1−1 . (16)

The multivariate t distribution based filter for system

L. Yan, C. Di, Q.M.J. Wu et al. / ISA Transactions 101 (2020) 160–169 y˜ y˜

ν0 − 2 E [˜yl|l−1 y˜ Tl|l−1 |Yl−1 ] ν0 ∫ ν0 − 2 ν0 − 2 = yˆ l|l−1 yˆ Tl|l−1 yl yT p(yl |Yl−1 )dyl − ν0 ∫ ∫ l ν0 ν0 − 2 = { yl yTl St(yl ; Hl xl , Rl , ν1 )dyl } ν0 ν0 − 2 × St(xl ; xˆ l|l−1 , Pl|l−1 , ν0 )dxl − yˆ l|l−1 yˆ Tl|l−1 ν 0 ∫ ν0 ν0 − 2 = Hl (xˆ l|l−1 xˆ Tl|l−1 + Pl|l−1 )HlT dxl ν0 ν0 − 2 ν0 − 2 ν1 (ν0 − 2) − yˆ l|l−1 yˆ Tl|l−1 + Rl ν0 (ν1 − 2)ν0 (ν0 − 2)ν1 = Hl Pl|l−1 HlT + Rl , ν0 (ν1 − 2)

Theorem 1 (Information Filter for Linear Dynamic Systems with Heavy-tailed Noises). The multivariate t distribution based information filter for system (6)–(8) is computed by

Pl|l−1 =

x˜ y˜

(21)

ν0 − 2 E [˜xl|l−1 y˜ Tl|l−1 |Yl−1 ] ν0 ∫ ∫ ν0 − 2 ν0 − 2 = xl yTl p(xl , yl |Yl−1 )dxl dyl − xˆ l|l−1 yˆ Tl|l−1 ν0 ∫ ν0 ∫ ν0 − 2 = xl { yTl p(yl |xl )dyl }p(xl |Yl−1 )dxl ν0 ν0 − 2 − xˆ l|l−1 yˆ Tl|l−1 ν0 ∫ ∫ ν0 − 2 xl { yTl St(yl ; Hl xl , Rl , ν1 )dyl } (22) = ν0 ν0 − 2 × St(xl ; xˆ l|l−1 , Pl|l−1 , ν0 )dxl − xˆ l|l−1 yˆ Tl|l−1 ν0 ν0 − 2 ν0 ν0 − 2 = xˆ l|l−1 yˆ Tl|l−1 (xˆ l|l−1 yˆ Tl|l−1 + Pl|l−1 )HlT − ν0 ν0 − 2 ν0 = Pl|l−1 HlT .

Pl|l−1 =

⎧ 1 ⎪ ⎪ ξˆl|l = a1 (ξˆl|l−1 + HlT R¯ − ⎪ l yl ) ⎪ l ⎪ ⎪ ⎪ 1 −1 −T ˆ ⎪ ⎪ξˆl|l−1 = [I − Al−1 (Al−1 + Q¯ l− −1 ) ]Fl−1 ξl−1|l−1 ⎪ ⎪ ⎪ ⎪ 1 ⎪ ⎪ Λl = a1 (Λl|l−1 + HlT R¯ − ⎪ l Hl ) l ⎪ ⎪ ⎪ ⎪ 1 −1 ⎪ Λl|l−1 = [I − Al−1 (Al−1 + Q¯ l− ⎪ −1 ) ]Al−1 ⎪ ⎪ ⎪ 2 (ν0 −2)(ν0 +∆ ) ⎪ ⎪ a = ν (ν +m −2)l ⎪ ⎪ 0 0 1 ⎨ l ¯Rl = ν1 (ν0 −2) Rl = bRl (ν1 −2)ν0 ⎪ ⎪ ⎪ ⎪ T −1 ⎪ ⎪ Al−1 = Fl− ⎪ −1 Λl−1 Fl−1 ⎪ ⎪ ⎪ ν (ν −2) ⎪ ⎪ Q¯ l−1 = (νw −02)ν Ql−1 ⎪ w ⎪ 0 ⎪ √ ⎪ ⎪ y˜ y˜ ⎪ T −1 ⎪ ⎪∆l = (yl − yˆ l|l−1 ) (Pl|l−1 ) (yl − yˆ l|l−1 ) ⎪ ⎪ ⎪ ⎪ y˜ y˜ 1 T ⎪ ¯ Pl|l−1 = Hl Λ− ⎪ l|l−1 Hl + Rl ⎪ ⎪ ⎪ ⎩ 1 ˆ yˆ l|l−1 = Hl Λ− l|l−1 ξl|l−1 ,

p(xl , yl |Yl−1 ) p(yl |Yl−1 )

= St(xl ; xˆ ′l|l , Pl′|l , ν0′ ),

(23)

{ Λl = Pl−|l 1 , Λl|l−1 = Pl−|l−1 1 ,

ξˆl|l = Pl−|l 1 xˆ l|l = Λl xˆ l|l , ξˆl|l−1 = Pl−|l−1 1 xˆ l|l−1 = Λl|l−1 xˆ l|l−1 .

Proof. From Lemmas 3 and 2, we have 1 Pl− |l =

=

1 al 1 al

y˜ y˜

(Pl|l−1 − Kl Pl|l−1 KlT )−1 1 T ¯ −1 (Pl− |l−1 + Hl Rl Hl ).

(33)

−1 1 Let Λl = Pl− |l and Λl|l−1 = Pl|l−1 , we have

1 al

1 (Λl|l−1 + HlT R¯ − l Hl ).

(34)

Let −1 −T −1 T −1 Al−1 = Fl− −1 Pl−1|l−1 Fl−1 = Fl−1 Λl−1 Fl−1 ,

(35)

Λl|l−1 = Pl−|l−1 1 = [Fl−1 Pl−1|l−1 FlT−1 + Q¯ l−1 ]−1

ν0′ = ν0 + m1 ,

(24)

xˆ ′l|l = xˆ l|l−1 + Kl [yl − yˆ l|l−1 ],

(25)

= Al−1 − Al−1 (I + Q¯ l−1 Al−1 )−1 Q¯ l−1 Al−1 = [I − Al−1 (Al−1 + Q¯ −1 )−1 ]Al−1 .

(26)

Let ξˆl|l−1 = Λl|l−1 xˆ l|l−1 , then from Eq. (17), we have

ν0 + ∆2l y˜ y˜ Pl′|l = (Pl|l−1 − Kl Pl|l−1 KlT ), ν0 + m1 x˜ y˜ y˜ y˜ Kl = Pl|l−1 (Pl|l−1 )−1 , √ y˜ y˜ ∆l = (yl − yˆ l|l−1 )T (Pl|l−1 )−1 (yl − yˆ l|l−1 ).

(27)

ξˆl|l−1 = Λl|l−1 Fl−1 xˆ l−1|l−1 .

(37)

Substitute (36) and (35) to (37), we have (28)

ξˆl|l−1 = Λl|l−1 Fl−1 xˆ l−1|l−1

matching approach is used, and the final state estimation and the

1 −1 ˆ l−1|l−1 = [I − Al−1 (Al−1 + Q¯ l− −1 ) ]Al−1 Fl−1 x − 1 = [I − Al−1 (Al−1 + Q¯ )−1 ]F −T Λl−1 F −1 Fl−1 xˆ l−1|l−1

estimation error scale matrix are given by

= [I − Al−1 (Al−1 + ¯

(ν0 − 2)ν0

Pl′|l =

(ν0 + ∆2l )(ν0 − 2)

ν0 (ν0 + m1 − 2)

l−1 1 −1 Ql− −1 )

]

l−1 l−1 T Fl− −1 l−1|l−1

ξˆ

(38)

.

Similar to Kalman filter, Pl|l and Kl in Lemma 3 can be rewritten

xˆ l|l = xˆ ′l|l = xˆ l|l−1 + Kl [yl − yˆ l|l−1 ], ′

(36)

l−1

From [23,28], to preserve the heavy-tailed property, the moment

(ν0 − 2)ν0′

(32)

from (17), and by the use of Lemma 2, we have

where xˆ ′l|l , Pl′|l and ν0′ are given by

Pl|l =

(31)

where m1 is the dimension of yl , and

Λl =

From Lemma 1, p(xl |Yl ) =

163

(29)

as Pl|l = al [I − Kl Hl ]Pl|l−1 ,

y˜ y˜

(Pl|l−1 − Kl Pl|l−1 KlT ). (30)

Kl =

1 al

1 Pl|l HlT R¯ − l .

(39) (40)

164

L. Yan, C. Di, Q.M.J. Wu et al. / ISA Transactions 101 (2020) 160–169

Corollary 1. The centralized fusion of linear dynamic system (1)– (5) by the use of the information filter presented in Theorem 1 can be generated by

Therefore, by the use of Eqs. (34), (39) and (40), we have

ξˆl|l = Λl xˆ l|l = Λl [ˆxl|l−1 + Kl (yl − yˆ l|l−1 )] = Λl [ˆxl|l−1 + = =

1 al 1 al

1 al

Λl

−1

1 HlT R− l (yl

¯

− yˆ l|l−1 )]

1 ˆ (Λl|l−1 + HlT R¯ − l Hl )xl|l−1 +

1 al

(41)

1 ˆ HlT R¯ − l (yl − Hl xl|l−1 )

1 [ξˆl|l−1 + HlT R¯ − l yl ].

The rest of Eq. (31) can be obtained from Eq. (17) directly.

⎧ a,T ¯ a,−1 a 1 ˆ ˆ ⎪ ⎪ ⎪ξc ,l|l = ac ,l (ξc ,l|l−1 + Hl Rl yl ) ⎪ ⎪ −1 −1 −T ⎪ ⎪ ⎨ξˆc ,l|l−1 = [I − Ac ,l−1 (Ac ,l−1 + Q¯ l−1 ) ]Fl−1 ξˆc ,l−1|l−1 a,T ¯ a,−1 a 1 Λc ,l = a (Λc ,l|l−1 + Hl Rl Hl ) ⎪ c ,l ⎪ ⎪ ⎪Λ ¯ −1 −1 ]Ac ,l−1 ⎪ c ,l|l−1 = [I − Ac ,l−1 (Ac ,l−1 + Ql−1 ) ⎪ ⎪ ⎩ −T −1 Ac ,l−1 = Fl−1 Λc ,l−1 Fl−1

,

(46)

where

4. The information fusion algorithms In this section, for simplicity, we consider the fusion algorithms under the assumption that νi = ν1 for all i = 2, 3, . . . , N. When the measurement noises have different dof of t distributions, the fusion algorithms will be introduced in Remark 3.

{ Λc ,l = Pc−,l1|l , Λc ,l|l−1 = Pc−,l1|l−1 ,

ξˆc ,l|l = Λc ,l xˆ c ,l|l , ξˆc ,l|l−1 = Λc ,l|l−1 xˆ c ,l|l−1 .

(47)

and where Hla , R¯ al , ac ,l , yal are computed by (42)–(44). 4.2. The distributed fusion algorithms

4.1. The centralized batch fusion Theorem 2. For the linear dynamic system (1)–(5), when νi = ν1 for all i = 2, 3, . . . , N, the state estimation by centralized fusion of sensors 1 to N can be computed by

⎧ ⎪ ⎪ xˆ c ,l|l−1 = Fl−1 xˆ c ,l−1|l−1 , ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ Pc ,l|l−1 = Fl−1 Pc ,l−1|l−1 FlT−1 + Q¯ l−1 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ yˆ c ,l|l−1 = Hla xˆ c ,l|l−1 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ xˆ c ,l|l = xˆ c ,l|l−1 + Kc ,l (yal − yˆ c ,l|l−1 ) ⎪ ⎪ ⎪ ⎪ ⎪ y˜ y˜ ⎪ ⎪ Pc ,l|l = ac ,l (Pc ,l|l−1 − Kc ,l Pc ,l|l−1 KcT,l ) ⎪ ⎪ ⎪ ⎪ ⎪ ⎨Kc ,l = Pcx˜,y˜l|l−1 (Pcy˜,y˜l|l−1 )−1

(42)

y˜ y˜ a ,T ⎪ Pc ,l|l−1 = Hla Pc ,l|l−1 Hl + R¯ al ⎪ ⎪ ⎪ ⎪ ⎪ x˜ y˜ a,T ⎪ ⎪ Pc ,l|l−1 = Pc ,l|l−1 Hl ⎪ ⎪ ⎪ 2 ⎪ (ν0 −2)(ν0 +∆c ,l ) ⎪ ⎪ ⎪ ⎪ac ,l = √ν0 (ν0 +m−2) ⎪ ⎪ ⎪ y˜ y˜ ⎪∆ = (ya − yˆ T −1 (ya − y ⎪ ˆ c ,l|l−1 ) c ,l c ,l|l−1 ) (Pc ,l|l−1 ) ⎪ l l ⎪ ⎪ ⎪ ⎪ ν (ν −2) ⎪ ⎪ Q¯ = w 0 Q ⎪ ⎪ l (νw −2)ν0 l ⎪ ⎪ ⎩¯a ν (ν −2) R = bRa = 1 0 Ra , l

l

(ν1 −2)ν0

∑N

l

where the subscript c denotes the centralized fusion, m = and

∑N

i=1

⎤ v1,l ⎢ v2,l ⎥ ⎥ a ⎢ ⎥ a ⎢ ⎥ yal = ⎢ ⎣ .. ⎦ , Hl = ⎣ .. ⎦ , vl = ⎣ .. ⎦ , . . . y N ,l HN ,l v N ,l y1,l ⎢ y2,l ⎥

⎡

⎤

H1,l ⎢ H2,l ⎥

⎡

⎤

Theorem 3 (The Optimal Distributed Fusion). For the linear dynamic system (1)–(5), when νi = ν1 for all i = 2, 3, . . . , N, the optimal distributed fusion of sensors 1 to N for state estimation can be computed by ⎧ 1ˆ ⎪ xˆ d,l|l = Λ− ⎪ d,l ξd,l|l ⎪ ⎪ ⎪ − ⎪ ⎪Pd,l|l = Λd,1l ⎪ ⎪ ⎪ N ⎪ ∑ ⎪ ⎪ ⎪ (ai,l Λi,l − Λd,l|l−1 )] Λd,l = a1 [Λd,l|l−1 + 1b ⎪ ⎪ d , l ⎪ ⎪ i=1 ⎪ ⎪ N ⎪ ⎪ ∑ ⎪ ⎪ ⎪ [ai,l ξˆi,l|l − ξˆd,l|l−1 ]) ξˆd,l|l = a1 (ξˆd,l|l−1 + 1b ⎪ d,l ⎪ ⎪ ⎪ i = 1 ⎪ ⎪ ⎪ ⎪ ¯ −1 −1 −T ˆ ˆ ⎪ ⎨ξd,l|l−1 = [I − Ad,l−1 (Ad,l−1 + Ql−1 ) ]Fl−1 ξd,l−1|l−1 1 −1 (48) Λd,l|l−1 = [I − Ad,l−1 (Ad,l−1 + Q¯ l− , −1 ) ]Ad,l−1 ⎪ ⎪ (ν0 −2)(ν0 +∆2d,l ) ⎪ ⎪ ⎪ ad,l = ν (ν +m−2) ⎪ 0 0 ⎪ ⎪ ⎪ T ⎪ Ad,l−1 = Fl− Λ F −1 ⎪ ⎪ ⎪ ]T ] [ [ −1 d,l−1 l−1 ⎪ ⎪ y˜ y˜ 2 −1 ⎪ T T T ⎪ y˜ T1,l y˜ T2,l · · · y˜ TN ,l ⎪∆d,l = y˜ 1,l y˜ 2,l · · · y˜ N ,l (Pd,l|l−1 ) ⎪ ⎪ ⎡ ⎤ ⎪ ⎪ ⎪ ⎪ H1,l Pd,l|l−1 H1T,l + bR1,l · · · H1,l Pd,l|l−1 HNT ,l ⎪ ⎪ ⎢ ⎥ ⎪ ⎪ ⎢ ⎥ . . y˜ y˜ .. ⎪ ⎪ . . Pd,l|l−1 = ⎢ ⎥ ⎪ . . . ⎪ ⎣ ⎦ ⎪ ⎪ ⎩ HN ,l Pd,l|l−1 H1T,l · · · HN ,l Pd,l|l−1 HNT ,l + bRN ,l

mi ,

⎡

p(vla ) = St(vla ; 0, Ral , ν1 ),

(43)

{

Ral = diag{R1,l , R2,l , . . . , RN ,l }.

(44)

Proof. From (43), (2) can be rewritten as yal = Hla xl + vla .

(45)

From the property of multivariate t distribution [32], we have (44). From Lemma 3, we have (42). From Theorems 1 and 2, we have the following corollary.

¯ where m = i=1 mi , b and Ql are the same as Lemma 3, the subscript d denotes the distributed fusion, the initial state meets (5). For i = 1, 2, . . . , N, the local state estimation and the local state prediction are computed by

⎧ −1 ⎪ ˆ ˆ ⎪ ⎪xi,l|l = Λi,l yi,l|l ⎪ ⎪ − 1 ⎪ Pi,l|l = Λi,l ⎪ ⎪ ⎪ ⎪ 1 ˆ T ¯ −1 ⎪ ˆ ⎪ ⎪ξi,l|l = ai,l (ξd,l|l−1 + Hi,l Ri,l yi,l ) ⎪ ⎪ ⎪ 1 ⎪Λi,l = 1 (Λd,l|l−1 + HiT,l R¯ − ⎪ i,l Hi,l ) ai,l ⎪ ⎪ 2 ⎪ (ν0 −2)(ν0 +∆ ) ⎨ ai,l = ν (ν +m −2)i,l i 0 0 ⎪ ¯ i,l = bRi,l ⎪ R ⎪ ⎪ √ ⎪ ⎪ ⎪∆i,l = y˜ T (P y˜ y˜ )−1 y˜ i,l ⎪ i,l i,l|l−1 ⎪ ⎪ ⎪ ⎪ y˜ y˜ 1 T ⎪ ¯ Pi,l|l−1 = Hi,l Λ− ⎪ d,l|l−1 Hi,l + Ri,l ⎪ ⎪ ⎪ − 1 ⎪yˆ i,l|l−1 = Hi,l Λ ˆ ⎪ ⎪ d,l|l−1 ξd,l|l−1 ⎪ ⎩ y˜ i,l = yi,l − yˆ i,l|l−1

(49)

L. Yan, C. Di, Q.M.J. Wu et al. / ISA Transactions 101 (2020) 160–169

Proof. From Theorem 1, by the use of the information filter to sensor i, the local state estimation can be computed by

⎧ −1 ˆ ⎪ ⎪ ⎪xˆ i,l|l = Λi,l ξi,l|l ⎪ ⎪ ⎪ 1 ⎪ Pi,l|l = Λ− ⎪ i,l ⎪ ⎪ ⎪ ⎪ 1 ⎪ ξˆi,l|l = a1 (ξˆi,l|l−1 + HiT,l R¯ − ⎪ i,l yi,l ) ⎪ i,l ⎪ ⎪ ⎪ 1 −1 −T ˆ ⎪ ⎪ξˆi,l|l−1 = [I − Ai,l−1 (Ai,l−1 + Q¯ l− ⎪ −1 ) ]Fl−1 ξi,l−1|l−1 ⎪ ⎪ ⎪ 1 ⎪ ⎪ Λi,l = a1 (Λi,l|l−1 + HiT,l R¯ − ⎪ i,l Hi,l ) i ,l ⎪ ⎪ ⎪ ⎪ 1 −1 ⎪ Λi,l|l−1 = [I − Ai,l−1 (Ai,l−1 + Q¯ l− ⎪ −1 ) ]Ai,l−1 ⎪ ⎪ 2 ) ⎪ ( ν − 2)( ν + ∆ 0 ⎨a = 0 i,l i ,l

yN ,l

=

=

ν0 (ν0 +mi −2)

(50)

⎪ ν (ν −2) ⎪ R¯ i,l = (ν1 −02)ν Ri,l ⎪ ⎪ 1 0 ⎪ ⎪ ⎪ −T −1 ⎪ ⎪ A = F Λ i , l − 1 l−1 i,l−1 Fl−1 ⎪ ⎪ ⎪ ⎪ ν (ν −2) ⎪ Q¯ l−1 = (νw −02)ν Ql−1 ⎪ ⎪ 0 ⎪ √w ⎪ ⎪ y˜ y˜ T ⎪ ⎪∆i,l = y˜ i,l (Pi,l|l−1 )−1 y˜ i,l ⎪ ⎪ ⎪ ⎪ y˜ y˜ ⎪ 1 T ¯ ⎪ Pi,l|l−1 = Hi,l Λ− ⎪ i,l|l−1 Hi,l + Ri,l ⎪ ⎪ ⎪ ⎪ 1 ˆ ⎪ yˆ i,l|l−1 = Hi,l Λ− ⎪ i,l|l−1 ξi,l|l−1 ⎪ ⎪ ⎩ y˜ i,l = yi,l − yˆ i,l|l−1

Λ c ,l =

=

=

1 ac , l 1 ac , l 1 ac , l 1 ac , l 1 ac , l

ac , l =

xˆ i,l|l−1 = Fl−1 xˆ d,l−1|l−1 , Pi,l|l−1 = Fl−1 Pd,l−1|l−1 FlT−1 + Q¯ l−1 ,

{

∆ c ,l =

(53)

ξˆd,l−1|l−1 = Λd,l−1 xˆ d,l−1|l−1 = Λc ,l−1 xˆ c ,l−1|l−1 = ξˆc ,l−1|l−1 .

(54)

Thus T −1 Ad,l−1 = Fl− −1 Λd,l−1 Fl−1 = Ac ,l−1 ,

(55)

1 −1 = [I − Ad,l−1 (Ad,l−1 + Q¯ l− (56) −1 ) ]Ad,l−1 = Λc ,l|l−1 , −1 −1 −T ˆ ¯ ˆ = [I − Ad,l−1 (Ad,l−1 + Ql−1 ) ]Fl−1 ξd,l−1|l−1 = ξc ,l|l−1 . (57)

a,T a,−1

(Λc ,l|l−1 + Hl R¯ l (Λd,l|l−1 +

[Λd,l|l−1 +

N 1∑

b

Hla )

1 HiT,l R¯ − i,l Hi,l )

(59)

i=1

N 1∑

b

(ai,l Λi,l − Λd,l|l−1 )].

i=1

,

(60)

y˜ y˜

(61)

(yal − yˆ c ,l|l−1 )T (Pc ,l|l−1 )−1 (yal − yˆ c ,l|l−1 ). a,T

y˜ y˜

+ R¯ al ⎤T ⎡

⎡

H1,l ⎢ H2,l ⎥

H1,l ⎢ H2,l ⎥

⎤

⎡

R1,l

0

⎢ ⎢ 0 ⎥ ⎥ ⎢ ⎢ = ⎣ . ⎦ Pd,l|l−1 ⎣ . ⎦ + b ⎢ ⎢ . .. .. ⎣ ..

R2,l

0

···

H N ,l

H1,l Pd,l|l−1 H1T,l + bR1,l

⎡ ⎢ =⎣ =

.. .

HN ,l Pd,l|l−1 H1T,l

··· .. . ···

··· ..

.

0

⎤

.. ⎥ . ⎥ ⎥ ⎥ ⎦

RN ,l

H1,l Pd,l|l−1 HNT ,l

.. .

⎤

HN ,l Pd,l|l−1 HNT ,l + bRN ,l

y˜ y˜ Pd,l|l−1 .

⎥ ⎦

(62) Hla xc ,l|l−1 ,

ˆ Substitute (62) and yˆ c ,l|l−1 = xˆ d,l|l−1 = xˆ c ,l|l−1 , Pd,l|l−1 = Pc ,l|l−1 , xˆ i,l|l−1 = xˆ d,l|l−1 = xˆ c ,l|l−1 to (61), we have y˜ y˜ ∆2c ,l = (yal − Hla xˆ c ,l|l−1 )T (Pc ,l|l−1 )−1 (yal − Hla xˆ c ,l|l−1 ) ⎛⎡ ⎤ ⎡ ⎤⎞T y1,l H1,l xˆ d,l|l−1 ⎜⎢ y2,l ⎥ ⎢ H2,l xˆ d,l|l−1 ⎥⎟ y˜ y˜ −1 ⎢ ⎥ ⎢ ⎥⎟ (P =⎜ .. ⎝⎣ .. ⎦ − ⎣ ⎦⎠ d,l|l−1 ) . . yN ,l HN ,l xˆ d,l|l−1 ⎛⎡ ⎤ ⎡ ⎤⎞ y1,l H1,l xˆ d,l|l−1 ⎜⎢ y2,l ⎥ ⎢ H2,l xˆ d,l|l−1 ⎥⎟ ⎢ ⎥ ⎢ ⎥⎟ ×⎜ .. ⎝⎣ .. ⎦ − ⎣ ⎦⎠ . .

HN ,l xˆ d,l|l−1

yN ,l yT1,l

[

= ˜ a,T a,−1 a Hl Rl yl )

¯

···

yTN ,l

˜

]

y˜ y˜

(Pd,l|l−1 )−1 y˜ T1,l

[

···

y˜ TN ,l

]T

= ∆2d,l , (63)

⎡

⎡ ⎤T R1,l H1,l ⎜ ⎢ H 1 ⎜ 1 ⎢ 2,l ⎥ ⎢ 0 ⎜ ⎢ ⎥ ˆ = ξd,l|l−1 + ⎣ . ⎦ ⎢ ⎢ . .. ac , l ⎜ b ⎣ .. ⎝ HN ,l

i=1

Pc ,l|l−1 = Hla Pc ,l|l−1 Hl

From Corollary 1, and by the use of (49), we have

⎛

(58)

N 1∑ (ai,l ξˆi,l|l − ξˆd,l|l−1 )], b

ν0 (ν0 + m − 2)

√

(52)

Λd,l−1 = Pd−,l1−1|l−1 = Pc−,l1−1|l−1 = Λc ,l−1 ,

(ξˆc ,l|l−1 +

1 HiT,l R− i,l yi,l )

i=1

HN ,l

Substitute (52) to (50), we have (49). In the sequel, we will prove this theorem by the use of the deduction method. Suppose xˆ d,l−1|l−1 = xˆ c ,l−1|l−1 and Pd,l−1|l−1 = Pc ,l−1|l−1 , we will show that xˆ d,l|l = xˆ c ,l|l and Pd,l|l = Pc ,l|l . Since xˆ d,l−1|l−1 = xˆ c ,l−1|l−1 and Pd,l−1|l−1 = Pc ,l−1|l−1 , we have

ac , l

b

(ν0 − 2)(ν0 + ∆2c ,l )

(51)

{ Λi,l|l−1 = Pi−,l|1l−1 = Pd−,l1|l−1 = Λd,l|l−1 , 1 ˆ i,l|l−1 = Pd−,l1|l−1 xˆ d,l|l−1 = ξˆd,l|l−1 . yˆ i,l|l−1 = Pi− ,l|l−1 x

ξˆc ,l|l =

[ξˆd,l|l−1 +

N 1∑

From Theorem 2, we have

which imply

1

(ξˆd,l|l−1 +

From Eq. (42),

To improve the accuracy and the robustness of the local state estimations, the fused state estimation is transmitted to the local, i.e., let

Λd,l|l−1 ξˆd,l|l−1

165

⎞ ⎡ ⎤ y1,l ⎟ ⎢ y2,l ⎥⎟ ⎢ ⎥ × ⎣ . ⎦⎟ .. ⎟ ⎠

0

0 R2,l

···

··· ..

.

0

⎤−1

.. ⎥ . ⎥ ⎥ ⎥ ⎦

RN ,l

≈

N ∑

y˜ Ti,l (Hi,l Pd,l|l−1 HiT,l + bRi,l )−1 y˜ i,l .

(64)

i=1

Substitute (63) to (60), we have ac ,l = ad,l .

(65)

166

L. Yan, C. Di, Q.M.J. Wu et al. / ISA Transactions 101 (2020) 160–169

Substitute (65) to (58) and (59), we have

ξˆc ,l|l = Λ c ,l =

1 ad,l 1 ad,l

N

[ξˆd,l|l−1 +

1∑

[Λd,l|l−1 +

b

(ai,l ξˆi,l|l − ξˆd,l|l−1 )] = ξˆd,l|l ,

(66)

i=1

N 1∑ (ai,l Λi,l − Λd,l|l−1 )] = Λd,l . b

(67)

i=1

The distributed fusion estimation is given by 1ˆ ˆ xˆ d,l|l = Λ− d,l ξd,l|l = xc ,l|l ,

Pd,l|l =

1 Λ− d,l

(68)

= Pc ,l|l .

(69)

The centralized fusion estimation and the distributed fusion estimation start from the same initial state x0 that meets (5), so, we have xˆ d,0|0 = xˆ 0|0 = xˆ c ,0|0 and Pd,0|0 = P0|0 = Pc ,0|0 . Similar to the above proof, it can be easily verified that xˆ d,1|1 = xˆ c ,1|1 and Pd,1|1 = Pc ,1|1 . This completes the proof. From the proof of Theorem 3, if ∆d,l is simplified by the use of (64) to replace (63), a suboptimal distributed fusion algorithm is obtained. In the sequel, we will give another suboptimal distributed fusion algorithm, which is easier to carry out in real applications. Corollary 2 (Suboptimal Distributed Fusion). For the linear dynamic system (1)–(5), when νi = ν1 for all i = 2, 3, . . . , N, the state estimation by suboptimal distributed fusion of sensors 1 to N can be computed by

⎧ 1 ˆ ⎪ xˆ sd,l|l = Λ− ⎪ sd,l ξsd,l|l ⎪ ⎪ ⎪ −1 ⎪ ⎪ ⎪Psd,l|l = Λsd,l ⎪ N ⎪ ∑ ⎪ ⎪ 1 1 ⎪ [ Λ + (ai,l Λi,l − Λi,l|l−1 )] Λ = ⎪ sd , l | l − 1 sd , l ⎪ asd,l b ⎪ ⎪ i=1 ⎪ ⎪ ⎪ N ⎪ ∑ ⎪ ⎪ 1 1 ˆ ˆ ⎪ ξ = ( ξ + [ai,l ξˆi,l|l − ξˆi,l|l−1 ]) ⎪ ⎨ sd,l|l asd,l sd,l|l−1 b i=1

1 −1 −T ˆ ⎪ ξˆsd,l|l−1 = [I − Asd,l−1 (Asd,l−1 + Q¯ l− ⎪ −1 ) ]Fl−1 ξsd,l−1|l−1 ⎪ ⎪ ⎪ ⎪Λsd,l|l−1 = [I − Asd,l−1 (Asd,l−1 + Q¯ −1 )−1 ]Asd,l−1 ⎪ l−1 ⎪ ⎪ ⎪ (ν0 −2)(ν0 +∆2sd,l ) ⎪ ⎪ a = sd , l ⎪ ν0 (ν0 +m−2) ⎪ ⎪ ⎪ −1 T ⎪ ⎪ A = Fl− sd , l − 1 −1 Λsd,l−1 Fl−1 ⎪ ⎪ ⎪ N ⎪ ∑ ⎪ ⎪ 2 ⎪ y˜ Ti,l (Hi,l Pi,l|l−1 HiT,l + bRi,l )−1 y˜ i,l , ⎩∆sd,l =

(70)

Remark 2. From Lemma 3 and Theorem 1, it can be easily verified that when ν0 , νw and ν1 approach to infinity, the algorithm given in Lemma 3 reduces to the classical Kalman filter, and formulas in Theorem 1 reduce to the classical information filter. At the same time, the algorithms given in Theorems 2, 3, and Corollary 2 reduce to the centralized batch fusion, the optimal distributed fusion, and the suboptimal distributed fusion algorithms of Gaussian driven systems, respectively. Therefore, the algorithms derived in this article based on t distribution are the generalization of the traditional ones that based on Gaussian Kalman filter. Actually, it is well known that the t distribution is quite similar to Gaussian distribution when the dof large enough. So, to better formulate the heavy-tailed noise, the dof of the t distribution should not be too large. To preserve the heavy-tailed property, the moment matching approach is used [23,28]. Remark 3. In Eq. (4), if νi ̸ = ν1 for i = 2, 3, . . . , N, we may use the moment matching approach to handle the corresponding information fusion estimation problem. To best preserve the heavy-tailed property, let ν ∗ = min{νi , i = 1, 2, . . . , N } [23]. By the use of the moment matching approach, i.e., approximate p(vi,l ) = St(vi,l ; 0, Ri,l , νi ) by p(vi∗,l ) = St(vi∗,l ; 0, R∗i,l , ν ∗ ), where (ν ∗ −2)ν

R∗i,l = (ν −2)ν ∗i Ri,l , i = 1, 2, . . . , N, then vi,l and vi∗,l have the same i mean and covariance. By the use of ν ∗ to replace ν1 , and R∗i,l to replace Ri,l in Theorem 2, Corollary 1, Theorem 3, and Corollary 2, we have the t distribution filter based centralized fusion, information filter based centralized fusion, the optimal distributed fusion, and the suboptimal distributed fusion algorithms for systems disturbed by heavy-tailed noises with different dofs, respectively. Remark 4. For the centralized fusion, the system parameters Fl , Ql , Hi,l , Ri,l , ν0 and the measurements yi,l for sensors i = 1, 2, . . . , N are sent to the fusion center. For the distributed fusion, both the measurements and the local estimations are sent to the fusion center. For the distributed fusion with feedback presented in Theorem 3, the fusion estimation is required to be transmitted to the local estimation too. So, the distributed fusion need much communication cost for the heavy-tailed noise systems. However, compared to the centralized fusion, the distributed fusion has the same estimation accuracy, but can improve the system’s robustness and the reliability of the estimation result in a hostile environment. Therefore, the distributed fusion is more practical for real applications. 5. Simulation

i=1

∑N

¯ where m = i=1 mi , b and Ql are the same as Lemma 3, the subscript sd denotes the suboptimal distributed fusion, the initial state meets (5). For i = 1, 2, . . . , N, the local state estimation and the local state prediction are computed by information filter of sensor i, i.e., Eq. (50). Remark 1. From Lemma 3, one can see that when the dof of the initial state, the process noise and the measurement noise are equal, the only difference between the t distribution based filter and Gaussian based classical Kalman filter lies in the computation of scale matrix Pl|l or the covariance of the state estimation error. This result is similar to many adaptive filters, such as, the strong-tracking-filter (STF), which is given based on Gaussian assumption but forced the alteration of Pl|l to improve the robustness of the classical Kalman filter under model inaccuracy [33]. The difference between this paper and the STF or the related work is that the parameter that determines the alteration of Pl|l is deduced under the formulation of t distribution.

We use an example to show the feasibility and the effectiveness of the presented algorithms. Considering a three-dimensional linear system with two sensors: xl+1 = Fl xl + wl ,

(71)

yl = Hi,l xl + vi,l , i = 1, 2, 3,

(72)

where

⎡

Ts2 2

⎤

1

Ts

Fl = ⎣0

1

Ts ⎦ ,

0

0

1

⎢

⎥

where, Ts is the sampling interval which is taken value as 0.01s. State xl = [sl s˙l s¨l ]T , where sl , s˙l and s¨l denote position, velocity and acceleration, respectively. H1 = [1 0 0], H2 = [0 1 0], H3 = [0 0 1]. The initial state is generated according to: p(x0 ) = St(x0 ; xˆ 0|0 , P0|0 , ν0 ),

(73)

L. Yan, C. Di, Q.M.J. Wu et al. / ISA Transactions 101 (2020) 160–169

167

Table 1 Average RMSEs with different filtering algorithms by single sensor. Algorithms

Position

Velocity

Acceleration

G-S1 H-S1 G-CF H-DF/CF H-SDF

0.4739 0.4470 0.4045 0.3687 0.3777

0.3815 0.3518 0.3511 0.3240 0.3340

0.3153 0.2997 0.3144 0.2927 0.2947

Table 2 Average CPU time per Monte Carlo run of different filtering algorithms.

Fig. 2. RMSEs of the position.

Fig. 3. RMSEs of the velocity.

where xˆ 0|0 = [10 0 0]T , P0|0 = diag{1, 1, 1} and ν0 = 3. The heavy-tailed noise wl and vi,l are generated according to: p(wl ) = St(wl ; 0, Q , νw ),

(74)

p(vi,l ) = St(vi,l ; 0, Ri , ν1 ), i = 1, 2, 3,

(75)

where Q = 0.16 × diag[1 1 1], R1 = 0.36, R2 = 0.25, R3 = 0.20 and νw = ν1 = 3. To evaluate performance of different algorithms, the root mean square error (RMSE) is employed as:

  L 1 ∑ i RMSE = √ (xl − xˆ il|l )2 , L

(76)

i=1

where xil and xˆ il|l denote the original signal and its estimate of the ith run, respectively. For each simulation, L = 1000 Monte Carlo simulations are run to get the RMSEs of the state estimates. Figs. 2 and 3 present the RMSE of the position and velocity by different algorithms, where ‘‘G-S1’’ and ‘‘G-CF’’ denote the RMSEs using the Gaussian information filter of Sensor 1 and the Gaussian information filter based centralized batch fusion that ν Q and regard the noises as Gaussian with covariances being ν− 2 ν R , respectively. ‘‘H-S1’’ and ‘‘H-S2’’ denote the RMSEs by the ν−2 i use of the t distribution based information filter of Sensor 1

Algorithms

H-S1

H-S2

H-S3

G-S1

G-S2

G-S3

CPU time (ms)

18.25

18.11

18.66

13.19

13.36

13.53

Algorithms

H-CF

H-DF

H-SDF

G-CF

CPU time (ms)

19.61

20.01

19.37

16.72

and Sensor 2, respectively. ‘‘H-CF’’, ‘‘H-DF’’ and ‘‘H-SDF’’ denote the RMSEs by using centralized fusion, distributed fusion and suboptimal distributed fusion, respectively. From Figs. 2 and 3, one can get that: (1) the proposed algorithms are more effective for estimation of both position and velocity than the Gaussian based algorithms when the noises are heavy-tailed; (2) the fusion results are better than the single sensor’s; (3) the H-DF and the H-CF have equivalent estimation accuracy, which are a little bit superior to the H-SDF. To better evaluate the proposed algorithms, the averaged RMSE by different algorithms are listed in Table 1, where H-Si denote the results by the use of the presented information filter of Sensor i and G-Si denote the results by the Gaussian information filter of Sensor i, i = 1, 2, 3, respectively. ‘‘G-CF’’ denotes the results by the Gaussian Kalman filter based centralized fusion. ‘‘H-CF’’, ‘‘H-DF’’ and ‘‘H-SDF’’ denote the results by the presented t distribution based centralized fusion, distributed fusion, and the suboptimal distributed fusion, respectively. One can find that the t distribution based algorithms have smaller averaged RMSEs than the Gaussian distribution based algorithms for each single sensor case and for the fusion case. Therefore, the effectiveness of the presented algorithms has been verified. For further evaluation of the performance of different algorithms, we show the CPU time of different algorithms in Table 2. One can see that the Gaussian distribution based algorithms have faster computation speed than the corresponding t distribution based algorithms, and this is consistent to our theory analysis that the computation of t distribution based algorithms is a little more complex than the Gaussian distribution based algorithms. As to different t distribution based fusion algorithms, the H-SDF shows shortest CPU time and followed by the H-CF and then the H-DF. The H-CF seems more valuable than the H-DF since they are equivalent in estimation accuracy but shows less CPU time. However, in practical applications, the measurements of different sensors can hardly be obtained at exactly the same time, so the H-DF seems to be more efficient in practical applications. The H-SDF is promising in practice since it runs fastest among all the t distribution based fusion algorithms, and has comparable estimation accuracy. According to the figures and tables in the simulation, all the heavy-tailed algorithms have higher estimation accuracy compared with the Gaussian Kalman filter based algorithms when dealing with the estimation problem of multisensor systems with heavy-tailed noises. But the Gaussian fusion algorithm has faster computing speed than the heavy-tailed fusion algorithms. Among three heavy-tailed fusion algorithms, the H-CF and the H-DF are shown to be equivalent and can obtain the most accurate estimations. The feedback is applied in the H-DF, which leads to an increasement of communication cost, thus results in longer

168

L. Yan, C. Di, Q.M.J. Wu et al. / ISA Transactions 101 (2020) 160–169

CPU time compared with the H-CF. The H-SDF is a simplification of the H-DF, which decreased the computation time with little sacrifices of accuracy.

6. Conclusions Based on multivariate t distribution, we derive the information filter, the centralized batch fusion and the distributed fusion algorithms for state estimation of heavy-tailed noise systems. The shown theoretical deductions and conducted simulations draw the following conclusions: (1) the presented t distribution based centralized fusion algorithm is effective, when the noise is heavytailed, it is superior to the classical Kalman filter based centralized fusion; (2) the t distribution based optimal distributed fusion algorithm is effective, and it is verified that they are equivalent to each other in the sense of LMSE; (3) the t distribution based suboptimal distributed fusion algorithm is effective, which is only a little worse than the optimal distributed fusion algorithm but has less computation complexity; (4) the traditional algorithms, including the classical Kalman filter, the information filter, the Kalman filter based centralized fusion and the optimal distributed fusion algorithms are the limit case of the corresponding algorithms presented in this paper on condition that the dof of t distribution tend to infinity. Thus, the proposed algorithms have potential values in a host of applications, such as aerospace, automation, and detection. There are still many interesting problems about the fusion estimation for the systems with heavy-tailed noises that can be studied, such as the fusion estimation for asynchronous multirate multisensor systems with heavy-tailed noises, the filter and the fusion estimation for the dynamic systems with correlated heavytailed noises and the distributed fusion estimation of nonlinear systems with heavy-tailed noises. These are our future work.

Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgment This work was supported by Beijing Natural Science Foundation, China under Grant 4202071 and in part by the National Natural Science Foundation of China under Grant 61703019.

References [1] Edward W, James L. Multisensor data fusion. Boston: Artech House; 1990. [2] Bar-Shalom Y, Li XR, Kirubarajan T. Estimation with Applications to Tracking and Navigation: Theory Algorithm and Software. Wiley-Interscience; 2001. [3] Dan S. Optimal state estimation: Kalman, H∞ , and nonlinear approaches. New York: John Wiley and Sons, Inc. Publication; 2006. [4] Sun SL. Optimal linear filters for discrete-time systems with randomly delayed and lost measurements with/without time stamps. IEEE Trans Automat Control 2013;58(6):1551–6. [5] Elliott RJ, van der Hoek J. Optimal linear estimation and data fusion. IEEE Trans Automat Control 2006;51(4):686–9.

[6] Yan LP, Li XR, Xia YQ, Fu MY. Optimal sequential and distributed fusion for state estimation in cross-correlated noise. Automatica 2013;49(12):3607–12. [7] Song E, Xu J, Zhu Y. Optimal distributed Kalman filtering fusion with singular covariances of filtering errors and measurement noises. IEEE Trans Automat Control 2014;59(5):1271–82. [8] Yan LP, Li XR, Xia YQ. Modeling and estimation of asynchronous multirate multisensor system with unreliable measurements. IEEE Trans Aerosp Electron Syst 2015;51(3):2012–26. [9] Hu J, Wang Z, Chen D, et al. Estimation, filtering and fusion for networked systems with network-induced phenomena: New progress and prospects. Inf Fusion 2016;31:65–75. [10] Qi H, Iyengar SS, Chakrabarty K. Distributed sensor networks–a review of recent research. J Franklin Inst B 2001;338(6):655–68. [11] Zhu J-W, Zhang W-A, Yu L, Zhang D. Robust distributed tracking control for linear multi–agent systems based on distributed intermediate estimator. J Franklin Inst B 2018;355(1):31–53. [12] Khaleghi B, Khamis A, Karray FO, et al. Multisensor data fusion: A review of the state-of-the-art. Inf Fusion 2013;14:28–44. [13] Svensen M, Bishop CM. Robust Bayesian mixture modeling. Neurocomputing 2005;64:235–52. [14] Niu W, Zhu J, Gu W, Chu J. Four statistical approaches for multisensor data fusion under non–Gaussian noise. In: International conference on control, automation and systems engineering. 2009, p. 27–30. [15] Li W, Jia M. Distributed consensus filtering for discrete-time nonlinear systems with non–Gaussian noise. Signal Process 2012;92:2464–70. [16] Shen C, Li J, Huang W. Robust centralized multi-sensor fusion using cubature information filter. In: The 30th chinese control and decision conference (2018 CCDC). 2018, p. 3297–302. [17] Wang Z, Luo X. A novel weight function–based robust iterative learning identification method for discrete box–jenkins models with student’s t–distribution noises. J Franklin Inst B 2017;354(18):8645–58. [18] Piché R, Särkkä S, Hartikainen J. Recursive outlier–robust filtering and smoothing for nonlinear systems using the multivariate student–t distribution. In: IEEE international workshop on machine learning for signal processing. Santander, Spain: IEEE; 2012, p. 1–6. [19] Roth M. On the multivariate t distribution. Technical report, Division of automatic control, Linkoping University; 2013. [20] Huang Y, Zhang Y, Li N, Chambers J. A robust Gaussian approximate filter for nonlinear systems with heavy tailed measurement noises. In: IEEE international conference on acoustics, speech and signal processing (ICASSP). 2016, p. 4209–13. [21] Huang Y, Zhang Y, Li N, Chambers J. A robust students t based cubature filter. In: 19th international conference on information fusion. Germany: Heidelberg; 2016, p. 9–16. [22] Amor N, Kahlaoui S, Chebbi S. Unscented particle filter using student–t distribution with non–Gaussian measurement noise. In: 2018 international conference on advanced systems and electric technologies. 2018, p. 34–8. [23] Roth M, O¨zkan E, Gustafsson F. A student’s filter for heavy tailed process and measurement noise. In: International conference on acoustics, speech and signal processing (ICASSP). 2013, p. 5770–4. [24] Baccelli F, Chatterjee A, Vishwanath S. Pairwise stochastic bounded confidence opinion dynamics: Heavy tails and stability. IEEE Trans Automat Control 2017;62(11):5678–93. [25] Kashima K, Aoyama H, Ohta Y. Stable process approach to analysis of systems under heavy-tailed noise: Modeling and stochastic linearizatione. IEEE Trans Automat Control 2019;64(4):1344–57. [26] Zhang L, Lan J, Li XR. A normal-Gamma filter for linear systems with heavy–tailed measurement noise. In: International conference on information fusion (FUSION). 2018. Xian, Shaanxi, p. 2552–9. [27] Agamennoni G, Nieto JI, Nebot EM. Approximate inference in statespace models with heavy–tailed noise. IEEE Trans Signal Process 2012;60(10):5024–37. [28] Huang Y, Zhang Y, Li N, Chambers J. Robust students t based nonlinear filter and smoother. IEEE Trans Aerosp Electron Syst 2016;52(5):2586–96. [29] Zhu H, Leung H, He Z. A variational Bayesian approach to robust sensor fusion based on student-t distribution. Inform Sci 2013;221:201–14. [30] Ristic B, Arulampalam S, Gordon N. Beyond the Kalman filter. Artech House; 2004. [31] Yan LP, Xia YQ, Liu B, Fu MY. Multisensor optimal estimation theory and its application. Beijing: The Science Publishing House; 2015. [32] Kotz S, Nadarajah S. Multivariate t-distributions and their applications. Cambridge University Press; 2004. [33] Xie XQ, Zhou DH, Jin YH. Strong tracking filter based adaptive generic model control. J Process Control 1999;9(4):337–50.

L. Yan, C. Di, Q.M.J. Wu et al. / ISA Transactions 101 (2020) 160–169 Liping Yan was born in Henan Province, China, in 1979. She received her B.S. degree and M.S. degree both in Mathematics from Henan University, Kaifeng city, Henan Province, P.R. China, in 2000 and 2003, respectively, and received her Ph.D. degree in Control Science and Engineering from Tsinghua University, Beijing, China, in 2007. From January 2007 to July 2009, she was a postdoctoral research associate in the Equipment Academy of Airforce. Since July 2009, she has been with the School of Automation, Beijing Institute of Technology, Beijing, first as an Assistant Professor, then, since 2011, as an Associate Professor. From March 2012 to March 2013, supported by the CSC, she was a visiting scholar in the University of New Orleans, New Orleans, LA, USA. From Sep. 2018 to Aug. 2019, supported by the CSC, she was a visiting scholar in the University of Windsor, Windsor, ON, Canada. She has co-authored five books and more than 60 journal and conference papers. Currently, she is an associate professor and Ph.D supervisor in BIT. Her research interests include multisensor data fusion, state estimation, image registration, intelligent navigation and integrated navigation, etc. Chenying Di was born in Shanxi Province, China, in 1994. She received her B.S. degree in Automation from Hebei University of Technology, Tianjin city, P.R. China, in 2017. She is currently working toward the Master degree in control science and engineering at the School of Automation in Beijing Institute of Technology, Beijing, China. Her research interests include multisensor data fusion, state estimation, etc.

Q. M. Jonathan Wu (M’92–SM’09) received the Ph.D. in electrical engineering from the University of Wales, Swansea, U.K., in 1990. From 1995 to 2005, he was affiliated with the National Research Council of Canada, Ottawa, ON, Canada, where he became a Senior Research Officer and a Group Leader. He is currently a Professor with the Department of Electrical and Computer Engineering, University of Windsor, Windsor, ON, Canada. He has published over 300 peer-reviewed papers in computer vision, image processing, intelligent systems, robotics, and integrated microsystems. His current research interests include 3-D computer vision, active video object tracking and extraction, interactive multimedia, sensor analysis and fusion, and visual sensor networks. Dr. Wu is the fellow of the Canadian Academy of Engineering, and holds the Tier 1 Canada Research Chair of Automotive Sensors and Information Systems. He is an Associate Editor of the IEEE Transactions on Neural Networks and Learning Systems, Cognitive Computation, and the International Journal of Robotics and Automation. He has served on technical program committees and international advisory committees for many prestigious conferences.

169

Yuanqing Xia was born in Anhui Province, China, in 1971 and graduated from the Department of Mathematics, Chuzhou University, Chuzhou, China, in 1991. He received his M.S. degree in Fundamental Mathematics from Anhui University, China, in 1998 and his Ph.D. degree in Control Theory and Engineering from Beijing University of Aeronautics and Astronautics, Beijing, China, in 2001. From 1991–1995, he was with Tongcheng Middle-School, Anhui, China, where he worked as a teacher. During January 2002 to November 2003, he was a postdoctoral research associate in the Institute of Systems Science, Academy of Mathematics and System Sciences, Chinese Academy of Sciences, Beijing, China, where he worked on navigation, guidance and control. From November 2003 to February 2004, he was with the National University of Singapore as a research fellow, where he worked on variable structure control. From February 2004 to February 2006, he was with the University of Glamorgan, Pontypridd, U.K., as a Research Fellow, where he worked on networked control systems. From February 2007 to June 2008, he was a Guest Professor with Innsbruck Medical University, Innsbruck, Austria, where he worked on biomedical signal processing. Since July 2004, he has been with the Department of Automation, Beijing Institute of Technology, Beijing, first as an Associate Professor, then, since 2008, as a Professor and in 2012, he was appointed as Xu Teli distinguished professor in Beijing Institute of Technology and obtained National Science Foundation for Distinguished Young Scholars of China. His current research interests are in the fields of networked information fusion and control systems, robust control, active disturbance rejection control and flight control. He has published 10 monographs in Springer and John Wiley, etc, and more than 100 papers in journals. He is a deputy editor of Journal of Beijing Institute of Technology, associate editor of Acta Automatica Sinica, Control Theory and Applications, International Journal of Innovative Computing, Information and Control. He obtained Second Award of Beijing Municipal Science and Technology (No.1) in 2010 and 2015, respectively, Second National Award for Science and Technology (No.2) in 2011, and second natural science award of The Ministry of Education (No.1) in 2012. Shida Liu was born in Mongolia province in 1988. He received the bachelor’s degree from Inner Mongolia University, Hohhot, China, in 2011, and the Ph.D. degree from Beijing Jiaotong University, Beijing, China, in 2017. He was a Post-Doctoral Researcher with the School of Automation Science and Electrical Engineering, Beihang University, Beijing from 2017 to 2019. Until now, he is a lecture in School of Electrical and Control Engineering, North China University of Technology, Beijing 10093, People’s Republic of China. His current research interests include intelligent traffic control, data-driven control, and artificial intelligence.