# System Model Bias Processing Approach for Regional Coordinated States Information Involved Filtering

Authors: Zebo Zhou; Yunlong Wu
Journal: Mathematical Problems in Engineering
 (2016)

---

## Abstract

In the Kalman filtering applications, the conventional dynamic model which connects the states information of two consecutive epochs by state transition matrix is usually predefined and assumed to be invariant. Aiming to improve the adaptability and accuracy of dynamic model, we propose multiple historical states involved filtering algorithm. An autoregressive model is used as the dynamic model which is subsequently combined with observation model for deriving the optimal window-recursive filter formulae in the sense of minimum mean square error principle. The corresponding test statistics characteristics of system residuals are discussed in details. The test statistics of regional predicted residuals are then constructed in a time-window for model bias testing with two hypotheses, that is, the null and alternative hypotheses. Based on the innovations test statistics, we develop a model bias processing procedure including bias detection, location identification, and state correction. Finally, the minimum detectable bias and bias-to-noise ratio are both computed for evaluating the internal and external reliability of overall system, respectively.


---

## Body

## 1. Introduction

The Kalman filter (KF) has been widely used in various fields, such as target tracking, navigation estimation, distributed network, and multiagent consensus problem [1–5]. Under the assumption of correct system functional and stochastic models, it generates a statistically optimal estimate of underlying system state by recursively operating on streams of noisy input data [6]. However, in real applications, the mathematical system model is not perfect such that the assumption is not always satisfied. For example, false modeled noise, improper system model, or unexpected sudden changes of state existed during filtering process [7], thus degrading the solution precision and even leading to divergence of filter.Many adaptive KF methods have been developed for dealing with abnormal system noise by compensating their influence into stochastic model where the corresponding variance of noise is inflated to reduce the contributions of dynamic model or observation model to the KF final solution [8, 9]. Such methods essentially reduce or even withdraw the inaccurate a priori knowledge of noise by adjusting the weights between the dynamic model and observation model. In fact, these methods are based on stochastic model, that is, variance of noises. A more straightforward strategy is to use their first moment information which can be referred to functional model modification. The unmodelled error is treated as a systematic error, namely, bias which is examined by system innovations based testing theory through generalized likelihood ratio method or the specific chi-square or consistency test [10, 11], and then the bias is further identified and adapted in the ways of data snooping [12, 13]. Another representative method is the interactive multiple model based Kalman filtering. It is in principle a hybrid state estimation filtering where a set of models must be selected in advance to capture the complex target dynamics [14]. Recently, a self-constructed trajectory based dynamic model has been constructed and successfully applied in navigation campaigns [2]. By using multiple states information, it well accommodates the dynamic model bias that cannot be adequately described in the predefined invariant dynamic model of conventional KF, for instance, accelerating, turning, and maneuverings versus constant velocity dynamic model. Substantially, it is an autoregressive (AR) model and this provides a new thought to establish a time-variant dynamic model. However, there is still potentially a considerable possibility of bias occurrence. Moreover, unexpected bias occurs during observation and thus also contaminates the observation model. It is therefore of great importance to develop the quality control strategy for dynamic system.Aiming to construct a general accurate dynamic model, we conduct an optimal window-recursive filter for processing AR-n dynamic model in the sense of minimum mean square error (MMSE) principle. Considering the test statistics characteristics of filtering, the regional innovation test statistics are constructed in a time-window for two hypotheses testing, that is, the null and alternative hypotheses. Based on the innovations test statistics, we further develop a model bias processing procedure including bias detection, location identification, and state correction. The minimum detectable bias and bias-to-noise ratio are both computed for the internal and external reliability analyses of overall system as well. This paper is organized as follows. Section 2 optimally derives the recursive formulae for multiple historical states involved filtering. The characteristics of local model test statistics are discussed in Section 3. The regional test statistics for system innovation are constructed in the two scenarios in Section 4. The system model bias processing procedure is developed in three stages in Section 5. Section 6 investigates the internal and external reliability of our proposed filter model. The concluding remarks are given in Section 7.


## 2. Principle of Multiple Historical States Involved Filtering

