Multisensor Estimation Fusion with Gaussian Process for Nonlinear Dynamic Systems

Abstract

The Gaussian process is gaining increasing importance in different areas such as signal processing, machine learning, robotics, control and aerospace and electronic systems, since it can represent unknown system functions by posterior probability. This paper investigates multisensor fusion in the setting of Gaussian process estimation for nonlinear dynamic systems. In order to overcome the difficulty caused by the unknown nonlinear system models, we associate the transition and measurement functions with the Gaussian process regression models, then the advantages of the non-parametric feature of the Gaussian process can be fully extracted for state estimation. Next, based on the Gaussian process filters, we propose two different fusion methods, centralized estimation fusion and distributed estimation fusion, to utilize the multisensor measurement information. Furthermore, the equivalence of the two proposed fusion methods is established by rigorous analysis. Finally, numerical examples for nonlinear target tracking systems demonstrate the equivalence and show that the multisensor estimation fusion performs better than the single sensor. Meanwhile, the proposed fusion methods outperform the convex combination method and the relaxed Chebyshev center covariance intersection fusion algorithm.

1. Introduction

Estimation in nonlinear systems is extremely important because almost all practical systems involve nonlinearities of one kind or another [1,2], such as target tracking, vehicle navigation, automatic control and robotics [3,4]. In the case of nonlinearities, estimation cannot be obtained analytically in general. Some methods based on exact parametric models and ideas of the Kalman filter (KF) have been developed. The Extended Kalman filter (EKF) [5] was the most common application to nonlinear systems, which simply linearizes all nonlinear functions via Taylor-series expansion and substitutes Jacobian matrices for the linear transformations into the KF equations. The unscented transformation was introduced to address the deficiencies of linearization, namely unscented Kalman filters (UKF) [1], which provided a more direct and explicit mechanism for transforming mean and covariance matrices. Under the Bayesian framework, particle filter (PF) was presented by constructing the posterior probability density function of the state based on all available information in Reference [6]. However, for nonlinear systems, these methods were usually assumed that the nonlinear relationships are known. Lack of modeling accuracy, including the identification of the noise and the model parameters, was inevitable [7]. In many applications, most real dynamic systems are difficult to obtain the exact models and system noises due to the complexity of the target motion environment, therefore, the parameterized estimation methods may be invalid.

To overcome the limitations of parametric models, researchers have recently employed so called non-parametric methods like the Gaussian process [8] to learn models for dynamic systems. More specifically, the functional representation needs to be learned from the training data before the filtering prediction and update step [9]. Gaussian processes have been increasingly attracting the interests in machine learning, signal processing, robotics, control and aerospace and electronic systems [10,11,12]. For example, Gaussian process models are used as the surrogate models for complex physics models in Reference [13]. The advantages stemmed from the fact that Gaussian processes take both the noise in the system and the uncertainty in the model into consideration [14]. In the context of modeling the dynamic systems, Gaussian processes can be used as prior over the transition function and measurement function. By analyzing the correlation among given training data, Gaussian process models can provide the posterior distributions over functions through the combination of the prior and the data. For the cases that ground truth data are unavailable or can only be determined approximately, Gaussian process latent variable models were developed in References [15,16] and they were extended to the setting of dynamical robotics systems in Reference [17]. In order to reduce the cubic complexity of Gaussian process training for the fixed number of training points, sparse Gaussian processes were developed (see e.g., [18,19,20,21,22,23]). K-optimality was used to improve the stability of the Gaussian process prediction in Reference [24].

So far, Gaussian process models have been applied successfully to massive nonlinear dynamic systems. Motivated by modeling human motion, learning nonlinear dynamical models with Gaussian process was investigated in Reference [16]. Many filtering methods, such as the Gaussian process extended Kalman filters (GP-EKF) [25], the Gaussian process unscented Kalman filters (GP-UKF) [26], the Gaussian process particle filters (GP-PF) [27], GP-BayesFilters [14] and the Gaussian process assumed density filters (GP-ADF) [7,10] were derived by incorporating the Gaussian process transition model and Gaussian process measurement model into the classic filters. GP-ADF was an efficient form of assumed density filtering (ADF) introduced in References [28,29,30] and propagated the full Gaussian density [7]. Although Gaussian processes have been around for decades, they mainly focus on single sensor systems.

With the rapid development of the sensor technology, computer science, communication technology and information technology, multisensor systems are widely used in military and civil fields [31,32,33,34,35]. Benefited from the application of multiple sensors, multisensor data fusion makes more comprehensive and accurate decision by integrating the available information from multiple sensors and has attracted lots of research interests. Generally speaking, the multisensor estimation fusion mainly contains centralized estimation fusion and distributed estimation fusion. Centralized estimation fusion is that a central processor receives all measurement data from sensors without processing and uses them to estimate the state [36,37]. In general, many filtering algorithms in single sensor system can be applied to the multisensor systems, since measurement value can be stacked and regarded as one measurement. On the other hand, distributed estimation fusion has several advantages in terms of reliability, survivability, communication bandwidth and computational burdens [38,39,40,41,42], which make it desirable in real applications such as surveillance, tracking and battle management. In the distributed setting, each local sensor deals with its own measurement data and transmits the local state estimation to the fusion center for the purpose of more accurate estimation fusion. So far, a variety of distributed fusion methods have been investigated for different occasions, such as References [37,40,43,44,45,46]. A convex combination method was given to fuse the local estimates in Reference [43]. For the unavailable cross correlation matrix problem, a relaxed Chebyshev center covariance intersection (RCC-CI) algorithm was also provided to fuse the local estimates in Reference [37]. A current review of distributed multisensor data fusion under unknown correlation can be seen in Reference [47]. However, these multisensor estimation fusion methods are mainly based on the exact dynamic models or known nonlinear functions. For the cases in which accurate parametric models are difficult to obtain, it is worth considering integrating Gaussian processes with multisensor estimation fusion to improve the system’s performance.

In this paper, we focus on multisensor estimation fusion including the centralized and distributed fusion methods with Gaussian process for nonlinear dynamic systems. Firstly, given the training data, we learn the dynamic models with the Gaussian process and derive multisensor estimation fusion methods based on the Gaussian process models for nonlinear dynamic systems, which can avoid inappropriate parametric models and improve predictive ability. Combining with the single senor GP-ADF and GP-UKF, respectively, the prediction step and update step of the multisensor estimation fusion are provided. In general, it is hard to analyze the performance of the non-parametric fusion methods. Since the Gaussian process fusion methods have analytic mean and covariance, we show that the distributed estimation fusion is equivalent to the centralized estimation fusion with the single sensor cross terms being full column rank. Numerical examples show that the equivalence is satisfied under the condition and the multisensor estimation fusion performs better than the single sensor. We also compare the proposed fusion methods with the RCC-CI [37] algorithm and the convex combination method [43]. In the Gaussian process model setting, the simulation results show that the multisensor estimation fusion methods outperform the RCC-CI algorithm and the convex combination method, as far as the estimation accuracy is concerned.

This article also extends our earlier work [48]. Compared with the conference paper, the main differences are as follows:

• Detailed proofs are provided.
• The enhancement of the multisensor fusion algorithm with GP-UKF is presented.
• Additional set of extensive experiments are carried out.
• The equivalence condition of Proposition 1 is analyzed.
• Comparison between GP-ADF fusion and GP-UKF is given.

The rest of the paper is organized as follows. In Section 2, the problem formulation and the Gaussian process are introduced. In Section 3, the centralized estimation fusion and the distributed estimation fusion methods are presented. In Section 4, simulations are provided to confirm our analysis. Some conclusions are drawn in Section 5.

2. Preliminaries

2.1. Problem Formulation

In this paper, we consider the state estimation problem of a nonlinear dynamic system with additive noise and N ($N\ge 2$) sensor measurements. The multisensor nonlinear dynamic system with state equation and N measurement equations (see Figure 1) is described as follows:

$$\begin{array}{cc}\hfill {\mathbf{x}}_k& =h\left({\mathbf{x}}_{k-1}\right)+{\mathbf{w}}_k,k=0,1,\dots ,\hfill \end{array}$$

(1)

$$\begin{array}{cc}\hfill {\mathbf{y}}_k^m& =g^m\left({\mathbf{x}}_k\right)+{\mathbf{v}}_k^m,m=1,\dots ,N,\hfill \end{array}$$

(2)

where ${\mathbf{x}}_k\in {\mathcal{R}}^D$ is the state of the dynamic system at time k and ${\mathbf{y}}_k^m\in {\mathcal{R}}^{l_m}$ is the measurement of the $m_{th}$ sensor at time k, $m=1,\dots ,N$. $h\left({\mathbf{x}}_k\right)$ is the transition function of the state ${\mathbf{x}}_k$ and $g^m\left({\mathbf{x}}_k\right)$ is the nonlinear measurement function of ${\mathbf{x}}_k$ at the $m_{th}$ sensor. ${\mathbf{w}}_k\sim N(\mathbf{0},{\mathbf{Q}}_k)$ is a Gaussian system noise and ${\mathbf{v}}_k^m\sim N(\mathbf{0},{\mathbf{R}}_k^m)$ is a Gaussian measurement noise, which is independent of each other.

Our goal is to estimate the state from all the available sensor measurement information. Gaussian processes are regarded as the prior of the transition function $h\left({\mathbf{x}}_{k-1}\right)$ and the sensor measurement function $g^m\left({\mathbf{x}}_k\right)$, $m=1,\dots ,N$, respectively. Then we make inference about the posterior distribution of the transition function and the sensor measurement function. The Gaussian process model represents a powerful tool to make Bayesian inference about functions [49].

2.2. Gaussian Processes

A Gaussian process is defined over functions, which is a generalization of the Gaussian probability distribution [8]. It means that we should consider inference directly in function space. Also, it gives a formal definition, “the Gaussian process is defined a collection of random variables, any finite number of which have a joint Gaussian distribution,” in [8].

