**Consensus-Based Distributed Filtering for GNSS** Consensus-Based Distributed Filtering for GNSS

DOI: 10.5772/intechopen.71138

Amir Khodabandeh, Peter J.G. Teunissen and Safoora Zaminpardaz Amir Khodabandeh, Peter J.G. Teunissen and Safoora Zaminpardaz

Additional information is available at the end of the chapter Additional information is available at the end of the chapter

http://dx.doi.org/10.5772/intechopen.71138

### Abstract

[6] Chong CY, Mori S. Optimal Fusion For Non-Zero Process Noise. Proceedings of the 16th

[7] Koch W, Govaers F. On accumulated state densities with applications to out-of-sequence measurement processing. In: IEEE Transactions on Aerospace and Electronic Systems.

[8] Reinhardt M. Linear Estimation in Interconnected Sensor Systems with Information

[9] Govaers F, Chong CY, Mori S, Koch W. Comparison of Augmented State Track Fusion Methods for Non-Full-Rate Communication. 18th International Conference on Informa-

[10] Pfaff F, Noack B, Hanebeck UD, Govaers F, Koch W. Information Form Distributed Kalman Filtering (IDKF) with Explicit Inputs. 20th International Conference on Informa-

[11] Koch W. Tracking and Sensor Data Fusion: Methodological Framework and Selected

International Conference on Information Fusion, Istanbul. 2013

Costraints. PhD thesis: KIT Scientific Publishing; 2014

tion Fusion, Washington, DC. 2015

tion Fusion, Xi'an, China. 2017

Applications. Springer; 2014:251

Vol. 47. 2011

272 Kalman Filters - Theory for Advanced Applications

Kalman filtering in its distributed information form is reviewed and applied to a network of receivers tracking Global Navigation Satellite Systems (GNSS). We show, by employing consensus-based data-fusion rules between GNSS receivers, how the consensus-based Kalman filter (CKF) of individual receivers can deliver GNSS parameter solutions that have a comparable precision performance as their network-derived, fusion center dependent counterparts. This is relevant as in the near future the proliferation of low-cost receivers will give rise to a significant increase in the number of GNSS users. With the CKF or other distributed filtering techniques, GNSS users can therefore achieve highprecision solutions without the need of relying on a centralized computing center.

Keywords: distributed filtering, consensus-based Kalman filter (CKF), global navigation satellite systems (GNSS), GNSS networks, GNSS ionospheric observables

### 1. Introduction

Kalman filtering in its decentralized and distributed forms has received increasing attention in the sensor network community and has been extensively studied in recent years, see e.g. [1–8]. While in the traditional centralized Kalman filter setup all sensor nodes have to send their measurements to a computing (fusion) center to obtain the state estimate, in the distributed filtering schemes the nodes only share limited information with their neighboring nodes (i.e. a subset of all other nodes) and yet obtain state estimates that are comparable to that of the centralized filter in a minimum-mean-squared-error sense. This particular feature of the distributed filters would potentially make the data communication between the nodes costeffective and develop the nodes' capacity to perform parallel computations.

Next to sensor networks, distributed filtering can therefore benefit several other applications such as formation flying of aerial vehicles [9], cooperative robotics [10] and disciplines that concern the Global Navigation Satellite Systems (GNSS). The latter is the topic of this present

© 2018 The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

© The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and eproduction in any medium, provided the original work is properly cited.

contribution. The GNSS have been proven to be an efficient tool for determination of time varying parameters that are of importance for Earth science disciplines like positioning, deformation, timing and atmosphere [11, 12]. Parameter estimation in GNSS often relies on the data processing of a network of receivers that collect measurements from visible GNSS satellites. In the context of sensor networks, GNSS network receivers therefore serve as sensor nodes, providing their data to a computing center, thereby computing network-based parameter solutions in a (near) real-time manner. In this contribution we intend to demonstrate how consensus algorithms [13] and the corresponding consensus-based Kalman filter (CKF), as a popular means for distributed filtering, can take an important role in GNSS applications for which a network of receivers are to be processed. Although each single receiver can run its own local filter to deliver GNSS-derived solutions, the precision of such single-receiver solutions is generally much lower than its network-derived counterparts, see e.g. [14, 15]. It will be shown, through a CKF setup, that single-receiver parameter solutions can achieve precision performances similar to that of their network-based versions, provided that a sufficient number of iterative communications between the neighboring receivers are established. The importance of such consensus-based single-receiver solutions is well appreciated in the light of the recent development of new GNSS constellations as well as the proliferation of low-cost massmarket receivers [16–18]. With the increase in the number and types of GNSS receivers, many more GNSS users can establish their own measurement setup to determine parameters that suit their needs. By taking recourse to the CKF or other distributed filtering techniques, GNSS users can therefore potentially deliver high-precision parameter solutions without the need of having a computing center.

means of the state-vectors x1, …, xt can be known, while their unknown realizations still need

on which set of observables prediction is based, we use the notation <sup>b</sup>xt∣<sup>τ</sup> as the predictor of xt

denoted by Eð Þ: , Cð Þ :; : and Dð Þ: , respectively. The capital Q is reserved for (co)variance matri-

To predict the state-vectors in an optimal sense, one often uses the minimum mean squared error (MMSE) principle as the optimality criterion, see e.g., [19–25]. In case no restrictions are placed on the class of predictors, the MMSE predictor <sup>b</sup>xt∣<sup>τ</sup> is given by the conditional mean

� �, known as the Best Predictor (BP). The BP is unbiased, but generally nonlinear, with

exemptions, for instance in the Gaussian case. In case xt and <sup>y</sup>½ � <sup>τ</sup> are jointly Gaussian, the BP becomes linear and identical to its linear counterpart, i.e. the Best Linear Predictor (BLP)

of a BLP is always uncorrelated with observables on which the BLP is based, i.e. C xt � <sup>b</sup>xt∣<sup>τ</sup>, � <sup>y</sup>½ � <sup>τ</sup> Þ ¼ 0. These two basic properties can be alternatively used to uniquely specify a BLP [26]. The Kalman filter is a recursive BP (Gaussian case) or a recursive BLP. A recursive predictor, say <sup>b</sup>xt∣<sup>t</sup>, can be obtained from the previous predictor <sup>b</sup>xt∣t�<sup>1</sup> and the newly collected observable

determination of temporally varying parameters. We now state the standard assumptions that

The dynamic model: The linear dynamic model, describing the time-evolution of the state-

for the time instances t, s ¼ 1, 2, …, with δt,s being the Kronecker delta. The nonsingular matrix Φt,t�<sup>1</sup> denotes the transition matrix and the random vector dt is the system noise. The system noise dt is thus assumed to have a zero mean, to be uncorrelated in time and to be uncorrelated

. Recursive prediction is thus very suitable for applications that require real-time

xt ¼ Φt,t�<sup>1</sup> xt�<sup>1</sup> þ dt, t ¼ 1, 2, … (2)

Eð Þ¼ x<sup>0</sup> x<sup>0</sup>∣<sup>0</sup>, Dð Þ¼ x<sup>0</sup> Qx0x<sup>0</sup> (3)

Eð Þ¼ dt 0, C dt ð Þ¼ ; ds St δt,s, C dt ð Þ¼ ; x<sup>0</sup> 0 (4)

<sup>y</sup>½ � <sup>τ</sup> <sup>y</sup>½ � <sup>τ</sup> <sup>y</sup>½ � <sup>τ</sup> � <sup>E</sup> <sup>y</sup>½ � <sup>τ</sup>

n o � � (1)

� � <sup>¼</sup> <sup>E</sup>ð Þ xt , and that (2) the prediction error

<sup>b</sup>xt∣<sup>τ</sup> <sup>¼</sup> E xð Þþ<sup>t</sup> Qxty½ � <sup>τ</sup> <sup>Q</sup>�<sup>1</sup>

. The expectation, covariance and dispersion operators are

. In the following, to show

275

Consensus-Based Distributed Filtering for GNSS http://dx.doi.org/10.5772/intechopen.71138

to be guessed (predicted) through observed realizations of y1, …, yt

<sup>1</sup> ; …; yT τ � �<sup>T</sup>

.

Eq. (1) implies that (1) the BLP is unbiased, i.e. E <sup>b</sup>xt∣<sup>τ</sup>

make the Kalman filter recursion feasible.

2.1. The Kalman filter standard assumptions

when based on <sup>y</sup>½ � <sup>τ</sup> <sup>¼</sup> yT

� � <sup>¼</sup> Qxty<sup>τ</sup>

ces. Thus C xt; y<sup>τ</sup>

E xtjy½ � <sup>τ</sup>

vector yt

with

and

vectors xt, is given as

The structure of this contribution is as follows. We first briefly review the principles of the standard Kalman filter and its information form in Section 2. The additivity property of the information filter that makes this filter particularly useful for distributed processing is also highlighted. In Section 3 we discuss average consensus rules on which the sensor nodes agree to fuse each other information. Different consensus protocols are discussed and a 'probabilistic' measure for the evaluation of their convergence rates is proposed. Section 4 is devoted to the CKF algorithmic steps. Its two time-scale nature is remarked and a three-step recursion for evaluating the consensus-based error variance matrix is developed. In Section 5 we apply the CKF theory to a small-scale network of GNSS receivers collecting ionospheric observables over time. Conducting a precision analysis, we compare the precision of the network-based ionospheric solutions with those of their single-receiver and consensus-based counterparts. It is shown how the CKF of each receiver responses to an increase in the number of iterative communications between the neighboring nodes. Concluding remarks and future outlook are provided in Section 6.

### 2. Kalman filtering

Consider a time series of observable random vectors y1,…, yt . The goal is to predict the unobservable random state-vectors x1, …, xt. By the term 'prediction', we mean that the observables y1, …, yt are used to estimate realizations of the random vectors x1, …, xt. Accordingly, the means of the state-vectors x1, …, xt can be known, while their unknown realizations still need to be guessed (predicted) through observed realizations of y1, …, yt . In the following, to show on which set of observables prediction is based, we use the notation <sup>b</sup>xt∣<sup>τ</sup> as the predictor of xt when based on <sup>y</sup>½ � <sup>τ</sup> <sup>¼</sup> yT <sup>1</sup> ; …; yT τ � �<sup>T</sup> . The expectation, covariance and dispersion operators are denoted by Eð Þ: , Cð Þ :; : and Dð Þ: , respectively. The capital Q is reserved for (co)variance matrices. Thus C xt; y<sup>τ</sup> � � <sup>¼</sup> Qxty<sup>τ</sup> .

#### 2.1. The Kalman filter standard assumptions

To predict the state-vectors in an optimal sense, one often uses the minimum mean squared error (MMSE) principle as the optimality criterion, see e.g., [19–25]. In case no restrictions are placed on the class of predictors, the MMSE predictor <sup>b</sup>xt∣<sup>τ</sup> is given by the conditional mean E xtjy½ � <sup>τ</sup> � �, known as the Best Predictor (BP). The BP is unbiased, but generally nonlinear, with exemptions, for instance in the Gaussian case. In case xt and <sup>y</sup>½ � <sup>τ</sup> are jointly Gaussian, the BP becomes linear and identical to its linear counterpart, i.e. the Best Linear Predictor (BLP)

$$
\widehat{\mathbf{x}}\_{t|\tau} = \mathbf{E}(\mathbf{x}\_t) + \mathbf{Q}\_{\mathbf{x}\_t y\_{[\tau]}} \mathbf{Q}\_{y\_{[\tau]} y\_{[\tau]}}^{-1} \left\{ y\_{[\tau]} - \mathbf{E} \left( y\_{[\tau]} \right) \right\} \tag{1}
$$

Eq. (1) implies that (1) the BLP is unbiased, i.e. E <sup>b</sup>xt∣<sup>τ</sup> � � <sup>¼</sup> <sup>E</sup>ð Þ xt , and that (2) the prediction error of a BLP is always uncorrelated with observables on which the BLP is based, i.e. C xt � <sup>b</sup>xt∣<sup>τ</sup>, � <sup>y</sup>½ � <sup>τ</sup> Þ ¼ 0. These two basic properties can be alternatively used to uniquely specify a BLP [26].

The Kalman filter is a recursive BP (Gaussian case) or a recursive BLP. A recursive predictor, say <sup>b</sup>xt∣<sup>t</sup>, can be obtained from the previous predictor <sup>b</sup>xt∣t�<sup>1</sup> and the newly collected observable vector yt . Recursive prediction is thus very suitable for applications that require real-time determination of temporally varying parameters. We now state the standard assumptions that make the Kalman filter recursion feasible.

The dynamic model: The linear dynamic model, describing the time-evolution of the statevectors xt, is given as

$$\mathbf{x}\_{t} = \Phi\_{t, t-1} \mathbf{x}\_{t-1} + d\_{t} \quad \text{ } t = 1, 2, \dots \tag{2}$$

with

contribution. The GNSS have been proven to be an efficient tool for determination of time varying parameters that are of importance for Earth science disciplines like positioning, deformation, timing and atmosphere [11, 12]. Parameter estimation in GNSS often relies on the data processing of a network of receivers that collect measurements from visible GNSS satellites. In the context of sensor networks, GNSS network receivers therefore serve as sensor nodes, providing their data to a computing center, thereby computing network-based parameter solutions in a (near) real-time manner. In this contribution we intend to demonstrate how consensus algorithms [13] and the corresponding consensus-based Kalman filter (CKF), as a popular means for distributed filtering, can take an important role in GNSS applications for which a network of receivers are to be processed. Although each single receiver can run its own local filter to deliver GNSS-derived solutions, the precision of such single-receiver solutions is generally much lower than its network-derived counterparts, see e.g. [14, 15]. It will be shown, through a CKF setup, that single-receiver parameter solutions can achieve precision performances similar to that of their network-based versions, provided that a sufficient number of iterative communications between the neighboring receivers are established. The importance of such consensus-based single-receiver solutions is well appreciated in the light of the recent development of new GNSS constellations as well as the proliferation of low-cost massmarket receivers [16–18]. With the increase in the number and types of GNSS receivers, many more GNSS users can establish their own measurement setup to determine parameters that suit their needs. By taking recourse to the CKF or other distributed filtering techniques, GNSS users can therefore potentially deliver high-precision parameter solutions without the need of

The structure of this contribution is as follows. We first briefly review the principles of the standard Kalman filter and its information form in Section 2. The additivity property of the information filter that makes this filter particularly useful for distributed processing is also highlighted. In Section 3 we discuss average consensus rules on which the sensor nodes agree to fuse each other information. Different consensus protocols are discussed and a 'probabilistic' measure for the evaluation of their convergence rates is proposed. Section 4 is devoted to the CKF algorithmic steps. Its two time-scale nature is remarked and a three-step recursion for evaluating the consensus-based error variance matrix is developed. In Section 5 we apply the CKF theory to a small-scale network of GNSS receivers collecting ionospheric observables over time. Conducting a precision analysis, we compare the precision of the network-based ionospheric solutions with those of their single-receiver and consensus-based counterparts. It is shown how the CKF of each receiver responses to an increase in the number of iterative communications between the neighboring nodes. Concluding remarks and future outlook are

unobservable random state-vectors x1, …, xt. By the term 'prediction', we mean that the observables y1, …, yt are used to estimate realizations of the random vectors x1, …, xt. Accordingly, the

. The goal is to predict the

having a computing center.

274 Kalman Filters - Theory for Advanced Applications

provided in Section 6.

2. Kalman filtering

Consider a time series of observable random vectors y1,…, yt

$$\mathbf{E}(\mathbf{x}\_0) = \mathbf{x}\_{0\mid 0} \quad \text{D}(\mathbf{x}\_0) = \mathbf{Q}\_{\mathbf{x}\_0 \mathbf{x}\_0} \tag{3}$$

and

$$\mathbf{E}(d\_t) = \mathbf{0}, \quad \mathbf{C}(d\_t, d\_s) = \mathbf{S}\_t \delta\_{t, s\prime} \quad \mathbf{C}(d\_t, \mathbf{x}\_0) = \mathbf{0} \tag{4}$$

for the time instances t, s ¼ 1, 2, …, with δt,s being the Kronecker delta. The nonsingular matrix Φt,t�<sup>1</sup> denotes the transition matrix and the random vector dt is the system noise. The system noise dt is thus assumed to have a zero mean, to be uncorrelated in time and to be uncorrelated with the initial state-vector x0. The transition matrix from epoch s to t is denoted as Φt,s. Thus Φ�<sup>1</sup> t,s ¼ Φs,t and Φt,t ¼ I (the identity matrix).

The measurement model: The link between the observables yt and the state-vectors xt is assumed given as

$$\mathbf{y}\_t = A\_t \mathbf{x}\_t + \varepsilon\_t \quad t = 1, 2, \dots \tag{5}$$

Measurement update: In the presence of new data yt

been met if the most recent data yt of <sup>y</sup>½ �<sup>t</sup> <sup>¼</sup> yT

the innovation of yt

since C xt � <sup>b</sup>xt∣t�<sup>1</sup>; yt

since C xt � <sup>b</sup>xt∣t�<sup>1</sup>; <sup>ε</sup><sup>t</sup>

� � <sup>¼</sup> 0.

gain matrix, is uniquely specified by

� � <sup>¼</sup> Pt∣t�<sup>1</sup>AT

The error variance matrix Pt∣<sup>t</sup> <sup>¼</sup> <sup>D</sup> xt � <sup>b</sup>xt∣<sup>t</sup>

agation law, together with C xt � <sup>b</sup>xt∣t�<sup>1</sup>; vt

2.3. A remark on the filter initialization

implies that yt <sup>¼</sup> Atbxt∣t�1, where Atbxt∣t�<sup>1</sup> is the BLP of yt

Kt <sup>¼</sup> Pt∣t�<sup>1</sup>AT

tion law to (10) gives the variance matrix of vt as follows

<sup>t</sup> Q�<sup>1</sup>

<sup>t</sup> and C vt; yt

� � <sup>¼</sup> Pt∣t�<sup>1</sup>A<sup>T</sup>

prediction error xt � <sup>b</sup>xt∣t�<sup>1</sup> and the measurement noise <sup>ε</sup>t, i.e. (cf. 5)

the BLP <sup>b</sup>xt∣<sup>t</sup>. Such a candidate fulfills the unbiasedness condition E <sup>b</sup>xt∣t�<sup>1</sup>

necessarily the zero-correlation condition, that is, Cðxt � <sup>b</sup>xt∣t�<sup>1</sup>, y½ �<sup>t</sup> <sup>Þ</sup> 6¼ 0. Note, however, that <sup>C</sup>ðxt � <sup>b</sup>xt∣t�<sup>1</sup>, y½ � <sup>t</sup>�<sup>1</sup> Þ ¼ 0. Thus the zero-correlation condition Cðxt � <sup>b</sup>xt∣t�<sup>1</sup>, y½ �<sup>t</sup> Þ ¼ 0 would have

data <sup>y</sup>½ � <sup>t</sup>�<sup>1</sup> , thereby fully predicted by <sup>y</sup>½ � <sup>t</sup>�<sup>1</sup> . Since an observable is its own best predictor, this

quantity vt <sup>¼</sup> yt � Atbxt∣t�<sup>1</sup> to be identically zero which is generally not the case. It is therefore the presence of vt that violates the zero-correlation condition. Note that vt is a function of the

that are both uncorrelated with <sup>y</sup>½ � <sup>t</sup>�<sup>1</sup> . Therefore, vt cannot be predicted by the previous data <sup>y</sup>½ � <sup>t</sup>�<sup>1</sup> , showing that vt contains truly new information. That is why vt is sometimes referred to as

linear function of vt to it. It reads <sup>b</sup>xt∣<sup>t</sup> <sup>¼</sup> <sup>b</sup>xt∣t�<sup>1</sup> <sup>þ</sup> Ktvt, with Kt being an unknown matrix to be chosen such that the zero-correlation condition is met. Such a matrix, known as the Kalman

� � <sup>¼</sup> Qvtvt

<sup>b</sup>xt∣<sup>t</sup> <sup>¼</sup> <sup>b</sup>xt∣t�<sup>1</sup> <sup>þ</sup> Ktvt, with Pt∣<sup>t</sup> <sup>¼</sup> Pt∣t�<sup>1</sup> � KtQvtvt

Qvtvt <sup>¼</sup> AtPt∣t�<sup>1</sup>A<sup>T</sup>

In the derivation of the Kalman filter one assumes the mean of the random initial state-vector x0, in Eq. (3), to be known, see e.g. [30–37]. This is because of the BLP structure (1) that needs knowledge of the means Eð Þ xt and Eðy½ � <sup>τ</sup> <sup>Þ</sup>. Since in many, if not most, applications the means of the state-vectors x1, …, xt are unknown, such derivation is therefore not appropriate. As shown in Ref. [38], one can do away with this need to have both the initial mean x<sup>0</sup>∣<sup>0</sup> and variance

vtvt <sup>⇔</sup> <sup>C</sup> xt � <sup>b</sup>xt∣t�<sup>1</sup> � Ktvt; yt

, see e.g. [27–29]. We now amend our earlier candidate <sup>b</sup>xt∣t�<sup>1</sup> by adding a

vt <sup>¼</sup> At xt � <sup>b</sup>xt∣t�<sup>1</sup>

½ � <sup>t</sup>�<sup>1</sup> ; yT t h i<sup>T</sup>

, one may yet offer <sup>b</sup>xt∣t�<sup>1</sup> as a candidate for

Consensus-Based Distributed Filtering for GNSS http://dx.doi.org/10.5772/intechopen.71138

would be a function of the previous

. But this would require the zero-mean

� � <sup>þ</sup> <sup>ε</sup>t, (10)

� � <sup>¼</sup> <sup>0</sup> (11)

. The measurement update reads then

<sup>t</sup> . Application of the covariance propaga-

<sup>t</sup> þ Rt (13)

<sup>t</sup> (12)

KT

� � follows by an application of the covariance prop-

� � <sup>¼</sup> <sup>E</sup>ð Þ xt , but not

277

with

$$\mathbf{E}(\varepsilon\_t) = \mathbf{0}, \quad \mathbf{C}(\varepsilon\_t, \varepsilon\_s) = \mathbf{R}\_t \delta\_{t, s\prime} \quad \mathbf{C}(\varepsilon\_t, \mathbf{x}\_0) = \mathbf{0}, \quad \mathbf{C}(\varepsilon\_t, d\_s) = \mathbf{0} \tag{6}$$

for t, s ¼ 1, 2, …, with At being the known design matrix. Thus the zero-mean measurement noise ε<sup>t</sup> is assumed to be uncorrelated in time and to be uncorrelated with the initial statevector x<sup>0</sup> and the system noise dt.

#### 2.2. The three-step recursion

Initialization: As the mean of x<sup>0</sup> is known, the best predictor of x<sup>0</sup> in the absence of data is the mean Eð Þ¼ x<sup>0</sup> x<sup>0</sup>∣0. Hence, the initialization is simply given by

$$
\widehat{\mathbf{x}}\_{0|0} = \mathbf{x}\_{0|0} \quad \text{ } P\_{0|0} = Q\_{\mathbf{x}\_0 \mathbf{x}\_0} \tag{7}
$$

