Unreliable networks with random parameter matrices and time-correlated noises: distributed estimation under deception attacks

1.   Introduction

Sensor networks have generated considerable research interest due to their wide range of practical applications and the rapid development of communication and information technologies (see, e.g., ^[1] and ^[2]). A widely-explored research topic in this domain is the distributed estimation problem, which arises in situations where the sensor nodes are spatially distributed according to a predetermined network topology and each sensor node can use data information from itself and its neighboring nodes to estimate the signal of interest. Each sensor node then acts as a local fusion centre, combining its own information with that obtained from adjacent nodes to improve the performance of the local estimators, based solely on its own measurement information. This collaborative signal estimation strategy offers many advantages, including ease of implementation, robustness, scalability and high reliability. In ^[3], the distributed estimation problem is addressed in the presence of random packet dropouts during data transmission. A distributed event-based filtering structure is proposed in ^[4] for a class of sensor networks with sensor saturations and cyber-attacks. Nonlinear systems over sensor networks, whose topologies are changeable subject to Round-Robin protocol in the finite horizon case, are considered in ^[5]. A recursive distributed filtering algorithm for networked systems with random parameter matrices and correlated noises is designed in ^[6]. In recent years, a large variety of distributed estimation algorithms have been designed for sensor networks exposed to various challenges and vulnerabilities, under a range of filter schemes. For example, a new distributed filtering strategy has been presented in ^[7] by fully taking both the prediction estimations and its own measurement innovation into consideration. In ^[8], a novel distributed filtering compensation algorithm is presented in terms of the available transmitted data. A new design of an innovation-based stealthy attack strategy against distributed state estimation over a sensor network is proposed in ^[9], while the state estimation problem in linear time-invariant systems by using a network of distributed observers with switching communication topology is studied in ^[10]. Some significant achievements in the field of distributed estimation for stochastic systems over sensor networks are reviewed in ^[11] and ^[12].

In the fields of physics, electronics and engineering, many applications can involve infinite-step colored measurement noises, particularly when the sampling frequency is high enough to make the noises significantly correlated over two or more consecutive sampling periods. In recent years, researchers have addressed the estimation problem under the assumption that the measurements are affected by infinite-step time-correlated channel noise, modeled as the output of a linear system driven by white noise. Two popular methods for dealing with this type of noise correlation are state augmentation, which is simple and direct but computationally expensive, and measurement differencing, which avoids the problem of increasing dimensions but requires two consecutive measurements to compute the difference. By referencing to the measurement differencing method, in ^[13] the time-correlation of the measurement noises is transformed into the cross-correlation between the equivalent measurement noise and the process noise. Convergence conditions of the optimal linear estimator are obtained in ^[14], by using a new measurement obtained from measurement differencing. In ^[15], by using the time-differencing approach, the available measurements are transformed into an equivalent set of observations that do not depend on the time-correlated noise. Alternative non-augmentation and non-differencing methods to address the state estimation problem, based on the direct estimation of the time-correlated additive noise, are described in ^[16] and ^[17].

Communication networks are typically subject to resource limitations, which can provoke network-related issues during signal measurement or transmission ^[18]. Some such issues, such as the presence of multiplicative noise, missing observations or fading measurements, can be described by introducing stochastic parameter matrices into the measurement equations. In recent years, considerable research has been conducted into estimation problems arising in systems with random parameter matrices. For example, in ^[19] the Tobit Kalman filtering problem is studied for a class of linear discrete-time system with random parameters, where the elements of both the system matrix and the measurement matrix are allowed to be random variables in order to reflect the reality. The distributed fusion estimation problem for networked systems whose multisensor measured outputs involve uncertainties modeled by random parameter matrices is investigated in ^[20]. In ^[21], the optimal linear filtering problem for linear discrete-time stochastic systems with random matrices in both the state and measurement equations is addressed. A Kalman-like recursive distributed optimal linear fusion predictor without feedback, for discrete-time linear stochastic systems with correlated random parameter matrices, is presented in ^[22]. In ^[23], the distributed fusion estimation problem is discussed in the presence of coupled noises, random delays and packet dropouts, for a class of uncertain systems, where the uncertainty in the measurement model is described by random parameter matrices.

In addition to the problem of random uncertainties in measurements and transmissions, a critical issue that cannot be ignored in any study of the estimation problem in networked systems is the possibility of suffering cyber-attacks. Security vulnerability is a common weakness that has been widely discussed in the literature; a bibliographical review of recent advances and challenges in this regard can be found in ^[24]. In particular, deception attacks have attracted significant research attention. This form of attack seeks to compromise data integrity by maliciously and randomly falsifying information. In this respect, ^[25] examined the centralized security-guaranteed filtering problem for linear time-invariant stochastic systems with multirate-sensor fusion subjected to deception attacks. In other approaches, the $H_\infty$-consensus filtering problem for discrete-time systems with multiplicative noises and deception attacks has been investigated by ^[26], and the distributed estimation problem in sensor networks with a specific topology structure under deception attacks has been addressed by ^[27,28,30]. More specifically, the distributed filtering problem in networked systems with fading measurements and multiplicative noises in both the signal and measurement equations is discussed in ^[27]. In ^[28], a positive system over a sensor network with simultaneous deception attacks and various network-induced constraints on sensor measurements is considered. The distributed secure state estimation problem is addressed in ^[29] for a class of general nonlinear systems over sensor networks under unknown deception attacks on innovations, while networked uncertain systems, containing uncertainties due to multiplicative and additive noises in the state and measurement equations, are considered in ^[30].

In view of the above considerations, the present study focuses on a class of networked systems, whose sensor nodes are distributed in space according to a fixed network topology. The main study goal is to address the least-squares linear distributed estimation problem from measurements affected by random parameter matrices and time-correlated additive noises, and simultaneously exposed to random deception attacks. For this purpose, a recursive algorithm for the distributed estimators is generated in two stages. First, each sensor node collects measurements from its neighbors to create intermediate least-squares linear estimators by an innovation approach. In the second stage, the intermediate estimators from neighboring sensors are combined to form distributed estimators through least-squares matrix-weighted linear combinations. A greater volume of information from different sensors is used in the second stage than in the first. This enhances the intermediate estimation performance and reduces disagreements among intermediate estimators from different sensors, by steering each distributed estimator closer to the global optimal linear estimator (hypothetically based on measurements from all network sensors).

This study makes the following main contributions. First, the consideration of random parameter matrices in the measurement equations provides a unified framework for handling common network-induced phenomena, such as multiplicative noise, missing observations and missing or fading measurements; thus, the proposed algorithm can be applied to a wide range of network systems with incomplete information. The model also integrates stochastic deception attacks, to which many networked systems are vulnerable. Second, both the random parameter matrices and the random variables modeling the deception attacks are time-varying, making it possible to consider general situations involving time-dependent network-induced phenomena and different random phenomena at different sensor nodes. Third, the study considers the infinite-step time correlation of the measurement noises and facilitates the direct estimation of time-correlated additive noise, without relying on the differencing method. Fourth, the covariance-based estimation technique employed does not require knowledge of the evolution model for the signal being estimated. Furthermore, and unlike most previous studies of distributed estimation, which usually obtain optimal linear estimators based on a given structure, we present optimal linear-distributed estimation algorithms under the mean squared error criterion, without requiring a particular structure for the estimators. Another advantage of the proposed distributed estimation scheme is that, while most other approaches incorporate an upper bound for the estimation error covariance, this paper derives an exact expression for the error covariance, which can be calculated offline, regardless of the specific measurement set to be processed. Finally, in comparison with the most closely related authors' previous studies on distributed estimation, the current system model, in contrast to the one in ^[6], incorporates deception attacks and infinite-step time-correlated noise. The system model in ^[27], as the system considered in this work, includes stochastic deception attacks, but, in contrast to the present study, white additive noises are considered. Furthermore, it is worth noting that the derivation of the distributed filtering algorithm in both ^[6] and ^[27] relies on the state-space model equations, whereas the algorithms proposed in this paper do not require explicit information about the state transition equation. Instead, they rely solely on the factorization of the state covariance matrix in a separable form. The superiority of the proposed distributed filter over the one in ^[27] in the presence of infinite-step time-correlated additive noises will be experimentally tested in a numerical simulation example.

