Distributed Kalman Filter

Abstract

The continuing trend toward connected sensors (“internet of things” and” ubiquitous computing”) drives a demand for powerful distributed estimation methodologies. In tracking applications, the distributed Kalman filter (DKF) provides an optimal solution under Kalman filter conditions. The optimal solution in terms of the estimation accuracy is also achieved by a centralized fusion algorithm, which receives all associated measurements. However, the centralized approach requires full communication of all measurements at each time step, whereas the DKF works at arbitrary communication rates since the calculation is fully distributed. A more recent methodology is based on ”accumulated state density” (ASD), which augments the states from multiple time instants to overcome spatial cross-correlations. This chapter explains the challenges in distributed tracking. Then, possible solutions are derived, which include the DKF and ASD approach.

1. Introduction

Modern tracking and surveillance system development is increasingly driving technological trends and algorithm developments toward networks of dislocated sensors. Besides classical target tracking, many other applications can be found, for instance, in robotics, manufacturing, health care, and financial economics. A multisensor network can exploit spatial diversity to compensate for the weak attributes of a single sensor such as limited field of view or high measurement noise. Also, heterogeneous sensors can reveal a more complete situational awareness and a more precise estimate of the underlying phenomena. Additionally, a sensor network is more robust against a single point of failure in comparison to a standalone system, if its architecture is chosen carefully.

Despite its unquestioned advantages, a multisensor network must cope with design-specific challenges, for instance, a full transmission of all the measurements to a fusion center can be hindered, when the communication links are unreliable or constrained in bandwidth or costs. A well-known approach to cope with the limited bandwidth data links is to apply data processing on the sensor sites to generate local track parameters that are transmitted to the fusion center. The latter then reconstructs the global track parameters by an application of a Track-to-Track Fusion (T2TF) scheme. Depending on the scenario constraints, this is a nontrivial task, since the local tracks are mutually correlated due to the common process noise. The first known solution in the literature to the T2TF problem proposed to apply an information filter-based multisensor fusion algorithm in [1], which later became famous as the “tracklet fusion.” However, the tracklet fusion also requires a transmission from each sensor after each time step in order to reconstruct the optimal solution.

This chapter presents the theory and the derivation of the distributed Kalman filter (DKF), which is an optimal solution of the T2TF problem under Kalman filter assumptions with respect to the mean squared error (MSE). Assuming Kalman conditions means that linear Gaussian models are provided for the motion model and all measurement models of the sensors. Moreover, it is assumed that measurement-to-track (at the sensors) and track-to-track (at the fusion center) association has been solved. The DKF approach requires, however, all remote sensor models to be known at each local sensor site, which is infeasible in most practical scenarios. Therefore, this chapter also presents a solution based on the accumulated state density (ASD), which is closely related to the DKF but does not require the measurement models to be known. Surveys that reflect the history of research in distributed estimation can be found, for instance, in [2, 3].

This chapter is structured as follows: Section 2 summarizes the problem formulation. A basic approach to T2TF is given in Section 3, where we present the least squares solution. Section 4 presents a simple fusion methodology, which is easy to compute and provides practical results for various applications. The reason why this approach is suboptimal is investigated in Section 5 by means of a recursive computation of the cross-covariances of the local tracks. In Section 6, the product representation of Gaussian probability densities is introduced, which is the basis for the derivation of the distributed Kalman filter in Section 7. An alternative derivation in terms of information parameters is provided in Section 8. Since the local measurement models are often unknown in practical applications, the distributed accumulated state density filter is introduced in Section 9. The chapter concludes with Section 10.

2. Problem formulation

Throughout this chapter, it is assumed that all S sensors produce their measurements $z_k^s\in {\mathbb{R}}^m$ at each time step $t_k$ of the same target with its true state $x_k\in {\mathbb{R}}^n$ in a synchronized way. The extension to the unsynchronized case is straightforward by means of standard Kalman filter predictions, and is therefore omitted for the sake of notational simplicity. The measurement equation for sensor $s\in \left\{1,\dots ,S\right\}$ is given by

$$z_k^s=H_k^s{x}_k+v_k^s$$E1

where $v_k^s\sim \mathrm{{\rm N}}\left(v_k^s;0,R_k^s\right)$ is the Gaussian distributed, zero-mean random variable, which models the noise of the sensing process. Therefore, the likelihood for a single measurement is fully described by the Gaussian

$$p\left(z_k^s|x_k\right)=\mathrm{{\rm N}}\left(z_k^s;H_k^s{x}_k,R_k^s\right)$$E2

Since the measurement processes across all sensors $Z_k=\left\{z_k^1,\dots ,z_k^S\right\}$ are mutually independent, the joint likelihood of all sensor data produced at time $t_k$ factorizes:

$$\begin{array}{cc}p\left(Z_k|x_k\right)& =\prod _{s=1}^{S}p\left(z_k^s|x_k\right)\\ & =\prod _{s=1}^S\mathrm{{\rm N}}\left(z_k^s;H_k^s{x}_k,R_k^s\right)\end{array}$$E3

The true state of the target itself is modeled as a time-variant stochastic process, where the transition from time $t_{k-1}$ to time $t_k$ is given by the following motion equation:

$$x_k=F_{k\mid k-1}{x}_{k-1}+w_k$$E4

where $w_k\sim \mathrm{{\rm N}}\left(w_k;0,Q_k\right)$ is the Gaussian distributed, zero-mean random variable to model the process noise of the system. Analogously to the likelihood, this provides the probability density function for the transition model:

$$p\left(x_k|x_{k-1}\right)=\mathrm{{\rm N}}\left(x_k;F_{k\mid k-1}{x}_{k-1},Q_k\right)$$E5

Based on the local processors, each sensor node produces a track at time $t_k$ in terms of an estimate $x_{k\mid k}^s$ and a corresponding estimation error covariance $P_{k\mid k}^s$. In a T2TF scheme, these parameters are the only information, which is transmitted to a fusion center (FC). It should be noted that the FC may also not be centralized, distinguished instance in the architecture, but each and every processing node can act as a FC. The introduction of a distinguished FC is only for clarification of different computation layers. An excellent overview of pros and cons of various data fusion layers can be found in [4].

The T2TF problem can now be stated as follows: compute a fused estimate $x_{k\mid k}$ of the state $x_k$ and a consistent error covariance $P_{k\mid k}$ by means of the local tracks and knowledge on their models:

$$x_{k\mid k}\leftarrow \left\{x_{k\mid k}^1,\dots ,x_{k\mid k}^S\right\}$$E6

3. Least squares estimate

In order to compute an estimate as a well-suited combination of the local tracks, it is useful to consider the joint likelihood function given by the following Gaussian:

$$p\left(x_{k\mid k}^1,\dots ,x_{k\mid k}^S|x_k\right)=\mathrm{N}\left(\left(\begin{array}{c}{x}_{k\mid k}^1\\ ⋮\\ x_{k\mid k}^S\end{array}\right);\left(\begin{array}{c}{x}_k\\ ⋮\\ x_k\end{array}\right),\left(\begin{array}{ccc}{P}_{k\mid k}^{1,1}& \cdots & P_{k\mid k}^{1,S}\\ ⋮& \ddots & ⋮\\ P_{k\mid k}^{S,1}& \cdots & P_{k\mid k}^{S,S}\end{array}\right)\right),$$E7

