Centralized Fusion Estimation in Networked Systems: Addressing Deception Attacks and Packet Dropouts with a Zero-Order Hold Approach

1. Introduction

With the advancement of technology and scientific progress, networked systems have become crucial in many research domains, including health management, target tracking, traffic control, environmental monitoring, fault diagnosis and security monitoring (see, e.g., [1−3]). Among these systems, multi-sensor systems in particular have garnered significant research attention due to their inherent characteristics such as flexibility, high performance, efficiency, and resource sharing. The literature extensively addresses state estimation problems for sensor networks facing uncertainties in both measurement devices and transmission processes. For instance, in [4], the state estimation problem for a class of spatial-temporal networks with time-varying delays is studied under event-triggered mechanisms and encoding-decoding schemes. In [5], using a prediction compensation strategy, a distributed fusion filtering algorithm is proposed for discrete-time linear network systems with stochastic disturbances in sensor measurements, and both random one-step delays and random packet dropouts (RPD) in transmission processes. The optimal, suboptimal and probability-dependent distributed Kalman filters are proposed in [6] for multi-sensor networked systems under stochastic communication protocols and correlated noises. Comprehensive reviews of the key findings and emerging challenges in this field are provided in [7−11].

Fusion estimation, which effectively integrates diverse sensor data to produce a more accurate estimate of the system state, is a critical research area in multi-sensor networked systems. Broadly speaking, fusion estimation methods can be divided into two main categories —centralized and distributed—, based on whether sensor measurements are transmitted to a fusion center (FC) or not. In centralized fusion (CF) estimation, individual sensors forward their observations to a FC where all incoming data is collectively processed. This approach relies on the FC to produce a single optimal state estimator utilizing measurements from all sensors across the network, thereby achieving optimal results in the least-squares (LS) sense. Representative contributions to this topic can be found in [12−15] and references therein. Despite their ability to yield optimal state estimators, CF estimators generally suffer from limited robustness and computational drawbacks. In contrast, distributed fusion estimation designs local estimators from pre-processed data at each sensor, which are then fused based on a specific fusion criterion to obtain distributed fusion estimators. This approach, which is widely used due to its reduced computational burden and improved robustness, is discussed in [16−18]. However, distributed fusion estimation demands that sensors have the computational ability to pre-process data. In situations where this capability is lacking, CF estimation is more suitable, as it processes all raw sensor data at a central location and provides globally optimal estimates.

Networked systems are typically subject to different stochastic uncertainties due to physical constraints, environmental complexities, changes in subsystem interconnections and random component failures or repairs. A common type of uncertainty in these environments is the presence of multiplicative stochastic disturbances in the system model (such as multiplicative noises in the system state), missing observations or degradation of measurements. These network-induced phenomena —which occur in diverse application fields, such as digital control of chemical processes, radar control, navigation systems and economic systems— can be described using stochastic parameter matrices within the system equations. In addition to these uncertainties, data transmission over communication networks may be affected by imperfect communication channels or network congestion, resulting in stochastic uncertainties, such as random delays and RPD in the transmitted measurements. A critical issue when dealing with RPD is how to compensate for lost measurements. A widely used compensation method involves employing the last successfully transmitted measurement when the actual one is unavailable —the zero-order hold strategy (ZOHS)— which effectively mitigates the impact of lost measurements. Numerous studies have focused on designing estimation algorithms for systems with random parameter matrices and transmission uncertainties, covering a broad range of network-induced uncertainties as discussed above. Relevant examples include estimation studies with random parameter matrices, random transmission delays, and packet losses in [19−24] and references therein.

Another major challenge is the presence of sequentially time-correlated measurement noises, which frequently occur in fields such as electronics and engineering. Extensive research has focused on addressing this issue, often assuming that measurements are affected by infinite-step time-correlated channel noises, modeled as the output of a linear system driven by white noises. Two well-known approaches for managing this noise correlation are state augmentation (which is straightforward but requires high computational cost) and measurement differencing (which avoids increased dimensionality but needs two consecutive measurements to compute the difference and generate a new measurement free of time-correlated noises). The design of fusion estimation algorithms for multi-sensor systems with time-correlated channel noises is a critical research area. In [24], centralized and distributed fusion filtering and smoothing algorithms are derived for multi-sensor systems with random parameter matrices using the measurement differencing method. However, when random delays or RPD occur, sensor measurements may not arrive at the processor in time, making the measurement differencing method ineffective. Additionally, from a computational viewpoint, developing novel non-augmentation methods to deal with time-correlated measurement noises remains a significant challenge. In this context, alternative approaches that do not rely on augmentation or differencing, but instead focus on the direct estimation of time-correlated additive noise, are discussed in [25] for stochastic uncertain systems over packet-dropping networks and in [26] for nonlinear stochastic uncertain systems with RPD compensations. Furthermore, even in the absence of transmission losses, direct estimation of time-correlated additive noises has been applied in [27] for networked uncertain systems where the stochastic uncertainties are characterized by white multiplicative noises, and in [28, 29] for systems with random parameter matrices and deception attacks.

In addition to the various stochastic uncertainties affecting both measurements and transmissions in networked systems, research on the estimation problem must take into account the high likelihood of cyber-attacks. Security vulnerabilities have been widely examined in the literature, with comprehensive reviews of recent advances and challenges provided in [30, 31]. In particular, random deception attacks (RDA), which aim to undermine data integrity by maliciously and randomly altering information, have received considerable research attention. The distributed filtering problem in sensor networks with specific communication topologies under RDA is explored in [32], considering systems with fading measurements and multiplicative noises in both signal and measurement equations. Further investigations into sensor networks with various network-induced constraints and simultaneous RDA on sensor measurements are presented in [33]. The distributed estimation problem has also been studied in [34] for a specific class of nonlinear systems under RDA. Additionally, networked uncertain systems affected by RDA –with uncertainties arising from both multiplicative and additive noises in the state and measurement equations– are examined in [35] to develop distributed optimal and self-tuning filtering solutions utilizing compressed data.

Despite the extensive literature on sensor network estimation, the above discussion leads us to conclude that several interesting challenges remain when the system dynamics include mixed complexities. These challenges serve as the main motivation for our research, which focuses on the design of recursive CF filtering and fixed-point smoothing algorithms for a class of multi-sensor networked uncertain systems with the following features:

(a) The state evolution equation includes both additive and multiplicative noises.

(b) The measured outputs contain stochastic uncertainties, that are depicted by considering that the entries of the measurement matrices at the different sensors are random variables (rv).

(c) The additive noise of the sensor measurement equations has infinite-step correlation, which, as usual, is described by a first-order auto-regressive model.

(d) In addition to the presence of sensor measurement uncertainties, RPD and RDA during transmission are also considered.

