State Estimation of Wireless Sensor Networks in the Presence of Data Packet Drops and Non-Gaussian Noise

The concept of correntropy, first proposed by Principe et al.[26], is a very practical method for evaluating the generalized similarity between random variables ${X,Y\in R}$ with the same dimensions. Here, correntropy is defined as

where E[.] is the expectation operator, ${\operatorname{V}(.)}$ represents the information potential, ${{\operatorname{F}_{XY}}(x,y)}$ represents the probability distribution function (PDF) with respect to variables ${X}$ and ${Y}$, and ${\kappa\left({\cdot{\rm{}},{\rm{}}\cdot}\right)}$ is the shift-invariant Mercer kernel. Here, we employ the Gaussian kernel, which is given as

where ${e=x-y}$ represents the error between elements ${x}$ and ${y}$, and ${\sigma>0}$ represents the kernel bandwidth (or kernel size) of the Gaussian kernel function.

However, only a limited amount of data related to the variables ${X}$ and ${Y}$ can be obtained in realistic scenarios, and the PDF ${{\operatorname{F}_{XY}}(x,y)}$ is usually unknown. Under these conditions, a sample estimator can be used to calculate the correntropy as follows:

where

and ${N}$ samples are employed to define ${{\operatorname{F}_{XY}}(x,y)}$.

Applying a Taylor expansion to the Gaussian kernel function yields the following:

We can infer from (5) that the correntropy is the weighted sum of all even-order moments of relation ${X-Y}$. Compared with other similarity measurement schemes, such as the mean-square error (MSE) criterion, correntropy contains all even-order moments, and is therefore useful for nonlinear and signal processing applications in non-Gaussian noise environments.

The state-space model of the DKF over the nodes of a WSN under conditions of data packet drops is the first problem that must be solved to estimate the state of a target node. Here, we adopt a powerful state-space model [19], which is described as follows.

Consider a WSN with ${N}$ sensors. The states and observations are as follows:

where ${{{\boldsymbol{x}}_{k}},{{\boldsymbol{x}}_{k-1}}\in{\mathbb{R}^{n\times 1}}}$ represent the states of the dynamical system at instants ${k}$ and ${k-1}$, ${\Omega}$ represents the set of all sensors in the network, ${{\boldsymbol{y}}_{k}^{i}\in{\mathbb{R}^{{m_{i}}\times 1}}({m_{i}}=\operatorname{rank}({{\boldsymbol{C}}_{i}})\forall i\in\Omega)}$ denotes the observations obtained by node ${i}$ at instant ${k}$, ${{{\boldsymbol{A}}_{k}}}$ and ${{{\boldsymbol{C}}^{i}}}$ represent the state transition matrix and observation matrix of the system, respectively, and ${{{\boldsymbol{q}}_{k}}}$ and ${{\boldsymbol{v}}_{k}^{i}}$ are the mutually uncorrelated process noise and measurement noise, respectively, with means of zero and covariances respectively given as follows:

We then define ${{\Omega_{i}}\subset\Omega}$ to represent the set of all neighboring sensors of the ${i}$th sensor. Therefore, set ${{\mho_{i}}={\Omega_{i}}\cup\{i\}}$ contains node ${i}$ itself and all its neighboring nodes. According to the properties of a WSN, the process whereby node ${i}$ receives the observations from its neighboring node ${j}$ at instant ${k}$ can be modeled as follows:

Here, the term ${\gamma_{k}^{i,j}}$ is used to indicate the existence of data packet drops over the WSN. We employ the Bernoulli model for modeling the data packet loss process and assume that the process has independent and identically distributed properties for ${\gamma_{k}^{i,j}}$ with the probability

and