Then-order autoregressive model are widely expressed by (1)
x
k
=
a
1
x
k
-
n
+
a
2
x
k
-
n
+
1
⋯
+
a
n
x
k
-
1
+
w
k
,where the subscript k is the time step; n denotes the number of involved epochs; x is the actual state with dimension t; a is the autoregressive coefficient; w is the zero-mean normal distributed process noise with the variance of Q
w. The intuitive representation of (1) in state space form can be given as(2)
x
k
x
k
-
1
⋯
x
k
-
n
+
1
=
a
n
a
n
-
1
⋯
a
1
1
1
×
t
0
1
×
t
0
1
×
t
0
1
×
t
⋯
⋯
⋯
⋯
0
1
×
t
0
1
×
t
1
1
×
t
0
1
×
t
x
k
-
1
x
k
-
2
⋯
x
k
-
n
+
w
k
.However, as n is getting large, the computational burden of such state augmentation will be accordingly increased. For this reason, we rewrite AR-n model (1) into matrix form without state augmentation: (3)
x
k
=
Φ
k
,
k
-
n
:
k
-
1
x
k
-
n
:
k
-
1
+
w
k
,where Φ
k
,
k
-
n
:
k
-
1
=
a
1
a
2
⋯
a
n is the transition matrix transforming the states of the previous n epochs x
k
-
n
:
k
-
1
=
x
k
-
n
T
x
k
-
n
+
1
T
⋯
x
k
-
1
T
T into the current one, x
k
-
n
:
k
-
1 is the vector that consists of the stacked state vectors from the epoch (
k
-
n
) to the epoch (
k
-
1
), and the corresponding covariance matrix is Q
x
^
k
-
n
:
k
-
1. On the other hand, the linear or linearized observation model reads(4)
l
k
=
A
k
x
k
+
ε
k
,where l is the measurement vector with dimension d; ε is the zero-mean normal distributed measurement noise with the variance of R; A is the design matrix. It is easy to compute the predicted state and its variance by(5)
x
-
k
=
Φ
k
,
k
-
n
:
k
-
1
x
^
k
-
n
:
k
-
1
,
(6)
Q
x
-
k
=
Φ
k
,
k
-
n
:
k
-
1
Q
x
^
k
-
n
:
k
-
1
Φ
k
,
k
-
n
:
k
-
1
T
+
Q
w
k
,where x
-
k and Q
x
-
k denote the predicted state and its corresponding variance, respectively; Q
x
^
k
-
n
:
k
-
1 is the variance of x
^
k
-
n
:
k
-
1. Let us start derivations with the following defined linear expression of estimated state:(7)
x
^
k
=
x
-
k
+
K
k
l
k
-
A
k
x
-
k
,where K
k is the so-called gain matrix to be determined. Denoting the a priori and posterior estimate errors x
~
k
,
k
-
1
=
x
-
k
-
x
k and x
~
k
=
x
^
k
-
x
k, respectively, it is rather easy to compute the estimated state variance with(8)
Q
x
^
k
=
E
x
~
k
x
~
k
T
=
E
x
~
k
,
k
-
1
-
K
k
A
k
x
~
k
,
k
-
1
+
ε
k
x
~
k
,
k
-
1
-
K
k
H
k
x
~
k
,
k
-
1
+
ε
k
T
=
Q
x
-
k
-
K
k
A
k
Q
x
-
k
-
Q
x
-
k
A
k
T
K
k
T
+
K
k
A
k
Q
x
-
k
A
k
T
K
k
T
+
K
k
R
k
K
k
T
,where E
(
·
) denotes the expectation operator. Alternative expression of (8) can be written as(9)
Q
x
^
k
=
I
-
K
k
A
k
Q
x
-
k
I
-
K
k
A
k
T
+
K
k
R
k
K
k
Tin which I is the identity matrix which has the same dimension with state vector. Next we will find out K
k in the sense of MMSE satisfying the following condition:(10)
K
k
=
arg
⁡
min
⁡
Tr
⁡
Q
x
^
kwith the fact of(11)
∂
Tr
⁡
Q
x
^
k
∂
K
k
=
-
2
I
-
K
k
H
k
Q
x
-
k
H
k
T
+
2
K
k
R
k
=
0
,where Tr
⁡
(
·
) denotes the trace of matrix. Thereby the unique solution of K
k is solved by (12)
K
k
=
Q
x
-
k
A
k
T
A
k
Q
x
-
k
A
k
T
+
R
k
-
1
.Furthermore, we also obtain the state estimate covariance:(13)
Q
x
^
k
=
I
-
K
k
A
k
Q
x
-
k
.Note that differing from conventional KF formulae derivation, the correlation between x
^
k and x
^
k
-
n
:
k
-
1 exists and should be rigorously considered when the time-window moves forward. Rewrite (8) into following form:(14)
x
^
k
=
I
-
K
k
A
k
x
-
k
+
K
k
l
k
.Inserting (5) into (14) yields(15)
x
^
k
=
I
-
K
k
A
k
Φ
k
,
k
-
n
:
k
-
1
x
^
k
-
n
:
k
-
1
+
K
k
l
k
.Noting that l
k is uncorrelated with x
^
k
-
n
:
k
-
1, the covariance matrix Q
x
^
k
x
^
k
-
n
:
k
-
1 between x
^
k and x
^
k
-
n
:
k
-
1 is then derived in terms of the error propagation law as(16)
Q
x
^
k
x
^
k
-
n
:
k
-
1
=
I
-
K
k
A
k
Φ
k
,
k
-
n
:
k
-
1
Q
x
^
k
-
n
:
k
-
1
.We further symbolize the variance of x
^
k
-
n
:
k in block matrix as(17)
Q
x
^
k
-
n
:
k
=
Q
x
^
k
-
n
:
k
-
1
Q
x
^
k
x
^
k
-
n
:
k
-
1
T
Q
x
^
k
x
^
k
-
n
:
k
-
1
Q
x
^
k
.Thus, with the window moving forward one step, the new epoch will be introduced and the first epoch removed. It is also rather easy to derive the filtering solution x
^
k
-
n
+
1
:
k
+
1 and Q
x
^
k
-
n
+
1
:
k
+
1. It should be pointed out that when the time-window length n
=
1, the multiple historical states involved filtering algorithm reduces to the KF. Its implementation procedure is summarized as follows:(i)
Initialize the historical state information ofx
^
k
-
n
:
k
-
1 and Q
x
^
k
-
n
:
k
-
1 of n epochs.
(ii)
Compute the predicted state vectorsx
-
k and Q
x
-
k by (5) and (6).
(iii)
Sequentially computeK
k, x
^
k, Q
x
^
k, and Q
x
^
k
x
^
k
-
n
:
k
-
1 with (12)~(14) and (16).
(iv)
Savex
^
k
-
n
+
1
:
k along with Q
x
^
k
-
n
+
1
:
k which is the submatrix of (17).
(v)
If no observation interruption occurs, update time epochk to k
+
1 and jump to (ii). Otherwise, restart the filter from (i) to reinitialize the information of multiple states.The above multiple states involved filtering generates optimal estimators of state only when the system noises are mutually uncorrelated and zero-mean normal distributed. However, such assumption is not always satisfied in real applications. Misspecifications in the system model and unexpected errors inevitably generate biased estimation results. It is therefore of great importance to correctly deal with these system model faults and outliers. In the next section, we will investigate the model bias processing theory to ensure the accuracy and reliability of overall system.


