Quadratic Filtering Algorithm Based on Covariances Using Correlated Uncertain Observations Coming from Different Sensors

Abstract

The least-squares quadratic estimation problem of signals from observations coming from multiple sensors is addressed when there is a nonzero probability that each observation does not contain the signal to be estimated. We assume that, at each sensor, the uncertainty about the signal being present or missing in the observation is modelled by correlated Bernoulli random variables, whose probabilities are not necessarily the same for all the sensors. A recursive algorithm is derived without requiring the knowledge of the signal state-space model but only the moments (up to the fourth-order ones) of the signal and observation noise, the uncertainty probabilities, and the correlation between the variables modelling the uncertainty. The estimators require the autocovariance and cross-covariance functions of the signal and their second-order powers in a semidegenerate kernel form. The recursive quadratic filtering algorithm is derived from a linear estimation algorithm for a suitably defined augmented system.

1. Introduction

In many real systems the signal to be estimated can be randomly missing in the observations due, for example, to intermittent failures in the observation mechanism, fading phenomena in propagation channels, target tracking, accidental loss of some measurements, or data inaccessibility during certain times. Usually, these situations are characterized by including in the observation equation not only an additive noise, but also a multiplicative noise consisting of a sequence of Bernoulli random variables taking the value one if the observation is state plus noise, or the value zero if it is only noise (uncertain observations). Since these models are appropriate in many practical situations with random failures in the transmission, the estimation problem in systems with uncertain observations has been widely studied in the literature under different hypotheses and approaches (see e.g., [1, 2] and references therein).

On the other hand, in some practical situations the state-space model of the signal is not available and another type of information must be processed for the estimation. In the last years, the estimation problem from uncertain observations has been investigated using covariance information, and algorithms with a simpler structure than those obtained when the state-space model is known have been derived (see, e.g., [3]).

Recently, the least-squares linear estimation problem using uncertain observations transmitted by multiple sensors, whose statistical properties are assumed not to be the same, has been studied by several authors under different approaches and hypotheses on the processes (see, e.g., [4, 5] for a state-space approach and [6, 7] for a covariance approach).

In this paper, using covariance information, recursive algorithms for the least-squares quadratic filtering problem from correlated uncertain observations coming from multiple sensors with different uncertainty characteristics are proposed. This paper extends the results in [6] in two directions: on the one hand, correlation at times 𝑘 and 𝑘+𝑟 between the random variables modelling the uncertainty in the observations is considered, and, on the other, the quadratic estimation problem is addressed. The quadratic estimation is also a new topic addressed in this paper over the results in Hermoso-Carazo et al. [7] which are also referred to observations with uncertainty modelled by Bernoulli variables correlated at times 𝑘 and 𝑘+𝑟 with arbitrary 𝑟, but coming from a single sensor. Furthermore, the current paper differs from [5] in the correlation model considered and in the information used to derive the algorithms (state-space model in [5] and covariance information in the current paper).

To address the quadratic estimation problem, augmented signal and observation vectors are introduced by assembling the original vectors with their second-order powers defined by the Kronecker product, thus obtaining a new augmented system and reducing the quadratic estimation problem in the original system to the linear estimation problem in the augmented system. By using an innovation approach, the linear estimator of the augmented signal based on the augmented observations is obtained, thus providing the required quadratic estimator.

The performance of the proposed filtering algorithms is illustrated by a numerical simulation example where the state of a first-order autoregressive model is estimated from uncertain observations coming from two sensors with different uncertainty characteristics correlated at times 𝑘 and 𝑘+𝑟, considering several values of 𝑟. The linear and quadratic estimation error covariance matrices are compared, showing the superiority of the quadratic estimators over the linear ones.

2. Observation Model and Hypotheses

The problem at hand is to determine the least-squares (LS) quadratic estimator of an 𝑛-dimensional discrete signal, 𝑧_𝑘, from noisy measurements coming from multiple sensors which may not contain the signal with different probabilities. In this section, we present the observation model and the hypotheses about the signal and noise processes involved.

Consider 𝑚 scalar sensors whose measurements at each sampling time, 𝑘, denoted by 𝑦^𝑖𝑘, may either contain the signal to be estimate, 𝑧_𝑘, or be only noise, 𝑣^𝑖𝑘; the uncertainty about the signal being present or missing in the observation is modelled by Bernoulli variables, 𝛾^𝑖𝑘. The observation model is thus described as follows:𝑦^𝑖𝑘=𝛾^𝑖𝑘𝐻^𝑖𝑘𝑧_𝑘+𝑣^𝑖𝑘,𝑘≥1,𝑖=1,…,𝑚.(2.1) If 𝛾^𝑖𝑘=1, then 𝑦^𝑖𝑘=𝐻^𝑖𝑘𝑧_𝑘+𝑣^𝑖𝑘 and the measurement coming from the 𝑖th sensor contains the signal; otherwise, if 𝛾^𝑖𝑘=0, then 𝑦^𝑖𝑘=𝑣^𝑖𝑘, which means that such measurement is only noise. Therefore, the variables {𝛾^𝑖𝑘;𝑘≥1} model the uncertainty of the observations coming from the 𝑖th sensor.

To simplify the notation, the observation equation (2.1) is rewritten in a compact form as follows:𝑦_𝑘=Υ_𝑘𝐻_𝑘𝑧_𝑘+𝑣_𝑘,𝑘≥1,(2.2) where 𝑦_𝑘=(𝑦^1𝑘,…,𝑦^𝑚𝑘)^𝑇,𝐻_𝑘=(𝐻^𝑘1𝑇,…,𝐻^𝑘𝑚𝑇)^𝑇,Υ_𝑘=Diag(𝛾^1𝑘,…,𝛾^𝑚𝑘), and 𝑣_𝑘=(𝑣^1𝑘,…,𝑣^𝑚𝑘)^𝑇.

It is known that if the signal 𝑧_𝑘 and the observations 𝑦₁,…,𝑦_𝑘 have finite second-order moments, the LS linear filter of 𝑧_𝑘 is the orthogonal projection of 𝑧_𝑘 on the space of 𝑛-dimensional random variables obtained as linear transformations of 𝑦₁,…,𝑦_𝑘. So, by defining the random vectors 𝑦^𝑖[2]=𝑦_𝑖⊗𝑦_𝑖 (⊗ denotes the Kronecker product [8]) and, if 𝐸[𝑦^𝑖[2]𝑇𝑦^𝑖[2]]<∞, the LS quadratic estimator of 𝑧_𝑘 based on the observations up to the sampling time 𝑘 is the orthogonal projection of 𝑧_𝑘 on the space of 𝑛-dimensional linear transformations of 𝑦₁,…,𝑦_𝑘 and their second-order powers 𝑦^1[2],…,𝑦^𝑘[2]. To guarantee the existence of the second-order moments of the vectors 𝑦^𝑖[2], the pertinent assumptions about the processes in (2.1) are now stated.

2.1. Hypotheses about the Model