We assume that variables ${\gamma_{k}^{i,j}}$ and ${\gamma_{k}^{j,i}}$ are mutually independent (i.e., ${\forall i\neq j}$), but the condition ${{p_{i,j}}={p_{j,i}}}$ is allowed, and ${\gamma_{k}^{i,j}}$ is independent of the process noise, measurement noise, and initial state of the dynamical system. According to (8), if sensor ${i}$ receives the observations from its neighboring node ${j}$ at instant ${k}$, then ${\gamma_{k}^{i,j}=1}$ and ${\gamma_{k}^{i,j}=0}$ when all the components of ${{{\boldsymbol{y}}_{k}^{j}}}$ are lost. We assume that node ${i}$ can receive all observations from itself at any time, so ${\gamma_{k}^{i,j}\equiv 1}$.

Based on the above analysis, the various factors of a DKF model that must take data packet drops into account can be expressed as

Here, ${{\text{vec}}\left\{\cdot\right\}}$, ${{\text{col}}\left\{\cdot\right\}}$, and ${\operatorname{diag}\left\{\cdot\right\}}$ denote the vectorization operation, the columnization operation, and the (block) diagonalization operation, respectively, and ${{{\boldsymbol{I}}_{{m_{j}}}}(j\in{\mho_{i}})}$ represents an identity matrix of dimension ${{m_{j}}\times{m_{j}}}$. According to the above discussion, the state-space model of the ${i}$th sensor in the case of data packet drops can be written in a compact form as

In addition, the covariance of ${{\boldsymbol{v}}_{k}^{{\mho_{i}}}}$ is

We first derive the stationary DMCKF in the case of data packet drops based on the model given by (12)-(14). The state prediction error of the algorithm can be written as

Combining this with the above state-space model yields the following augmented model:

where ${{{{\boldsymbol{I}}_{n}}}}$ represents an ${n\times n}$ identity matrix and ${{\boldsymbol{g}}_{k}^{i}}$ represents the augmented noise vector containing the state and measurement errors of the dynamical system, which is defined as

Assuming that the covariance matrix of the augmented vector ${{\text{E}}[{\boldsymbol{g}}_{k}^{i}{({\boldsymbol{g}}_{k}^{i})^{\text{T}}}]}$ is positive definite yields the following:

with

Here, ${{\boldsymbol{B}}_{k}^{i}}$ and ${{({\boldsymbol{B}}_{k}^{i})^{\text{T}}}}$ can be obtained by the Cholesky decomposition of ${{\text{E}}[{\boldsymbol{g}}_{k}^{i}{({\boldsymbol{g}}_{k}^{i})^{\text{T}}}]}$, and matrix ${{\boldsymbol{D}}_{p;k}^{{\mho_{i}}}}$ represents the expectation of ${{{\boldsymbol{D}}_{\gamma;k}^{{\mho_{i}}}}}$. Left multiplying each term of (17) by ${{({\boldsymbol{B}}_{k}^{i})^{-1}}}$ yields the following:

where

According to the above derivation, we propose the following cost function based on the MC criterion:

where ${d_{k}^{i;h}}$ represents the ${h}$th element of ${{\boldsymbol{D}}_{k}^{i}}$, ${{\boldsymbol{w}}_{k}^{i;h}}$ represents the ${h}$th row of ${{\boldsymbol{W}}_{k}^{i}}$, and

denotes the number of elements of ${{\boldsymbol{D}}_{k}^{i}}$. Then, the objective function of the optimal state estimation ${{\boldsymbol{x}}_{k}^{i}}$ based on the MC criterion is

where ${e_{k}^{i;h}=d_{k}^{i;h}-{\boldsymbol{w}}_{k}^{i;h}{\boldsymbol{x}}_{k}^{i}}$ and ${\left({h=1,2,\cdots L}\right)}$ is the ${h}$th element of ${e_{k}^{i}}$. Finally, the optimal state estimation of ${{\boldsymbol{\hat{x}}}_{k}^{i}}$ is achieved by maximizing the information potential ${\operatorname{V}(.)}$. To that end, we set the gradient of the cost function ${{{J_{L}}({\boldsymbol{x}}_{k}^{i})}}$ regarding ${{{\boldsymbol{x}}_{k}^{i}}}$ to zero:

and obtaining the optimal state of ${{{\boldsymbol{x}}_{k}^{i}}}$ is relatively easy:

Because ${e_{k}^{i;h}=d_{k}^{i;h}-{\boldsymbol{w}}_{k}^{i;h}{\boldsymbol{x}}_{k}^{i}}$ is a function of ${{\boldsymbol{x}}_{k}^{i}}$, the optimal state estimation in (27) is a fixed-point iterative equation of ${{\boldsymbol{x}}_{k}^{i}}$, and can be rewritten in the following form:

where

According to the above derivation, the iterative equation in (28) can be written as follows:

Here, ${{({\boldsymbol{\hat{x}}}_{k}^{i})_{t+1}}}$ represents the result of ${{\boldsymbol{x}}_{k}^{i}}$ at the fixed-point iteration ${t+1}$, and (27) can be further written in the form of matrix multiplication:

where

According to (22), (31), and (32), we obtain the following:

We then apply the following matrix inversion lemma [28]

with the identifications

and obtain the expression (36) given below.

According to (22), (31), and (32), we obtain the following:

Substituting formulas (36) and (37) into (31) yields

where

Remark 1. Equation (38) is the optimal solution of ${{\boldsymbol{x}}_{k}^{i}}$, and it depends on the prior estimate ${{\boldsymbol{\hat{x}}}_{k|k-1}^{i}}$ and available observation information ${{{\boldsymbol{s}}_{k}^{{\mho_{i}}}}}$. The observation information obtainable by the ${i}$th node is determined by the matrix ${{{\boldsymbol{D}}_{\gamma;k}^{{\mho_{i}}}}}$.

According to the above derivations, we summarize the steps of the proposed DMCKF-DPD algorithm as follows.

1. 1.

   Initialize the values of the ${\sigma}$, ${\varepsilon}$ (a small positive number), ${{\boldsymbol{\hat{x}}}_{0|0}^{i}}$, and ${{\boldsymbol{P}}_{0|0}^{i}}$.

2. 2.

   Apply (19) to calculate ${{\boldsymbol{B}}_{k}^{i}}$, and calculate ${{\boldsymbol{\hat{x}}}_{k|k-1}^{i}}$ and ${{\boldsymbol{P}}_{k|k-1}^{i}}$ as follows:

   $$\begin{split}\begin{gathered}{\boldsymbol{x}}_{k}^{i}={\boldsymbol{Ax}}_{k-1}^{i},\hfill\\ {\boldsymbol{P}}_{k|k-1}^{i}={\boldsymbol{AP}}_{k-1|k-1}^{i}{{\boldsymbol{A}}^{\text{T}}}+{\boldsymbol{Q}}.\hfill\\ \end{gathered}\end{split}$$

   (40)

3. 3.

   Update the parameters ${{\boldsymbol{D}}_{k}^{i}={\boldsymbol{W}}_{k}^{i}{\boldsymbol{x}}_{k}^{i}+{\boldsymbol{e}}_{k}^{i}}$ in the new system.