where $P_{k\mid k}^{s,s}=P_{k\mid k}^s$ are the track covariances on the block-diagonal entries and $P_{k\mid k}^{i,j}\triangleq \text{cov}\left[x_{k\mid k}^i,x_{k\mid k}^j|x_k\right]=P_{k\mid k}^{j,iT}$ are the cross-covariances on the off-diagonal entries of the joint error covariance.

Since the joint likelihood from above is proportional to an exponential function:

$$p\left(x_{k\mid k}^1,\dots ,x_{k\mid k}^S|x_k\right)\propto exp\left\{-\frac{1}{2}{\left(\left(\begin{array}{c}{x}_{k\mid k}^1\\ ⋮\\ x_{k\mid k}^S\end{array}\right)-\left(\begin{array}{c}{x}_k\\ ⋮\\ x_k\end{array}\right)\right)}^T{\left(\begin{array}{ccc}{P}_{k\mid k}^{1,1}& \cdots & P_{k\mid k}^{1,S}\\ ⋮& \ddots & ⋮\\ P_{k\mid k}^{S,1}& \cdots & P_{k\mid k}^{S,S}\end{array}\right)}^{-1}\left(\left(\begin{array}{c}{x}_{k\mid k}^1\\ ⋮\\ x_{k\mid k}^S\end{array}\right)-\left(\begin{array}{c}{x}_k\\ ⋮\\ x_k\end{array}\right)\right)\right\}$$E8

the maximum likelihood (ML) estimate can be computed in terms of the least squares:

$$x_{k\mid k}=arg{min}_{x_k}\left\{{\left(x_{k\mid k}^{1:S}-I^{1:S}{x}_k\right)}^T{\left(P_{k\mid k}^{1:S}\right)}^{-1}\left(x_{k\mid k}^{1:S}-I^{1:S}{x}_k\right)\right\},$$E9

where $x_{k\mid k}^{1:S}={\left(x_{k\mid k}^{1T},\dots ,x_{k\mid k}^{ST}\right)}^T$, $I^{1:S}={\left(I,\dots ,I\right)}^T,$ and $P_{k\mid k}^{1:S}$ for the joint error covariance have been introduced for notational simplicity. A closed form solution of the ML estimates can be obtained by setting the gradient with respect to the state to zero:

$$0={\nabla }_{x_k}{\left(x_{k\mid k}^{1:S}-I^{1:S}{x}_k\right)}^T{\left(P_{k\mid k}^{1:S}\right)}^{-1}\left(x_{k\mid k}^{1:S}-I^{1:S}{x}_k\right)=2{\left(I^{1:S}\right)}^T{\left(P_{k\mid k}^{1:S}\right)}^{-1}\left(x_{k\mid k}^{1:S}-I^{1:S}{x}_k\right),$$E10

Therefore, the ML estimate is given by:

$$x_{k\mid k}={\left({\left(I^{1:S}\right)}^T{\left(P_{k\mid k}^{1:S}\right)}^{-1}\left(I^{1:S}\right)\right)}^{-1}{\left(I^{1:S}\right)}^T{\left(P_{k\mid k}^{1:S}\right)}^{-1}{x}_{k\mid k}^{1:S}.$$E11

For information fusion applications, it is also important to have a consistent estimate of the squared error, in other words, we need to compute the corresponding error covariance:

$$\begin{array}{ll}\mathrm{cov}\left[x_{k\mid k}\text{|}{x}_k\right]\hfill & =\text{E}\left[\left(x_k-x_{k\mid k}\right){\left(x_k-x_{k\mid k}\right)}^{\text{T}}\right]\hfill \\ \hfill & \triangleq \text{E}\left[{\left(x_k-x_{k\mid k}\right)}^2\right]\hfill \\ \hfill & \text{E}\left[{\left(x_k-{\left({\left(I^{1:S}\right)}^{\text{T}}{\left(P_{k\mid k}^{1:S}\right)}^{-1}\left(I^{1:S}\right)\right)}^{-1}{\left(I^{1:S}\right)}^{\text{T}}{\left(P_{k\mid k}^{1:S}\right)}^{-1}{x}_{k\mid k}^{1:S}\right)}^2\right]\hfill \\ \hfill & \text{E}[({\left({\left(I^{1:S}\right)}^{\text{T}}{\left(P_{k\mid k}^{1:S}\right)}^{-1}\left(I^{1:S}\right)\right)}^{-1}(\left({\left(I^{1:S}\right)}^{\text{T}}{\left(P_{\left(k\text{|}k\right)}^{1:S}\right)}^{-1}\left(I^{1:S}\right)\right)x_k\hfill \\ \hfill & {-{\left(I^{1:S}\right)}^{\text{T}}{\left(P_{k\mid k}^{1:S}\right)}^{-1}{x}_{k\mid k}^{1:S}))}^2]\hfill \\ \hfill & ={\left({\left(I^{1:S}\right)}^{\text{T}}{\left(P_{k\mid k}^{1:S}\right)}^{-1}\left(I^{1:S}\right)\right)}^{-1}{\left(I^{1:S}\right)}^{\text{T}}{\left(P_{k\mid k}^{1:S}\right)}^{-1}\text{E}\left[{\left(\left(I^{1:S}\right)x_k-x_{k\mid k}^{1:S}\right)}^2\right]\hfill \\ \hfill & \begin{array}{c}{\left(P_{k\mid k}^{1:S}\right)}^{-1}\left(I^{1:S}\right){\left({\left(I^{1:S}\right)}^{\text{T}}{\left(P_{k\mid k}^{1:S}\right)}^{-1}\left(I^{1:S}\right)\right)}^{-1}\end{array}\hfill \\ \hfill & ={\left({\left(I^{1:S}\right)}^{\text{T}}{\left(P_{k\mid k}^{1:S}\right)}^{-1}\left(I^{1:S}\right)\right)}^{-1}.\hfill \end{array}$$E12

The last equation holds due to the fact that $\mathrm{E}\left[{\left(\left(I^{1:S}\right)x_k-x_{k\mid k}^{1:S}\right)}^2\right]=P_{k\mid k}^{1:S}$ is the joint covariance. Concluding the derivations from above equation, one can obtain:

$$\begin{array}{cc}{x}_{k\mid k}& =P_{k\mid k}{\left(I^{1:S}\right)}^{\mathrm{T}}{\left(P_{k\mid k}^{1:S}\right)}^{-1}{x}_{k\mid k}^{1:S},\\ P_{k\mid k}& ={\left({\left(I^{1:S}\right)}^{\mathrm{T}}{\left(P_{k\mid k}^{1:S}\right)}^{-1}\left(I^{1:S}\right)\right)}^{-1}\end{array}$$E13

4. Naïve fusion

It is obvious that for the ML estimate, it is assumed that the cross-covariances $P_{k\mid k}^{i,j},i,j\in \left\{1,\dots ,S\right\}$ are known. Since this might not be given in practical scenarios, a simple approximation is to assume them to be zero. This approach is called Naïve fusion. It implies that the joint error covariance is given in block-diagonal form:

$$P_{k\mid k}^{1:S}=\left(\begin{array}{ccc}{P}_{k\mid k}^{1,1}& & \\ & \ddots & \\ & & P_{k\mid k}^{S,S}\end{array}\right)$$E14