The primary challenges of the addressed problem are: 1) how to compensate for lost measurements to improve the estimation accuracy under RPD; 2) how to effectively obtain recursive estimators when infinite-sept correlation of the sensor noises and RPD in transmission occur simultaneously, a challenge that cannot be addressed by subtracting two successive observations due to the random loss of data packets; 3) how to adequately address the risks posed by RDA, which can lead to significant estimation errors if not properly handled. The most outstanding contributions achieved in this work are closely related to the strategies used to overcome these challenges; namely, the use of the ZOHS to counterbalance the effect of RPD, the use of direct estimation of measurement noise to deal with time-correlation, and the incorporation of the effect of RDA into the estimation algorithm based on the knowledge of their probability of success, without needing to know whether a particular random attack was successful or not. In contrast to the above commonly used methods for dealing with the estimation problem in systems with time-correlated additive noise, we propose a novel approach that combines the estimators of the state process and the measurement noise to effectively handle the correlation of such processes. Another significant contribution of the current research is the use of the state and noise mean and covariance functions, instead of the explicit state and noise evolution models, together with an innovation-based methodology. This allows the algorithms to be easily derived while retaining their precision and simple structure.

The remainder of this paper is structured as follows: Section 2 introduces the system model, which accounts for random parameter matrices and time-correlated noises in sensor outputs, along with simultaneous RDA and RPD in transmissions. This section also outlines the assumptions made regarding the stochastic processes involved. In Section 3, the CF estimation problem is formulated, with detailed derivations of the filtering and fixed-point smoothing algorithms. Section 4 presents an illustrative example that demonstrates the effectiveness of the proposed estimation algorithms and examines how both RDA and RPD probabilities affect estimation accuracy. Finally, Section 5 concludes with a summary of the key findings.

Notation and abbreviations. The following standard notation is used throughout the paper. \( \mathbb{R}^{n} \) and \( \mathbb{R}^{m\times n} \) denote the sets of \( n \)-dimensional real vectors and \( m\times n \) real matrices, respectively. For a matrix \( A \), \( A^T \) and \( A^{-1} \) denote its transpose and inverse, respectively. \( Diag(a_1,\ldots, a_m) \) is the diagonal matrix with entries \( a_1, \ldots, a_m \), and \( (A_1 \ | \ A_2) \) represents the partitioned matrix with submatrices \( A_1 \) and \( A_2 \). \( I_{n} \) and \( \mathbf{1}_{n} \) denote the \( n\times n \) identity and all-ones matrices, respectively, and \( 0_{n\times m} \) is the \( n\times m \) all-zero matrix. \( \otimes \) and \( \circ \) denote the Kronecker and Hadamard products of matrices and \( \oplus \) represent the direct sum of matrices. If the dimensions of a vector or a matrix are not explicitly stated, they are assumed to be compatible with algebraic operations. For simplicity, we write \( G_{k}=G_{k,k} \) for any function \( G_{k,h} \), depending on time instants \( k \) and \( h \), when \( h=k \). \( \mathbb{E}[\cdot] \) denotes the mathematical expectation and \( P(B) \) is the probability of event \( B \). \( \delta_{k,h} \) denotes the Kronecker delta function. Finally, the following abbreviations are used:

2. Problem Formulation and System Description

This paper addresses the LS linear estimation problem for multi-sensor systems experiencing various stochastic uncertainties in sensor measurements. It further considers that the transmissions are simultaneously subject to RPD and RDA, which randomly alter the real measurements by replacing them with noises. The design of LS linear estimators is based on the CF architecture, in which sensor measurements are transmitted through different communication channels to a FC, where they are merged to obtain the required estimators. Utilizing an innovation approach and the ZOHS —which compensates for missing current measurements at the FC by using the last received measurement— covariance-based recursive algorithms for CF filtering and fixed-point smoothing estimators are derived. The system model and the assumptions necessary to address the LS linear estimation problem are detailed in the following subsections.

2.1. State and Measurement Models

Consider a class of networked multi-sensor systems with stochastic uncertainties characterized by random parameter matrices and multiplicative and additive noises, described by the state and measurement equations:

where \( \mathcal{M}=\{1,2,\cdots, m\} \) denotes the set of sensors, \( x_{k}\in \mathbb{R}^{n_x} \) is the state vector and \( z_{i,k}\in \mathbb{R}^{n_z} \) is the measured output collected by the \( i \)th sensor. The matrices \( \{A_k\}_{k\geq 0} \), \( \{\breve{A}_k\}_{k\geq 0} \) are known time-varying matrices.

All processes in equations (1) and (2) are assumed to be second-order and the following assumptions apply:

(I) The initial state vector \( x_{0} \) and the noise processes \( \{\alpha_k\}_{k\geq 0} \) and \( \{w_{k}\}_{k\geq 0} \) are mutually independent and verify:

● The initial state \( x_{0} \) is a zero-mean random vector with known covariance matrix, \( K^x_{0}=\mathbb{E}[x_0x_0^T] \).

● The multiplicative noise \( \{\alpha_k\}_{k\geq 0} \) is a zero-mean scalar white process with known variances, \( \sigma^2_{\alpha,k}=\mathbb{E}[\alpha^2_k], \ k\geq 0. \)

● The additive noise \( \{w_{k}\}_{k\geq 0} \) is a zero-mean white process with known covariance matrices, \( K^w_{k}=\mathbb{E}[w_{k}w^T_{k}], \ k\geq 0 \).

Remark 1. From equation (1) and the above assumptions, it is easy to prove that the state process \( \{x_k\}_{k\geq 0} \) has zero mean and its covariance matrices \( K^x_s= \mathbb{E}[x_sx_s^T] \) can be recursively obtained as:

Then, by denoting \( \mathbb{A}_k=\prod_{h=0}^{k-1}A_{h} \) and \( \mathbb{B}^T_s=\mathbb{A}_s^{-1}K^x_s \), the state covariance function \( K^x_{k,s}=\mathbb{E}[x_{k}x^T_{s}] \) can be expressed as: \( K^x_{k,s}=\mathbb{A}_k \mathbb{B}^T_s,\ \ s\leq k. \)

(II)\( \{H_{i,k}\}_{k\geq 1}, \ i\in \mathcal{M} \), are independent sequences of independent random parameter matrices whose entries \( h_{i,pq}(k) \), \( p=1,\cdots, n_z \) and \( q=1,\cdots, n_x \), have known means, \( \mathbb{E}[h_{i,pq}(k)] \). The expectations \( \mathbb{E}[h_{i,pq}(k)h_{i,p' q'}(k)\big] \) are also assumed to be known, for \( p,p'=1,\cdots, n_z \) and \( q, q'=1,\cdots, n_x \). The mean matrices are denoted by \( \overline{H}_{i,k} = \mathbb{E}[H_{i,k}] \), with \( \mathbb{E}[H_{i,k}]=\big(\mathbb{E}[h_{i,pq}(k)]\big)_{n_z\times n_x} \).

