Distributed hybrid consensus–based square-root cubature quadrature information filter and its application to maneuvering target tracking

Abstract

To handle nonlinear filtering problems with networked sensors in a distributed manner, a novel distributed hybrid consensus–based square-root cubature quadrature information filter is proposed. The proposed hybrid consensus–based square-root cubature quadrature information filter exploits fifth-order spherical simplex-radial quadrature rule to tackle system nonlinearities and incorporates a novel measurement update strategy into the hybrid consensus filtering framework, which takes the predicted measurement error into account and hence produces more accurate estimates. In addition, the proposed hybrid consensus–based square-root cubature quadrature information filter inherits the complementary positive features of both consensus on information and consensus on measurements methods and avoids sensitive matrix operations such as square-root decompositions and inversion of covariances, which is beneficial for numerical stability. Stability analysis with respect to consensus, convergence, and consistency for the proposed hybrid consensus–based square-root cubature quadrature information filter is also developed. The effectiveness of the proposed hybrid consensus–based square-root cubature quadrature information filter is validated through a maneuvering target tracking scenario. The simulation results show that the proposed hybrid consensus–based square-root cubature quadrature information filter outperforms the existing algorithms at the expense of a slight increase in computational cost.

Introduction

Distributed nonlinear filtering (DNF) algorithms have attracted enormous popularity in recent years. Unlike the centralized solutions, no fusion center is required for distributed filtering, and each node only communicates with its immediate neighbors. Due to its scalability, low communication cost, and robustness to sensor failures, the consensus-based strategies have been extensively studied in the distributed filtering community. In the consensus-based estimation framework, through iteratively exchanging estimated information with its immediate neighbors, each node is able to achieve consensus on final estimates and exhibit comparable performance to its centralized counterpart.

In the linear setting, the well-known Kalman consensus filter (KCF)^1,2 updates the posterior estimate with average consensus performed on the prior estimates from neighboring nodes. KCF works well under the condition that all nodes can observe the target of interest. However, it performs poorly in the presence of naive nodes, where the target is neither observed by the node itself nor observed by its immediate neighbors. ³ For nonlinear filtering problems, based on the information weighted consensus filter (IWCF)^3,4 and extended Kalman filter, the extended IWCF (EIWCF) ⁵ was proposed. However, EIWCF suffers from low accuracy and poor stability due to neglect of higher order terms in the Taylor expansion. In such a case, the unscented Kaman filter (UKF), which is a kind of deterministic sampling method, is a better choice. ⁶ By applying consensus on local estimates from neighboring nodes, Li et al. ⁷ proposed a weighted average consensus-based unscented Kaman filter (WACUKF) to estimate the true state, and its estimation error has been proved to be bounded in mean square sense. However, since WACUKF directly performs consensus on local estimates, the global optimality cannot be achieved. With two parallel consensus protocols performed on both prior distributions and likelihoods, the consensus on likelihoods and priors approach based on UKF (CLCP-UKF) was proposed to deal with nonlinear filtering problems with naive nodes, which guarantees both estimation accuracy and filtering stability. ⁸ To solve filtering problems with modeling uncertainties, a distributed multiple model UKF is derived, which computes the mode-conditioned estimates and mode probabilities with consensus strategies. ⁹ In the above-mentioned algorithms, UKF acts as the underlying filter to solve system nonlinearities. However, an appropriate scaling parameter is required for UKF to achieve acceptable performance. A negative weight in UKF may cause the loss of positive definiteness for the covariance matrix, especially in high dimensional systems, which may make UKF stop its operation.^10,11

Similar to UKF, the cubature Kalman filter (CKF) is also a deterministic sampling method to settle nonlinearities, but there are no requirements for scaling parameters. ¹⁰ Only dimension of the state vector is needed to compute the cubature points, which can be carried out in advance, and hence eases the burden on curse of dimensionality. In the light of consensus on information (CI), ¹² the distributed cubature information filter (DCIF) was developed. ¹³ However, DCIF ignores the problem of information redundancy. By embedding CKF into the IWCF framework, the cubature information consensus filter (Cub-ICF) was developed for space object tracking with space-based object (SBO) measurements, which effectively overcomes the drawbacks of DCIF. ¹⁴ To avoid square-root operations on covariance matrices and improve numerical stability, the square-root version cubature IWCF (SCICF) was proposed. ¹⁵ However, in actual applications, only a small number of consensus iterations is allowed, and hence, IWCF-based algorithms may lose the property of consistency. ¹⁶ To solve this problem, the hybrid consensus–based cubature Kalman filter (HCCKF) was developed by combining CI and consensus on measurements (CM) together, which takes both accuracy and consistency into account. ¹⁷ To tackle the problem of measurement uncertainty, the variational Bayesian consensus CKF (VB-CCKF) was proposed, where VB inference ¹⁸ is exploited to iteratively approximate sufficient statistics of the measurement noise covariance. ¹⁹ However, the algorithms discussed above are all designed based on the statistical linear error propagation methodology, ²⁰ which is exploited to compute the pseudo measurement matrix. Since the predicted measurement errors are ignored, the corresponding information filters are not equal to its Kalman filtering counterparts.^21,22 With approximate compensation for the linearization errors, the hybrid CMCI based on CKF (HCMCI-CKF) was developed for collaborative state estimation, which effectively improves the estimation accuracy. ²³

Motivated by the above research, a distributed hybrid consensus–based square-root cubature quadrature information filter (SRCQIF-HC) is proposed. The main contributions include the following: (1) the high-degree cubature quadrature rule is combined with the square-root information filter to develop the proposed SRCQIF-HC, which takes the linearization error of measurement equation into consideration; (2) consensus and convergence of the proposed SRCQIF-HC to the centralized scheme is discussed and analyzed; and (3) stability analysis in terms of consistency of local estimates is facilitated by introducing pseudo transition and measurement matrices. The remainder of this article is organized as follows. Section “Problem formulation” formulates the problem of DNF with networked sensors. A novel Gaussian approximation information filtering paradigm with consideration for the predicted measurement errors is developed in section “A novel Gaussian approximation information filtering paradigm.” In section “High-degree spherical simplex-radial quadrature rule,” the high-degree spherical simplex-radial quadrature rule based on Mysovskikh’s method is derived. The proposed SRCQIF-HC and analysis with regard to consensus, convergence, and consistency of local estimates are presented in section “Distributed SRCQIF-HC.” Simulation experiments are given in section “Simulation experiments.”

Problem formulation

Consider a sensor network consisting of $N$ nodes in the surveillance area. In the network, there are two types of nodes: communication nodes and sensor nodes. Communication nodes are only capable of processing and exchanging local information with its neighbors, while sensor nodes are also able to sense information from the environment. It is worth mentioning that communication nodes serve as relays of information among distant nodes, which is beneficial for improving network connectivity as well as reducing the requirements of transmission power of each node, and in turn extending lifetime of the sensor network.

The communication topology in the sensor network is denoted by an undirected graph $\mathcal{G}=(\mathcal{N},\mathcal{E})$ . Here, $\mathcal{N}=\mathcal{S}\cup \mathcal{C}=\{1,2,\dots ,N\}$ is the set of nodes, where $\mathcal{S}$ and $\mathcal{C}$ represent sensor nodes and communication nodes, respectively. $\mathcal{E}$ is the set of edges such that $(i,j)\in \mathcal{E}$ if node $i$ is able to broadcast or receive information from node $j$ , and vice versa. The set of neighbors of node $i$ is denoted by ${\mathcal{N}}_i=\{j|(i,j)\in \mathcal{E},\forall j\ne i\}$ . Further, ${\mathcal{J}}_i={\mathcal{N}}_i\cup \left\{i\right\}$ denotes the in-neighbors (including node $i$ itself). If the target of interest is not observed by the in-neighbors of node $i$ (including node $i$ itself), then node $i$ is called naive node.^3,24–26

The DNF problem with $N_s$ sensors can be formulated as follows. Consider a discrete-time nonlinear dynamic system

$${\mathit{x}}_k=\mathit{f}\left({\mathit{x}}_{k-1}\right)+{\mathit{w}}_{k-1}$$

(1)

and a set of sensor nodes $\mathcal{S}$ with local measurements

$${\mathit{z}}_{i,k}={\mathit{h}}_i\left({\mathit{x}}_k\right)+{\mathit{v}}_{i,k},i\in \mathcal{S}$$

(2)

where ${\mathit{x}}_k\in {\mathbb{R}}^{n_x}$ and ${\mathit{z}}_{i,k}\in {\mathbb{R}}^{n_z}$ are the state and measurement at time instant $k$ and $\mathit{f}(·)$ and ${\mathit{h}}_i(·)$ are the nonlinear transition and measurement functions, respectively. The white Gaussian noises ${\mathit{w}}_{k-1}$ and ${\mathit{v}}_{i,k}$ are assumed to be mutually uncorrelated with covariances ${\mathit{Q}}_{k-1}=\mathbb{E}\{{\mathit{w}}_{k-1}{\mathit{w}}_{k-1}^{\mathrm{T}}\}$ and ${\mathit{R}}_{i,k}=\mathbb{E}\{{\mathit{v}}_{i,k}{\mathit{v}}_{i,k}^{\mathrm{T}}\}$ .

The goal in a distributed filtering paradigm is to obtain, at each time instant $k$ and for each node $i\in \mathcal{N}$ , an estimate ${\hat{\mathit{x}}}_{i,k|k}$ of the unknown state ${\mathit{x}}_k$ , which is achieved with only information sharing with its neighbors ${\mathcal{N}}_i$ . It is also desirable for the estimate ${\hat{\mathit{x}}}_{i,k|k}$ to be comparable to the centralized estimate, which is achieved based on all available measurements over the entire sensor network. It is also essential for each node to achieve consensus on local estimates such that the global decision can be carried out independently at each local node.

A novel Gaussian approximation information filtering paradigm

For the recursive Bayesian filtering in Gaussian domain, the posterior estimate is obtained through two steps: time update and measurement update, which are carried out based on the Chapman–Kolmogorov equation ²⁷ and Bayes rule, respectively. Both steps involve Gaussian weighted integrals to capture the first two moments of the state. However, it is intractable to analytically compute these integrals. Fortunately, some efficient quadrature rules are available to eventually derive different Gaussian approximation filters, such as UKF based on unscented transformation^6,20 and CKF based on spherical-radial cubature rule.^10,28,29 With resort to the statistical linear error propagation methodology,^15,29 the information form of Gaussian approximation filter can be developed, which is useful for multiple sensor estimation. However, the conventional nonlinear information filters, such as the UIF^20,30 and CIF,^21,29 are all designed based on the extended information filter architecture, where the predicted measurement errors are not taken into account in the measurement update. It has been pointed out that these information filters are not equivalent to their Kalman filter counterparts and, hence, are less accurate than the corresponding Kalman filters.^21,22 Here, the statistical linear regression method ²² is introduced to linearize the measurement equation, where the linearization effect is compensated in the linearization noise. This drives the resulting information filter, which is strictly equivalent to its Kalman filter counterpart.

The integral with respect to the standard Gaussian weighting function can be approximated by Arasaratnam and Haykin ¹⁰

$$\underset{{\mathbb{R}}^n}{\int }\mathit{g}\left(\mathit{x}\right)N(\mathit{x};0,\mathit{I})\mathrm{d}\mathit{x}\approx \sum _{s=1}^{N_p}{w}_s\mathit{g}\left({\mathit{\xi }}_s\right)$$

(3)