4. 4.

   Utilize the fixed-point iterative algorithm to estimate the state of system as follows:

   $$\begin{split}\begin{gathered}{\left({{\boldsymbol{\hat{x}}}_{k|k}^{i}}\right)_{t+1}}\hfill\\ ={\boldsymbol{\hat{x}}}_{k|k-1}^{i}+{\boldsymbol{\bar{K}}}_{k}^{i}\left[{{\boldsymbol{s}}_{k}^{{\mho_{i}}}-{\boldsymbol{D}}_{\gamma;k}^{{\mho_{i}}}{{\boldsymbol{C}}^{{\mho_{i}}}}{\boldsymbol{\hat{x}}}_{k|k-1}^{i}}\right],\hfill\\ \end{gathered}\end{split}$$

   (41)

   where

   $$\begin{split}\begin{gathered}{\boldsymbol{\tilde{K}}}_{k}^{i}={\boldsymbol{\tilde{P}}}_{k|k-1}^{i}{\left({{{\boldsymbol{C}}^{{\mho_{i}}}}}\right)^{\text{T}}}{\left({{\boldsymbol{D}}_{\gamma;k}^{{\mho_{i}}}}\right)^{\text{T}}}\times\hfill\\ \left[{{\boldsymbol{D}}_{\gamma;k}^{{\mho_{i}}}{{\boldsymbol{C}}^{{\mho_{i}}}}{\boldsymbol{\bar{P}}}_{k|k-1}^{i}{{\left({{{\boldsymbol{C}}^{{\mho_{i}}}}}\right)}^{\text{T}}}{{\left({{\boldsymbol{D}}_{\gamma;k}^{{\mho_{i}}}}\right)}^{\text{T}}}+{\boldsymbol{\tilde{R}}}_{k}^{i}}\right],\hfill\\ \end{gathered}\end{split}$$

   (42)

   $$\begin{split}{\boldsymbol{\tilde{P}}}_{k|k-1}^{i}={\boldsymbol{B}}_{P\left({k|k-1}\right)}^{i}{\left({{\boldsymbol{\tilde{\Lambda}}}_{x;k}^{i}}\right)^{-1}}{\left({{\boldsymbol{B}}_{P\left({k|k-1}\right)}^{i}}\right)^{\text{T}}},\end{split}$$

   (43)

   $$\begin{split}{\boldsymbol{\tilde{R}}}_{k}^{i}={\boldsymbol{B}}_{R;k}^{i}{\left({{\boldsymbol{\tilde{\Lambda}}}_{y;k}^{i}}\right)^{-1}}{\left({{\boldsymbol{B}}_{R;k}^{i}}\right)^{\text{T}}},\end{split}$$

   (44)

   $$\begin{split}{\boldsymbol{\tilde{\Lambda}}}_{x;k}^{i}=\operatorname{diag}\left[{{\operatorname{G}_{\sigma}}\left({\tilde{e}_{k}^{i;1}}\right),\cdots,{\operatorname{G}_{\sigma}}\left({\tilde{e}_{k}^{i;n}}\right)}\right],\end{split}$$

   (45)

   $$\begin{split}{\boldsymbol{\tilde{\Lambda}}}_{y;k}^{i}=\operatorname{diag}\left[{{\operatorname{G}_{\sigma}}\left({\tilde{e}_{k}^{i;n+1}}\right),\cdots,{\operatorname{G}_{\sigma}}\left({\tilde{e}_{k}^{i;n+{m_{i;a}}}}\right)}\right],\end{split}$$

   (46)

   $$\begin{split}\tilde{e}_{k}^{i;h}=d_{k}^{i;h}-{\boldsymbol{w}}_{k}^{i;h}{\left({{\boldsymbol{\hat{x}}}_{k|k}^{i;h}}\right)_{t}}.\end{split}$$

   (47)

5. 5.

   Compare the state values obtained by the fixed-point iterative algorithm at ${t}$ and ${t+1}$. If these two state estimates satisfy the following criterion:

   $$\begin{split}\frac{{\left\|{{{\left({{\boldsymbol{\hat{x}}}_{k|k}^{i}}\right)}_{t+1}}-{{\left({{\boldsymbol{\hat{x}}}_{k|k}^{i}}\right)}_{t}}}\right\|}}{{{{\left({{\boldsymbol{\hat{x}}}_{k|k}^{i}}\right)}_{t}}}}\leqslant\varepsilon,\end{split}$$

   (48)

   where ${\varepsilon}$ is the termination condition of the fixed-point iterative algorithm, set the optimal estimate of ${{\boldsymbol{\hat{x}}}_{k|k}^{i}}$ equal to ${{({\boldsymbol{\hat{x}}}_{k|k}^{i})_{t+1}}}$, and continue to step 6). Otherwise, set ${t\leftarrow t+1}$ and return to step 4).