Remark 2. From the independence hypotheses of assumption (Ⅱ), for arbitrary \( i,j\in\mathcal{M} \) and any deterministic matrix \( R=\big(r_{pq}\big)_{n_x\times n_x} \), it is clear that \( \mathbb{E}\big[H_{i,k}RH^T_{j,s}\big]=\overline{H}_{i,k}R\overline{H}^T_{j,s} \) for \( j\neq i \) or \( s\neq k \). When \( j= i \) and \( s=k \), the entries of the matrix of \( \mathbb{E}\big[H_{i,k}RH^T_{i,k}\big] \) are computed by

(III) The measurement noises \( \{v_{i,k}\}_{k\geq 0} \), \( i\in \mathcal{M} \), are time-correlated processes generated by:

where \( C_{i,k} \in \mathbb{R}^{n_z\times n_z} \), \( i\in \mathcal{M} \), are non-singular known deterministic matrices. The initial vectors \( v_{i,0}, \ i\in {\cal M}, \) and the noises \( \{u_{i,k}\}_{k\geq 0}, \ i\in {\cal M}, \) are independent and they verify:

● The initial noises \( v_{i,0}, \ i\in {\cal M}, \) are zero-mean random vectors with known covariance and cross-covariance matrices, \( K^{v}_{ij,0}=\mathbb{E}[v_{i,0} v^T_{j,0}], \ i,j\in {\cal M}. \)

● The noises \( \{u_{i,k}\}_{k\geq 0}, \ i\in {\cal M}, \) are zero-mean white processes, independent of each other at different times, with known covariance and cross-covariance functions \( K^{u}_{ij,k}=\mathbb{E}[u_{i,k}u^T_{j,k}], \ k\geq 0; \ i,j\in {\cal M} \).

Remark 3. From this assumption, it is evident that the measurement noises \( \{v_{i,k}\}_{k\geq 0}, \ i\in \mathcal{M}, \) have zero mean and their cross-covariance matrices \( K^{v}_{ij,k,s}=\mathbb{E}[v_{i,k}v^T_{j,s}] \), \( i, j\in \mathcal{M}, \) can be expressed as \( K^v_{ij,k,s}=C_{i,k-1}\cdots C_{i,s}K^{v}_{ij,s}, \ s\leq k, \) where the matrices \( K^{v}_{ij,s}= \mathbb{E}[v_{i,s}v^T_{j,s}] \) are recursively obtained as:

2.2. Measurements Subject to Random Deception Attacks and Random Packet Losses

Consider the scenario in which the sensor measurements \( \{z_{i,k}\}_{k\geq 1}, \ i\in \mathcal{M}, \) are transmitted to the remote FC via individual communication channels subject to both RDA and RPD during data transmission. Binary indicator rv, \( \lambda_{i,k} \in \{0,1\} \) and \( \gamma_{i,k} \in \{0,1\} \), are incorporated into the observation model to represent these random phenomena. More specifically:

* Random deception attacks: The sensor measurements subject to RDA, denoted as \( \vec{z}_{i,k} \), are described by

where \( \lambda_{i,k} \in \{0,1\} \) is the successful attack indicator for the \( i \)th sensor at time \( k \) and \( \breve{z}_{i,k}=- z_{i,k} + \varepsilon_{i,k} \) is the RDA signal introduced by the attacker to nullify the actual measurement \( z_{i,k} \) and replace it with deceptive information represented by noise \( \varepsilon_{i,k} \). Expression (4) for the attacked measurement outputs can then be equivalently rewritten as follows:

\( \ast \) Random packet dropouts: Consider that data packets are simultaneously subject to random losses during transmission to the FC. Let \( \gamma_{i,k} \in \{0,1\} \) be the packet arrival indicator of the \( i \)th sensor measurement at the FC at time \( k \). To mitigate the effect of these packet losses, the ZOHS is adopted; specifically, the most recently received data packet is considered to counterbalance packet losses. Thus, denoting \( {y}_{i,k} \) the measurements processed for estimation under ZOHS, the following equation holds:

The next assumption is made regarding the processes involved in equations (5) and (6):

(IV) The attack noises \( \{\varepsilon_{i,k}\}_{k\geq 1} \) and the sequences \( \{\lambda_{i,k}\}_{k\geq 1} \) and \( \{\gamma_{i,k}\}_{k\geq 1} \) satisfy:

● \( \{\varepsilon_{i,k}\}_{k\geq 1}, \ i\in {\cal M}, \) are zero-mean white processes, independent of each other at different times, with known covariance and cross-covariance functions, \( K^{\varepsilon}_{ij,k}=\mathbb{E}[\varepsilon_{i,k}\varepsilon^T_{j,k}], \ k\geq 1; \ i,j\in {\cal M} \).

● Both \( \{\lambda_{i,k}\}_{k\geq 1}, \ i\in {\cal M}, \) and \( \{\gamma_{i,k}\}_{k\geq 1}, \ i\in {\cal M}, \) are independent sequences of independent Bernoulli rv with known probabilities \( P(\lambda_{i,k}=1)=\overline{\lambda}_{i,k} \) and \( P(\gamma_{i,k}=1)=\overline{\gamma}_{i,k} \).

Remark 4. From the above independence assumption, it is clear that \( \mathbb{E}[\lambda_{i,k}\gamma_{j,s}]=0 \), \( \forall i,j\in \mathcal{M} \), \( \forall k,s\geq 1 \), and

Finally, the following independence assumption is required:

(V) For each \( i\in {\cal M} \), the state process, \( \{x_k\}_{k\geq 0} \), the random parameter matrices, \( \{H_{i,k}\}_{k\geq 1} \), the noises, \( \{v_{i,k}\}_{k\geq 1} \) and \( \{\varepsilon_{i,k}\}_{k\geq 1} \), and the sequences \( \{\lambda_{i,k}\}_{k\geq 1} \) and \( \{\gamma_{i,k}\}_{k\geq 1} \), are mutually independent.

3. Centralized Fusion Estimators

In this section, the proposed methodology to address the CF LS linear estimation problem for the system model presented in Section 2 is described, and recursive filtering and smoothing algorithms are designed. Actually, using an innovation approach, we derive recursive algorithms for the LS linear estimators \( \widehat{x}_{k/N} \) of the state, \( x_k \), based on the measurements \( \{ y_{i,1}, \cdots, y_{i,N}, \ i\in \mathcal{M}\} \) defined by (6), using the CF methodology. Specifically, we start by deriving a recursive algorithm for the LS linear filter, \( \widehat{x}_{k/k} \), based on the measurements \( \{ y_{i,1}, \cdots, y_{i,k}, \ i\in \mathcal{M}\} \). Then, for each fixed \( k \), when successive measurements \( y_{i,k+1}, y_{i,k+2}, \cdots, y_{i,k+L}, \ i\in \mathcal{M} \) become available, the filtering estimator will be recursively updated to obtain the fixed-point smoothing estimators \( \widehat{x}_{k/k+L}, \ L\geq 1 \).