where $n$ is the dimension of variable $\mathit{x}$ and ${\mathit{\xi }}_s$ and $w_s$ are selected quadrature points and the associated weights, respectively. $N_p$ is the total number of quadrature points. $\mathit{g}(\mathit{x})$ is a nonlinear function. $N(\mathit{x};0,\mathit{I})$ is the probability density function (pdf) of standard multi-variate Gaussian distribution.

For a more general Gaussian distribution $N(\mathit{x};\hat{\mathit{x}},\mathit{P})$ , let $\mathit{x}=\sqrt{\mathit{P}}\mathit{y}+\hat{\mathit{x}}$ , where $\sqrt{\mathit{P}}$ is obtained by Cholesky decomposition such that $\mathit{P}=\sqrt{\mathit{P}}{\sqrt{\mathit{P}}}^{\mathrm{T}}$ . Then, $(\mathit{x}-\hat{\mathit{x}}{)}^{\mathrm{T}}{\mathit{P}}^{-1}(\mathit{x}-\hat{\mathit{x}})={\mathit{y}}^{\mathrm{T}}\mathit{y}$ and $\mathrm{d}\mathit{x}=|\sqrt{\mathit{P}}|\mathrm{d}\mathit{y}=\sqrt{|\mathit{P}|}\mathrm{d}\mathit{y}$ . The desired integral can be written as

$$\begin{array}{c}\multicolumn{1}{c}{\underset{{\mathbb{R}}^n}{\int }\mathit{g}\left(\mathit{x}\right)N(\mathit{x};\hat{\mathit{x}},\mathit{P})\mathrm{d}\mathit{x}={\left(2\pi \right)}^{\frac{-n}{2}}{\left|\mathit{P}\right|}^{\frac{-1}{2}}}\\ \multicolumn{1}{c}{\underset{{\mathbb{R}}^n}{\int }\mathit{g}(\sqrt{\mathit{P}}\mathit{y}+\hat{\mathit{x}})\exp \left(\frac{-{\mathit{y}}^{\mathrm{T}}\mathit{y}}{2}\right)\sqrt{\left|\mathit{P}\right|}\mathrm{d}\mathit{y}}\\ \multicolumn{1}{c}{=\underset{{\mathbb{R}}^n}{\int }\mathit{g}(\sqrt{\mathit{P}}\mathit{y}+\hat{\mathit{x}})N(\mathit{y};0,\mathit{I})\mathrm{d}\mathit{y}}\\ \multicolumn{1}{c}{\approx \sum _{s=1}^{N_p}{w}_s\mathit{g}(\sqrt{\mathit{P}}{\mathit{\xi }}_s+\hat{\mathit{x}})}\end{array}$$

(4)

According to the numerical integration approximations in equations (3) and (4), the Gaussian approximation information filtering paradigm can be summarized as follows. The time update step is identical with the conventional information filter. The predicted estimate is

$$\begin{array}{c}\multicolumn{1}{c}{{\hat{\mathit{x}}}_{k|k-1}=\underset{{\mathbb{R}}^{n_x}}{\int }\mathit{f}\left({\mathit{x}}_{k-1}\right)N({\mathit{x}}_{k-1};{\hat{\mathit{x}}}_{k-1|k-1},{\mathit{P}}_{k-1|k-1})\mathrm{d}{\mathit{x}}_{k-1}}\\ \multicolumn{1}{c}{=\sum _{s=1}^{N_p}{w}_s\mathit{f}\left({\mathit{\chi }}_{s,k-1|k-1}\right)}\end{array}$$

(5)

where ${\mathit{\chi }}_{s,k-1|k-1}=\sqrt{{\mathit{P}}_{k-1|k-1}}{\mathit{\xi }}_s+{\hat{\mathit{x}}}_{k-1|k-1}$ . The associated covariance ${\mathit{P}}_{k|k-1}$ is computed by

$$\begin{array}{c}\multicolumn{1}{c}{{\mathit{P}}_{k|k-1}=\mathbb{E}\left\{({\mathit{x}}_k-{\hat{\mathit{x}}}_{k|k-1}){({\mathit{x}}_k-{\hat{\mathit{x}}}_{k|k-1})}^{\mathrm{T}}|{\mathit{Z}}_{k-1}\right\}}\\ \multicolumn{1}{c}{=\underset{{\mathbb{R}}^{n_x}}{\int }\mathit{f}\left({\mathit{x}}_{k-1}\right){\mathit{f}}^{\mathrm{T}}\left({\mathit{x}}_{k-1}\right)\times N({\mathit{x}}_{k-1};{\hat{\mathit{x}}}_{k-1|k-1},{\mathit{P}}_{k-1|k-1})\mathrm{d}{\mathit{x}}_{k-1}-{\hat{\mathit{x}}}_{k|k-1}{\hat{\mathit{x}}}_{k|k-1}^{\mathrm{T}}+{\mathit{Q}}_{k-1}}\\ \multicolumn{1}{c}{=\sum _{s=1}^{N_p}{w}_s\mathit{f}\left({\mathit{\chi }}_{s,k-1|k-1}\right){\mathit{f}}^{\mathrm{T}}\left({\mathit{\chi }}_{s,k-1|k-1}\right)-{\hat{\mathit{x}}}_{k|k-1}{\hat{\mathit{x}}}_{k|k-1}^{\mathrm{T}}+{\mathit{Q}}_{k-1}}\end{array}$$

(6)

Based on the statistical linear regression method, nonlinear measurement equation (2) can be transformed into the following linear form^22,31

$${\mathit{z}}_{i,k}={\mathcal{H}}_{i,k}{\mathit{x}}_k+{\mathit{b}}_{i,k}+{\overline{\mathit{v}}}_{i,k}$$

(7)

where

$${\mathcal{H}}_{i,k}={\mathit{P}}_{i,\mathrm{xz},k|k-1}^{\mathrm{T}}{\mathit{P}}_{k|k-1}^{-1}$$

(8)

$${\mathit{b}}_{i,k}={\hat{\mathit{z}}}_{i,k|k-1}-{\mathcal{H}}_{i,k}{\hat{\mathit{x}}}_{k|k-1}$$

(9)

and ${\overline{\mathit{v}}}_{i,k}$ is a zero-mean Gaussian noise with covariance ${\overline{\mathit{R}}}_{i,k}={\mathit{R}}_{i,k}+{\mathit{P}}_{i,\mathrm{zz},k|k-1}-{\mathcal{H}}_{i,k}{\mathit{P}}_{k|k-1}{\mathcal{H}}_{i,k}^{\mathrm{T}}$ . The predicted measurement ${\hat{\mathit{z}}}_{i,k|k-1}$ is computed by

$$\begin{array}{c}\multicolumn{1}{c}{{\hat{\mathit{z}}}_{i,k|k-1}=\underset{{\mathbb{R}}^{n_x}}{\int }{\mathit{h}}_i\left({\mathit{x}}_k\right)N({\mathit{x}}_k;{\hat{\mathit{x}}}_{k|k-1},{\mathit{P}}_{k|k-1})\mathrm{d}{\mathit{x}}_k}\\ \multicolumn{1}{c}{=\sum _{s=1}^{N_p}{w}_s{\mathit{h}}_i\left({\mathit{\chi }}_{s,k|k-1}\right)}\end{array}$$

(10)

where ${\mathit{\chi }}_{s,k|k-1}=\sqrt{{\mathit{P}}_{k|k-1}}{\mathit{\xi }}_s+{\hat{\mathit{x}}}_{k|k-1}$ . The ${\mathit{P}}_{i,\mathrm{zz},k|k-1}$ involved in ${\overline{\mathit{R}}}_{i,k}$ is computed by

$$\begin{array}{c}\multicolumn{1}{c}{{\mathit{P}}_{i,\mathrm{zz},k|k-1}=\underset{{\mathbb{R}}^{n_x}}{\int }{\mathit{h}}_i\left({\mathit{x}}_k\right){\mathit{h}}_i^{\mathrm{T}}\left({\mathit{x}}_k\right)N({\mathit{x}}_k;{\hat{\mathit{x}}}_{k|k-1},{\mathit{P}}_{k|k-1})\mathrm{d}{\mathit{x}}_k}\\ \multicolumn{1}{c}{-{\hat{\mathit{z}}}_{k|k-1}{\hat{\mathit{z}}}_{k|k-1}^{\mathrm{T}}}\\ \multicolumn{1}{c}{=\sum _{s=1}^{N_p}{w}_s{\mathit{h}}_i\left({\mathit{\chi }}_{s,k|k-1}\right){\mathit{h}}_i^{\mathrm{T}}\left({\mathit{\chi }}_{s,k|k-1}\right)-{\hat{\mathit{z}}}_{k|k-1}{\hat{\mathit{z}}}_{k|k-1}^{\mathrm{T}}}\end{array}$$

(11)

and the cross-covariance ${\mathit{P}}_{i,\mathrm{xz},k|k-1}$ is computed by

$$\begin{array}{c}\multicolumn{1}{c}{{\mathit{P}}_{i,\mathrm{xz},k|k-1}=\underset{{\mathbb{R}}^{n_x}}{\int }{\mathit{x}}_k{\mathit{h}}_i^{\mathrm{T}}\left({\mathit{x}}_k\right)N({\mathit{x}}_k;{\hat{\mathit{x}}}_{k|k-1},{\mathit{P}}_{k|k-1})\mathrm{d}{\mathit{x}}_k}\\ \multicolumn{1}{c}{-{\hat{\mathit{x}}}_{k|k-1}{\hat{\mathit{z}}}_{i,k|k-1}^{\mathrm{T}}}\\ \multicolumn{1}{c}{=\sum _{s=1}^{N_p}{w}_s{\mathit{\chi }}_{s,k|k-1}{\mathit{h}}_i^{\mathrm{T}}\left({\mathit{\chi }}_{s,k|k-1}\right)-{\hat{\mathit{x}}}_{k|k-1}{\hat{\mathit{z}}}_{i,k|k-1}^{\mathrm{T}}}\end{array}$$

(12)

Therefore, the measurement update for node $i$ in the information form can be expressed as follows

$${\hat{\mathit{y}}}_{k|k}={\hat{\mathit{y}}}_{k|k-1}+\underset{{\mathit{u}}_{i,k}}{\underbrace{{\mathcal{H}}_{i,k}^{\mathrm{T}}{\overline{\mathit{R}}}_{i,k}^{-1}{\stackrel{⌣}{\mathit{z}}}_{i,k}}}$$

(13)

$${\mathit{Y}}_{k|k}={\mathit{Y}}_{k|k-1}+\underset{U_{i,k}}{\underbrace{{\mathcal{H}}_{i,k}^{\mathrm{T}}{\overline{\mathit{R}}}_{i,k}^{-1}{\mathcal{H}}_{i,k}}}$$

(14)

where ${\hat{\mathit{y}}}_{k|k-1}={\mathit{P}}_{k|k-1}^{-1}{\hat{\mathit{x}}}_{k|k-1}$ , ${\mathit{Y}}_{k|k-1}={\mathit{P}}_{k|k-1}^{-1}$ , ${\stackrel{⌣}{\mathit{z}}}_{i,k}={\mathit{z}}_{i,k}-{\widehat{\mathit{z}}}_{i,k|k-1}+{\mathcal{H}}_{i,k}{\widehat{\mathit{x}}}_{i,k|k-1}$ . ${\mathit{u}}_{i,k}$ and ${\mathit{U}}_{i,k}$ are, respectively, the information contribution vector and information contribution matrix.