Similar to a Gaussian distribution, knowledge of both mean and covariance function can specify a Gaussian process. The mean function $m\left(\mathbf{x}\right)$ and covariance function (also called a kernel) $k(\mathbf{x},{\mathbf{x}}^{\prime })$ of a Gaussian process $f\left(\mathbf{x}\right)$ are defined as follows,

$$\begin{array}{cc}\hfill m\left(\mathbf{x}\right)& =\mathbb{E}\left[f\right(\mathbf{x}\left)\right],\hfill \\ \hfill k(\mathbf{x},{\mathbf{x}}^{\prime })& =\mathbb{E}\left[(f\left(\mathbf{x}\right)-m\left(\mathbf{x}\right))(f\left({\mathbf{x}}^{\prime }\right)-m\left({\mathbf{x}}^{\prime }\right))\right],\hfill \end{array}$$

where $\mathbb{E}\left[·\right]$ denotes the expectation. Thus, we write the Gaussian process as

$$\begin{array}{c}\hfill f\left(\mathbf{x}\right)\sim \mathcal{GP}(m\left(\mathbf{x}\right),k(\mathbf{x},{\mathbf{x}}^{\prime })).\end{array}$$

Unless stated otherwise, a prior mean function is usually assumed to be 0. The choice of the covariance functions depends on the application [14]. In this paper, we employ the most commonly-used squared exponential (SE) kernel in machine learning, which is the prototypical stationary covariance function and useful for modeling particularly smooth functions. It is defined as

$$\begin{array}{cc}\hfill & \mathrm{Cov}(f\left(\mathbf{x}\right),f\left({\mathbf{x}}^{\prime }\right))=k(\mathbf{x},{\mathbf{x}}^{\prime })\hfill \\ \hfill & ={\alpha }^2\exp \{-\frac{1}{2}{(\mathbf{x}-{\mathbf{x}}^{\prime })}^{\mathrm{T}}{\mathbf{\Lambda }}^{-1}(\mathbf{x}-{\mathbf{x}}^{\prime })\},\hfill \end{array}$$

(3)

where parameter ${\alpha }^2$ represents the variance of the function f that controls the uncertainty of predictions in areas of less training sample density and parameter $\mathbf{\Lambda }$ is a diagonal matrix of the characteristic length-scales of the SE kernel. Other commonly employed kernel functions can be seen in Reference [8].

A Gaussian process implies a distribution over functions based upon the obtained training data. Assume that we have obtained a set of training data $\mathcal{X}=\{\mathbf{X},\mathbf{y}\}$. $\mathbf{X}$ and $\mathbf{y}$ are made up of multiple samples drawn from the following standard Gaussian process regression model

$$\begin{array}{c}\hfill y=f\left(\mathbf{x}\right)+\epsilon ,\epsilon \sim \mathcal{N}(0,{\sigma }_{\epsilon }^2),\end{array}$$

(4)

where $f:{\mathcal{R}}^D\to \mathcal{R}$ and $f\left(\mathbf{x}\right)\sim \mathcal{GP}$, $\mathcal{N}$ denotes the normalized Gaussian probability density function. Note that Gaussian process regression uses the fact that any finite set of training data and testing data of a Gaussian process are jointly Gaussian distributed. Let $\mathit{\theta }=\{{\alpha }^2,\mathbf{\Lambda },{\sigma }_{\epsilon }^2\}$, called the hyper-parameters of the Gaussian process. Using the evidence maximization method [8,50], we obtain a point estimate

$$\begin{array}{cc}\hfill \widehat{\mathit{\theta }}& =arg\underset{\mathit{\theta }}{max}\log p(\mathbf{y}\mid \mathbf{X},\mathit{\theta }),\hfill \end{array}$$

from the known training data [51]. We can solve the optimization problem with numerical optimization techniques such as conjugate gradient ascent [8,14]. Next, we infer the posterior distribution over function value $f\left({\mathbf{x}}_{\ast }\right)$ from training data for each inputs ${\mathbf{x}}_{\ast }\in {\mathcal{R}}^D$. In addition, the test input ${\mathbf{x}}_{\ast }$ can be categorized into two classes, relying on whether it is uncertain or not. The corresponding predictive distribution over $f\left({\mathbf{x}}_{\ast }\right)$ will be presented as follows.

(1) Predictive Distribution over Univariate Function

Suppose that the training data is $\mathcal{X}=(\mathbf{X},\mathbf{y})$=$\{({\mathbf{x}}_i,y_i):i=1,\dots ,n\}$, where ${\mathbf{x}}_i$ is a D-dimensional input vector, $y_i$ is a scalar output, and n is the number of training data. Next, two cases, deterministic inputs and uncertain inputs, are considered.

● Deterministic Inputs

Conditioned on the training data and the deterministic test input ${\mathbf{x}}_{\ast }$, the predictive distribution over $f\left({\mathbf{x}}_{\ast }\right)$ is a Gaussian with the mean

$$\begin{array}{c}\hfill m_f\left({\mathbf{x}}_{\ast }\right)=\mathbb{E}\left[f_{\ast }\right]={\mathbf{k}}_{\ast }^{\mathrm{T}}{(\mathbf{K}+{\sigma }_{\epsilon }^2\mathbf{I})}^{-1}y={\mathbf{k}}_{\ast }^{\mathrm{T}}\mathit{\beta },\end{array}$$

(5)

and the variance

$$\begin{array}{c}\hfill {\sigma }_f^2\left({\mathbf{x}}_{\ast }\right)={\mathrm{Var}}_f\left[f_{\ast }\right]=k_{\ast \ast }-{\mathbf{k}}_{\ast }^{\mathrm{T}}{(\mathbf{K}+{\sigma }_{\epsilon }^2\mathbf{I})}^{-1}{\mathbf{k}}_{\ast },\end{array}$$

(6)

where ${\mathrm{Var}}_f\left[·\right]$ represents the variance with respect to f,

$$\begin{array}{cc}\hfill {\mathbf{k}}_{\ast }& :=k(\mathbf{X},{\mathbf{x}}_{\ast }),\hfill \\ \hfill k_{\ast \ast }& :=k({\mathbf{x}}_{\ast },{\mathbf{x}}_{\ast }),\hfill \\ \hfill \mathit{\beta }& :={(\mathbf{K}+{\sigma }_{\epsilon }^2\mathbf{I})}^{-1}\mathbf{y},\hfill \end{array}$$

and $\mathbf{K}$ is the kernel matrix with each element satisfying $K_{ij}=k({\mathbf{x}}_i,{\mathbf{x}}_j)$. The uncertainty of prediction is characterized by the variance ${\sigma }_f^2\left({\mathbf{x}}_{\ast }\right)$. More details can be seen in Reference [8].

● Uncertain Inputs

When the test input is uncertain, namely ${\mathbf{x}}_{\ast }$ has a probability distribution, the prediction problem is relatively difficult. Let us consider the predictive distribution of $f\left({\mathbf{x}}_{\ast }\right)$ for the uncertain input ${\mathbf{x}}_{\ast }\sim \mathcal{N}(\mathit{\mu },\mathbf{\Sigma })$. According to the results in Reference [7] (or reviewing References [52,53]), we introduce the predictive distribution over the function value as

$$\begin{array}{c}\hfill p\left(f\left({\mathbf{x}}_{\ast }\right)\right|\mathit{\mu },\mathbf{\Sigma })=\int p\left(f\left({\mathbf{x}}_{\ast }\right)\right|{\mathbf{x}}_{\ast })p\left({\mathbf{x}}_{\ast }\right|\mathit{\mu },\mathbf{\Sigma })d{\mathbf{x}}_{\ast },\end{array}$$

(7)

where the mean and variance of the distribution $p\left(f\left({\mathbf{x}}_{\ast }\right)\right|{\mathbf{x}}_{\ast })$ are provided with Equations (5) and (6), respectively. Based on the conditional expectation and variance formulae, it yields the mean ${\mu }_{\ast }$ and variance ${\sigma }_{\ast }^2$ of the distribution $p\left(f\left({\mathbf{x}}_{\ast }\right)\right|\mathit{\mu },\mathbf{\Sigma })$ in closed form. In particular, we have

$$\begin{array}{cc}\hfill {\mu }_{\ast }& ={\mathbb{E}}_{{\mathbf{x}}_{\ast }}\left[{\mathbb{E}}_f\left(f\left({\mathbf{x}}_{\ast }\right)\right|{\mathbf{x}}_{\ast })\right|\mathit{\mu },\mathbf{\Sigma }]\hfill \\ \hfill & ={\mathbb{E}}_{{\mathbf{x}}_{\ast }}\left[m_f\left({\mathbf{x}}_{\ast }\right)\right|\mathit{\mu },\mathbf{\Sigma }]\hfill \\ \hfill & =\int m_f\left({\mathbf{x}}_{\ast }\right)\mathcal{N}\left({\mathbf{x}}_{\ast }\right|\mathit{\mu },\mathbf{\Sigma })d{\mathbf{x}}_{\ast }\hfill \\ \hfill & ={\mathit{\beta }}^{\mathrm{T}}\mathbf{l},\hfill \end{array}$$

(8)

and

$$\begin{array}{cc}\hfill {\sigma }_{\ast }^2& ={\mathbb{E}}_{{\mathbf{x}}_{\ast }}\left[m_f{\left({\mathbf{x}}_{\ast }\right)}^2\right|\mathit{\mu },\mathbf{\Sigma }]+{\mathbb{E}}_{{\mathbf{x}}_{\ast }}\left[{\sigma }_f^2\left({\mathbf{x}}_{\ast }\right)\right|\mathit{\mu },\mathbf{\Sigma }]\hfill \\ \hfill & -{\mathbb{E}}_{{\mathbf{x}}_{\ast }}{\left[m_f\left({\mathbf{x}}_{\ast }\right)\right|\mathit{\mu },\mathbf{\Sigma }]}^2\hfill \\ \hfill & ={\mathit{\beta }}^{\mathrm{T}}\tilde{\mathbf{L}}\mathit{\beta }+{\alpha }^2-tr\left({(\mathbf{K}+{\sigma }_{\epsilon }^2\mathbf{I})}^{-1}\tilde{\mathbf{L}}\right)-{\left({\mu }_{\ast }\right)}^2.\hfill \end{array}$$