As a direct consequence of the matrix inversion lemma (see Appendix 12.1), the inverse can be obtained in closed form solution:

$${\left(P_{k\mid k}^{1:S}\right)}^{-1}=\left(\begin{array}{ccc}{\left(P_{k\mid k}^{1,1}\right)}^{-1}& & \\ & \ddots & \\ & & {\left(P_{k\mid k}^{S,S}\right)}^{-1}\end{array}\right)$$E15

Filling into the maximum likelihood formulas directly yields.

$$\begin{array}{cc}{x}_{k\mid k}& =P_{k\mid k}\sum _{s=1}^S{\left(P_{k\mid k}^s\right)}^{-1}{x}_{k\mid k}^s,\\ P_{k\mid k}& ={\left(\sum _{s=1}^S{\left(P_{k\mid k}^s\right)}^{-1}\right)}^{-1}\end{array}$$E16

Thus, by means of a simple approximation of the ML estimate, we have obtained a first practical fusion rule for the FC, which we call convex combination due to its structure for further usage. The fusion scheme as a whole can be outlined schematically as in the flowing Figure 1. Each sensor node $s$ processes its produced sensor data by means of a local filter, which results in a track in terms of an estimate together with an error covariance. These parameters are transmitted to the fusion center, which applies the convex combination to compute the fused result.

Figure 1.

Fusion scheme of the Naive Fusion approach.

5. What makes the Naïve fusion naïve?

For the Naïve fusion, we have assumed that the cross-covariances vanish. It is worth to be aware of the structure of the cross-covariances to see the conditions whether this holds or does not hold. This can be achieved by a recursive computation of the posterior cross-covariance $P_{k\mid k}^{i,j}$ of two sensors with indices $i$ and $j$, which process their data by means of local Kalman filters. At the beginning of the estimation process at $t_0$ tracks are not yet correlated, that is, $P_{0\mid 0}^{i,j}=0$, due to the fact that initial measurements are mutually uncorrelated. A recursive computation can be achieved by a prediction-filtering cycle.

5.1. Cross-covariance prediction

For the prediction step, it is assumed that a previous posterior cross-covariance $P_{k\mid k-1}^{i,j}$ has been computed. The prior parameters are obtained by means of the motion model:

$$\begin{array}{ll}{P}_{k\mid k-1}^{i,j}\hfill & =\text{E}\left[\left(x_k-x_{k\mid k-1}^i\right){\left(x_k-x_{k\mid k-1}^j\right)}^{\text{T}}\right]\hfill \\ \hfill & =\text{E}\left[\left(F_{k\mid k-1}{x}_{k-1}+w_k-F_{k\mid k-1}{x}_{k-1\mid k-1}^i\right){\left(F_{k\mid k-1}{x}_{k-1}+w_k-F_{k\mid k-1}{x}_{k-1\mid k-1}^j\right)}^{\text{T}}\right]\hfill \\ \hfill & =F_{k\mid k-1}\text{E}\left[\left(x_{k-1}-x_{k-1\mid k-1}^i\right){\left(x_{k-1}-x_{k-1\mid k-1}^j\right)}^{\text{T}}\right]F_{k\mid k-1}^{\text{T}}\hfill \\ \hfill & \begin{array}{c}-F_{k\mid k-1}\text{E}\left[\left(x_{k-1}-x_{k-1\mid k-1}^i\right)w_k{}^{\text{T}}\right]-\text{E}\left[w_k{\left(x_{k-1}-x_{k-1\mid k-1}^j\right)}^{\text{T}}\right]F_{k\mid k-1}^{\text{T}}+\text{E}\left[w_k{w}_k{}^{\text{T}}\right]\end{array}\hfill \\ \hfill & \begin{array}{c}=F_{k\mid k-1}{P}_{k-1\mid k-1}^{i,j}{F}_{k\mid k-1}^{\text{T}}+Q_k,\end{array}\hfill \end{array}$$E17

where the last equality holds due to the fact that the estimation errors at time $t_{k-1}$ of both sensor processors are uncorrelated to the process noise $w_k$.

5.2. Cross-covariance filtering

In the filtering step, both sensors compute their posterior parameters based on the produced measurements $z_k^i$ and $z_k^j,$ respectively. It is assumed that the local processors have applied local Kalman filter update steps with Kalman gains $W_{k\mid k-1}^i$ and $W_{k\mid k-1}^j$. The cross-covariance posterior matrix can be obtained by a straightforward computation:

$$\begin{array}{ll}{P}_{k\mid k}^{i,j}\hfill & =\text{E}\left[\left(x_k-x_{k\mid k}^i\right){\left(x_k-x_{k\mid k}^j\right)}^{\text{T}}\right]\hfill \\ \hfill & \begin{array}{c}=\text{E}\left[\left(x_k-x_{k\mid k-1}^i-W_{k\mid k-1}^i\left(z_k^i-H_k^i{x}_{k\mid k-1}^i\right)\right){\left(x_k-x_{k\mid k-1}^j-W_{k\mid k-1}^j\left(z_k^j-H_k^j{x}_{k\mid k-1}^j\right)\right)}^{\text{T}}\right]\end{array}\hfill \\ \hfill & \begin{array}{c}=\text{E}[\left(x_k-x_{k\mid k-1}^i-W_{k\mid k-1}^i\left(H_k^i{x}_k+v_k^i-H_k^i{x}_{k\mid k-1}^i\right)\right)(x_k-x_{k\mid k-1}^j-W_{k\mid k-1}^j\end{array}\hfill \\ \hfill & {\left(H_k^j{x}_k+v_k^j-H_k^j{x}_{k\mid k-1}^j\right))}^{\text{T}}]\hfill \\ \hfill & \begin{array}{c}=\text{E}\left[\left(\left(I-W_{k\mid k-1}^i{H}_k^i\right)\left(x_k-x_{k\mid k-1}^i\right)-W_{k\mid k-1}^i{v}_k^i\right){\left(\left(I-W_{k\mid k-1}^j{H}_k^j\right)\left(x_k-x_{k\mid k-1}^j\right)-W_{k\mid k-1}^j{v}_k^j\right)}^{\text{T}}\right]\end{array}\hfill \\ \hfill & \begin{array}{c}=\left(I-W_{k\mid k-1}^i{H}_k^i\right)\text{E}\left[\left(x_k-x_{k\mid k-1}^i\right){\left(x_k-x_{k\mid k-1}^j\right)}^{\text{T}}\right]{\left(I-W_{k\mid k-1}^j{H}_k^i\right)}^{\text{T}}+W_{k\mid k-1}^i\text{E}\left[v_k^i{v}_k^j{}^{\text{T}}\right]W_{k\mid k-1}^{j\text{T}}\end{array}\hfill \\ \hfill & \begin{array}{c}=\left(I-W_{k\mid k-1}^i{H}_k^i\right)P_{k\mid k-1}^{i,j}{\left(I-W_{k\mid k-1}^i{H}_k^i\right)}^{\text{T}}\end{array}\hfill \end{array}$$E18

For these equations, we have used the fact that the prior estimate error $\left(x_k-x_{k\mid k-1}^i\right)$ is independent of the measurement noise $v_k^i$, and that $v_k^i$ and $v_k^j$ are mutually independent.