The rest of this paper is structured as follows. The networked system model with random parameter matrices, time-correlated noises and deception attacks is presented in Section 2, together with the assumptions required of the stochastic processes involved. The distributed estimation problem is then formulated in Section 3. After this, the distributed estimators are derived in two steps, described in Section 4 and Section 5, respectively. Finally, Section 6 provides an illustrative example highlighting the effectiveness of the proposed distributed estimation algorithms and the main conclusions drawn are summarized in Section 7.

Notation and abbreviations

The mathematical notation and abbreviations used in this paper are detailed in the following table. If not explicitly stated, all vector and matrix dimensions are assumed to be compatible with algebraic operations.

2.   Model description and assumptions

Consider a second-order $n_x$-dimensional discrete-time signal process $\left\{x_k\right\}_{k\geq 1}$, measured by different sensors which are spatially distributed according to a fixed network topology. Specifically, the sensor network is represented by a digraph $\mathcal{G} = (\mathcal{N}, \mathcal{E}, \mathcal{D})$ of order $m$, where $\mathcal{N} = \left\{1, \ldots, m\right\}$ denotes the set of sensor nodes, $\mathcal{E}\subseteq \mathcal{N}\times\mathcal{N}$ is the set of edges connecting different nodes and $(j, i)\in \mathcal{E}$ means that the sensor node $i$ receives the information from node $j$. These link relations among sensors are specified by the adjacency matrix, $\mathcal{D} = \big(d_{ij}\big)_{m\times m}$, with $d_{ij} = 1$ for $(i, j)\in \mathcal{E}$ and $d_{ij} = 0$ otherwise; since any sensor receives its own information, it is clear that $d_{ii} = 1$. For each node $i\in \mathcal{N}$, the set of adjacent nodes plus the node itself is denoted by $\mathcal{N}_{i} = \left\{j\in \mathcal{N}: d_{ji} = 1 \right\}$; therefore $\mathcal{N}_{i}$, which we term the neighborhood of node $i$, is the set of sensor nodes that transmit their information to node $i$.

The network sensors provide noisy measurements of the signal with multiplicative perturbations described by random parameter matrices, and the measurement noises in each sensor are assumed to be sequentially correlated. In our study context, deception attacks may be launched by potential adversaries to replace these measurements by deception noises before they are processed.

We now address the distributed filtering and smoothing problems at each node $i\in {\cal N}$, based on the available signal information in that node, i.e. the measurements derived from all its neighbor nodes $j\in {\cal N}_i$. The estimation is addressed under the LS approach, using only covariance information on the processes involved in the model. Therefore, assumptions (A1)–(A6), listed below, are set to guarantee the existence and knowledge of the first and second-order moments of the signal to be estimated and those of the observations on which the estimation is based. Regarding the signal process, the following assumption –which is key to the recursivity of the estimation algorithms– is required.

(A1) (On the signal). The signal process $\left\{x_k\right\}_{k\geq 1}$ is a zero mean second-order process whose covariance function is expressed in a separable form as follows:

where the factors $\mathcal{A}_k, \mathcal{B}_k$ are $n_x\times N$-dimensional known matrices.

In the following, we describe the measurement model and specify the assumptions made for the random measurement matrices and additive noises.

2.1.   Noisy measurements with random parameter matrices

As previously indicated, at any time $k\geq 1$, the signal $x_k$ is measured by all the network sensors, which provide local outputs perturbed by random parameter matrices and time-correlated additive measurement noises. To describe this situation, consider the following measurement model:

where $\breve{y}_k^{(i)}\in {\Bbb R^{n_{y}}}$ is the measurement of the signal provided by the $i$-th sensor at time $k$, $C_k^{(i)}$ is a random parameter matrix and $v_k^{(i)}$ is the time-correlated measurement noise, which is assumed to be generated by a white process, $\left\{\xi^{(i)}_{k}\right\}_{k\geq 0}$, from an initial noise $v_0^{(i)}$:

where $H^{(i)}_{k}$ are non-singular known deterministic matrices.

To address the estimation problem using covariance information, the following assumptions about the first and second-order moments of the random parameter matrices and the measurement noises are required.

(A2) (On the random parameter matrices). $\left\{C_k^{(i)}\right\}_{k\geq 1}, \ \ i\in {\cal N},$ with $C_k^{(i)} = \Big(c^{(i)}_{_{pq}}(k)\Big)_{n_y\times n_x}$, are independent sequences of independent random parameter matrices whose entries have known first and second-order moments. The mean matrices are denoted by $\overline{C}_k^{(i)} = E\left[C_k^{(i)}\right]$, and their $(p, q)$-entries are given by $E\left[c^{(i)}_{_{pq}}(k)\right]$.

Remark 1. The existence of second-order moments guarantees that of $E\left[C^{(i)}_{k}GC^{(j)T}_{k}\right]$, for any random matrix ${G}$ with mean $\overline{G}$. Moreover, if $G$ is independent of the matrices $C^{(i)}_{k}$ and $C^{(i)}_{k}$, the $(p, q)$-entry of this expectation is given by $\sum_{a = 1}^{n_x}\sum_{b = 1}^{n_x}E\left[c^{(i)}_{_{pa}}(k)c^{(j)}_{_{qb}}(k)\right]\overline{G}_{ab}$.

(A3) (On the measurement noises). The measurement noises $\left\{v_k^{(i)}\right\}_{k\geq 1}, \ i\in {\cal N},$ are time-correlated sequences as described in (2.2), where:

● The initial vectors $v_0^{(i)}, \ i\in {\cal N},$ have zero-mean and known cross-covariance matrices, $\Sigma^{v^{(ij)}}_{0} = E\left[v^{(i)}_0 v^{(j)T}_0\right], \ i, j\in {\cal N}.$

● The white processes $\left\{\xi^{(i)}_{k}\right\}_{k\geq 0}, \ i\in {\cal N},$ are assumed to be independent of each other at different times (consequently, $E\left[\xi^{(i)}_k\xi^{(j)T}_l\right] = 0, \ l\neq k; \ i, j\in {\cal N}$). Their covariance and cross-covariance functions are denoted by $\Sigma^{\xi^{(ij)}}_{k} = E\left[\xi^{(i)}_k\xi^{(j)T}_k\right], \ k\geq 0; \ i, j\in {\cal N}$.

● $\left\{\xi^{(i)}_{k}\right\}_{k\geq 0}, \ i\in{\cal N},$ are independent of the initial vectors $v_0^{(i)}, \ i\in{\cal N}$, and, consequently, $E\left[v^{(i)}_0\xi^{(j)T}_k\right] = 0, \ \ k\geq 0; \ i, j \in{\cal N}.$

Remark 2. As a consequence of (A3), the measurement noises, $\left\{v_k^{(i)}\right\}_{k\geq 1}, \ i\in {\cal N},$ are zero-mean second-order processes with cross-covariance matrices $\Sigma_{k, l}^{v^{(ij)}} = E\left[v^{(i)}_k v^{(j)T}_l\right] = H_{k-1}^{(i)}\cdots H_l^{(i)}\Sigma_{l}^{v^{(ij)}}, \ 1\leq l < k$, and $\Sigma_{k}^{v^{(ij)}}$ is recursively obtained from the relation $\Sigma_k^{v^{(ij)}} = H_{k-1}^{(i)} \Sigma_{k-1}^{v^{(ij)}} H_{k-1}^{(j)T}+ \Sigma^{\xi^{(ij)}}_{k-1}, \ \ k\geq 1,$ with initial condition $\Sigma^{v^{(ij)}}_{0}$.

2.2.   Random deception attack model

In practice, sensor networks are often exposed to attacks from potential adversaries, seeking to modify or deteriorate the real measurements. In this situation, the observations to be processed for estimation may differ from the actual measurements, and so the mathematical model for the measurement outputs after the attacks must be specified.

In this paper, we assume that the measurements are subject to deception attacks that, if successful, will neutralize them and insert deceptive information. Such attacks may or may not be successful, and this uncertainty is incorporated into the model via Bernoulli random variables. Therefore, the mathematical model that describes the potentially attacked measurements to be processed in the signal estimation, which will be denoted by $y^{(i)}_k$, is formulated as:

where $\lambda^{(i)}_k$ are Bernoulli random variables modeling the randomness on the success ($\lambda^{(i)}_k = 1$) or failure ($\lambda^{(i)}_k = 0$) of the attacks, and $\dot{y}^{(i)}_k = -\breve{y}^{(i)}_k + w^{(i)}_k$ is the signal inserted by the attacker, neutralizing the actual measurement, $\breve{y}^{(i)}_k$, and replacing it with a deceptive noise represented by $w^{(i)}_k$.