(9)

Note that $\mathbf{l}$ in Equation (8) is from $\mathbf{l}={[l_1,\dots ,l_n]}^{\mathrm{T}}$,

$$\begin{array}{cc}\hfill l_j& =\int k_f({\mathbf{x}}_j,{\mathbf{x}}_{\ast })p\left({\mathbf{x}}_{\ast }\right)d{\mathbf{x}}_{\ast }={\alpha }^2{|\mathbf{\Sigma }{\mathbf{\Lambda }}^{-1}+\mathbf{I}|}^{-\frac{1}{2}}\hfill \\ \hfill & \times \exp \{-\frac{1}{2}{({\mathbf{x}}_j-\mathit{\mu })}^{\mathrm{T}}{(\mathbf{\Sigma }+\mathbf{\Lambda })}^{-1}({\mathbf{x}}_j-\mathit{\mu })\}.\hfill \end{array}$$

$tr\left(·\right)$ represents the trace in Equation (9) and we have

$$\begin{array}{cc}\hfill {\tilde{L}}_{ij}& =\frac{k_f({\mathbf{x}}_i,\mathit{\mu })k_f({\mathbf{x}}_j,\mathit{\mu })}{|2\mathbf{\Sigma }{\mathbf{\Lambda }}^{-1}{+\mathbf{I}|}^{\frac{1}{2}}}\hfill \\ \hfill & \times \exp \{{({\tilde{\mathbf{z}}}_{ij}-\mathit{\mu })}^{\mathrm{T}}{(\mathbf{\Sigma }+\frac{1}{2}\mathbf{\Lambda })}^{-1}\mathbf{\Sigma }{\mathbf{\Lambda }}^{-1}({\tilde{\mathbf{z}}}_{ij}-\mathit{\mu })\},\hfill \end{array}$$

and

$$\begin{array}{c}\hfill {\tilde{\mathbf{z}}}_{ij}:=\frac{1}{2}({\mathbf{x}}_i+{\mathbf{x}}_j).\end{array}$$

In general, the predictive distribution in Equation (7) cannot be analytically calculated since it leads to a non-Gaussian distribution, if a Gaussian distribution is mapped through a nonlinear function. Using moment-matching method, the distribution $p\left(f\left({\mathbf{x}}_{\ast }\right)\right|\mathit{\mu },\mathbf{\Sigma })$ can be approximated by the Gaussian distribution $\mathcal{N}({\mu }_{\ast },{\sigma }_{\ast }^2)$.

(2) Predictive Distribution over Multivariate Function

For the model (4), let us turn to the multiple output case that $f:{\mathcal{R}}^D\to {\mathcal{R}}^E$, $f\sim \mathcal{GP}$. Following the results in References [7,51], we need to train E Gaussian process regression models independently. The ath model is learned from the training data $[X,{\mathbf{y}}^a]$, ${\mathbf{y}}^a={[y_1^a,\dots ,y_n^a]}^{\mathrm{T}}$, $a=1,\dots ,E$, where $y_j^a$ is the ath element of $y_j$. It implies that any two targeted dimensions are conditionally independent given the input.

For any deterministic input ${\mathbf{x}}_{\ast }$, the mean and variance of each target dimension are obtained by Equations (5) and (6). The predictive mean of $f\left({\mathbf{x}}_{\ast }\right)$ is a vector stacked by E targeted dimension mean, and the corresponding covariance is a diagonal matrix whose diagonal element is the variance of each targeted dimension.

With the uncertain input ${\mathbf{x}}_{\ast }\sim \mathcal{N}(\mathit{\mu },\mathbf{\Sigma })$, the predictive mean ${\mathit{\mu }}_{\ast }$ of $f\left({\mathbf{x}}_{\ast }\right)$ also is the stacked E predictive mean ${\mathit{\mu }}_{\ast }^a$ given by Equation (8), $a=1,\dots ,E$. Unlike predicting at deterministic inputs, targeted dimensions are dependent due to the uncertain input. Denote $f_a^{\ast }=f_a\left({\mathbf{x}}_{\ast }\right)$ and the corresponding predictive covariance matrix is given by

$$\begin{array}{cc}\hfill {\mathbf{\Sigma }}_{\ast }& \mid \mathit{\mu },\mathbf{\Sigma }\hfill \\ \hfill & =\left[\begin{array}{ccc}\mathrm{Var}\left[f_1^{\ast }\right|\mathit{\mu },\mathbf{\Sigma }]& \cdots & \mathrm{Cov}[f_1^{\ast },f_E^{\ast }|\mathit{\mu },\mathbf{\Sigma }]\\ ⋮& \ddots & ⋮\\ \mathrm{Cov}[f_E^{\ast },f_1^{\ast }|\mathit{\mu },\mathbf{\Sigma }]& \cdots & \mathrm{Var}\left[f_E^{\ast }\right|\mathit{\mu },\mathbf{\Sigma }]\end{array}\right].\hfill \end{array}$$

(10)

It is obvious that the predictive covariance is no longer a diagonal matrix. Using Equation (9), we can obtain the diagonal element of covariance matrix (10). For each $a,b\in \{1,\dots ,E\}$, the off-diagonal elements satisfy

$$\begin{array}{c}\hfill \mathrm{Cov}[f_a^{\ast },f_b^{\ast }|\mathit{\mu },\mathbf{\Sigma }]={\mathbb{E}}_{f,{\mathbf{x}}_{\ast }}\left[f_a\left({\mathbf{x}}_{\ast }\right)f_b\left({\mathbf{x}}_{\ast }\right)\right|\mathit{\mu },\mathbf{\Sigma }]-{\mu }_{\ast }^a{\mu }_{\ast }^b.\end{array}$$

Given ${\mathbf{x}}_{\ast }$, $f_a\left({\mathbf{x}}_{\ast }\right)$ and $f_b\left({\mathbf{x}}_{\ast }\right)$ are independent, then we have

$$\begin{array}{cc}\hfill & {\mathbb{E}}_{f,{\mathbf{x}}_{\ast }}\left[f_a^{\ast }{f}_b^{\ast }\right|\mathit{\mu },\mathbf{\Sigma }]\hfill \\ \hfill =& \int \int f_a^{\ast }{f}_b^{\ast }p(f_a,f_b|{\mathbf{x}}_{\ast })p\left({\mathbf{x}}_{\ast }\right|\mathit{\mu },\mathbf{\Sigma })dfd{\mathbf{x}}_{\ast }\hfill \\ \hfill =& {\mathit{\beta }}_a^{\mathrm{T}}\int k_f^a(\mathbf{X},{\mathbf{x}}_{\ast })k_f^b({\mathbf{x}}_{\ast },\mathbf{X})p\left({\mathbf{x}}_{\ast }\right|\mathit{\mu },\mathbf{\Sigma })d{\mathbf{x}}_{\ast }{\mathit{\beta }}_b,\hfill \end{array}$$

(11)

where ${\mathit{\beta }}_a,{\mathit{\beta }}_b$ can be similarly obtained in Equation (5). For notational simplicity, define

$$\begin{array}{cc}\hfill \mathbf{L}& :=\int k_f^a(\mathbf{X},{\mathbf{x}}_{\ast })k_f^b({\mathbf{x}}_{\ast },\mathbf{X})p\left({\mathbf{x}}_{\ast }\right|\mathit{\mu },\mathbf{\Sigma })d{\mathbf{x}}_{\ast },\hfill \end{array}$$

where

$$\begin{array}{cc}\hfill L_{ij}& ={\alpha }_a^2{\alpha }_b^2{|({\mathbf{\Lambda }}_a^{-1}+{\mathbf{\Lambda }}_b^{-1})\mathbf{\Sigma }+\mathbf{I}|}^{-\frac{1}{2}}\hfill \\ \hfill & \times \exp \{-\frac{1}{2}{({\mathbf{x}}_i-{\mathbf{x}}_j)}^{\mathrm{T}}{({\mathbf{\Lambda }}_a+{\mathbf{\Lambda }}_b)}^{-1}({\mathbf{x}}_i-{\mathbf{x}}_j)\}\hfill \\ \hfill & \times \exp \{-\frac{1}{2}{({\mathbf{z}}_{ij}-\mathit{\mu })}^{\mathrm{T}}{\mathbf{R}}^{-1}({\mathbf{z}}_{ij}-\mathit{\mu })\},\hfill \\ \hfill \mathbf{R}& :={({\mathbf{\Lambda }}_a^{-1}+{\mathbf{\Lambda }}_b^{-1})}^{-1}+\mathbf{\Sigma },\hfill \\ \hfill {\mathbf{z}}_{ij}& :={\mathbf{\Lambda }}_b{({\mathbf{\Lambda }}_a+{\mathbf{\Lambda }}_b)}^{-1}{\mathbf{x}}_i+{\mathbf{\Lambda }}_a{({\mathbf{\Lambda }}_a+{\mathbf{\Lambda }}_b)}^{-1}{\mathbf{x}}_j.\hfill \end{array}$$

Thus, we also approximate the distribution $p\left(f\left({\mathbf{x}}_{\ast }\right)\right|\mathit{\mu },\mathbf{\Sigma })$ with the Gaussian distribution $\mathcal{N}({\mathit{\mu }}_{\ast },{\mathbf{\Sigma }}_{\ast })$. More details can be found in Reference [51].

3. Multisensor Estimation Fusion

In this section, an essential prerequisite is that the transition function (1) and the measurement functions (2) are either not known or no longer accessible. Thus, we model the latent functions with Gaussian processes. For the N-sensor dynamic systems (1) and (2), suppose that we have obtained some training data of the target state and sensor measurement. In the following, we discuss the centralized and distributed estimation fusion methods to estimate the state from a sequence of noisy sensor measurements.