Concluding the calculations from this section, we have obtained a recursive solution for the cross-covariances:

$$\begin{array}{cc}{P}_{0\mid 0}^{i,j}& =0\\ P_{k\mid k-1}^{i,j}& =F_{k\mid k-1}{P}_{k-1\mid k-1}^{i,j}{F}_{k\mid k-1}^{\mathrm{T}}+Q_k\\ P_{k\mid k}^{i,j}& =\left(I-W_{k\mid k-1}^i{H}_k^i\right)P_{k\mid k-1}^{i,j}{\left(I-W_{k\mid k-1}^i{H}_k^i\right)}^{\mathrm{T}}\end{array}$$E19

One can see that the cross-covariances are zero, if and only if the process noise covariance $Q_k$ vanishes. In other words, if the tracks refer to a deterministically moving target and all sensors do local Kalman filtering, then the convex combination equations yield the optimal fusion result. If this is not the case, as maybe in most practical applications, then the Naïve fusion method is an approximation and its degree of approximation depends primarily on the level of the process noise.

6. Gaussian product representation

The basic concept of the distributed Kalman filter is to make the local parameters stochastically independent, even if process noise is present. This is achieved by a product representation, which directly follows from the fact that.

$$\begin{array}{c}exp\left\{-\frac{1}{2}{\left(\left(\begin{array}{c}{x}_{k\mid k}^1\\ ⋮\\ x_{k\mid k}^S\end{array}\right)-\left(\begin{array}{c}{x}_k\\ ⋮\\ x_k\end{array}\right)\right)}^{\mathrm{T}}\left(\begin{array}{ccc}{\left(P_{k\mid k}^1\right)}^{-1}& & \\ & \ddots & \\ & & {\left(P_{k\mid k}^S\right)}^{-1}\end{array}\right)\left(\left(\begin{array}{c}{x}_{k\mid k}^1\\ ⋮\\ x_{k\mid k}^S\end{array}\right)-\left(\begin{array}{c}{x}_k\\ ⋮\\ x_k\end{array}\right)\right)\right\}\\ =exp\left\{-\frac{1}{2}\sum _{s=1}^S{\left(x_{k\mid k}^s-x_k\right)}^{\mathrm{T}}{\left(P_{k\mid k}^s\right)}^{-1}\left(x_{k\mid k}^s-x_k\right)\right\}\propto \prod _{s=1}^{S}N\left(x_k;x_{k\mid k}^s,P_{k\mid k}^s\right)\end{array}$$E20

Thus, the Gaussian product representation is equivalent to uncorrelated track parameters for each processing node. It should be noted that the product representation is not normalized, that is, the integral for $S>1$ is not unity. This, however, is not of practical relevance since the fused estimate density is a Gaussian and the parameters of which are provided by the convex combination.

7. Derivation of the distributed Kalman filter

For the DKF, we are going to modify the local processing scheme for each sensor in order to have the product representation hold at each instant of time. Then, when the fusion center receives the parameters from all sensors, the convex combination can be applied to compute the optimal global estimate. Note that the convex combination does not consider a local prior of the fusion center; therefore, the result will be independent from previous transmissions. This can be of great benefit, if communication channels with unreliable links have to be considered, since the full information on the target state is distributed in the sensor network. However, for completeness, it should also be noted that the modified local parameters are not optimal anymore in a local sense. One could say that local optimality is given up for the sake of global optimality [5].

In the following sections, the derivation of a prediction-filtering recursion of the DKF is discussed.

7.1. DKF prediction

For the prediction, it is assumed that the previous posterior is given in product representation:

$$p\left(x_{k-1}|Z^{k-1}\right)\propto \prod _{s=1}^{S}N\left(x_{k-1};x_{k-1\mid k-1}^s,P_{k-1\mid k-1}^s\right)$$E21

For notational simplicity, we have conditioned the posterior on the full data set $Z^{k-1}=\left\{Z_1^{k-1},\dots ,Z_s^{k-1}\right\}$, where $Z_s^{k-1}\left\{Z_1^s,\dots ,Z_{k-1}^s\right\}$ are from all sensors and all time steps up to time $t_{k-1}$, of which the local tracks are sufficient statistics. At the initialization phase, that is, $k-1=0$, the product representation is trivial, since the initial estimates can be based on first measurements, which are mutually independent.

To derive a closed form solution for the prediction of product representation, it is required to globalize the covariances $P_{k-1\mid k-1}^s$ of the local processing nodes at first. This process changes the local parameters such that the same previous posterior density is factorized; however, the local covariances will be unified to a single ${\stackrel{\sim }{P}}_{k-1\mid k-1}$. Rigorously speaking, this matrix does not represent a meaningful covariance in the sense of an expected estimation error squared anymore. Still, the fused result will be optimal since the global density is not changed during this process. Thus, if we set

$$p\left(x_{k-1}|Z^{k-1}\right)\propto \prod _{s=1}^{S}N\left(x_{k-1};{\stackrel{\sim }{x}}_{k-1\mid k-1}^s,{\stackrel{\sim }{P}}_{k-1\mid k-1}\right)$$E22

where

$$\begin{array}{cc}{{\displaystyle \stackrel{\sim }{x}}}_{k-1\mid k-1}^s& =S{P}_{k-1\mid k-1}{\left(P_{k-1\mid k-1}^s\right)}^{-1}{x}_{k-1\mid k-1}^s\\ {{\displaystyle \stackrel{\sim }{P}}}_{k-1\mid k-1}& =S{P}_{k-1\mid k-1}\\ P_{k-1\mid k-1}& ={\left(\sum _{s=1}^S{\left(P_{k-1\mid k-1}^s\right)}^{-1}\right)}^{-1}\end{array}$$E23

then the same fused density will be obtained, which is easily verified by means of the convex combination. It should be noted that the remote error covariances $P_{k-1\mid k-1}^s$ are required to compute the globalized covariance matrix ${\stackrel{\sim }{P}}_{k-1\mid k-1}$. Since Kalman filter conditions are assumed, they can be computed without data transmission, as they do not depend on the local sensor measurements. Therefore, it is sufficient to be aware of the remote sensor models.

The prediction formulas can now be obtained by a marginalization of the joint density of the current and the last time step:

$$\begin{array}{cc}p\left(x_k|Z^{k-1}\right)& =\int \mathrm{d}{x}_{k-1}p\left(x_k,x_{k-1}|Z^{k-1}\right)\\ & =\int \mathrm{d}{x}_{k-1}p\left(x_k|x_{k-1},Z^{k-1}\right)p\left(x_{k-1}|Z^{k-1}\right)\\ & =\int \mathrm{d}{x}_{k-1}p\left(x_k|x_{k-1}\right)p\left(x_{k-1}|Z^{k-1}\right)\end{array}$$E24

The last equality holds due to the Markov property of the system. Filling in our linear Gaussian transition model and the previous posterior yields

$$p\left(x_k|Z^{k-1}\right)\propto \int d{x}_{k-1}\mathrm{{\rm N}}\left(x_k;F_{k\mid k-1}{x}_{k-1},Q_k\right)\prod _{s=1}^{S}N\left(x_{k-1};{\stackrel{\sim }{x}}_{k-1\mid k-1}^s,{\stackrel{\sim }{P}}_{k-1\mid k-1}\right)$$E25