As for multiple sensor nodes, the centralized estimate can be obtained by summation of ${\mathit{u}}_{i,k}$ and ${\mathit{U}}_{i,k}$ over all sensors

$${\hat{\mathit{y}}}_{k|k}={\hat{\mathit{y}}}_{k|k-1}+\sum _{i\in \mathcal{S}}{\mathit{u}}_{i,k}$$

(15)

$${\mathit{Y}}_{k|k}={\mathit{Y}}_{k|k-1}+\sum _{i\in \mathcal{S}}{\mathit{U}}_{i,k}$$

(16)

The distributed filtering algorithms are designed to achieve comparable performance to the centralized results in equations (15) and (16), which will be shown in the later section. It should be noted that in the conventional nonlinear information filter, ${\mathit{R}}_{i,k}$ is used in equations (13) and (14) instead of ${\overline{\mathit{R}}}_{i,k}$ , where the linearization error ${\mathit{P}}_{i,\mathrm{zz},k|k-1}-{\mathcal{H}}_{i,k}{\mathit{P}}_{k|k-1}{\mathcal{H}}_{i,k}^{\mathrm{T}}$ is neglected. Therefore, the estimation accuracy will be degraded, which will be illustrated in the simulation experiments. In addition, it has been proved that the proposed information filter is strictly equivalent to its centralized Kalman filter counterpart, which is an approximation derived from the recursive Bayesian nonlinear filtering framework.

High-degree spherical simplex-radial quadrature rule

It can be seen from section “A novel Gaussian approximation information filtering paradigm” that the choice of quadrature points and associated weights is of fundamental significance for the estimation accuracy. The basic problem is how to choose ${\mathit{\xi }}_s$ and $w_s$ to achieve higher approximation accuracy in equation (3) than the existing schemes such as unscented and cubature transformations.

Consider the Gaussian weighted integral in equation (3), let $\mathit{x}=r\mathit{s}$ with ${\mathit{s}}^{\mathrm{T}}\mathit{s}=1$ , such that $r=\sqrt{{\mathit{x}}^{\mathrm{T}}\mathit{x}}$ for $r\in [0,+\infty )$ . Then, in the spherical-radial coordinate system, the integral in equation (3) becomes

$$\begin{array}{c}\multicolumn{1}{c}{\underset{{\mathbb{R}}^n}{\int }\mathit{g}\left(\mathit{x}\right)N(\mathit{x};0,\mathit{I})\mathrm{d}\mathit{x}={\left(2\pi \right)}^{\frac{-n}{2}}\underset{{\mathbb{R}}^n}{\int }\mathit{g}\left(\mathit{x}\right)\exp \left(\frac{-{\mathit{x}}^{\mathrm{T}}\mathit{x}}{2}\right)\mathrm{d}\mathit{x}}\\ \multicolumn{1}{c}{={\pi }^{\frac{-n}{2}}\underset{{\mathbb{R}}^n}{\int }\mathit{g}\left(\sqrt{2}\mathit{x}\right)\exp (-{\mathit{x}}^{\mathrm{T}}\mathit{x})\mathrm{d}\mathit{x}}\\ \multicolumn{1}{c}{={\pi }^{\frac{-n}{2}}\underset{0}{\overset{\infty }{\int }}\underset{U_n}{\int }\mathit{g}\left(\sqrt{2}r\mathit{s}\right)r^{n-1}}\\ \multicolumn{1}{c}{\exp (-r^2)\mathrm{d}\sigma \left(\mathit{s}\right)\mathrm{d}r}\end{array}$$

(17)

where $U_n=\{\mathit{s}\in {\mathbb{R}}^n:{\mathit{s}}^{\mathrm{T}}\mathit{s}=1\}$ and $\sigma (·)$ is the spherical surface measure on $U_n$ . Note that the second step in equation (17) is developed with the linear transformation $\mathit{x}=\sqrt{2}\mathit{y}$ . Then, equation (17) can be decomposed into the spherical integral $\mathit{S}(r)={\int }_{U_n}\mathit{g}(\sqrt{2}r\mathit{s})\mathrm{d}\sigma (\mathit{s})$ and the radial integral $\mathit{I}(\mathit{g})={\int }_0^{\infty }\mathit{S}(r)r^{n-1}\exp (-r^2)\mathrm{d}r$ .

The spherical integral $\mathit{S}(r)$ can be effectively solved by Genz’s method with arbitrary accuracy. ³² As higher degree spherical rules are applied, the approximation is closer to the true value, but it also requires more sampling points and hence increases computational burden for the system. To make a balance between the accuracy and computational burden, here the Mysovskikh’s method ³³ is chosen for higher accuracy but less sampling points than other methods. The third-degree spherical simplex-radial rule is built as

$${\mathit{S}}_3\left(r\right)=\frac{A_n}{2(n+1)}\sum _{j=1}^{n+1}(\mathit{g}\left(\sqrt{2}r{\overline{\mathit{s}}}_j\right)+\mathit{g}(-\sqrt{2}r{\overline{\mathit{s}}}_j))$$

(18)

where $A_n=2{\mathrm{\Gamma }}^n(1/2)/\mathrm{\Gamma }(n/2)=2\sqrt{{\pi }^n}/\mathrm{\Gamma }(n/2)$ is the surface area of the unit sphere. $\mathrm{\Gamma }(z)$ is the gamma function defined by $\mathrm{\Gamma }(z)={\int }_0^{\infty }\exp (-\lambda ){\lambda }^{z-1}d\lambda$ . The points ${\overline{s}}_j$ are determined by

$${\overline{\mathit{s}}}_j={[{\overline{s}}_{j,1},{\overline{s}}_{j,2},\dots ,{\overline{s}}_{j,n}]}^{\mathrm{T}},j=1,2,\dots ,n+1$$

with