That the initial error variance matrix <sup>P</sup><sup>0</sup>∣<sup>0</sup> <sup>¼</sup> <sup>D</sup> <sup>x</sup><sup>0</sup> � <sup>b</sup>x<sup>0</sup>∣<sup>0</sup> � � is identical to the variance matrix Qx0x<sup>0</sup> follows from the equality D x<sup>0</sup> � x<sup>0</sup>∣<sup>0</sup> � � <sup>¼</sup> <sup>D</sup>ð Þ <sup>x</sup><sup>0</sup> .

Time update: Let us choose <sup>Φ</sup>t,t�<sup>1</sup>bxt�1∣t�<sup>1</sup> as a candidate for the BLP <sup>b</sup>xt∣t�1. According to Eq. (1), the candidate would be the BLP if it fulfills two conditions: (1) it must be unbiased and (2) it must have a prediction error uncorrelated with the previous data <sup>y</sup>½ � <sup>t</sup>�<sup>1</sup> . The first condition, i.e. EðΦt,t�<sup>1</sup>bxt�1∣t�<sup>1</sup>Þ ¼ <sup>E</sup>ð Þ xt , follows from Eq. (2) and the equalities E <sup>b</sup>xt�1∣t�<sup>1</sup> � � <sup>¼</sup> <sup>E</sup>ð Þ xt�<sup>1</sup> and <sup>E</sup>ð Þ¼ dt 0. The second condition, i.e. Cðxt � <sup>Φ</sup>t,t�<sup>1</sup>bxt�1∣t�<sup>1</sup>; <sup>y</sup>½ � <sup>t</sup>�<sup>1</sup> Þ ¼ 0, follows from the fact that the prediction error xt � <sup>Φ</sup>t,t�<sup>1</sup>bxt�1∣t�<sup>1</sup> is a function of the previous BLP prediction error xt�<sup>1</sup> � <sup>b</sup>xt�1∣t�<sup>1</sup> and the system noise dt, i.e. (cf. 2)

$$\mathbf{x}\_{t} - \Phi\_{t, t-1} \widehat{\mathbf{x}}\_{t-1|t-1} = \Phi\_{t, t-1} \left( \mathbf{x}\_{t-1} - \widehat{\mathbf{x}}\_{t-1|t-1} \right) + d\_{t\prime} \tag{8}$$

that are both uncorrelated with the previous data <sup>y</sup>½ � <sup>t</sup>�<sup>1</sup> . Hence, the time update is given by

$$
\widehat{\mathbf{x}}\_{t|t-1} = \Phi\_{t, t-1} \widehat{\mathbf{x}}\_{t-1|t-1} \quad P\_{t|t-1} = \Phi\_{t, t-1} P\_{t-1|t-1} \Phi\_{t, t-1}^T + \mathbf{S}\_t \tag{9}
$$

The error variance matrix Pt∣t�<sup>1</sup> <sup>¼</sup> <sup>D</sup> xt � <sup>b</sup>xt∣t�<sup>1</sup> � � follows by applying the covariance propagation law to (8), together with C xt�<sup>1</sup> � <sup>b</sup>xt�1∣t�<sup>1</sup>; dt � � <sup>¼</sup> 0.

Measurement update: In the presence of new data yt , one may yet offer <sup>b</sup>xt∣t�<sup>1</sup> as a candidate for the BLP <sup>b</sup>xt∣<sup>t</sup>. Such a candidate fulfills the unbiasedness condition E <sup>b</sup>xt∣t�<sup>1</sup> � � <sup>¼</sup> <sup>E</sup>ð Þ xt , but not necessarily the zero-correlation condition, that is, Cðxt � <sup>b</sup>xt∣t�<sup>1</sup>, y½ �<sup>t</sup> <sup>Þ</sup> 6¼ 0. Note, however, that <sup>C</sup>ðxt � <sup>b</sup>xt∣t�<sup>1</sup>, y½ � <sup>t</sup>�<sup>1</sup> Þ ¼ 0. Thus the zero-correlation condition Cðxt � <sup>b</sup>xt∣t�<sup>1</sup>, y½ �<sup>t</sup> Þ ¼ 0 would have been met if the most recent data yt of <sup>y</sup>½ �<sup>t</sup> <sup>¼</sup> <sup>y</sup><sup>T</sup> ½ � <sup>t</sup>�<sup>1</sup> ; yT t h i<sup>T</sup> would be a function of the previous data <sup>y</sup>½ � <sup>t</sup>�<sup>1</sup> , thereby fully predicted by <sup>y</sup>½ � <sup>t</sup>�<sup>1</sup> . Since an observable is its own best predictor, this implies that yt <sup>¼</sup> Atbxt∣t�1, where Atbxt∣t�<sup>1</sup> is the BLP of yt . But this would require the zero-mean quantity vt <sup>¼</sup> yt � Atbxt∣t�<sup>1</sup> to be identically zero which is generally not the case. It is therefore the presence of vt that violates the zero-correlation condition. Note that vt is a function of the prediction error xt � <sup>b</sup>xt∣t�<sup>1</sup> and the measurement noise <sup>ε</sup>t, i.e. (cf. 5)

$$
\sigma\_t = A\_t \left( \mathbf{x}\_t - \widehat{\mathbf{x}}\_{t|t-1} \right) + \varepsilon\_t \tag{10}
$$

that are both uncorrelated with <sup>y</sup>½ � <sup>t</sup>�<sup>1</sup> . Therefore, vt cannot be predicted by the previous data <sup>y</sup>½ � <sup>t</sup>�<sup>1</sup> , showing that vt contains truly new information. That is why vt is sometimes referred to as the innovation of yt , see e.g. [27–29]. We now amend our earlier candidate <sup>b</sup>xt∣t�<sup>1</sup> by adding a linear function of vt to it. It reads <sup>b</sup>xt∣<sup>t</sup> <sup>¼</sup> <sup>b</sup>xt∣t�<sup>1</sup> <sup>þ</sup> Ktvt, with Kt being an unknown matrix to be chosen such that the zero-correlation condition is met. Such a matrix, known as the Kalman gain matrix, is uniquely specified by

$$\mathbf{K}\_{t} = P\_{t|t-1} \mathbf{A}\_{t}^{T} \mathbf{Q}\_{v\_{t}v\_{t}}^{-1} \quad \Leftrightarrow \quad \mathbf{C} \left(\mathbf{x}\_{t} - \widehat{\mathbf{x}}\_{t|t-1} - \mathbf{K}\_{t} \boldsymbol{\upsilon}\_{t}, \boldsymbol{y}\_{t}\right) = \mathbf{0} \tag{11}$$

since C xt � <sup>b</sup>xt∣t�<sup>1</sup>; yt � � <sup>¼</sup> Pt∣t�<sup>1</sup>AT <sup>t</sup> and C vt; yt � � <sup>¼</sup> Qvtvt . The measurement update reads then

$$
\widehat{\boldsymbol{\mathfrak{X}}}\_{t|t} = \widehat{\boldsymbol{\mathfrak{X}}}\_{t|t-1} + \boldsymbol{\mathcal{K}}\_{t} \boldsymbol{\mathfrak{v}}\_{t\prime} \quad \text{with} \quad \boldsymbol{P}\_{t|t} = \boldsymbol{P}\_{t|t-1} - \boldsymbol{\mathcal{K}}\_{t} \boldsymbol{Q}\_{\boldsymbol{\mathfrak{v}}\_{t}\boldsymbol{\mathfrak{v}}\_{t}} \boldsymbol{\mathsf{K}}\_{t}^{T} \tag{12}
$$

The error variance matrix Pt∣<sup>t</sup> <sup>¼</sup> <sup>D</sup> xt � <sup>b</sup>xt∣<sup>t</sup> � � follows by an application of the covariance propagation law, together with C xt � <sup>b</sup>xt∣t�<sup>1</sup>; vt � � <sup>¼</sup> Pt∣t�<sup>1</sup>A<sup>T</sup> <sup>t</sup> . Application of the covariance propagation law to (10) gives the variance matrix of vt as follows

$$Q\_{v\_l v\_l} = A\_t P\_{t|t-1} A\_t^T + R\_t \tag{13}$$

since C xt � <sup>b</sup>xt∣t�<sup>1</sup>; <sup>ε</sup><sup>t</sup> � � <sup>¼</sup> 0.

with the initial state-vector x0. The transition matrix from epoch s to t is denoted as Φt,s. Thus

The measurement model: The link between the observables yt and the state-vectors xt is assumed

for t, s ¼ 1, 2, …, with At being the known design matrix. Thus the zero-mean measurement noise ε<sup>t</sup> is assumed to be uncorrelated in time and to be uncorrelated with the initial state-

Initialization: As the mean of x<sup>0</sup> is known, the best predictor of x<sup>0</sup> in the absence of data is the

Time update: Let us choose <sup>Φ</sup>t,t�<sup>1</sup>bxt�1∣t�<sup>1</sup> as a candidate for the BLP <sup>b</sup>xt∣t�1. According to Eq. (1), the candidate would be the BLP if it fulfills two conditions: (1) it must be unbiased and (2) it must have a prediction error uncorrelated with the previous data <sup>y</sup>½ � <sup>t</sup>�<sup>1</sup> . The first condition,

<sup>E</sup>ð Þ¼ dt 0. The second condition, i.e. Cðxt � <sup>Φ</sup>t,t�<sup>1</sup>bxt�1∣t�<sup>1</sup>; <sup>y</sup>½ � <sup>t</sup>�<sup>1</sup> Þ ¼ 0, follows from the fact that the prediction error xt � <sup>Φ</sup>t,t�<sup>1</sup>bxt�1∣t�<sup>1</sup> is a function of the previous BLP prediction error

xt � <sup>Φ</sup>t,t�<sup>1</sup>bxt�1∣t�<sup>1</sup> <sup>¼</sup> <sup>Φ</sup>t,t�<sup>1</sup> xt�<sup>1</sup> � <sup>b</sup>xt�1∣t�<sup>1</sup>

<sup>b</sup>xt∣t�<sup>1</sup> <sup>¼</sup> <sup>Φ</sup>t,t�<sup>1</sup>bxt�1∣t�<sup>1</sup>, Pt∣t�<sup>1</sup> <sup>¼</sup> <sup>Φ</sup>t,t�<sup>1</sup> Pt�1∣t�<sup>1</sup> <sup>Φ</sup><sup>T</sup>

� � <sup>¼</sup> 0.

that are both uncorrelated with the previous data <sup>y</sup>½ � <sup>t</sup>�<sup>1</sup> . Hence, the time update is given by

mean Eð Þ¼ x<sup>0</sup> x<sup>0</sup>∣0. Hence, the initialization is simply given by

� � <sup>¼</sup> <sup>D</sup>ð Þ <sup>x</sup><sup>0</sup> .

i.e. EðΦt,t�<sup>1</sup>bxt�1∣t�<sup>1</sup>Þ ¼ <sup>E</sup>ð Þ xt , follows from Eq. (2) and the equalities E <sup>b</sup>xt�1∣t�<sup>1</sup>

That the initial error variance matrix <sup>P</sup><sup>0</sup>∣<sup>0</sup> <sup>¼</sup> <sup>D</sup> <sup>x</sup><sup>0</sup> � <sup>b</sup>x<sup>0</sup>∣<sup>0</sup>

xt�<sup>1</sup> � <sup>b</sup>xt�1∣t�<sup>1</sup> and the system noise dt, i.e. (cf. 2)

The error variance matrix Pt∣t�<sup>1</sup> <sup>¼</sup> <sup>D</sup> xt � <sup>b</sup>xt∣t�<sup>1</sup>

tion law to (8), together with C xt�<sup>1</sup> � <sup>b</sup>xt�1∣t�<sup>1</sup>; dt

yt ¼ Atxt þ εt, t ¼ 1, 2, … (5)

<sup>b</sup>x<sup>0</sup>∣<sup>0</sup> <sup>¼</sup> <sup>x</sup><sup>0</sup>∣<sup>0</sup>, P<sup>0</sup>∣<sup>0</sup> <sup>¼</sup> Qx0x<sup>0</sup> (7)

� � is identical to the variance matrix Qx0x<sup>0</sup>

� � <sup>þ</sup> dt, (8)

� � follows by applying the covariance propaga-

� � <sup>¼</sup> <sup>E</sup>ð Þ xt�<sup>1</sup> and

t,t�<sup>1</sup> <sup>þ</sup> St (9)

Eð Þ¼ ε<sup>t</sup> 0, C ε<sup>t</sup> ð Þ¼ ; ε<sup>s</sup> Rt δt,s, C ε<sup>t</sup> ð Þ¼ ; x<sup>0</sup> 0, C ε<sup>t</sup> ð Þ¼ ; ds 0 (6)

Φ�<sup>1</sup>

given as

with

t,s ¼ Φs,t and Φt,t ¼ I (the identity matrix).

276 Kalman Filters - Theory for Advanced Applications

vector x<sup>0</sup> and the system noise dt.

follows from the equality D x<sup>0</sup> � x<sup>0</sup>∣<sup>0</sup>

2.2. The three-step recursion

#### 2.3. A remark on the filter initialization

In the derivation of the Kalman filter one assumes the mean of the random initial state-vector x0, in Eq. (3), to be known, see e.g. [30–37]. This is because of the BLP structure (1) that needs knowledge of the means Eð Þ xt and Eðy½ � <sup>τ</sup> <sup>Þ</sup>. Since in many, if not most, applications the means of the state-vectors x1, …, xt are unknown, such derivation is therefore not appropriate. As shown in Ref. [38], one can do away with this need to have both the initial mean x<sup>0</sup>∣<sup>0</sup> and variance matrix Qx0x<sup>0</sup> , given in Eq. (3), known. The corresponding three-step recursion would then follow the Best Linear Unbiased Prediction (BLUP) principle and not that of the BLP. The BLUP is also a MMSE predictor, but within a more restrictive class of predictors. It replaces the means Eð Þ xt and Eðy½ � <sup>τ</sup> <sup>Þ</sup> by their corresponding Best Linear Unbiased Estimators (BLUEs). Within such BLUP recursion, the initialization Eq. (7) is revised and takes place at time instance t ¼ 1 in the presence of the data y1. Provided that matrix A<sup>1</sup> is of full column rank, the predictor <sup>b</sup>x<sup>1</sup>∣<sup>1</sup> follows from solving the normal equations

$$N\_1 \hat{\mathbf{x}}\_{1|1} = r\_1 \quad \text{with} \quad N\_1 = A\_1^T R\_1^{-1} A\_1 \quad r\_1 = A\_1^T R\_1^{-1} y\_1 \tag{14}$$

Thus

$$
\widehat{\mathbf{x}}\_{1|1} = \mathbf{N}\_1^{-1} \mathbf{r}\_1 \quad \text{and} \quad P\_{1|1} = \mathbf{N}\_1^{-1} \tag{15}
$$

information vector it∣<sup>τ</sup> : <sup>¼</sup> <sup>P</sup>�<sup>1</sup>

Time-update : <sup>ð</sup>Φt,t�<sup>1</sup> Pt�1∣t�<sup>1</sup> <sup>Φ</sup><sup>T</sup>

Measurement-update : <sup>ð</sup>Pt∣t�<sup>1</sup> � Pt∣t�<sup>1</sup>A<sup>T</sup>

<sup>t</sup>�1∣t�<sup>1</sup>Φ<sup>t</sup>�1,t.

<sup>t</sup>�1,t <sup>P</sup>�<sup>1</sup>

matrix-inversion equalities

where Mt <sup>¼</sup> <sup>Φ</sup><sup>T</sup>

vector it∣<sup>t</sup> and matrix I<sup>t</sup>∣<sup>t</sup>.

yt

<sup>t</sup>∣<sup>τ</sup> <sup>b</sup>xt∣<sup>τ</sup> and information matrix <sup>I</sup><sup>t</sup>∣<sup>τ</sup> : <sup>¼</sup> <sup>P</sup>�<sup>1</sup>

�<sup>1</sup> <sup>¼</sup> <sup>P</sup>�<sup>1</sup>

�<sup>1</sup> <sup>¼</sup> Mt � MtðMt <sup>þ</sup> <sup>S</sup>�<sup>1</sup>

<sup>t</sup>∣t�<sup>1</sup> <sup>þ</sup> <sup>A</sup><sup>T</sup>

<sup>t</sup> R�<sup>1</sup> <sup>t</sup> At

Consensus-Based Distributed Filtering for GNSS http://dx.doi.org/10.5772/intechopen.71138

Given the definition above, the information filter recursion concerning the time-evolution of it∣<sup>t</sup> and I<sup>t</sup>∣<sup>t</sup> would then follow from the recursion Eqs. (15), (9) and (12), along with the following

t,t�<sup>1</sup> <sup>þ</sup> St<sup>Þ</sup>

The algorithmic steps of the information filter are presented in Figure 1. In the absence of data, the filter is initialized by the zero information i<sup>1</sup>∣<sup>0</sup> ¼ 0 and I<sup>1</sup>∣<sup>0</sup> ¼ 0. In the presence of the data

, the corresponding normal matrix Nt and right-hand-side vector rt are added to the time update information it∣t�<sup>1</sup> and I<sup>t</sup>∣t�<sup>1</sup> to obtain the measurement update information it∣<sup>t</sup> and I<sup>t</sup>∣<sup>t</sup>.

Figure 1. Algorithmic steps of the information filter recursion concerning the time-evolution of the information

AtPt∣t�<sup>1</sup>Þ

<sup>t</sup> Qvtvt

<sup>t</sup>∣<sup>τ</sup> (17)

279

(18)

<sup>t</sup> Þ �1 Mt

The above error variance matrix P<sup>1</sup>∣<sup>1</sup> is thus not dependent on the variance matrix of x1, i.e. Qx1x<sup>1</sup> <sup>¼</sup> <sup>Φ</sup>1, <sup>0</sup> Qx0x<sup>0</sup> <sup>Φ</sup><sup>T</sup> <sup>1</sup>,<sup>0</sup> þ S1. This is, however, not the case with the variance matrix of the predictor <sup>b</sup>x<sup>1</sup>∣<sup>1</sup> itself, i.e. <sup>Q</sup>bx<sup>1</sup>∣<sup>1</sup>bx<sup>1</sup>∣<sup>1</sup> <sup>¼</sup> <sup>D</sup> <sup>b</sup>x<sup>1</sup>∣<sup>1</sup> � �. This variance matrix is given by [38]

$$Q\_{\widehat{\mathcal{X}}\_{1|1}\widehat{\mathcal{X}}\_{1|1}} = Q\_{x\_1x\_1} + P\_{1|1} \tag{16}$$

showing that <sup>P</sup><sup>1</sup>∣<sup>1</sup> 6¼ <sup>Q</sup>bx<sup>1</sup>∣<sup>1</sup>bx<sup>1</sup>∣<sup>1</sup> . Matrices Pt∣<sup>t</sup> and <sup>Q</sup>bxt∣<sup>t</sup>bxt∣<sup>t</sup> (t ¼ 1, 2, …) are used for two different purposes. The error variance matrix Pt∣<sup>t</sup> <sup>¼</sup> <sup>D</sup> xt � <sup>b</sup>xt∣<sup>t</sup> � � is a measure of 'closeness' of <sup>b</sup>xt∣<sup>t</sup> to its target random vector xt, thereby meant to describe the 'quality' of prediction, i.e. precision of the prediction error xt � <sup>b</sup>xt∣<sup>t</sup> � �. The variance matrix <sup>Q</sup>bxt∣<sup>t</sup>bxt∣<sup>t</sup> <sup>¼</sup> <sup>D</sup> <sup>b</sup>xt∣<sup>t</sup> � � however, is a measure of closeness of <sup>b</sup>xt∣<sup>t</sup> to the nonrandom vector Eð Þ xt , as D <sup>b</sup>xt∣<sup>t</sup> � � <sup>¼</sup> <sup>D</sup>ðEð Þ� xt <sup>b</sup>xt∣<sup>t</sup>Þ. Thus <sup>Q</sup>bxt∣<sup>t</sup>bxt∣<sup>t</sup> does not describe the quality of prediction, but instead the precision of the predictor <sup>b</sup>xt∣<sup>t</sup>.

The MMSE of the BLUP recursion is never smaller than that of the Kalman filter, as the Kalman filter makes use of additional information, namely, the known mean x<sup>0</sup>∣<sup>0</sup> and variance matrix Qx0x<sup>0</sup> . When the stated information is available, the BLUP recursion is shown to encompass the Kalman filter as a special case [39]. In the following we therefore assume that the means of the state-vectors x1, …, xt are unknown, a situation that often applies to GNSS applications.

#### 2.4. Filtering in information form

The three-step recursion presented in Eqs. (7), (9) and (12) concerns the time-evolution of the predictor <sup>b</sup>xt∣<sup>t</sup> and the error variance matrix Pt∣<sup>t</sup>. As shown in Eq. (15), both <sup>P</sup><sup>1</sup>∣<sup>1</sup> and <sup>b</sup>x<sup>1</sup>∣<sup>1</sup> can be determined by the normal matrix <sup>N</sup><sup>1</sup> <sup>¼</sup> <sup>P</sup>�<sup>1</sup> <sup>1</sup>∣<sup>1</sup> and the right-hand-side vector <sup>r</sup><sup>1</sup> <sup>¼</sup> <sup>P</sup>�<sup>1</sup> <sup>1</sup>∣<sup>1</sup>bx<sup>1</sup>∣1. One can therefore alternatively develop recursion concerning the time-evolution of P�<sup>1</sup> <sup>t</sup>∣<sup>t</sup> and P�<sup>1</sup> <sup>t</sup>∣<sup>t</sup> <sup>b</sup>xt∣<sup>t</sup>. From a computational point of view, such recursion is found to be very suitable when the inverse-variance or information matrices S�<sup>1</sup> <sup>t</sup> and R�<sup>1</sup> <sup>t</sup> serve as input rather than the variance matrices St and Rt. To that end, one may define [34]

information vector it∣<sup>τ</sup> : <sup>¼</sup> <sup>P</sup>�<sup>1</sup> <sup>t</sup>∣<sup>τ</sup> <sup>b</sup>xt∣<sup>τ</sup> and information matrix <sup>I</sup><sup>t</sup>∣<sup>τ</sup> : <sup>¼</sup> <sup>P</sup>�<sup>1</sup> <sup>t</sup>∣<sup>τ</sup> (17)

Given the definition above, the information filter recursion concerning the time-evolution of it∣<sup>t</sup> and I<sup>t</sup>∣<sup>t</sup> would then follow from the recursion Eqs. (15), (9) and (12), along with the following matrix-inversion equalities

$$\begin{array}{llll}\text{Time-update}: & \left(\boldsymbol{\Phi}\_{t,t-1}\boldsymbol{P}\_{t-1|t-1}\boldsymbol{\Phi}\_{t,t-1}^{T} + \boldsymbol{S}\_{t}\right)^{-1} &=& M\_{t} - M\_{t}(M\_{t} + \boldsymbol{S}\_{t}^{-1})^{-1}M\_{t} \\\\ \text{Measurement-update}: & \left(\boldsymbol{P}\_{t|t-1} - \boldsymbol{P}\_{t|t-1}\boldsymbol{A}\_{t}^{T}\boldsymbol{Q}\_{t,t\_{t}}A\_{t}\boldsymbol{P}\_{t|t-1}\right)^{-1} &=& P\_{t|t-1}^{-1} + A\_{t}^{T}R\_{t}^{-1}A\_{t} \end{array} \tag{18}$$