## 3. Test Statistics Characteristics of Local Filter Test

According to (3) and (4), three error sources exist in the filter solutions, that is, the system process noise w
k, the observation noise ε
k, and the errors of x
^
k
-
n
:
k
-
1. Based on this point of view, the system models can be reformulated as the following residual equations for pseudo-observations [7]: (18)
v
x
k
=
z
x
k
-
x
^
k
-
w
^
k
,
(19)
v
w
k
=
z
w
k
-
w
^
k
,
(20)
v
l
k
=
z
l
k
-
A
k
x
^
k
,where z
x
k
=
Φ
k
,
k
-
n
:
k
-
1
x
^
k
-
n
:
k
-
1, z
w
k
=
0, and z
l
k
=
l
k with their pseudo-observation variance matrices by(21)
Q
z
x
k
=
Φ
k
,
k
-
n
:
k
-
1
Q
x
^
k
-
n
:
k
-
1
Φ
k
,
k
-
n
:
k
-
1
T
,
Q
z
w
k
=
Q
w
k
,
Q
z
l
k
=
R
k
.Defining a new observation vector z
k
=
z
x
k
T
z
w
k
T
z
l
k
T
T, we obtain its corresponding partitioned matrix Q
z
k
=
Q
z
x
k
Q
z
w
k
Q
z
l
k. By denoting the new residual as v
z
k
=
v
x
k
T
v
w
k
T
v
l
k
T
T, then the following test statistics widely used in system diagnosing of local system model can be constructed as (22)
v
z
k
T
Q
z
k
-
1
v
z
k
~
χ
2
r
,
0
,where r is the redundancy of z
k. In practice, the following alternate chi-square test is more popular: (23)
v
-
k
T
Q
v
-
k
-
1
v
-
k
~
χ
2
r
,
0
,where v
-
k is the system innovation defined as the difference between raw observation and predicted one derived from the predicted state:(24)
v
-
k
=
l
k
-
A
k
x
-
k
.Evidently, it mixes up all the error sources in (18)~(20). Further inserting (5) into (24), the innovation vector arrives at(25)
v
-
k
=
l
k
-
A
k
Φ
k
,
k
-
n
:
k
-
1
x
^
k
-
n
:
k
-
1
.Accordingly, its variance matrix is derived by(26)
Q
v
-
k
=
R
k
+
A
k
Φ
k
,
k
-
n
:
k
-
1
Q
x
^
k
-
n
:
k
-
1
Φ
k
,
k
-
n
:
k
-
1
T
+
Q
w
k
A
k
T
.It should be pointed out that the two test statistics in (22) and (23) have been proved to be equivalent [15]. Therefore, we use the distribution of system innovation vector to test the mathematical model of filter.


## 4. Regional Test Statistics for System Innovation