(H1)The 𝑛×1 signal process {𝑧_𝑘;𝑘≥1} has zero mean, and its autocovariance function, 𝐾^{𝑧𝑘,𝑠}, as well as the autocovariance function of the second-order powers, 𝐾^𝑧[2]_𝑘,𝑠, is expressed in a semidegenerate kernel form 𝐾^{𝑧𝑘,𝑠}=𝐴_𝑘𝐵^𝑇𝑠,𝑠≤𝑘,𝐾^𝑧[2]_𝑘,𝑠=𝑎_𝑘𝑏^𝑇𝑠,𝑠≤𝑘,(2.3) where the 𝑛×𝑀 matrix functions 𝐴, 𝐵, and the 𝑛²×𝐿 matrix functions 𝑎, 𝑏, are known. Moreover, it is assumed that the covariance function of the signal and its second-order powers, 𝐾^𝑧𝑧[2]_𝑘,𝑠, can also be expressed as 𝐾^𝑧𝑧[2]_𝑘,𝑠=𝛼^𝑘𝜀𝛽^𝑇𝑠,𝑠≤𝑘,_𝑘𝛿^𝑇𝑠,𝑘≤𝑠,(2.4) where 𝛼, 𝛽, 𝜀, and 𝛿 are 𝑛×𝑁, 𝑛²×𝑁, 𝑛×𝑃, and 𝑛²×𝑃 known matrix functions, respectively. (H2)For 𝑖=1,…,𝑚, the sensor additive noises, {𝑣^𝑖𝑘;𝑘≥1}, are zero mean white processes, and their moments, up to the fourth one, are known; we will denote 𝑅_𝑘=Cov[𝑣_𝑘], 𝑅^𝑘(3)=Cov[𝑣_𝑘,𝑣^𝑘[2]] and 𝑅^𝑘(4)=Cov[𝑣^𝑘[2]]. (H3)For 𝑖=1,…,𝑚, the noises {𝛾^𝑖𝑘;𝑘≥1} are sequences of Bernoulli random variables with 𝑃[𝛾^𝑖𝑘=1]=𝑝^𝑖𝑘; the variables 𝛾^𝑖𝑘 and 𝛾^𝑖𝑠 are independent for |𝑘−𝑠|≥2, and Cov[𝛾^𝑖𝑘,𝛾^𝑖𝑘+1] are assumed to be known. (H4)The signal process, {𝑧_𝑘;𝑘≥1}, and the noise processes, {𝛾_𝑘;𝑘>1} and {𝑣_𝑘;𝑘≥1}, where 𝛾_𝑘=(𝛾^1𝑘,…,𝛾^𝑚𝑘)^𝑇, are mutually independent.

3. Augmented System

Given the observation model (2.1) with assumptions (H1)–(H4), the problem is to find the LS quadratic filter of the signal 𝑧_𝑘, which will be denoted 𝑧^{𝑞𝑘/𝑘}. The technique used to obtain this estimator consists of augmenting the signal and observation by assembling the original vectors and their second-order powers, 𝒵_𝑘=(𝑧^𝑇𝑘,𝑧^𝑘[2]𝑇)^𝑇,𝒴_𝑘=(𝑦^𝑇𝑘,𝑦^𝑘[2]𝑇)^𝑇, and deriving the estimator 𝑧^{𝑞𝑘/𝑘} as the vector constituted by the first 𝑛 entries of the LS linear filter of 𝒵_𝑘 based on 𝒴_𝑘.

To obtain this linear estimator, the first- and second-order statistical properties of the augmented vectors 𝒵_𝑘 and 𝒴_𝑘 are now analyzed.

3.1. Properties of the Augmented Vectors

By using the Kronecker product properties and denoting 𝐷^𝛾𝑘=Diag(Υ_𝑘,Υ^𝑘[2]), ℋ_𝑘=Diag(𝐻_𝑘,𝐻^𝑘[2]) and 𝒱_𝑘=𝑣^𝑘𝐼_𝑚²+𝐾_𝑚2Υ_𝑘𝐻_𝑘𝑧_𝑘⊗𝑣_𝑘+𝑣^𝑘[2],(3.1) (𝐼_𝑚2 is the 𝑚²×𝑚² identity matrix and 𝐾_𝑚2 is the 𝑚²×𝑚² commutation matrix, [8]) the following model with uncertain observations is obtained, 𝒴_𝑘=𝐷^𝛾𝑘ℋ_𝑘𝒵_𝑘+𝒱_𝑘,𝑘≥1.(3.2) It should be noted that the signal, 𝒵_𝑘, and the noise, 𝒱_𝑘, in this new model have non-zero mean. Nevertheless, this handicap can be overcome by considering the centered augmented vectors 𝑍_𝑘=𝒵_𝑘−𝐸[𝒵_𝑘] and 𝑌_𝑘=𝒴_𝑘−𝐸[𝒴_𝑘] which, taking into account that 𝐸[𝐷^𝛾𝑘ℋ_𝑘𝑍_𝑘]=𝐸[𝐷^𝛾𝑘]ℋ_𝑘𝐸[𝑍_𝑘], satisfy𝑌_𝑘=𝐷^𝛾𝑘ℋ_𝑘𝑍_𝑘+𝑉_𝑘,𝑘≥1,(3.3) where 𝑉_𝑘=𝑣^𝑘𝐼_𝑚²+𝐾_𝑚2Υ_𝑘𝐻_𝑘𝑧_𝑘⊗𝑣_𝑘+𝑣^𝑘[2]−vec𝑅_𝑘+𝐷^𝛾𝑘−𝐷^𝑝𝑘ℋ_𝑘𝐸𝒵_𝑘,(3.4) being 𝐷^𝑝𝑘=𝐸[𝐷^𝛾𝑘] and vec the operator that vectorizes a matrix [8].

Note that the LS linear estimator of 𝒵_𝑘 based on 𝒴₁,…,𝒴_𝑘 is obtained from the LS linear estimator of 𝑍_𝑘 based on 𝑌₁,…,𝑌_𝑘, just adding the mean vector 𝐸[𝒵_𝑘]=(0^𝑇𝑛×1,(vec(𝐴_𝑘𝐵^𝑇𝑘))^𝑇)^𝑇. Hence, since the first 𝑛 components of 𝐸[𝒵_𝑘] are zero, the required quadratic estimator 𝑧^{𝑞𝑘/𝑘} is just the vector constituted by the first 𝑛 entries of the LS linear filter of 𝑍_𝑘. Henceforth, these centered vectors will be referred to as the augmented signal and observation vectors, respectively.

The signal and noise processes {𝑍_𝑘;𝑘≥1} and {𝑉_𝑘;𝑘≥1} involved in model (3.3) are zero mean. In the following propositions the second-order statistical properties of these processes are established.

Proposition 3.1. If the signal process {𝑧_𝑘;𝑘≥1} satisfies (H1), the autocovariance function of the augmented signal process {𝑍_𝑘;𝑘≥1} can be expressed in a semidegenerate kernel form, namely, 𝐾^{𝑍𝑘,𝑠}=𝐸𝑍_𝑘𝑍^𝑇𝑠=𝒜_𝑘ℬ^𝑇𝑠,𝑠≤𝑘,(3.5) where 𝒜_𝑘=𝐴^𝑘𝛼_𝑘0_𝑛×𝑃0_𝑛×𝐿ℬ_𝑘0_𝑛2×𝑀0_𝑛2×𝑁𝛿_𝑘𝑎_𝑘,=𝐵^𝑘0_𝑛×𝑁𝜀_𝑘0_𝑛×𝐿0_𝑛2×𝑀𝛽_𝑘0_𝑛2×𝑃𝑏_𝑘.(3.6)