The section is organized as follows. In Subsection 3.1, the stacked observation model is introduced to address the estimation problem under the CF architecture. Subsection 3.2 discusses the implementation of the innovation technique and its integration with the proposed estimation method. Finally, Subsection 3.3 details the derivation of the filtering and smoothing algorithms.

3.1. Stacked Observation Model: Properties

At the FC, the stacked vectors \( y_k=\big({y}^{T}_{1,k}, \cdots, {y}^{T}_{m,k}\big)^T \) with all the measurement information used for estimation are considered. Next, we establish the mathematical model for these vectors and examine the properties that will be necessary to address the LS linear estimation problem.

The following notation is used for the stacked vectors and matrices:

From equation (6), the mathematical model for the stacked vectors \( y_k \) is:

where, from (5), the vectors \( \vec{z}_k \) of measurements subject to RDA are:

and, from (2), the vectors \( z_k \) are given by:

Note that, from (3), the additive noise \( \{v_k\}_{k\geq 0} \) is a time-correlated sequence, satisfying

with \( C_k=C_{i,k}\oplus \cdots\oplus C_{m,k} \) and \( u_k=\big(u^{T}_{1,k}, \cdots, u^{T}_{m,k}\big)^T \).

Once we have established the model for the stacked vectors, we must derive the LS linear estimators \( \widehat{x}_{k/N} \) of the state, \( x_k \), based on the measurements \( \{y_1, \cdots, y_N\} \) given by (7). For this purpose, the statistical properties of the processes in the observation model (7)-(10) need to be studied. These properties are summarized in the following proposition, whose proof is straightforward from assumptions (I)-(V).

Proposition 1. The processes involved in equations (7)-(10) satisfy the following properties:

\( {\bf (a)} \)\( \left\{H_k\right\}_{k\geq 1} \) is a sequence of independent random parameter matrices with know means \( \overline{H}_k= \big(\overline{H}_{1,k}^T \ \cdots \ \overline{H}_{m,k}^T\big)^T \) and, for any deterministic matrix \( R \), the expectations \( \mathbb{E}\left[H_{k}RH_{s}^{T}\right]=\Big(\mathbb{E}[H_{i,k}RH^{T}_{j,s}]\Big)_{i,j\in \mathcal{M}} \) are also known. Specifically, their entries \( \mathbb{E}[H_{i,k}RH^{T}_{j,s}], \ i,j\in \mathcal{M} \), are defined in Remark 2.

\( {\bf (b)} \) The measurement noise, \( \left\{v_k\right\}_{k\geq 1} \), is a zero-mean time-correlated sequence, with covariance matrices \( K^v_{k,s}= \mathbb{E}[v_kv_s^T] \) given by \( K^v_{k,s}=C_{k-1}\cdots C_sK_{s}^{v}, \ s<k \), where \( K^v_{s}= \mathbb{E}[v_kv_s^T] \) is recursively obtained as: \( K^v_s=C_{s-1}K^v_{s-1} C_{s-1}^T+ K^u_{s-1}, \ \ s\geq 1, \) with \( K^{u}_{s}=\left(K^u_{ij,s}\right)_{i,j\in {\cal M}} \) and initial condition \( K^v_{0}=\left(K^v_{ij,0}\right)_{i,j\in {\cal M}} \). Moreover, denoting \( \mathbb{C}_k=\prod_{h=0}^{k-1}C_{h} \) and \( \mathbb{D}^T_s=\mathbb{C}_s^{-1}K^v_s \), the noise covariance function \( K^v_{k,s} \) can be factorized as: \( K^v_{k,s}=\mathbb{C}_k \mathbb{D}^T_s,\ \ s\leq k. \)

\( {\bf (c)} \) The attack noise \( \{\varepsilon_k\}_{ k\geq1} \) is a zero-mean white process with covariance matrices \( K^{\varepsilon}_{k}= \mathbb{E}[\varepsilon_k \varepsilon^T_k]= \big(K^{w}_{ij,k}\big)_{i,j\in \mathcal{M}} \).

\( {\bf (d)} \)\( \{\Lambda_{k}\}_{ k\geq1} \) and \( \{\Gamma_{k}\}_{ k\geq1} \) are sequences of diagonal independent random matrices with know means

Denoting \( \lambda_{k}= \big(\lambda_{1,k},\cdots, \lambda_{m,k}\big)^T\otimes \mathbf{1}_{n_z} \) and \( \gamma_{k}= \big(\gamma_{1,k},\cdots, \gamma_{m,k}\big)^T\otimes \mathbf{1}_{n_z} \), \( k\geq 1, \) the matrices \( K^{\beta}_{k}= \mathbb{E}\big[\beta_{k} \beta^T_{k}\big] \), \( K^{\mathbf{1}-\beta}_{k}= \mathbb{E}\big[(\mathbf{1}_{mn_z}-\beta_{k})(\mathbf{1}_{mn_z}- \beta_{k})^{T}\big] \) and \( K^{\beta(\mathbf{1}-\beta)}_{k}= \mathbb{E}\big[\beta_{k}(\mathbf{1}_{mn_z}- \beta_{k})^{T}\big] \), with \( \beta_k=\lambda_k, \gamma_k \), are known and their entries can be computed taking into account Remark 4.

\( {\bf (e)} \) The state process \( \{x_k\}_{k\geq 1} \), the random parameter matrices, \( \{H_k\}_{k\geq 1} \), the noises, \( \left\{v_k\right\}_{k\geq 1} \) and \( \left\{\varepsilon_k\right\}_{ k\geq 1} \), and the diagonal random matrices \( \{\Lambda_k\}_{k\geq 1} \) and \( \{\Gamma_{k}\}_{ k\geq1} \), are all mutually independent.

The covariance and cross-covariance functions of the observation processes involved in equations (7)-(9) will be used in the filtering algorithm. The expressions for such functions are provided in the following proposition, whose proof follows easily from the aforementioned properties.

Proposition 2. The sequences \( \left\{\vec{z}_k\right\}_{k\geq 1} \) and \( \left\{y_k\right\}_{k\geq 1} \) are zero-mean second-order processes whose covariance and cross-covariance functions, \( K^{\vec{z}}_{k,s}= \mathbb{E}[\vec{z}_k\vec{z}_s^T] \), \( K^y_{k,s}= \mathbb{E}[y_ky_s^T] \) and \( K^{\vec{z}y}_{k,s}= \mathbb{E}[\vec{z}_ky_s^T] \), are given by the following expressions:

● \( K^{\vec{z}}_{k}=K^{\mathbf{1}-\lambda}_{k}\circ \big(\mathbb{E}[H_k \mathbb{A}_k \mathbb{B}_k^T H_k^T ]+ {\mathbb C}_k {\mathbb D}_k^T\big) + K^{\lambda}_{k}\circ K^{\varepsilon}_{k}, \ \ k\geq 1. \)

● \( K^{\vec{z}}_{k,s}=(I_{mn_z}-\overline{\Lambda}_{k})\big(\overline{H}_k \mathbb{A}_k \mathbb{B}_s^T \overline{H}_s^T + {\mathbb C}_k {\mathbb D}_s^T\big)(I_{mn_z}-\overline{\Lambda}_{s}), \ \ s\leq k-1. \)

● \( K^{\vec{z}y}_{k,s}=K^{\vec{z}}_{k,s}\overline{\Gamma}_{s}+ K^{\vec{z}y}_{k,s-1}(I_{mn_z}-\overline{\Gamma}_{s}), \ \ \ 2 \leq s \leq k-1; \ \ \ K^{\vec{z}y}_{k,1}= K^{\vec{z}}_{k,1}\overline{\Gamma}_{1}. \)

● \( K^y_{k}=K^{\gamma}_{k}\circ K^{\vec{z}}_{k}+K^{(\mathbf{1}-\gamma)}_{k}\circ K^{y}_{k-1} +K^{\gamma(\mathbf{1}-\gamma)}_{k}\circ K^{\vec{z}y}_{k,k-1}+\big(K^{\gamma(\mathbf{1}-\gamma)}_{k}\circ K^{\vec{z}y}_{k,k-1}\big)^T , \ \ k\geq 2;\ \ \ \ K^y_1=K^{\gamma}_{1}\circ K^{\vec{z}}_{1}. \)

● \( K^y_{k,k-1}=\overline{\Gamma}_{k}K_{k,k-1}^{\vec{z}y}+(I_{mn_z}-\overline{\Gamma}_{k})K_{k-1}^{y}, \ \ k\geq 2. \)

Remark 5. For the sake of notational simplicity in the derivation of the filtering algorithm, let us consider a new process \( \left\{g_k\right\}_{k\geq 1} \) defined as \( g_{k}=y_k-(I_{mn_z}-\overline{\Gamma}_{k})y_{k-1}, \ k\geq 2; \ \ g_{1}=y_1. \) This process has zero mean and their covariance matrices, \( K^g_k=\mathbb{E}[g_kg_k^T], \ k\geq 1 \), are given by:

The statistical properties of the processes involved in the stacked model, established in propositions 1 and 2, ensure the existence of the LS linear estimators, the design of which is detailed in the next subsection.

3.2. Innovation Approach to the LS Estimation Problem

Given the observations \( \{y_1,\cdots ,y_N\} \), the goal is to derive recursive algorithms for the CF LS linear estimators, \( \widehat{x}_{k/N} \), of the state \( x_k \) based on these observations. As the observations are generally nonorthogonal vectors, the Gram-Schmidt orthogonalization procedure is employed to transform the observation process, \( \left\{y_k\right\}_{k\geq 1} \), into a white process \( \left\{\mu_k\right\}_{k\geq 1} \), referred to as the innovation process. The estimators can then be expressed as a linear combination of the innovations and the orthogonality of this new process greatly simplifies the derivation of the algorithms compared to using the observations directly.

Specifically, replacing the observation process with the innovation process enables the LS linear estimate \( \widehat{\xi}_{k/N} \)(of an arbitrary second-order vector \( \xi_k \) based on the observations \( \left\{y_h, \ h\leq N\right\} \)) to be expressed as the following linear combination of the innovations \( \left\{\mu_h, \ h\leq N\right\} \):

where \( {\Pi}_h=\mathbb{E}[\mu_{h}\mu_{h}^{T}] \) denotes the innovation covariance matrix.

The innovation at time \( k \) is defined as \( \mu_{k}=y_{k}-\widehat{y}_{k/k-1} \), where \( \widehat{y}_{k/k-1} \) is the LS linear one-step predictor of the observation \( y_{k} \). From the Orthogonal Projection Lemma (OPL), equations (7)-(9) and the properties established in Proposition 1, we have that \( \widehat{y}_{k/k-1}=\overline{\Gamma}_k\widehat{\vec{z}}_{k/k-1}+(I_{mn_z}-\overline{\Gamma}_k) y_{k-1}, \ \ k\geq2, \) with \( \widehat{\vec{z}}_{k/k-1}=(I_{mn_z}-\overline{\Lambda}_k)\widehat{z}_{k/k-1} \) and \( \widehat{z}_{k/k-1}= \overline{H}_k\widehat{x}_{k/k-1}+\widehat{v}_{k/k-1} \).

Thus, the one-step predictor \( \widehat{y}_{k/k-1} \) is expressed in terms of the state predictor, \( \widehat{x}_{k/k-1} \), and the noise predictor, \( \widehat{v}_{k/k-1} \); specifically, \( \widehat{y}_{k/k-1} \) is given by:

For simplicity, we denote \( \widehat{\mathbf{\Psi}}_{k/k-1}=\left(\hskip -0.1cm \begin{array}{c}\widehat{x}_{k/k-1}\\ \widehat{v}_{k/k-1}\end{array}\hskip -0.1cm\right) \) the one-step predictor of the vector \( \mathbf{\Psi}_{k}=\left(\hskip -0.1cm \begin{array}{c}x_k\\ v_k\end{array}\hskip -0.1cm\right) \). Then, expression (12) can be equivalently rewritten as:

Therefore, the determination of the innovation requires the linear one-step predictor \( \widehat{\mathbf{\Psi}}_{k/k-1} \) which, together with expression (13) provides the foundation for deriving the CF filtering and fixed-point smoothing algorithms presented in the following subsection.

3.3. Centralized Fusion Filtering and Fixed-Point Smoothing Algorithms

In this section, we first present a recursive algorithm for the CF LS linear filter, \( \widehat{x}_{k/k} \), of the state \( x_k \) based on the observations \( y_1, \cdots, y_k \). Then, in order to obtain the smoothers at the fixed-point \( k \), \( \widehat{x}_{k/k+L} \), \( L\geq 1 \), the filter \( \widehat{x}_{k/k} \) will be recursively updated as successive observations \( y_{k+1}, y_{k+2}, \cdots, y_{k+L} \) become available. The performance of the CF estimators is measured by the magnitude of the estimation errors, \( \widetilde{x}_{k/k+L}=x_k-\widehat{x}_{k/k+L} \) and, more specifically, by their covariance matrices, \( K^{\tilde{x}}_{k/k+L}=\mathbb{E}\big[\widetilde{x}_{k/k+L}\widetilde{x}_{k/k+L}^T\big] \), \( L\geq 0 \). Both proposed algorithms also provide recursive formulas for the estimation error covariance matrices, thereby offering a measure of the estimation accuracy under the LS optimality criterion.