Therefore, $\lambda^{(i)}_k = 0$ means that the attack against the $i$-th sensor at time $k$ has failed and the measurement to be processed is the actual one ($y^{(i)}_k = \breve{y}^{(i)}_k$), while $\lambda^{(i)}_k = 1$ means that the attack was successful and the processed measurement is the deceptive one ($y^{(i)}_k = w^{(i)}_k$). An equivalent expression for the attacked measurement outputs (2.3) is given by:

The following assumptions are made regarding the processes involved in these observation equations.

(A4) (On the success of attacks). The processes $\left\{\lambda^{(i)}_k\right\}_{k\geq 1}, \ i\in {\cal N},$ are independent sequences of independent Bernoulli random variables with known success probabilities, $P(\lambda^{(i)}_k = 1) = \overline{\lambda}^{(i)}_k$. As a consequence, the first and second-order moments of these variables are:

(A5) (On the deception noises). The noises inserted by successful attacks, $\left\{w^{(i)}_k\right\}_{k\geq 1}, \ i\in {\cal N},$ consist of independent white processes with known covariance and cross-covariance matrices, $\Sigma^{w^{(ij)}}_{k} = E\left[w_k^{(i)}w^{(j)T}_k\right], \ k\geq 1; \ \ i, j\in {\cal N}.$

Remark 3. In this paper, the attack probabilities of success and the covariances and cross-covariances of the noises of the attacks are assumed to be known. If they were unknown, they should be identified before applying the proposed algorithms. To deal with this issue, a distributed self-tuning filtering algorithm is proposed in ^[30], based on the identification of unknown characteristics through the sample zero-order and first-order correlation functions of the observations.

Finally, the following assumption on the processes involved in the described model is required to address the estimation problem.

(A6) (Mutual independence). For each $i\in {\cal N}$, the signal process, $\left\{x_k\right\}_{k\geq 1}$, the random parameter matrices, $\left\{C^{(i)}_k\right\}_{k\geq 1}$, the measurement noise, $\left\{v^{(i)}_k\right\}_{k\geq 1}$, and the processes involved in the attacks, $\left\{\lambda^{(i)}_k\right\}_{k\geq 1}$ and $\left\{w^{(i)}_k\right\}_{ k\geq 1}$, are mutually independent.

3.   Estimation problem formulation

In formulating the estimation problem, our aim is to use the distributed fusion method to obtain estimators of the signal at each sensor node, $i\in {\cal N}$, based on the information available at that node. In accordance with the network topology described in Section 2, this information comes not only from the sensor itself but also from all the neighboring ones that transmit their information to it. Therefore, in the $i$-th sensor, the distributed fusion estimator of the signal $x_k$ based on the information up to time $h$, denoted by $\widehat{x}^{D(i)}_{k/h}$, is obtained by using the local information from the sensor together with that from the nodes in its neighborhood, $\mathcal{N}_{i} = \left\{j\in \mathcal{N}: d_{ji} = 1 \right\}$.

The proposed distributed estimators are derived in two steps for each $i\in \mathcal{N}$: a) derive the intermediate LS linear estimators; b) fuse the neighboring intermediate estimators. In the first step, the LS linear estimator, $\widehat{x}^{(i)}_{k/h}$, is obtained from the potentially attacked observations (2.4) coming from all the sensor nodes $j\in \mathcal{N}_{i}$. Then, in the second step, the distributed estimator $\widehat{x}^{D(i)}_{k/h}$ is determined as the LS matrix-weighted linear combination of the neighboring intermediate estimators, $\widehat{x}^{(j)}_{k/h}, \ j\in \mathcal{N}_{i}$.

In order to unify the derivation of the intermediate estimators in all sensors, we jointly consider all the network output information at each sampling time $k\geq 1$, which is described by the gathered vector $y_k = \left(y_k^{(1)T}, \ldots, y_k^{(m)T}\right)^T$; the outputs corresponding to the sensors in ${\cal N}_i$ are then extracted from this vector to be processed in each sensor $i\in {\cal N}$. Therefore, the intermediate estimator $\widehat{x}^{(i)}_{k/h},$ for each $i\in {\cal N}$, is based on the measurements described by $Y_l^{(i)} = D^{(i)}_yy_l, \ l = 1, \ldots, h$, where $D^{(i)}_y$ is the matrix obtained by removing the all-zero rows of $Diag(d_{1i}, \ldots, d_{mi})\otimes I_{n_y}.$

Next, we specify the gathered measurement model and the statistical properties of the processes involved –derived from our prior assumptions about the local measurements–.

3.1.   Gathered measurement model and properties

To describe the gathered observations to be processed in the estimation, the components of the local observation models in Section 2 are stacked as follows:

From this, the $mn_y$-dimensional observation vector $y_k$, which gathers all the sensor measurements at time $k$, is expressed by:

where $\breve{y}_k$ is the actual measurement vector gathered before the attacks, given by:

The following properties of the processes involved in (3.1)-(3.2) are directly obtained from the corresponding assumptions stated in Section 2.

● $\left\{C_k\right\}_{k\geq 1}$ is a sequence of independent second-order random parameter matrices with means $\overline{C}_k = \left(\overline{C}_k^{(1)T}\big| \cdots \big| \ \overline{C}_k^{(m)T}\right)^T$. Moreover, for any first-order random matrix, ${G}$, $\exists \ E\left[C_{k}GC^{T}_{k}\right] = \left(E\left[C^{(i)}_{k}GC^{(j)T}_{k}\right]\right)_{i, j\in {\cal N}}$, with $(i, j)$-components as defined in Remark 1.

● The measurement noise, $\left\{v_k\right\}_{k\geq 1}$, is a zero-mean second-order time-correlated sequence. From Remark 2, the corresponding covariance matrices, $\Sigma_{k, l}^{v} = E\left[v_{k}v^{T}_{l}\right]$, are given by $\Sigma_{k, l}^{v} = H_{k-1}\cdots H_l\Sigma_{l}^{v}, \ 1\leq l < k$, and $\Sigma_{k}^{v}$ is recursively obtained from the relation $\Sigma_k^{v} = H_{k-1} \Sigma_{k-1}^{v} H_{k-1}^{T}+ \Sigma^{\xi}_{k-1}, \ \ k\geq 1,$ with initial condition $\Sigma^v_{0} = \left(\Sigma_0^{v^{(ij)}}\right)_{i, j\in {\cal N}}$ and $\Sigma^{\xi}_{k} = \left(\Sigma_k^{\xi^{(ij)}}\right)_{i, j\in {\cal N}}$.

Remark 4. The non-singularity of the matrices $H_k$ allows us to factorize the noise covariance matrices in a similar way to those of the signal in (A1); namely:

where $\mathcal{H}_k = H_{k-1} \cdots H_{0}$ and $\mathcal{F}_k = \Sigma_k^v \mathcal{H}_k^{-T}, \ k\geq 1.$

● $\left\{\Lambda_k\right\}_{k\geq 1}$ are diagonal independent matrices with means $\overline{\Lambda}_k = Diag\left(\overline{\lambda}_{k}^{(1)}, \ldots, \overline{\lambda}_{k}^{(m)}\right)\otimes I_{n_y}, \ \ k\geq 1$. For the purpose of further developments, if $G$ is a random matrix independent of $\left\{\Lambda_k\right\}_{k\geq 1}$, these matrices operate as follows:

where $K^{\lambda}_{k} = \left(E\left[\lambda^{(i)}_{k}\lambda_{k}^{(j)}\right]\right)_{i, j\in {\cal N}}\otimes \mathbf{1}_{n_y\times n_y}, \ \ K^{\mathbf{1}-\lambda}_{k} = \left(E\left[\left(1-\lambda^{(i)}_{k}\right)\left(1-\lambda_{k}^{(j)}\right)\right]\right)_{i, j\in {\cal N}}\otimes \mathbf{1}_{n_y\times n_y},$ and the entries $E\left[\lambda^{(i)}_{k}\lambda_{k}^{(j)}\right]$ are given in (A4).

● The deception noise, $\left\{w_k\right\}_{k\geq 1},$ is a white process whose covariance matrices, $\Sigma^{w}_{k} = E\left[w_kw^{T}_k\right]$, are obtained from (A5), $\Sigma^{w}_{k} = \left(\Sigma_k^{w^{(ij)}}\right)_{i, j\in {\cal N}}, \ k\geq 1.$

Finally, assumption (A6) implies the independence of the processes involved in the gathered model (3.1)-(3.2):