The aforementioned local test statistics can be potentially used to detect sudden failures or outliers in system model. However, there are two drawbacks: (i) the local testing is too insensitive to detection of the slowly growing sensor errors such as gyroscope drift and accelerometer zero-bias; (ii) only the current observation information is included in the test statistics, thus degrading the reliability of testing results. Therefore, in this section, we introduce a sliding window based test statistics which contain information of multiple system innovations rather than one, thereby enhancing system reliability in the aspect of model fault detection. The predicted residual vectors ofm epochs can be expressed as(27)
v
-
k
-
m
+
1
:
k
=
v
-
k
-
m
+
1
T
v
-
k
-
m
+
1
T
⋯
v
-
k
T
T
.The following null and alternative hypotheses are considered:(28)
H
0
:
v
-
k
-
m
+
1
:
k
~
N
0
,
Q
v
-
k
-
m
+
1
:
k
,
H
α
:
v
-
k
-
m
+
1
:
k
~
N
Δ
v
-
k
-
m
+
1
:
k
,
Q
v
-
k
-
m
+
1
:
k
,where Δ
v
-
k
-
m
+
1
:
k is an unknown residual bias vector with dimension ∑
i
=
k
-
m
+
1
k
d
i. Define the following relationship between Δ
v
-
k
-
m
+
1
:
k and unknown model bias vector b: (29)
Δ
v
-
k
-
m
+
1
:
k
=
C
v
-
k
-
m
+
1
:
k
b
,where C
v
-
k
-
m
+
1
:
k
=
C
v
-
k
-
m
+
1
T
C
v
-
k
-
m
+
2
T
⋯
C
v
-
k
T
T is a projection matrix which can be interpreted as the transformation from system model bias into innovation bias. Here two typical model errors scenarios are considered.Case A.
A bias error vectorb occurs in the dynamic model at epoch s, and the dynamic model at arbitrary epoch is expressed as(30)
x
i
=
Φ
i
,
i
-
n
:
i
-
1
x
i
-
1
+
B
i
b
+
w
i
,where B
i is the coefficient matrix denoting how the bias vector b enters into the dynamics. For convenience, the transition matrix Φ
i
,
i
-
n
:
i
-
1 is denoted as Φ
i in the latter expressions.Then (30) can accommodate the unmodeled error in dynamic model, for example, maneuverings, accelerating for vehicles with constant velocity. The matrix C
v
-
i is conducted by(31)
C
v
-
i
=
-
A
i
M
i
,
i
=
s
,
s
+
1
,
…
,
k
,
M
i
+
1
=
B
i
+
Φ
i
+
1
∑
p
=
1
i
-
s
+
1
δ
i
-
s
-
p
+
2
⊗
I
t
×
t
-
K
s
+
p
-
1
A
s
+
p
-
1
B
s
+
p
-
1
,
M
s
=
B
s
,where ⊗ denotes the Kronecker product; δ
p
=
δ
p
,
1
δ
p
,
2
⋯
δ
p
,
n
T and the elements satisfy the Dirac function:(32)
δ
p
,
q
=
1
,
p
=
q
,
δ
p
,
q
=
0
,
p
≠
q
.Case B.
A bias error vectorb occurs in the observational model at epoch s, and the observational model at arbitrary epoch is expressed as(33)
l
i
=
A
i
x
i
+
H
i
b
+
ε
i
,where H
i is the coefficient matrix of observation bias vector. Thereby the matrix C
v
-
i is computed by(34)
C
v
-
i
=
H
i
-
A
i
N
i
,
i
=
s
,
s
+
1
,
…
,
k
,
N
i
+
1
=
Φ
i
+
1
G
i
+
1
,
N
s
=
0
,
G
i
=
∑
p
=
1
i
-
s
δ
p
⊗
Λ
i
-
p
,
Λ
p
=
K
s
H
s
,
p
=
s
,
Λ
p
=
Φ
p
-
1
-
K
p
-
1
A
p
-
1
Φ
p
-
1
∑
q
=
1
p
-
s
δ
q
⊗
Λ
p
-
q
+
K
p
H
p
,
p
>
s
.


### 4.1. Generalized Likelihood Ratio Test Statistics