By means of a simple algebraic manipulation, it is possible to factorize the transition kernel Gaussian up to proportionality:

$$\begin{array}{c}\mathrm{{\rm N}}\left(x_k;F_{k\mid k-1}{x}_{k-1},Q_k\right)\propto exp\left\{-\frac{1}{2}{\left(F_{k\mid k-1}{x}_{k-1}-x_k\right)}^{\mathrm{T}}{\left(Q_k\right)}^{-1}\left(F_{k\mid k-1}{x}_{k-1}-x_k\right)\right\}\\ =exp\left\{-\frac{1}{2}{\left(F_{k\mid k-1}{x}_{k-1}-x_k\right)}^{\mathrm{T}}S{\left({\mathit{SQ}}_k\right)}^{-1}\left(F_{k\mid k-1}{x}_{k-1}-x_k\right)\right\}\\ \propto {\left[\mathrm{{\rm N}}\left(x_k;F_{k\mid k-1}{x}_{k-1},{\mathit{SQ}}_k\right)\right]}^S\end{array}$$E26

Thus, we can factorize the integration term of the global prior completely:

$$p\left(x_k|Z^{k-1}\right)\propto \int d{x}_{k-1}\prod _{s=1}^S{\rm N}\left(x_k;F_{k\mid k-1}{x}_{k-1},{\mathit{SQ}}_k\right)N\left(x_{k-1};{\stackrel{\sim }{x}}_{k-1\mid k-1}^s,{\stackrel{\sim }{P}}_{k-1\mid k-1}\right)$$E27

An application of the product formula (Section 12.2 in the appendix) yields:

$$p\left(x_k|Z^{k-1}\right)\propto \prod _{s=1}^S\mathrm{{\rm N}}\left(x_k;x_{k\mid k-1}^s,{\stackrel{\sim }{P}}_{k\mid k-1}^s\right)\int d{x}_{k-1}\prod _{s=1}^{S}N\left(x_{k-1};y^s,Y\right),$$E28

where

$$\begin{array}{cc}{x}_{k\mid k-1}^s& =F_{k\mid k-1}{{\displaystyle \stackrel{\sim }{x}}}_{k-1\mid k-1}^s\\ {{\displaystyle \stackrel{\sim }{P}}}_{k\mid k-1}^s& =F_{k\mid k-1}{{\displaystyle \stackrel{\sim }{P}}}_{k-1\mid k-1}{F_{k\mid k-1}}^{\mathrm{T}}+{\mathit{SQ}}_k\\ y^s& =Y\left({\left({{\displaystyle \stackrel{\sim }{P}}}_{k\mid k-1}^s\right)}^{-1}{x}_{k\mid k-1}^s+{F_{k\mid k-1}}^{\mathrm{T}}{\left(S{Q}_k\right)}^{-1}{x}_k\right)\\ Y& ={\left({\left({{\displaystyle \stackrel{\sim }{P}}}_{k\mid k-1}^s\right)}^{-1}+{F_{k\mid k-1}}^{\mathrm{T}}{\left(S{Q}_k\right)}^{-1}{F}_{k\mid k-1}\right)}^{-1}.\end{array}$$E29

At this point, we have derived factorized prediction formulas for the DKF prediction, and to our knowledge, the remaining integral is part of the normalization constant. This, however, is not trivial, since the parameter $y^s$ depends on $x_k$. A successive application of the product formula yields the desired result:

$$\begin{array}{c}\int \mathrm{d}{x}_{k-1}\prod _{s=1}^{S}N\left(x_{k-1}-Y{F_{k\mid k-1}}^{\mathrm{T}}{\left(S{Q}_k\right)}^{-1}{x}_k;Y{\left({\stackrel{\sim }{P}}_{k\mid k-1}^s\right)}^{-1}{x}_{k\mid k-1}^s,Y\right)\\ \propto \int \mathrm{d}{x}_{k-1}N\left(x_{k-1}-Y{F_{k\mid k-1}}^{\mathrm{T}}{\left(S{Q}_k\right)}^{-1}{x}_k;{\displaystyle \overline{\overline{y}}},{\displaystyle \overline{\overline{Y}}}\right)=1,\end{array}$$E30

for some auxiliary variables ${\displaystyle \overline{\overline{y}}}$ and ${\displaystyle \overline{\overline{Y}}}$. All factors, which are independent of $x_k$ have been omitted. This proves that the integral is independent of the state variable $x_k$. Concluding the derivations from above, we have derived the prediction formulas of the local estimation parameters as:

$$\begin{array}{cc}{{\displaystyle \stackrel{\sim }{x}}}_{k-1\mid k-1}^s& =S{P}_{k-1\mid k-1}{\left(P_{k-1\mid k-1}^s\right)}^{-1}{x}_{k-1\mid k-1}^s\\ {{\displaystyle \stackrel{\sim }{P}}}_{k-1\mid k-1}& =S{P}_{k-1\mid k-1}\\ x_{k\mid k-1}^s& =F_{k\mid k-1}{{\displaystyle \stackrel{\sim }{x}}}_{k-1\mid k-1}^s\\ {{\displaystyle \stackrel{\sim }{P}}}_{k\mid k-1}^s& =F_{k\mid k-1}{{\displaystyle \stackrel{\sim }{P}}}_{k-1\mid k-1}{F_{k\mid k-1}}^{\mathrm{T}}+{\mathit{SQ}}_k.\end{array}$$

7.2. DKF filtering

Let $Z_k$ denotes the set of measurements produced by all sensors at time $t_k$. The posterior density can be inferred by using the Bayes theorem:

$$p\left(x_k|Z^k\right)=\frac{p\left(Z_k|x_k\right)p\left(x_k|Z^{k-1}\right)}{p\left(Z_k|Z^{k-1}\right)}$$E31

Due to the mutual independence of the measurement noises, the joint likelihood function is given by:

$$p\left(Z_k|x_k\right)=\prod _{s=1}^{S}p\left(Z_k^s|x_k\right)$$E32

This is particularly useful for the structure of the product representation used for the DKF. Filling in the linear Gaussian models and neglecting the normalization constant in the denominator directly yields:

$$p\left(x_k|Z^k\right)\propto \prod _{s=1}^S\mathrm{{\rm N}}\left(Z_k;H_k^s{x}_k,R_k^s\right)\mathrm{{\rm N}}\left(x_k;x_{k\mid k-1}^s,{\stackrel{\sim }{P}}_{k\mid k-1}^s\right)$$E33

Thus, the product formula again can be applied to compute the posterior parameters:

$$p\left(x_k|Z^k\right)\propto \prod _{s=1}^S\mathrm{{\rm N}}\left(x_k;x_{k\mid k}^s,P_{k\mid k}^s\right)$$E34

where

$$\begin{array}{cc}{x}_{k\mid k}^s& =x_{k\mid k-1}^s+W_{k\mid k-1}^s{\nu }_k^s\\ P_{k\mid k}^s& =P_{k\mid k-1}^s-W_{k\mid k-1}^s{S}_k^s{W}_{k\mid k-1}^{s\mathrm{T}}\\ {\nu }_k^s& =z_k^s-H_k^s{x}_{k\mid k-1}^s\\ W_{k\mid k-1}^s& =P_{k\mid k-1}^s{H}_k^{s\mathrm{T}}{\left(S_k^s\right)}^{-1}\\ S_k^s& =H_k^s{P}_{k\mid k-1}^s{H}_k^{s\mathrm{T}}+R_k^s.\end{array}$$E35