$${\overline{s}}_{j,i}=\{\begin{array}{cc}-\sqrt{\frac{n+1}{n(n-i+2)(n-i+1)}}\hfill & 1\le i<j\hfill \\ \multicolumn{1}{c}{\sqrt{\frac{(n+1)(n-j+1)}{n(n-j+2)}}}& i=j,i\le n\hfill \\ \multicolumn{1}{c}{0}& j<i\le n\hfill \end{array}$$

In a similar way, the fifth-degree spherical simplex-radial rule is constructed by

$$\begin{array}{c}\multicolumn{1}{c}{{\mathit{S}}_5\left(r\right)=w_{s,1}\sum _{j=1}^{n+1}(\mathit{g}\left(\sqrt{2}r{\overline{\mathit{s}}}_j\right)+\mathit{g}(-\sqrt{2}r{\overline{\mathit{s}}}_j))}\\ \multicolumn{1}{c}{+w_{s,2}\sum _{j=1}^{n(n+1)/2}(\mathit{g}\left(\sqrt{2}r{\tilde{\mathit{s}}}_j\right)+\mathit{g}(-\sqrt{2}r{\tilde{\mathit{s}}}_j))}\end{array}$$

(19)

where the points set of ${\tilde{\mathit{s}}}_j$ is determined by

$$\left\{{\tilde{\mathit{s}}}_j\right\}=\{\sqrt{\frac{n}{2(n-1)}}({\overline{s}}_k+{\overline{s}}_l):k<l;k,l=1,2,\dots ,n+1\}$$

with associated weights

$${\overline{w}}_{s,1}=\frac{(7-n)n{A}_n}{2{(n+1)}^2(n+2)},{\overline{w}}_{s,2}=\frac{2{(n-1)}^2{A}_n}{n{(n+1)}^2(n+2)}$$

Note that the number of points in set $\{{\tilde{\mathit{s}}}_j\}$ is $C_{n+1}^2=n(n+1)/2$ , and thus, the total number of points is $N_s=2\times (n+1)+2\times C_{n+1}^2=n^2+3n+2$ , while the number of points for the fifth-degree spherical rule with Genz’s method is $2n^2$ .^28,34 For high dimensional weighted integrals with $n\ge 4$ , the Mysovskikh’s method requires less points and, hence, shows better performance in terms of computation cost.

The radial integral $\mathit{I}(\mathit{g})$ is computed by the moment-matching method, which can achieve arbitrary degree radial rule. The idea is to find the points $r_i$ and the associated weights $w_{r,i}$ which satisfy the following moment equation

$$\underset{0}{\overset{\infty }{\int }}\mathit{S}\left(r\right)r^{n-1}\exp (-r^2)\mathrm{d}r=\sum _{i=1}^{N_r}{w}_{r,i}\mathit{S}\left(r_i\right)$$

(20)

where $\mathit{S}(r)=r^l$ is a monomial in $r$ with an even integer $l$ . Then, the left-hand side of equation (20) reduces to $\mathrm{\Gamma }((n+1)/2)/2$ . Due to the fact that the spherical-simplex rule and the spherical simplex-radial rule are both symmetric, only the even-degree monomials need to be matched. To acquire a (2m + 1)th-degree radial rule, it is necessary for equation (20) to be exact with $l=0,2,\dots ,2m$ , which consists of $(m+1)$ equations. ²⁸

For the third-degree radial rule $(N_r=1)$ , only one quadrature point $\{r_1,w_{r,1}\}$ is enough to match equation (20) with $l=0$ and $l=2$ . To determine $\{r_1,w_{r,1}\}$ , the following equations should be satisfied

$$\{\begin{array}{c}{w}_{r,1}{r}_1^0=\frac{1}{2}\mathrm{\Gamma }\left(\frac{n}{2}\right)\hfill \\ w_{r,1}{r}_1^2=\frac{1}{2}\mathrm{\Gamma }(\frac{n}{2}+1)=\frac{n}{4}\mathrm{\Gamma }\left(\frac{n}{2}\right)\hfill \end{array}$$

(21)

where the last equality holds for the fact that $\mathrm{\Gamma }(z+1)=z\mathrm{\Gamma }(z)$ . The point and associated weight are given by solving equation (21) as

$$\{\begin{array}{c}{r}_1=\sqrt{\frac{1}{2}n}\hfill \\ w_{r,1}=\frac{\mathrm{\Gamma }\left(\frac{n}{2}\right)}{2}\hfill \end{array}$$

(22)

Similarly, the fifth-degree radial rule with two quadrature points is determined as

$$\underset{0}{\overset{\infty }{\int }}\mathit{S}\left(r\right)r^{n-1}\exp (-r^2)\mathrm{d}r=w_{r,1}\mathit{S}\left(0\right)+w_{r,2}\mathit{S}\left(\sqrt{\frac{n}{2}+1}\right)$$

(23)

where $w_{r,1}=\mathrm{\Gamma }(n/2)/(n+2)$ and $w_{r,2}=n\mathrm{\Gamma }(n/2)/(2n+4)$ . The interested reader is referred to Jia et al. ²⁸ for more details.

Exploiting equations (18) and (22) in equation (17), the third-degree spherical simplex-radial rule is obtained as

$$\begin{array}{c}\underset{{\mathbb{R}}^n}{\int }\mathit{g}\left(\mathit{x}\right)N(\mathit{x};0,\mathit{I})\mathrm{d}\mathit{x}\hfill \\ =\frac{w_{r,1}}{{\pi }^{n/2}}·\frac{A_n}{2(n+1)}\sum _{j=1}^{n+1}(\mathit{g}\left(\sqrt{2}{r}_1{\overline{\mathit{s}}}_{\mathit{j}}\right)+\mathit{g}(-\sqrt{2}{r}_1{\overline{\mathit{s}}}_j))\hfill \\ =\frac{1}{2(n+1)}\sum _{j=1}^{n+1}(\mathit{g}\left(\sqrt{n}{\overline{\mathit{s}}}_j\right)+\mathit{g}(-\sqrt{n}{\overline{\mathit{s}}}_j))\hfill \end{array}$$

(24)

Similarly, substituting equations (19) and (23) into equation (17), the fifth-degree spherical simplex-radial rule is formulated as

$$\begin{array}{c}\underset{{\mathbb{R}}^n}{\int }\mathit{g}\left(\mathit{x}\right)N(\mathit{x};0,\mathit{I})\mathrm{d}\mathit{x}={\pi }^{\frac{-n}{2}}(w_{r,1}{\mathit{S}}_5\left(0\right)+w_{r,2}{\mathit{S}}_5\left(\sqrt{\left(\frac{n}{2}\right)+1}\right))\hfill \\ =\frac{2}{n+2}\mathit{g}\left(0\right)+\frac{n^2(7-n)}{2{(n+1)}^2{(n+2)}^2}\hfill \\ \sum _{j=1}^{n+1}(\mathit{g}(\sqrt{n+2}·{\overline{\mathit{s}}}_j)+\mathit{g}(-\sqrt{n+2}·{\overline{s}}_j))\hfill \\ +\frac{2{(n-1)}^2}{{(n+1)}^2{(n+2)}^2}\hfill \\ \sum _{j=1}^{n(n+1)/2}(\mathit{g}(\sqrt{n+2}·{\tilde{\mathit{s}}}_j)+\mathit{g}(-\sqrt{n+2}·{\tilde{\mathit{s}}}_j))\hfill \end{array}$$

(25)

Therefore, spherical simplex-radial quadrature points ${\mathit{\xi }}_s$ and the associated weights $w_s$ in equation (3) are determined.

Distributed SRCQIF-HC

As is discussed in Battistelli and Chisci, ¹² the CI method can guarantee stability in the mean square error boundedness sense for any number of consensus iterations (even a single iteration). However, it adopts a conservative fusion rule based on the assumption that the correlation between the estimates from different nodes is completely unknown. As a result, the novel information is unavoidably underestimated such that the estimation accuracy may be deteriorated. In contrast, the CM method only fuses the novel information and hence avoids any conservative assumption on the correlation between estimates. The multiplication by some scalar weights in CM can result in improved estimation performance. However, it ensures stability only with an adequate number of consensus iterations.^8,16 Undoubtedly, this will be challenging when communication resources are constrained, which is often the case in actual applications. To integrate the benefits and neutralize the drawbacks of both approaches, the HCMCI method is developed to perform consensus on both prior and novel information. ¹⁶ In addition, it can be seen from section “A novel Gaussian approximation information filtering paradigm” that the generation of quadrature points involves square-root decomposition of covariance matrix, which is a numerically sensitive operation requiring the covariance to be positive definite. However, due to limited word length and rounding errors in an embedded system, the positive definiteness may be destroyed and, thus, render the system unstable. To further enhance the numerical stability and accomplish the nonlinear filtering task in a distributed manner, a distributed SRCQIF-HC is proposed.

Square-root cubature quadrature information filter based on hybrid consensus

Assume, at time instant $k$ , the prior information vector ${\hat{\mathit{y}}}_{i,k|k-1}$ and square-root factor ${\mathit{S}}_{i,y,k|k-1}$ with ${\mathit{Y}}_{i,k|k-1}={\mathit{S}}_{i,y,k|k-1}{\mathit{S}}_{i,y,k|k-1}^{\mathrm{T}}$ are known. Since ${\mathit{P}}_{i,k|k-1}={\mathit{Y}}_{i,k|k-1}^{-1}={\mathit{S}}_{i,k|k-1}{\mathit{S}}_{i,k|k-1}^{\mathrm{T}}$ , one has ${\mathit{S}}_{i,k|k-1}={\mathit{S}}_{i,y,k|k-1}^{-\mathrm{T}}$ , and hence, the quadrature points become ${\mathit{\chi }}_{i,s,k|k-1}={\mathit{S}}_{i,y,k|k-1}^{-\mathrm{T}}{\mathit{\xi }}_s+{\hat{\mathit{x}}}_{i,k|k-1}$ with ${\hat{\mathit{x}}}_{i,k|k-1}={\mathit{S}}_{i,y,k|k-1}^{-\mathrm{T}}{\mathit{S}}_{i,y,k|k-1}^{-1}{\hat{\mathit{y}}}_{k|k-1}$ . Then, the pseudo measurement matrix ${\mathcal{H}}_{i,k}$ can be computed from equations (8) and (12) as follows

$$\begin{array}{cc}\multicolumn{1}{c}{{\mathcal{H}}_{i,k}}& ={\mathit{P}}_{i,\mathrm{xz},k|k-1}^{\mathrm{T}}{\mathit{Y}}_{i,k|k-1}\hfill \\ \multicolumn{1}{c}{}& ={\mathcal{Z}}_{i,k|k-1}{\mathcal{X}}_{i,k|k-1}^{\mathrm{T}}{\mathit{S}}_{i,y,k|k-1}{\mathit{S}}_{i,y,k|k-1}^{\mathrm{T}}\hfill \end{array}$$

(26)

where

$$\begin{array}{l}{\mathcal{X}}_{i,k|k-1}=\left[\begin{array}{cccc}\sqrt{w_1}\left({\mathit{\chi }}_{i,1,k|k-1}-{\widehat{\mathit{x}}}_{i,k|k-1}\right)& \sqrt{w_2}\left({\mathit{\chi }}_{i,2,k|k-1}-{\widehat{\mathit{x}}}_{i,k|k-1}\right)& \cdots & \sqrt{w_{N_p}}\left({\mathit{\chi }}_{i,N_p,k|k-1}-{\widehat{\mathit{x}}}_{i,k|k-1}\right)\end{array}\right]\\ {\mathcal{Z}}_{i,k|k-1}=\left[\begin{array}{cccc}\sqrt{w_1}\left({\mathit{{\rm Z}}}_{i,1,k|k-1}-{\widehat{\mathit{z}}}_{i,k|k-1}\right)& \sqrt{w_2}\left({\mathit{{\rm Z}}}_{i,2,k|k-1}-{\widehat{\mathit{z}}}_{i,k|k-1}\right)& \cdots & \sqrt{w_{N_p}}\left({\mathit{{\rm Z}}}_{i,N_p,k|k-1}-{\widehat{\mathit{z}}}_{i,k|k-1}\right)\end{array}\right]\\ {\mathit{{\rm Z}}}_{i,s,k|k-1}={\mathit{h}}_i\left({\mathit{\chi }}_{i,s,k|k-1}\right)\end{array}$$

and ${\hat{z}}_{i,k|k-1}$ can be computed by equation (10).

According to equation (8), one has

$${\mathcal{H}}_{i,k}{\mathit{P}}_{i,k|k-1}{\mathcal{H}}_{i,k}^{\mathrm{T}}={\mathcal{H}}_{i,k}{\mathit{P}}_{i,\mathrm{xz},k|k-1}={\mathit{P}}_{i,\mathrm{xz},k|k-1}^{\mathrm{T}}{\mathcal{H}}_{i,k}^{\mathrm{T}}$$

(27)

Exploiting equations (26) and (27) in the computation of ${\overline{\mathit{R}}}_{i,k}$ , one has

$$\begin{array}{cc}\multicolumn{1}{c}{{\overline{\mathit{R}}}_{i,k}}& ={\mathit{R}}_{i,k}+{\mathit{P}}_{i,\mathrm{zz},k|k-1}-{\mathcal{H}}_{i,k}{\mathit{P}}_{i,k|k-1}{\mathcal{H}}_{i,k}^{\mathrm{T}}\hfill \\ \multicolumn{1}{c}{}& ={\mathit{R}}_{i,k}+{\mathit{P}}_{i,\mathrm{zz},k|k-1}-{\mathcal{H}}_{i,k}{\mathit{P}}_{i,\mathrm{xz},k|k-1}-{\mathit{P}}_{i,\mathrm{xz},k|k-1}^{\mathrm{T}}{\mathcal{H}}_{i,k}^{\mathrm{T}}+{\mathcal{H}}_{i,k}{\mathit{P}}_{i,k|k-1}{\mathcal{H}}_{i,k}^{\mathrm{T}}\hfill \\ \multicolumn{1}{c}{}& ={\mathit{S}}_{i,R,k}{\mathit{S}}_{i,R,k}^{\mathrm{T}}+{\mathcal{Z}}_{i,k|k-1}{\mathcal{Z}}_{i,k|k-1}^{\mathrm{T}}-{\mathcal{H}}_{i,k}{\mathcal{X}}_{i,k|k-1}{\mathcal{Z}}_{i,k|k-1}^{\mathrm{T}}-{\mathcal{Z}}_{i,k|k-1}{\mathcal{X}}_{i,k|k-1}^{\mathrm{T}}{\mathcal{H}}_{i,k}^{\mathrm{T}}+{\mathcal{H}}_{i,k}{\mathcal{X}}_{i,k|k-1}{\mathcal{X}}_{i,k|k-1}^{\mathrm{T}}{\mathcal{H}}_{i,k}^{\mathrm{T}}\hfill \\ \multicolumn{1}{c}{}& ={\mathit{S}}_{i,R,k}{\mathit{S}}_{i,R,k}^{\mathrm{T}}+({\mathcal{Z}}_{i,k|k-1}-{\mathcal{H}}_{i,k}{\mathcal{X}}_{i,k|k-1}){({\mathcal{Z}}_{i,k|k-1}-{\mathcal{H}}_{i,k}{\mathcal{X}}_{i,k|k-1})}^{\mathrm{T}}\hfill \\ \multicolumn{1}{c}{}& =[{\mathcal{Z}}_{i,k|k-1}-{\mathcal{H}}_{i,k}{\mathcal{X}}_{i,k|k-1}{\mathit{S}}_{i,R,k}]{[{\mathcal{Z}}_{i,k|k-1}-{\mathcal{H}}_{i,k}{\mathcal{X}}_{i,k|k-1}{\mathit{S}}_{i,R,k}]}^{\mathrm{T}}\hfill \end{array}$$

(28)

Therefore, the square-root of ${\overline{\mathit{R}}}_{i,k}$ can be written as

$${\mathit{S}}_{i,\overline{R},k}=\mathit{Tria}\left(\left[\begin{array}{cc}{\mathcal{Z}}_{i,k|k-1}-{\mathcal{H}}_{i,k}{\mathcal{X}}_{i,k|k-1}\hfill & {\mathit{S}}_{i,R,k}\hfill \end{array}\right]\right)$$

(29)

Here, $\mathit{S}=\mathit{Tria}(\mathit{A})$ denotes the triangularization decomposition with $\mathit{S}=[\mathrm{QR}({\mathit{A}}^{\mathrm{T}}){]}^{\mathrm{T}}$ , where $\mathrm{QR}$ represents the QR decomposition.^10,35

Substituting equation (29) into equations (13) and (14), one has

$${\mathit{u}}_{i,k}={\mathcal{H}}_{i,k}^{\mathrm{T}}{\mathit{S}}_{i,\overline{R},k}^{-\mathrm{T}}{\mathit{S}}_{i,\overline{R},k}^{-1}({\mathit{z}}_{i,k}-{\hat{\mathit{z}}}_{i,k|k-1}+{\mathcal{H}}_{i,k}{\hat{\mathit{x}}}_{i,k|k-1})$$

(30)

$${\mathit{S}}_{i,u,k}={\mathcal{H}}_{i,k}^{\mathrm{T}}{\mathit{S}}_{i,\overline{R},k}^{-\mathrm{T}}$$

(31)

Performing consensus on the prior information, one has

$${\hat{\mathit{y}}}_{i,k|k-1}^l={\hat{\mathit{y}}}_{i,k|k-1}^{l-1}+\sum _{j\in {\mathcal{N}}_i}{\pi }_{i,j}({\hat{\mathit{y}}}_{j,k|k-1}^{l-1}-{\hat{\mathit{y}}}_{i,k|k-1}^{l-1})$$

(32)

$${\mathit{Y}}_{i,k|k-1}^l={\mathit{Y}}_{i,k|k-1}^{l-1}+\sum _{j\in {\mathcal{N}}_i}{\pi }_{i,j}({\mathit{Y}}_{j,k|k-1}^{l-1}-{\mathit{Y}}_{i,k|k-1}^{l-1})$$

(33)

where $l$ denotes the consensus iterations with initialization of ${\hat{\mathit{y}}}_{i,k|k-1}^0={\hat{\mathit{y}}}_{i,k|k-1}$ , ${\mathit{Y}}_{i,k|k-1}^0={\mathit{Y}}_{i,k|k-1}$ . ${\pi }_{i,j}$ is the consensus weight satisfying $\sum _{j\in {\mathcal{J}}_i}{\pi }_{i,j}=1,{\pi }_{i,j}\ge 0$ . The square-root form of equation (33) can be written as

$$\begin{array}{cc}\multicolumn{1}{c}{{\mathit{Y}}_{i,k|k-1}^l}& ={\pi }_{i,i}{\mathit{S}}_{i,y,k|k-1}^{l-1}{\left({\mathit{S}}_{i,y,k|k-1}^{l-1}\right)}^{\mathrm{T}}\hfill \\ \multicolumn{1}{c}{}& +\sum _{j\in {\mathcal{N}}_i}{\pi }_{i,j}{\mathit{S}}_{j,y,k|k-1}^{l-1}{\left({\mathit{S}}_{j,y,k|k-1}^{l-1}\right)}^{\mathrm{T}}\hfill \\ \multicolumn{1}{c}{}& ={\mathit{S}}_{i,y,k|k-1}^l{\left({\mathit{S}}_{i,y,k|k-1}^l\right)}^{\mathrm{T}}\hfill \end{array}$$

(34)

where

$${\mathit{S}}_{i,y,k|k-1}^l=\mathrm{Tria}\left(\left[\begin{array}{cccc}\sqrt{{\pi }_{i,i}}{\mathit{S}}_{i,y,k|k-1}^{l-1}\hfill & \sqrt{{\pi }_{i,j_1}}{\mathit{S}}_{j_1,y,k|k-1}^{l-1}\hfill & \cdots \hfill & \sqrt{{\pi }_{i,j_{\left|{\mathcal{N}}_i\right|}}}{\mathit{S}}_{j_{{}_{\left|{\mathcal{N}}_i\right|}},y,k|k-1}^{l-1}\hfill \end{array}\right]\right)$$

(35)

Here, ${\mathit{S}}_{i,y,k|k-1}^0={\mathit{S}}_{i,y,k|k-1}$ , $j_h\in {\mathcal{N}}_i,h=1,2,\dots ,|{\mathcal{N}}_i|$ is the index of neighbors of node $i$ . The notation $|•|$ means the cardinality of a set.

Similarly, performing consensus on the novel information with initialization of ${\mathit{u}}_{i,k}^0={\mathit{u}}_{i,k}$ , ${\mathit{U}}_{i,k}^0={\mathit{U}}_{i,k}$ , one has

$${\mathit{u}}_{i,k}^l={\mathit{u}}_{i,k}^{l-1}+\sum _{j\in {\mathcal{N}}_i}{\pi }_{i,j}({\mathit{u}}_{j,k}^{l-1}-{\mathit{u}}_{i,k}^{l-1})$$

(36)

$${\mathit{U}}_{i,k}^l={\mathit{U}}_{i,k}^{l-1}+\sum _{j\in {\mathcal{N}}_i}{\pi }_{i,j}({\mathit{U}}_{j,k}^{l-1}-{\mathit{U}}_{i,k}^{l-1})$$

(37)

Substituting equation (31) into (37) yields

$$\begin{array}{cc}\multicolumn{1}{c}{{\mathit{U}}_{i,k}^l}& ={\pi }_{i,i}{\mathit{S}}_{i,u,k}^{l-1}{\left({\mathit{S}}_{i,u,k}^{l-1}\right)}^{\mathrm{T}}+\sum _{j\in {\mathcal{N}}_i}{\pi }_{i,j}{\mathit{S}}_{j,u,k}^{l-1}{\left({\mathit{S}}_{j,u,k}^{l-1}\right)}^{\mathrm{T}}\hfill \\ \multicolumn{1}{c}{}& ={\mathit{S}}_{i,u,k}^l{\left({\mathit{S}}_{i,u,k}^l\right)}^{\mathrm{T}}\hfill \end{array}$$

(38)

where

$$\begin{array}{c}{\mathit{S}}_{i,u,k}^l=\hfill \\ \mathit{Tria}\left(\left[\begin{array}{cccc}\sqrt{{\pi }_{i,i}}{\mathit{S}}_{i,u,k}^{l-1}\hfill & \sqrt{{\pi }_{i,j_1}}{\mathit{S}}_{j_1,u,k}^{l-1}\hfill & \cdots \hfill & \sqrt{{\pi }_{i,j_{\left|{\mathcal{N}}_i\right|}}}{\mathit{S}}_{j_{\left|{\mathcal{N}}_i\right|},u,k}^{l-1}\hfill \end{array}\right]\right)\hfill \end{array}$$

(39)

Here, ${\mathit{S}}_{i,u,k}^0={\mathit{S}}_{i,u,k}$ . The posterior estimate is updated by

$${\hat{\mathit{y}}}_{i,k|k}={\hat{\mathit{y}}}_{i,k|k-1}^L+{\lambda }_{i,k}{\mathit{u}}_{i,k}^L$$

(40)

$${\mathit{Y}}_{i,k|k}={\mathit{Y}}_{i,k|k-1}^L+{\lambda }_{i,k}{\mathit{U}}_{i,k}^L$$

(41)

where ${\lambda }_{i,k}$ is the node-independent scalar weights.

Due to the fact that ${\mathit{Y}}_{i,k|k}={\mathit{S}}_{i,y,k|k-1}^L{\left({\mathit{S}}_{i,y,k|k-1}^L\right)}^{\mathrm{T}}+{\mathit{\lambda }}_{i,k}{\mathit{S}}_{i,u,k}^L{\left({\mathit{S}}_{i,u,k}^L\right)}^{\mathrm{T}}$ , the result in equation (41) can be rewritten as

$${\mathit{S}}_{i,y,k|k}=\mathit{Tria}\left(\left[\begin{array}{cc}{\mathit{S}}_{i,y,k|k-1}^L\hfill & \sqrt{{\lambda }_{i,k}}{\mathit{S}}_{i,u,k}^L\hfill \end{array}\right]\right)$$

(42)

where $L$ is the total number of consensus iterations.

The prior estimate ${\hat{\mathit{x}}}_{i,k+1|k}$ at time instant $k+1$ can be computed by equation (5), and the corresponding square-root of ${\mathit{P}}_{i,k+1|k}$ can be computed by

$${\mathit{S}}_{i,k+1|k}=\mathit{Tria}\left(\left[\begin{array}{cc}{\mathcal{X}}_{i,k+1|k}^{*}\hfill & {\mathit{S}}_{Q,k}\hfill \end{array}\right]\right)$$

(43)

where ${\mathit{S}}_{Q,k}$ is the square-root of ${\mathit{Q}}_k$ satisfying ${\mathit{Q}}_k={\mathit{S}}_{Q,k}{\mathit{S}}_{Q,k}^{\mathrm{T}}$ and

$${\mathcal{X}}_{i,k+1|k}^{*}=\left[\begin{array}{cccc}\sqrt{w_1}({\mathit{X}}_{i,1,k+1|k}-{\hat{\mathit{x}}}_{i,k+1|k})\hfill & \sqrt{w_2}({\mathit{X}}_{i,2,k+1|k}-{\hat{\mathit{x}}}_{i,k+1|k})\hfill & \cdots \hfill & \sqrt{w_{N_p}}({\mathit{X}}_{i,N_p,k+1|k}-{\hat{\mathit{x}}}_{i,k+1|k})\hfill \end{array}\right]$$

(44)

Here, ${\mathit{X}}_{i,s,k+1|k}=\mathit{f}({\mathit{S}}_{i,y,k|k}^{-\mathrm{T}}{\mathit{\xi }}_s+{\mathit{S}}_{i,y,k|k}^{-\mathrm{T}}{\mathit{S}}_{i,y,k|k}^{-1}{\hat{\mathit{y}}}_{i,k|k})$ . Due to the fact that ${\mathit{Y}}_{i,k+1|k}={\mathit{P}}_{i,k+1|k}^{-1}$ , one has

$$\{\begin{array}{c}{\mathit{S}}_{i,y,k+1|k}={\mathit{S}}_{i,k+1|k}^{-\mathrm{T}}\hfill \\ \multicolumn{1}{c}{{\hat{\mathit{y}}}_{i,k+1|k}={\mathit{S}}_{i,y,k+1|k}{\mathit{S}}_{i,y,k+1|k}^{\mathrm{T}}{\hat{\mathit{x}}}_{i,k+1|k}}\end{array}$$

(45)

Up to now, the filtering recursion from time instant $k$ to $k+1$ is completed. The detailed steps of the proposed SRCQIF-HC are summarized in Algorithm 1.

Algorithm 1. Distributed SRCQIF-HC in node $i$ at time instant $k+1$ .

Input: The prior information message $\{{\hat{\mathit{y}}}_{i,k|k-1},{\mathit{S}}_{i,y,k|k-1}\}$ at time instant k, the consensus weights ${\pi }_{\mathrm{ij}}$ , and the total number of consensus iterations L. Step 1. Compute the square-root ${\mathit{S}}_{i,\overline{R},k}$ of the total linearized measurement noise covariance matrix by equation (29). Step 2. Compute the information vector ${\mathit{u}}_{i,k}$ and the square-root information matrix ${\mathit{S}}_{i,u,k}$ by equations (30) and (31). If node $i$ is a communication node $i\in \mathcal{C}$ , ${\mathit{u}}_{i,k}=0$ , ${\mathit{S}}_{i,u,k}=0$ . Step 3. Perform consensus iterations   1) Initialization: ${\hat{\mathit{y}}}_{i,k|k-1}^0={\hat{\mathit{y}}}_{i,k|k-1}$ , ${\mathit{S}}_{i,y,k|k-1}^0={\mathit{S}}_{i,y,k|k-1}$ , ${\mathit{u}}_{i,k}^0={\mathit{u}}_{i,k}$ , ${\mathit{S}}_{i,u,k}^0={\mathit{S}}_{i,u,k}$ .   2) Consensus for $l=1,2,\dots ,L$ do   a) Broadcast information message $\{{\hat{\mathit{y}}}_{i,k|k-1}^{l-1},{\mathit{S}}_{i,y,k|k-1}^{l-1},{\mathit{u}}_{i,k}^{l-1},{\mathit{S}}_{i,u,k}^{l-1}\}$ to its immediate neighbors $j\in {\mathcal{N}}_i$ ;   b) Receive information message $\{{\hat{\mathit{y}}}_{j,k|k-1}^{l-1},{\mathit{S}}_{j,y,k|k-1}^{l-1},{\mathit{u}}_{j,k}^{l-1},{\mathit{S}}_{j,u,k}^{l-1}\}$ from its immediate neighbors $j\in {\mathcal{N}}_i$ ;   c) Fuse the prior information quantities ${\hat{\mathit{y}}}_{i,k|k-1}^{l-1}$ and ${\mathit{S}}_{i,y,k|k-1}^{l-1}$ by equations (32) and (35), and in parallel fuse the novel information quantities ${\mathit{u}}_{i,k}^{l-1}$ and ${\mathit{S}}_{i,u,k}^{l-1}$ by equations (36) and (39). end for Step 4. Update the posterior information vector ${\hat{\mathit{y}}}_{i,k|k}$ and the corresponding square-root information matrix ${\mathit{S}}_{i,y,k|k}$ by equations (40) and (42). Step 5. Compute the prior information vector ${\hat{\mathit{y}}}_{i,k+1|k}$ at time instant $k+1$ and the corresponding square-root information matrix ${\mathit{S}}_{i,y,k+1|k}$ by equation (45). Output: The prior information message $\{{\hat{\mathit{y}}}_{i,k+1|k},{\mathit{S}}_{i,y,k+1|k}\}$ at time instant $k+1$ .