3.1. Centralized Estimation Fusion

First, we consider the centralized estimation fusion method. To facilitate fusion, we stack the multisensor measurement information as follows

$$\begin{array}{cc}\hfill & {\mathbf{y}}_k={[{\left({\mathbf{y}}_k^1\right)}^{\mathrm{T}},\dots ,{\left({\mathbf{y}}_k^N\right)}^{\mathrm{T}}]}^{\mathrm{T}},\hfill \\ \hfill & g\left({\mathbf{x}}_k\right)={[{\left(g^1\left({\mathbf{x}}_k\right)\right)}^{\mathrm{T}},\dots ,{\left(g^N\left({\mathbf{x}}_k\right)\right)}^{\mathrm{T}}]}^{\mathrm{T}},\hfill \\ \hfill & {\mathbf{v}}_k={[{\left({\mathbf{v}}_k^1\right)}^{\mathrm{T}},\dots ,{\left({\mathbf{v}}_k^N\right)}^{\mathrm{T}}]}^{\mathrm{T}}.\hfill \end{array}$$

(12)

Thus the dynamic system (1) and (2) can be rewritten as

$$\begin{array}{cc}\hfill {\mathbf{x}}_k& =h\left({\mathbf{x}}_{k-1}\right)+{\mathbf{w}}_k,\hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill {\mathbf{y}}_k& =g\left({\mathbf{x}}_k\right)+{\mathbf{v}}_k,\hfill \end{array}$$

(14)

where the covariance of the noise ${\mathbf{v}}_k$ is given by

$$\begin{array}{cc}\hfill \mathrm{Cov}\left({\mathbf{v}}_k\right)=& {\mathbf{\Sigma }}_k=diag({\mathbf{\Sigma }}_k^1,\dots ,{\mathbf{\Sigma }}_k^N),\hfill \\ \hfill \mathrm{Cov}\left({\mathbf{v}}_k^m\right)=& {\mathbf{\Sigma }}_k^m,m=1,\dots ,N.\hfill \end{array}$$

Considering the Gaussian process dynamic system setup, two Gaussian process models can be trained using evidence maximization. ${\mathcal{GP}}_h$ models the mapping ${\mathbf{x}}_{k-1}\mapsto {\mathbf{x}}_k,{\mathcal{R}}^D\to {\mathcal{R}}^D$ and ${\mathcal{GP}}_g$ models the mapping ${\mathbf{x}}_k\mapsto {\mathbf{y}}_k,{\mathcal{R}}^D\to {\mathcal{R}}^{N\ast E}$.

As we know, the key of the estimation is to recursively infer the posterior distribution over the state ${\mathbf{x}}_k$ of the dynamic system (13)–(14), which are based on all the sensor measurements ${\mathbf{y}}_{1:k}={\left\{{\mathbf{y}}_{\tau }\right\}}_{\tau =1}^k$, including the prediction step and the update step. In particular, the prediction step uses the posterior distribution from the previous time step to produce a predictive distribution of the state at the current time step. In the update step, the current predictive distribution is combined with current measurement information to refine the posterior distribution at the current time step. Next, centralized estimation fusion methods with GP-ADF and GP-UKF are presented, respectively.

(1) Fusion with GP-ADF

Note that some approximation methods are required to obtain an analytic posterior distribution for the nonlinear dynamic systems. In particular, for the Gaussian process dynamic system, it is easy to make some Gaussian approximation to the posterior distribution. GP-ADF exploits the fact that the true moments of the Gaussian process predictive distribution can be computed in closed form. The predictive distribution is approximated by a Gaussian with the exact predictive mean and covariance [7]. Based on the following lemma, we present the centralized estimation fusion method with GP-ADF.

Lemma 1.

For the Gaussian process dynamic system (13)–(14), the posterior distribution of the state ${\mathbf{x}}_k$ can be approximated by the Gaussian distribution $\mathcal{N}({\mathit{\mu }}_k^e,{\mathbf{C}}_k^e)$, in which the mean ${\mathit{\mu }}_k^e$ and covariance ${\mathbf{C}}_k^e$ are given as [7]:

$$\begin{array}{cc}\hfill & {\mathit{\mu }}_k^e={\mathit{\mu }}_k^p+{\mathbf{C}}_{{\mathbf{x}}_k{\mathbf{y}}_k}{\left({\mathbf{C}}_k^y\right)}^{-1}({\mathbf{y}}_k-{\mathit{\mu }}_k^y),\hfill \end{array}$$

(15)

$$\begin{array}{cc}\hfill & {\mathbf{C}}_k^e={\mathbf{C}}_k^p-{\mathbf{C}}_{{\mathbf{x}}_k{\mathbf{y}}_k}{\left({\mathbf{C}}_k^y\right)}^{-1}{\mathbf{C}}_{{\mathbf{x}}_k{\mathbf{y}}_k}^{\mathrm{T}},\hfill \end{array}$$

(16)

where

$$\begin{array}{cc}\hfill & {\mathit{\mu }}_k^p=\mathbb{E}\left[{\mathbf{x}}_k\right|{\mathbf{y}}_{1:k-1}],\hfill \\ \hfill & {\mathbf{C}}_k^p=\mathrm{Cov}\left({\mathbf{x}}_k\right|{\mathbf{y}}_{1:k-1}),\hfill \\ \hfill & {\mathit{\mu }}_k^y=\mathbb{E}\left[{\mathbf{y}}_k\right|{\mathbf{y}}_{1:k-1}],\hfill \\ \hfill & {\mathbf{C}}_k^y=\mathrm{Cov}\left({\mathbf{y}}_k\right|{\mathbf{y}}_{1:k-1}),\hfill \\ \hfill & {\mathbf{C}}_{x_k{y}_k}=\mathrm{Cov}({\mathbf{x}}_k,{\mathbf{y}}_k|{\mathbf{y}}_{1:k-1}).\hfill \end{array}$$

(17)

Remark 1.

Since GP-ADF algorithm [7] is a state estimation method in the single sensor Gaussian process dynamic system, it is applicable to the multisensor system with all the sensor measurements stacked into a vector. Therefore, Equation (15) and Equation (16) are called the centralized estimation fusion method with GP-ADF here.

Then, let us present the prediction step and update step of the above fusion method in detail.

● Prediction Step

First, this step needs to use the previous posterior distribution $p\left({\mathbf{x}}_{k-1}\right|{\mathbf{y}}_{1:k-1})\approx \mathcal{N}({\mathit{\mu }}_{k-1}^e,{\mathbf{C}}_{k-1}^e)$ to produce a current predictive distribution $p\left({\mathbf{x}}_k\right|{\mathbf{y}}_{1:k-1})$. We write the predictive distribution as

$$p\left({\mathbf{x}}_k\right|{\mathbf{y}}_{1:k-1})=\int p\left({\mathbf{x}}_k\right|{\mathbf{x}}_{k-1})p\left({\mathbf{x}}_{k-1}\right|{\mathbf{y}}_{1:k-1})d{\mathbf{x}}_{k-1},$$

and $p\left({\mathbf{x}}_k\right|{\mathbf{x}}_{k-1})$ is a Gaussian distribution due to $h\sim \mathcal{GP}$ and ${\mathbf{w}}_k\sim \mathcal{N}(0,{\mathbf{\Sigma }}_{{\mathbf{w}}_k})$. From this, an analogy with Equation (7) yields the mean ${\mathit{\mu }}_k^p$ and the variance ${\mathbf{C}}_k^p$ of the state ${\mathbf{x}}_k$ conditioned on the measurement ${\mathbf{y}}_{1:k-1}$. In particular,

$$\begin{array}{cc}\hfill {\mathit{\mu }}_k^p& =\mathbb{E}\left({\mathbf{x}}_k\right|{\mathbf{y}}_{1:k-1})\hfill \\ \hfill & =\mathbb{E}\left[h\left({\mathbf{x}}_{k-1}\right)\right|{\mathbf{y}}_{1:k-1}]\hfill \\ \hfill & ={\mathbb{E}}_{{\mathbf{x}}_{k-1}}\left[{\mathbb{E}}_h\left(h\left({\mathbf{x}}_{k-1}\right)\right|{\mathbf{x}}_{k-1})\right|{\mathbf{y}}_{1:k-1}].\hfill \end{array}$$

(18)

${\mathit{\mu }}_k^p$ is a D-dimension mean vector and the computation of every target dimension can be seen in Equation (8). The corresponding covariance matrix

$$\begin{array}{cc}\hfill {\mathbf{C}}_k^p& =\mathrm{Cov}\left(h\left({\mathbf{x}}_{k-1}\right)\right|{\mathbf{y}}_{1:k-1})+\mathrm{Cov}\left({\mathbf{w}}_k\right)\hfill \\ \hfill & ={\mathbf{\Sigma }}_h+{\mathbf{\Sigma }}_{w_k},\hfill \end{array}$$

(19)

where ${\mathbf{\Sigma }}_{w_k}$ is the covariance matrix of transition noise obtained by the evidence maximization method, and the computation of ${\mathbf{\Sigma }}_h$ can be seen in Equation (10).

● Update Step

Next, we approximate the joint distribution $p({\mathbf{x}}_k,{\mathbf{y}}_k|{\mathbf{y}}_{1:k-1})$ with a joint Gaussian $\mathcal{N}({\mathit{\mu }}_k,{\mathbf{C}}_k)$, where

$${\mathit{\mu }}_k=\left[\begin{array}{c}{\mathit{\mu }}_k^p\\ {\mathit{\mu }}_k^y\end{array}\right],{\mathbf{C}}_k=\left[\begin{array}{cc}{\mathbf{C}}_k^p& {\mathbf{C}}_{{\mathbf{x}}_k{\mathbf{y}}_k}\\ {\mathbf{C}}_{{\mathbf{x}}_k{\mathbf{y}}_k}^{\mathrm{T}}& {\mathbf{C}}_k^y\end{array}\right].$$