Proof. It is immediate from hypothesis (H1) about the covariance functions of the signal and its second-order moments.

Proposition 3.2. Under (H1)–(H4), the noise {𝑉_𝑘;𝑘≥1} is a sequence of random vectors with covariance matrices, 𝐸[𝑉_𝑘𝑉^𝑇𝑠]=𝑅^{𝑉𝑘,𝑠}, given by 𝑅^{𝑉𝑘,𝑠}=⎧⎪⎪⎨⎪⎪⎩C𝑅o_𝑘v+Γ𝐶_𝑘𝑜𝑣Γ_𝑘∘ℋ_𝑘𝐸𝒵_𝑘𝐸𝒵^𝑇𝑘ℋ^𝑇𝑘,𝑠=𝑘0,,Γ_𝑘+1∘ℋ_𝑘𝐸𝒵_𝑘𝐸𝒵^𝑇𝑘+1ℋ^𝑇𝑘+1,𝑠||𝑘=−𝑘𝑠||+≠10,1,(3.7) where ∘ denotes the Hadamard product and Γ_𝑘=𝛾^𝑘𝛾^𝑘[2],𝑅_𝑘=𝑅^𝑘𝑅^𝑘(3)𝑅^𝑘(3)𝑇𝑅^𝑘22(3.8) with 𝑅^𝑘22=𝐼_𝑚2+𝐾_𝑚2𝐸𝛾_𝑘𝛾^𝑇𝑘∘𝐻_𝑘𝐴_𝑘𝐵^𝑇𝑘𝐻^𝑇𝑘⊗𝑅_𝑘𝐼_𝑚²+𝐾_𝑚2+𝑅^𝑘(4).(3.9) Moreover, {𝑉_𝑘;𝑘≥1} is uncorrelated with the processes {𝑍_𝑘;𝑘≥1} and {𝐷^𝛾𝑘ℋ_𝑘𝑍_𝑘;𝑘≥1}.

Proof. It is obvious that 𝐸[𝑉_𝑘]=0,forall𝑘≥1. On the other hand, since 𝑉_𝑘=𝒱_𝑘−𝐸[𝒱_𝑘]+(𝐷^𝛾𝑘−𝐷^𝑝𝑘)ℋ_𝑘𝐸[𝒵_𝑘], and {𝑧_𝑘;𝑘≥1}, {𝑣_𝑘;𝑘≥1} and {𝛾_𝑘;𝑘≥1} are mutually independent, it is easy to see that 𝐸[(𝒱_𝑘−𝐸[𝒱_𝑘])((𝐷^𝛾𝑠−𝐷^𝑝𝑠)ℋ_𝑠𝐸[𝒵_𝑠])^𝑇]=0,forall𝑘,𝑠 and hence 𝐸𝑉_𝑘𝑉^𝑇𝑠=Cov𝒱_𝑘,𝒱_𝑠+𝐸𝐷^𝛾𝑘−𝐷^𝑝𝑘ℋ_𝑘𝐸𝒵_𝑘𝐷^𝛾𝑠−𝐷^𝑝𝑠ℋ_𝑠𝐸𝒵_𝑠^𝑇.(3.10) Firstly, we prove that Cov𝒱_𝑘,𝒱_𝑠=𝑅^11𝑘,𝑠𝑅^12𝑘,𝑠𝑅^{12𝑇𝑘,𝑠}𝑅^22𝑘,𝑠=𝑅^𝑘𝑅^𝑘(3)𝑅^𝑘(3)𝑇𝑅^𝑘22𝛿^𝑘,𝑠,(3.11) where 𝛿 denotes the Kronecker delta function.
Indeed, since {𝑣_𝑘;𝑘≥1} is a zero mean white sequence with covariances 𝑅_𝑘,forall𝑘≥1, it is clear that 𝑅^11𝑘,𝑠=𝑅_𝑘𝛿_𝑘,𝑠. Moreover, from the mutual independence, the Kronecker and Hadamard products properties lead to 𝑅𝑅^12𝑘,𝑠_22𝑘,𝑠=Υ^𝑝𝑠𝐻_𝑠𝐸𝑧_𝑠^𝑇⊗𝑅_𝑘𝛿_𝑘,𝑠𝐼^𝑚2+𝐾_𝑚2+𝑅^𝑘(3)𝛿_𝑘,𝑠=𝑅^𝑘(3)𝛿_𝑘,𝑠,=𝐼_𝑚2+𝐾_𝑚2𝐸𝛾_𝑘𝛾^𝑇𝑠∘𝐻_𝑘𝐴_𝑘𝐵^𝑇𝑠𝐻^𝑇𝑠⊗𝑅_𝑘𝛿_𝑘,𝑠𝐼_𝑚²+𝐾_𝑚2+𝑅^𝑘(4)𝛿_𝑘,𝑠=𝑅^𝑘22𝛿_𝑘,𝑠.(3.12)
On the other hand, 𝐸𝐷^𝛾𝑘=−C𝐷ov^𝑝𝑘Γℋ_𝑘𝑘𝐸𝒵_𝑘𝐷^𝛾𝑠−𝐷^𝑝𝑠ℋ_𝑠𝐸𝒵_𝑠^𝑇,,Γ_𝑠∘ℋ_𝑘𝐸𝒵_𝑘𝐸𝒵^𝑇𝑠ℋ^𝑇𝑠(3.13) and since Cov[Γ_𝑘,Γ_𝑠]=0 for |𝑘−𝑠|≥2, the covariance matrices 𝑅^{𝑉𝑘,𝑠} are obtained.
The uncorrelation between {𝑉_𝑘;𝑘≥1} and the processes {𝑍_𝑘;𝑘≥1} and {𝐷^𝛾𝑘ℋ_𝑘𝑍_𝑘;𝑘≥1} is derived in a similar way, taking into account that {𝑧_𝑘;𝑘≥0}, {𝑣_𝑘;𝑘≥1}, and {𝛾_𝑘;𝑘≥1} are mutually independent and using the Kronecker and Hadamard products properties.

4. Quadratic Filtering Algorithm

As indicated above, to obtain the LS quadratic estimators of the signal 𝑧_𝑘 based on observations (2.1), we consider the LS linear estimators of the augmented signal, 𝑍_𝑘, based on the augmented observations (3.3). As known, the LS linear filter of 𝑍_𝑘 is the orthogonal projection of the vector 𝑍_𝑘 onto ℒ(𝑌₁,…,𝑌_𝑘), the linear space spanned by {𝑌₁,…,𝑌_𝑘}; so the Orthogonal Projection Lemma (OPL) states that the estimator, 𝑍_𝑘/𝑘, is the only linear combination satisfying the orthogonality property 𝐸𝑍_𝑘−𝑍_𝑘/𝑘𝑌^𝑇𝑠=0,𝑠=1,…,𝑘.(4.1)

Due to the fact that the observations are generally nonorthogonal vectors, we will use an innovation approach, consisting of transforming the observation process {𝑌_𝑘;𝑘≥1} into an equivalent process (innovation process) of orthogonal vectors {𝜈_𝑘;𝑘≥1}, equivalent in the sense that each set {𝜈₁,…,𝜈_𝑘} spans the same linear subspace as {𝑌₁,…,𝑌_𝑘}; that is, ℒ(𝜈₁,…,𝜈_𝑘)=ℒ(𝑌₁,…,𝑌_𝑘).