Consensus and convergence

Consensus and convergence are of great importance for the consensus-based algorithms. Here, a brief discussion about the two issues is given below. Before analysis of the proposed SRCQIF-HC, the following lemma is necessary.

Lemma 1

If the consensus matrix $\mathit{\Pi }$ , with (i, j)th element being ${\pi }_{i,j},i,j\in \mathcal{N}$ , is doubly stochastic and primitive, then^2,7,36

$$\underset{x\to \infty }{lim}{\mathit{\Pi }}^l=\frac{1}{N}11^{\mathrm{T}}$$

(46)

where ${\mathit{\Pi }}^l$ is the lth power of $\mathit{\Pi }$ and 1 is a column vector with all 1 entries.

Theorem 1

If $\mathit{\Pi }$ is doubly stochastic and primitive, the information package $\{{\hat{\mathit{y}}}_{i,k|k},{\mathit{S}}_{i,y,k|k}\}$ obtained by the proposed SRCQIF-HC for each node $i\in \mathcal{N}$ will achieve average consensus as $L\to \infty$ .

Proof

It is possible to write the prior information vector ${\hat{\mathit{y}}}_{i,k|k-1}^l$ for all nodes in a collective form as ${\hat{\mathit{y}}}_{k|k-1}^l\stackrel{\mathrm{\Delta }}{=}\mathrm{col}({\hat{\mathit{y}}}_{i,k|k-1}^l,i\in \mathcal{N})$ . Then, equation (32) can be formulated as