Recursive centralized fusion filtering algorithm. Under assumptions (I)-(V), the CF LS linear filtering estimator, \( \widehat{x}_{k/k} \), and the corresponding error covariance matrix, \( K^{\tilde{x}}_{k/k} \), are given by:

where the vectors \( \mathbf{e}_k \) and the matrices \( K^{\mathbf{e}}_k= \mathbb{E}\big[\mathbf{e}_k\mathbf{e}^{T}_k\big] \) are recursively obtained from:

and the matrices \( \mathbf{\Phi}_k= \mathbb{E}\big[\mathbf{e}_k\mu_k^{T}\big] \) satisfy:

The innovations, \( \mu_k \), and their covariance matrices, \( {\Pi}_k \), are calculated by:

with \( K^g_k \) given in Remark 5.

Proof. To derive the innovation from (13), we need to calculate the predictor \( \widehat{\mathbf{\Psi}}_{k/k-1} \). Given that the filter of the state vector is \( \widehat{x}_{k/k}=\left(I_{n_x} \ \ 0_{n_x\times mn_z}\right)\widehat{\mathbf{\Psi}}_{k/k} \) (which consists of the first \( n_x \) entries of \( \widehat{\mathbf{\Psi}}_{k/k} \)), we must obtain both the prediction and filtering estimators \( \widehat{\mathbf{\Psi}}_{k/s}, \ s\leq k. \) Using expression (11) for the LS linear estimators, we write

so the coefficients \( \mathcal{B}_{k,h}=\mathbb{E}\big[\mathbf{\Psi}_{k}\mu^{T}_{h}\big]=\mathbb{E}\big[\mathbf{\Psi}_{k}y^{T}_{h}\big]-\mathbb{E}\big[\mathbf{\Psi}_{k}\widehat{y}^{T}_{h/h-1}\big], \ 1 \leq h \leq k \), need to be computed.

− On the one hand, by using equations (7)-(9), Remark 1 and the properties outlined in Proposition 1, we have:

− On the other hand, by employing (13) for \( \widehat{y}_{h/h-1} \) in conjunction with (21) for the predictor \( \widehat{\mathbf{\Psi}}_{h/h-1} \), we obtain:

Consequently,

ensuring that \( {\mathcal B}_{k,h}=(\mathbb{A}_k\oplus \mathbb{C}_k)\mathbf{\Phi}_h, \ 1 \leq h \leq k, \) where

Now, by defining the vectors \( \mathbf{e}_k= {\displaystyle\sum\limits_{h=1}^{k}}\mathbf{\Phi}_h\Pi^{-1}_h\mu_h, \ k\geq 1 \), we have that \( \widehat{\mathbf{\Psi}}_{k/s}=(\mathbb{A}_k\oplus\mathbb{C}_k)\ \mathbf{e}_s, \ s\leq k. \)

After these preliminary steps, the equations of the filtering algorithm are proven as follows:

● Since \( \widehat{x}_{k/k}=\left(I_{n_x} \ \ 0_{n_x\times mn_z}\right)\widehat{\mathbf{\Psi}}_{k/k} \), using that \( \widehat{\mathbf{\Psi}}_{k/k}=(\mathbb{A}_k\oplus\mathbb{C}_k)\ \mathbf{e}_{k-1} \), expression (14) is clear.

● Taking into account that \( K^{\mathbf{e}}_k= \mathbb{E}\big[\mathbf{e}_k\mathbf{e}^{T}_k\big] \), along with Remark 1 and (14), the expression for the filtering error covariance matrices \( K^{\tilde{x}}_{k/k}=\mathbb{E}[x_kx_k^T]- \mathbb{E}[\widehat{x}_{k/k}\widehat{x}^{T}_{k/k}] \) given in (15) is readily obtained.

● From their definition, the vectors \( \mathbf{e}_k \) clearly satisfy (16).

● The recursive formula (17) for the covariance matrices \( K^{\mathbf{e}}_k= \mathbb{E}\big[\mathbf{e}_k\mathbf{e}_k^{T}\big] \) is derived directly from (16), taking into account that the innovation is a white process.

● Using again the definition of the vectors \( \mathbf{e}_k \) and taking into account that the innovations form a white process, it is clear that \( K^{\mathbf{e}}_k= {\displaystyle\sum\limits_{h=1}^{k}}\mathbf{\Phi}_h\Pi^{-1}_h\mathbf{\Phi}^{T}_h \). This expression for \( K^{\mathbf{e}}_k \) along with (22), immediately leads to (18).

● By substituting \( \widehat{\mathbf{\Psi}}_{k/k-1}=(\mathbb{A}_k\oplus\mathbb{C}_k)\ \mathbf{e}_{k-1} \) in (13), the expression for the innovation in (19) is directly obtained.

● To derive expression (20), we use the notation \( g_{k}=y_k- (I_{mn_z}-\overline{\Gamma}_k) y_{k-1} \), introduced in Remark 5, to express the innovation as \( \mu_k= g_k-\overline{\Gamma}_k\widehat{\vec{z}}_{k/k-1} \). Given that \( \mathbb{E}\big[g_k\widehat{\vec{z}}^T_{k/k-1}\big]\overline{\Gamma}_k = \overline{\Gamma}_k\mathbb{E}\big[\vec{z}_{k}\widehat{\vec{z}}^T_{k/k-1}\big]\overline{\Gamma}_k \) and, by the OPL, \( \mathbb{E}\big[\vec{z}_{k}\widehat{\vec{z}}^T_{k/k-1}\big]= \mathbb{E}\big[\widehat{\vec{z}}_{k/k-1}\widehat{\vec{z}}^T_{k/k-1}\big] \), we obtain \( \Pi_k=\mathbb{E}\big[g_kg_k^T\big]-\overline{\Gamma}_k \mathbb{E}\big[\widehat{\vec{z}}_{k/k-1} \widehat{\vec{z}}_{k/k-1}^T\big]\overline{\Gamma}_k \). Thus, expression (20) is derived straightforwardly, taking into account that \( K^g_k=\mathbb{E}\big[g_kg_k^T\big] \) and \( \widehat{\vec{z}}_{k/k-1}=(I_{mn_z}-\overline{\Lambda}_k)\left(\ \overline{H}_k\mathbb{A}_k \ | \ \mathbb{C}_k\right)\mathbf{e}_{k-1} \).

Recursive centralized fixed-point smoothing algorithm. Under assumptions (I)-(V), starting at any fixed sampling time \( k\geq 1 \) with the filter, \( \widehat{x}_{k/k} \), and its error covariance matrix, \( K^{\tilde{x}}_{k/k} \), as initial conditions, the CF LS linear smoothing estimators, \( \widehat{x}_{k/k+L}, \ L\geq 1 \), and their error covariance matrices, \( K^{\tilde{x}}_{k/k+L}, \ L\geq 1 \), are computed using the following recursive relations:

In these relations, the matrices \( G_{k,k+L}=\mathbb{E}[x_k\mu^{T}_{k+L}] \) satisfy;

and \( \mathbf{F}_{k,k+L}=\mathbb{E}[x_k\mathbf{e}^{T}_{k+L}] \) is recursively obtained from:

Proof. The formulas of the fixed-point smoothing algorithm are proven in the following steps:

● From the general expression (11), it is evident that, at a fixed sampling time \( k\geq 1 \), the smoother \( \widehat{x}_{k/k+L}, \ L\geq 1 \), can be recursively calculated by (23), starting from the linear filter, \( \widehat{x}_{k/k} \), as initial condition.

● Using (23), the smoothing errors can be expressed as \( \widetilde{x}_{k/k+L}=\widetilde{x}_{k/k+L-1}- G_{k,k+L}\Pi^{-1}_{k+L}\mu_{k+L} \), and (24) for the error covariance matrices, \( K^{\tilde{x}}_{k/k+L} \), is clear.

● Expression (25) for the coefficients

is obtained by denoting \( \mathbf{F}_{k,k+L}=\mathbb{E}[x_k\mathbf{e}^{T}_{k+L}] \) and taking into account that:

– \( \mathbb{E}[x_{k}\vec{z}^{T}_{k+L}]=\mathbb{B}_k\mathbb{A}^T_{k+L}\overline{H}^T_{k+L}(I_{mn_z}-\overline{\Lambda}_{k+L}) \)

– \( \mathbb{E}\big[x_k\widehat{\vec{z}}_{k+L/k+L-1}^T\big]=\mathbb{E}\big[ x_k\mathbf{e}^{T}_{k+L-1}\big]\left( \overline{H}_{k+L}\mathbb{A}_{k+L} \ | \ \mathbb{C}_{k+L} \right)^T(I_{mn_z}-\overline{\Lambda}_{k+L}). \)

● Finally, the recursion (26) for the matrices \( \mathbf{F}_{k,k+L} \) is easily derived using (16) for \( \mathbf{e}_{k+L} \). The initial condition is obtained by applying the LPO, from which \( \mathbb{E}\left[ x_k\mathbf{e}^{T}_{k}\right]=\mathbb{E}\left[ \widehat{x}_{k/k}\mathbf{e}^{T}_{k}\right] \), and using (14) for \( \widehat{x}_{k/k} \).

Remark 6. It is important to note that the state evolution equation (1) has not been directly used in the design of the proposed LS filtering and smoothing algorithms. Instead, it was utilized only to derive the state covariance matrix in a separable form (Remark 1), which serves as the foundation for the algorithms derivation, together with the factorization of the time-correlated noise covariance (Proposition 1). This makes our methodology applicable to a wide variety of models, as long as they meet this covariance factorization property. The structure of the proposed algorithms, designed on the basis of covariance information, differs from that of conventional LS linear estimation algorithms based on the state-space model, but both approaches provide optimal linear estimators. Nevertheless, although the proposed filtering and smoothing algorithms have attractive recursive structures due to the factorization property of both state and noise covariance functions, the presence of mixed uncertainties causes additional difficulties in deriving simple formulas for the innovation covariance matrix, \( \Pi_k \) (expression (20)). The definition of the vectors \( \widehat{\mathbf{\Psi}}_{k/k-1}=\left(\hskip -0.1cm \begin{array}{c}\widehat{x}_{k/k-1}\\ \widehat{v}_{k/k-1}\end{array}\hskip -0.1cm\right) \), which combine the estimators of the state process and the measurement noise, overcomes this issue and makes the derivation of the algorithm significantly simpler.

4. Numerical Simulation Example

In this section, a four-sensor system with random uncertainties is considered to show the validity of the designed estimation algorithms and assess the impact of different random uncertainties on the estimation accuracy. Specifically, the following model is considered:

State equation: Described by (1), with \( A_k=0.9, \ \breve{A}_k=0.05 \). The initial state, \( x_0 \), and the rv in the sequences \( \left\{ \alpha_k\right\}_{k \geq 1} \) and \( \left\{ w_k\right\}_{k \geq 1} \) obey the standard normal distribution.

Sensor measurement equations: Described by (2), with \( \mathcal{M}=\{1,2,3,4\} \), \( H_{1,k}= 0.9\vartheta_{1,k},\ H_{2,k}= 0.8\vartheta_{2,k}, \ H_{3,k}= 0.7\vartheta_{3,k} \) and \( H_{4,k}= 0.8\vartheta_{4,k}, \) where \( \{\vartheta_{i,k}\}_{k\geq 1} \), \( i\in\mathcal{M} \), are sequences of independent rv with the following distributions:

− \( \vartheta_{1,k} \) and \( \vartheta_{2,k} \) have uniform distributions over the intervals \( [0.1, 0.8 ] \) and \( [0.2, 0.9 ] \), respectively.

− \( \vartheta_{3,k} \) are discrete rv with \( P(\vartheta_{3,k}=0)=0.1, \ P(\vartheta_{3,k}=0.5)=0.5 \), \( P(\vartheta_{3,k}=1)=0.4, \ \forall k\geq 1. \)

− \( \vartheta_{4,k} \) have Bernoulli distributions with \( P(\theta_{4,k}=1)=\overline{\theta}_{4}, \ \forall k\geq 1. \)

The noises \( \{v_{k}^{(i)}\}_{k\geq 1} \), \( i\in \mathcal{M} \), are defined in (3), with \( C_{1,k}=C_{4,k}=0.7 \) and \( C_{2,k}=C_{3,k}=0.6 \), \( \forall k\geq 1 \). The noise \( u_{i,k} \) satisfies \( u_{i,k}=a_{i}\xi_{k} \), with \( a_1=a_3=0.25 \), \( a_2=a_4=0.5 \) and the rv in the white sequence \( \{\xi_{k}\}_{k\geq 0} \) obey the standard normal distribution. For \( i\in \mathcal{M} \), \( v_{i,0}=\nu_0 \) and \( \nu_0 \) is also a standard normal rv.

Random deception attacks and packet losses: The random sequences \( \{\varepsilon_{i,k}\}_{k\geq 1} \) in (5) are defined as \( \varepsilon_{i,k}=\varepsilon_i\zeta_{k} \), where \( \varepsilon_1=0.25 \), \( \varepsilon_2=0.5 \), \( \varepsilon_3=\varepsilon_4=0.75 \), and the rv in the white sequence \( \{\zeta_{k}\}_{ k\geq 1} \) obey the standard normal distribution; so, \( K^{\varepsilon}_{ij,k}= \varepsilon_i\varepsilon_j, \ i,j\in \mathcal{M}. \)\( \big\{\lambda_{i,k}\big\}_{ k\geq1} \) and \( \big\{\gamma_{i,k}\big\}_{ k\geq1} \), are white sequences of Bernoulli rv with \( P(\lambda_{i,k}=1)=\overline{\lambda}, \ \ P(\gamma_{i,k}=1)=\overline{\gamma}, \ \ \ k\geq 1, \ \ i\in \mathcal{M} \).