Thus far, the alternative hypothesis has been specified. Next we introduce the generalized likelihood ratio test statistic for testingH
0 against H
α as follows [13]:(35)
T
k
-
m
+
1
:
k
=
v
-
k
-
m
+
1
:
k
T
Q
v
-
k
-
m
+
1
:
k
-
1
C
v
-
k
-
m
+
1
:
k
C
v
-
k
-
m
+
1
:
k
T
Q
v
-
k
-
m
+
1
:
k
-
1
C
v
-
k
-
m
+
1
:
k
-
1
C
v
-
k
-
m
+
1
:
k
T
Q
v
-
k
-
m
+
1
:
k
-
1
v
-
k
-
m
+
1
:
k
.Reject H
0 if (36)
T
k
-
m
+
1
:
k
>
χ
2
∑
i
=
k
-
m
+
1
k
d
i
,
0
,where Q
v
-
k
-
m
+
1
:
k is the covariance matrix of v
-
k
-
m
+
1
:
k. As aforementioned, the bias occurs at epoch s, and then (35) can be further rewritten as (37)
T
s
:
k
=
v
-
s
:
k
T
Q
v
-
s
:
k
-
1
C
v
-
s
:
k
C
v
-
s
:
k
T
Q
v
-
s
:
k
-
1
C
v
-
s
:
k
-
1
C
v
-
s
:
k
T
Q
v
-
s
:
k
-
1
v
-
s
:
k
,where the innovation residual vectors from s to k are calculated by(38)
v
-
s
=
l
s
-
A
s
Φ
s
x
^
s
-
n
:
s
-
1
v
-
s
+
1
=
l
s
+
1
-
A
s
+
1
Φ
s
+
1
x
^
s
-
n
+
1
:
s
⋮
v
-
k
=
l
k
-
A
k
Φ
k
x
^
k
-
n
:
k
-
1
.Evidently, the innovation residuals vectors in (38) are strongly correlated with each other. Thus, Q
v
-
s
:
k is not a diagonal matrix and its covariance should be rigorously considered. Let us reformulate the above arbitrary innovation residual vector as(39)
v
-
q
=
∑
i
=
s
q
O
q
,
i
l
i
+
G
q
x
^
s
-
n
:
s
-
1
,
q
=
s
,
…
,
k
.Defining Ψ
p
=
0
n
-
1
t
×
1
I
n
-
1
t
×
n
-
1
t
Φ
p, then G
q and O
q
,
i are accordingly specified as(40)
G
q
=
-
A
q
Φ
q
,
q
=
s
,
G
q
=
-
A
q
Φ
q
∏
p
=
1
q
-
k
-
m
+
1
Ψ
q
-
p
-
K
q
-
p
A
q
-
p
Φ
q
-
p
,
s
<
q
≤
k
,
O
q
,
i
=
I
d
q
×
d
q
,
i
=
q
,
O
q
,
i
=
A
q
Φ
q
K
i
,
i
=
q
-
1
,
O
q
,
i
=
A
q
Φ
q
∏
p
=
1
q
-
i
+
1
Ψ
q
-
p
-
K
q
-
p
A
q
-
p
Φ
q
-
p
K
i
,
i
<
q
-
1
.Inserting (40) into (39), and reformulating (38) into matrix form, we obtain(41)
v
-
s
:
k
=
F
s
:
k
U
s
:
k
,where (42)
F
s
:
k
=
G
s
O
s
,
s
0
0
0
0
G
s
+
1
O
s
+
1
,
s
O
s
+
1
,
s
+
1
0
0
0
G
s
+
2
O
s
+
2
,
s
O
s
+
2
,
s
+
1
O
s
+
2
,
s
+
2
0
0
⋯
⋯
⋯
⋯
⋯
⋯
G
k
O
k
,
s
O
k
,
s
+
1
O
k
,
s
+
2
⋯
O
k
,
k
,
U
s
:
k
=
l
s
l
s
+
1
⋯
l
k
x
^
s
-
n
:
s
-
1
.Therefore, in terms of the error propagation law, the covariance matrix Q
v
-
s
:
k is derived by(43)
Q
v
-
s
:
k
=
F
s
:
k
Q
U
s
:
k
F
s
:
k
T
,where Q
U
s
:
k is the covariance of U
s
:
k. For the independence of subvectors, Q
U
s
:
k is written in a form of block diagonal matrix as(44)
Q
U
s
:
k
=
R
s
0
0
0
0
0
R
s
+
1
0
0
0
0
0
⋯
0
0
0
0
0
R
k
0
0
0
0
0
Q
x
^
s
-
n
:
s
-
1
.


## 4.1. Generalized Likelihood Ratio Test Statistics