where Mt <sup>¼</sup> <sup>Φ</sup><sup>T</sup> <sup>t</sup>�1,t <sup>P</sup>�<sup>1</sup> <sup>t</sup>�1∣t�<sup>1</sup>Φ<sup>t</sup>�1,t.

matrix Qx0x<sup>0</sup> , given in Eq. (3), known. The corresponding three-step recursion would then follow the Best Linear Unbiased Prediction (BLUP) principle and not that of the BLP. The BLUP is also a MMSE predictor, but within a more restrictive class of predictors. It replaces the means Eð Þ xt and Eðy½ � <sup>τ</sup> <sup>Þ</sup> by their corresponding Best Linear Unbiased Estimators (BLUEs). Within such BLUP recursion, the initialization Eq. (7) is revised and takes place at time instance t ¼ 1 in the presence of the data y1. Provided that matrix A<sup>1</sup> is of full column rank,

1R�<sup>1</sup>

<sup>1</sup> <sup>r</sup>1, and <sup>P</sup><sup>1</sup>∣<sup>1</sup> <sup>¼</sup> <sup>N</sup>�<sup>1</sup>

<sup>1</sup>,<sup>0</sup> þ S1. This is, however, not the case with the variance matrix of the

� �. This variance matrix is given by [38]

The above error variance matrix P<sup>1</sup>∣<sup>1</sup> is thus not dependent on the variance matrix of x1, i.e.

target random vector xt, thereby meant to describe the 'quality' of prediction, i.e. precision of

The MMSE of the BLUP recursion is never smaller than that of the Kalman filter, as the Kalman filter makes use of additional information, namely, the known mean x<sup>0</sup>∣<sup>0</sup> and variance matrix Qx0x<sup>0</sup> . When the stated information is available, the BLUP recursion is shown to encompass the Kalman filter as a special case [39]. In the following we therefore assume that the means of the state-vectors x1, …, xt are unknown, a situation that often applies to GNSS applications.

The three-step recursion presented in Eqs. (7), (9) and (12) concerns the time-evolution of the predictor <sup>b</sup>xt∣<sup>t</sup> and the error variance matrix Pt∣<sup>t</sup>. As shown in Eq. (15), both <sup>P</sup><sup>1</sup>∣<sup>1</sup> and <sup>b</sup>x<sup>1</sup>∣<sup>1</sup> can be

From a computational point of view, such recursion is found to be very suitable when the

<sup>t</sup> and R�<sup>1</sup>

can therefore alternatively develop recursion concerning the time-evolution of P�<sup>1</sup>

<sup>1</sup> <sup>A</sup>1, r<sup>1</sup> <sup>¼</sup> <sup>A</sup><sup>T</sup>

1R�<sup>1</sup>

¼ Qx1x<sup>1</sup> þ P<sup>1</sup>∣<sup>1</sup> (16)

� � is a measure of 'closeness' of <sup>b</sup>xt∣<sup>t</sup> to its

� � <sup>¼</sup> <sup>D</sup>ðEð Þ� xt <sup>b</sup>xt∣<sup>t</sup>Þ. Thus <sup>Q</sup>bxt∣<sup>t</sup>bxt∣<sup>t</sup>

<sup>¼</sup> <sup>D</sup> <sup>b</sup>xt∣<sup>t</sup>

<sup>1</sup>∣<sup>1</sup> and the right-hand-side vector <sup>r</sup><sup>1</sup> <sup>¼</sup> <sup>P</sup>�<sup>1</sup>

<sup>t</sup> serve as input rather than the variance

(t ¼ 1, 2, …) are used for two different

� � however, is a measure of

does

<sup>1</sup>∣<sup>1</sup>bx<sup>1</sup>∣1. One

<sup>t</sup>∣<sup>t</sup> <sup>b</sup>xt∣<sup>t</sup>.

<sup>t</sup>∣<sup>t</sup> and P�<sup>1</sup>

<sup>1</sup> y<sup>1</sup> (14)

<sup>1</sup> (15)

the predictor <sup>b</sup>x<sup>1</sup>∣<sup>1</sup> follows from solving the normal equations

Thus

Qx1x<sup>1</sup> <sup>¼</sup> <sup>Φ</sup>1, <sup>0</sup> Qx0x<sup>0</sup> <sup>Φ</sup><sup>T</sup>

predictor <sup>b</sup>x<sup>1</sup>∣<sup>1</sup> itself, i.e. <sup>Q</sup>bx<sup>1</sup>∣<sup>1</sup>bx<sup>1</sup>∣<sup>1</sup>

278 Kalman Filters - Theory for Advanced Applications

showing that <sup>P</sup><sup>1</sup>∣<sup>1</sup> 6¼ <sup>Q</sup>bx<sup>1</sup>∣<sup>1</sup>bx<sup>1</sup>∣<sup>1</sup>

the prediction error xt � <sup>b</sup>xt∣<sup>t</sup>

2.4. Filtering in information form

determined by the normal matrix <sup>N</sup><sup>1</sup> <sup>¼</sup> <sup>P</sup>�<sup>1</sup>

inverse-variance or information matrices S�<sup>1</sup>

matrices St and Rt. To that end, one may define [34]

<sup>N</sup><sup>1</sup>bx<sup>1</sup>∣<sup>1</sup> <sup>¼</sup> <sup>r</sup>1, with <sup>N</sup><sup>1</sup> <sup>¼</sup> <sup>A</sup><sup>T</sup>

<sup>b</sup>x<sup>1</sup>∣<sup>1</sup> <sup>¼</sup> <sup>N</sup>�<sup>1</sup>

<sup>¼</sup> <sup>D</sup> <sup>b</sup>x<sup>1</sup>∣<sup>1</sup>

purposes. The error variance matrix Pt∣<sup>t</sup> <sup>¼</sup> <sup>D</sup> xt � <sup>b</sup>xt∣<sup>t</sup>

closeness of <sup>b</sup>xt∣<sup>t</sup> to the nonrandom vector Eð Þ xt , as D <sup>b</sup>xt∣<sup>t</sup>

<sup>Q</sup>bx<sup>1</sup>∣<sup>1</sup>bx<sup>1</sup>∣<sup>1</sup>

� �. The variance matrix <sup>Q</sup>bxt∣<sup>t</sup>bxt∣<sup>t</sup>

not describe the quality of prediction, but instead the precision of the predictor <sup>b</sup>xt∣<sup>t</sup>.

. Matrices Pt∣<sup>t</sup> and <sup>Q</sup>bxt∣<sup>t</sup>bxt∣<sup>t</sup>

The algorithmic steps of the information filter are presented in Figure 1. In the absence of data, the filter is initialized by the zero information i<sup>1</sup>∣<sup>0</sup> ¼ 0 and I<sup>1</sup>∣<sup>0</sup> ¼ 0. In the presence of the data yt , the corresponding normal matrix Nt and right-hand-side vector rt are added to the time update information it∣t�<sup>1</sup> and I<sup>t</sup>∣t�<sup>1</sup> to obtain the measurement update information it∣<sup>t</sup> and I<sup>t</sup>∣<sup>t</sup>.


Figure 1. Algorithmic steps of the information filter recursion concerning the time-evolution of the information vector it∣<sup>t</sup> and matrix I<sup>t</sup>∣<sup>t</sup>.

The transition matrix <sup>Φ</sup>t,t�<sup>1</sup> and inverse-variance matrix <sup>S</sup>�<sup>1</sup> <sup>t</sup> would then be used to time update the previous information it�1∣t�<sup>1</sup> and I<sup>t</sup>�1∣t�1.

Singular matrix St: In the first expression of Eq. (18) one assumes the variance matrix St to be nonsingular and invertible. There are, however, situations where some of the elements of the state-vector xt are nonrandom, i.e., the corresponding system noise is identically zero. As a consequence, the variance matrix St becomes singular and the inverse-matrix S�<sup>1</sup> <sup>t</sup> does not exist. An example of such concerns the presence of the GNSS carrier-phase ambiguities in the filter state-vector which are treated constant in time. In such cases the information time update in Figure 1 must be generalized so as to accommodate singular variance matrices St. Let S~<sup>t</sup> be an invertible sub-matrix of St that has the same rank as that of St. Then there exists a fullcolumn rank matrix Ht such that

$$\mathcal{S}\_t = H\_t \tilde{\mathcal{S}}\_t H\_t^T \tag{19}$$

yt ¼ Atxt þ ε<sup>t</sup> )

independently from one another. This yields

C εi,t; εj,t

vector rt <sup>¼</sup> AT

where

<sup>t</sup> R�<sup>1</sup>

with the average quantities

y1,t ⋮ yi,t ⋮ yn,t

network of n sensor nodes where each collects its own observable vector yi,t

With the extra assumption Eq. (23), the normal matrix Nt <sup>¼</sup> <sup>A</sup><sup>T</sup>

<sup>t</sup> yt can then be, respectively, expressed as

i,tR�<sup>1</sup>

i¼1

Nt <sup>¼</sup> <sup>X</sup><sup>n</sup> i¼1

Ni,t <sup>¼</sup> <sup>A</sup><sup>T</sup>

ually added to the information states I<sup>t</sup>∣t�<sup>1</sup> and it∣t�1, that is

<sup>I</sup><sup>t</sup>∣<sup>t</sup> <sup>¼</sup> <sup>I</sup><sup>t</sup>∣t�<sup>1</sup> <sup>þ</sup>X<sup>n</sup>

Nt <sup>¼</sup> <sup>1</sup> n Xn i¼1

¼

Accordingly, the observable vector yt is partitioned into n sub-vectors yi,t (i ¼ 1,…, n), each having its own design matrix Ai,t and measurement noise vector εi,t. One can think of a

determine a common state-vector xt. Let us further assume that the nodes collect observables

Thus the measurement noise vectors εi,t (i ¼ 1, …, n) are assumed to be mutually uncorrelated.

Ni,t, and rt <sup>¼</sup> <sup>X</sup><sup>n</sup>

i,t Ai,t, and ri,t <sup>¼</sup> AT

Ni,t, it∣<sup>t</sup> <sup>¼</sup> it∣t�<sup>1</sup> <sup>þ</sup>X<sup>n</sup>

According to Eq. (24), the measurement information of each node, say Ni,t and ri,t, is individ-

Now consider the situation where each node runs its own local information filter, thus having its own information states Ii,t∣<sup>t</sup> and ii,t∣<sup>t</sup> (i ¼ 1, …, n). The task is to recursively update the local states Ii,t∣<sup>t</sup> and ii,t∣<sup>t</sup> in a way that they remain equal to their central counterparts I<sup>t</sup>∣<sup>t</sup> and it∣<sup>t</sup> given in Eq. (26). Suppose that such equalities hold at the time update, i.e. Ii,t∣t�<sup>1</sup> ¼ I<sup>t</sup>∣t�<sup>1</sup> and ii,t∣t�<sup>1</sup> ¼ it∣t�1. Given the number of contributing nodes n, each node just needs to be provided

Ni,t, and rt <sup>¼</sup> <sup>1</sup>

n Xn i¼1

A1,t ⋮ Ai,t ⋮ An,t

xt þ

� � <sup>¼</sup> Ri,t <sup>δ</sup>i,j, for i, j <sup>¼</sup> <sup>1</sup>,…, n, and <sup>t</sup> <sup>¼</sup> <sup>1</sup>, <sup>2</sup>, … (23)

i¼1

i,t R�<sup>1</sup>

i¼1

<sup>t</sup> R�<sup>1</sup>

ε1,t ⋮ εi,t ⋮ εn,t

, t ¼ 1, 2, … (22)

Consensus-Based Distributed Filtering for GNSS http://dx.doi.org/10.5772/intechopen.71138

, but aiming to

281

<sup>t</sup> At and right-hand-side

ri,t (24)

i,t yi,t (25)

ri,t (26)

ri,t (27)

Matrix Ht can be, for instance, structured by the columns of the identity matrix I corresponding to the columns of St on which the sub-matrix S~<sup>t</sup> is positioned. The special case

$$\mathbf{S}\_{t} = \begin{bmatrix} \tilde{\mathbf{S}}\_{t} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} \end{bmatrix} = \begin{bmatrix} I \\ \mathbf{0} \end{bmatrix} \tilde{\mathbf{S}}\_{t} \begin{bmatrix} I \\ \mathbf{0} \end{bmatrix}^{T} \quad \Rightarrow \quad H\_{t} = \begin{bmatrix} I \\ \mathbf{0} \end{bmatrix} \tag{20}$$

shows an example of the representation (19). With Eq. (19), a generalization of the time update (Figure 1) can be shown to be given by

$$\mathcal{Z}\_{t|t-1} = M\_t - M\_t H\_t \left[ H\_t^T M\_t H\_t + \tilde{\mathcal{S}}\_t^{-1} \right]^{-1} H\_t^T M\_t \tag{21}$$

Thus instead of S�<sup>1</sup> <sup>t</sup> , the inverse-matrix <sup>S</sup>~�<sup>1</sup> <sup>t</sup> and Ht are assumed available.

#### 2.5. Additivity property of the information measurement update

As stated previously, the information filter delivers outcomes equivalent to those of the Kalman filter recursion. Thus any particular preference for the information filter must be attributed to the computational effort required for obtaining the outcomes. For instance, if handling matrix inversion requires low computational complexity when working with the input inverse-matrices S�<sup>1</sup> t and R�<sup>1</sup> <sup>t</sup> , the information filter appears to be more suitable. In this subsection we will highlight yet another property of the information filter that makes this recursion particularly useful for distributed processing.

As shown in Figure 1, the information measurement update is additive in the sense that the measurement information Nt and rt is added to the information states I<sup>t</sup>∣t�<sup>1</sup> and it∣t�1. We now make a start to show how such additivity property lends itself to distributed filtering. Let the measurement model Eq. (5) be partitioned as

#### Consensus-Based Distributed Filtering for GNSS http://dx.doi.org/10.5772/intechopen.71138 281

$$\mathbf{y}\_t = A\_t \mathbf{x}\_t + \boldsymbol{\varepsilon}\_t \quad \Rightarrow \quad \begin{bmatrix} \mathbf{y}\_{1,t} \\ \vdots \\ \mathbf{y}\_{i,t} \\ \vdots \\ \vdots \\ \mathbf{y}\_{n,t} \end{bmatrix} = \begin{bmatrix} A\_{1,t} \\ \vdots \\ A\_{i,t} \\ \vdots \\ \vdots \\ A\_{n,t} \end{bmatrix} \mathbf{x}\_t + \begin{bmatrix} \boldsymbol{\varepsilon}\_{1,t} \\ \vdots \\ \boldsymbol{\varepsilon}\_{i,t} \\ \vdots \\ \boldsymbol{\varepsilon}\_{n,t} \end{bmatrix}, \quad t = 1, 2, \dots \tag{22}$$

Accordingly, the observable vector yt is partitioned into n sub-vectors yi,t (i ¼ 1,…, n), each having its own design matrix Ai,t and measurement noise vector εi,t. One can think of a network of n sensor nodes where each collects its own observable vector yi,t , but aiming to determine a common state-vector xt. Let us further assume that the nodes collect observables independently from one another. This yields

$$\mathbb{C}\{\varepsilon\_{i,t},\varepsilon\_{j,t}\} = \mathbb{R}\_{i,t}\delta\_{i,j} \quad \text{for} \quad i,j = 1,\ldots,n, \quad \text{and} \quad t = 1,2,\ldots \tag{23}$$

Thus the measurement noise vectors εi,t (i ¼ 1, …, n) are assumed to be mutually uncorrelated. With the extra assumption Eq. (23), the normal matrix Nt <sup>¼</sup> <sup>A</sup><sup>T</sup> <sup>t</sup> R�<sup>1</sup> <sup>t</sup> At and right-hand-side vector rt <sup>¼</sup> AT <sup>t</sup> R�<sup>1</sup> <sup>t</sup> yt can then be, respectively, expressed as

$$N\_t = \sum\_{i=1}^{n} N\_{i,t} \quad \text{and} \quad r\_t = \sum\_{i=1}^{n} r\_{i,t} \tag{24}$$

where

The transition matrix <sup>Φ</sup>t,t�<sup>1</sup> and inverse-variance matrix <sup>S</sup>�<sup>1</sup>

St <sup>¼</sup> <sup>S</sup>~<sup>t</sup> <sup>0</sup> 0 0 " #

<sup>t</sup> , the inverse-matrix <sup>S</sup>~�<sup>1</sup>

the measurement model Eq. (5) be partitioned as

2.5. Additivity property of the information measurement update

Singular matrix St: In the first expression of Eq. (18) one assumes the variance matrix St to be nonsingular and invertible. There are, however, situations where some of the elements of the state-vector xt are nonrandom, i.e., the corresponding system noise is identically zero. As a

exist. An example of such concerns the presence of the GNSS carrier-phase ambiguities in the filter state-vector which are treated constant in time. In such cases the information time update in Figure 1 must be generalized so as to accommodate singular variance matrices St. Let S~<sup>t</sup> be an invertible sub-matrix of St that has the same rank as that of St. Then there exists a full-

St <sup>¼</sup> Ht <sup>S</sup>~<sup>t</sup> <sup>H</sup><sup>T</sup>

Matrix Ht can be, for instance, structured by the columns of the identity matrix I corresponding to the columns of St on which the sub-matrix S~<sup>t</sup> is positioned. The special case

shows an example of the representation (19). With Eq. (19), a generalization of the time update

As stated previously, the information filter delivers outcomes equivalent to those of the Kalman filter recursion. Thus any particular preference for the information filter must be attributed to the computational effort required for obtaining the outcomes. For instance, if handling matrix inversion requires low computational complexity when working with the input inverse-matrices S�<sup>1</sup>

As shown in Figure 1, the information measurement update is additive in the sense that the measurement information Nt and rt is added to the information states I<sup>t</sup>∣t�<sup>1</sup> and it∣t�1. We now make a start to show how such additivity property lends itself to distributed filtering. Let

<sup>t</sup> , the information filter appears to be more suitable. In this subsection we will highlight yet another property of the information filter that makes this recursion particularly useful for

<sup>t</sup> MtHt <sup>þ</sup> <sup>S</sup>~�<sup>1</sup> t

<sup>t</sup> and Ht are assumed available.

h i�<sup>1</sup>

<sup>¼</sup> <sup>I</sup> 0 � � S~t I 0 � �<sup>T</sup>

<sup>I</sup><sup>t</sup>∣t�<sup>1</sup> <sup>¼</sup> Mt � MtHt <sup>H</sup><sup>T</sup>

consequence, the variance matrix St becomes singular and the inverse-matrix S�<sup>1</sup>

update the previous information it�1∣t�<sup>1</sup> and I<sup>t</sup>�1∣t�1.

column rank matrix Ht such that

280 Kalman Filters - Theory for Advanced Applications

(Figure 1) can be shown to be given by

Thus instead of S�<sup>1</sup>

distributed processing.

and R�<sup>1</sup>

<sup>t</sup> would then be used to time

<sup>t</sup> (19)

, (20)

<sup>t</sup> Mt (21)

t

) Ht <sup>¼</sup> <sup>I</sup>

0 � �

H<sup>T</sup>

<sup>t</sup> does not

$$N\_{i,t} = A\_{i,t}^T R\_{i,t}^{-1} A\_{i,t} \quad \text{and} \quad r\_{i,t} = A\_{i,t}^T R\_{i,t}^{-1} y\_{i,t} \tag{25}$$

According to Eq. (24), the measurement information of each node, say Ni,t and ri,t, is individually added to the information states I<sup>t</sup>∣t�<sup>1</sup> and it∣t�1, that is

$$\mathcal{Z}\_{t|t} = \mathcal{Z}\_{t|t-1} + \sum\_{i=1}^{n} \mathcal{N}\_{i,t} \qquad \mathcal{i}\_{t|t} = \mathcal{i}\_{t|t-1} + \sum\_{i=1}^{n} r\_{i,t} \tag{26}$$

Now consider the situation where each node runs its own local information filter, thus having its own information states Ii,t∣<sup>t</sup> and ii,t∣<sup>t</sup> (i ¼ 1, …, n). The task is to recursively update the local states Ii,t∣<sup>t</sup> and ii,t∣<sup>t</sup> in a way that they remain equal to their central counterparts I<sup>t</sup>∣<sup>t</sup> and it∣<sup>t</sup> given in Eq. (26). Suppose that such equalities hold at the time update, i.e. Ii,t∣t�<sup>1</sup> ¼ I<sup>t</sup>∣t�<sup>1</sup> and ii,t∣t�<sup>1</sup> ¼ it∣t�1. Given the number of contributing nodes n, each node just needs to be provided with the average quantities

$$\overline{N\_t} = \frac{1}{n} \sum\_{i=1}^{n} N\_{i,t} \quad \text{and} \quad \overline{r}\_t = \frac{1}{n} \sum\_{i=1}^{n} r\_{i,t} \tag{27}$$

The local states Ii,t∣t�<sup>1</sup> and ii,t∣t�<sup>1</sup> would then be measurement updated as (cf. 26)

$$\mathcal{L}\_{i,t|t} = \mathcal{L}\_{i,t|t-1} + nN\_{t\prime} \qquad \qquad i\_{i,t|t} = \quad i\_{i,t|t-1} + n\overline{r}\_t \tag{28}$$

links though, i.e. the edge set E<sup>k</sup> depends on k. As in Figure 2 (b), the number of links between the nodes can be different for different sessions k ¼ 1,…, ko. Each session represents a graph that may not be connected. In a 'connected' graph, every vertex is linked to all other vertices at least through one path. In order for information to flow from each node to all other nodes, the

Figure 2. Communications graphs of 20 sensor nodes. The edges represent two-way communication links between the nodes. (a) Network with 49 links. (b) Network with different numbers of links over four sessions: 7 links in session 1 (R), 6

> ko k¼1

is therefore assumed to be connected. We define the neighbors of node i as those to which the node i has direct links. For every session k, they are collected in the set N i, <sup>k</sup> ¼ f g jj ð Þ j; i ∈E<sup>k</sup> . For instance for network (a) of Figure 2, we have only one session, i.e. ko ¼ 1, in which N <sup>2</sup>,<sup>1</sup> ¼ f g 1; 3; 4; 5 represents the neighbors of node 2. In case of network (b) however, we have different links over four sessions, i.e. ko ¼ 4. In this case, the neighbors of node 2 are given by four sets: N <sup>2</sup>,<sup>1</sup> ¼ f g5 in session 1 (red), N <sup>2</sup>, <sup>2</sup> ¼ fg in session 2 (yellow), N <sup>2</sup>, <sup>3</sup> ¼ f g4 in session 3

Given the right-hand-vector ri,t, suppose that node i aims to obtain the average rt for which all other vectors rj,t (∀j 6¼ i) are required to be available (cf. 27). But the node i only has access to those of its neighbors, i.e. the vectors rj,t (j∈ N i, <sup>k</sup>). For the first session k ¼ 1, it would then

<sup>j</sup> <sup>∈</sup>f g <sup>i</sup>;<sup>N</sup> i,<sup>1</sup>

as an approximation of rt, where the scalars wijð Þ1 (j ∈f g i; N i, <sup>1</sup> ) denote the corresponding weights at session k ¼ 1. Now assume that all other nodes j 6¼ i agree to apply the fusion rule

seem to be reasonable to compute a weighted-average of the available vectors, i.e.

ri,tð Þ¼ <sup>1</sup> <sup>X</sup>

E<sup>k</sup> ð Þ ko : a finite number (29)