Again, we have omitted the factors, which are independent of $x_k$.

8. Information filter formulation of the DKF

In [6], an elegant derivation of the DKF formulas was published based on the information filter (IF). The IF uses information matrices, which are inverted covariances, and information states, which are information matrices multiplied with states. The optimal, centralized update formulas for $S$ sensors based on the predicted information parameters ${\left(P_{k\mid k-1}\right)}^{-1}$ and ${\left(P_{k\mid k-1}\right)}^{-1}{x}_{k\mid k-1}$ are given by:

$$\begin{array}{cc}{\left(P_{k\mid k}\right)}^{-1}{x}_{k\mid k}& ={\left(P_{k\mid k-1}\right)}^{-1}{x}_{k\mid k-1}+\sum _{s=1}^S{i}_k^s\\ {\left(P_{k\mid k}\right)}^{-1}& ={\left(P_{k\mid k-1}\right)}^{-1}+\sum _{s=1}^S{I}_k^s\end{array}$$E36

where $i_k^s=H_k^{s\mathrm{T}}{\left(R_k^s\right)}^{-1}{z}_k^s$ and $I_k^s=H_k^{s\mathrm{T}}{\left(R_k^s\right)}^{-1}{H}_k^s$ are the local information contribution from the current measurements at time $t_k$, which were received by the FC. If we want to distribute the computation to $S$ nodes, we will have them uncorrelated as in the DKF previously. Since the fused estimate is obtained via the convex combination, the local information parameters are summed up:

$$\begin{array}{cc}{\left(P_{k\mid k}\right)}^{-1}{x}_{k\mid k}& =\sum _{s=1}^S{\left(P_{k\mid k}^s\right)}^{-1}{x}_{k\mid k}^s\\ {\left(P_{k\mid k}\right)}^{-1}& =\sum _{s=1}^S{\left(P_{k\mid k}^s\right)}^{-1}\end{array}$$E37

This summation structure can be used to provide a closed prediction-filtering cycle.

8.1. Information DKF prediction

The prediction of the state is easier than a direct transition of the information parameters. Based on the fused estimate, one can obtain.

$$\begin{array}{cc}{x}_{k\mid k-1}& =F_{k\mid k-1}{x}_{k-1\mid k-1}\\ & =F_{k\mid k-1}{P}_{k-1\mid k-1}\sum _{s=1}^S{\left(P_{k-1\mid k-1}^s\right)}^{-1}{x}_{k-1\mid k-1}^s\\ & =\sum _{s=1}^S{F}_{k\mid k-1}{P}_{k-1\mid k-1}{\left(P_{k-1\mid k-1}^s\right)}^{-1}{x}_{k-1\mid k-1}^s.\end{array}$$E38

Thus, we have given the local predicted state parameters as:

$$x_{k\mid k-1}^s=F_{k\mid k-1}{P}_{k-1\mid k-1}{\left(P_{k\mid k}^s\right)}^{-1}{x}_{k\mid k}^s$$E39

Analogously, one obtains for the prior covariance:

$$\begin{array}{cc}{P}_{k\mid k-1}& =F_{k\mid k-1}{P}_{k-1\mid k-1}{F}_{k\mid k-1}^{\mathrm{T}}+Q_k\\ & =\frac{1}{S}{F}_{k\mid k-1}{\mathit{SP}}_{k-1\mid k-1}{F}_{k\mid k-1}^{\mathrm{T}}+{\mathit{SQ}}_k\end{array}$$E40

Thus, if we set $P_{k\mid k-1}^s=S{P}_{k\mid k-1}=F_{k\mid k-1}{\mathit{SP}}_{k-1\mid k-1}{F}_{k\mid k-1}^{\mathrm{T}}+{\mathit{SQ}}_k$, then the convex combination yields the exact global fused covariance.

8.2. Information DKF filtering

For the filtering, it is assumed that each sensor has computed its local information contribution parameter $i_k^s$ and $I_k^s$ from its own sensor model and, in addition, the information matrix contributions $I_k^l$ from all remote sensors $l$ by using the individual sensor models which are again assumed to be known. Then, the information state is updated via.

$$\begin{array}{cc}{\left(P_{k\mid k}\right)}^{-1}{x}_{k\mid k}& ={\left(P_{k\mid k-1}\right)}^{-1}{x}_{k\mid k-1}+\sum _{s=1}^S{i}_k^s\\ & =\sum _{s=1}^S{\left(P_{k\mid k-1}^s\right)}^{-1}{x}_{k\mid k-1}^s+i_k^s.\end{array}$$E41

As a direct consequence, the updated parameters of the local processors follow the standard IF filtering equations:

$${\left(P_{k\mid k}^s\right)}^{-1}{x}_{k\mid k}^s={\left(P_{k\mid k-1}^s\right)}^{-1}{x}_{k\mid k-1}^s+i_k^s.$$E42

For the globalized information matrix, the remote information parameters from the sensor models are used:

$$\begin{array}{cc}{\left(P_{k\mid k}\right)}^{-1}& ={\left(P_{k\mid k-1}\right)}^{-1}+\sum _{s=1}^S{I}_k^s\\ & =\sum _{s=1}^S{\left(P_{k\mid k-1}^s\right)}^{-1}+I_k^s,\\ {\left(P_{k\mid k}^s\right)}^{-1}& ={\left(P_{k\mid k-1}^s\right)}^{-1}+I_k^s\end{array}$$E43

It is important to note that the local processing nodes compute both the local pseudo information matrix ${\left(P_{k\mid k}^s\right)}^{-1}$ and the global fused error covariance $P_{k\mid k}$, which is required for the prediction to the next time step.

9. Distributed accumulated state density filter

The DKF from Sections 7 and 8 can be considered a big step toward distributed state estimation, tracking, and information inference. However, in practical applications, the exact solution is often hindered by the fact that the exact remote sensor model parameters are unknown and can only be approximated based on local state estimates. The good news is that there is another exact solution based on the accumulated state density (ASD). The distributed ASD equations turn the spatial correlations into temporal correlations of successive states. Originally, the ASD equations were introduced to solve the out-of-sequence problem, which handles delayed transmissions of measurements into an ongoing fusion process in an optimal manner. Therefore, the temporal correlations can well be coped with the ASD approach.

At first, let us introduce the ASD state $x_{k:n}$ as

$$x_{k:n}={\left(x_k^T,\dots ,x_n^T\right)}^T$$E44

where $t_k$ refers to the current time of the filtering process and $t_n$ refers to the initialization time. The ASD approach now considers the conditional joint density $p\left(x_{k:n}|z^k\right)$, that is, the posterior of the full trajectory of the target between $t_n$ and $t_k$. In particular, the individual state densities of a single instant of time can be obtained via marginalization. Also, it should be noted that the Rauch-Tung-Striebel (RTS) smoothing equations are inherently integrated in the ASD posterior, since all states are conditioned on the full set of measurements up to time $t_k$ [7].