Thus far, the alternative hypothesis has been specified. Next we introduce the generalized likelihood ratio test statistic for testingH
0 against H
α as follows [13]:(35)
T
k
-
m
+
1
:
k
=
v
-
k
-
m
+
1
:
k
T
Q
v
-
k
-
m
+
1
:
k
-
1
C
v
-
k
-
m
+
1
:
k
C
v
-
k
-
m
+
1
:
k
T
Q
v
-
k
-
m
+
1
:
k
-
1
C
v
-
k
-
m
+
1
:
k
-
1
C
v
-
k
-
m
+
1
:
k
T
Q
v
-
k
-
m
+
1
:
k
-
1
v
-
k
-
m
+
1
:
k
.Reject H
0 if (36)
T
k
-
m
+
1
:
k
>
χ
2
∑
i
=
k
-
m
+
1
k
d
i
,
0
,where Q
v
-
k
-
m
+
1
:
k is the covariance matrix of v
-
k
-
m
+
1
:
k. As aforementioned, the bias occurs at epoch s, and then (35) can be further rewritten as (37)
T
s
:
k
=
v
-
s
:
k
T
Q
v
-
s
:
k
-
1
C
v
-
s
:
k
C
v
-
s
:
k
T
Q
v
-
s
:
k
-
1
C
v
-
s
:
k
-
1
C
v
-
s
:
k
T
Q
v
-
s
:
k
-
1
v
-
s
:
k
,where the innovation residual vectors from s to k are calculated by(38)
v
-
s
=
l
s
-
A
s
Φ
s
x
^
s
-
n
:
s
-
1
v
-
s
+
1
=
l
s
+
1
-
A
s
+
1
Φ
s
+
1
x
^
s
-
n
+
1
:
s
⋮
v
-
k
=
l
k
-
A
k
Φ
k
x
^
k
-
n
:
k
-
1
.Evidently, the innovation residuals vectors in (38) are strongly correlated with each other. Thus, Q
v
-
s
:
k is not a diagonal matrix and its covariance should be rigorously considered. Let us reformulate the above arbitrary innovation residual vector as(39)
v
-
q
=
∑
i
=
s
q
O
q
,
i
l
i
+
G
q
x
^
s
-
n
:
s
-
1
,
q
=
s
,
…
,
k
.Defining Ψ
p
=
0
n
-
1
t
×
1
I
n
-
1
t
×
n
-
1
t
Φ
p, then G
q and O
q
,
i are accordingly specified as(40)
G
q
=
-
A
q
Φ
q
,
q
=
s
,
G
q
=
-
A
q
Φ
q
∏
p
=
1
q
-
k
-
m
+
1
Ψ
q
-
p
-
K
q
-
p
A
q
-
p
Φ
q
-
p
,
s
<
q
≤
k
,
O
q
,
i
=
I
d
q
×
d
q
,
i
=
q
,
O
q
,
i
=
A
q
Φ
q
K
i
,
i
=
q
-
1
,
O
q
,
i
=
A
q
Φ
q
∏
p
=
1
q
-
i
+
1
Ψ
q
-
p
-
K
q
-
p
A
q
-
p
Φ
q
-
p
K
i
,
i
<
q
-
1
.Inserting (40) into (39), and reformulating (38) into matrix form, we obtain(41)
v
-
s
:
k
=
F
s
:
k
U
s
:
k
,where (42)
F
s
:
k
=
G
s
O
s
,
s
0
0
0
0
G
s
+
1
O
s
+
1
,
s
O
s
+
1
,
s
+
1
0
0
0
G
s
+
2
O
s
+
2
,
s
O
s
+
2
,
s
+
1
O
s
+
2
,
s
+
2
0
0
⋯
⋯
⋯
⋯
⋯
⋯
G
k
O
k
,
s
O
k
,
s
+
1
O
k
,
s
+
2
⋯
O
k
,
k
,
U
s
:
k
=
l
s
l
s
+
1
⋯
l
k
x
^
s
-
n
:
s
-
1
.Therefore, in terms of the error propagation law, the covariance matrix Q
v
-
s
:
k is derived by(43)
Q
v
-
s
:
k
=
F
s
:
k
Q
U
s
:
k
F
s
:
k
T
,where Q
U
s
:
k is the covariance of U
s
:
k. For the independence of subvectors, Q
U
s
:
k is written in a form of block diagonal matrix as(44)
Q
U
s
:
k
=
R
s
0
0
0
0
0
R
s
+
1
0
0
0
0
0
⋯
0
0
0
0
0
R
k
0
0
0
0
0
Q
x
^
s
-
n
:
s
-
1
.


## 5. The System Model Bias Processing Procedure