● The signal process $\left\{x_k\right\}_{k\geq 1}$, the random parameter matrices, $\left\{C_k\right\}_{k\geq 1}$, the measurement noise, $\left\{v_k\right\}_{k\geq 1}$, and the processes $\left\{\Lambda_k\right\}_{k\geq 1}$ and $\left\{w_k\right\}_{ k\geq 1}$, are all mutually independent.

The above properties guarantee that the observations $\left\{y_k\right\}_{k\geq 1}$ constitute a zero-mean second-order process, whose covariance matrices, $\Sigma^{y}_{k} = E\left[y_ky^{T}_k\right]$, are given by:

3.2.   Innovation approach to the LS estimation

In this paper, the intermediate LS linear estimation problem is addressed by an innovation approach, according to which the observation process at the $i$-th sensor, $\left\{Y_k^{(i)}\right\}_{k\geq 1}$, with $Y_k^{(i)} = D^{(i)}_yy_k$, is replaced by a white process $\left\{\eta_k^{(i)}\right\}_{k\geq 1}$ – innovation process–, such that the LS linear estimator of an arbitrary zero-mean second-order vector $a_k$ based on the observations $\left\{Y_1^{(i)}, \ldots, Y_h^{(i)}\right\}$ can be obtained by the following linear combination of the corresponding innovations, $\left\{\eta_1^{(i)}, \ldots, \eta_h^{(i)}\right\}$:

where $\Sigma^{\eta^{(i)}}_k = E\left[\eta_{k}^{(i)}\eta_{k}^{(i)T}\right]$ denotes the innovation covariance.

The innovation at time $k$ is defined as $\eta_{k}^{(i)} = Y_{k}^{(i)}-\widehat{Y}_{k/k-1}^{(i)},$ where $\widehat{Y}_{k/k-1}^{(i)}$ is the LS linear one-stage predictor of the observation $Y_{k}^{(i)} = D^{(i)}_yy_k$, which is clearly given by $\widehat{Y}_{k/k-1}^{(i)} = D^{(i)}_y\widehat{y}_{k/k-1}^{(i)}$. Hence, the innovations are written as:

where, according to (3.5), expressions (3.1)-(3.2) and the independence properties on the model, the predictor of $y_k$ is expressed in terms of those for the signal and noise as follows:

Expressions (3.6) and (3.7) provide the basis for obtaining the intermediate estimation algorithms presented in the following section.

4.   Intermediate LS linear estimation

In this section, the LS linear filtering and fixed-point smoothing estimators of the signal $x_k$ are obtained in each sensor according to the observations available at that sensor. More precisely, for each $i\in {\cal N},$ we first derive a recursive algorithm for the intermediate LS linear filter, $\widehat{x}^{(i)}_{k/k}$, based on the observations $Y_l^{(i)} = D^{(i)}_yy_l, \ l = 1, \ldots, k.$ The filtering estimator of $x_k$ is then recursively updated with successive observations, $Y_l^{(i)} = D^{(i)}_yy_l, \ l = k+1, \ldots, h,$ in order to obtain the intermediate LS linear smoothing estimators, $\widehat{x}^{(i)}_{k/h}, \ h > k$.

The intermediate filtering and fixed-point smoothing algorithms are presented in Theorems 1 and 2, respectively. Both theorems also include recursive formulas for the error covariance matrices $\Sigma_{k/h}^{\widetilde{x}^{(i)}} = E\left[\widetilde{x}^{(i)}_{k/h}\widetilde{x}^{(i)T}_{k/h}\right]$, with $\widetilde{x}^{(i)}_{k/h} = x_k-\widehat{x}_{k/h}^{(i)}$, thus providing the natural measure of the estimation accuracy under the LS criterion.

Theorem 1. For any $i\in {\cal N}$, the intermediate LS linear filtering estimators, $\widehat{x}^{(i)}_{k/k}$, and the error covariance matrices, $\Sigma_{k/k}^{\widetilde{x}^{(i)}}$, are calculated as:

The vectors $\mathbf{z}^{(i)}_k$ and their covariance matrices, $\Sigma^{\mathbf{z}^{(i)}}_k$, satisfy the following recursive formulas:

where the gain coefficients, ${\cal Z}^{(i)}_k$, are:

and the innovations, $\eta^{(i)}_k$, and their covariance matrices, $\Sigma^{\eta^{(i)}}_k$, are given by:

The matrices $\mathcal{H}_k$ and $\mathcal{F}_k$ are the noise covariance factors defined in (3.3) and the observation covariance, $\Sigma^{y}_{k}$, is given in (3.4).

Proof. See Appendix A.

The key of the recursive algorithm in Theorem 1 is the intermediate filter factorization (4.1), in which the vectors $\mathbf{z}_k^{(i)}$ can be recursively obtained as indicated in (4.3). Equation (4.2) provides the filtering error covariance matrices, which are needed to evaluate the performance of the filters; since these covariances do not depend of the observations, the filter performance can be measured before the observation set is available. In view of (4.1) and (4.2), the problem is focused on obtaining the vectors $\mathbf{z}_k^{(i)}$ and their covariances, $\Sigma_k^{\mathbf{z}^{(i)}}$, which are derived from the innovations and their covariance matrices by the recursive equations (4.3) and (4.4), respectively. Specifically, by starting from $\mathbf{z}_{k-1}^{(i)}$, $\Sigma_{k-1}^{\mathbf{z}^{(i)}}$ and the new observation $y_k^{(i)}$, the gain, innovation and innovation covariance are calculated by equations (4.5), (4.6) and (4.7), respectively; then, equations (4.3) and (4.4) are used for obtaining $\mathbf{z}_k^{(i)}$ and $\Sigma_k^{\mathbf{z}^{(i)}}$.

Theorem 2. For any $i\in {\cal N}$, the intermediate LS linear smoothing estimators, $\widehat{x}^{(i)}_{k/h}, \ h > k$, and the error covariance matrices, $\Sigma_{k/h}^{\widetilde{x}^{(i)}}$, are recursively calculated from the filter, $\widehat{x}^{(i)}_{k/k}$, and the error covariance $\Sigma_{k/k}^{\widetilde{x}^{(i)}}$, respectively, given in Theorem 1:

where:

and

Proof. See Appendix B.

By starting from the filter, equation (4.8) explains how the estimators of $x_k$ are updated when successive observations are available. Besides the innovations –which are given in Theorem 1–, for the update it is necessary to know the gain coefficients ${\cal X}^{(i)}_{k, h}$, which are obtained from equation (4.10) and require, in turn, the cross-covariance matrices $\Sigma^{x\mathbf{z}^{(i)}}_{k, h}$, that are given by equation (4.11). Finally, equation (4.9) provides the smoothing error covariances in a recursive way.

Remark 5. By replacing $D^{(i)}_y$ by the matrix obtained taking $d_{ji} = 0, \ j\neq i$, the LS local linear filter and smoothers, based only on the measurements from the $i$-th sensor, can be obtained from the formulas in Theorems 1 and 2, respectively. Furthermore, the hypothetical global optimal linear estimators –based on measurements from all sensors– can be calculated by replacing $D^{(i)}_y$ by the identity matrix.

5.   Distributed fusion estimation

The second step of the distributed fusion method is performed in order to derive the distributed filter and fixed-point smoother of the signal, $x_k$. As indicated above, for each $i\in \mathcal{N}$, the distributed estimators, $\widehat{x}^{D(i)}_{k/h}, \ h\geq k,$ are derived as the matrix-weighted linear combination of the neighboring intermediate estimators, $\widehat{x}^{(j)}_{k/h}, \ j\in \mathcal{N}_i,$ using the mean squared error as an optimality criterion.

Specifically, the distributed estimators are expressed as $\widehat{x}^{D(i)}_{k/h} = F^{(i)}_{k/h}\widehat{X}^{(i)}_{k/h}, \ h\geq k,$ for each $i\in \mathcal{N},$ where $\widehat{X}^{(i)}_{k/h}$ is the vector obtained by assembling all the intermediate estimators of the sensors connected to the $i$-th one, and where $F^{(i)}_{k/h}$ is a matrix to be determined in order to minimize the mean squared error, $E\left[\left(x_k-F_{k/h}^{(i)}\widehat X^{(i)}_{k/h}\right)^T \left(x_k-F_{k/h}^{(i)}\widehat X^{(i)}_{k/h}\right)\right].$

As it is known, the solution to this problem is given by:

and, consequently, both matrices, $E\left[x_{k}\widehat X^{(i)T}_{k/h}\right]$ and $E\left[\widehat X^{(i)}_{k/h}\widehat X^{(i)T}_{k/h}\right]$, must be determined for each $i\in {\cal N}$.

The entries of the above matrices are extracted from $E\left[x_{k}\widehat X^{T}_{k/h}\right]$ and $E\left[\widehat X_{k/h}\widehat X^{T}_{k/h}\right]$, respectively, where $\widehat{X}_{k/h} = \left(\widehat{x}^{(1)T}_{k/h}, \dots, \widehat{x}^{(m)T}_{k/ h}\right)^T$ denotes the vector stacking all the intermediate estimators in the network. For this purpose, we write $\widehat{X}^{(i)}_{k/h} = D^{(i)}_x\widehat{X}_{k/h}$, where $D^{(i)}_x$ is the matrix obtained by removing the all-zero rows of $Diag(d_{1i}, \dots, d_{mi})\otimes I_{n_x}$. The optimal distributed filter and the smoothers of the signal $x_k$ are then given by:

Hence, both matrices, $E\left[x_{k}\widehat X^{T}_{k/h}\right]$ and $E\left[\widehat X_{k/h}\widehat X^{T}_{k/h}\right]$ must be determined. Taking into account the OPL, each intermediate estimator $\widehat{x}^{(r)}_{k/h}$ is uncorrelated with the estimation error; therefore, the entries of $E\left[x_{k}\widehat X^{T}_{k/h}\right]$ can be rewritten as $E\left[x_k\widehat{x}^{(r)T}_{k/h}\right] = E\left[\widehat{x}^{(r)T}_{k/h}\widehat{x}^{(r)T}_{k/h}\right], \ r\in{\cal N}$ and then, denoting the cross-covariance between intermediate estimators as $\Sigma^{\widehat{x}^{(rs)}}_{k/h} = E\left[\widehat x^{(r)}_{k/h}\widehat x^{(s)T}_{k/h} \right], \ r, s\in{\cal N}$, the matrices involved in (5.1) are given by $E\left[\widehat X_{k/h}\widehat X^{T}_{k/h}\right] = \left(\Sigma^{\widehat{x}^{(rs)}}_{k/h}\right)_{r, s\in {\cal N}}$ and $E\left[x_{k}\widehat X^{T}_{k/h}\right] = \left(\Sigma^{\widehat{x}^{(1)}}_{k/h}\ \big| \cdots \big| \ \Sigma^{\widehat{x}^{(m)}}_{k/h} \right).$

The algorithms needed to obtain the cross-covariances between the intermediate filters and smoothers, $\Sigma^{\widehat{x}^{(rs)}}_{k/h}, \ h\geq k, \ r, s\in{\cal N},$ are stated in the following lemmas.

Lemma 1 presents the expression of the cross-covariances between the filters, $\Sigma^{\widehat{x}^{(rs)}}_{k/k}, \ r, s\in {\cal N},$ which is directly obtained from (4.1). Thus, it depends on the cross-covariances $\Sigma^{\mathbf{z}^{(rs)}}_k = E\left[\mathbf{z}^{(r)}_{k}\mathbf{z}^{(s)T}_{k}\right]$; a recursive algorithm is then stated for $\Sigma^{\mathbf{z}^{(rs)}}_k$, which in turn depends on the innovation cross-covariance matrices, $\Sigma^{\eta^{(rs)}}_{k} = E\left[\eta_{k}^{(r)}\eta_{k}^{(s)T}\right]$.

Lemma 1. For $r, s \in \mathcal{N},$ the cross-covariance matrices between the intermediate filters, $\Sigma^{\widehat{x}^{(rs)}}_{k/k} = E\left[\widehat{x}^{(r)}_{k/k}\widehat{x}^{(s)T}_{k/k}\right],$ are computed by:

where $\Sigma^{\mathbf{z}^{(rs)}}_k$ satisfies the following recursive relation:

The matrices $\mathcal{Z}^{(rs)}_{k-1, k} = E\left[\mathbf{z}^{(r)}_{k-1}\eta_k^{(s)T}\right]$ are given by:

and the innovation cross-covariance matrices, $\Sigma^{\eta^{(rs)}}_{k}$, satisfy:

Proof. See Appendix C.

In the second lemma, for any $r, s\in {\cal N}$, the cross-covariances between the intermediate smoothers, $\Sigma^{\widehat{x}^{(rs)}}_{k/h}, \ h > k,$ are obtained from that of the filter, $\Sigma^{\widehat{x}^{(rs)}}_{k/k}$, by a recursive algorithm.

Lemma 2. For $r, s \in \mathcal{N},$ the cross-covariance matrices between the intermediate fixed-point smoothers, $\Sigma^{\widehat{x}^{(rs)}}_{k/h} = E\left[\widehat{x}^{(r)}_{k/h}\widehat{x}^{(s)T}_{k/h}\right], \ h > k,$ are recursively obtained from:

with initial condition $\Sigma^{\widehat{x}^{(rs)}}_{k/k}$, given in Lemma 1.

The matrices $\widehat{{\cal X}}^{(rs)}_{k, h} = E\left[\widehat{x}^{(r)}_{k/h-1}\eta_{h}^{(s)T}\right]$ are derived as:

where $\Sigma^{x\mathbf{z}^{(r)}}_{k, h}$ is given in Theorem 2 and $\Sigma^{\widehat{x}\mathbf{z}^{(rs)}}_{k, h}$ is recursively computed by:

with initial condition $\Sigma^{\widehat{x}\mathbf{z}^{(rs)}}_{k} = (\mathcal{A}_k \ |\ 0)\Sigma_{k}^{\mathbf{z}^{(rs)}}, \ k\geq 1$.

Proof. See Appendix D.

The above results for the distributed filters and fixed-point smoothers are presented in the following theorem, which also includes the formula to calculate their estimation error covariance matrices, $\Sigma_{k/h}^{D(i)} = E\left[\left(x_k-\widehat{x}^{D(i)}_{k/h}\right)\left(x_k-\widehat{x}^{D(i)}_{k/h}\right)^T \right]$.

Theorem 3. For each $i\in \mathcal{N},$ the distributed filter and fixed-point smoothers, $\widehat{x}^{D(i)}_{k/h}, \ h\geq k$, are given by:

where $\widehat{X}_{k/h}$ is the stacked vector of intermediate estimators and $\Sigma^{\widehat{x}^{(rs)}}_{k/h}$ are the cross-covariances between them, given in Lemma 1 and Lemma 2 for the filter and smoothers, respectively. The estimation error covariance matrices, $\Sigma_{k/h}^{D(i)}$, are given by:

Proof. Expression (5.9) is merely a rewrite of (5.1). To derive (5.10), the OPL must be used to express $\Sigma_{k/h}^{D(i)} = E\left[x_kx_k^T \right]-E\left[\widehat{x}^{D(i)}_{k/h}\widehat{x}^{D(i)T}_{k/h} \right]$, and the formula is straightforwardly achieved by using the signal covariance factorization, (A1), together with (5.9). □

6.   Numerical simulation example

The applicability and effectiveness of the proposed distributed estimators are illustrated by the sensor network shown in Figure 1, where the topology is represented by a digraph $\mathcal{G} = (\mathcal{N}, \mathcal{E}, \mathcal{D})$, with the set of nodes $\mathcal{N} = \left\{1, 2, 3, 4, 5\right\}$. The elements in the adjacency matrix, $\mathcal{D}$, are $d_{ji} = 1$ when sensor $i$ can obtain information from sensor $j$; otherwise, $d_{ji} = 0$:

6.1.   Signal process and measurement model

6.1.1.   Signal process

Consider a two-dimensional stochastic signal, $x_k\hskip -0.1cm = \hskip -0.1cm \left(\hskip -0.1cm \begin{array}{c}x_{1k}\\ x_{2k}\end{array}\hskip -0.1cm \right)$, whose first and second component, $x_{1k}$ and $x_{2k}$, denote the position and velocity of a target, respectively. Moreover, assume that the signal evolution is the same as in ^[30]; namely:

where