Note that we are not aiming at approximating the distribution $p\left({\mathbf{x}}_k\right|{\mathbf{y}}_{1:k})$ directly. To obtain the joint approximation Gaussian distribution, we use the Gaussian distribution $\mathcal{N}({\mathit{\mu }}_k^y,{\mathbf{C}}_k^y)$ to approximate the distribution $p\left({\mathbf{y}}_k\right|{\mathbf{y}}_{1:k-1})$, which can be done in a same way to the prediction step due to

$$p\left({\mathbf{y}}_k\right|{\mathbf{y}}_{1:k-1})=\int p\left({\mathbf{y}}_k\right|{\mathbf{x}}_k)p\left({\mathbf{x}}_k\right|{\mathbf{y}}_{1:k-1})d{\mathbf{x}}_k.$$

On the other hand, the cross term satisfies

$$\begin{array}{cc}\hfill {\mathbf{C}}_{{\mathbf{x}}_k,{\mathbf{y}}_k}& ={\mathbb{E}}_{{\mathbf{x}}_k,g}\left[{\mathbf{x}}_k{\left({\mathbf{y}}_k\right)}^{\mathrm{T}}\right|{\mathbf{y}}_{1:k-1}]-{\mathit{\mu }}_k^p{\left({\mathit{\mu }}_k^y\right)}^{\mathrm{T}}.\hfill \end{array}$$

For the unknown term ${\mathbb{E}}_{{\mathbf{x}}_k,g}\left[{\mathbf{x}}_k{\left({\mathbf{y}}_k\right)}^{\mathrm{T}}\right|{\mathbf{y}}_{1:k-1}]$, we hav

$$\begin{array}{cc}\hfill & {\mathbb{E}}_{{\mathbf{x}}_k,g_{\tilde{a}}}\left[{\mathbf{x}}_k{\mathbf{y}}_k^{\tilde{a}}\right|{\mathbf{y}}_{1:k-1}]={\mathbb{E}}_{{\mathbf{x}}_k,g_{\tilde{a}}}\left[{\mathbf{x}}_k{g}_{\tilde{a}}\left({\mathbf{x}}_k\right)\right|{\mathbf{y}}_{1:k-1}]\hfill \\ \hfill =& \int {\mathbf{x}}_k\left[\sum _{i=1}^n{\beta }_i^{\tilde{a}}{k}_{g^{\tilde{a}}}({\mathbf{x}}_k,{\mathbf{x}}_i)\right]p\left({\mathbf{x}}_k\right)d{\mathbf{x}}_k\hfill \\ \hfill =& \sum _{i=1}^n{\beta }_i^{\tilde{a}}\int {\mathbf{x}}_k{c}_1\mathcal{N}\left({\mathbf{x}}_k\right|{\mathbf{x}}_i,{\Lambda }_{\tilde{a}})\mathcal{N}\left({\mathbf{x}}_k\right|{\mathit{\mu }}_k^p,{\mathbf{C}}_k^p)d{\mathbf{x}}_k\hfill \\ \hfill =& c_1{\left(c_2\right)}^{-1}\sum _{i=1}^n{\beta }_i^{\tilde{a}}\psi (x_i,{\mathit{\mu }}_k^p),\hfill \end{array}$$

(20)

where $\tilde{a}=1,2,\dots ,N\ast E$, $c_1^{-1}$ is the normalization constant of the unnormalized SE kernel, $\psi ({\mathbf{x}}_i,{\mathit{\mu }}_k^p)$ is the mean of a new Gaussian distribution that is the product of two Gaussian density functions and $c_2^{-1}$ is the normalization constant of the new Gaussian distribution. Consequently, from the joint distribution $p({\mathbf{x}}_k,{\mathbf{y}}_k|{\mathbf{y}}_{1:k-1})$, we obtain the posterior distribution $p\left({\mathbf{x}}_k\right|{\mathbf{y}}_{1:k})\approx \mathcal{N}({\mathit{\mu }}_k^e,{\mathbf{C}}_k^e)$ with the mean and variance as Equations (15) and (16), respectively.

Remark 2.

From Equations (15) and (16), we can see that the centralized fusion method with GP-ADF is similar to the unified optimal linear estimation fusion [40,54]. The difference is that the computational way of ${\mathit{\mu }}_k^p,{\mathbf{C}}_k^p,{\mathit{\mu }}_k^y$,${\mathbf{C}}_k^y$ and ${\mathbf{C}}_{{\mathbf{x}}_k{\mathbf{y}}_k}$ and they are mainly based on the training data due to the nonparametric dynamic system model.

(2) Fusion with GP-UKF

Following the single sensor GP-UKF (see e.g., References [14,26]), we present the prediction step and update step of the GP-UKF fusion.

● Prediction Step

At time $k-1$, based on the unscented transform [1], we obtain ${\mathcal{X}}_{k-1}$ containing $2D+1$ points through the mean ${\mathit{\mu }}_{k-1}^e$ and covariance ${\mathbf{C}}_{k-1}^e$. Using the GP prediction model, the transformed set is given by

$$\begin{array}{c}\hfill {\overline{\mathcal{X}}}_k^{\left[i\right]}={\mathcal{GP}}_h^{\left[\mu \right]}\left({\mathcal{X}}_{k-1}^{\left[i\right]}\right),fori=1,\dots ,2D+1,\end{array}$$

(21)

and the process noise is computed as

$$\begin{array}{c}\hfill {\mathbf{Q}}_k={\mathcal{GP}}_h^{[\Sigma ]}\left({\mathit{\mu }}_{k-1}^e\right),\end{array}$$

(22)

where ${\mathcal{GP}}_h^{\left[\mu \right]}$ and ${\mathcal{GP}}_h^{[\Sigma ]}$ are the mean and covariance models with respect to the Gaussian process ${\mathcal{GP}}_h$. The predictive mean and covariance are

$$\begin{array}{cc}\hfill {\mathit{\mu }}_k^p& =\sum _{i=1}^{2D+1}{\mathbf{W}}^{\left[i\right]}{\overline{\mathcal{X}}}_k^{\left[i\right]},\hfill \end{array}$$

(23)

$$\begin{array}{cc}\hfill {\mathbf{C}}_k^p& =\sum _{i=1}^{2D+1}{\mathbf{W}}^{\left[i\right]}({\overline{\mathcal{X}}}_k^{\left[i\right]}-{\mathit{\mu }}_k^p){({\overline{\mathcal{X}}}_k^{\left[i\right]}-{\mathit{\mu }}_k^p)}^{\mathrm{T}}+{\mathbf{Q}}_k,\hfill \end{array}$$

(24)

where ${\mathbf{W}}^{\left[i\right]}$ is the weight generated in the unscented transform. Using the predictive mean ${\mathit{\mu }}_k^p$ and covariance ${\mathbf{C}}_k^p$, we obtain ${\widehat{\mathcal{X}}}_k^{\left[i\right]}$ through the unscented transform. Thus, based on the GP observation model,

$$\begin{array}{c}\hfill {\widehat{\mathcal{Y}}}_k^{\left[i\right]}={\mathcal{GP}}_g^{\left[\mu \right]}\left({\widehat{\mathcal{X}}}_k^{\left[i\right]}\right)fori=1,\dots ,2D+1,\end{array}$$

(25)

and the observation noise matrix is determined by

$$\begin{array}{c}\hfill {\mathbf{R}}_k={\mathcal{GP}}_g^{[\Sigma ]}\left({\widehat{\mathcal{X}}}_k^{\left[i\right]}\right).\end{array}$$

(26)

Then the predicted observation and innovation covariance are calculated by

$$\begin{array}{cc}\hfill {\widehat{\mathbf{y}}}_k& =\sum _{i=1}^{2D+1}{\mathbf{W}}^{\left[i\right]}{\widehat{\mathcal{Y}}}_k^{\left[i\right]},\hfill \end{array}$$

(27)

$$\begin{array}{cc}\hfill {\mathbf{S}}_k& =\sum _{i=1}^{2D+1}{\mathbf{W}}^{\left[i\right]}({\widehat{\mathcal{Y}}}_k^{\left[i\right]}-{\widehat{\mathbf{y}}}_k){({\widehat{\mathcal{Y}}}_k^{\left[i\right]}-{\widehat{\mathbf{y}}}_k)}^{\mathrm{T}}+{\mathbf{R}}_k.\hfill \end{array}$$

(28)

Meanwhile, the cross covariance is given by

$$\begin{array}{c}\hfill {\mathbf{C}}_{{\mathbf{x}}_k{\mathbf{y}}_k}=\sum _{i=1}^{2D+1}{\mathbf{W}}^{\left[i\right]}({\widehat{\mathcal{X}}}_k^{\left[i\right]}-{\mathit{\mu }}_k^p){({\widehat{\mathcal{Y}}}_k^{\left[i\right]}-{\widehat{\mathbf{y}}}_k)}^{\mathrm{T}}\end{array}$$

(29)

● Update Step

Finally, the update can be obtained by the following equations:

$$\begin{array}{cc}\hfill {\mathit{\mu }}_k^e& ={\mathit{\mu }}_k^p+{\mathbf{C}}_{{\mathbf{x}}_k{\mathbf{y}}_k}{\mathbf{S}}_k^{-1}({\mathbf{y}}_k-{\widehat{\mathbf{y}}}_k),\hfill \end{array}$$

(30)

$$\begin{array}{cc}\hfill {\mathbf{C}}_k^e& ={\mathbf{C}}_k^p-{\mathbf{C}}_{{\mathbf{x}}_k{\mathbf{y}}_k}{\mathbf{S}}_k^{-1}{\mathbf{C}}_{{\mathbf{x}}_k{\mathbf{y}}_k}^{\mathrm{T}}.\hfill \end{array}$$

(31)

Remark 3.