Performance analysis. Assuming \( \overline{\vartheta}_4=\overline{\lambda}=\overline{\gamma}=0.5 \), Figure 1 displays the error variances of the proposed estimators and shows that the smoothers outperform the filters. It is also clear from this figure that the accuracy of the smoothers improves as \( L \) increases and this fact is particularly notable for \( L\leq 6 \), beyond which the difference in accuracy becomes practically negligible.

Figure 1. Filtering and smoothing estimation error variances for \( \overline{\vartheta}_4=\overline{\lambda} = \overline{\gamma}=0.5 \).

Influence of missing probability, \( 1-\overline{\vartheta}_4 \), at sensor 4. Considering again \( \overline{\lambda}=\overline{\gamma}=0.5 \) and different values of \( \overline{\vartheta}_4 \), Figure 2 presents the filtering and smoothing (\( L=1,3 \)) error variances. From this figure, we infer that the probability \( \overline{\vartheta}_4 \) of the state being present in the measured outputs of sensor 4 significantly influences the estimation precision. More specifically, since the error variances decrease as \( \overline{\vartheta}_4 \) increases, we conclude that the proposed estimators perform better when the missing probability, \( 1-\overline{\vartheta}_4 \), decreases. As in Figure 1, it is also clear that the filtering error variances are greater than those of the smoothers and also that the smoothing estimators with lag \( L=3 \) outperform those with lag \( L=1 \).

Figure 2. Filtering and smoothing error variances for \( \overline{\vartheta}_4 = 0.3, 0.5, 0.7 \) and \( 0.9 \).

Effect of successful deception attacks probability \( \overline{\lambda} \). For \( \overline{\vartheta}_4=\overline{\gamma}=0.5 \), the filtering and smoothing error variances are analyzed as \( \overline{\lambda} \) varies from \( 0.1 \) to \( 0.9 \). The filtering error variances over all iterations are plotted in Figure 3(a), which clearly shows that the filter performance deteriorates as \( \overline{\lambda} \) increases (an analogous conclusion is drawn for the smoother performance). Figure 3(b), displays the estimation error variances at \( k = 50 \), versus \( \overline{\lambda} \). Despite showing only the results in a particular iteration, this figure provides better visibility of the increasing trend of the error variances with increasing probability \( \overline{\lambda} \).

Figure 3. (a) Filtering error variances for \( \overline{\lambda}=0.1 \) to 0.9; (b) Filtering and smoothing error variances at \( k = 50 \), versus \( \overline{\lambda} \).

Impact of transmission loss probability \( 1-\overline{\gamma} \). Considering \( \overline{\vartheta}_4= \overline{\lambda}=0.5 \), we analyze the influence of the loss probability, \( 1-\overline{\gamma} \), on the estimation accuracy, by comparing the filtering and smoothing error variances as \( \overline{\gamma} \) varies from \( 0.1 \) to \( 0.9 \). From Figure 4(a) we observe that the filtering error variances decrease as \( \overline{\gamma} \) increases; hence, better estimations are obtained when the transmission loss probability \( 1-\overline{\gamma} \) is low. This fact is better reflected in Figure 4(b), which presents the filtering and smoothing error variances at \( k = 50 \) versus \( \overline{\gamma} \).

Figure 4. (a) Filtering error variances for \( \overline{\gamma}=0.1 \) to 0.9; (b) Filtering and smoothing error variances at \( k = 50 \), versus \( \overline{\gamma} \).

Comparative analysis with CF filters in [24] and [29]. Finally, the designed CF filter is compared with the filters in [24] and [29]. From two thousand independent simulations, the comparison is made on the basis of the mean-squared error (MSE), whose empirical value at time \( k \) is calculated by:

where \( x_{k}^{(s)} \) and \( \widehat{x}_{k/k}^{(s)} \) denote the simulated state and filtering estimates, respectively, at time \( k \) for the \( s \)-th simulation run. Assuming again \( \overline{\vartheta}_4 =\overline{\lambda} = \overline{\gamma}=0.5 \), the results are depicted in Figure 5, which reveals that the MSE values of the proposed CF filtering estimates are consistently lower than those of the other two CF filtering estimates. This outcome was expected, as the proposed filter accounts for both RPD and RDA, while the filter in [29] does not account for random packet losses and the filter in [24] overlooks both random attacks and losses during data transmissions.

Figure 5. MSE comparison of the proposed CF filter and the filters in [24] and [29].

5. Conclusions

In this paper, we have designed recursive algorithms for the CF filtering and fixed-point smoothing problems in networked systems operating under the dual challenges of RDA and RPD, in the presence of different perturbations including multiplicative noise and time-correlated additive noise. The use of the process mean and covariance functions rather than the explicit state evolution model, simplifies the derivation of the algorithms without sacrificing accuracy. The impact of packet losses is mitigated by using ZOHS as a compensation methodology and time-correlated noises are not handled by measurement differencing, but by direct estimation of the noises.

Through numerical simulations, we have demonstrated that our fusion estimation scheme not only meets theoretical expectations but also outperforms some existing filters in practical scenarios. The superior performance of our approach underscores its potential for application in real-world networked systems where resilience to RDA and RPD is critical. The empirical results have also shown that, as the probability of deception attacks increases, the estimation accuracy decreases. It would be valuable to explore this inverse relationship theoretically in future research, by conducting a rigorous monotonicity analysis of the error covariance with respect to the attack probability. Future work will focus on extending this framework to more complex scenarios and further enhance its robustness and applicability. More specifically, new challenges and future research topics include:

● Extending the proposed study to deal with quantization effects to cover systems with limited network bandwidth, where it is common for sensor measurements to be quantized before transmission over the communication network.

● Considering different types of random attacks (replay attacks, DoS attacks, etc.).

● Incorporating energy-saving communication protocols to reduce communication burden and enhance transmission efficiency (event-triggered mechanisms, random access protocol, etc.).

Author Contributions: Raquel Caballero-Águila: Conceptualization, Methodology, Formal analysis, Investigation, Writing —original draft, Writing —review & editing, Software, Funding acquisition. Josefa Linares-Pérez: Conceptualization, Methodology, Formal analysis, Investigation, Writing —original draft, Writing—review & editing, Software, Funding acquisition.

Funding: Grant PID2021-124486NB-I00 funded by MICIU/AEI/ 10.13039/501100011033 and ERDF/EU.

Data Availability Statement: Not applicable.

Conflict of Interest: The authors declare no conflict of interest.