A recursive computation of the ASD posterior can be achieved by using the Bayes theorem:

$$p\left(x_{k:n}\right|Z^k)=\frac{p\left(Z_k|x_{k:n}\right)p\left(x_{k:n}|Z^{k-1}\right)}{p\left(Z_k|Z^{k-1}\right)}$$E45

Since the measurements conditioned on the whole trajectory only depend on the state at time $t_k$, the joint likelihood function is given by:

$$p\left(Z_k|x_{k:n}\right)=p\left(Z_k|x_k\right)=\prod _{s=1}^{S}p\left(z_k^s|x_k\right)$$E46

The second factor can be reformulated as follows:

$$\begin{array}{cc}p\left(x_{k:n}|Z^{k-1}\right)& =p\left(x_k|x_{k-1:n},Z^{k-1}\right)p\left(x_{k-1:n}|Z^{k-1}\right)\\ & =p\left(x_k|x_{k-1}\right)p\left(x_{k-1:n}|Z^{k-1}\right)\end{array}$$E47

where we have used the Markov property of the system in the last equation. This recursive representation can now be repeated on the term $p\left(x_{k-1:n}|Z^{k-1}\right)$. A successive application of this procedure yields

$$\begin{array}{c}p(x_{k:n}\left|Z^k\right)\propto \prod _{s=1}^S\left\{p\left(z_k^s|x_k\right)\right\}p\left(x_k|x_{k-1}\right)\prod _{s=1}^S\left\{p\left(z_{k-1}^s|x_{k-1}\right)\right\}p\left(x_{k-1}|x_{k-2}\right)\bullet \bullet \bullet \\ \prod _{s=1}^S\left\{p\left(z_{n+1}^s|x_{n+1}\right)\right\}p\left(x_{n+1}|x_n\right)\bullet p\left(x_n|Z^n\right)\end{array}$$E48

where we have neglected the normalization constant in the denominator. Filling in our Gaussian models and using the factorization of the transition model from above equation yields

$$p(x_{k:n}\left|Z^k\right)\propto \prod _{l=n+1}^k\prod _{s=1}^S\left\{N\left(z_l^s;H_l^s{x}_l,R_l^s\right)N\left(x_l;F_{l\mid l-1}{x}_{l-1},{\mathit{SQ}}_l\right)\right\}\bullet p\left(x_n|Z^n\right)$$E49

Since the initial density usually is based on a first measurement, we can assume that it factorizes into independent local track starts:

$$p\left(x_n|Z^n\right)\propto \prod _{s=1}^{S}N\left(x_n;x_{n\mid n}^s,P_{n\mid n}^s\right)$$E50

When the posterior is fully factorized in the number of sensors and in the time steps, each processing node can compute the resulting ASD Gaussian with mean $x_{k:n\mid k}^s$ and covariance matrix $P_{k:n\mid k}^s$ [7]:

$$p(x_{k:n}\left|z^k\right)\propto \prod _{s=1}^{S}N\left(x_{k:n};x_{k:n\mid k}^s,P_{k:n\mid k}^s\right)$$E51

where the parameters are given by:

$$\begin{array}{c}{x}_{k:n\mid k}^s={\left(x_{k\mid k}^{s\mathrm{T}},\dots ,x_{n\mid k}^{s\mathrm{T}}\right)}^{\mathrm{T}},\hfill \\ P_{k:n\mid k}^s=\left(\begin{array}{ccccc}{P}_{k\mid k}^s& P_{k\mid k}^s{W}_{k-1\mid k}^{\mathrm{s}\mathrm{T}}& P_{k\mid k}^s{W}_{k-2\mid k}^{\mathrm{s}\mathrm{T}}& \cdots & P_{k\mid k}^s{W}_{n\mid k}^{\mathrm{s}\mathrm{T}}\\ W_{k-1\mid k}^s{P}_{k\mid k}^s& P_{k-1\mid k}^s& P_{k-1\mid k}^s{W}_{k-2\mid k-1}^{\mathrm{s}\mathrm{T}}& \ast & P_{k-1\mid k}^s{W}_{n\mid k-1}^{\mathrm{s}\mathrm{T}}\\ W_{k-2\mid k}^s{P}_{k\mid k}^s& W_{k-2\mid k-1}^s{P}_{k-1\mid k}^s& P_{k-2\mid k}^s& \ast & P_{k-2\mid k}^s{W}_{n\mid k-2}^{\mathrm{s}\mathrm{T}}\\ ⋮& \ast & \ast & \ast & ⋮\\ W_{n\mid k}^s{P}_{k\mid k}^s& W_{n\mid k-1}^s{P}_{k-1\mid k}^s& \cdots & W_{n\mid +1}^s{P}_{n+1\mid k}^s& P_{n\mid k}^{s.}\end{array}\right)\end{array}$$

We have used a short notation such that $x_{l\mid k}=\mathrm{E}\left[x_l|z^k\right]$ are the smoothed state estimates and $P_{l\mid k}=cov\left[x_l|z^k\right]$ are the covariances, respectively, which result from the Rauch-Tung-Striebel equations. Also, the combined retrodiction gain matrices are known from the RTS smoother:

$$\begin{array}{cc}{W}_{l\mid k}& =\prod _{i=l}^{k-1}{W}_{i\mid i+1}\\ W_{i\mid i+1}& =P_{i\mid i}{F}_{i+1\mid i}^{\mathrm{T}}{\left(P_{i+1\mid i}\right)}^{-1}\end{array}$$E52

Thus, when the FC receives the local ASD parameters, the optimal fused estimate can be obtained via the convex combination:

$$\begin{array}{cc}{x}_{k:n\mid k}& =P_{k:n\mid k}\sum _{s=1}^S{\left(P_{k:n\mid k}^s\right)}^{-1}{x}_{k:n\mid k}^s,\\ P_{k:n\mid k}& ={\left(\sum _{s=1}^S{\left(P_{k:n\mid k}^s\right)}^{-1}\right)}^{-1}\end{array}$$E53

For a continuous state estimation process, it is convenient to formulate the distributed ASD solution in terms of a prediction-filtering cycle.

9.1. Distributed ASD prediction

For the prediction step, it is assumed that the local processing node $s$ has computed the previous filtering parameters $x_{k-1:n\mid k-1}^s$ and $P_{k-1:n\mid k-1}^s$, which refer to time $t_{k-1}$. Then, the prior parameters are given by:

$$\begin{array}{c}{x}_{k:n\mid k-1}^s={\left(x_{k\mid k-1}^{s\mathrm{T}},x_{k-1:n\mid k-1}^{s\mathrm{T}}\right)}^{\mathrm{T}}\hfill \\ P_{k:n\mid k-1}^s=\left(\begin{array}{cc}{P}_{k\mid k-1}^s& P_{k-1:n\mid k-1}^s{W}_{k-1:n}^{s\mathrm{T}}\\ W_{k-1:n}^s{P}_{k-1:n\mid k-1}^s& P_{k-1:n\mid k-1}^s\end{array}\right)\\ x_{k\mid k-1}^s=F_{k\mid k-1}{x}_{k-1\mid k-1}^s\hfill \end{array}$$E54