6. 6.

   Update the covariance matrix as follows:

   $$\begin{split}\begin{gathered}{\boldsymbol{P}}_{k|k}^{i}=\left({{\boldsymbol{I}}-{\boldsymbol{\tilde{K}}}_{k}^{i}{\boldsymbol{D}}_{\gamma;k}^{{\mho_{i}}}{{\boldsymbol{C}}^{{\mho_{i}}}}}\right){\boldsymbol{P}}_{k|k-1}^{i}\hfill\\ \left[{{\boldsymbol{I}}-{\boldsymbol{\tilde{K}}}_{k}^{i}{{\left({{{\boldsymbol{C}}^{{\mho_{i}}}}}\right)}^{\text{T}}}{{\left({{\boldsymbol{D}}_{\gamma;k}^{{\mho_{i}}}}\right)}^{\text{T}}}}\right]+{\boldsymbol{\tilde{K}}}_{k}^{i}{\boldsymbol{R}}_{k}^{{\mho_{i}}}{\left({{\boldsymbol{\tilde{K}}}_{k}^{i}}\right)^{\text{T}}},\hfill\\ \end{gathered}\end{split}$$

   (49)

   set ${k\leftarrow k+1}$, and return to step 2).

Remark 2. The kernel bandwidth ${\sigma}$ is an important parameter in the MC criterion. Broadly speaking, the convergence rate of the algorithm increases with increasing ${\sigma}$, but the accuracy of the algorithm decreases with increasing ${\sigma}$.

The computational complexity of the proposed DMCKF-DPD algorithm can be compared with that of the stationary DKF [19] in terms of the equations and operations employed by the two algorithms, which are listed in Table 1.

The stationary DKF algorithm mainly includes (13), (14), and (40) cited in the present paper, and (5) and (6) cited within[19]. Therefore, we can evaluate the computational complexity of the stationary DKF as

The proposed algorithm mainly includes (13) and (14), (40), (2), and (41)-(47). Accordingly, the computational complexity of the DMCKF-DPD algorithm can be defined based on an average number of fixed-point iterative algorithm iterations ${T}$ as

We can infer from this discussion that the computational complexity of the DMCKF-DPD algorithm is moderate compared with that of the stationary DKF, provided that the value of ${T}$ is small, which is indeed the case, as will be demonstrated later in Section IV.

A full analysis of the convergence behavior of the DMCKF-DPD algorithm based on the fixed-point iterative algorithm is very complicated. Therefore, we give only a sufficient condition that ensures the convergence of the fixed-point iterative algorithm. However, a detailed proof process is not presented here because the convergence condition is similar to the analysis presented in an earlier work[28], which can be consulted for additional details.

Theorem 1. First, we assume the conditions ${{\beta_{i}}>{\zeta_{i}}=\frac{{\sqrt{n}\sum\limits_{h=1}^{L}{{{\left\|{{{\left({{\boldsymbol{w}}_{k}^{i;h}}\right)}^{\text{T}}}}\right\|}_{1}}\left|{d_{k}^{i;h}}\right|}}}{{{\lambda_{\min}}\sum\limits_{h=1}^{L}{{{\left({{\boldsymbol{w}}_{k}^{i;h}}\right)}^{\text{T}}}{\boldsymbol{w}}_{k}^{i;h}}}}}$ and ${{\sigma_{i}}\geq\max{\rm{\{}}\sigma_{i}^{*},\sigma_{i}^{\diamondsuit}{\rm{\}}}}$. Here, ${\sigma_{i}^{*}}$ is the optimal result of the equation ${\phi({\sigma_{i}})={\beta_{i}}}$, where