Note that GP-ADF propagates the full Gaussian density by exploiting specific properties of Gaussian process models [7]. GP-UKF maps the sigma points through the Gaussian process models instead of the parametric functions and the distributions are described by finite-sample approximations. In contrast to GP-UKF, GP-ADF is consistent and moment preserving [7]. In addition GP-EKF requires linearization of the nonlinear prediction and observation models in order to propagate the state and observation, respectively [14]. GP-PF needs to perform one Gaussian process mean and variance computation per particle and has very high computational cost. For more details about these single sensor filtering methods, it can be seen in [7,14], and so forth. Due to the linearization loss of GP-EKF and high computational cost of GP-PF, we do not consider them and only introduce the multisensor estimation fusion with GP-ADF and GP-UKF in this paper.

3.2. Distributed Estimation Fusion

In this subsection, we mainly present the distributed estimation fusion with GP-ADF. For the dynamic systems (1) and (2), assume that there is a prior on initial state ${\mathbf{x}}_0$ and each local sensor sends its estimation to the fusion center. Thus, the posterior distribution of $m_{th}$ local state estimation is a Gaussian distribution $\mathcal{N}({\mathit{\mu }}_k^{em},{\mathbf{C}}_k^{em})$, where ${\mathit{\mu }}_k^{em}$ and ${\mathbf{C}}_k^{em}$ are calculated as follows

$$\begin{array}{cc}\hfill & {\mathit{\mu }}_k^{em}={\mathit{\mu }}_k^{pm}+{\mathbf{C}}_{{\mathbf{x}}_k{\mathbf{y}}_k^m}{\left({\mathbf{C}}_k^{ym}\right)}^{-1}({\mathbf{y}}_k^m-{\mathit{\mu }}_k^{ym}),\hfill \end{array}$$

(32)

$$\begin{array}{cc}\hfill & {\mathbf{C}}_k^{em}={\mathbf{C}}_k^{pm}-{\mathbf{C}}_{{\mathbf{x}}_k{\mathbf{y}}_k^m}{\left({\mathbf{C}}_k^{ym}\right)}^{-1}{\mathbf{C}}_{{\mathbf{x}}_k{\mathbf{y}}_k^m}^{\mathrm{T}},\hfill \end{array}$$

(33)

with

$$\begin{array}{cc}\hfill & {\mathit{\mu }}_k^{pm}=\mathbb{E}\left[{\mathbf{x}}_k\right|{\mathbf{y}}_{1:k-1}^m],\hfill \\ \hfill & {\mathbf{C}}_k^{pm}=\mathrm{Cov}\left({\mathbf{x}}_k\right|{\mathbf{y}}_{1:k-1}^m),\hfill \\ \hfill & {\mathit{\mu }}_k^{ym}=\mathbb{E}\left[{\mathbf{y}}_k^m\right|{\mathbf{y}}_{1:k-1}^m],\hfill \\ \hfill & {\mathbf{C}}_k^{ym}=\mathrm{Cov}\left({\mathbf{y}}_k^m\right|{\mathbf{y}}_{1:k-1}^m),\hfill \\ \hfill & {\mathbf{C}}_{x_k{y}_k^m}=\mathrm{Cov}({\mathbf{x}}_k,{\mathbf{y}}_k^m|{\mathbf{y}}_{1:k-1}^m).\hfill \end{array}$$

(34)

Then, the local posterior mean and covariance are transmitted to the fusion center to yield the posterior distribution of global state estimation.

In general, the centralized estimation fusion method with directly fusing raw data has better fusion performance than the distributed estimation fusion method based on the processed data. However, in some cases, the distributed estimation fusion is equivalent to the centralized estimation fusion. In particular, in what follows, we propose a distributed estimation fusion method and prove that the distributed estimation fusion method is equivalent to the above centralized estimation fusion method when all cross terms ${\mathbf{C}}_{{\mathbf{x}}_k{\mathbf{y}}_k^m}$, $m=1,\dots ,N$, have full column rank.

Proposition 1.

Assume that all cross terms ${\mathbf{C}}_{{\mathbf{x}}_k{\mathbf{y}}_k^m}$, $m=1,\dots ,N$ are full column rank, the distributed estimation fusion is equivalent to the centralized estimation fusion as follows

$$\begin{array}{cc}\hfill {\mathit{\mu }}_k^e& ={\mathit{\mu }}_k^p+{\mathbf{C}}_{{\mathbf{x}}_k{\mathbf{y}}_k}\sum _{m=1}^N{\left({\mathbf{C}}_k^y\right)}^{-1}\left(m\right)({\mathit{\mu }}_k^{ym}-{\mathit{\mu }}_k^y\left(m\right))\hfill \\ \hfill & +{\mathbf{C}}_{{\mathbf{x}}_k{\mathbf{y}}_k}\sum _{m=1}^N{\left({\mathbf{C}}_k^y\right)}^{-1}\left(m\right){\mathbf{C}}_k^{ym}{\mathbf{C}}_{{\mathbf{x}}_k{\mathbf{y}}_k^m}^{+}({\mathit{\mu }}_k^{em}-{\mathit{\mu }}_k^{pm}),\hfill \end{array}$$

(35)

where

$$\begin{array}{cc}\hfill & {\mathit{\mu }}_k^y=({\mathit{\mu }}_k^y\left(1\right),\dots ,{\mathit{\mu }}_k^y\left(N\right)),\hfill \\ \hfill & {\mathit{\mu }}_k^y\left(m\right)=\mathbb{E}\left({\mathbf{y}}_k^m\right|{\mathbf{y}}_{1:k-1}),\hfill \end{array}$$

and

$${\left({\mathbf{C}}_k^y\right)}^{-1}=[{\left({\mathbf{C}}_k^y\right)}^{-1}\left(1\right),{\left({\mathbf{C}}_k^y\right)}^{-1}\left(2\right),\dots ,{\left({\mathbf{C}}_k^y\right)}^{-1}\left(N\right)]$$

is an appropriate partition of matrix ${\left({\mathbf{C}}_k^y\right)}^{-1}$ such that Equation (35) holds. The superscript “+” stands for pseudoinverse.

Proof. 

According to Equation (15), the mean of centralized fused state is given by

$$\begin{array}{cc}\hfill {\mathit{\mu }}_k^e& ={\mathit{\mu }}_k^p+{\mathbf{C}}_{{\mathbf{x}}_k{\mathbf{y}}_k}{\left({\mathbf{C}}_k^y\right)}^{-1}({\mathbf{y}}_k-{\mathit{\mu }}_k^y)\hfill \end{array}$$

(36)

$$\begin{array}{cc}\hfill & ={\mathit{\mu }}_k^p-{\mathbf{C}}_{{\mathbf{x}}_k{\mathbf{y}}_k}{\left({\mathbf{C}}_k^y\right)}^{-1}{\mathit{\mu }}_k^y+{\mathbf{C}}_{{\mathbf{x}}_k{\mathbf{y}}_k}{\left({\mathbf{C}}_k^y\right)}^{-1}{\mathbf{y}}_k\hfill \end{array}$$

(37)

$$\begin{array}{cc}\hfill & ={\mathit{\mu }}_k^p-{\mathbf{C}}_{{\mathbf{x}}_k{\mathbf{y}}_k}{\left({\mathbf{C}}_k^y\right)}^{-1}{\mathit{\mu }}_k^y+{\mathbf{C}}_{{\mathbf{x}}_k{\mathbf{y}}_k}\sum _{m=1}^N{\left({\mathbf{C}}_k^y\right)}^{-1}\left(m\right){\mathbf{y}}_k^m.\hfill \end{array}$$

(38)

Based on the assumption that all $C_{x_k{y}_k^m}$ have full column rank, we obtain

$$\begin{array}{c}\hfill {\mathbf{C}}_{{\mathbf{x}}_k{\mathbf{y}}_k^m}^{+}{\mathbf{C}}_{{\mathbf{x}}_k{\mathbf{y}}_k^m}=\mathbf{I},m=1,\dots ,N.\end{array}$$

(39)

Thus, combining Equations (37) and (38) and Equation (39) yields

$$\begin{array}{cc}\hfill & {\mathbf{C}}_{{\mathbf{x}}_k{\mathbf{y}}_k}{\left({\mathbf{C}}_k^y\right)}^{-1}{\mathbf{y}}_k\hfill \\ \hfill =& {\mathbf{C}}_{{\mathbf{x}}_k{\mathbf{y}}_k}\sum _{m=1}^N{\left({\mathbf{C}}_k^y\right)}^{-1}\left(m\right){\mathbf{y}}_k^m\hfill \\ \hfill =& {\mathbf{C}}_{{\mathbf{x}}_k{\mathbf{y}}_k}\sum _{m=1}^N{\left({\mathbf{C}}_k^y\right)}^{-1}\left(m\right){\mathbf{C}}_k^{ym}{\mathbf{C}}_{{\mathbf{x}}_k{\mathbf{y}}_k^m}^{+}{\mathbf{C}}_{{\mathbf{x}}_k{\mathbf{y}}_k^m}{\left({\mathbf{C}}_k^{ym}\right)}^{-1}{\mathbf{y}}_k^m.\hfill \end{array}$$

(40)

In order to obtain the centralized estimation mean, that is, fuse the mean of local sensor state estimation at the fusion center, we use Equations (32) and (40) to eliminate ${\mathbf{y}}_k$ from Equation (36). From Equation (32), we have

$$\begin{array}{cc}\hfill & {\mathbf{C}}_{{\mathbf{x}}_k{\mathbf{y}}_k^m}{\left({\mathbf{C}}_k^{ym}\right)}^{-1}{\mathbf{y}}_k^m\hfill \\ \hfill =& {\mathit{\mu }}_k^{em}-{\mathit{\mu }}_k^{pm}+{\mathbf{C}}_{{\mathbf{x}}_k{\mathbf{y}}_k^m}{\left({\mathbf{C}}_k^{ym}\right)}^{-1}{\mathit{\mu }}_k^{ym}.\hfill \end{array}$$

(41)