The innovation process is constructed by the Gram-Schmidt orthogonalization procedure, using an inductive reasoning. Starting with 𝜈₁=𝑌₁, the projection of the next observation, 𝑌₂, onto ℒ(𝜈₁) is given by 𝑌_2/1=𝐸[𝑌₂𝜈^𝑇1](𝐸[𝜈₁𝜈^𝑇1])⁻¹𝜈₁; then, the vector 𝜈₂=𝑌₂−𝑌_2/1 is orthogonal to 𝜈₁, and clearly ℒ(𝜈₁,𝜈₂)=ℒ(𝑌₁,𝑌₂). Let {𝜈₁,…,𝜈_𝑘−1} be the set of orthogonal vectors satisfying ℒ(𝜈₁,…,𝜈_𝑘−1)=ℒ(𝑌₁,…,𝑌_𝑘−1); if now we have an additional observation 𝑌_𝑘, we project it onto ℒ(𝜈₁,…,𝜈_𝑘−1) and the orthogonality property allows us to find the projection by separately projecting onto each of the previous orthogonal vectors; that is, 𝑌_{𝑘/𝑘−1}=^𝑘−1_𝑗=1𝐸𝑌_𝑘𝜈^𝑇𝑗𝐸𝜈_𝑗𝜈^𝑇𝑗⁻¹𝜈_𝑗;(4.2) so the next vector, 𝜈_𝑘=𝑌_𝑘−𝑌_{𝑘/𝑘−1}, is orthogonal to the previous ones and ℒ(𝜈₁,…,𝜈_𝑘)=ℒ(𝑌₁,…,𝑌_𝑘).

Note that the projection 𝑌_{𝑘/𝑘−1} is the part of the observation 𝑌_𝑘 that is determined by knowledge of {𝑌₁,…,𝑌_𝑘−1}; thus the remainder vector 𝜈_𝑘=𝑌_𝑘−𝑌_{𝑘/𝑘−1} can be regarded as the “new information” or the “innovation" provided by 𝑌_𝑘 and the process {𝜈_𝑘;𝑘≥1} as the innovation process associated with {𝑌_𝑘;𝑘≥1}. The causal and causally invertible linear relation existing between the observation and innovation processes makes the innovation process unique.

Next, taking into account that the innovations constitute a white process, we derive a general expression for the LS linear estimator of the augmented signal, 𝑍_𝑘, based on {𝑌₁,…,𝑌_𝐿}, which will be denoted by 𝑍_𝑘/𝐿. Replacing {𝑌₁,…,𝑌_𝐿} by the equivalent set of orthogonal vectors {𝜈₁,…,𝜈_𝐿}, the signal estimator is 𝑍_𝑘/𝐿=^𝐿_𝑗=1ℎ_𝑘,𝑗𝜈_𝑗,(4.3) where the impulse-response function ℎ_𝑘,𝑗,𝑗=1,…,𝐿 is calculated from the orthogonality property, 𝐸𝑍_𝑘−𝑍_𝑘/𝐿𝜈^𝑇𝑠=0,𝑠≤𝐿,(4.4) which leads to the Wiener-Hopf equation𝐸𝑍_𝑘𝜈^𝑇𝑠=^𝐿_𝑗=1ℎ_𝑘,𝑗𝐸𝜈_𝑗𝜈^𝑇𝑠,𝑠≤𝐿.(4.5) Due to the whiteness of the innovation process, 𝐸[𝜈_𝑗𝜈^𝑇𝑠]=0 for 𝑗≠𝑠 and the Wiener-Hopf equation is expressed as 𝐸𝑍_𝑘𝜈^𝑇𝑠=ℎ_𝑘,𝑠𝐸𝜈_𝑠𝜈^𝑇𝑠,𝑠≤𝐿;(4.6) consequently, ℎ_𝑘,𝑠=𝐸𝑍_𝑘𝜈^𝑇𝑠𝐸𝜈_𝑠𝜈^𝑇𝑠⁻¹,𝑠≤𝐿,(4.7) and, therefore, the following general expression for the LS linear filter of the augmented signal is obtained𝑍_𝑘/𝐿=^𝐿_𝑖=1𝑆_𝑘,𝑖Π^𝑖−1𝜈_𝑖,(4.8) where 𝑆_𝑘,𝑖=𝐸[𝑍_𝑘𝜈^𝑇𝑖] and Π_𝑖=𝐸[𝜈_𝑖𝜈^𝑇𝑖].

Using the properties of the processes involved in (3.3), as established in Propositions 3.1 and 3.2, and expression (4.8) for the filter, we derive a recursive algorithm for the linear filtering estimators, 𝑍_𝑘/𝑘, of the augmented signal 𝑍_𝑘. As indicated above, the first 𝑛 entries of these estimators provide the required quadratic filter of the original signal 𝑧_𝑘.

Theorem 4.1. The quadratic filter, 𝑧^{𝑞𝑘/𝑘}, of the original signal 𝑧_𝑘 is given by 𝑧^{𝑞𝑘/𝑘}=Θ𝑍_𝑘/𝑘,𝑘≥1,(4.9) where Θ is the operator which extracts the first 𝑛 entries of 𝑍_𝑘/𝑘, the linear filter of the augmented signal 𝑍_𝑘, which is obtained by 𝑍_𝑘/𝑘=𝒜_𝑘𝑂_𝑘,𝑘≥1,(4.10) where the vectors 𝑂_𝑘 are recursively calculated from 𝑂_𝑘=𝑂_𝑘−1+𝐽_𝑘Π^𝑘−1𝜈_𝑘,𝑘≥1;𝑂₀=0.(4.11) The innovation, 𝜈_𝑘, satisfies 𝜈_𝑘=𝑌_𝑘−𝐷^𝑝𝑘ℋ_𝑘𝒜_𝑘𝑂_𝑘−1−Ξ_{𝑘,𝑘−1}𝜈_𝑘−1,𝑘≥2;𝜈₁=𝑌₁,(4.12) with Ξ_{𝑘,𝑘−1}=CovΓ_𝑘,Γ_𝑘−1∘ℋ_𝑘𝒜_𝑘ℬ^𝑇𝑘−1ℋ^𝑇𝑘−1+𝑅^{𝑉𝑘,𝑘−1}Π^−1𝑘−1,𝑘≥2,(4.13) and Π_𝑘, the covariance matrix of the innovation, verifies Π_𝑘=𝐸Γ_𝑘Γ^𝑇𝑘∘ℋ_𝑘𝒜_𝑘ℬ^𝑇𝑘ℋ^𝑇𝑘−𝐷^𝑝𝑘ℋ_𝑘𝒜_𝑘𝑟_𝑘−1𝒜^𝑇𝑘ℋ^𝑇𝑘𝐷^p_𝑘−Ξ_{𝑘,𝑘−1}Π_𝑘−1Ξ^{𝑇𝑘,𝑘−1}Π−𝐷^𝑝𝑘ℋ_𝑘𝒜_𝑘𝐽_𝑘−1Ξ^{𝑇𝑘,𝑘−1}−Ξ_{𝑘,𝑘−1}𝐽^𝑇𝑘−1𝒜^𝑇𝑘ℋ^𝑇𝑘𝐷^𝑝𝑘+𝑅^{𝑉𝑘,𝑘},𝑘≥2,₁=𝐸Γ₁Γ^𝑇1∘ℋ₁𝒜₁ℬ^𝑇1ℋ^𝑇1+𝑅^𝑉1,1.(4.14) The matrix function 𝐽 is given by 𝐽𝐽_𝑘=ℬ^𝑇𝑘−𝑟_𝑘−1𝒜^𝑇𝑘ℋ^𝑇𝑘𝐷^𝑝𝑘−𝐽_𝑘−1Ξ^{𝑇𝑘,𝑘−1},𝑘≥2,₁=ℬ^𝑇1ℋ^𝑇1𝐷^𝑝1,(4.15) where 𝑟_𝑘 is recursively obtained from 𝑟_𝑘=𝑟_𝑘−1+𝐽_𝑘Π^𝑘−1𝐽^𝑇𝑘,𝑘≥1;𝑟₀=0.(4.16)