In this section, based on the testing principle of previous section, we propose a system model bias processing procedure including three steps, that is, bias detection, location identification, and state correction.Step 1 (bias detection).
First we implement the regional testing with test statisticT to examine whether the bias exists in the time-window of filter. Therefore, the overall system model testing with (35) is carried out based on following H
0 and H
α, respectively,(45)
H
0
:
T
k
-
m
+
1
:
k
~
χ
2
∑
i
=
k
-
m
+
1
k
d
i
,
0
,
H
α
:
T
k
-
m
+
1
:
k
~
χ
2
∑
i
=
k
-
m
+
1
k
d
i
,
μ
,where μ is the noncentrality parameter and can be computed by(46)
μ
=
b
T
C
v
-
k
-
m
+
1
:
k
T
Q
v
-
k
-
m
+
1
:
k
-
1
C
v
-
k
-
m
+
1
:
k
b
.
H
0 will be accepted if T
k
-
m
+
1
:
k
<
χ
2
∑
j
=
k
-
m
+
1
k
d
j
,
0; that is, no bias occurs during the filtering process. Otherwise, we will accept the alternative hypothesis H
α.Step 2 (time and location identification).
If the null hypothesisH
0 is rejected, we need to locate the starting time and position of model bias. Since the time when the bias occurs is still unknown, the test statistics screening strategy k is proposed here to find the most likely starting time of bias from the beginning of the time-window k
-
m
+
1 to the current epoch. The starting time and location of bias can be identified as follows:(i)
Similar to the data snooping method, we construct alternative hypotheses for potential blunders through making the elements in coefficient matrixB or H of model bias equal 1 in order: that is, 0
⋯
1
j
⋯
0
T, where j ranges from 1 to dimension of B or H.
(ii)
Compute all the candidates of test statisticsT
i
:
k
i
=
k
-
m
+
1
,
k
-
m
+
2
,
…
,
k with the specified coefficient matrix B or H in (i) in the time-window.
(iii)
Find the maximum value ofT
s
:
k
(
k
-
m
+
1
≤
s
≤
k
) with the corresponding B or H and the epoch index s and the element order j are the most likely starting time and location of bias if the alternative hypothesis H
α is true.Step 3 (state correction).
Once the starting time and location of bias are identified, a bias correction step further needs to be undertaken to make the null hypothesis accepted. It can be fulfilled by additional bias parameterization, which treats the bias vector as a new parameter to be estimated. Accordingly, the dynamics of bias vector is assumed as a constant vector:(47)
b
i
=
b
i
-
1
,
i
=
s
+
1
,
s
+
2
…
,
k
.With this operation, the identified alternative hypothesis turns into a new null hypothesis. Based on the two-step filtering theory of bias treatment [12], the final optimum estimate x
^ which exempts the bias influences can be expressed by two components as(48)
x
^
i
=
x
~
i
+
Δ
x
^
i
,where x
~ is the bias-free estimate calculated as if no bias exists in system model; Δ
x is the estimate correction term induced by bias vector during the recursive filtering and it is given by(49)
Δ
x
^
i
=
V
x
i
b
^
i
.The two terms on the right hand side of (49) can be computed by (50)
b
^
i
=
I
-
K
b
i
A
i
b
^
i
-
1
+
K
b
i
v
-
~
i
,
V
x
i
=
Y
x
i
-
K
~
x
i
S
i
,where the bias-free predicted residual v
-
~
i is as follows:(51)
v
-
~
i
=
l
i
-
A
i
Φ
i
x
~
i
-
n
:
i
-
1
.The introduced intermedia matrix is(52)
Y
x
i
=
Φ
i
-
1
S
x
i
-
1
+
B
i
-
1and the gain matrices of bias-free state x and bias state b are(53)
K
~
x
i
=
Q
x
-
~
k
A
k
T
A
k
Q
x
-
~
k
A
k
T
+
R
k
-
1
,
(54)
K
b
i
=
J
i
+
1
V
x
i
T
A
i
T
+
H
i
R
i
-
1in which Q
x
-
~
k
=
Φ
k
Q
x
~
k
-
n
:
k
-
1
Φ
k
T
+
Q
w
k and J
i
+
1 is recursively computed by means of(55)
J
i
+
1
=
J
i
-
J
i
S
i
T
A
i
Q
x
-
~
i
A
i
T
+
R
i
+
S
i
J
i
S
i
T
-
1
S
i
J
iwith (56)
S
i
=
A
i
Y
x
i
+
H
i
.Note that Q
x
-
~
k in (53) is the variance of bias-free estimated state. The initial information of Y introduced at epoch s is given as Y
x
s
=
0. Up to now, the bias-corrected state estimation under the alternative hypothesis H
α can be recursively realized.


## 6. Reliability Analysis

When we implement the overall system innovation model testing, to ensure the probability of success in bias detection, we should control the probabilities of committing both type I errorα and type II error β which correspond to the probabilities of making wrong decisions during test [16, 17]. According to Neyman-Pearson principle, if the size of type I error α is fixed, type II error β is naturally associated with the alternative hypothesis H
α. In this section, with the given sizes α and β, we will have insight into the system reliability, which is also known as internal reliability and external reliability.


### 6.1. Minimal Detectable Bias

The minimal detectable bias (MDB) is considered as a measure of the internal reliability, which indicates a lower bound of the bias that can be successfully detected. It has a strong relation with the alternative hypothesis and size of type II errorβ, the probability of accepting null hypothesis when H
α is true. It is worthy of pointing out that one should maintain a protection against the size of type II error β. It depends on the parameters involved in H
α, that is, the chosen level of significance α, the noncentrality parameter μ, and the degrees of freedom d. More straightforwardly, type II error manifests the power of statistical test γ through the following function:(57)
γ
=
1
-
β
f
α
,
μ
,
d
,where the subscript f is the function symbol. Another expression of (57) reads(58)
μ
0
=
μ
f
1
α
0
,
β
0
,
d
0
=
μ
f
2
α
0
,
γ
0
,
d
0
.Obviously, given the protected power of statistical test γ
0, α
0, and d
0
=
∑
j
=
s
k
d
j against type II error, we obtain the corresponding noncentrality parameter μ
0. Consider the fixed noncentrality parameter μ
0 in (58) and relation with system innovation information of alternative hypothesis in (46); accordingly the MDB is then calculated by(59)
b
M
D
B
=
μ
0
C
v
-
k
-
m
+
1
:
k
T
Q
v
-
k
-
m
+
1
:
k
-
1
C
v
-
k
-
m
+
1
:
k
.