Consensus-Based Distributed Filtering for GNSS http://dx.doi.org/10.5772/intechopen.71138 283

wijð Þ1 rj,t (30)

union of the graphs G<sup>k</sup> (k ¼ 1,…, ko), i.e.

(green) and N <sup>2</sup>, <sup>4</sup> ¼ fg in session 4 (blue).

3.2. Consensus protocols

G ¼ ð Þ V; E with E ¼ ⋃

links in session 2 (Y), 8 links in session 3 (G) and 7 links in session 4 (B).

that are equal to the central states I<sup>t</sup>∣<sup>t</sup> and it∣<sup>t</sup>, respectively. In this way one has multiple distributed local filters i ¼ 1, …, n, where each recursively delivers results identical to those of a central filter.

To compute the average quantities Nt and rt, node i may need to receive all other information Nj,t and rj,t (j 6¼ i). In other words, node i would require direct connections to all other nodes j 6¼ i, a situation that makes data communication and processing power very expensive (particularly for a large number of nodes). In the following cheaper ways of evaluating the averages Nt and rt are discussed.

#### 3. Average consensus rules

In the previous section, we briefly discussed the potential applicability of the information filter as a tool for handling the measurement model Eq. (22) in a distributed manner. With the representation Eq. (28) however, one may be inclined to conclude that such applicability is limited to the case where the nodes i ¼ 1, …, n, have 'direct' communication connections to one another in order to receive/send their measurement information Ni,t and ri,t (i ¼ 1, …, n).

Instead of having direct connections, the idea is now to relax such a stringent requirement by assuming that the nodes are linked to each other at least through a 'path' so that information can flow from each node to all other nodes. It is therefore assumed that each node along the path plays the role of an agent transferring information to other nodes. To reach the averages Nt and rt, the nodes would then agree on specific 'fusion rules' or consensus protocols, see e.g. [6, 8, 40]. Note that each node exchanges information with neighboring nodes (i.e. those to which the node has direct connections) and not the entire nodes. Therefore, a repeated application of the consensus protocols is required to be carried out. The notion is made precise below.

#### 3.1. Communication graphs

The way the nodes interact with each other to transfer information is referred to as the interaction topology between the nodes. The interaction topology is often described by a directed graph whose vertices and edges, respectively, represent the nodes and communication links [4]. The interaction topology may also undergo a finite number of changes over sessions k ¼ 1,…, ko. In case of one-way links, the directions of the edges face toward the receiving nodes (vertices). Here we assume that the communication links between the nodes are two-way, thus having undirected (or bidirectional) graphs. Examples of such representing a network of 20 nodes with their interaction links are shown in Figure 2. Let an undirected graph at session k be denoted by G<sup>k</sup> ¼ ð Þ V; E<sup>k</sup> where V ¼ f g 1; …; n is the vertex set and E<sup>k</sup> ⊂ f g ð Þj i; j i; j ∈V is the edge set. We assume that the nodes remain unchanged over time, that is why the subscript k is omitted for V. This is generally not the case with their interaction

Figure 2. Communications graphs of 20 sensor nodes. The edges represent two-way communication links between the nodes. (a) Network with 49 links. (b) Network with different numbers of links over four sessions: 7 links in session 1 (R), 6 links in session 2 (Y), 8 links in session 3 (G) and 7 links in session 4 (B).

links though, i.e. the edge set E<sup>k</sup> depends on k. As in Figure 2 (b), the number of links between the nodes can be different for different sessions k ¼ 1,…, ko. Each session represents a graph that may not be connected. In a 'connected' graph, every vertex is linked to all other vertices at least through one path. In order for information to flow from each node to all other nodes, the union of the graphs G<sup>k</sup> (k ¼ 1,…, ko), i.e.

$$\mathcal{G} = (\mathcal{V}, \mathcal{E}) \quad \text{with} \quad \mathcal{E} = \bigcup\_{k=1}^{k\_o} \mathcal{E}\_k \quad (k\_o \text{ : a finite number}) \tag{29}$$

is therefore assumed to be connected. We define the neighbors of node i as those to which the node i has direct links. For every session k, they are collected in the set N i, <sup>k</sup> ¼ f g jj ð Þ j; i ∈E<sup>k</sup> . For instance for network (a) of Figure 2, we have only one session, i.e. ko ¼ 1, in which N <sup>2</sup>,<sup>1</sup> ¼ f g 1; 3; 4; 5 represents the neighbors of node 2. In case of network (b) however, we have different links over four sessions, i.e. ko ¼ 4. In this case, the neighbors of node 2 are given by four sets: N <sup>2</sup>,<sup>1</sup> ¼ f g5 in session 1 (red), N <sup>2</sup>, <sup>2</sup> ¼ fg in session 2 (yellow), N <sup>2</sup>, <sup>3</sup> ¼ f g4 in session 3 (green) and N <sup>2</sup>, <sup>4</sup> ¼ fg in session 4 (blue).

#### 3.2. Consensus protocols

The local states Ii,t∣t�<sup>1</sup> and ii,t∣t�<sup>1</sup> would then be measurement updated as (cf. 26)

a central filter.

averages Nt and rt are discussed.

282 Kalman Filters - Theory for Advanced Applications

3. Average consensus rules

3.1. Communication graphs

that are equal to the central states I<sup>t</sup>∣<sup>t</sup> and it∣<sup>t</sup>, respectively. In this way one has multiple distributed local filters i ¼ 1, …, n, where each recursively delivers results identical to those of

To compute the average quantities Nt and rt, node i may need to receive all other information Nj,t and rj,t (j 6¼ i). In other words, node i would require direct connections to all other nodes j 6¼ i, a situation that makes data communication and processing power very expensive (particularly for a large number of nodes). In the following cheaper ways of evaluating the

In the previous section, we briefly discussed the potential applicability of the information filter as a tool for handling the measurement model Eq. (22) in a distributed manner. With the representation Eq. (28) however, one may be inclined to conclude that such applicability is limited to the case where the nodes i ¼ 1, …, n, have 'direct' communication connections to one another in order to receive/send their measurement information Ni,t and ri,t (i ¼ 1, …, n).

Instead of having direct connections, the idea is now to relax such a stringent requirement by assuming that the nodes are linked to each other at least through a 'path' so that information can flow from each node to all other nodes. It is therefore assumed that each node along the path plays the role of an agent transferring information to other nodes. To reach the averages Nt and rt, the nodes would then agree on specific 'fusion rules' or consensus protocols, see e.g. [6, 8, 40]. Note that each node exchanges information with neighboring nodes (i.e. those to which the node has direct connections) and not the entire nodes. Therefore, a repeated application of the consensus protocols is required to be carried out. The notion is made precise below.

The way the nodes interact with each other to transfer information is referred to as the interaction topology between the nodes. The interaction topology is often described by a directed graph whose vertices and edges, respectively, represent the nodes and communication links [4]. The interaction topology may also undergo a finite number of changes over sessions k ¼ 1,…, ko. In case of one-way links, the directions of the edges face toward the receiving nodes (vertices). Here we assume that the communication links between the nodes are two-way, thus having undirected (or bidirectional) graphs. Examples of such representing a network of 20 nodes with their interaction links are shown in Figure 2. Let an undirected graph at session k be denoted by G<sup>k</sup> ¼ ð Þ V; E<sup>k</sup> where V ¼ f g 1; …; n is the vertex set and E<sup>k</sup> ⊂ f g ð Þj i; j i; j ∈V is the edge set. We assume that the nodes remain unchanged over time, that is why the subscript k is omitted for V. This is generally not the case with their interaction

Ii,t∣<sup>t</sup> ¼ Ii,t∣t�<sup>1</sup> þ nNt, ii,t∣<sup>t</sup> ¼ ii,t∣t�<sup>1</sup> þ nrt (28)

Given the right-hand-vector ri,t, suppose that node i aims to obtain the average rt for which all other vectors rj,t (∀j 6¼ i) are required to be available (cf. 27). But the node i only has access to those of its neighbors, i.e. the vectors rj,t (j∈ N i, <sup>k</sup>). For the first session k ¼ 1, it would then seem to be reasonable to compute a weighted-average of the available vectors, i.e.

$$r\_{i,t}(1) = \sum\_{j \in \{i, \mathcal{N}\_{i,1}\}} w\_{i\bar{\jmath}}(1) r\_{j,t} \tag{30}$$

as an approximation of rt, where the scalars wijð Þ1 (j ∈f g i; N i, <sup>1</sup> ) denote the corresponding weights at session k ¼ 1. Now assume that all other nodes j 6¼ i agree to apply the fusion rule Eq. (30) to those of their own neighbors. Thus the neighboring nodes j∈ N i, <sup>k</sup> also have their own weighted-averages rj,tð Þ1 . But they may have access to those to which the node i has no direct links. In other words, the weighted-averages rj,tð Þ1 (j ∈ N i, 1) contain information on the nodes to which the node i has no access. For the next session k ¼ 2, it is therefore reasonable for the node i to repeat the fusion rule Eq. (30), but now over the new vectors rj,tð Þ1 (j∈ f g i; N i, <sup>2</sup> ), aiming to improve on the earlier approximation ri,tð Þ1 . This yields the following iterative computations

$$r\_{i,t}(k) = \sum\_{j \in \{i, N\_{i,k}\}} w\_{ij}(k) r\_{j,t}(k-1), \quad k = 1, 2, \dots \tag{31}$$

matrices W kð Þ (k ¼ 1, …, ko) have bounded nonnegative entries with positive diagonals, i.e.

Examples of such consensus protocols are given in Table 1. As shown, the weights form a symmetric weight matrix W kð Þ, i.e. wjið Þ¼ k wijð Þk . In all protocols presented, self-weights

and 2 belong to the class of 'maximum-degree' weights, while those of Protocol 3 are referred to as 'Metropolis' weights [8]. The weights of Protocols 1 and 3 are driven by the degrees

(a) of Figure 2 we have dg1ð Þ¼ 1 4 as node 1 has 4 neighbors, while dg14ð Þ¼ 1 7 as node 14 has 7 neighbors. Protocol 4 is only applicable to networks like (b) in Figure 2, i.e. when each node has at most one neighbor at a session [4]. In this case, each node exchanges its information to just one neighbor at a session. Thus for two neighboring nodes i and j we have wiið Þ¼ k wjjð Þ¼ k wijð Þ¼ k 0:5, each averaging ri,tð Þ k � 1 and rj,tð Þ k � 1 to obtain ri,tð Þ¼ k rj,tð Þk . To provide insight into the applicability of the protocols given in Table 1, we apply them to the networks of Figure 2. Twenty values (scalars), say ri (i ¼ 1, …, 20), are generated whose average is equal to 5, i.e. r ¼ 5. Each value is assigned to its corresponding node. For network (a), Protocols 1, 2 and 3 are separately applied, whereas Protocol 4 is only applied to network (b). The corresponding results, up to 30 iterations, are presented in Figure 3. As shown, the iterated values rið Þk (i ¼ 1,…, 20) get closer to their average (i.e. r ¼ 5), the more the number

Figure 3 shows that the states ri,tð Þk (i ¼ 1, …, n) converge to their average rt, but with different rates. The convergence rate depends on the initial states ri,tð Þ¼ 0 ri,t and on the consensus protocol employed. From the figure it seems that the convergence rates of Protocols 1 and 3 are about the same, higher than those of Protocols 2 and 4. Note that the stated results are

max <sup>u</sup> <sup>∈</sup> f g <sup>1</sup>;…;<sup>n</sup>

<sup>j</sup>¼<sup>1</sup> wijð Þ¼ <sup>k</sup> 1 and

285

<sup>j</sup>¼<sup>1</sup> wijð Þ¼ <sup>k</sup> 1 is satisfied. The weights of Protocols 1

� � <sup>w</sup>iið Þ <sup>k</sup>

<sup>n</sup> <sup>1</sup> �<sup>X</sup>

<sup>1</sup>þmax dgið Þ<sup>k</sup> ;dgj f g ð Þ<sup>k</sup> <sup>1</sup> �<sup>X</sup>

<sup>2</sup> <sup>1</sup> �<sup>X</sup>

dgu f g ð Þ<sup>k</sup> <sup>1</sup> �<sup>X</sup>

ð Þ¼ k #N i, <sup>k</sup>. Protocol 4 is only applicable when each node has

u6¼i wiuð Þk

u6¼i wiuð Þk

u6¼i wiuð Þk

u6¼i wiuð Þk

ð Þ¼ k N i, <sup>k</sup>. For instance, in network

Consensus-Based Distributed Filtering for GNSS http://dx.doi.org/10.5772/intechopen.71138

wijð Þ<sup>k</sup> <sup>≥</sup> 0 and wiið Þ<sup>k</sup> <sup>&</sup>gt; 0, having row- and column-sums equal to one, i.e. <sup>P</sup><sup>n</sup>

<sup>i</sup>¼<sup>1</sup> wijð Þ¼ <sup>k</sup> 1 (i, j <sup>¼</sup> <sup>1</sup>, …, n), see e.g. [3, 5, 40, 42, 43].

(number of neighbors) of nodes i ¼ 1, …, n, denoted by dgi

wiið Þ<sup>k</sup> are chosen so that the condition <sup>P</sup><sup>n</sup>

of iterative communications.

3.3. On convergence of consensus states

Protocols <sup>w</sup>ijð Þ <sup>k</sup> <sup>i</sup>;<sup>j</sup> � �∈E<sup>k</sup>

Protocol 1 <sup>1</sup>

Protocol 3 <sup>1</sup>

The degree (number of neighbors) of node i is denoted by dgi

Table 1. Examples of average-consensus protocols forming the weights wijð Þk in Eq. (31).

Protocol 2 <sup>1</sup>

Protocol 4 <sup>1</sup>

otherwise wijð Þ¼ k 0

at most one neighbor at a session.

P<sup>n</sup>

with rj,tð Þ0 : ¼ rj,t. Choosing a set of weights wijð Þk , the nodes i ¼ 1, …, n agree on the consensus protocol (31) to iteratively fuse their information vectors ri,tð Þk . Here and in the following, we use the letter 'k' for the 'session number' k ¼ 1, …, ko (cf. 29) and for the 'number of iterative communications' k ¼ 1, …, kn (cf. 34). The maximum iteration kn is assumed to be not smaller than the maximum session number ko, i.e. kn ≥ ko.

The question that now comes to the fore is how to choose the weights wijð Þk such that the approximation ri,tð Þk gets close to rt through the iteration Eq. (31). More precisely, the stated iteration becomes favorable if ri,tð Þ!k rt when k ! ∞ for all nodes i ¼ 1, …, n. To address this question, we use a multivariate formulation. Let p be the size of the vectors ri,t (i ¼ 1, …, n). We define the higher-dimensioned vector <sup>r</sup> <sup>¼</sup> <sup>r</sup><sup>T</sup> 1,t;…;r<sup>T</sup> n,t h i<sup>T</sup> . The multivariate version of Eq. (31) reads then

$$\begin{bmatrix} r\_{1,t}(k) \\ \vdots \\ r\_{i,t}(k) \\ \vdots \\ r\_{n,t}(k) \end{bmatrix} = \begin{bmatrix} w\_{11}(k)I\_p & \dots & w\_{1i}(k)I\_p & \dots & w\_{1n}(k)I\_p \\ \vdots & \ddots & \vdots & \ddots & \vdots \\ w\_{i1}(k)I\_p & \dots & w\_{ii}(k)I\_p & \dots & w\_{in}(k)I\_p \\ \vdots & \ddots & \vdots & \ddots & \vdots \\ w\_{n1}(k)I\_p & \dots & w\_{in}(k)I\_p & \dots & w\_{nn}(k)I\_p \end{bmatrix} \begin{bmatrix} r\_{1,t}(k-1) \\ \vdots \\ r\_{i,t}(k-1) \\ \vdots \\ r\_{n,t}(k-1) \end{bmatrix}, \quad k = 1,2,\dots \tag{32}$$

or

$$r(k) = \left[\mathcal{W}(k) \otimes I\_{\mathcal{P}}\right] r(k-1), \quad k = 1, 2, \dots \tag{33}$$

The n � n weight matrix W kð Þ is structured by wijð Þk (j∈ f g i; N i, <sup>k</sup> ) and wijð Þ¼ k 0.

(j∉f g i; N i, <sup>k</sup> ). The symbol ⊗ is the Kronecker matrix product [41]. According to Eq. (33), after kn iterations the most recent iterated vector r kð Þ<sup>n</sup> is linked to the initial vector rð Þ0 by Qkn <sup>k</sup>¼<sup>1</sup> W kð Þ <sup>⊗</sup> Ip h irð Þ<sup>0</sup> . Thus the vectors ri,tð Þ<sup>k</sup> (<sup>i</sup> <sup>¼</sup> <sup>1</sup>, …, n) converge to rt when