Thus, substituting Equation (41) into Equation (40) and recalling Equations (36)–(38), we have

$$\begin{array}{cc}\hfill {\mathit{\mu }}_k^e& ={\mathit{\mu }}_k^p-{\mathbf{C}}_{{\mathbf{x}}_k{\mathbf{y}}_k}{\left({\mathbf{C}}_k^y\right)}^{-1}{\mathit{\mu }}_k^y\hfill \\ \hfill & +{\mathbf{C}}_{{\mathbf{x}}_k{\mathbf{y}}_k}\sum _{m=1}^N{\left({\mathbf{C}}_k^y\right)}^{-1}\left(m\right){\mathbf{C}}_k^{ym}{\mathbf{C}}_{{\mathbf{x}}_k{\mathbf{y}}_k^m}^{+}\hfill \\ \hfill & \times ({\mathit{\mu }}_k^{em}-{\mathit{\mu }}_k^{pm}+{\mathbf{C}}_{{\mathbf{x}}_k{\mathbf{y}}_k^m}{\left({\mathbf{C}}_k^{ym}\right)}^{-1}{\mathit{\mu }}_k^{ym})\hfill \\ \hfill & ={\mathit{\mu }}_k^p+{\mathbf{C}}_{{\mathbf{x}}_k{\mathbf{y}}_k}\sum _{m=1}^N{\left({\mathbf{C}}_k^y\right)}^{-1}\left(m\right)({\mathit{\mu }}_k^{ym}-{\mathit{\mu }}_k^y\left(m\right))\hfill \\ \hfill & +{\mathbf{C}}_{{\mathbf{x}}_k{\mathbf{y}}_k}\sum _{m=1}^N{\left({\mathbf{C}}_k^y\right)}^{-1}\left(m\right){\mathbf{C}}_k^{ym}{\mathbf{C}}_{{\mathbf{x}}_k{\mathbf{y}}_k^m}^{+}({\mathit{\mu }}_k^{em}-{\mathit{\mu }}_k^{pm}).\hfill \end{array}$$

Therefore, we obtain the state mean of the distributed estimation fusion. □

Remark 4.

The proposed distributed estimation fusion formula is equivalent to the centralized estimation fusion, therefore, the two fusion methods have the same fusion performance. However, distributed estimation fusion is more desirable in real applications in terms of reliability, survivability, communication bandwidth and computational burdens. From Equation (35), we can see that ${\mathit{\mu }}_k^p$, ${\mathbf{C}}_k^p$, ${\mathit{\mu }}_k^y$, ${\mathbf{C}}_k^y$, ${\mathbf{C}}_{{\mathbf{x}}_k{\mathbf{y}}_k}$ can be calculated in advance and the terms ${\mathit{\mu }}_k^{pm}$, ${\mathit{\mu }}_k^{ym}$, ${\mathit{\mu }}_k^{em}$, ${\mathbf{C}}_k^{ym}$, ${\mathbf{C}}_{{\mathbf{x}}_k{\mathbf{y}}_k^m}$ need to be transmitted by local sensors in real time.

Remark 5.

From the update step (30) of the centralized estimation fusion with GP-UKF, the distributed estimation fusion method (35) is also suitable for the GP-UKF. Since the detailed description and proof for the distributed GP-UKF fusion are similar to Proposition 1, we omit the results. Note that the distributed GP-UKF fusion is different from the distributed GP-ADF fusion in the prediction step.

4. Numerical Examples

In this section, we assess the performance of our fusion methods with two different examples.

4.1. 1D Example

We consider the nonlinear dynamic system with two sensor measurements as follows:

$$\begin{array}{cc}\hfill & {\mathbf{x}}_{k+1}=0.5{\mathbf{x}}_k+\frac{25{\mathbf{x}}_k}{1+{\mathbf{x}}_k^2}+{\mathbf{w}}_k,{\mathbf{w}}_k\sim \mathcal{N}(0,{\sigma }^2)\hfill \end{array}$$

(42)

$$\begin{array}{cc}\hfill & {\mathbf{y}}_k^m=5sin(2{\mathbf{x}}_k+{\mathbf{z}}^m)+{\mathbf{v}}_k^m,{\mathbf{v}}_k^m\sim \mathcal{N}(0,0.1^2),m=1,2,\hfill \end{array}$$

(43)

where ${\mathbf{z}}^m,m=1,2$ are given by ${\mathbf{z}}^1=0$ and ${\mathbf{z}}^2=\frac{\pi }{4}$, respectively. This nonlinear system is popularly used in many papers (see References [7,10]) to measure the performance of nonlinear filters. We regard the first 200 samples as the training data and use the rest 50 samples called testing data to test the fusion methods. Assume that the ${\mathbf{x}}_{-200}$ is a Gaussian distribution with prior mean ${\mu }_{-200}=0$ and variance ${\sigma }_{-200}^2=1.5^2$. The estimation performance of the single sensor and multisensor fusion is evaluated by the average Mahalanobis distances that are also used in Reference [7]. The Mahalanobis distance, which is defined between the ground truth and the filtered mean, is as follows:

$$\begin{array}{c}\hfill \mathrm{Maha}=\sqrt{{({\mathbf{x}}_k-{\widehat{\mathbf{x}}}_k)}^{\mathrm{T}}{\left({\mathbf{C}}_k^e\right)}^{-1}({\mathbf{x}}_k-{\widehat{\mathbf{x}}}_k)},\end{array}$$

(44)

where ${\mathbf{C}}_k^e$ is the estimated covariance. For the Mahalanobis distance, lower values indicate better performance [7]. Figure 2 and Figure 3 show the average Mahalanobis distances of the single filter with GP-ADF and GP-UKF, and the fusion methods after 1000 independent runs. Figure 4 shows the average Mahalanobis distances for different system noises.

According to the average Mahalanobis distances in Figure 2 and Figure 3, we can see that the estimation performance of Sensor 2 is better than that of Sensor 1. Furthermore, the multisensor estimation fusion methods perform better than the single sensor filters with GP-ADF and GP-UKF, respectively. Thus the effectiveness and advantages of the multisensor estimation fusion with Gaussian process can be observed. In addition, GP-ADF fusion performs better than GP-UKF fusion for the system noise ${\sigma }^2=1$ and GP-UKF fusion performs better than GP-ADF fusion for the system noise ${\sigma }^2=2$. It means that GP-ADF fusion and GP-UKF fusion have their comparative advantages for different system noises. It suggests that GP-ADF fusion may be a better choice for the small system noise and GP-UKF fusion may be worth considering for the relatively large system noise. From Figure 4, we can see that the proposed fusion methods enjoy better performance for different system noises. It demonstrates the consistence of our fusion methods for different system noises.

4.2. 2D Nonlinear Dynamic System

In this subsection, we elaborate the performance of the fusion methods with GP-ADF and GP-UKF, and we consider the example with constant turn motion model in target tracking. The root mean square error (RMSE) is used as the estimation performance measure. It is defined as follows:

$$\begin{array}{c}\hfill {\mathrm{RMSE}}_k=\sqrt{\frac{1}{M}\sum _{j=1}^M{\parallel {\mathbf{x}}_k^j-{\widehat{\mathbf{x}}}_k^j\parallel }_2^2},\end{array}$$

(45)

where M is the total number of simulation runs, ${\mathbf{x}}_k^j$ is the true simulated state for the jth simulation, and ${\widehat{\mathbf{x}}}_k^j$ is the estimated state value for the jth simulation at time k, $j=1,\dots ,M$. The lower the value of the RMSE is, the better the performance of corresponding method is.

4.2.1. Experimental Setup

We consider the constant turn motion model with two sensors as follows:

$$\begin{array}{cc}\hfill & {\mathbf{x}}_{k+1}=\mathbf{F}{\mathbf{x}}_k+{\mathbf{w}}_k,{\mathbf{w}}_k\sim \mathcal{N}(0,{\mathbf{Q}}_k),\hfill \\ \hfill & {\mathbf{y}}_k^m=\left[\begin{array}{c}\sqrt{{({\mathbf{x}}_k\left(1\right)-{\mathbf{z}}^m\left(1\right))}^2+{({\mathbf{x}}_k\left(2\right)-{\mathbf{z}}^m\left(2\right))}^2}\\ arctan\frac{{\mathbf{x}}_k\left(2\right)-{\mathbf{z}}^m\left(2\right)}{{\mathbf{x}}_k\left(1\right)-{\mathbf{z}}^m\left(1\right)}\end{array}\right]+{\mathbf{v}}_k^m,\hfill \end{array}$$

(46)

$$\begin{array}{cc}\hfill & {\mathbf{v}}_k^m\sim \mathcal{N}(0,{\mathbf{R}}_k),m=1,2,\hfill \end{array}$$

(47)

where the state transition matrix $\mathbf{F}$ is given by

$$\mathbf{F}=\left[\begin{array}{cc}\cos \Omega & \sin \Omega \\ -\sin \Omega & \cos \Omega \end{array}\right]$$

with angular velocity $\Omega$, the covariance matrix of transition noise satisfies

$${\mathbf{Q}}_k=\left[\begin{array}{cc}5& 0\\ 0& 5\end{array}\right],$$

and the covariance matrix of measurement noise is

$${\mathbf{R}}_k=\left[\begin{array}{cc}{5}^2& 0\\ 0& {\left(\frac{0.1\ast \pi }{180}\right)}^2\end{array}\right].$$

The sensor positions are given by ${\mathbf{z}}^1={[-200,200]}^{\mathrm{T}}$ and ${\mathbf{z}}^2={[200,-200]}^{\mathrm{T}}$, respectively.