### 6.2. Bias-to-Noise Ratio

Besides internal reliability mentioned above, the influences of the model bias on the final state solutions also should be evaluated since the model biases are not always detectable. As an important external reliability criterion, the bias-to-noise ratio is an efficient tool for quantizing the influences of model bias of overall system. A scalar squared bias-to-noise ratio is defined as [13](60)
c
B
N
R
=
Δ
x
^
T
Q
x
~
-
1
Δ
x
^
,where Δ
x
^ is the influences term of model bias on state in (48) for the alternative hypothesis H
α and Q
x
~ is the covariance matrix of bias-free state for the null hypothesis. Further inserting (13) and (49) into (60) yields(61)
c
B
N
R
=
b
^
T
V
x
T
I
-
K
~
k
A
k
Q
x
-
~
k
-
1
V
x
b
^
T
.Evidently, a large value of c
B
N
R indicates the significant influences of model bias b on state estimation, while a small value of c
B
N
R exhibits a good external reliability of system. It is therefore of importance to design a reasonable system model structure for minimizing the influences of model bias on solutions.


## 6.1. Minimal Detectable Bias

The minimal detectable bias (MDB) is considered as a measure of the internal reliability, which indicates a lower bound of the bias that can be successfully detected. It has a strong relation with the alternative hypothesis and size of type II errorβ, the probability of accepting null hypothesis when H
α is true. It is worthy of pointing out that one should maintain a protection against the size of type II error β. It depends on the parameters involved in H
α, that is, the chosen level of significance α, the noncentrality parameter μ, and the degrees of freedom d. More straightforwardly, type II error manifests the power of statistical test γ through the following function:(57)
γ
=
1
-
β
f
α
,
μ
,
d
,where the subscript f is the function symbol. Another expression of (57) reads(58)
μ
0
=
μ
f
1
α
0
,
β
0
,
d
0
=
μ
f
2
α
0
,
γ
0
,
d
0
.Obviously, given the protected power of statistical test γ
0, α
0, and d
0
=
∑
j
=
s
k
d
j against type II error, we obtain the corresponding noncentrality parameter μ
0. Consider the fixed noncentrality parameter μ
0 in (58) and relation with system innovation information of alternative hypothesis in (46); accordingly the MDB is then calculated by(59)
b
M
D
B
=
μ
0
C
v
-
k
-
m
+
1
:
k
T
Q
v
-
k
-
m
+
1
:
k
-
1
C
v
-
k
-
m
+
1
:
k
.


## 6.2. Bias-to-Noise Ratio

Besides internal reliability mentioned above, the influences of the model bias on the final state solutions also should be evaluated since the model biases are not always detectable. As an important external reliability criterion, the bias-to-noise ratio is an efficient tool for quantizing the influences of model bias of overall system. A scalar squared bias-to-noise ratio is defined as [13](60)
c
B
N
R
=
Δ
x
^
T
Q
x
~
-
1
Δ
x
^
,where Δ
x
^ is the influences term of model bias on state in (48) for the alternative hypothesis H
α and Q
x
~ is the covariance matrix of bias-free state for the null hypothesis. Further inserting (13) and (49) into (60) yields(61)
c
B
N
R
=
b
^
T
V
x
T
I
-
K
~
k
A
k
Q
x
-
~
k
-
1
V
x
b
^
T
.Evidently, a large value of c
B
N
R indicates the significant influences of model bias b on state estimation, while a small value of c
B
N
R exhibits a good external reliability of system. It is therefore of importance to design a reasonable system model structure for minimizing the influences of model bias on solutions.


## 7. Concluding Remarks

This paper proposes a model bias processing approach for the autoregressive (AR) dynamic model involved dynamic system. The optimal time-window recursive filtering formulae are derived in the case of AR-n dynamic model in the sense of minimum mean square error (MMSE) principle. The test statistics of regional system innovation, that is, predicted residuals under the null and alternative hypotheses, are constructed for model bias testing specified in two bias occurrence scenarios. Based on the system innovation test statistics, we develop a model bias processing procedure including bias detection, location identification, and state correction. The minimum detectable bias and bias-to-noise ratio are both computed and comprehensively analyzed for evaluating the internal reliability and external reliability, respectively, of overall system. In the future, since the recursive time-window will inevitably generate the strong correlations between innovations and state estimates, the multiple hypotheses will be further constructed considering the probability of wrong decisions on the bias location when the null hypothesis is rejected. In addition, the computations for the regional testing and processing of model bias are time-consuming; thus, an efficient computation strategy is also desired in the future work.



---

*Source: 4860785-2016-03-15.xml*