$$L(k\_n) \coloneqq \prod\_{k=1}^{k\_n} W(k) \to \frac{1}{n} \varepsilon\_n e\_{n'}^T \quad \text{as} \quad k\_n \to \infty \tag{34}$$

where the n-vector en contains ones. If the condition Eq. (34) is met, the set of nodes 1f g ;…; n can asymptotically reach average consensus [4]. It can be shown that (34) holds if the weight matrices W kð Þ (k ¼ 1, …, ko) have bounded nonnegative entries with positive diagonals, i.e. wijð Þ<sup>k</sup> <sup>≥</sup> 0 and wiið Þ<sup>k</sup> <sup>&</sup>gt; 0, having row- and column-sums equal to one, i.e. <sup>P</sup><sup>n</sup> <sup>j</sup>¼<sup>1</sup> wijð Þ¼ <sup>k</sup> 1 and P<sup>n</sup> <sup>i</sup>¼<sup>1</sup> wijð Þ¼ <sup>k</sup> 1 (i, j <sup>¼</sup> <sup>1</sup>, …, n), see e.g. [3, 5, 40, 42, 43].

Examples of such consensus protocols are given in Table 1. As shown, the weights form a symmetric weight matrix W kð Þ, i.e. wjið Þ¼ k wijð Þk . In all protocols presented, self-weights wiið Þ<sup>k</sup> are chosen so that the condition <sup>P</sup><sup>n</sup> <sup>j</sup>¼<sup>1</sup> wijð Þ¼ <sup>k</sup> 1 is satisfied. The weights of Protocols 1 and 2 belong to the class of 'maximum-degree' weights, while those of Protocol 3 are referred to as 'Metropolis' weights [8]. The weights of Protocols 1 and 3 are driven by the degrees (number of neighbors) of nodes i ¼ 1, …, n, denoted by dgi ð Þ¼ k N i, <sup>k</sup>. For instance, in network (a) of Figure 2 we have dg1ð Þ¼ 1 4 as node 1 has 4 neighbors, while dg14ð Þ¼ 1 7 as node 14 has 7 neighbors. Protocol 4 is only applicable to networks like (b) in Figure 2, i.e. when each node has at most one neighbor at a session [4]. In this case, each node exchanges its information to just one neighbor at a session. Thus for two neighboring nodes i and j we have wiið Þ¼ k wjjð Þ¼ k wijð Þ¼ k 0:5, each averaging ri,tð Þ k � 1 and rj,tð Þ k � 1 to obtain ri,tð Þ¼ k rj,tð Þk .

To provide insight into the applicability of the protocols given in Table 1, we apply them to the networks of Figure 2. Twenty values (scalars), say ri (i ¼ 1, …, 20), are generated whose average is equal to 5, i.e. r ¼ 5. Each value is assigned to its corresponding node. For network (a), Protocols 1, 2 and 3 are separately applied, whereas Protocol 4 is only applied to network (b). The corresponding results, up to 30 iterations, are presented in Figure 3. As shown, the iterated values rið Þk (i ¼ 1,…, 20) get closer to their average (i.e. r ¼ 5), the more the number of iterative communications.

#### 3.3. On convergence of consensus states

Eq. (30) to those of their own neighbors. Thus the neighboring nodes j∈ N i, <sup>k</sup> also have their own weighted-averages rj,tð Þ1 . But they may have access to those to which the node i has no direct links. In other words, the weighted-averages rj,tð Þ1 (j ∈ N i, 1) contain information on the nodes to which the node i has no access. For the next session k ¼ 2, it is therefore reasonable for the node i to repeat the fusion rule Eq. (30), but now over the new vectors rj,tð Þ1 (j∈ f g i; N i, <sup>2</sup> ), aiming to improve on the earlier approximation ri,tð Þ1 . This yields the following iterative computations

with rj,tð Þ0 : ¼ rj,t. Choosing a set of weights wijð Þk , the nodes i ¼ 1, …, n agree on the consensus protocol (31) to iteratively fuse their information vectors ri,tð Þk . Here and in the following, we use the letter 'k' for the 'session number' k ¼ 1, …, ko (cf. 29) and for the 'number of iterative communications' k ¼ 1, …, kn (cf. 34). The maximum iteration kn is assumed to be not smaller

The question that now comes to the fore is how to choose the weights wijð Þk such that the approximation ri,tð Þk gets close to rt through the iteration Eq. (31). More precisely, the stated iteration becomes favorable if ri,tð Þ!k rt when k ! ∞ for all nodes i ¼ 1, …, n. To address this question, we use a multivariate formulation. Let p be the size of the vectors ri,t (i ¼ 1, …, n). We

> 1,t;…;r<sup>T</sup> n,t h i<sup>T</sup>

r1,tð Þ k � 1 ⋮ ri,tð Þ k � 1 ⋮ rn,tð Þ k � 1

� �r kð Þ � <sup>1</sup> , k <sup>¼</sup> <sup>1</sup>, <sup>2</sup>,… (33)

wijð Þk rj,tð Þ k � 1 , k ¼ 1, 2, … (31)

. The multivariate version of Eq. (31)

, k ¼ 1, 2, … (32)

<sup>n</sup> , as kn ! ∞ (34)

ri,tð Þ¼ <sup>k</sup> <sup>X</sup>

than the maximum session number ko, i.e. kn ≥ ko.

define the higher-dimensioned vector <sup>r</sup> <sup>¼</sup> <sup>r</sup><sup>T</sup>

reads then

or

Qkn

<sup>k</sup>¼<sup>1</sup> W kð Þ <sup>⊗</sup> Ip h i

r1,tð Þk ⋮ ri,tð Þk ⋮ rn,tð Þk

284 Kalman Filters - Theory for Advanced Applications

<sup>j</sup><sup>∈</sup> f g <sup>i</sup>;<sup>N</sup> i, <sup>k</sup>

w11ð Þk Ip … w1ið Þk Ip … w1nð Þk Ip ⋮⋱⋮⋱⋮ wi1ð Þk Ip … wiið Þk Ip … winð Þk Ip ⋮⋱⋮⋱⋮ wn1ð Þk Ip … wnið Þk Ip … wnnð Þk Ip

r kð Þ¼ W kð Þ ⊗ Ip

L kð Þ<sup>n</sup> <sup>≔</sup><sup>Y</sup> kn

k¼1

The n � n weight matrix W kð Þ is structured by wijð Þk (j∈ f g i; N i, <sup>k</sup> ) and wijð Þ¼ k 0.

(j∉f g i; N i, <sup>k</sup> ). The symbol ⊗ is the Kronecker matrix product [41]. According to Eq. (33), after kn iterations the most recent iterated vector r kð Þ<sup>n</sup> is linked to the initial vector rð Þ0 by

rð Þ0 . Thus the vectors ri,tð Þk (i ¼ 1, …, n) converge to rt when

where the n-vector en contains ones. If the condition Eq. (34) is met, the set of nodes 1f g ;…; n can asymptotically reach average consensus [4]. It can be shown that (34) holds if the weight

W kð Þ! <sup>1</sup> n ene T Figure 3 shows that the states ri,tð Þk (i ¼ 1, …, n) converge to their average rt, but with different rates. The convergence rate depends on the initial states ri,tð Þ¼ 0 ri,t and on the consensus protocol employed. From the figure it seems that the convergence rates of Protocols 1 and 3 are about the same, higher than those of Protocols 2 and 4. Note that the stated results are


The degree (number of neighbors) of node i is denoted by dgi ð Þ¼ k #N i, <sup>k</sup>. Protocol 4 is only applicable when each node has at most one neighbor at a session.

Table 1. Examples of average-consensus protocols forming the weights wijð Þk in Eq. (31).

Figure 3. Performances of protocols 1, 2 and 3 (network (a) of Figure 2) and protocols 4 (network (b)) in delivering the average of 20 values (scalars). The iterated values get closer to their average (i.e. r ¼ 5), the more the number of iterative communications.

obtained on the basis of specific 'realizations' of ri,t (i ¼ 1,…, n). Consider the states ri,t to be random vectors. In that case, can the results be still representative for judging the convergence performances of the protocols? To answer this question, let us define the difference vectors dri,tð Þ¼ <sup>k</sup> ri,tð Þ� <sup>k</sup> rt that are collected in the higher-dimensioned vector dr kð Þ¼ dr<sup>T</sup> <sup>1</sup>,tð Þk ;…; h dr<sup>T</sup> n,t ð Þ� <sup>k</sup> <sup>T</sup>. The more the number of iterations, the smaller the norm of dr kð Þ becomes. According to Eq. (33), after kn iterations the difference vector drð Þ kn is linked to <sup>r</sup> <sup>¼</sup> <sup>r</sup><sup>T</sup> 1,t;…;r<sup>T</sup> n,t h i<sup>T</sup> through

$$dr(k\_n) = \left[ \left( L(k\_n) - \frac{1}{n} e\_n e\_n^T \right) \otimes I\_p \right] r \tag{35}$$

with the eigenvalues λ<sup>1</sup> ≤… ≤ λ<sup>n</sup>�<sup>1</sup> < λ<sup>n</sup> ¼ 1, and the corresponding orthogonal unit eigenvec-

The above equation shows that the entries of the variance matrix (36) are largely driven by the second largest eigenvalue of W, i.e. λ<sup>n</sup>�1. The smaller the scalar ∣λ<sup>n</sup>�<sup>1</sup>∣, the faster the quantity

<sup>n</sup>�<sup>1</sup> tends to zero, as kn ! <sup>∞</sup>. The scalar <sup>∣</sup>λ<sup>n</sup>�<sup>1</sup><sup>∣</sup> is thus often used as a measure to judge the convergence performances of the protocols [7]. For the networks of Figure 2, ∣λ<sup>n</sup>�<sup>1</sup>∣ of Protocols 1, 2 and 3 are about 0.92, 0.97, 0.91, respectively. As Protocol 3 has the smallest ∣λ<sup>n</sup>�<sup>1</sup>∣, it is therefore expected to have the best performance. Note, in Protocol 4, that the weight matrix W kð Þ varies in every session, the performance of which cannot be judged by a single eigenvalue λ<sup>n</sup>�1. One can therefore think of another means of measuring the convergence performance. Due to the randomness of the information vectors ri,t (i ¼ 1, …, n), one may propose

<sup>i</sup> k k dri,tð Þ kn <sup>Q</sup> <sup>≤</sup> <sup>q</sup>

to the probability that the maximum-norm of the difference vectors dri,tð Þ¼ kn ri,tð Þ� kn rt (i ¼ 1, …, n) is not larger than a given positive scalar q for a fixed number of iterations kn. The higher the probability Eq. (39), the better the performance of a protocol. For the scalar case

<sup>i</sup> <sup>j</sup>dri,tð Þj kn <sup>≤</sup> <sup>q</sup><sup>σ</sup>

which is the probability that the absolute differences ∣dri,tð Þ kn ∣ (i ¼ 1, …, n) are not larger than q times the standard-deviation σ. For the networks of Figure 2, 100,000 normally-distributed

for Protocols 1, 2, 3 and 4 are presented in Figure 4. The stated probability is plotted as a function of q for three numbers of iterative communications kn ¼ 10, 20 and 30. As shown, Protocol 3 gives rise to highest probabilities, while Protocol 2 delivers lowest probabilities. After 10 iterations, the probability of having absolute differences smaller than one-fifth of the standard-deviation σ (i.e. q ¼ 0:2) is about 80% for Protocol 1, whereas it is less than 5% for Protocol 2. After 30 iterations, the stated probability increases to 80% for Protocol 2, but close

Figure 4 demonstrates that the convergence performance of Protocol 4 is clearly better than that of Protocol 2, as it delivers higher probabilities (for the networks of Figure 2). Such a conclusion however, cannot be made on the basis of the results of Figure 3. This shows that

Prob max

L kð Þ¼ <sup>n</sup> <sup>W</sup>kn . Substitution into Eq. (36), together with Eq. (37), gives finally

Xn�1 i¼1

λ<sup>2</sup>kn <sup>i</sup> uiu<sup>T</sup> i

!

Dð Þ¼ dr kð Þ<sup>n</sup>

Prob max

to evaluate the convergence rates of the protocols, where k k dri,t <sup>2</sup>

<sup>n</sup> <sup>p</sup> ð Þen. By a repeated application of the protocol <sup>W</sup>, we get

⊗ Q ≤ λ<sup>2</sup>kn

� �, qð Þ <sup>&</sup>gt; <sup>0</sup> (39)

<sup>Q</sup> : <sup>¼</sup> dr<sup>T</sup> i,t Q�<sup>1</sup>

� � (40)

<sup>T</sup> are simulated to evaluate the probability (40). The results

dri,t. Eq. (39) refers

<sup>n</sup>�<sup>1</sup> ½ � In <sup>⊗</sup> <sup>Q</sup> (38)

Consensus-Based Distributed Filtering for GNSS http://dx.doi.org/10.5772/intechopen.71138 287

tors <sup>u</sup>1, …, un�<sup>1</sup>, un <sup>¼</sup> <sup>1</sup><sup>=</sup> ffiffiffi

'probabilistic' measures such as

<sup>Q</sup> <sup>¼</sup> <sup>σ</sup>2, Eq. (39) is reduced to

vectors as samples of r ¼ ½ � r1;…;r<sup>20</sup>

to 100% for Protocols 1 and 3.

λ<sup>2</sup>kn

Now let the initial states ri,t have the same mean and the same variance matrix Dð Þ¼ ri,t Q (i ¼ 1, …, n), but mutually uncorrelated. An application of the covariance propagation law to (35), together with L kð Þ<sup>n</sup> en ¼ en, gives

$$\mathbf{D}(dr(k\_n)) = \left(L^2(k\_n) - \frac{1}{n}\varepsilon\_n e\_n^T\right) \otimes \mathbb{Q} \tag{36}$$

Thus the closer the squared matrix L<sup>2</sup> ð Þ kn to 1ð Þ <sup>=</sup><sup>n</sup> ene<sup>T</sup> <sup>n</sup> , the smaller the variance matrix Eq. (36) becomes. In the limit when kn ! ∞, the stated variance matrix tends to zero. This is what one would expect, since dr kð Þ!<sup>n</sup> 0. Under the conditions stated in Eq. (34), matrices W kð Þ have λ<sup>n</sup> ¼ 1 as the largest absolute value of their eigenvalues [42]. A symmetric weight matrix W can then be expressed in its spectral form as

$$\mathcal{W} = \sum\_{i=1}^{n-1} \lambda\_i \mu\_i \mu\_i^T + \frac{1}{n} \mathbf{c}\_n \mathbf{e}\_n^T \tag{37}$$

with the eigenvalues λ<sup>1</sup> ≤… ≤ λ<sup>n</sup>�<sup>1</sup> < λ<sup>n</sup> ¼ 1, and the corresponding orthogonal unit eigenvectors <sup>u</sup>1, …, un�<sup>1</sup>, un <sup>¼</sup> <sup>1</sup><sup>=</sup> ffiffiffi <sup>n</sup> <sup>p</sup> ð Þen. By a repeated application of the protocol <sup>W</sup>, we get L kð Þ¼ <sup>n</sup> <sup>W</sup>kn . Substitution into Eq. (36), together with Eq. (37), gives finally

$$\mathbf{D}(dr(k\_n)) = \left(\sum\_{i=1}^{n-1} \lambda\_i^{2k\_n} u\_i u\_i^T\right) \otimes \mathbf{Q} \le \lambda\_{n-1}^{2k\_n} \left[I\_n \otimes \mathbf{Q}\right] \tag{38}$$

The above equation shows that the entries of the variance matrix (36) are largely driven by the second largest eigenvalue of W, i.e. λ<sup>n</sup>�1. The smaller the scalar ∣λ<sup>n</sup>�<sup>1</sup>∣, the faster the quantity λ<sup>2</sup>kn <sup>n</sup>�<sup>1</sup> tends to zero, as kn ! <sup>∞</sup>. The scalar <sup>∣</sup>λ<sup>n</sup>�<sup>1</sup><sup>∣</sup> is thus often used as a measure to judge the convergence performances of the protocols [7]. For the networks of Figure 2, ∣λ<sup>n</sup>�<sup>1</sup>∣ of Protocols 1, 2 and 3 are about 0.92, 0.97, 0.91, respectively. As Protocol 3 has the smallest ∣λ<sup>n</sup>�<sup>1</sup>∣, it is therefore expected to have the best performance. Note, in Protocol 4, that the weight matrix W kð Þ varies in every session, the performance of which cannot be judged by a single eigenvalue λ<sup>n</sup>�1. One can therefore think of another means of measuring the convergence performance. Due to the randomness of the information vectors ri,t (i ¼ 1, …, n), one may propose 'probabilistic' measures such as

$$\text{Prob}\left(\max\_{i} \|dr\_{i,t}(k\_n)\|\_Q \le q\right), \quad (q > 0) \tag{39}$$

to evaluate the convergence rates of the protocols, where k k dri,t <sup>2</sup> <sup>Q</sup> : <sup>¼</sup> dr<sup>T</sup> i,t Q�<sup>1</sup> dri,t. Eq. (39) refers to the probability that the maximum-norm of the difference vectors dri,tð Þ¼ kn ri,tð Þ� kn rt (i ¼ 1, …, n) is not larger than a given positive scalar q for a fixed number of iterations kn. The higher the probability Eq. (39), the better the performance of a protocol. For the scalar case <sup>Q</sup> <sup>¼</sup> <sup>σ</sup>2, Eq. (39) is reduced to

obtained on the basis of specific 'realizations' of ri,t (i ¼ 1,…, n). Consider the states ri,t to be random vectors. In that case, can the results be still representative for judging the convergence performances of the protocols? To answer this question, let us define the difference vectors

Figure 3. Performances of protocols 1, 2 and 3 (network (a) of Figure 2) and protocols 4 (network (b)) in delivering the average of 20 values (scalars). The iterated values get closer to their average (i.e. r ¼ 5), the more the number of iterative

ð Þ� <sup>k</sup> <sup>T</sup>. The more the number of iterations, the smaller the norm of dr kð Þ becomes. According

1 n ene T n

� �

⊗ Ip

� �

Now let the initial states ri,t have the same mean and the same variance matrix Dð Þ¼ ri,t Q (i ¼ 1, …, n), but mutually uncorrelated. An application of the covariance propagation law to

> ð Þ� kn 1 n ene T n

ð Þ kn to 1ð Þ <sup>=</sup><sup>n</sup> ene<sup>T</sup>

becomes. In the limit when kn ! ∞, the stated variance matrix tends to zero. This is what one would expect, since dr kð Þ!<sup>n</sup> 0. Under the conditions stated in Eq. (34), matrices W kð Þ have λ<sup>n</sup> ¼ 1 as the largest absolute value of their eigenvalues [42]. A symmetric weight matrix W

> λ<sup>i</sup> uiu<sup>T</sup> <sup>i</sup> þ 1 n ene T

� �

<sup>1</sup>,tð Þk ;…;

through

h

1,t;…;r<sup>T</sup> n,t h i<sup>T</sup>

r (35)

⊗ Q (36)

<sup>n</sup> , the smaller the variance matrix Eq. (36)

<sup>n</sup> (37)

dri,tð Þ¼ <sup>k</sup> ri,tð Þ� <sup>k</sup> rt that are collected in the higher-dimensioned vector dr kð Þ¼ dr<sup>T</sup>

to Eq. (33), after kn iterations the difference vector drð Þ kn is linked to <sup>r</sup> <sup>¼</sup> <sup>r</sup><sup>T</sup>

dr kð Þ¼ <sup>n</sup> L kð Þ� <sup>n</sup>

<sup>D</sup>ð Þ¼ dr kð Þ<sup>n</sup> <sup>L</sup><sup>2</sup>

<sup>W</sup> <sup>¼</sup> <sup>X</sup><sup>n</sup>�<sup>1</sup> i¼1

dr<sup>T</sup> n,t

communications.

286 Kalman Filters - Theory for Advanced Applications

(35), together with L kð Þ<sup>n</sup> en ¼ en, gives

Thus the closer the squared matrix L<sup>2</sup>

can then be expressed in its spectral form as

$$\text{Prob}\left(\max\_{i}|dr\_{i,t}(k\_n)| \le q\sigma\right) \tag{40}$$

which is the probability that the absolute differences ∣dri,tð Þ kn ∣ (i ¼ 1, …, n) are not larger than q times the standard-deviation σ. For the networks of Figure 2, 100,000 normally-distributed vectors as samples of r ¼ ½ � r1;…;r<sup>20</sup> <sup>T</sup> are simulated to evaluate the probability (40). The results for Protocols 1, 2, 3 and 4 are presented in Figure 4. The stated probability is plotted as a function of q for three numbers of iterative communications kn ¼ 10, 20 and 30. As shown, Protocol 3 gives rise to highest probabilities, while Protocol 2 delivers lowest probabilities. After 10 iterations, the probability of having absolute differences smaller than one-fifth of the standard-deviation σ (i.e. q ¼ 0:2) is about 80% for Protocol 1, whereas it is less than 5% for Protocol 2. After 30 iterations, the stated probability increases to 80% for Protocol 2, but close to 100% for Protocols 1 and 3.

Figure 4 demonstrates that the convergence performance of Protocol 4 is clearly better than that of Protocol 2, as it delivers higher probabilities (for the networks of Figure 2). Such a conclusion however, cannot be made on the basis of the results of Figure 3. This shows that

than the sampling rate Δ so as to be able to incorporate consensus protocols into the information filter setup. Such a consensus-based Kalman filter (CKF) would thus generally be of a two time-scale nature [2], the data sampling time-scale t ¼ 1, 2, …, versus the data sending timescale k ¼ 1, …, kn. The CKF is a suitable tool for handling real-time data processing in a distributed manner for the applications in which the state-vectors xt (t ¼ 1, 2,…) change rather slowly over time (i.e. Δ can take large values) and/or for the cases where the sensor nodes

Figure 5. The two time-scale nature of a consensus-based Kalman filter (CKF): The data sampling time-scale t ¼ 1, 2, …, versus the data sending time-scale k ¼ 1,…, kn. The sending rate δ must be reasonably smaller than the sampling rate Δ so

Under the assumption δ ≤ Δ, the CKF recursion follows from the Kalman filter recursion by considering an extra step, namely, the 'consensus update'. The algorithmic steps of the CKF in information form are presented in Figure 6. Compare the recursion with that of the information filter given in Figure 1. Similar to the information filter, the CKF at node i is initialized by

local normal matrix Ni,t and right-hand-side vector ri,t to send them to its neighboring nodes j∈ N i, <sup>k</sup> (k ¼ 1, …, kn). In the consensus update, iterative communications between the neighboring nodes f g i; N i, <sup>k</sup> are carried out to approximate the averages Nt and rt by Ni,tð Þ kn and ri,tð Þ kn , respectively. After a finite number of communications kn, the consensus states Ni,tð Þ kn and ri,tð Þ kn are, respectively, added to the time update information Ii,t∣t�<sup>1</sup> and ii,t∣t�<sup>1</sup> to obtain their measurement update version Ii,t∣<sup>t</sup> and ii,t∣<sup>t</sup> at node i (cf. 28). The time update goes along

With the consensus-based information filter, presented in Figure 6, it is therefore feasible to develop multiple distributed filters, all running in parallel over time. By taking recourse to an average-consensus protocol, not all the nodes are needed to be directly linked, thereby allowing non-neighboring nodes to also benefit from information states of each other. The price one has to pay for such an attractive feature of the CKF is that the local predictors

will have a poorer precision performance than that of their central counterpart <sup>b</sup>xt∣<sup>t</sup>. This is due to the fact that the consensus states Ni,tð Þ kn and ri,tð Þ kn ð Þ i ¼ 1;…; n are just approximations of

i,t∣<sup>t</sup> ii,t∣t, i ¼ 1, …, n, (41)

, node i computes its

Consensus-Based Distributed Filtering for GNSS http://dx.doi.org/10.5772/intechopen.71138 289

transfer their data rather quickly (i.e. δ can take small values).

as to be able to incorporate consensus protocols into the CKF setup.

the same lines as that of the information filter.

4.2. Time evolution of the CKF error covariances

the zero information Ii,1∣<sup>0</sup> ¼ 0 and ii,<sup>1</sup>∣<sup>0</sup> ¼ 0. In the presence of the data yi,t

<sup>b</sup>xi,t∣<sup>t</sup> <sup>¼</sup> <sup>I</sup>�<sup>1</sup>

Figure 4. Probability Eq. (40) as a function of q for three numbers of iterative communications kn ¼ 10, 20 and 30 (protocols 1 (P.1), 2 (P.2), 3 (P.3) and 4 (P.4)). It refers to the probability that the absolute differences ∣ri,tð Þ� kn rt∣ (i ¼ 1,…, 20) are not larger than q times the standard-deviation σ (cf. Figure 2).

results obtained on the basis of specific 'realizations' of ri,t (i ¼ 1, …, n) are not necessarily representative.

### 4. Consensus-based Kalman filters

#### 4.1. Two time-scale approach

In Section 2.5 we discussed how the additivity property of the measurement update Eq. (26) offers possibilities for developing multiple distributed local filters i ¼ 1, …, n, each delivering local states Ii,t∣<sup>t</sup> and ii,t∣<sup>t</sup> equal to their central counterparts I<sup>t</sup>∣<sup>t</sup> and it∣<sup>t</sup>. In doing so, each node has to evaluate the averages Nt and rt at every time instance t. Since in practice the nodes do not necessarily have direct connections to each other, options such as the consensus-based fusion rules (cf. Section 3) can alternatively be employed to 'approximate' Nt and rt. As illustrated in Figures 3 and 4, such consensus-based approximation requires a number of iterative communications between the nodes in order to reach the averages Nt and rt. The stated iterative communications clearly require some time to be carried out and must take place during every time interval ½ � t; t þ 1 (see Figure 5). We distinguish between the sampling rate Δ and the sending rate δ. The sampling rate refers to the frequency with which the node i collects its observables yi,t (t ¼ 1, 2, …), while the sending rate refers to the frequency with which the node i sends/receives information Nj,tð Þk and rj,tð Þk (k ¼ 1, …, kn) to/from its neighboring nodes. As shown in Figure 5, the sending rate δ should therefore be reasonably smaller

Figure 5. The two time-scale nature of a consensus-based Kalman filter (CKF): The data sampling time-scale t ¼ 1, 2, …, versus the data sending time-scale k ¼ 1,…, kn. The sending rate δ must be reasonably smaller than the sampling rate Δ so as to be able to incorporate consensus protocols into the CKF setup.

than the sampling rate Δ so as to be able to incorporate consensus protocols into the information filter setup. Such a consensus-based Kalman filter (CKF) would thus generally be of a two time-scale nature [2], the data sampling time-scale t ¼ 1, 2, …, versus the data sending timescale k ¼ 1, …, kn. The CKF is a suitable tool for handling real-time data processing in a distributed manner for the applications in which the state-vectors xt (t ¼ 1, 2,…) change rather slowly over time (i.e. Δ can take large values) and/or for the cases where the sensor nodes transfer their data rather quickly (i.e. δ can take small values).

Under the assumption δ ≤ Δ, the CKF recursion follows from the Kalman filter recursion by considering an extra step, namely, the 'consensus update'. The algorithmic steps of the CKF in information form are presented in Figure 6. Compare the recursion with that of the information filter given in Figure 1. Similar to the information filter, the CKF at node i is initialized by the zero information Ii,1∣<sup>0</sup> ¼ 0 and ii,<sup>1</sup>∣<sup>0</sup> ¼ 0. In the presence of the data yi,t , node i computes its local normal matrix Ni,t and right-hand-side vector ri,t to send them to its neighboring nodes j∈ N i, <sup>k</sup> (k ¼ 1, …, kn). In the consensus update, iterative communications between the neighboring nodes f g i; N i, <sup>k</sup> are carried out to approximate the averages Nt and rt by Ni,tð Þ kn and ri,tð Þ kn , respectively. After a finite number of communications kn, the consensus states Ni,tð Þ kn and ri,tð Þ kn are, respectively, added to the time update information Ii,t∣t�<sup>1</sup> and ii,t∣t�<sup>1</sup> to obtain their measurement update version Ii,t∣<sup>t</sup> and ii,t∣<sup>t</sup> at node i (cf. 28). The time update goes along the same lines as that of the information filter.

#### 4.2. Time evolution of the CKF error covariances

results obtained on the basis of specific 'realizations' of ri,t (i ¼ 1, …, n) are not necessarily

Figure 4. Probability Eq. (40) as a function of q for three numbers of iterative communications kn ¼ 10, 20 and 30 (protocols 1 (P.1), 2 (P.2), 3 (P.3) and 4 (P.4)). It refers to the probability that the absolute differences ∣ri,tð Þ� kn rt∣ (i ¼ 1,…, 20)

In Section 2.5 we discussed how the additivity property of the measurement update Eq. (26) offers possibilities for developing multiple distributed local filters i ¼ 1, …, n, each delivering local states Ii,t∣<sup>t</sup> and ii,t∣<sup>t</sup> equal to their central counterparts I<sup>t</sup>∣<sup>t</sup> and it∣<sup>t</sup>. In doing so, each node has to evaluate the averages Nt and rt at every time instance t. Since in practice the nodes do not necessarily have direct connections to each other, options such as the consensus-based fusion rules (cf. Section 3) can alternatively be employed to 'approximate' Nt and rt. As illustrated in Figures 3 and 4, such consensus-based approximation requires a number of iterative communications between the nodes in order to reach the averages Nt and rt. The stated iterative communications clearly require some time to be carried out and must take place during every time interval ½ � t; t þ 1 (see Figure 5). We distinguish between the sampling rate Δ and the sending rate δ. The sampling rate refers to the frequency with which the node i collects its observables yi,t (t ¼ 1, 2, …), while the sending rate refers to the frequency with which the node i sends/receives information Nj,tð Þk and rj,tð Þk (k ¼ 1, …, kn) to/from its neighboring nodes. As shown in Figure 5, the sending rate δ should therefore be reasonably smaller

representative.

4. Consensus-based Kalman filters

are not larger than q times the standard-deviation σ (cf. Figure 2).

4.1. Two time-scale approach

288 Kalman Filters - Theory for Advanced Applications

With the consensus-based information filter, presented in Figure 6, it is therefore feasible to develop multiple distributed filters, all running in parallel over time. By taking recourse to an average-consensus protocol, not all the nodes are needed to be directly linked, thereby allowing non-neighboring nodes to also benefit from information states of each other. The price one has to pay for such an attractive feature of the CKF is that the local predictors

$$
\widehat{\mathfrak{X}}\_{i,t|t} = \mathcal{T}\_{i,t|t}^{-1} \quad \mathfrak{i}\_{i,t|t} \quad \mathfrak{i} = \mathbf{1}, \ldots, n,\tag{41}
$$

will have a poorer precision performance than that of their central counterpart <sup>b</sup>xt∣<sup>t</sup>. This is due to the fact that the consensus states Ni,tð Þ kn and ri,tð Þ kn ð Þ i ¼ 1;…; n are just approximations of

xt � <sup>b</sup>xi,t∣<sup>t</sup> <sup>¼</sup> <sup>I</sup>�<sup>1</sup>

entries of the product matrix L kð Þ<sup>n</sup> in Eq. (34), one obtains4

lij rj,t � Nj,t xt

8 < :

shown, the node i would need an extra input, i.e., the term P<sup>n</sup>

time instance <sup>t</sup>, the stated variance matrix reads Pt∣<sup>t</sup> <sup>¼</sup> ð Þ <sup>1</sup>=tn <sup>N</sup>�<sup>1</sup>

Xn j¼1

<sup>j</sup>¼<sup>1</sup> lij <sup>¼</sup> 1, give

in which the scalar α is given by

matrix Pi,t∣<sup>t</sup> as follows (Figure 7)

Ni,tð Þ¼ kn

j¼1

law to (42) results in the error variance matrix

Pi,t∣<sup>t</sup> <sup>¼</sup> <sup>I</sup>�<sup>1</sup>

Note that the terms xt � <sup>b</sup>xi,t∣t�<sup>1</sup>

ri,tð Þ� kn Ni,tð Þ kn xt <sup>½</sup> � ¼ <sup>X</sup><sup>n</sup>

that is not necessarily equal to I�<sup>1</sup>

since D rj,t � Nj,t xt

together with P<sup>n</sup>

i,t∣<sup>t</sup> <sup>I</sup>i,t∣t�<sup>1</sup> xt � <sup>b</sup>xi,t∣t�<sup>1</sup>

� � � n ri,tð Þ� kn Ni,tð Þ kn xt ½ � � � (42)

j¼1 l 2

Consensus-Based Distributed Filtering for GNSS http://dx.doi.org/10.5772/intechopen.71138

i,t∣<sup>t</sup> (44)

ijNj,t in order to be able to

. We now compare Pt∣<sup>t</sup> with its

ij Nj,t ¼ αnN (45)

ij (46)

ij Nj,t (43)

291

� � and ri,tð Þ� kn Ni,tð Þ kn xt ½ � are uncorrelated. With lij as the

� �, ) <sup>D</sup> ri,tð Þ� kn Ni,tð Þ kn xt ð½ �Þ ¼ <sup>X</sup><sup>n</sup>