We use the transition function (46) and measurement functions (47) to simulate a group of data. The values of angular velocity $\Omega$ are given by $\Omega =\frac{2\pi }{90}\mathrm{rad}/\mathrm{s}$ and $\Omega =0.5\mathrm{rad}/\mathrm{s}$, respectively, which correspond to the small turn motion model and large turn motion model. The first $\kappa$ samples are regarded as the training data and the rest 50 samples called testing data are used to test the fusion methods. The ${\{{\mathbf{x}}_{\tau },{\mathbf{y}}_{\tau }^1,{\mathbf{y}}_{\tau }^2\}}_{\tau =-\kappa }^{-1}$ are the true simulated state and two-sensor observations used to train the models. Assume that the ${\mathbf{x}}_{-\kappa }$ is a Gaussian distribution with the prior mean ${\mathit{\mu }}_{-\kappa }={[0,0]}^{\mathrm{T}}$ and the identity matrix variance ${\mathbf{S}}_{-\kappa }=\mathbf{I}$. The initial value of filtering for each sensor is set as the Cartesian coordinate transformation value based on the Polar coordinate observation with a random perturbation for each dimension. We use the random perturbation $N(10,1)$ for the $\kappa =300$ training data case and $N(20,1)$ for the $\kappa =30,60$ training data cases, respectively. The initial value of fusion filtering is the mean of the initial values of all sensors. Based on the available observation information, the outputs of GP-ADF and GP-UKF for each single sensor are the mean and covariance terms ${\mathit{\mu }}_k^{pm}$, ${\mathbf{C}}_k^{pm}$, ${\mathit{\mu }}_k^{em}$, ${\mathbf{C}}_k^{em}$, ${\mathit{\mu }}_k^{ym}$,${\mathbf{C}}_k^{ym}$, ${\mathbf{C}}_{{\mathbf{x}}_k{\mathbf{y}}_k^m}$. The distributed fusion methods, the RCC-CI algorithm [37] and the convex combination method [43], are also used as a comparison with our distributed fusion method. The RCC-CI algorithm and the convex combination method directly use the estimated mean and covariance matrix to fuse the estimated results. Our methods take full advantages of the results of the single sensor Gaussian process filters such as the cross term ${\mathbf{C}}_{{\mathbf{x}}_k{\mathbf{y}}_k^m}$ which is easily obtained. Note that all of the fusion methods are based on the local estimates that are distilled from the measurements. Thus, the comparison is fair from the available observation information. The simulation results of the RCC-CI algorithm are based on the YALMIP toolbox [55]. We compare the fusion methods with the RCC-CI algorithm and the convex combination method by RMSE after $M=1000$ independent simulation runs. Meanwhile, we also compare the computation time of the multisensor estimation fusion methods with the RCC-CI algorithm and the convex combination method. The ratio of full rank, defined as the percentage of full column rank of the cross term ${\mathbf{C}}_{{\mathbf{x}}_k{\mathbf{y}}_k^m}$ for the single sensor filtering at every time step after $M=1000$ independent runs, is used to test the condition of equivalence. If the ratio of full rank is less than $100\%$, the condition is broken. The simulations are done under Matlab (MathWorks, Inc., Natick, MA, USA) R2018b with ThinkPad W540.

Figure 5 describes the ratios of full column rank in the single sensor case for testing data with $\kappa =300$ training data and $\Omega =\frac{2\pi }{90},0.5\mathrm{rad}/\mathrm{s}$. The RMSEs for $\Omega =\frac{2\pi }{90},0.5\mathrm{rad}/\mathrm{s}$ in the case of $\kappa =300$ training data are depicted in Figure 6 and Figure 7, respectively. The average computation time of these fusion methods for $\Omega =\frac{2\pi }{90},0.5\mathrm{rad}/\mathrm{s}$ in the case of $\kappa =300$ training data after 1000 independent simulation runs are depicted in Figure 8 and Figure 9, respectively. The diamond line represents the RMSE of the centralized estimation fusion with GP-ADF (CMF-GP-ADF) and the star line represents the distributed estimation fusion with GP-ADF (DMF-GP-ADF). The circle line represents the RMSE of the centralized estimation fusion with GP-UKF (CMF-GP-UKF) and the x line represents the distributed estimation fusion with GP-UKF (DMF-GP-UKF). The upper triangle line represents the RMSE of the RCC-CI fusion algorithm with GP-ADF (RCC-CI-GP-ADF) and the lower triangle line represents the RMSE of the RCC-CI fusion algorithm with GP-UKF (RCC-CI-GP-UKF). The square line represents the RMSE of the convex combination method, that is, covariance weight, with GP-ADF (CW-GP-ADF) and the + line represents the RMSE of the convex combination method with GP-UKF (CW-GP-UKF). At the same time, solid line is used for the single sensor GP-ADF filter and dash-dotted line is used for the single GP-UKF filter. Similarly, the ratios of full column rank in the single sensor case for testing data with $\kappa =30,60$ training data and $\Omega =\frac{2\pi }{90}\mathrm{rad}/\mathrm{s}$ are depicted in Figure 10 and Figure 11, respectively. The RMSEs for $\Omega =\frac{2\pi }{90}\mathrm{rad}/\mathrm{s}$ in the case of $\kappa =30,60$ training data are depicted in Figure 12 and Figure 13, respectively. The average computation time of these fusion methods for testing data with $\kappa =30,60$ training data and $\Omega =\frac{2\pi }{90}\mathrm{rad}/\mathrm{s}$ after 1000 independent simulation runs are depicted in Figure 14 and Figure 15, respectively.

4.2.2. Experimental Analysis

We divide the experimental analysis into two cases, that is, the equivalence condition is satisfied or not. Some phenomena and analysis are described as follows:

Case 1: the equivalence condition is satisfied.

• From Figure 5, we can see that the ratios of full rank are all equal to $100\%$. It means that the cross terms of the single sensor filters are all full column rank for the $\kappa =300$ training data case and thus the equivalence condition of centralized fusion and distributed fusion is satisfied in Proposition 1. Meanwhile, from Figure 6 and Figure 7, the RMSE of the distributed estimation fusion is the same as that of the centralized estimation fusion based on GP-ADF and GP-UKF, respectively. It demonstrates the equivalence between the centralized estimation fusion and the distributed estimation fusion in Proposition 1 under the condition of full column rank.
• From Figure 6 and Figure 7, it can be seen that the RMSE of the multisensor estimation fusion method is lower than that of the single sensor filtering. It shows that the multisensor fusion improves the estimation accuracy.
• In addition, the RMSE of multisensor estimation fusion method is lower than that of the RCC-CI algorithm and the convex combination method. It implies the effectiveness of the fusion methods based on Gaussian processes. The possible reason is that our estimation fusion methods extract more extra correlation information, and the RCC-CI algorithm and the convex combination method only use the local estimates with mean and covariance.
• Compared Figure 6 with Figure 7, we can find that our methods do well in the different angular velocity cases, i.e., the small turn motion model and large turn motion model. It confirms that our fusion methods have fine applicability with better performance.
• From Figure 8 and Figure 9, we can find that the computation time of the distributed estimation fusion is less than that of the centralized estimation fusion. It demonstrates the superiority of the distributed estimation fusion under the same fusion performance. The computation time of the proposed fusion methods is much less than that of the RCC-CI algorithm. The possible reason is due to solve a optimization problem for the RCC-CI algorithm. The convex combination method takes the least computation time, since it directly uses the weight combination with covariance.

Case 2: the equivalence condition is not satisfied.

• From Figure 10 and Figure 11, it can be seen that the ratios of full rank are both less than $100\%$ for the single sensor GP-UKF case with $\kappa =30$ and $\kappa =60$ training data. Thus, the equivalence between the centralized and distributed estimation fusion with GP-UKF is broken in Figure 12 and Figure 13, respectively. The reason may be that the Gaussian models are relatively inaccurate with less training data, which can be known from the comparison with Figure 6, Figure 12 and Figure 13 for the same methods. At the same time, the finite-sample approximation of GP-UKF seriously depends on the Gaussian process models and the computation way of the cross terms is the sum about rank-one matrices for GP-UKF. However, the equivalence is still satisfied for the GP-ADF fusion. It implies GP-ADF is more stable than GP-UKF, which is also referred to in Reference [7].
• We can also see that the performance of GP-UKF fusion is better than that of GP-ADF fusion with $\kappa =300$ training data from Figure 6 and a little worse with $\kappa =30,60$ training data from Figure 12 and Figure 13. Meanwhile, from Figure 8, Figure 14 and Figure 15, the average computation time of GP-UKF fusion is less than that of GP-ADF fusion with $\kappa =300$ training data and is contrary with $\kappa =30,60$ training data. It may inspire us that the GP-UKF fusion is suitable for the enough training data case and GP-ADF fusion does well in the small number of training data case for the turn motion systems.

In a word, distributed estimation fusion with GP-ADF is more stable, and has a better performance in the case of the small number of training data. If we have enough training data, GP-UKF fusion may be the better choice.

5. Conclusions

In the context of this paper, the property matters if the real system does not closely follow idealized models or a parametric model cannot easily be determined for nonlinear systems. A non-parametric method, the Gaussian process, is introduced to learn models from training data, since it takes both model uncertainty and sensor measurement noise into account. In order to estimate the state of the multisensor nonlinear dynamic systems, we have used the Gaussian process models as the prior of the transition and measurement function of the dynamic system. Then the transition function and measurement function have been trained with Gaussian processes, respectively. Based on the Gaussian process models for all available local sensor measurements, we have developed two fusion methods, centralized estimation fusion and distributed estimation fusion, with GP-ADF and GP-UKF, respectively. Taking full advantage of the nature of the Gaussian process, the equivalence between centralized estimation fusion and distributed estimation fusion has been derived under mild conditions. Simulations show that the equivalence is satisfied under given conditions and the multisensor estimation fusion performs better than the single sensor filters. Compared with the RCC-CI algorithm, the multisensor estimation fusion methods not only have higher accuracy, but also require less computation time. The estimation performance of multisensor estimation fusion methods is also better than that of the convex combination method. Future work may involve state initialization, Gaussian process latent variable models without ground truth states, multiple model fusion with different Gaussian processes for maneuvering target tracking and validate the methods with the real sensor data.