The initial signal, $x_0$, is a two-dimensional standard Gaussian vector. The multiplicative noises $\big\{\alpha^n_{k}\big\}_{k\geq 0}$, $\big\{\beta^n_{k}\big\}_{k\geq 0}$, $n = 1, 2$, are white Gaussian scalar processes, with variances $Var\big(\alpha^n_{k}\big) = 0.16$, $Var\big(\beta^n_{k}\big) = 0.11$, $n = 1, 2$. The additive noise $\big\{\mu_k\big\}_{k\geq 0}$ is a zero-mean Gaussian scalar process with variance 0.5. All of these noise sequences and the initial signal are assumed to be mutually independent. The signal covariance function can then be expressed in a separable form as $E\Big[x_kx_l^T\Big] = \mathcal{A}_k \mathcal{B}^T_l$, $l\leq k$, with $\mathcal{A}_k = A^k$ and $\mathcal{B}^T_l = A^{-l}\Sigma^x_l$, where $\Sigma^x_l = E\Big[x_lx_l^T\Big]$ is recursively obtained by:

6.1.2.   Sensor measurements with multiplicative perturbations and time-correlated noises

Scalar measurements of the signal, simultaneously disturbed by fading effects and both additive and multiplicative noises, are provided by the five sensor nodes, according to model (2.1), where the processes involved have the following characteristics:

● The multiplicative perturbations are described by $C^{(i)}_k = \theta^{(i)}_k\Big(\dot{C}^{(i)}+\rho^{(i)}_k\ddot{C}^{(i)}\Big), \ i\in \mathcal{N}$, where:

- $\dot{C}^{(1)} = (0.8, \ 0.9)$, $\dot{C}^{(2)} = (0.9, \ 0.7)$, $\dot{C}^{(3)} = (0.6, \ 0.7)$, $\dot{C}^{(4)} = (0.7, \ 0.8)$, $\dot{C}^{(5)} = (0.9, \ 0.5)$, $\ddot{C}^{(i)} = (1, \ 0), \ i = 1, 2,$ and $\ddot{C}^{(i)} = (0, \ 1), \ i = 3, 4, 5$.

- $\big\{\theta^{(i)}_k\big\}_{k\geq 1}$ and $\big\{\rho^{(i)}_k\big\}_{k\geq 1}$, $i\in \mathcal{N}$, are independent sequences of independent and identically distributed random variables; namely:

* $\theta^{(1)}_k$ and $\theta^{(2)}_k$ are random variables uniformly distributed over $[0.3, 0.7]$ and $[0.2, 0.8]$, respectively, representing a continuous random fading effect on the measurements from sensors 1 and 2.

* $\theta^{(3)}_k$ and $\theta^{(4)}_k$ are discrete random variables, representing a discrete random fading effect on the measurements from sensors 3 and 4, according to the following probability mass functions:

* $\theta^{(5)}_k$ are Bernoulli random variables, where $P\big(\theta^{(5)}_k = 1\big) = \overline{\theta}$, represents the randomly missing measurement phenomenon in sensor 5.

* The multiplicative components $\rho^{(i)}_k$, $i\in \mathcal{N},$ are standard Gaussian variables.

● The time-correlated noises $\big\{v_{k}^{(i)}\big\}_{k\geq 1}$, $i\in \mathcal{N}$, are defined by (2.2), with $H^{(1)} = H^{(3)} = H^{(5)} = 0.8$ and $H^{(2)} = H^{(4)} = 0.7$. These noises are generated by the processes $\big\{\xi^{(i)}_{k}\big\}_{k\geq 0}$, defined by $\xi_{k}^{(i)} = a^{(i)}\xi_{k}$, $i\in \mathcal{N}$, where $\big\{\xi_{k}\big\}_{k\geq 0}$ is a standard Gaussian white process, $a^{(1)} = a^{(3)} = a^{(5)} = 0.5$, $a^{(2)} = a^{(4)} = 0.25$, and the initial conditions are $v_{0}^{(i)} = v_0$, for $i\in \mathcal{N}$, and where $v_0$ is a standard Gaussian variable.

Under these assumptions, the covariance matrices of the stacked noise $v_k = \left(v^{(1)}_k, \ldots, v^{(5)}_k\right)^T$ can always be factorized as $\Sigma_{k, l}^v = \mathcal{H}_k \mathcal{F}_l^T, \ l\leq k,$ where $\mathcal{H}_k = H^k$ and $\mathcal{F}^T_l = H^{-l}\Sigma^v_l$, being $H = Diag(0.8, 0.7, 0.8, 0.7, 0.8)$ and where $\Sigma^v_l$ is recursively obtained by $\Sigma_l^{v} = H \Sigma_{l-1}^{v} H^{T}+\Sigma^{\xi}_{l-1}, \ \ l\geq 1,$ with $\Sigma^{\xi}_{l} = \left(a^{(i)}a^{(j)}\right)_{i, j\in{\cal N}}$ and initial condition $\Sigma^v_{0} = \mathbf{1}_{5\times 5}$.

6.1.3.   Sensor measurements subject to deception attacks

In accordance with the theoretical model, let us suppose that the measurements at each sensor are subject to deception attacks and that the attacked measurement outputs are given by (2.4), where:

● The noises of the false data injection attacks are defined as $w_{k}^{(i)} = b^{(i)}w_{k}$, for $i\in \mathcal{N}$, where $b^{(1)} = b^{(2)} = b^{(5)} = 0.5$, $b^{(3)} = b^{(4)} = 0.75$, and $\big\{w_{k}\big\}_{ k\geq 1}$ is a standard Gaussian white process. Clearly, these attack noises are correlated and $\Sigma^{w^{(ij)}}_{k} = b^{(i)}b^{(j)}, \ i, j \in \mathcal{N}.$

● The status of the attacks is described by mutually independent sequences of independent and identically distributed Bernoulli random variables, $\big\{\lambda^{(i)}_k\big\}_{ k\geq1}$, $i\in \mathcal{N}$, with known probabilities $P\big(\lambda^{(i)}_k = 1\big) = \overline{\lambda}.$

6.2.   Applicability of the proposed distributed algorithms

To illustrate the effectiveness of the proposed distributed filtering and fixed-point smoothing algorithms and to quantify the estimation accuracy obtained, the estimation error variances of the first and second signal components (position and velocity) were calculated at every sensor node $i\in \mathcal{N}$. First, the local, intermediate, distributed and global estimators were compared, for fixed values of the probabilities $\overline{\theta}$ (probability that the signal is present in the measured outputs of sensor 5) and $\overline{\lambda}$ (probability of a successful attack). Different values of the probabilities $\overline{\theta}$ and $\overline{\lambda}$ were then considered to highlight the effects of the missing measurement and attack phenomena, respectively, on the performance of the proposed distributed estimators, analyzing how these probabilities influence the estimation error variances for both signal components.

6.2.1.   Comparison of local, intermediate and distributed estimators

Considering the same value of 0.5 for the probabilities $\overline{\theta}$ and $\overline{\lambda}$, Figure 2 depicts –for the first signal component– the error variances of the local filters (obtained using only the measurements from the sensor itself) and those of the proposed intermediate filters, $\widehat{x}^{(i)}_{k/k}$, and distributed filters, $\widehat{x}^{D(i)}_{k/k}$, and smoothers, $\widehat{x}^{D(i)}_{k/k+L}$ (with lag $L = 1, 2, 3, 4, 5$), at every sensor node $i\in \mathcal{N}.$ On the one hand, this figure shows that the error variances corresponding to the intermediate filters are significantly less than those of the local filters and that the distributed filters outperform the intermediate ones. In addition, it is apparent that the distributed smoothing error variances are less than the filtering ones and, also, that at each fixed-point $k$, the fixed-point smoothers become more accurate as the number of available observations, $k+L$, increases. As expected, the improvement is smaller as $L$ increases; indeed, in this example, the improvement is practically imperceptible for $L\geq 5$. Similar results are obtained for the second signal component.

6.2.2.   Comparison of intermediate, distributed and global optimal linear estimators

Close inspection of Figure 2 reveals little difference in the values of the distributed estimation error variances over the five nodes. For a better evaluation of these differences, Figure 3 shows, for the first and second signal components, the error variances of the intermediate filters together with those of the distributed filters and smoothers (for $L = 1, 2$) at the different nodes, as well as those of the global optimal linear filtering and smoothing estimators, based on the set of measurements obtained from the five nodes of the network. A desirable property for distributed estimators over sensor networks is that the discrepancies between different nodes should be as small as possible; indeed, as we can see in Figure 3, the proposed distributed estimators considerably reduce the disagreements among intermediate estimators from different sensors. Thus, the distance between the error variances of the distributed estimators at the different nodes is fairly small and close to the global optimal error variances. In consequence, not only do the proposed distributed estimators present only slight discrepancies among the sensors, but they also provide a very similar level of performance to that of the global optimal estimators. Moreover, the proximity between the error variances of the global optimal filtering and smoothing estimators and those of the corresponding proposed distributed estimators shows that the latter estimators perform well.