j¼1 l 2 ij Nj,t

i,t∣<sup>t</sup> (see the following discussion on Eqs. (47) and (48)).

<sup>j</sup>¼<sup>1</sup> <sup>l</sup> 2

9 = ; I�<sup>1</sup>

� � <sup>¼</sup> Nj,t. With this in mind, an application of the covariance propagation

i,t∣<sup>t</sup> <sup>I</sup>i,t∣t�<sup>1</sup> Pi,t∣t�<sup>1</sup> <sup>I</sup>i,t∣t�<sup>1</sup> <sup>þ</sup> <sup>n</sup><sup>2</sup> <sup>X</sup><sup>n</sup>

In Figure 7 we present the three-step recursion of the error variance matrix Pi,t∣<sup>t</sup> (for node i). As

compute Pi,t∣<sup>t</sup>. In practice however, such additional information is absent in the CKF setup. This means that the node i does not have enough information to evaluate the error variance matrix Pi,t∣<sup>t</sup>. Despite such restriction, it will be shown in Section 5 how the recursion of Pi,t∣<sup>t</sup> conveys useful information about the performance of the local filters i ¼ 1, …, n, thereby allowing one to a-priori design and analyze sensor networks with different numbers of iterative communications. To better appreciate the recursion given in Figure 7, let us consider a special case where a stationary state-vector xt is to be predicted over time. Thus Φt,t�<sup>1</sup> ¼ I and St ¼ 0 ð Þ t ¼ 1; 2;… . Moreover, we assume that all nodes deliver the same normal matrices Ni,t ¼ N ið Þ ¼ 1;…; n . The central error variance matrix Pt∣<sup>t</sup> would then simply follow by inverting the sum of all normal matrices over n nodes and t time instances. Collecting observables up to and including

consensus-based local counterpart at node i, i.e. Pi,t∣<sup>t</sup>. The aforementioned assumptions,

lij Nj,t <sup>¼</sup> N, and <sup>n</sup><sup>2</sup>X<sup>n</sup>

Xn j¼1 l 2

α : ¼ n

Substitution into the stated recursion provides us with the time-evolution of the error variance

j¼1 l 2

$$t\_{\mathbb{J}\mathbf{1}|\mathbb{C}} = 0, \qquad \tau\_{\mathbb{J}\mathbf{1}|\mathbb{O}} = 0.$$

$$r\_{i,t} = \mathcal{A}\_{i,t}^{\top} \mathcal{R}\_{i,t}^{-1} y\_{i,t} \quad \mathcal{A}\_{i,t}$$

$$r\_{i,t} = \mathcal{A}\_{i,t}^{\top} \mathcal{R}\_{i,t}^{-1} y\_{i,t} \quad \quad \mathcal{N}\_{i,t} = \mathcal{A}\_{i,t}^{\top} \mathcal{R}\_{i,t}^{-1} \mathcal{A}\_{i,t}$$

$$r\_{i,t}(0) \stackrel{\sim}{\rightarrow} r\_{i,t} \quad \quad \quad \mathcal{N}\_{i,t}(0) \stackrel{\sim}{\rightarrow} \mathcal{N}\_{i,t}$$

$$\left[\mathbb{N}\_{\vec{\imath},\sharp}(k+1), r\_{\vec{\imath},\sharp}(k+1)\right] - \sum\_{j\in\{\vec{\imath},\imath\_{\imath}\mathbb{N}\_{\imath},\dagger\}} \text{tr}\_{\vec{\imath},\sharp}(k) [\mathbb{N}\_{j,\sharp}(k), r\_{j,\sharp}(k)]$$

$$\mathfrak{i}\_{\sharp,\sharp\vert\mathfrak{t}} = \mathfrak{i}\_{\sharp,\sharp\vert\mathfrak{t}} \quad \mathfrak{1} \quad \mathfrak{n} \ r\_{\sharp,\sharp}(\mathfrak{k}\_{\flat}) , \qquad \mathfrak{I}\_{\sharp,\sharp\vert\mathfrak{t}} = \mathfrak{I}\_{\sharp\sharp\vert\mathfrak{t}} \quad \mathfrak{1} \quad \mathfrak{n} \ \mathfrak{N}\_{\sharp,\sharp}(\mathfrak{k}\_{\flat})$$

$$\begin{array}{rcl} \stackrel{\scriptstyle{\scriptstyle{\mathcal{U}}}\stackrel{\scriptstyle{\scriptstyle{\mathcal{U}}}\stackrel{\scriptstyle{\scriptstyle{\mathcal{U}}}}{\scriptstyle{\mathcal{E}}}}{\scriptstyle{\mathcal{E}}}\Phi\_{\textit{i},\textit{\ell}-1},\ \mathcal{S}\_{\textit{i}}^{-1} \\\ M\_{\textit{i},\textit{t}} & = & \Phi\_{\textit{i},\textit{t}-1}^{-\widetilde{\mathcal{X}}}\mathcal{X}\_{\textit{i},\textit{t}-1}\Phi\_{\textit{i},\textit{t}-1}^{-1} \\\ T\_{\textit{i},\textit{t}|\textit{t}-1} & = & \mathcal{M}\_{\textit{i},\textit{t}} \quad \mathcal{M}\_{\textit{i},\textit{t}}(\mathcal{M}\_{\textit{i}\boldsymbol{t}} + \mathcal{S}\_{\textit{i}}^{-1}) \quad ^{\scriptstyle 1}\mathcal{M}\_{\textit{i},\textit{t}} \\\ t\_{\textit{i},\textit{t}} & = & \mathcal{X}\_{\textit{i},\textit{t}|\textit{t}-1}\Phi\_{\textit{i},\textit{t}-1}\mathcal{X}\_{\textit{i},\textit{t}-1}^{-1} \mathcal{A}\_{\textit{i},\textit{t}-1,\textit{t}-1} \\\ \text{co }\mathsf{S}\textbf{t}\textbf{c}\textbf{p} \; 2 \end{array}$$

Figure 6. Algorithmic steps of the CKF in information form concerning the time-evolution of the local information vector ii,t∣<sup>t</sup> and matrix Ii,t∣<sup>t</sup> of node i.

the averages Nt and rt. Although they reach the stated averages as kn ! ∞, one of course always comes up with a finite number of communications kn. As a consequence, while the inverse-matrix I�<sup>1</sup> <sup>t</sup>∣<sup>t</sup> represents the error variance matrix Pt∣<sup>t</sup> <sup>¼</sup> <sup>D</sup> xt � <sup>b</sup>xt∣<sup>t</sup> � � (cf. 17), the inversematrices I�<sup>1</sup> i,t∣<sup>t</sup> ð Þ <sup>i</sup> <sup>¼</sup> <sup>1</sup>;…; <sup>n</sup> do not represent the error variance matrices Pi,t∣<sup>t</sup> <sup>¼</sup> <sup>D</sup> xt � <sup>b</sup>xi,t∣<sup>t</sup> � �. To see this, consider the local prediction errors xt � <sup>b</sup>xi,t∣<sup>t</sup> � � which can be expressed as (Figure 6)

$$\mathbf{x}\_t - \widehat{\mathbf{x}}\_{i,t|t} = \mathcal{T}\_{i,t|t}^{-1} \left\{ \mathcal{T}\_{i,t|t-1} \left( \mathbf{x}\_t - \widehat{\mathbf{x}}\_{i,t|t-1} \right) - n \left[ r\_{i,t}(\mathbf{k}\_n) - \mathcal{N}\_{i,t}(\mathbf{k}\_n) \mathbf{x}\_t \right] \right\} \tag{42}$$

Note that the terms xt � <sup>b</sup>xi,t∣t�<sup>1</sup> � � and ri,tð Þ� kn Ni,tð Þ kn xt ½ � are uncorrelated. With lij as the entries of the product matrix L kð Þ<sup>n</sup> in Eq. (34), one obtains4

$$\mathbb{E}\left[r\_{i,t}(k\_n) - \mathcal{N}\_{i,t}(k\_n)\mathbf{x}\_t\right] = \sum\_{j=1}^n l\_{\vec{\eta}} (r\_{\bar{\eta},t} - \mathcal{N}\_{\bar{\mu},t}\mathbf{x}\_t)\_{\prime} \quad \Rightarrow \quad \mathbf{D}([r\_{i,t}(k\_n) - \mathcal{N}\_{i,t}(k\_n)\mathbf{x}\_t]) = \sum\_{j=1}^n l\_{\vec{\eta}}^2 \mathcal{N}\_{\bar{\mu},t} \tag{43}$$

since D rj,t � Nj,t xt � � <sup>¼</sup> Nj,t. With this in mind, an application of the covariance propagation law to (42) results in the error variance matrix

$$P\_{i,t|t} = \mathcal{Z}\_{i,t|t}^{-1} \left\{ \mathcal{Z}\_{i,t|t-1} P\_{i,t|t-1} \mathcal{Z}\_{i,t|t-1} + n^2 \sum\_{j=1}^n l\_{ij}^2 N\_{j,t} \right\} \mathcal{Z}\_{i,t|t}^{-1} \tag{44}$$

that is not necessarily equal to I�<sup>1</sup> i,t∣<sup>t</sup> (see the following discussion on Eqs. (47) and (48)).

In Figure 7 we present the three-step recursion of the error variance matrix Pi,t∣<sup>t</sup> (for node i). As shown, the node i would need an extra input, i.e., the term P<sup>n</sup> <sup>j</sup>¼<sup>1</sup> <sup>l</sup> 2 ijNj,t in order to be able to compute Pi,t∣<sup>t</sup>. In practice however, such additional information is absent in the CKF setup. This means that the node i does not have enough information to evaluate the error variance matrix Pi,t∣<sup>t</sup>. Despite such restriction, it will be shown in Section 5 how the recursion of Pi,t∣<sup>t</sup> conveys useful information about the performance of the local filters i ¼ 1, …, n, thereby allowing one to a-priori design and analyze sensor networks with different numbers of iterative communications.

To better appreciate the recursion given in Figure 7, let us consider a special case where a stationary state-vector xt is to be predicted over time. Thus Φt,t�<sup>1</sup> ¼ I and St ¼ 0 ð Þ t ¼ 1; 2;… . Moreover, we assume that all nodes deliver the same normal matrices Ni,t ¼ N ið Þ ¼ 1;…; n .

The central error variance matrix Pt∣<sup>t</sup> would then simply follow by inverting the sum of all normal matrices over n nodes and t time instances. Collecting observables up to and including time instance <sup>t</sup>, the stated variance matrix reads Pt∣<sup>t</sup> <sup>¼</sup> ð Þ <sup>1</sup>=tn <sup>N</sup>�<sup>1</sup> . We now compare Pt∣<sup>t</sup> with its consensus-based local counterpart at node i, i.e. Pi,t∣<sup>t</sup>. The aforementioned assumptions, together with P<sup>n</sup> <sup>j</sup>¼<sup>1</sup> lij <sup>¼</sup> 1, give

$$N\_{i,t}(k\_n) = \sum\_{j=1}^{n} l\_{\vec{\eta}} N\_{j,t} = N\_\prime \quad \text{and} \quad n^2 \sum\_{j=1}^{n} l\_{\vec{\eta}}^2 N\_{j,t} = \alpha nN \tag{45}$$

in which the scalar α is given by

the averages Nt and rt. Although they reach the stated averages as kn ! ∞, one of course always comes up with a finite number of communications kn. As a consequence, while the

Figure 6. Algorithmic steps of the CKF in information form concerning the time-evolution of the local information

i,t∣<sup>t</sup> ð Þ <sup>i</sup> <sup>¼</sup> <sup>1</sup>;…; <sup>n</sup> do not represent the error variance matrices Pi,t∣<sup>t</sup> <sup>¼</sup> <sup>D</sup> xt � <sup>b</sup>xi,t∣<sup>t</sup>

� � (cf. 17), the inverse-

� � which can be expressed as (Figure 6)

� �. To

<sup>t</sup>∣<sup>t</sup> represents the error variance matrix Pt∣<sup>t</sup> <sup>¼</sup> <sup>D</sup> xt � <sup>b</sup>xt∣<sup>t</sup>

inverse-matrix I�<sup>1</sup>

vector ii,t∣<sup>t</sup> and matrix Ii,t∣<sup>t</sup> of node i.

290 Kalman Filters - Theory for Advanced Applications

see this, consider the local prediction errors xt � <sup>b</sup>xi,t∣<sup>t</sup>

matrices I�<sup>1</sup>

$$\alpha := n \sum\_{j=1}^{n} l\_{ij}^{2} \tag{46}$$

Substitution into the stated recursion provides us with the time-evolution of the error variance matrix Pi,t∣<sup>t</sup> as follows (Figure 7)

α ¼ e T <sup>n</sup> en � � l

Eq. (34), this can be realized if L kð Þ!<sup>n</sup> ð Þ <sup>1</sup>=<sup>n</sup> ene<sup>T</sup>

5. Applications to GNSS

actual error variance matrices of the predictors <sup>b</sup>xi,t∣<sup>t</sup>, i.e. <sup>I</sup>�<sup>1</sup>

as l

solutions.

about 8 km).

Tl � � ≥ l

Ten <sup>¼</sup> 1. Thus scalar <sup>α</sup> is never smaller than 1, i.e. Pi,t∣<sup>t</sup> <sup>≥</sup> Pt∣<sup>t</sup>, showing that the performance of the consensus-based predictor <sup>b</sup>xi,t∣<sup>t</sup> is never better than that of its central version <sup>b</sup>xt∣<sup>t</sup>. The lower-bound Eq. (48) is reached when l ¼ ð Þ 1=n en, i.e. when lij ¼ 1=n (j ¼ 1, …, n). According to

be required to be reasonably large. The conclusion reads therefore that the local filters at nodes i ¼ 1, …, n, generate information matrices Ii,t∣<sup>t</sup>, the inverse of which are different from the

The purpose of this section is to demonstrate how the CKF theory, discussed in Section 4, can play a pivotal role in applications for which the GNSS measurements of a network of receivers are to be processed in a real-time manner. In a GNSS network setup, each receiver serves as a sensor node for receiving observables from visible GNSS satellites to determine a range of different parameters such as positions and velocities in an Earth-centered Earth-fixed coordinate system, atmospheric delays, timing and instrumental biases, see e.g. [11, 12]. As the observation equations of the receivers have satellite specific parameters in common, the receivers' observables are often integrated through a computing (fusion) center to provide network-derived parameter solutions that are more precise than their single-receiver versions. Now the idea is to deliver GNSS parameter solutions without the need of having a computing center, such that their precision performance is still comparable to that of network-derived

As previously discussed, consensus-based algorithms and in particular the CKF can be employed to process network data in a distributed filtering scheme, i.e. no computing center is required. In order to illustrate such applicability, we simulate a network of 13 GNSS receivers located in Perth, Western Australia (Figure 8). As shown in the figure, each node (white circle) represents a receiver having data links (red lines) to its neighbors with interstation distances up to 4 km. We therefore assume that the receivers receive each other data within the ranges not longer than 4 km. For instance, receiver 1 is directly connected to receivers 2 and 6, but not to receiver 3 (the inter-station distance between receivers 1 and 3 is

Although the GNSS observables contain information on various positioning and nonpositioning parameters, here we restrict ourselves to ionospheric observables of the GPS pseudorange measurements only [44]. One should however bear in mind that such restriction is made just for the sake of presentation and illustration of the theory discussed in Sections 3 and 4.

5.1. GNSS ionospheric observables: Dynamic and measurement models

Ten

� �<sup>2</sup> <sup>¼</sup> <sup>1</sup> (48)

Consensus-Based Distributed Filtering for GNSS http://dx.doi.org/10.5772/intechopen.71138 293

<sup>n</sup> , for which the number of iterations kn might

i,t∣<sup>t</sup> 6¼ Pi,t∣<sup>t</sup>.

$$\mathcal{I}\_{\text{i},\text{l}\vert 0} = 0, \qquad F\_{\text{i},\text{l}} = \mathbf{0}$$

$$\begin{array}{cccc} & & & & \mathbf{M}\_{\boldsymbol{\upbeta},\boldsymbol{\upbeta}}^{\text{aug}}(\mathbf{N}\_{\boldsymbol{\upbeta},\boldsymbol{\upbeta}}(\mathbf{k}\_{\boldsymbol{\upbeta}}),\mathbf{N}\_{\boldsymbol{\upbeta}}^{\text{aug}},\mathbf{N}\_{\boldsymbol{\upbeta},\boldsymbol{\upbeta}}^{\text{aug}})\\ \mathbf{Z}\_{\boldsymbol{\upbeta},\boldsymbol{\upbeta}} & = & \mathbf{Z}\_{\boldsymbol{\upbeta},\boldsymbol{\upbeta}}^{\text{aug}}(\mathbf{I}\_{\boldsymbol{\upbeta}}+n\mathbf{N}\_{\boldsymbol{\upbeta},\boldsymbol{\upbeta}}^{\text{aug}}(\mathbf{k}\_{\boldsymbol{\upbeta}})\\ \mathbf{P}\_{\boldsymbol{\upbeta},\boldsymbol{\upbeta}}^{\text{aug}} & = & \mathbf{Z}\_{\boldsymbol{\upbeta},\boldsymbol{\upbeta}}^{-1}(\mathbf{F}\_{\boldsymbol{\upbeta},\boldsymbol{\upbeta}}+n^{2}\sum\_{j=1}^{n}l\_{j,j}^{2}\mathbf{N}\_{j,\boldsymbol{\upbeta}}^{\text{aug}}\mathbf{Z}\_{j,\boldsymbol{\upbeta}}^{-1}) \end{array}$$

$$\begin{array}{rcl} & \stackrel{sim\_{\mathit{m}}\mathit{m}\_{\mathit{i}}}{\rightarrow} & \Phi\_{\mathit{i},\mathit{i}-\mathit{i},\mathit{e}} \; S\_{\mathit{i}}^{-1} \\ & \stackrel{P\_{\mathit{i},\mathit{i}}\mathit{i}-1}{\rightarrow} & \Phi\_{\mathit{i},\mathit{i}-1} \; \Phi\_{\mathit{i},\mathit{i}-1} \; \Phi\_{\mathit{i},\mathit{i}-1}^{\mathit{T}} - S\_{\mathit{i}} \\ & \begin{array}{rcl} \mathcal{M}\_{\mathit{i},\mathit{i}} & = & \Phi\_{\mathit{i},\mathit{i}-\mathit{i}} \; \mathcal{Z}\_{\mathit{i},\mathit{i}-1} \; \! \! \! \! \! \! \! \! \!} \mathcal{M}\_{\mathit{i},\mathit{i}-1} \\ & \mathcal{Z}\_{\mathit{i},\mathit{i},\mathit{i}-1} & = & \mathcal{M}\_{\mathit{i},\mathit{i}} \; \mathcal{M}\_{\mathit{i},\mathit{i}} \; (\mathcal{M}\_{\mathit{i},\mathit{i}} - \mathcal{S}\_{\mathit{i}}^{-1})^{-1} \mathcal{M}\_{\mathit{i},\mathit{i}} \\ & \mathcal{T}\_{\mathit{i},\mathit{i}} & = & \mathcal{T}\_{\mathit{i},\mathit{i}-1} \; \mathcal{P}\_{\mathit{i},\mathit{i}} \; \! \! \! \! \! \! \end{array} \\ & \text{to to } \mathit{Step 2} \end{array}$$

Figure 7. The three-step recursion of the error variance matrix Pi,t∣<sup>t</sup> <sup>¼</sup> <sup>D</sup> xt � <sup>b</sup>xi,t∣<sup>t</sup> � � for node i. The extra term P<sup>n</sup> <sup>j</sup>¼<sup>1</sup> <sup>l</sup> 2 ij Nj,t would be required to compute Pi,t∣<sup>t</sup>. The entries of the product matrix L kð Þ<sup>n</sup> in Eq. (34) are denoted by lij.

$$\begin{aligned} \mathcal{Z}\_{i,1|0} &= 0, \quad F\_{i,1} = 0 \\ &\downarrow \\ \mathcal{Z}\_{i,1|1} &= nN, \quad \quad \Rightarrow \quad P\_{i,1|1} = a \, \frac{1}{n}N^{-1} \\ &\downarrow \\ \mathcal{Z}\_{i,2|1} &= \mathcal{Z}\_{i,1|1} \quad F\_{i,2} = a \, nN \\ &\vdots & \quad \vdots & \quad \vdots \\ \mathcal{Z}\_{i,t|t-1} &= \mathcal{Z}\_{i,t-1|t-1} \quad F\_{i,t} = a \, (t-1)nN \quad \Rightarrow \quad P\_{i,t|t-1} = P\_{i,t-1|t-1} \\ &\downarrow \\ \mathcal{Z}\_{i,tt} &= tnN, \quad \Rightarrow \quad P\_{i,tt|t} = a \, \frac{1}{tn}N^{-1} \end{aligned} \tag{47}$$

This shows that the consensus-based error variance matrix Pi,t∣<sup>t</sup> is α times its central counterpart Pt∣<sup>t</sup> <sup>¼</sup> ð Þ <sup>1</sup>=tn <sup>N</sup>�<sup>1</sup> . With the vector l≔ li1;…; l ½ � in <sup>T</sup>, application of the Cauchy-Schwarz inequality gives the lower-bound

$$\alpha = \left(e\_n^T e\_n\right) \left(l^T l\right) \ge \left(l^T e\_n\right)^2 = 1\tag{48}$$

as l Ten <sup>¼</sup> 1. Thus scalar <sup>α</sup> is never smaller than 1, i.e. Pi,t∣<sup>t</sup> <sup>≥</sup> Pt∣<sup>t</sup>, showing that the performance of the consensus-based predictor <sup>b</sup>xi,t∣<sup>t</sup> is never better than that of its central version <sup>b</sup>xt∣<sup>t</sup>. The lower-bound Eq. (48) is reached when l ¼ ð Þ 1=n en, i.e. when lij ¼ 1=n (j ¼ 1, …, n). According to Eq. (34), this can be realized if L kð Þ!<sup>n</sup> ð Þ <sup>1</sup>=<sup>n</sup> ene<sup>T</sup> <sup>n</sup> , for which the number of iterations kn might be required to be reasonably large. The conclusion reads therefore that the local filters at nodes i ¼ 1, …, n, generate information matrices Ii,t∣<sup>t</sup>, the inverse of which are different from the actual error variance matrices of the predictors <sup>b</sup>xi,t∣<sup>t</sup>, i.e. <sup>I</sup>�<sup>1</sup> i,t∣<sup>t</sup> 6¼ Pi,t∣<sup>t</sup>.

### 5. Applications to GNSS

Ii, <sup>1</sup>∣<sup>0</sup> ¼ 0, Fi, <sup>1</sup> ¼ 0 ↓

would be required to compute Pi,t∣<sup>t</sup>. The entries of the product matrix L kð Þ<sup>n</sup> in Eq. (34) are denoted by lij.

Figure 7. The three-step recursion of the error variance matrix Pi,t∣<sup>t</sup> <sup>¼</sup> <sup>D</sup> xt � <sup>b</sup>xi,t∣<sup>t</sup>

↓

↓

. With the vector l≔ li1;…; l ½ � in

part Pt∣<sup>t</sup> <sup>¼</sup> ð Þ <sup>1</sup>=tn <sup>N</sup>�<sup>1</sup>

inequality gives the lower-bound

292 Kalman Filters - Theory for Advanced Applications

Ii, <sup>1</sup>∣<sup>1</sup> ¼ nN, ) Pi, <sup>1</sup>∣<sup>1</sup> ¼ α