$$\begin{array}{c}{P}_{k\mid k-1}^s=F_{k\mid k-1}{P}_{k-1\mid k-1}^s{F}_{k\mid k-1}^{\mathrm{T}}+{\mathit{SQ}}_k\\ W_{k-1:n}^s=\left(\begin{array}{c}{W}_{k-1\mid k}^s\\ W_{k-2\mid k}^s\\ ⋮\\ W_{n\mid k}^s\end{array}\right)\hfill \end{array}$$

9.2. Distributed ASD filtering

Since the prior is factorized in form of a product representation and the current measurements from time $t_k$ are mutually uncorrelated, local Kalman filters can be applied to obtain the posterior parameters:

$$x_{k:n\mid k}^s=x_{k:n\mid k-1}^s+W_{k:n\mid k}^s\left(z_k^s-H_k^s{\Pi }_k{x}_{k:n\mid k-1}^s\right)$$

$$P_{k:n\mid k}^s=P_{k:n\mid k-1}^s-W_{k:n\mid k}^s{S}_k^s{W}_{k:n\mid k}^{s\mathrm{T}}$$

$$S_k^s=H_k^s{\Pi }_k{P}_{k:n\mid k-1}^s{\Pi }_k^{\mathrm{T}}{H}_k^{s\mathrm{T}}+R_k^s$$

$$W_{k:n\mid k}^s=P_{k:n\mid k-1}^s{\Pi }_k^T{H}_k^{sT}{\left(S_k^s\right)}^{-1}$$E55

10. Conclusion

In this chapter, we have introduced the least squares solution to the track-to-track fusion problem, where cross-covariances of the track estimation errors are required. Neglecting the cross-covariances has led us to the Naïve fusion, a simple but powerful fusion algorithm for practical applications. By recursive computation of the cross-covariances, we have seen that they primarily depend on the process noise of the state transition kernel. Since a centralized computation of the cross-covariances is infeasible in practical applications, more sophisticated solutions are required for optimal fusion results. The distributed Kalman filter, which uses the product representation to keep the local parameters decorrelated achieved this. However, this approach only works, if the local processors know all measurement models at each time step. Then, the distributed accumulated state density filter uses the temporal correlations to factorize the global posterior density. This approach does not require remote sensor models and is therefore, well suited for extensions with measurement ambiguity or nonlinear measurement functions.

In the study, one can find more extensions based on the distributed Kalman filter to overcome the lack of knowledge on the remote sensor models. In [8] and the references therein, a debiasing matrix is introduced to compensate for globally biased gain matrices of the local filters. An application of the tracklet fusion based on the distributed accumulated state density filter can be found in [9]. Then, in [10], the information filter formulation of the distributed Kalman filter also was extended to scenarios with input information on the transition process.

A. Appendix

A.1. Matrix inversion lemma

Let $A,B,C$, and $D$ be matrices of a block matrix such that $A$ and $D$ are invertable, and also such that the Schur-complements $D-C{A}^{-1}B$ and $A-B{D}^{-1}C$ have full rank. Then, the inversion of the block matrix is given by:

$${\left(\begin{array}{cc}A& B\\ C& D\end{array}\right)}^{-1}=\left(\begin{array}{cc}{\left(A-B{D}^{-1}C\right)}^{-1}& -{\left(A-B{D}^{-1}C\right)}^{-1}B{D}^{-1}\\ -D^{-1}C{\left(A-B{D}^{-1}C\right)}^{-1}& D^{-1}+D^{-1}C{\left(A-B{D}^{-1}C\right)}^{-1}B{D}^{-1}\end{array}\right)$$

$$=\left(\begin{array}{cc}{A}^{-1}+A^{-1}B{\left(D-C{A}^{-1}B\right)}^{-1}C{A}^{-1}& -A^{-1}B{\left(D-C{A}^{-1}B\right)}^{-1}\\ -{\left(D-C{A}^{-1}B\right)}^{-1}C{A}^{-1}& {\left(D-C{A}^{-1}B\right)}^{-1}\end{array}\right).$$E56

In particular, it holds that

$$\begin{array}{c}{\left(A-B{D}^{-1}C\right)}^{-1}=A^{-1}+A^{-1}B{\left(D-C{A}^{-1}B\right)}^{-1}C{A}^{-1}-{\left(A-B{D}^{-1}C\right)}^{-1}B{D}^{-1}=\hfill \\ -A^{-1}B{\left(D-C{A}^{-1}B\right)}^{-1}-D^{-1}C{\left(A-B{D}^{-1}C\right)}^{-1}=-{\left(D-C{A}^{-1}B\right)}^{-1}C{A}^{-1}{D}^{-1}\hfill \\ +D^{-1}C{\left(A-B{D}^{-1}C\right)}^{-1}B{D}^{-1}={\left(D-C{A}^{-1}B\right)}^{-1}\hfill \end{array}$$E57

Proof. Let the inverted block matrix be given by submatrices $E,F,G,$ and $H$. By definition, it holds that

$$\left(\begin{array}{cc}A& B\\ C& D\end{array}\right)\left(\begin{array}{cc}E& F\\ G& H\end{array}\right)=\left(\begin{array}{cc}I& O\\ O& I\end{array}\right)\text{and}\left(\begin{array}{cc}E& F\\ G& H\end{array}\right)\left(\begin{array}{cc}A& B\\ C& D\end{array}\right)=\left(\begin{array}{cc}I& O\\ O& I\end{array}\right),$$

where $I$ and $O$ are the identity and zero matrix, respectively. A matrix multiplication of the first and second equality yields two blocks of equations:

$$\begin{array}{cc}\mathit{AE}+\mathit{BG}& =\mathit{I}\\ \mathit{AF}+\mathit{BH}& =\mathit{O}\\ \mathit{CE}+\mathit{DG}& =\mathit{O}\\ \mathit{CF}+\mathit{DH}& =I\end{array}$$E58

which we call block A and

$$\begin{array}{cc}\mathit{EA}+\mathit{FC}& =\mathit{I}\\ \mathit{EB}+\mathit{FD}& =\mathit{O}\\ \mathit{GA}+\mathit{HC}& =\mathit{O}\\ \mathit{GB}+\mathit{HD}& =I\end{array}$$E59

which we call block B. Resolving block A for $E,F,G,$ and $H$ yields the first version of the inverted block matrix, whereas resolving block B yields the second version.

A.2. Product formula for Gaussian densities

For two Gaussian distributed random variables $x$ and $z$, it holds that

$$N\left(x;y,P\right)N\left(z;\mathit{Hx},R\right)=N\left(x;\overline{y},\overline{P}\right)N\left(z;\overline{z},S\right)$$E60

where

$$\begin{array}{cc}\overline{y}& =\left\{\begin{array}{c}y+\mathit{W\nu }\\ P^{-1}y+H^{\mathrm{T}}{R}^{-1}z\end{array}\right)\\ \overline{P}& =\left\{\begin{array}{c}P-\mathit{WS}{W}^{\mathrm{T}}\\ {\left(P^{-1}+H^{\mathrm{T}}{R}^{-1}H\right)}^{-1}\end{array}\right)\\ S& =\mathit{HP}{H}^{\mathrm{T}}+\mathit{R}\\ \mathit{W}& =P{H}^{\mathrm{T}}{S}^{-1}\\ \nu & =z-\mathit{Hy}\end{array}$$E61

A proof can be found, for instance, in [11].