Proof. We start by obtaining an explicit formula for the innovations, 𝜈_𝑘=𝑌_𝑘−𝑌_{𝑘/𝑘−1}, or, equivalently, for the one-stage predictor of 𝑌_𝑘, which by denoting 𝑇_𝑘,𝑖=𝐸[𝑌_𝑘𝜈^𝑇𝑖] is given by 𝑌_{𝑘/𝑘−1}=^𝑘−1_𝑖=1𝑇_𝑘,𝑖Π^𝑖−1𝜈_𝑖,𝑘≥2;𝑌_1/0=0.(4.17) Using the hypotheses on the model, it is deduced that 𝑇_𝑘,𝑖=𝐷^𝑝𝑘ℋ_𝑘𝑆_𝑘,𝑖,𝑖<𝑘−1,(4.18) and hence, 𝑌_{𝑘/𝑘−1}=𝐷^𝑝𝑘ℋ_𝑘𝑍_{𝑘/𝑘−1}+𝑇_{𝑘,𝑘−1}−𝐷^𝑝𝑘ℋ_𝑘𝑆_{𝑘,𝑘−1}Π^−1𝑘−1𝜈_𝑘−1.(4.19) Using again the hypotheses on the model, we obtain 𝑇_{𝑘,𝑘−1}−𝐷^𝑝𝑘ℋ_𝑘𝑆_{𝑘,𝑘−1}=CovΓ_𝑘,Γ_𝑘−1∘ℋ_𝑘𝒜_𝑘ℬ^𝑇𝑘−1ℋ^𝑇𝑘−1+𝑅^{𝑉𝑘,𝑘−1};(4.20) consequently, 𝜈_𝑘=𝑌_𝑘−𝐷^𝑝𝑘ℋ_𝑘𝑍_{𝑘/𝑘−1}−Ξ_{𝑘,𝑘−1}𝜈_𝑘−1,𝑘≥2;𝜈₁=𝑌₁,(4.21) with Ξ_{𝑘,𝑘−1} given by expression (4.13).
Next, expression (4.10) for the filter 𝑍_𝑘/𝑘 is derived. For this purpose, taking into account expression (4.8), we obtain formulas to calculate the coefficients 𝑆_𝑘,𝑖=𝐸[𝑍_𝑘𝜈^𝑇𝑖],𝑖≤𝑘. From the hypotheses on the model, replacing 𝜈_𝑖 by its expression in (4.21), and using (4.8) for 𝑍_{𝑖/𝑖−1}, we have 𝑆_𝑘,𝑖=𝒜_𝑘ℬ^𝑇𝑖ℋ^𝑇𝑖𝐷^𝑝𝑖−^𝑖−1_𝑗=1𝑆𝑆_𝑘,𝑗Π^𝑗−1𝑆^{𝑇𝑖,𝑗}ℋ^𝑇𝑖𝐷^𝑝𝑖−𝑆_{𝑘,𝑖−1}Ξ^{𝑇𝑖,𝑖−1},2≤𝑖≤𝑘;_𝑘,1=𝒜₁ℬ^𝑇1ℋ^𝑇1𝐷^𝑝1(4.22) or, equivalently, 𝑆_𝑘,𝑖=𝐴_𝑘𝐽_𝑖,𝑖≤𝑘,(4.23) where 𝐽 is a function satisfying 𝐽_𝑖=ℬ^𝑇𝑖ℋ^𝑇𝑖𝐷^𝑝𝑖−^𝑖−1_𝑗=1𝐽𝐽_𝑗Π^𝑗−1𝑆^{𝑇𝑖,𝑗}ℋ^𝑇𝑖𝐷^𝑝𝑖−𝐽_𝑖−1Ξ^{𝑇𝑖,𝑖−1},2≤𝑖≤𝑘;₁=ℬ^𝑇1ℋ^𝑇1𝐷^𝑝1.(4.24) Then, from (4.8) and (4.23), expression (4.10) for the filter is deduced, where 𝑂_𝑘, defined by 𝑂_𝑘=^𝑘_𝑖=1𝐽_𝑖Π^𝑖−1𝜈_𝑖,𝑘≥1;𝑂₀=0(4.25) satisfies the recursive relation (4.11). Analogously, it is obtained that the one-stage predictor of the signal is given by 𝑍_{𝑘/𝑘−1}=𝒜_𝑘𝑂_𝑘−1, which substituted into (4.21) leads to formula (4.12) for the innovation.
Expression (4.15) for 𝐽_𝑘 is derived making 𝑖=𝑘 in (4.24), using (4.23), and defining the function 𝑟_𝑘=𝐸[𝑂_𝑘𝑂^𝑇𝑘]=∑^𝑘𝑗=1𝐽_𝑗Π^𝑗−1𝐽^𝑇𝑗,𝑘≥1;𝑟₀=0. From this definition, the recursive relation (4.16) is also immediately derived.
Finally, we obtain expression (4.14) for the innovation covariance matrix; from the hypotheses on the model, expression (4.12), and the definition of 𝑟_𝑘, the following equation is obtained: Π_𝑘=𝐸Γ_𝑘Γ^𝑇𝑘∘ℋ_𝑘𝒜_𝑘ℬ^𝑇𝑘ℋ^𝑇𝑘+𝑅^{𝑉𝑘,𝑘}−𝐷^𝑝𝑘ℋ_𝑘𝒜_𝑘𝑟_𝑘−1𝒜^𝑇𝑘ℋ^𝑇𝑘𝐷^𝑝𝑘−Ξ_{𝑘,𝑘−1}Π_𝑘−1Ξ^{𝑇𝑘,𝑘−1}−𝐷^𝑝𝑘ℋ_𝑘𝒜_𝑘𝐸𝑂_𝑘−1𝜈^𝑇𝑘−1Ξ^{𝑇𝑘,𝑘−1}−Ξ_{𝑘,𝑘−1}𝐸𝜈_𝑘−1𝑂^𝑇𝑘−1𝒜^𝑇𝑘ℋ^𝑇𝑘𝐷^𝑝𝑘,𝑘≥2.(4.26) So, expression (4.14) for Π_𝑘 is deduced taking into account that 𝐸[𝑂_𝑘−1𝜈^𝑇𝑘−1]=𝐽_𝑘−1, which follows from (4.11) using that the vector 𝑂_𝑘−2 is orthogonal to 𝜈_𝑘−1.