Ii,2∣<sup>1</sup> ¼ Ii, <sup>1</sup>∣<sup>1</sup>, Fi,<sup>2</sup> ¼ αnN ) Pi, <sup>2</sup>∣<sup>1</sup> ¼ Pi,<sup>1</sup>∣<sup>1</sup>

Ii,t∣<sup>t</sup> ¼ tnN, ) Pi,t∣<sup>t</sup> ¼ α

Ii,t∣t�<sup>1</sup> ¼ Ii,t�1∣t�<sup>1</sup>, Fi,t ¼ αð Þ t � 1 nN ) Pi,t∣t�<sup>1</sup> ¼ Pi,t�1∣t�<sup>1</sup>

This shows that the consensus-based error variance matrix Pi,t∣<sup>t</sup> is α times its central counter-

⋮ ⋮⋮

1 n N�<sup>1</sup>

� � for node i. The extra term P<sup>n</sup>

(47)

<sup>j</sup>¼<sup>1</sup> <sup>l</sup> 2 ij Nj,t

1 tn N�<sup>1</sup>

<sup>T</sup>, application of the Cauchy-Schwarz

The purpose of this section is to demonstrate how the CKF theory, discussed in Section 4, can play a pivotal role in applications for which the GNSS measurements of a network of receivers are to be processed in a real-time manner. In a GNSS network setup, each receiver serves as a sensor node for receiving observables from visible GNSS satellites to determine a range of different parameters such as positions and velocities in an Earth-centered Earth-fixed coordinate system, atmospheric delays, timing and instrumental biases, see e.g. [11, 12]. As the observation equations of the receivers have satellite specific parameters in common, the receivers' observables are often integrated through a computing (fusion) center to provide network-derived parameter solutions that are more precise than their single-receiver versions. Now the idea is to deliver GNSS parameter solutions without the need of having a computing center, such that their precision performance is still comparable to that of network-derived solutions.

As previously discussed, consensus-based algorithms and in particular the CKF can be employed to process network data in a distributed filtering scheme, i.e. no computing center is required. In order to illustrate such applicability, we simulate a network of 13 GNSS receivers located in Perth, Western Australia (Figure 8). As shown in the figure, each node (white circle) represents a receiver having data links (red lines) to its neighbors with interstation distances up to 4 km. We therefore assume that the receivers receive each other data within the ranges not longer than 4 km. For instance, receiver 1 is directly connected to receivers 2 and 6, but not to receiver 3 (the inter-station distance between receivers 1 and 3 is about 8 km).

#### 5.1. GNSS ionospheric observables: Dynamic and measurement models

Although the GNSS observables contain information on various positioning and nonpositioning parameters, here we restrict ourselves to ionospheric observables of the GPS pseudorange measurements only [44]. One should however bear in mind that such restriction is made just for the sake of presentation and illustration of the theory discussed in Sections 3 and 4.

i � s (see Figure 9). They are computed with respect to those of the reference IPP at time

are computed based on the mean Earth's radius 6378:137 km and height of layer 450 km. The

<sup>0</sup>:<sup>02</sup> <sup>þ</sup> sin <sup>θ</sup><sup>s</sup>

σ is set to σ ≈ 65:6 cm as the zenith-referenced standard-deviation of the GPS 'geometry-free'

Suppose that m number of satellites s ¼ 1, …, m, are tracked by the network receivers i ¼ 1, …, n ¼ 13, during the observational campaign. The state-vector sought is structured as

Thus the state-vector xt contains three TEC parameters νo,t, νϕ,t, νψ,t and ð Þ m � 1 between-

3 7

parameters νo,t, νϕ,t, νψ,t is captured by a random-walk process. The corresponding zero-mean

Figure 9. Longitude (ψ) and latitude (ϕ) of an ionospheric piercing point (IPP) corresponding to a receiver-to-satellite

do dϕ dψ

2 6 4 p1 <sup>t</sup> ; b p2 <sup>t</sup> ;…; b pm t

<sup>t</sup> (s 6¼ p). The dynamic model is assumed to be given by (cf. 2, 4 and 21)

<sup>5</sup>, and <sup>b</sup>

<sup>t</sup> are assumed constant in time, while the temporal behavior of the TEC

ps <sup>t</sup> ¼ b ps

h i<sup>T</sup>

xt ¼ νo,t; νϕ,t; νψ,t; b

νo,t�<sup>1</sup> νϕ,t�<sup>1</sup> νψ,t�<sup>1</sup>

2 6 4 � � <sup>¼</sup> <sup>1</sup>:02<sup>2</sup>

i,t are assumed to be mutually uncorrelated with the dispersion (cf. 6)

i,t

i,t denotes the zenith angle of the IPPs. These angles

� � � � <sup>2</sup> <sup>σ</sup><sup>2</sup> (51)

i,t is the satellite elevation angle. The scalar

Consensus-Based Distributed Filtering for GNSS http://dx.doi.org/10.5772/intechopen.71138

<sup>t</sup>�<sup>1</sup> <sup>s</sup> 6¼ <sup>p</sup> (53)

(52)

295

. The angle zs

D ε<sup>s</sup> i,t

forming the variance matrices Rt in Eq. (6), where θ<sup>s</sup>

instance t, i.e. ψo,t and ϕo,t

pseudo-range measurements [48].

ps

νo,t νϕ,t νψ,t

2 6 4

ps

line-of-sight. The distance scales are exaggerated.

measurement noises ε<sup>s</sup>

satellite DCBs b

Thus the DCBs b

Figure 8. A network of 13 GNSS receivers simulated over Perth, Western Australia. Each node (white circle) represents a receiver tracking GNSS satellites. The receivers have data links to their neighbors with inter-station distances up to 4 km. The data links are shown by red lines.

Would one, for instance, make use of the very precise carrier-phase measurements and/or formulate a multi-GNSS measurement setup, solutions of higher precision levels are therefore expected.

Let the scalar ys i,t denote the pseudo-range ionospheric observable that the receiver i collects from satellite s at time instance t. The corresponding measurement model, formed by the between-satellite differences y ps i,t ≔ys i,t � y p i,t ð Þ s 6¼ p , reads (cf. 5)

$$y\_{i,t}^{ps} = \left\{ a\_{i,t;o}^{ps} \nu\_{o,t} + a\_{i,t;\phi}^{ps} \nu\_{\phi,t} + a\_{i,t;\psi}^{ps} \nu\_{\psi,t} \right\} - b\_t^{ps} + \varepsilon\_{i,t}^{ps} \tag{49}$$

where the term within f g: refers to the first-order slant ionospheric delays, and b ps <sup>t</sup> denotes the between-satellite differential code biases (DCBs). We use a regional single-layer model [45, 46] to represent the slant ionospheric delays in terms of 1) νo,t as the vertical total electron content (TEC), 2) νϕ,t and 3) νψ,t as the south-to-north and west-to-east spatial gradient of νo,t, respectively. The corresponding known coefficients follow from Ref. [47]

$$a\_{i,t;o}^{s} = \frac{1}{\cos\left(z\_{i,t}^{s}\right)}, \quad d\_{i,t;\phi}^{s} = \frac{1}{\cos\left(z\_{i,t}^{s}\right)} \left(\phi\_{i,t}^{s} - \phi\_{o,t}\right), \quad a\_{i,t;\psi}^{s} = \frac{1}{\cos\left(z\_{i,t}^{s}\right)} \cos\left(\phi\_{i,t}^{s}\right) \left(\psi\_{i,t}^{s} - \psi\_{o,t}\right) \tag{50}$$

with ð Þ: ps≔ð Þ: <sup>s</sup> � ð Þ: <sup>p</sup> . The angles ψ<sup>s</sup> i,t and <sup>ϕ</sup><sup>s</sup> i,t , respectively, denote the longitude and latitude of the ionospheric piercing points (IPPs) corresponding to the receiver-to-satellite line-of-sight i � s (see Figure 9). They are computed with respect to those of the reference IPP at time instance t, i.e. ψo,t and ϕo,t . The angle zs i,t denotes the zenith angle of the IPPs. These angles are computed based on the mean Earth's radius 6378:137 km and height of layer 450 km. The measurement noises ε<sup>s</sup> i,t are assumed to be mutually uncorrelated with the dispersion (cf. 6)

$$\mathcal{D}\left(\varepsilon\_{i,t}^{s}\right) = \frac{1.02^2}{\left(0.02 + \sin\left(\theta\_{i,t}^{s}\right)\right)^2}\sigma^2 \tag{51}$$

forming the variance matrices Rt in Eq. (6), where θ<sup>s</sup> i,t is the satellite elevation angle. The scalar σ is set to σ ≈ 65:6 cm as the zenith-referenced standard-deviation of the GPS 'geometry-free' pseudo-range measurements [48].

Suppose that m number of satellites s ¼ 1, …, m, are tracked by the network receivers i ¼ 1, …, n ¼ 13, during the observational campaign. The state-vector sought is structured as

$$\mathbf{x}\_t = \begin{bmatrix} \nu\_{o,t}, \nu\_{\phi,t}, \nu\_{\psi,t}, \boldsymbol{b}\_t^{p1}, \boldsymbol{b}\_t^{p2}, \dots, \boldsymbol{b}\_t^{pm} \end{bmatrix}^T \tag{52}$$

Thus the state-vector xt contains three TEC parameters νo,t, νϕ,t, νψ,t and ð Þ m � 1 betweensatellite DCBs b ps <sup>t</sup> (s 6¼ p). The dynamic model is assumed to be given by (cf. 2, 4 and 21)

Would one, for instance, make use of the very precise carrier-phase measurements and/or formulate a multi-GNSS measurement setup, solutions of higher precision levels are therefore

Figure 8. A network of 13 GNSS receivers simulated over Perth, Western Australia. Each node (white circle) represents a receiver tracking GNSS satellites. The receivers have data links to their neighbors with inter-station distances up to 4 km.

from satellite s at time instance t. The corresponding measurement model, formed by the

i,t;<sup>ϕ</sup> νϕ,t þ a

between-satellite differential code biases (DCBs). We use a regional single-layer model [45, 46] to represent the slant ionospheric delays in terms of 1) νo,t as the vertical total electron content (TEC), 2) νϕ,t and 3) νψ,t as the south-to-north and west-to-east spatial gradient of νo,t, respec-

> i,t � ϕo,t � �

the ionospheric piercing points (IPPs) corresponding to the receiver-to-satellite line-of-sight

n o

ps

where the term within f g: refers to the first-order slant ionospheric delays, and b

i,t and <sup>ϕ</sup><sup>s</sup> i,t

i,t denote the pseudo-range ionospheric observable that the receiver i collects

i,t ð Þ s 6¼ p , reads (cf. 5)

ps i,t;<sup>ψ</sup> νψ,t

, as

� b ps <sup>t</sup> þ ε ps

i,t;<sup>ψ</sup> <sup>¼</sup> <sup>1</sup>

cos zs i,t � � cos <sup>ϕ</sup><sup>s</sup>

, respectively, denote the longitude and latitude of

i,t (49)

ps

i,t � �

ψs i,t � ψo,t � �

<sup>t</sup> denotes the

(50)

expected.

as

i,t; <sup>o</sup> <sup>¼</sup> <sup>1</sup>

cos zs i,t � � , as

with ð Þ: ps≔ð Þ: <sup>s</sup> � ð Þ: <sup>p</sup>

Let the scalar ys

between-satellite differences y

The data links are shown by red lines.

294 Kalman Filters - Theory for Advanced Applications

ps i,t ≔ys i,t � y p

tively. The corresponding known coefficients follow from Ref. [47]

cos zs i,t � � <sup>ϕ</sup><sup>s</sup>

i,t;<sup>ϕ</sup> <sup>¼</sup> <sup>1</sup>

. The angles ψ<sup>s</sup>

y ps i,t ¼ a ps i,t; <sup>o</sup> νo,t þ a

$$
\begin{bmatrix} \boldsymbol{\nu}\_{\boldsymbol{\phi},t} \\ \boldsymbol{\nu}\_{\boldsymbol{\phi},t} \\ \boldsymbol{\nu}\_{\boldsymbol{\psi},t} \end{bmatrix} = \begin{bmatrix} \boldsymbol{\nu}\_{\boldsymbol{\phi},t-1} \\ \boldsymbol{\nu}\_{\boldsymbol{\phi},t-1} \\ \boldsymbol{\nu}\_{\boldsymbol{\psi},t-1} \end{bmatrix} + \begin{bmatrix} \boldsymbol{d}\_{\boldsymbol{\phi}} \\ \boldsymbol{d}\_{\boldsymbol{\phi}} \\ \boldsymbol{d}\_{\boldsymbol{\psi}} \end{bmatrix}, \quad \text{and} \quad \boldsymbol{b}\_{t}^{\rm ps} = \boldsymbol{b}\_{t-1}^{\rm ps} \quad s \neq p \tag{53}
$$

Thus the DCBs b ps <sup>t</sup> are assumed constant in time, while the temporal behavior of the TEC parameters νo,t, νϕ,t, νψ,t is captured by a random-walk process. The corresponding zero-mean

Figure 9. Longitude (ψ) and latitude (ϕ) of an ionospheric piercing point (IPP) corresponding to a receiver-to-satellite line-of-sight. The distance scales are exaggerated.

process noises are assumed to be mutually uncorrelated, having the standard-deviations <sup>σ</sup>do <sup>¼</sup> 1 mm/ ffiffiffiffiffiffi <sup>s</sup>ec <sup>p</sup> and <sup>σ</sup><sup>d</sup><sup>ϕ</sup> <sup>¼</sup> <sup>σ</sup>dpsi <sup>¼</sup> 5 mm/rad/ ffiffiffiffiffiffi sec p [49].

5.3. Central (network-based) versus local (single-receiver) solutions

of the local TEC solutions are ffiffiffiffiffi

5.4. Role of CKF in improving local solutions

Before discussing the precision performance of the CKF solutions, we first compare the networkbased (central) TEC solutions with the solutions that are obtained by the data of one singlereceiver only (referred to as the local solutions). At the filter initialization, the standard-deviations

square-root of the number of nodes). This is because of the fact that each of the 13 network receivers independently provides equally precise solutions. In that case, the central solution follows then by averaging all the 13 local solutions. Due to the common dynamic model Eq. (53) however, the local solutions become correlated over time. After the filter initialization, the central solution would therefore not follow the average of its local versions. The standard-deviation results, after one hour of the filter initialization, are presented in Figure 11. Only the results of the receiver 1 are shown as local solutions (in red). As shown, the standard-deviations get stable over time as the filters reach their steady-state. On the right panel of the figure, the local-to-central standard-deviation ratios are also presented. In case of the vertical TECs νo,t, the ratios vary from 1.5 to 3. For the horizontal gradients νϕ,t and νψ,t, the ratios are about 2 and 2.5, respectively.

With the results of Figure 11, we observed that the central TEC solutions considerably outperform their local counterparts in the sense of delivering more precise outcomes, i.e. the local-to-central standard-deviation ratios are considerably larger than 1. We now employ the

Figure 11. Left: Standard-deviation of the central (green) and local (red) solutions of the TEC parameters νo,t (top), νϕ,t (middle) and νψ,t (bottom) as functions of time. Right: The corresponding local-to-central standard-deviation ratios.