$$\begin{array}{cc}\multicolumn{1}{c}{{\hat{\mathit{y}}}_{k|k-1}^l}& =(\mathit{\Pi }\otimes I_{n_x}){\hat{\mathit{y}}}_{k|k-1}^{l-1}\hfill \\ \multicolumn{1}{c}{}& =(\mathit{\Pi }\otimes I_{n_x})\cdots (\mathit{\Pi }\otimes I_{n_x})(\mathit{\Pi }\otimes I_{n_x}){\hat{\mathit{y}}}_{k|k-1}^0\hfill \\ \multicolumn{1}{c}{}& =({\mathit{\Pi }}^l\otimes I_{n_x}){\hat{\mathit{y}}}_{k|k-1}^0\hfill \end{array}$$

(47)

Since $\mathit{\Pi }$ is doubly stochastic and primitive, in view of Lemma 1, $\underset{l\to 0}{lim}{\mathit{\Pi }}^l={11}^{\mathrm{T}}/N$ . Then one has

$$\underset{l\to \infty }{lim}{\hat{\mathit{y}}}_{k|k-1}^l=\frac{1}{N}({11}^{\mathrm{T}}\otimes I_{n_x}){\hat{\mathit{y}}}_{k|k-1}^0$$

(48)

That is, for any node $i\in \mathcal{N}$

$$\begin{array}{c}\multicolumn{1}{c}{\underset{l\to \infty }{lim}{\hat{\mathit{y}}}_{i,k|k-1}^l=\frac{1}{N}\sum _{j\in \mathcal{N}}{\hat{\mathit{y}}}_{j,k|k-1}^0}\\ \multicolumn{1}{c}{={\hat{\mathit{y}}}_{k|k-1}^{*}}\end{array}$$

(49)

Since ${\mathit{S}}_{i,y,k|k-1}^l$ is the square-root of ${\mathit{Y}}_{i,y,k|k-1}^l$ , equation (35) is mathematically equivalent to equation (33). In the similar way, for any node $i\in \mathcal{N}$ , one has

$$\begin{array}{cc}\multicolumn{1}{c}{\underset{l\to \infty }{lim}{\mathit{Y}}_{i,y,k|k-1}^l}& =\frac{1}{N}\sum _{j\in \mathcal{N}}{\mathit{Y}}_{j,k|k-1}^0\hfill \\ \multicolumn{1}{c}{}& ={\mathit{Y}}_{k|k-1}^{*}\hfill \end{array}$$

(50)

Therefore, ${\mathit{S}}_{i,y,k|k-1}^l$ converges to ${\mathit{S}}_{y,k|k-1}^{*}$ as $l\to \infty$ . With the similar technique, when $L$ goes to infinity, one can conclude that ${\mathit{u}}_{i,k}^L$ and ${\mathit{S}}_{i,u,k}^L$ will achieve convergence to ${\mathit{u}}_k^{*}$ and ${\mathit{S}}_{u,k}^{*}$ , respectively. It should be noted that the scalar weights ${\lambda }_{i,k}$ are node independent, and for each node, the same ${\lambda }_k^{*}$ is applied. Therefore, the posterior information vector ${\hat{\mathit{y}}}_{i,k|k}$ and the square-root information matrix ${\mathit{S}}_{i,y,k|k}$ will be of consensus as $L\to \infty$