To conclude, as a measure of the estimation accuracy, we have calculated the filtering error covariance matrices, Σ_𝑘/𝑘=𝐸[𝑍_𝑘𝑍^𝑇𝑘]−𝐸[𝑍_𝑘/𝑘𝑍^{𝑇𝑘/𝑘}], which clearly are obtained by Σ_𝑘/𝑘=𝒜_𝑘[ℬ^𝑇𝑘−𝑟_𝑘𝒜^𝑇𝑘],𝑘≥1.

5. Generalization to Correlation at Times 𝑘 and 𝑠, with |𝑘−𝑠|=𝑟

The observation model considered in Section 2 assumes that the uncertainty is modeled by Bernoulli variables correlated at consecutive sampling times, but independent otherwise. In this section, such model is generalized by assuming correlation in the uncertainty at times 𝑘 and 𝑠 differing 𝑟 units of time. Specifically, hypothesis (H3) is replaced by the following one.(H3')For 𝑖=1,…,𝑚, the noises {𝛾^𝑖𝑘;𝑘≥1} are sequences of Bernoulli random variables with 𝑃[𝛾^𝑖𝑘=1]=𝑝^𝑖𝑘. For 𝑖,𝑗=1,…,𝑚, the variables 𝛾^𝑖𝑘 and 𝛾^𝑗𝑠 are assumed to be independent for |𝑘−𝑠|≠𝑟, and Cov[𝛾^𝑖𝑘,𝛾^𝑗𝑠] are known for |𝑘−𝑠|=𝑟.

This correlation model allows us to consider certain situations where the signal cannot be missing at 𝑟+1 consecutive observations.

Similar considerations to those made in Section 3 for the case of consecutive sampling times lead now to the following expression for the covariance matrices of the noise {𝑉_𝑘;𝑘≥1}: 𝑅^{𝑉𝑘,𝑠}=⎧⎪⎪⎨⎪⎪⎩𝑅Co_𝑘v+ΓC_𝑘ovΓ_𝑘∘ℋ_𝑘𝐸𝒵_𝑘𝐸𝒵^𝑇𝑘ℋ^𝑇𝑘,𝑠=𝑘,0,,Γ_𝑠∘ℋ_𝑘𝐸𝒵_𝑘𝐸𝒵^𝑇𝑠ℋ^𝑇𝑠,||||𝑘𝑘−−𝑠𝑠||||=≠0𝑟,,𝑟.(5.1) Then, performing the same steps as in the proof of Theorem 4.1, the following algorithm is deduced.

Theorem 5.1. The quadratic filter, 𝑧^{𝑞𝑘/𝑘}, of the original signal 𝑧_𝑘 is given by 𝑧^{𝑞𝑘/𝑘}=Θ𝑍_𝑘/𝑘,𝑘≥1,(5.2) where Θ is the operator which extracts the first 𝑛 entries of 𝑍_𝑘/𝑘, the linear filter of the augmented signal 𝑍_𝑘, which is obtained by 𝑍_𝑘/𝑘=𝒜_𝑘𝑂_𝑘,𝑘≥1,(5.3) where the vectors 𝑂_𝑘 are recursively calculated from 𝑂_𝑘=𝑂_𝑘−1+𝐽_𝑘Π^𝑘−1𝜈_𝑘,𝑘≥1;𝑂₀=0.(5.4) The innovation, 𝜈_𝑘, satisfies 𝜈𝜈_1𝑘=𝑌₁,𝜈=𝑌_𝑘−𝐷^𝑝𝑘ℋ_𝑘𝒜_𝑘𝑂_𝑘−1,2≤𝑘≤𝑟,_𝑘=𝑌_𝑘−𝐷^𝑝𝑘ℋ_𝑘𝒜_𝑘𝑂_𝑘−1−Ξ_{𝑘,𝑘−𝑟}𝜈^𝑘−𝑟−^𝑘−1_𝑖=_{𝑘−𝑟+1}𝑇^{𝑇𝑖,𝑘−𝑟}Π^𝑖−1𝜈_𝑖,𝑘>𝑟(5.5) with Ξ_{𝑘,𝑘−𝑟}=CovΓ_𝑘,Γ_𝑘−𝑟∘ℋ_𝑘𝒜_𝑘ℬ^{𝑇𝑘−𝑟}ℋ^{𝑇𝑘−𝑟}+𝑅^{𝑉𝑘,𝑘−𝑟}Π^{−1𝑘−𝑟}𝑇_{𝑖,𝑘−𝑟}=𝐷^𝑝𝑖ℋ_𝑖𝑆_{𝑖,𝑘−𝑟}−Ξ_{𝑖,𝑖−𝑟}𝑇^{𝑇𝑘−𝑟,𝑖−𝑟},𝑖=𝑘−𝑟+1,…,𝑘−1.(5.6) The covariance matrix of the innovation, Π_𝑘, verifies ΠΠ_1𝑘=𝐸Γ₁Γ^𝑇1∘ℋ₁𝒜₁ℬ^𝑇1ℋ^𝑇1+𝑅^𝑉1,1,Π=CovΓ_𝑘∘ℋ_𝑘𝒜_𝑘ℬ^𝑇𝑘ℋ^𝑇𝑘+𝑅^{𝑉𝑘,𝑘}+𝐷^𝑝𝑘ℋ_𝑘𝒜_𝑘𝐽_𝑘,2≤𝑘≤𝑟,_𝑘=CovΓ_𝑘∘ℋ_𝑘𝒜_𝑘ℬ^𝑇𝑘ℋ^𝑇𝑘+𝑅^{𝑉𝑘,𝑘}+𝐷^𝑝𝑘ℋ_𝑘𝒜_𝑘𝐽_𝑘−Ξ_{𝑘,𝑘−𝑟}Π^𝑘−𝑟+^𝑘−1_𝑖=_{𝑘−𝑟+1}𝑇^{𝑇𝑖,𝑘−𝑟}Π^𝑖−1𝑇_{𝑖,𝑘−𝑟}Ξ^{𝑇𝑘,𝑘−𝑟}−𝐷^𝑝𝑘ℋ_𝑘ℬ_𝑘−𝒜_𝑘𝑟_𝑘−1−𝐽^𝑇𝑘𝒜^𝑇𝑘ℋ^𝑇𝑘𝐷^𝑝𝑘,𝑘>𝑟.(5.7) The matrix function 𝐽 is given by 𝐽𝐽_1𝑘=𝑝₁ℬ^𝑇1ℋ^𝑇1,𝐽=ℬ^𝑇𝑘−𝑟_𝑘−1𝒜^𝑇𝑘ℋ^𝑇𝑘𝐷^𝑝𝑘,2≤𝑘≤𝑟,_𝑘=ℬ^𝑇𝑘−𝑟_𝑘−1𝒜^𝑇𝑘ℋ^𝑇𝑘𝐷^𝑝𝑘−𝐽_𝑘−𝑟−^𝑘−1_𝑗=_{𝑘−𝑟+1}𝐽_𝑗Π^𝑗−1𝑇_{𝑗,𝑘−𝑟}Ξ^{𝑇𝑘,𝑘−𝑟},𝑘>𝑟,(5.8) where 𝑟_𝑘 is recursively obtained from 𝑟_𝑘=𝑟_𝑘−1+𝐽_𝑘Π^𝑘−1𝐽^𝑇𝑘,𝑘≥1;𝑟₀=0.(5.9)

Proof. It is analogous to that of Theorem 4.1, taking into account that, in this case, the one-stage predictor of 𝑌_𝑘 satisfies 𝑌_{𝑘/𝑘−1}=𝐷^𝑝𝑘ℋ_𝑘𝑍_{𝑘/𝑘−1}+^𝑘−1_𝑖=_𝑘−𝑟𝑇_𝑘,𝑖−𝐷^𝑝𝑘ℋ_𝑘𝑆_𝑘,𝑖Π^𝑖−1𝜈_𝑖,(5.10) and, from the model hypotheses, 𝑇_{𝑘,𝑘−𝑟}−𝐷^𝑝𝑘ℋ_𝑘𝑆_{𝑘,𝑘−𝑟}=CovΓ_𝑘,Γ_𝑘−𝑟∘ℋ_𝑘𝒜_𝑘ℬ^{𝑇𝑘−𝑟}ℋ^{𝑇𝑘−𝑟}+𝑅^{𝑉𝑘,𝑘−𝑟},(5.11) and for 𝑖>𝑘−𝑟, 𝑇_𝑘,𝑖−𝐷^𝑝𝑘ℋ_𝑘𝑆_𝑘,𝑖=−CovΓ_𝑘,Γ_𝑘−𝑟∘ℋ_𝑘𝒜_𝑘ℬ^{𝑇𝑘−𝑟}ℋ^{𝑇𝑘−𝑟}+𝑅^{𝑉𝑘,𝑘−𝑟}Π^{−1𝑘−𝑟}𝑇^{𝑇𝑖,𝑘−𝑟}.(5.12)

6. Numerical Simulation Example

To illustrate the application of the proposed filtering algorithm a numerical simulation example is shown now. To check the effectiveness of the proposed quadratic filter, we ran a program in MATLAB which, at each iteration, simulates the signal and the observed values and provides the linear and quadratic filtering estimates, as well as the corresponding error covariance matrices.

For the simulations, this program has been applied to a scalar signal {𝑧_𝑘;𝑘≥1}, generated by the following first-order autoregressive model: 𝑧_𝑘=0.95𝑧_𝑘−1+𝑤_𝑘−1,𝑘≥1,(6.1) where the initial state is a zero mean Gaussian variable with Var[𝑧₀]=1 and {𝑤_𝑘;𝑘≥0} is a zero mean white Gaussian noise with Var[𝑤_𝑘]=1.

The autocovariance functions of the signal and their second-order powers are given in a semidegenerate kernel form, specifically, 𝐾^{𝑧𝑘,𝑠}=1.025641×0.95^𝑘−𝑠,𝐾^𝑧2_𝑘,𝑠=2.1038795×0.95^{2(𝑘−𝑠)},𝑠≤𝑘(6.2)

and the covariance function of the signal and their second-order powers is given by 𝐾^𝑧𝑧2_𝑘,𝑠=0,forall𝑠,𝑘. According to hypothesis (H1), the functions constituting these covariance functions can be defined as follows:𝐴_𝑘=1.025641×0.95^𝑘,𝐵_𝑘=𝛼_𝑘0.95^−𝑘,𝑎_𝑘=2.1038795×0.95^2𝑘,𝑏_𝑘=0.95^−2𝑘,=𝛽_𝑘=𝜀_𝑘=𝛿_𝑘=0.(6.3)

Consider two sensors whose measurements, according to our theoretical study, are perturbed by sequences of Bernoulli random variables {𝛾^𝑘(𝑖);𝑘≥1}, 𝑖=1,2, and by additive white noises, {𝑣^𝑘(𝑖);𝑘≥1}, 𝑖=1,2; that is, 𝑦^𝑖𝑘=𝛾^𝑘(𝑖)𝑧_𝑘+𝑣^𝑘(𝑖),𝑘≥1,𝑖=1,2.(6.4)

Assume that the additive noises are independent and have the following probability distributions: 𝑃𝑣^𝑘(2)=−8=18,𝑃𝑣^𝑘(1)=8𝑃𝑣7=78,∀𝑘≥1,(6.5)_𝑘(2)=1=1518,𝑃𝑣^𝑘(2)=−3=218,𝑃𝑣^𝑘(2)=−9=118,∀𝑘≥1.(6.6)

Now, in accordance with the proposed uncertain observation model, we assume that the uncertainty at any time 𝑘>𝑟 is correlated with the uncertainty at time 𝑠 only if |𝑘−𝑠|=𝑟, but independent otherwise.

To model the uncertainty in this way, we can consider two independent sequences of independent Bernoulli random variables, {𝜃^𝑘(𝑖);𝑘≥1}, 𝑖=1,2, and define 𝛾^𝑘(𝑖)=1−𝜃^𝑘(𝑖)(1−𝜃^{(𝑖)𝑘+𝑟}),𝑘≥1. It is obvious that the variables 𝛾^𝑘(𝑖) take the value zero if and only if 𝜃^𝑘(𝑖)=1 and 𝜃^{(𝑖)𝑘+𝑟}=0; otherwise, 𝛾^𝑘(𝑖)=1; therefore, {𝛾^𝑘(𝑖);𝑘≥1} are Bernoulli variables with 𝑃[𝛾^𝑘(𝑖)=1]=1−𝑃[𝜃^𝑘(𝑖)=1]𝑃[𝜃^{(𝑖)𝑘+𝑟}=0]. Note that 𝛾^𝑘(𝑖)=0 implies 𝜃^{(𝑖)𝑘+𝑟}=0 and, consequently, 𝛾^{(𝑖)𝑘+𝑟}=1; this fact implies that if the signal is missing at time 𝑘, it is assured that, at time 𝑘+𝑟, the observation contains the signal; therefore, the signal cannot be missing in 𝑟+1 consecutive observations.

For the application, we have assumed that the variables 𝜃^𝑘(𝑖) in each sensor have the same distribution; that is, 𝑃[𝜃^𝑘(𝑖)=1]=𝜃_𝑖 independent of 𝑘. So, in each sensor, the probability that the observation contains the signal, 𝑝_𝑖=𝑃[𝛾^𝑘(𝑖)=1]=1−𝜃_𝑖(1−𝜃_𝑖), is constant for all the observations.

Since {𝜃^𝑘(1);𝑘≥1} is independent of {𝜃^𝑘(2);𝑘≥1}, also {𝛾^𝑘(1);𝑘≥1} is independent of {𝛾^𝑘(2);𝑘≥1} and hence, forall𝑘,𝑠, Cov[𝛾^𝑘(𝑖),𝛾^𝑠(𝑗)]=0 for 𝑖≠𝑗. For fixed 𝑖=1,2, the variance of 𝛾^𝑘(𝑖) is Var[𝛾^𝑘(𝑖)]=𝑝_𝑖(1−𝑝_𝑖), and the correlation between two different variables 𝛾^𝑘(𝑖) and 𝛾^𝑠(𝑖) is obtained as follows.(I)For |𝑘−𝑠|≠0,𝑟 the variables (𝜃^𝑘(𝑖),𝜃^{(𝑖)𝑘+𝑟}) and (𝜃^𝑠(𝑖),𝜃^{(𝑖)𝑠+𝑟}) are independent, thus being the variables 𝛾^𝑘(𝑖)=1−𝜃^𝑘(𝑖)(1−𝜃^{(𝑖)𝑘+𝑟}) and 𝛾^𝑠(𝑖)=1−𝜃^𝑠(𝑖)(1−𝜃^{(𝑖)𝑠+𝑟}) also independent and hence, uncorrelated; that is, Cov[𝛾^𝑘(𝑖),𝛾^𝑠(𝑖)]=𝐸[𝛾^𝑘(𝑖)𝛾^𝑠(𝑖)]−𝐸[𝛾^𝑘(𝑖)]𝐸[𝛾^𝑠(𝑖)]=0.(II)For |𝑘−𝑠|=𝑟, 𝐸[𝜃^𝑘(𝑖)(1−𝜃^{(𝑖)𝑘+𝑟})𝜃^𝑠(𝑖)(1−𝜃^{(𝑖)𝑠+𝑟})]=0 since 𝑘=𝑠+𝑟 or 𝑠=𝑘+𝑟, and 𝜃^(𝑖)(1−𝜃^(𝑖))=0. Then, 𝐸[𝛾^𝑘(𝑖)𝛾^𝑠(𝑖)]=1−𝐸[𝜃^𝑘(𝑖)(1−𝜃^{(𝑖)𝑘+𝑟})]−𝐸[𝜃^𝑠(𝑖)(1−𝜃^{(𝑖)𝑠+𝑟})]=1−2𝜃_𝑖(1−𝜃_𝑖)=2𝑝_𝑖−1, which implies that Cov[𝛾^𝑘(𝑖),𝛾^𝑠(𝑖)]=𝐸[𝛾^𝑘(𝑖)𝛾^𝑠(𝑖)]−𝐸[𝛾^𝑘(𝑖)]𝐸[𝛾^𝑠(𝑖)]=2𝑝_𝑖−1−𝑝^2𝑖=−(1−𝑝_𝑖)².

Summarizing, the correlation function of {𝛾^𝑘(𝑖);𝑘≥1} is given by Cov𝛾^𝑘(𝑖),𝛾^𝑠(𝑖)=⎧⎪⎪⎨⎪⎪⎩𝑝0−,11−−𝑝𝑝_𝑖𝑖,||||𝑘𝑘−−𝑠𝑠||||=≠𝑟0,^𝑖,²,||𝑘−𝑠||=𝑟,(6.7)

and, hence, the measurements above described are in accordance with the proposed correlation model.

To analyze the performance of the proposed estimators, the linear and quadratic filtering error variances have been calculated for different values of 𝑟 and also for different 𝜃₁ and 𝜃₂, which provide different values of the probabilities 𝑝₁ and 𝑝₂. Since 𝑝_𝑖 are the same if the value 1−𝜃_𝑖 is considered instead of 𝜃_𝑖, only the case 𝜃_𝑖≤0.5 is examined here (note that, in such case, 𝑝_𝑖 is a decreasing function of 𝜃_𝑖); more specifically, the values 𝜃_𝑖=0.1,0.2,0.3,0.4,0.5 (which lead to 𝑝_𝑖=0.91,0.84,0.78,0.76,0.75, resp.) have been used.

First, considering 𝑟=4, the linear and quadratic filtering error variances are calculated for the values 𝜃₁=𝜃₂=0.1, 𝜃₁=0.2 and 𝜃₂=0.3, 𝜃₁=0.4 and 𝜃₂=0.5. Figure 1 shows the results obtained; for all the values of 𝜃_𝑖, the error variances corresponding to the quadratic filter are always considerably less than those of the linear filter, thus confirming the superiority of the quadratic filter over the linear one in the estimation accuracy. Also, from this figure it is gathered that, as 𝜃₁ or 𝜃₂ increase (which means that the probability that the signal is present in the observations coming from the corresponding sensor decreases), the filtering error variances become greater and, hence, worse estimations are obtained.

Figure 1 

Linear and quadratic filtering error variances for 𝑟=4 and different values of 𝜃1 and 𝜃2.

Next, we compare the performance of the linear and quadratic filtering estimators for the values 𝜃_𝑖=0.1,0.2,0.3,0.4,0.5; since the linear and quadratic filtering error variances show insignificant variation from the 5th iteration onwards only the error variances at a specific iteration are considered.

In Figure 2 the linear and quadratic filtering error variances at 𝑘=50 are displayed versus 𝜃₁ (for constant values of 𝜃₂), and, in Figure 3, these variances are shown versus 𝜃₂ (for constant values of 𝜃₁). From these figures it is gathered that, as 𝜃₁ or 𝜃₂ decrease (and, consequently, the probability that the signal is not present in the observations coming from the corresponding sensor, 1−𝑝_𝑖, decreases), the filtering error variances become smaller and, hence, better estimations are obtained. Note that this improvement is more significant for small values of 𝜃₁ or 𝜃₂, that is, when the probability that the signal is present in the observations coming from one of the sensors is large. On the other hand, both figures show that, for all the values of 𝜃₁ and 𝜃₂, the error variances corresponding to the quadratic filter are always considerably less than those of the linear filter, confirming again the superiority of the quadratic filter over the linear one.

(a)

(b)

Figure 2 

Linear and quadratic filtering error variances at 𝑘=50 versus 𝜃1, for constant values of 𝜃2=0.1,0.2,0.3,0.4,0.5.

(a)

(b)

Figure 3 

Linear and quadratic filtering error variances at 𝑘=50 versus 𝜃2, for constant values of 𝜃2=0.1,0.2,0.3,0.4,0.5.

Finally, for 𝜃₁=𝜃₂=0.5 (these values produce the maximum value for the probability that the signal is not present in the observations coming from both sensors) and considering different values of 𝑟, specifically 𝑟=1,…,16, the error variances at 𝑘=50 for the linear and quadratic filters are displayed in Figure 4. From this figure it is deduced that the performance of the estimators improves when the values of 𝑟 are smaller and, hence, a greater distance between the correlated variables produces worse estimations (in the sense of the mean squared error). As expected, this figure also shows that the estimation accuracy of the quadratic filters is superior to that of the linear filters and also that the error variances show insignificant variation when the values of 𝑟 are greater.

(a)

(b)

Figure 4 

Linear and quadratic filtering error variances at 𝑘=50 versus 𝑟, when 𝜃1=𝜃2=0.5.

7. Conclusion

A recursive quadratic filtering algorithm is proposed from correlated uncertain observations coming from multiple sensors with different uncertainty characteristics. This is a realistic assumption in situations concerning sensor data that are transmitted over communication networks where, generally, multiple sensors with different properties are involved. The uncertainty in each sensor is modelled by a sequence of Bernoulli random variables which are correlated at times 𝑘 and 𝑘+𝑟. A real application of such observation model arises, for example, in signal transmission problems where a failure in one of the sensors at time 𝑘 is detected and the old sensor is replaced at time 𝑘+𝑟, thus avoiding the possibility of missing signal in 𝑟+1 consecutive observations.

Using covariance information, the algorithm is derived by applying the innovation technique to suitably defined augmented signal and observation vectors, and the LS quadratic estimator of the signal is obtained from the LS linear estimator of the augmented signal based on the augmented observations.

The performance of the proposed filtering algorithm is illustrated by a numerical simulation example where the state of a first-order autoregressive model is estimated from uncertain observations coming from two sensors with different uncertainty characteristics correlated at times 𝑘 and 𝑘+𝑟, considering several values of 𝑟.

Acknowledgment

This paper is supported by Ministerio de Educación y Ciencia (Grant no. MTM2008-05567) and Junta de Andalucía (Grant no. P07-FQM-02701).