<sup>13</sup> <sup>p</sup> <sup>≈</sup> <sup>3</sup>:6 times larger than those of the central TEC solutions (i.e.

Consensus-Based Distributed Filtering for GNSS http://dx.doi.org/10.5772/intechopen.71138 297

#### 5.2. Observational campaign

The network receivers i ¼ 1, …, n (n ¼ 13), shown in Figure 8, are assumed to track GPS satellites over 16 hours from 8:00 to 24:00 Perth local time, on 02-06-2016. The observation sampling rate is set to Δ ¼ 1 minute. Thus the number of observational epochs (time instances) is 960. As to the data sending rate δ (cf. 5), we assume three different sending rates δ ¼ 5, 10 and 15 seconds. Thus the number of iterative communications between the neighboring receivers takes the values kn ¼ 4, 6 and 12. The consensus protocol 3 (Table 1) is applied to the CKF of each receiver.

As the satellites revolve around the Earth, not all of which are simultaneously visible to the 'small-scale' network of Figure 8. Their visibility over time is shown in Figure 10 (left panel) in which the satellites with elevation angles smaller than 10 degrees are excluded. There are 31 GPS satellites (i.e. m ¼ 31), with PRN 4 absent (PRN refers to the satellite identifier). PRN 22 has the maximum duration of visibility, while PRN 21 has the minimum duration of visibility. Note also that PRNs 2, 6, 16, 17, 19, 26 and 32 disappear (set) and reappear (re-rise). That is why their visibility is shown via two separate time intervals. Figure 10 (right panel) shows the trajectories of the ionospheric pierce points on the ionospheric single layer that are made by receiver-to-satellite line-of-sight paths. It is the spatial distribution of these points that drives the coefficients a ps i,t; <sup>o</sup>, a ps i,t;ϕ, a ps i,t;<sup>ψ</sup> in Eq. (49).

In the following we present precision analyses on the measurement update solutions of xt in Eq. (52), given the network and satellite configurations shown in Figures 8 and 10, respectively. Throughout the text, PRN 10 is chosen as the pivot satellite p (cf. (49)). By the term 'standard-deviation', we mean the square-root of prediction errors' variance.

Figure 10. Left: GPS satellites visibility over time, viewed from Perth, Western Australia. Right: Trajectories of the corresponding IPPs made by receiver-to-satellite line-of-sight paths. The satellites are indicated by different colors.

#### 5.3. Central (network-based) versus local (single-receiver) solutions

Before discussing the precision performance of the CKF solutions, we first compare the networkbased (central) TEC solutions with the solutions that are obtained by the data of one singlereceiver only (referred to as the local solutions). At the filter initialization, the standard-deviations of the local TEC solutions are ffiffiffiffiffi <sup>13</sup> <sup>p</sup> <sup>≈</sup> <sup>3</sup>:6 times larger than those of the central TEC solutions (i.e. square-root of the number of nodes). This is because of the fact that each of the 13 network receivers independently provides equally precise solutions. In that case, the central solution follows then by averaging all the 13 local solutions. Due to the common dynamic model Eq. (53) however, the local solutions become correlated over time. After the filter initialization, the central solution would therefore not follow the average of its local versions. The standard-deviation results, after one hour of the filter initialization, are presented in Figure 11. Only the results of the receiver 1 are shown as local solutions (in red). As shown, the standard-deviations get stable over time as the filters reach their steady-state. On the right panel of the figure, the local-to-central standard-deviation ratios are also presented. In case of the vertical TECs νo,t, the ratios vary from 1.5 to 3. For the horizontal gradients νϕ,t and νψ,t, the ratios are about 2 and 2.5, respectively.

#### 5.4. Role of CKF in improving local solutions

process noises are assumed to be mutually uncorrelated, having the standard-deviations

The network receivers i ¼ 1, …, n (n ¼ 13), shown in Figure 8, are assumed to track GPS satellites over 16 hours from 8:00 to 24:00 Perth local time, on 02-06-2016. The observation sampling rate is set to Δ ¼ 1 minute. Thus the number of observational epochs (time instances) is 960. As to the data sending rate δ (cf. 5), we assume three different sending rates δ ¼ 5, 10 and 15 seconds. Thus the number of iterative communications between the neighboring receivers takes the values kn ¼ 4, 6 and 12. The consensus protocol 3 (Table 1) is applied to the CKF of each receiver. As the satellites revolve around the Earth, not all of which are simultaneously visible to the 'small-scale' network of Figure 8. Their visibility over time is shown in Figure 10 (left panel) in which the satellites with elevation angles smaller than 10 degrees are excluded. There are 31 GPS satellites (i.e. m ¼ 31), with PRN 4 absent (PRN refers to the satellite identifier). PRN 22 has the maximum duration of visibility, while PRN 21 has the minimum duration of visibility. Note also that PRNs 2, 6, 16, 17, 19, 26 and 32 disappear (set) and reappear (re-rise). That is why their visibility is shown via two separate time intervals. Figure 10 (right panel) shows the trajectories of the ionospheric pierce points on the ionospheric single layer that are made by receiver-to-satellite line-of-sight paths. It is the spatial distribution of these points that drives

In the following we present precision analyses on the measurement update solutions of xt in Eq. (52), given the network and satellite configurations shown in Figures 8 and 10, respectively. Throughout the text, PRN 10 is chosen as the pivot satellite p (cf. (49)). By the term

Figure 10. Left: GPS satellites visibility over time, viewed from Perth, Western Australia. Right: Trajectories of the corresponding IPPs made by receiver-to-satellite line-of-sight paths. The satellites are indicated by different colors.

sec p [49].

<sup>s</sup>ec <sup>p</sup> and <sup>σ</sup><sup>d</sup><sup>ϕ</sup> <sup>¼</sup> <sup>σ</sup>dpsi <sup>¼</sup> 5 mm/rad/ ffiffiffiffiffiffi

i,t;<sup>ψ</sup> in Eq. (49).

'standard-deviation', we mean the square-root of prediction errors' variance.

<sup>σ</sup>do <sup>¼</sup> 1 mm/ ffiffiffiffiffiffi

the coefficients a

ps i,t; <sup>o</sup>, a ps i,t;ϕ, a ps

5.2. Observational campaign

296 Kalman Filters - Theory for Advanced Applications

With the results of Figure 11, we observed that the central TEC solutions considerably outperform their local counterparts in the sense of delivering more precise outcomes, i.e. the local-to-central standard-deviation ratios are considerably larger than 1. We now employ the

Figure 11. Left: Standard-deviation of the central (green) and local (red) solutions of the TEC parameters νo,t (top), νϕ,t (middle) and νψ,t (bottom) as functions of time. Right: The corresponding local-to-central standard-deviation ratios.

CKF for each node (receiver) i ¼ 1, …, 13, to improve the local solutions' precision performance via consensus-based iterative communications between the receivers. In doing so, we make use of the three-step recursion given in Figure 7 to evaluate the error variance matrices Pi,t∣<sup>t</sup> (i ¼ 1, …, 13), thereby computing the CKF-to-central standard-deviation ratios. The stated ratios are presented in Figure 12 for two different data sending rates δ ¼ 15 seconds (left panel) and δ ¼ 5 seconds (right panel). In both cases, the CKF-to-central standard-deviation ratios are smaller than their local-to-central versions shown in Figure 11 (right panel), illustrating that employing the CKF does indeed improve the local solutions' precision. Since more iterative communications take place for δ ¼ 5, the corresponding ratios are very close to 1. In that case, the CKF of each receiver is expected to have a similar precision performance to that of the central (network-based) filter. For the case δ ¼ 15 however, the CKF performance of each receiver does very much depend on the number of the receiver's neighbors. This is because of the fact that only 4 iterative communications between the receivers take place (i.e. kn ¼ 4). The receivers with the minimum number of neighbors, i.e. receivers 1, 3 and 13 (Figure 8), have the worst precision performance as the corresponding ratios take largest values. On the other hand, the receivers with the maximum number of neighbors, i.e. receivers 4, 7, 9 and 8, have the best performance as the corresponding ratios are close to 1.

visible, the smaller the standard-deviation is expected. We now consider the required time to have between-satellite DCBs solutions with standard-deviation smaller than 0.5 nanoseconds. Because of the stated difference in the standard-deviations, each between-satellite DCB corresponds to a different required time. For the central filter, the minimum value of such required time is 7 minutes, with the 25th percentile as 12, median as 38, 75th percentile as 63 and the maximum as 84 minutes. Thus after 84 minutes of the filter initialization, all central DCB solutions have standard-deviations smaller than 0.5 nanoseconds. Such percentiles can be represented by a 'boxplot'. We compute the stated percentiles for all the CKF solutions and compare their boxplots with the central one in Figure 13. The results are presented for three

Consensus-Based Distributed Filtering for GNSS http://dx.doi.org/10.5772/intechopen.71138 299

Figure 13. Boxplots of the required time (minutes) to have between-satellite DCBs solutions with standard-deviation smaller than 0.5 nanoseconds. The performance of the CKF of each node (receiver) is compared to that of the central filter (Cen.). In each boxplot, the horizontal lines from bottom to top show the minimum (black), 25th percentile (blue), median (green), 75th percentile (blue) and maximum (black) of the stated required time. The data sending rate is set to Top: δ ¼ 15 seconds (i.E. 4 iterative communications), Middle: δ ¼ 10 seconds (i.E. 6 iterative communications), and Bottom: δ ¼ 5

seconds (i.E. 12 iterative communications).

Next to the solutions of the TEC parameters νo,t, νϕ,t and νψ,t, we also analyze CKF solutions of the between-satellite DCBs b ps <sup>t</sup> (s 6¼ p) in Eq. (52). Because of the difference in the satellites visibility over time (cf. Figure 10), the DCBs' standard-deviations are quite distinct and very much depend on the duration of the satellites visibility. The more a pair of satellites p � s are

Figure 12. Time-series of the CKF-to-central standard-deviation ratios corresponding to the TEC parameters νo,t (top), νϕ,t (middle) and νψ,t (bottom). Left: The data sending rate is set to δ ¼ 15 seconds (i.E. 4 iterative communications). Right: The data sending rate is set to δ ¼ 5 seconds (i.E. 12 iterative communications). The results of the nodes (receivers) are indicated by different colors.

visible, the smaller the standard-deviation is expected. We now consider the required time to have between-satellite DCBs solutions with standard-deviation smaller than 0.5 nanoseconds. Because of the stated difference in the standard-deviations, each between-satellite DCB corresponds to a different required time. For the central filter, the minimum value of such required time is 7 minutes, with the 25th percentile as 12, median as 38, 75th percentile as 63 and the maximum as 84 minutes. Thus after 84 minutes of the filter initialization, all central DCB solutions have standard-deviations smaller than 0.5 nanoseconds. Such percentiles can be represented by a 'boxplot'. We compute the stated percentiles for all the CKF solutions and compare their boxplots with the central one in Figure 13. The results are presented for three

CKF for each node (receiver) i ¼ 1, …, 13, to improve the local solutions' precision performance via consensus-based iterative communications between the receivers. In doing so, we make use of the three-step recursion given in Figure 7 to evaluate the error variance matrices Pi,t∣<sup>t</sup> (i ¼ 1, …, 13), thereby computing the CKF-to-central standard-deviation ratios. The stated ratios are presented in Figure 12 for two different data sending rates δ ¼ 15 seconds (left panel) and δ ¼ 5 seconds (right panel). In both cases, the CKF-to-central standard-deviation ratios are smaller than their local-to-central versions shown in Figure 11 (right panel), illustrating that employing the CKF does indeed improve the local solutions' precision. Since more iterative communications take place for δ ¼ 5, the corresponding ratios are very close to 1. In that case, the CKF of each receiver is expected to have a similar precision performance to that of the central (network-based) filter. For the case δ ¼ 15 however, the CKF performance of each receiver does very much depend on the number of the receiver's neighbors. This is because of the fact that only 4 iterative communications between the receivers take place (i.e. kn ¼ 4). The receivers with the minimum number of neighbors, i.e. receivers 1, 3 and 13 (Figure 8), have the worst precision performance as the corresponding ratios take largest values. On the other hand, the receivers with the maximum number of neighbors, i.e. receivers 4, 7, 9 and 8, have

Next to the solutions of the TEC parameters νo,t, νϕ,t and νψ,t, we also analyze CKF solutions of

visibility over time (cf. Figure 10), the DCBs' standard-deviations are quite distinct and very much depend on the duration of the satellites visibility. The more a pair of satellites p � s are

Figure 12. Time-series of the CKF-to-central standard-deviation ratios corresponding to the TEC parameters νo,t (top), νϕ,t (middle) and νψ,t (bottom). Left: The data sending rate is set to δ ¼ 15 seconds (i.E. 4 iterative communications). Right: The data sending rate is set to δ ¼ 5 seconds (i.E. 12 iterative communications). The results of the nodes (receivers) are

<sup>t</sup> (s 6¼ p) in Eq. (52). Because of the difference in the satellites

the best performance as the corresponding ratios are close to 1.

ps

the between-satellite DCBs b

298 Kalman Filters - Theory for Advanced Applications

indicated by different colors.

Figure 13. Boxplots of the required time (minutes) to have between-satellite DCBs solutions with standard-deviation smaller than 0.5 nanoseconds. The performance of the CKF of each node (receiver) is compared to that of the central filter (Cen.). In each boxplot, the horizontal lines from bottom to top show the minimum (black), 25th percentile (blue), median (green), 75th percentile (blue) and maximum (black) of the stated required time. The data sending rate is set to Top: δ ¼ 15 seconds (i.E. 4 iterative communications), Middle: δ ¼ 10 seconds (i.E. 6 iterative communications), and Bottom: δ ¼ 5 seconds (i.E. 12 iterative communications).

different data sending rates δ ¼ 15 seconds (top), δ ¼ 10 seconds (middle) and δ ¼ 5 seconds (bottom). As shown, the more the number of iterative communications, the more similar the boxplots becomes, i.e. the nodes (receivers) are reaching consensus. Similar to the TEC solutions, the DCB precision performance of the CKF corresponding to the receivers 4, 7, 9 and 8 is almost similar to that of the central one, irrespective of the number of iterative communications. This follows from the fact that the stated receivers have the maximum number of neighbors (Figure 8), thus efficiently approximating the averages Nt and rt in Eq. (28) after a few iterations. On the other hand, the receivers with the minimum number of neighbors require more number of iterative communications in order for their CKF precision performance to get similar to that of the central filter.

Author details

Amir Khodabandeh<sup>1</sup>

References

\*, Peter J.G. Teunissen1,2 and Safoora Zaminpardaz1

Consensus-Based Distributed Filtering for GNSS http://dx.doi.org/10.5772/intechopen.71138 301

[1] Cattivelli FS, Sayed AH. Diffusion strategies for distributed Kalman filtering and smooth-

[2] Das S, Moura JMF. Consensus+innovations distributed Kalman filter with optimized

[3] Jadbabaie A, Lin J, Stephen Morse A. Coordination of groups of mobile autonomous agents using nearest neighbor rules. IEEE Transactions on Automatic Control. 2003;48(6):988-1001

[4] Kingston DB, Beard RW. Discrete-time average-consensus under switching network topologies. In: American Control Conference, 2006. IEEE; 2006. pp. 3551-3556

[5] Moreau L. Stability of multiagent systems with time-dependent communication links.

[6] Olfati-Saber R, Murray RM. Consensus problems in networks of agents with switching topology and time-delays. IEEE Transactions on Automatic Control. 2004;49(9):1520-1533

[7] Scherber DS, Papadopoulos HC. Locally constructed algorithms for distributed computations in ad-hoc networks. In: Proceedings of the 3rd International Symposium on

[8] Xiao L, Boyd S, Lall S. A scheme for robust distributed sensor fusion based on average consensus. In: Proceedings of the 4th International Symposium on Information Proces-

[9] Rigatos GG. Distributed filtering over sensor networks for autonomous navigation of

[10] Sugar TG, Kumar V. Control of cooperating mobile manipulators. IEEE Transactions on

[11] Hofmann-Wellenhof B, Lichtenegger H, Wasle E. GNSS: Global Navigation Satellite

[12] Teunissen PJG, Montenbruck O, editors. Springer Handbook of Global Navigation Satel-

Systems: GPS, Glonass, Galileo, and More. New York: Springer; 2008

ing. IEEE Transactions on Automatic Control. 2010;55(9):2069-2084

gains. IEEE Transactions on Signal Processing. 2017;65(2):467-481

IEEE Transactions on Automatic Control. 2005;50(2):169-182

Information Processing in Sensor Networks. ACM; 2004. pp. 11-19

sing in Sensor Networks. IEEE Press; 2005. pp. 63-70

UAVs. Intelligent Service Robotics. 2012;5(3):179-198

Robotics and Automation. 2002;18(1):94-103

lite Systems. Switzerland: Springer; 2017

\*Address all correspondence to: amir.khodabandeh@curtin.edu.au

1 Curtin University of Technology, Australia

2 Delft University of Technology, The Netherlands

### 6. Concluding remarks and future outlook

In this contribution we reviewed Kalman filtering in its information form and showed how the additive measurement update (28) can be realized by employing average-consensus rules, even when not all nodes are directly connected, thus allowing the sensor nodes to develop their own distributed filters. The nodes are assumed linked to each other at least through a 'path' so that information can flow from each node to all other nodes. Under this assumption, average-consensus protocols can deliver consensus states ½ � Ni,tð Þ kn ;ri,tð Þ kn as an approximation of the averages Nt;rt in Eq. (28) at every time instance <sup>t</sup> <sup>¼</sup> <sup>1</sup>, <sup>2</sup>, …, thus allowing one to establish a CKF recursion at every node i ¼ 1, …, n. To improve the stated approximation, the neighboring nodes have to establish a number of iterative data communications to transfer and receive their consensus states. This makes the CKF implementation applicable only for the applications in which the state-vectors change rather slowly over time (i.e. the sampling rate Δ can take large values) and/or for the cases where the sensor nodes transfer their data rather quickly (i.e. the sending rate δ can take small values).

We developed a three-step recursion of the CKF error variance matrix (Figure 7). This recursion conveys useful information about the precision performance of the local filters i ¼ 1, …, n, thereby enabling one to a-priori design and analyze sensor networks with different numbers of iterative communications. As an illustrative example, we applied the stated recursion to a small-scale network of GNSS receivers and showed the role taken by the CKF in improving the precision of the solutions at each single receiver. In near future the proliferation of low-cost receivers will give rise to an increase in the number of GNSS users. Employing the CKF or other distributed filtering techniques, GNSS users can therefore potentially deliver highprecision parameter solutions without the need of having a computing center.

### Acknowledgements

The second author is the recipient of an Australian Research Council Federation Fellowship (project number FF0883188). This support is gratefully acknowledged.

### Author details

different data sending rates δ ¼ 15 seconds (top), δ ¼ 10 seconds (middle) and δ ¼ 5 seconds (bottom). As shown, the more the number of iterative communications, the more similar the boxplots becomes, i.e. the nodes (receivers) are reaching consensus. Similar to the TEC solutions, the DCB precision performance of the CKF corresponding to the receivers 4, 7, 9 and 8 is almost similar to that of the central one, irrespective of the number of iterative communications. This follows from the fact that the stated receivers have the maximum number of neighbors (Figure 8), thus efficiently approximating the averages Nt and rt in Eq. (28) after a few iterations. On the other hand, the receivers with the minimum number of neighbors require more number of iterative communications in order for their CKF precision performance to get

In this contribution we reviewed Kalman filtering in its information form and showed how the additive measurement update (28) can be realized by employing average-consensus rules, even when not all nodes are directly connected, thus allowing the sensor nodes to develop their own distributed filters. The nodes are assumed linked to each other at least through a 'path' so that information can flow from each node to all other nodes. Under this assumption, average-consensus protocols can deliver consensus states ½ � Ni,tð Þ kn ;ri,tð Þ kn as an approxima-

establish a CKF recursion at every node i ¼ 1, …, n. To improve the stated approximation, the neighboring nodes have to establish a number of iterative data communications to transfer and receive their consensus states. This makes the CKF implementation applicable only for the applications in which the state-vectors change rather slowly over time (i.e. the sampling rate Δ can take large values) and/or for the cases where the sensor nodes transfer their data rather

We developed a three-step recursion of the CKF error variance matrix (Figure 7). This recursion conveys useful information about the precision performance of the local filters i ¼ 1, …, n, thereby enabling one to a-priori design and analyze sensor networks with different numbers of iterative communications. As an illustrative example, we applied the stated recursion to a small-scale network of GNSS receivers and showed the role taken by the CKF in improving the precision of the solutions at each single receiver. In near future the proliferation of low-cost receivers will give rise to an increase in the number of GNSS users. Employing the CKF or other distributed filtering techniques, GNSS users can therefore potentially deliver high-

The second author is the recipient of an Australian Research Council Federation Fellowship

precision parameter solutions without the need of having a computing center.

(project number FF0883188). This support is gratefully acknowledged.

in Eq. (28) at every time instance <sup>t</sup> <sup>¼</sup> <sup>1</sup>, <sup>2</sup>, …, thus allowing one to

similar to that of the central filter.

300 Kalman Filters - Theory for Advanced Applications

tion of the averages Nt;rt

Acknowledgements

6. Concluding remarks and future outlook

quickly (i.e. the sending rate δ can take small values).

Amir Khodabandeh<sup>1</sup> \*, Peter J.G. Teunissen1,2 and Safoora Zaminpardaz1

\*Address all correspondence to: amir.khodabandeh@curtin.edu.au


### References


[13] Casbeer DW, Beard R. Distributed information filtering using consensus filters. In: American Control Conference, 2009. ACC'09. IEEE; 2009. pp. 1882-1887

[31] Bar-Shalom Y, Li XR. Estimation and Tracking- Principles, Techniques, and Software. Vol.

Consensus-Based Distributed Filtering for GNSS http://dx.doi.org/10.5772/intechopen.71138 303

[32] Grewal MS, Andrews AP. Kalman Filtering; Theory and Practice Using MATLAB. 3rd ed.

[34] Maybeck PS. Stochastic Models, Estimation, and Control. Vol. 1. Academic Press; 1979.

[35] Simon D. Optimal State Estimation: Kalman, H [Infinity] and Nonlinear Approaches.

[36] Sorenson HW. Kalman filtering techniques. In: Leondes CT editor. Advances in Control

[37] Stark H, Woods JW. Probability, Random Processes, and Estimation Theory for Engi-

[38] Teunissen PJG, Khodabandeh A. BLUE, BLUP and the Kalman filter: Some new results.

[39] Khodabandeh A, Teunissen PJG. A recursive linear MMSE filter for dynamic systems with unknown state vector means. GEM - International Journal on Geomathematics. 2014;

[40] Ren W, Beard RW. Consensus of information under dynamically changing interaction topologies. In: American Control Conference, 2004. Proceedings of the 2004, Vol. 6. IEEE;

[41] Henderson HV, Pukelsheim F, Searle SR. On the history of the Kronecker product. Linear

[43] Wolfowitz J. Products of indecomposable, aperiodic, stochastic matrices. Proceedings of

[44] Blewitt G. An automatic editing algorithm for GPS data. Geophysical Research Letters.

[45] Mannucci AJ, Wilson BD, Yuan DN, Ho CH, Lindqwister UJ, Runge TF. A global mapping technique for GPS-derived ionospheric total electron content measurements. Radio

[46] Schaer S. Mapping and predicting the Earth's ionosphere using the global positioning

[47] Brunini C, Azpilicueta FJ. Accuracy assessment of the GPS-based slant total electron

[42] Horn RA, Johnson CR. Matrix Analysis. New York: Cambridge UP; 1985

system. PhD thesis. Bern, Switzerland: University of Bern; 1999

content. Journal of Geodesy. 2009;83(8):773-785

the American Mathematical Society. 1963;14(5):733-737

[33] Kailath T, Sayed AH, Hassibi B. Linear Estimation. New Jersey: Prentice-Hall; 2000

1993. Norwood, MA: Artech House, Inc; 1993

New Jersey: John Wiley and Sons; 2008

New Jersey: John Wiley and Sons; 2006

Journal of Geodesy. 2013;87(5):461-473

and Multilinear Algebra. 1983;14(2):113-120

Systems Theory and Applications, Vol. 3; 1966. pp. 219-292

neers. Englewood Cliffs, New Jersey: Prentice-Hall; 1986

Republished 1994

5(1):17-31

2004. pp. 4939-4944

1990;17(3):199-202

Science. 1998;33(3):565-582


[13] Casbeer DW, Beard R. Distributed information filtering using consensus filters. In: Amer-

[14] Khodabandeh A, Teunissen PJG. An analytical study of PPP-RTK corrections: Precision,

[15] Li W, Nadarajah N, Teunissen PJG, Khodabandeh A, Chai Y. Array-aided single-frequency state-space RTK with combined GPS, Galileo, IRNSS, and QZSS L5/E5a observations.

[16] Li X, Ge M, Dai X, Ren X, Fritsche M, Wickert J, Schuh H. Accuracy and reliability of multi-GNSS real-time precise positioning: GPS, GLONASS, BeiDou, and Galileo. Journal

[17] Odolinski R, Teunissen PJG. Single-frequency, dual-GNSS versus dual-frequency, single-GNSS: A low-cost and high-grade receivers GPS-BDS RTK analysis. Journal of Geodesy.

[18] Zaminpardaz S, Teunissen PJG, Nadarajah N. GLONASS CDMA L3 Ambiguity Resolu-

[20] Candy JV. Signal Processing: Model Based Approach. New Jersey: McGraw-Hill, Inc.; 1986

[22] Gibbs BP. Advanced Kalman Filtering, Least-Squares and Modeling: A Practical Hand-

[23] Jazwinski AH. Stochastic Processes and Filtering Theory. Maryland: Dover Publications;

[24] Kailath T. Lectures on Wiener and Kalman Filtering. Number 140. Vienna: Springer; 1981 [25] Kalman RE. A new approach to linear filtering and prediction problems. Journal of Basic

[26] Teunissen PJG. Best prediction in linear models with mixed integer/real unknowns:

[27] Bode HW, Shannon CE. A simplified derivation of linear least square smoothing and

[28] Kailath T. An innovations approach to least-squares estimation–part I: Linear filtering in additive white noise. Automatic Control, IEEE Transactions on. 1968;13(6):646-655 [29] Zadeh LA, Ragazzini JR. An extension of Wiener's theory of prediction. Journal of

[30] Anderson BDO, Moore JB. Optimal Filtering. Vol. 11. New Jersey: Prentice-hall Engle-

tion and Positioning. Berlin, Heidelberg: GPS Solutions; 2016. pp. 1-15

[19] Brammer K, Siffling G. Kalman-Bucy Filters. Berlin: Artech House; 1989

Theory and application. Journal of Geodesy. 2007;81(12):759-780

prediction theory. Proceedings of the IRE. 1950;38(4):417-425

[21] Gelb A. Applied Optimal Estimation. London: MIT press; 1974

ican Control Conference, 2009. ACC'09. IEEE; 2009. pp. 1882-1887

Journal of Surveying Engineering. 2017;143(4):04017006

of Geodesy. 2015;89(6):607-635

302 Kalman Filters - Theory for Advanced Applications

book. New Jersey: Wiley; 2011

Engineering. 1960;82(1):35-45

Applied Physics. 1950;21(7):645-655

wood Cliffs; 1979

1991

2016;90(11):1255-1278

correlation and user-impact. Journal of Geodesy. 2015;89(11):1109-1132


[48] Khodabandeh A, Teunissen PJG. Array-aided multifrequency GNSS Ionospheric sensing: Estimability and precision analysis. IEEE Transactions on Geoscience and Remote Sens-

[49] Julien O, Macabiau C, Issler JL. Ionospheric delay estimation strategies using Galileo E5 signals only. In: GNSS 2009, 22nd International Technical Meeting of The Satellite Divi-

sion of the Institute of Navigation, Savannah; 2009. pp. 3128-3141

ing. 2016;54(10):5895-5913

304 Kalman Filters - Theory for Advanced Applications