$$\{\begin{array}{c}{\hat{\mathit{y}}}_{i,k|k}={\hat{\mathit{y}}}_{k|k-1}^{*}+{\lambda }_k^{*}{\mathit{u}}_k^{*}\hfill \\ \multicolumn{1}{c}{{\mathit{S}}_{i,y,k|k}=\mathit{Tria}\left(\left[\begin{array}{cc}{\mathit{S}}_{y,k|k-1}^{*}\hfill & \sqrt{{\lambda }_k^{*}}{\mathit{S}}_{u,k}^{*}\hfill \end{array}\right]\right)}\end{array}$$

(51)

The proof is now completed.

It should be noted that after consensus, the estimation accuracy depends on the choice of the scalar weights ${\lambda }_{i,k}$ . The following Theorem 2 ensures the convergence of the proposed SRCQIF-HC to its centralized counterpart with ${\lambda }_{i,k}=N$ .

Theorem 2

If $\mathit{\Pi }$ is doubly stochastic and primitive, and ${\lambda }_{i,k}=N$ , the proposed SRCQIF-HC will converge to its centralized counterpart shown in equations (15) and (16) as $L$ goes to $\infty$ .

Proof

In the light of Lemma 1, one has $\underset{L\to \infty }{lim}{\pi }_{i,j}^L=1/N$ , where ${\pi }_{i,j}^L$ is the (i, j)th element of ${\mathit{\Pi }}^L$ . Then, as $L\to \infty$ , one has

$${\hat{\mathit{y}}}_{i,k|k-1}^L=\frac{1}{N}\sum _{j\in \mathcal{N}}{\hat{\mathit{y}}}_{j,k|k-1},{\mathit{Y}}_{i,k|k-1}^L=\frac{1}{N}\sum _{j\in \mathcal{N}}{\mathit{Y}}_{j,k|k-1}$$

(52)

and

$${\mathit{u}}_{i,k}^L=\frac{1}{N}\sum _{j\in \mathcal{N}}{\mathit{u}}_{j,k},{\mathit{U}}_{i,k}^L=\frac{1}{N}\sum _{j\in \mathcal{N}}{\mathit{U}}_{j,k}$$

(53)

Therefore, it is achievable to write equations (40) and (41) as

$$\{\begin{array}{c}{\hat{\mathit{y}}}_{i,k|k}=\frac{1}{N}\sum _{j\in \mathcal{N}}{\hat{\mathit{y}}}_{j,k|k-1}+\frac{{\lambda }_{i,k}}{N}\sum _{j\in \mathcal{N}}{\mathit{u}}_{j,k}\hfill \\ {\mathit{Y}}_{i,k|k}=\frac{1}{N}\sum _{j\in \mathcal{N}}{\mathit{Y}}_{j,k|k-1}+\frac{{\lambda }_{i,k}}{N}\sum _{j\in \mathcal{N}}{\mathit{U}}_{j,k}\hfill \end{array}$$

(54)

Suppose each node enjoys identical prior estimates, that is ${\hat{\mathit{y}}}_{j,k|k-1}={\hat{\mathit{y}}}_{k|k-1}$ , ${\mathit{Y}}_{j,k|k-1}={\mathit{Y}}_{k|k-1}$ . Due to choice of scalar weights ${\lambda }_{i,k}=N$ , one has

$${\hat{\mathit{y}}}_{i,k|k}={\hat{\mathit{y}}}_{k|k-1}+\sum _{j\in \mathcal{N}}{\mathit{u}}_{j,k},{\mathit{Y}}_{i,k|k}={\mathit{Y}}_{k|k-1}+\sum _{j\in \mathcal{N}}{\mathit{U}}_{j,k}$$

(55)

For communication nodes, there is no novel information available, and hence, ${\mathit{u}}_{j,k}=0,{\mathit{U}}_{j,k}=0$ for $j\in \mathcal{C}$ . Then, equation (55) is of the same form as equations (15) and (16) in the centralized scheme.

It is shown in Theorem 2 that the proposed SRCQIF-HC can achieve the performance of its centralized counterpart with an infinite number of consensus iterations. However, in practical applications, only a finite, even small number of consensus iterations are accessible. In such a case, the choice with ${\lambda }_{i,k}=N$ may cause some nodes to overestimate the novel information,^3,23,37 which should be avoided to preserve consistency of the estimates. An alternative solution can be exploited to compute a normalization factor in a distributed manner, which is capable of improving the filtering performance as well as maintaining consistency of local estimates. In a sensor network, not all nodes can observe the target of interest, and hence, it is essential to consider whether a node is equipped with observability or not. It is feasible to determine ${\lambda }_{i,k}$ via consensus as follows ¹⁶

$${\gamma }_{i,k}^{l+1}=\sum _{j\in {\mathcal{N}}_i}{\pi }_{i,j}{\gamma }_{j,k}^l,l=0,1,\dots ,L-1$$

(56)

where ${\gamma }_{i,k}^0=1$ if $i\in \mathcal{S}$ and ${\gamma }_{i,k}^0=0$ otherwise. The choice

$${\lambda }_{i,k}=\{\begin{array}{cc}\frac{1}{{\gamma }_{i,k}^L}\hfill & \mathrm{if}{\gamma }_{i,k}^L\ne 0\hfill \\ \multicolumn{1}{c}{1}& \mathrm{otherwise}\hfill \end{array}$$

(57)

ensures that the novel information for each node is never overestimated. Actually, the weight ${\lambda }_{i,k}{\pi }_{i,j}^L={\pi }_{i,j}^L/{\gamma }_{i,k}^L$ , which is exploited to multiply ${\mathit{u}}_{j,k}$ and ${\mathit{U}}_{j,k}$ in equations (40) and (41), is guaranteed to be no more than 1 by construction of ${\gamma }_{i,k}^L=\sum _{j\in \mathcal{S}}{\pi }_{i,j}^L$ . This choice will produce satisfactory performance even a small $L$ is provided. ³⁸ In addition, ${\pi }_{i,j}$ is also a crucial parameter for the convergence rate of consensus algorithms. A poor choice may lead to consensus on local information with more iterations, which is not always practical in real-time applications. Here, the Metropolis weights ³⁹ are adopted to ensure fast consensus and convergence condition for $\mathit{\Pi }$ required in Lemma 1

$${\pi }_{\mathrm{ij}}=\{\begin{array}{cc}{(1+max\{d_i,d_j\})}^{-1}\hfill & \mathrm{if}\{i,j\}\in \mathcal{E}\hfill \\ \multicolumn{1}{c}{1-\sum _{j\in {\mathcal{N}}_i}{\pi }_{\mathrm{ij}}}& \mathrm{if}i=j\hfill \\ \multicolumn{1}{c}{0}& \mathrm{otherwise}\hfill \end{array}$$

Consistency of local estimates

The stability of the proposed SRCQIF-HC is analyzed in terms of consistency. The square-root filter is mathematically equivalent to its standard counterpart, and it is used to overcome the numerical problems resulting from finite word length processor. It has been proved in Simon ⁴⁰ that the stability of square-root filter is exactly equivalent to that of its standard counterpart. To discuss the consistency analysis of the proposed SRCQIF-HC, the statistical linear error propagation methodology^20,29 is applied here to linearize the nonlinear system.

The cross-covariance between the previous estimate and the current prediction can be written as

$$\begin{array}{cc}\multicolumn{1}{c}{{\mathit{P}}_{i,x_{k-1|k-1},x_{k|k-1}}}& =\mathbb{E}\left\{({\mathit{x}}_{k-1}-{\hat{\mathit{x}}}_{i,k-1|k-1}){({\mathit{x}}_k-{\hat{\mathit{x}}}_{i,k|k-1})}^T\right\}\hfill \\ \multicolumn{1}{c}{}& \approx \mathbb{E}\left\{(x_{k-1}-{\hat{\mathit{x}}}_{i,k-1|k-1}){[{\mathcal{F}}_{i,k-1}({\mathit{x}}_{k-1}-{\hat{\mathit{x}}}_{i,k-1|k-1})+{\mathit{w}}_{k-1}]}^{\mathrm{T}}\right\}\hfill \\ \multicolumn{1}{c}{}& ={\mathit{P}}_{i,k-1|k-1}{\mathcal{F}}_{i,k-1}^{\mathrm{T}}\hfill \end{array}$$

(58)

Then, the pseudo transition matrix is

$${\mathcal{F}}_{i,k-1}={\mathit{P}}_{i,x_{k-1|k-1},x_{k|k-1}}^{\mathrm{T}}{\mathit{Y}}_{i,k-1|k-1}$$

(59)

where ${\mathit{P}}_{i,x_{k-1|k-1},x_{k|k-1}}$ can be approximated by

$$\begin{array}{c}{\mathit{P}}_{i,x_{k-1|k-1},x_{k|k-1}}=\hfill \\ \sum _{s=1}^{N_p}{w}_s({\mathit{\chi }}_{i,s,k-1|k-1}-{\hat{\mathit{x}}}_{i,k-1|k-1}){(\mathit{f}\left({\mathit{\chi }}_{i,s,k-1|k-1}\right)-{\hat{\mathit{x}}}_{i,k|k-1})}^{\mathrm{T}}\hfill \end{array}$$

(60)

The pseudo measurement matrix ${\mathcal{H}}_{i,k}$ can be computed by equation (26). Therefore, the system equations (1) and (2) can be rewritten as

$$\{\begin{array}{c}{\mathit{x}}_k={\mathit{\alpha }}_{i,k-1}{\mathcal{F}}_{i,k-1}{\mathit{x}}_{k-1}+{\mathit{w}}_{k-1}\hfill \\ \multicolumn{1}{c}{{\mathit{z}}_{i,k}={\mathit{\beta }}_{i,k}{\mathcal{H}}_{i,k}{\mathit{x}}_k+{\mathit{v}}_{i,k}}\end{array}$$

(61)

where the unknown instrumental matrix ${\mathit{\alpha }}_{i,k-1}=\mathrm{diag}\left({\mathit{\alpha }}_{i,k-1}^1,\dots ,{\mathit{\alpha }}_{i,k-1}^{n_x}\right)$ and ${\mathit{\beta }}_{i,k}=\mathrm{diag}({\mathit{\beta }}_{i,k}^1,\dots ,{\mathit{\beta }}_{i,k}^{n_z})$ are introduced here to neutralize the possible approximation errors.

According to equation (61), ${\mathit{Y}}_{i,k|k-1}$ , ${\mathit{u}}_{i,k}$ , and ${\mathit{U}}_{i,k}$ can be, respectively, written as

$$\{\begin{array}{l}{\mathit{Y}}_{i,k|k-1}={\left({\mathit{\alpha }}_{i,k-1}{\mathcal{F}}_{i,k-1}{\mathit{Y}}_{i,k-1|k-1}^{-1}{\mathcal{F}}_{i,k-1}^{\mathrm{T}}{\mathit{\alpha }}_{i,k-1}^{\mathrm{T}}+{\mathit{Q}}_{k-1}\right)}^{-1}\hfill \\ {\mathit{u}}_{i,k}={\left({\mathit{\beta }}_{i,k}{\mathcal{H}}_{i,k}\right)}^{\mathrm{T}}{\mathit{R}}_{i,k}^{-1}{\mathit{z}}_{i,k}\hfill \\ {\mathit{U}}_{i,k}={\left({\mathit{\beta }}_{i,k}{\mathcal{H}}_{i,k}\right)}^{\mathrm{T}}{\mathit{R}}_{i,k}^{-1}{\mathit{\beta }}_{i,k}{\mathcal{H}}_{i,k}\hfill \end{array}$$

(62)

For simplicity, ${\mathit{Y}}_{i,k|k-1}$ can be rewritten as

$${\mathit{Y}}_{i,k|k-1}={\Psi }_{k-1}\left({\mathit{Y}}_{i,k-1|k-1}\right)$$

(63)

where ${\Psi }_{k-1}\left({\mathit{Y}}_{i,k-1|k-1}\right)={\left({\mathit{\alpha }}_{i,k-1}{\mathcal{F}}_{i,k-1}{\mathit{Y}}_{i,k-1|k-1}^{-1}{\mathcal{F}}_{i,k-1}^{\mathrm{T}}{\mathit{\alpha }}_{i,k-1}^{\mathrm{T}}+{\mathit{Q}}_{k-1}\right)}^{-1}$ .

Definition 1

Consider a random vector $\mathit{x}$ .^12,13,41 Let $\hat{\mathit{x}}$ and $\mathit{P}$ be, respectively, the estimate of $\mathit{x}$ and the estimate of the corresponding error covariance. Then the pair $(\hat{\mathit{x}},\mathit{P})$ is said to be consistent if $\mathbb{E}\left\{\left(\widehat{\mathit{x}}-\mathit{x}\right){\left(\widehat{\mathit{x}}-\mathit{x}\right)}^{\mathrm{T}}\right\}\le \mathit{P}$ .

It indicates in Definition 1 that the true error covariance is always upper bounded (in the positive definite sense) by the estimated error covariance. From the perspective of information filter, the information pair $(\hat{\mathit{y}},\mathit{Y})=({\mathit{P}}^{-1}\hat{\mathit{x}},{\mathit{P}}^{-1})$ is consistent if $\mathit{Y}\le (\mathbb{E}\{({\mathit{Y}}^{-1}\hat{\mathit{y}}-\mathit{x})({\mathit{Y}}^{-1}\hat{\mathit{y}}-\mathit{x}{)}^{\mathrm{T}}\}{)}^{-1}$ .

Lemma 2

For any positive semidefinite matrices ${\mathit{Y}}_1\le {\mathit{Y}}_2$ , ¹² one has $0\le {\Psi }_k({\mathit{Y}}_1)\le {\Psi }_k({\mathit{Y}}_2)$ .

Theorem 3

If the initial estimate $({\hat{\mathit{y}}}_{i,0|0},{\mathit{Y}}_{i,0|0})$ for $i\in \mathcal{N}$ is consistent, the information pair $({\hat{\mathit{y}}}_{i,k|k},{\mathit{Y}}_{i,k|k})$ obtained by the proposed SRCQIF-HC is consistent in that

$${\mathit{Y}}_{i,k|k}\le {\left(\mathbb{E}\left\{({\hat{\mathit{x}}}_{i,k|k}-{\mathit{x}}_k){({\hat{\mathit{x}}}_{i,k|k}-{\mathit{x}}_k)}^{\mathrm{T}}\right\}\right)}^{-1}$$

(64)

with ${\hat{\mathit{x}}}_{i,k|k}={\mathit{Y}}_{i,k|k}^{-1}{\hat{\mathit{y}}}_{i,k|k}$ .

Proof

See the Appendix 1.

Note that it is easy for the initial estimate of each node to be consistent. In the worst case where no prior information is available, the initial values can be chosen as ${\hat{\mathit{y}}}_{i,0|0}=0$ , ${\mathit{Y}}_{i,0|0}=0$ for any $i\in \mathcal{N}$ . It is shown in Theorem 3 that the proposed SRCQIF-HC avoids unaware reuse of the identical information and, hence, preserves consistency of local estimates.

Simulation experiments

In this section, a typical nonlinear air traffic control scenario^10,17 is considered to test the performance of the proposed SRCQIF-HC. The nonlinear dynamic is modeled by

${\mathit{x}}_k=\left(\begin{array}{ccccc}1\hfill & \frac{\sin \mathrm{\Omega }T}{\Omega }\hfill & 0\hfill & -\left(\frac{1-\cos \Omega T}{\Omega }\right)\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& \cos \Omega T\hfill & 0\hfill & -\sin \Omega T\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& \frac{1-\cos \Omega T}{\Omega }\hfill & 1\hfill & \frac{\sin \Omega T}{\Omega }\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& \sin \Omega T\hfill & 0\hfill & \cos \Omega T\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 1\hfill \end{array}\right){\mathit{x}}_{k-1}+{\mathit{w}}_{k-1}$

where ${\mathit{x}}_k=[x_k,{\stackrel{·}{x}}_k,y_k,{\stackrel{·}{y}}_k,\Omega {]}^{\mathrm{T}}$ is the state vector, $(x_k,y_k)$ is the target position, and $({\stackrel{·}{x}}_k,{\stackrel{·}{y}}_k)$ denotes the velocity. $\Omega$ is the unknown turn rate. $T=1\mathrm{s}$ is the sampling interval. ${\mathit{w}}_{k-1}$ is the zero-mean Gaussian noise with covariance $\mathit{Q}=\mathrm{diag}(q_1\mathit{M},q_1\mathit{M},q_{2}T)$ . Same as Arasaratnam and Haykin, ¹⁰ the parameters are set as follows: $q_1=0.1{\mathrm{m}}^2/{\mathrm{s}}^3$ , $q_2={10}^{-6}/{\mathrm{s}}^3$ , and

$$\mathit{M}=\left(\begin{array}{cc}\frac{T^3}{3}\hfill & \frac{T^2}{2}\hfill \\ \multicolumn{1}{c}{\frac{T^2}{2}}& T\hfill \end{array}\right)$$

Each simulation runs $K=100$ steps. The initial estimate is randomly generated from $N({\mathit{x}}_0;{\hat{\mathit{x}}}_0,{\mathit{P}}_0)$ , where $x_0=[1000\mathrm{m},300\mathrm{m}/\mathrm{s},1000\mathrm{m},0\mathrm{m}/\mathrm{s},-3°/\mathrm{s}]$ and ${\mathit{P}}_0=\mathrm{diag}(100{\mathrm{m}}^2,25{\mathrm{m}}^2/{\mathrm{s}}^2,100{\mathrm{m}}^2,25{\mathrm{m}}^2/{\mathrm{s}}^2,100\mathrm{m}\mathrm{ra}{\mathrm{d}}^2/{\mathrm{s}}^2)$ are the same as that in Arasaratnam and Haykin. ¹⁰

Similar to Chen et al., ¹⁷ there are 10 sensors randomly located at the [−4000 m, 4000 m] × [−4000 m, 4000 m] area. The network consists of four sensor nodes and six communication nodes. The communication topology is shown in Figure 1. The measurement equation is modeled by

$${\mathit{z}}_{i,k}=\left[\begin{array}{c}\sqrt{{(x_k-x_i)}^2+{(y_k-y_i)}^2}\hfill \\ \multicolumn{1}{c}{a\tan \left(\frac{y_k-y_i}{x_k-x_i}\right)}\end{array}\right]+{\mathit{v}}_{i,k}$$

Figure 1.

An illustrated sensor network with sparse connectivity.

where $(x_i,y_i)$ denotes the sensor position. $v_{i,k}$ is the zero-mean Gaussian noise with covariance ${\mathit{R}}_{i,k}=\mathrm{diag}({\sigma }_r^2,{\sigma }_{\theta }^2)$ , where ${\sigma }_r=10\mathrm{m}$ and ${\sigma }_{\theta }=\sqrt{10}\mathrm{mrad}$ are the same as that in Wang et al. ²³

We carry out $M=100$ independent experiments, and the total number of consensus iterations is set to $L=10$ . The Cub-ICF, ¹⁴ HCCKF, ¹⁷ HCMCI-UKF, ⁸ HCMCI-CKF, ²³ SCICF, ¹⁵ the proposed third-degree SRCQIF-HC (SRTCQIF-HC), and fifth-degree SRCQIF-HC (SRFCQIF-HC) are considered in this simulation for comparison. The Cramér–Rao lower bound (CRLB)^42,43 is introduced here as an indication of the best possible performance that can be achieved. The suggested optimal scaling parameter $\kappa =3-n=-2$ for UKF leads to negative weights, which may halt its operation, ¹⁰ Here, we set $\kappa =1$ as used in Jia et al. ²⁸

Performance metrics

Several standard performance metrics are used here to compare the accuracy, consistency, consensus, and computational cost of the distributed algorithms under discussion.

1. Root mean square error (RMSE): The RMSE indicates the accuracy of an estimation algorithm. ²³

2. Normalized estimation error squared (NEES): The NEES is a measure to check for filter consistency. ⁴⁴

3. Averaged consensus estimate error (ACEE): The ACEE suggests the estimation bias between different nodes. ⁴⁵

4. Computational cost: The computational cost is measured by the average running time (RT) of one-step filtering over all time steps. ¹⁵

The RMSE behaviors in position, velocity, and turn rate with a single consensus iteration are shown in Figures 2–4. It can be seen that the SCICF and Cub-ICF are less accurate than other distributed algorithms discussed in the scenario. HCCKF and HCMCI-UKF, which are based on the conventional information filtering paradigm, are less accurate than HCMCI-CKF, SRTCQIF-HC, and SRFCQIF-HC, which update the posterior estimates with the novel scheme presented in section “A novel Gaussian approximation information filtering paradigm.” In addition, the proposed SRTCQIF-HC and SRFCQIF-HC exhibits higher estimation accuracy than HCMCI-CKF. To further investigate the estimation accuracy with different consensus iterations, the average RMSE in position is shown in Figure 5. As consensus iteration increases, the estimation accuracy of all algorithms is improved. In specific, HCCKF and HCMCI-UKF have comparable performance but are much better than SCICF and Cub-ICF. The performance of HCMCI-CKF and SRTCQIF-HC are also analogous to each other but much less accurate than that of SRFCQIF-HC derived from higher quadrature rules, which approaches closer to the CRLB with increasing consensus iterations. To facilitate a more comprehensive analysis on estimation accuracy, the average RMSEs in position, velocity, and turn rate are given in Tables 1–3, which indicate that the proposed SRFCQIF-HC achieves higher accuracy than the algorithms compared.

Figure 2.

The RMSE in position with a single iteration.

Figure 3.

The RMSE in velocity with a single iteration.

Figure 4.

The RMSE in turn rate with a single iteration.

Figure 5.

The average RMSE in position with different number of consensus iterations.

Table 1.

The average RMSE in position for different algorithms.

RMSE (m)	Cub-ICF	HCCKF	HCMCI-CKF	HCMCI-UKF	SCICF	SRTCQIF-HC	SRFCQIF-HC	CRLB
$L=1$	15.60	14.88	13.40	14.65	18.63	13.31	12.81	8.75
$L=3$	13.40	12.73	11.61	12.54	16.47	11.36	10.85	8.75
$L=5$	12.44	11.95	10.83	11.75	15.60	10.57	10.06	8.75
$L=10$	11.27	10.99	9.87	10.79	14.38	9.74	9.23	8.75

RMSE: root mean square error; Cub-ICF: cubature information consensus filter; HCCKF: hybrid consensus–based cubature Kalman filter; CKF: cubature Kalman filter; UKF: unscented Kaman filter; SCICF: square-root version cubature IWCF; SRTCQIF-HC: third-degree SRCQIF-HC; SRFCQIF-HC: fifth-degree SRCQIF-HC; CRLB: Cramér–Rao lower bound; IWCF: information weighted consensus filter; SRCQIF-HC: hybrid consensus–based square-root cubature quadrature information filter.

Table 2.

The average RMSE in velocity for different algorithms.

RMSE (m/s)	Cub-ICF	HCCKF	HCMCI-CKF	HCMCI-UKF	SCICF	SRTCQIF-HC	SRFCQIF-HC	CRLB
$L=1$	10.41	9.71	9.21	9.76	11.49	9.21	8.83	7.31
$L=3$	9.70	9.12	8.61	9.07	11.00	8.59	8.21	7.31
$L=5$	9.29	8.88	8.37	8.80	10.80	8.34	7.96	7.31
$L=10$	8.75	8.56	8.04	8.51	10.50	8.06	7.68	7.31

RMSE: root mean square error; Cub-ICF: cubature information consensus filter; HCCKF: hybrid consensus–based cubature Kalman filter; CKF: cubature Kalman filter; UKF: unscented Kaman filter; SCICF: square-root version cubature IWCF; SRTCQIF-HC: third-degree SRCQIF-HC; SRFCQIF-HC: fifth-degree SRCQIF-HC; CRLB: Cramér–Rao lower bound; IWCF: information weighted consensus filter; SRCQIF-HC: hybrid consensus–based square-root cubature quadrature information filter.

Table 3.

The average RMSE in turn rate for different algorithms.

RMSE (deg/s)	Cub-ICF	HCCKF	HCMCI-CKF	HCMCI-UKF	SCICF	SRTCQIF-HC	SRFCQIF-HC	CRLB
$L=1$	1.327	1.267	1.237	1.266	1.343	1.237	1.232	1.196
$L=3$	1.287	1.250	1.219	1.248	1.331	1.219	1.214	1.196
$L=5$	1.270	1.242	1.212	1.241	1.327	1.212	1.207	1.196
$L=10$	1.248	1.233	1.204	1.232	1.320	1.204	1.199	1.196

RMSE: root mean square error; Cub-ICF: cubature information consensus filter; HCCKF: hybrid consensus–based cubature Kalman filter; CKF: cubature Kalman filter; UKF: unscented Kaman filter; SCICF: square-root version cubature IWCF; SRTCQIF-HC: third-degree SRCQIF-HC; SRFCQIF-HC: fifth-degree SRCQIF-HC; CRLB: Cramér–Rao lower bound; IWCF: information weighted consensus filter; SRCQIF-HC: hybrid consensus–based square-root cubature quadrature information filter.

Figure 6 shows the average NEES curves of different algorithms with 95% confidence interval. The curve of Cub-ICF lies much higher than the concentration region, which shows poor consistency of local estimates. On the other hand, most points of the NEES curves of the remaining algorithms lie either below or within the concentration region, which demonstrates improved consistency. The ACEE of different algorithms is shown in Figure 7. With increasing consensus iterations, better consensus on local estimates is shown through lower ACEE. Moreover, the proposed SRFCQIF-HC has much lower ACEE than all the other algorithms, which indicates that SRFCQIF-HC achieves the fastest consensus rate among the discussed algorithms.

Figure 6.

The NEES of different algorithms with a single iteration.

Figure 7.

The ACEE with different consensus iterations.

The computational cost measured by average RT for all the algorithms is shown in Table 4. As higher cubature quadrature rules are exploited, more sampling points are required during the recursion of state estimation. Although it introduces extra computational burden on the system, the proposed SRFCQIF-HC possesses outstanding performance in the aspects of estimation accuracy, consistency of local estimates, and consensus on estimates across the entire sensor network.

Table 4.

Comparison of the running time for a single filtering step.

RT (×10−3 s)	Cub-ICF	HCCKF	HCMCI-CKF	HCMCI-UKF	SCICF	SRTCQIF-HC	SRFCQIF-HC
$L=1$	1.80	1.92	1.96	2.02	2.30	2.80	4.34
$L=3$	2.06	2.48	2.50	2.58	3.16	4.52	5.97
$L=5$	2.33	3.04	3.06	3.12	3.98	6.21	7.66
$L=10$	3.01	4.37	4.39	4.50	6.09	10.4	11.9

RT: running time; Cub-ICF: cubature information consensus filter; HCCKF: hybrid consensus–based cubature Kalman filter; CKF: cubature Kalman filter; UKF: unscented Kaman filter; SCICF: square-root version cubature IWCF; SRTCQIF-HC: third-degree SRCQIF-HC; SRFCQIF-HC: fifth-degree SRCQIF-HC; IWCF: information weighted consensus filter; SRCQIF-HC: hybrid consensus–based square-root cubature quadrature information filter.

Conclusion

In this article, to tackle nonlinear state estimation problems with networked sensors, a novel distributed SRCQIF-HC is proposed. The proposed SRCQIF-HC preserves the positive features, including consistency and accuracy, of both CI and CM algorithms and is mathematically equivalent to its centralized counterpart with resort to the novel measurement information update strategy presented in section “A novel Gaussian approximation information filtering paradigm.” Stability analysis with regard to consensus on local estimates, convergence to its centralized counterpart, and consistency is also developed. The simulation results indicate that the proposed SRCQIF-HC achieves an improved performance than the existing algorithms at the cost of a little more computational burden. Here, both the process and measurement noises are assumed to be Gaussian; in future work, we will extend the proposed algorithm to be applicable to nonlinear filtering under non-Gaussian process or measurement noises.