6.2.3.   Effect of the probability of missing measurements on estimation accuracy

Assuming, as above, that the attack probability is $\overline{\lambda} = 0.5$, we now evaluate the performance of the proposed distributed estimators with respect to the missing measurements phenomenon in sensor node 5. To do so, the distributed filtering and smoothing ($L = 2$) error variances for the first signal component are plotted in Figure 4 for different values of the probability $\overline{\theta}$ (namely, $\overline{\theta} = 0.3, 0.5, 0.7$ and $0.9$). For these values, the distributed estimation error variances at sensor node 5 are plotted in the left-side panel of Figure 4. This figure shows that the performance of the distributed estimators is indeed influenced by these probabilities and, as expected, that the distributed estimation error variances decrease as the probability $\overline{\theta}$ increases. Hence, both the filtering and the smoothing distributed estimators achieve better estimation accuracy when $1-\overline{\theta}$, the probability of missing measurements, decreases, as this means that further information about the signal is available. Analogous results were obtained for sensor nodes 2, 3 and 4, since these nodes all utilize information from node 5; in fact, node 4 uses the measurements from node 5 to obtain its intermediate estimators, which are subsequently sent to nodes 2 and 3 to construct the distributed estimators at these nodes. In node 1, however, the distributed error variances do not change since the estimators in this node do not use the measurements from node 5. Since the behavior of the distributed error variances is analogous in all the iterations, for a better visualization of the decreasing trend of the error variances as the probability $\overline{\theta}$ increases at all sensor nodes –except sensor 1, in which it remains constant–, the right-side panel of Figure 4 displays the distributed error variance at the iteration $k = 100$. As in the previous figures, Figure 4 also shows that the error variances corresponding to the smoothers are less than those of the filters. Similar results are obtained for the second signal component and therefore the same conclusions are drawn.

6.2.4.   Effect of the probability of a successful attack on estimation accuracy

Next, our aim is to examine the association between the probability of a successful attack and the performance of the estimators. For this purpose, we compare the distributed filtering error variances for different values of this probability (namely, $\overline{\lambda} = 0.1$ to 0.9). The distributed filtering error variances at sensor node 1, corresponding to the second signal component, are plotted in the left-side panel of Figure 5. Here, as expected, these error variances rise in line with $\overline{\lambda}$. Furthermore, this increase is more pronounced for higher values of $\overline{\lambda}$. Similar results are obtained in all nodes, as shown in the right-side panel of Figure 5, which displays the distributed filtering and smoothing error variances at $k = 100$ versus $\overline{\lambda}$ in the five sensor nodes. The discrepancies between the different nodes are negligible and, as shown in all the other figures, the smoother with lag $L = 2$ outperforms the one with $L = 1$ which, in turn, outperforms the filter. Similar results are inferred for the first signal component, and therefore the same conclusions are drawn.

6.2.5.   Comparative analysis with the distributed filter ^[27]

Our final aim in this section is to show the superior performance of the proposed distributed filter in the presence of infinite-step time-correlated additive noises. For this purpose, we conduct a comparative analysis between the distributed filter proposed in this paper and the one proposed in ^[27] for networked systems with fading measurements, multiplicative noises in both the signal and measurements equations and stochastic deception attacks, but without infinite-step time correlation of the measurement noises.

At every sensor node and for each of the two components of the distributed filtering estimates, the comparison is made on the basis of the empirical values of the mean-squared error at each time instant, which are calculated from two thousand independent simulations by

where, for each sampling time $k$ and for the $s$-th simulation run, $x_{a, k}^{(s)}$ denotes the $a$-th component of the simulated signal, and $\widehat{x}_{a, k/k}^{D(i, s)}$ is the $a$-th component of the distributed filter calculated in the $i$-th sensor node.

Assuming again the same fixed value 0.5 for the probabilities $\overline{\theta}$ and $\overline{\lambda}$, the results are displayed in Figure 6, which shows that for all the sensor nodes and for both the first and second components, the empirical mean-squared error values of the proposed distributed filtering estimates are lower than those of the distributed filtering estimates in ^[27]. Note that the proposed filter was indeed expected to outperform the filter in ^[27], since the latter does not take into account the infinite-step time correlation of the sensor measurement noises.

7.   Conclusions

In this paper, we investigate the distributed estimation problem –including filtering and fixed-point smoothing– in networked systems whose sensor nodes are spatially distributed according to a predetermined network topology, represented by a directed graph. Random parameter matrices and stochastic deception attacks are incorporated into the measurement model. Thus, a broad theoretical framework is provided with which to address general stochastic multi-sensor systems with different network-induced uncertainties. The presence of time-correlated additive noise in the observation model is handled by a non-augmentation method, based on the direct estimation of the noise. For every sensor node, the proposed distributed estimation algorithm runs in two phases. The first yields an intermediate least-squares linear estimator using its own local measurements and those received from its neighboring nodes. In the second phase, the own-sensor intermediate estimator is combined with those calculated by its neighboring nodes to obtain the desired distributed estimator as the minimum mean squared error matrix-weighted linear combination of the intermediate estimators. The proposed estimation strategy does not rely on explicit information provided by the signal evolution equation, but rather on the factorization of the signal and time-correlated noise covariance matrices in a separable form. As a result, whether or not the signal evolution model is completely known, the proposed distributed estimation technique can be used to estimate a wide class of stochastic signals, including those whose evolution is affected by multiplicative noises. The simulation experiment performed shows that the theoretical system model we present covers common random imperfections, such as the presence of multiplicative noise, missing observations and fading effects. The numerical results obtained were used to examine the influence of two degrading effects on estimation performance: a) the probability of missing measurements; b) the probability of successful attacks. Comparative analysis of the estimation error variances shows that the proposed distributed estimators outperform the intermediate ones and reduce the disagreements between different sensors by bringing each distributed estimator closer to the global optimal linear estimator, based on the full set of measurements of the entire network. Finally, in the presence of infinite-step colored measurement noises, the proposed estimators are shown to outperform the distributed estimators in the authors' previous work ^[27].

A challenging topic for future studies is the derivation of self-tuning estimation algorithms for the case when the attack probabilities and/or the covariances and cross-covariances of the attack noises are unknown. It would also be interesting to consider the possibility that attack noise is not stochastic but a constant or time-varying deterministic sequence, as well as the scenario of random packet dropouts in the transmissions among sensor nodes.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This research was suported by the "Ministerio de Ciencia e Innovación, Agencia Estatal de Investigación" of Spain and the European Regional Development Fund [grant number PID2021-124486NB-I00].

Conflict of interest

The authors declare there is no conflict of interest.

Appendix A. Proof of Theorem 1

According to expression (3.5) for the LS linear estimators, the coefficients ${{\cal X}}^{(i)}_{k, l} = E\left[x_k\eta^{(i)T}_l\right], \ 1\leq l\leq k,$ with $\eta^{(i)}_l$ given from (3.6) and (3.7), must be calculated for the signal filter $\widehat{x}^{(i)}_{k/k}$. For this purpose, we use expressions (3.1) and (3.2) for $y_l$ which, together with (A1) on the signal covariance and the independence properties, easily lead to:

Now, by denoting ${\mathcal V}^{(i)}_{k, l} = E\left[v_k\eta^{(i)T}_{l}\right]$ and using (3.5) for the one-stage predictors $\widehat{x}^{(i)T}_{l/l-1}$ and $\widehat{v}^{(i)T}_{l/l-1}$, we obtain:

and ${{\cal X}}^{(i)}_{k, l}$ is then written as:

In a similar way, using now the factorization (3.3) for the noise covariance, the coefficients ${\mathcal V}^{(i)}_{k, l}$ are also written as:

These expressions guarantee the following factorizations:

where

and, hence, the augmented matrices ${\cal Z}^{(i)}_k = \Big({\cal Z}^{x(i)T}_k \big| \ {\cal Z}^{v(i)T}_k\big)^T$ satisfy:

From these preliminaries, we can now derive the expressions of the theorem at hand. Expression (4.1) for the signal filter comes directly from (3.5) using the factorization (7.1) for ${\mathcal X}^{(i)}_{k, l}$, and defining $\mathbf{z}^{(i)}_k\equiv {\sum_{l = 1}^{k}}{\cal Z}^{(i)}_l\Big(\Sigma_l^{\eta^{(i)}}\Big)^{-1}\eta^{(i)}_l$ which clearly satisfies (4.3). Now, using the OPL –which guarantees the uncorrelation between the estimation error and the observations–, we express $\Sigma_{k/k}^{\widetilde{x}^{(i)}} = E\left[x_kx_k^T\right]- E\left[\widehat{x}^{(i)}_{k/k}\widehat{x}^{(i)T}_{k/k}\right]$, and (4.2) is easily obtained from (4.1), by denoting $\Sigma^{\mathbf{z}^{(i)}}_k = E\left[\mathbf{z}^{(i)}_k\mathbf{z}_k^{(i)T}\right]$. Recursive formula (4.4) for these matrices is derived by from the fact that the innovation process is white, which leads to $\Sigma^{\mathbf{z}^{(i)}}_k = {\sum_{l = 1}^{k}}{\cal Z}^{(i)}_l\Big(\Sigma_l^{\eta^{(i)}}\Big)^{-1}{\cal Z}^{(i)T}_l$; now, (4.5) is merely expression (7.2) for ${\cal Z}^{(i)}_k$ rewritten in terms of $\Sigma^{\mathbf{z}^{(i)}}_{k-1}$.

Again, taking again into account the factorizations (7.1), the one-stage predictors of the signal and noise are given by $\widehat{x}^{(i)}_{k/k-1} = ({\cal A}_k \ | \ 0)\mathbf{z}^{(i)}_{k-1}$ and $\widehat{v}^{(i)}_{k/k-1} = (0 \ | \ \mathcal{H}_k)\mathbf{z}^{(i)}_{k-1}$, respectively, and so (3.7) can be rewritten as:

Expression (4.6) for the innovation is obtained immediately from (3.6), and (4.7) is deduced by again using the OPL to express the innovation covariance matrix as $\Sigma^{\eta^{(i)}}_k = D^{(i)}_y \left(\Sigma_k^{y}-E\left[\widehat {y}^{(i)}_{k/k-1} \widehat {y}^{(i)T}_{k/k-1}\right]\right)D^{(i)T}_y.$ □

Appendix B. Proof of Theorem 2

Recursive formula (4.8) is an immediate consequence of the general expression (3.5). From this, the estimation error is expressed as $\widetilde{x}^{(i)}_{k/h} = \widetilde{x}^{(i)}_{k/h-1}-{\cal X}^{(i)}_{k, h}\Big(\Sigma_h^{\eta^{(i)}}\Big)^{-1}\eta^{(i)}_{h}$, and (4.9) is derived straightforwardly taking into account that, from the OPL, $E\left[\widetilde{x}^{(i)}_{k/h-1}\eta^{(i)T}_{h}\right] = E\left[x_k\eta^{(i)T}_{h}\right] = {\cal X}^{(i)}_{k, h}$.

Since, for $h > k$, $E\left[x_ky^{T}_{h}\right] = {\cal B}_k{\cal A}_h^{T}\overline{C}_h^T\left(I_{mn_y}-\overline{\Lambda}_h\right)$, expression (4.10) for the gain coefficients, ${\cal X}^{(i)}_{k, h} = E\left[x_k\eta^{(i)T}_{h}\right]$, is obtained using (4.6) for $\eta^{(i)}_{h}$ and denoting $\Sigma^{x\mathbf{z}^{(i)}}_{k, h} = E\left[x_k\mathbf{z}^{(i)T}_{h}\right], \ h\geq k$.

Finally, recursive formula (4.11) is deduced from (4.3) for $\mathbf{z}^{(i)}_{h}$; its initial condition is another consequence of the OPL, from which $E\left[x_k\mathbf{z}^{(i)T}_{k}\right] = E\left[\widehat{x}^{(i)}_{k/k}\mathbf{z}^{(i)T}_{k}\right]$, and it is obtained simply by using (4.1) for $\widehat{x}^{(i)}_{k/k}$.

□

Appendix C. Proof of Lemma 1

Clearly, expression (5.2) for the cross-covariance matrices is obtained from (4.1) for $\widehat{x}^{(r)}_{k/k}$ and $\widehat{x}^{(s)}_{k/k}$, being $\Sigma^{\mathbf{z}^{(rs)}}_k = E\left[\mathbf{z}^{(r)}_{k}\mathbf{z}_k^{(s)T}\right]$. Expression (5.3) for $\Sigma^{\mathbf{z}^{(rs)}}_k$ is derived by using (4.3) for $\mathbf{z}^{(r)}_k$ and $\mathbf{z}^{(s)}_k,$ and denoting $\mathcal{Z}^{(rs)}_{k-1, k} = E\left[\mathbf{z}^{(r)}_{k-1}\eta_k^{(s)T}\right]$.

From (4.6) for $\eta^{(s)}_k$, we have $\mathcal{Z}^{(rs)}_{k-1, k} = \left(E[\mathbf{z}^{(r)}_{k-1}y_{k}^{T}]-\Sigma^{\mathbf{z}^{(rs)}}_{k-1}\left(\overline{C}_k{\cal A}_k \ | \ \mathcal{H}_k\right)^T\left(I_{mn_y}-\overline{\Lambda}_k\right)\right)D^{(s)T}_{y}.$ Now, using the OPL, we express $E\left[\mathbf{z}_{k-1}^{(r)}y_k^T\right] = E\left[\mathbf{z}_{k-1}^{(r)}\widehat{y}^{(r)T}_{k/k-1}\right]$, and from (7.3) for $\widehat{y}^{(r)}_{k/k-1}$, the following identity holds:

and (5.4) is obtained.

To end the proof, we again use (4.6) for $\eta^{(r)}_k$ together with the definition $\mathcal{Z}^{(rs)}_{k-1, k} = E\left[\mathbf{z}^{(r)}_{k-1}\eta_k^{(s)T}\right]$ to express the innovation cross-covariance as $\Sigma^{\eta^{(rs)}}_{k} = D^{(r)}_y\left(E\left[y_k\eta_{k}^{(s)T}\right]-\left(I_{mn_y}-\overline{\Lambda}_{k}\right)\left(\overline{C}_k{\cal A}_k \ | \ \mathcal{H}_k\right){\cal Z}^{(rs)}_{k-1, k}\right)$. Then, from (4.6) for $\eta^{(s)}_k$, and taking into account (13.1) for $E\left[y_k\mathbf{z}^{(s)T}_{k-1}\right]$, it is clear that:

and expression (5.5) is concluded simply by using (5.3) for $\mathcal{Z}^{(rs)}_{k-1, k}$. □

Appendix D. Proof of Lemma 2

Expression (5.6) for $\Sigma_{k/h}^{\widehat{x}^{(rs)}}, \ h > k,$ is directly obtained from (4.8), by denoting $\widehat{{\cal X}}^{(rs)}_{k, h} = E\left[\widehat{x}^{(r)}_{k/h-1}\eta_{h}^{(s)T}\right].$ In order to derive (5.7), we use (4.6) for $\eta_{h}^{(s)}$, and denote $\Sigma^{\widehat{x}\mathbf{z}^{(rs)}}_{k, h-1} = E\left[\widehat{x}^{(r)}_{k/h-1}\mathbf{z}_{h-1}^{(s)T}\right]$ to obtain $\widehat{{\cal X}}^{(rs)}_{k, h} = \left(E\left[\widehat{x}^{(r)}_{k/h-1}y_h^T\right]-\Sigma^{\widehat{x}\mathbf{z}^{(rs)}}_{k, h-1}\left(\overline{C}_h{\cal A}_h \ | \ \mathcal{H}_h\right)^T\left(I_{mn_y}-\overline{\Lambda}_h\right)\right)D^{(s)T}_{y}$. Then, taking into account that –from the OPL– $E\left[\widehat{x}^{(r)}_{k/h-1}y_h^T\right] = E\left[x_k\widehat{y}_{h/h-1}^{(r)T}\right]$, we need only use (7.3) for $\widehat{y}_{h/h-1}^{(r)}$.

Finally, expression (5.8) of $\Sigma^{\widehat{x}\mathbf{z}^{(rs)}}_{k, h}, \ h > k,$ is obtained straightforwardly simply by using (4.8) for $\widehat{x}^{(r)}_{k/h}$ and (4.3) for $\mathbf{z}_h^{(s)}$. Its initial condition is derived directly from (4.1) and $\Sigma_{k}^{\mathbf{z}^{(rs)}} = E\left[\mathbf{z}^{(r)}_k \mathbf{z}^{(s)T}_k \right].$ □