Adaptive Filter as Efficient Tool for Data Assimilation under Uncertainties

Abstract

In this contribution, the problem of data assimilation as state estimation for dynamical systems under uncertainties is addressed. This emphasize is put on high-dimensional systems context. Major difficulties in the design of data assimilation algorithms is a concern for computational resources (computational power and memory) and uncertainties (system parameters, statistics of model, and observational errors). The idea of the adaptive filter will be given in detail to see how it is possible to overcome uncertainties as well as to explain the main principle and tools for implementation of the adaptive filter for complex dynamical systems. Simple numerical examples are given to illustrate the principal differences of the AF with the Kalman filter and other methods. The simulation results are presented to compare the performance of the adaptive filter with the Kalman filter.

1. Introduction

In this chapter, the adaptive filter (AF) is considered as a computational device that yields estimates of the system state by minimizing recursively (in time) the error between the predicted output of the device and its observed signal in real time. As the main objective of the AF is to produce estimates of the state in high-dimensional systems (HdSs), we shall focus the attention on the mathematical form of the AF in a state-space form as that used in the Kalman filter (KF) [1]. In this chapter, the HdS is referred to as a system whose state dimension is of order $\mathrm{O}\left({10}^7\right)-\mathrm{O}\left({10}^8\right)$.

The assimilation problem in this chapter is formulated as a standard filtering problem. For simplicity, let the dynamical system be described by the equation

$$x\left(t+1\right)=\Phi x\left(t\right)+w\left(t\right),x\left(0\right)=x_0,t=0,1,\dots$$E1

where $x\left(t\right)$ is the system state at the $t$ time instant. At each time instant $t$, we are given the observation for the system output

$$z\left(t+1\right)=Hx\left(t+1\right)+v\left(t+1\right),t=0,1,\dots$$E2

In (1) and (2), $w\left(t\right)$ is the model error (ME), $v\left(t\right)$is the observation error (ObE), and $\Phi$ represents the system dynamics. In general, the system (1) and (2) may be nonlinear with $\Phi x=f\left(x\right)$, $\mathit{Hx}=h\left(x\right)$. The filtering problem for a partial observed dynamical system (1) and (2) is to obtain the best possible estimate for the state $x\left(t\right)$at each instant $t$, given the set of observations $Z\left(1:t\right)=\left[z\left(1\right),\dots ,z\left(t\right)\right]$.

There exist different techniques to solve estimation problems. The simplest approach is related to linear estimator [2], since it requires only first two moments. Linear estimation is frequently used in practice when there is a limitation in computational complexity. Among others, the widely used methods are maximum likelihood, least squares, method of moments, the Bayesian estimation, minimum mean square error (MMSE), etc. For more details, see [3].

There are limitations of optimal filters. In practice, the difficulties are numerous: the statistics of signals which may not be available or cannot be accurately estimated; there may not be available time for statistical estimation (real-time); the signals and systems may be non-stationary; memory required and computational load may be prohibitive. All these difficulties become insurmountable, especially for HdSs.

In order to deal with real-time applications, the AFs appear to be a valuable tool in solving estimation problems when there is no time for statistical estimation and when we are dealing with non-stationary signals and/or systems environment. They can operate satisfactorily in unknown and possibly time-varying environments without user intervention. They improve their performance during operation by learning statistics from current signal observations. Finally, they can track variations in the signal operating environment [4].

It is well-known that the MMSE estimator in the class of Borel measurable (with respect to (wrt) $Z\left(1:t\right)$) functions is given by the conditional mean

$$\widehat{x}\left(t\right)=E\left[x\left(t\right)/Z\left(1:t\right)\right]$$E3

Under standard conditions, related to the noise sequences $w\left(t\right),v\left(t\right)$ (Gaussian i.i.d.—identically independent (temporal) distributed), the estimate (3)$\widehat{x}\left(t\right)$ for $x\left(t\right)$ can be obtained from the KF in the recursive form

$$\widehat{x}\left(t+1\right)=\widehat{x}\left(t+1/t\right)+\mathit{K\zeta }\left(t+1\right),$$E4

$$\widehat{x}\left(t+1/t\right)=\Phi \widehat{x}\left(t\right),\zeta \left(t+1\right)=z\left(t+1\right)-\widehat{x}\left(t+1/t\right)$$E5

$$K\left(t\right)=M\left(t\right)H^T{\left[\mathit{HM}\left(t\right)H^T+R\right]}^{-1}$$E6

$$M\left(t+1\right)=\Phi P\left(t\right){\Phi }^T+Q$$E7

$$P\left(t\right)=\left[I-K\left(t\right)H\right]M\left(t\right){\left[I-K\left(t\right)H\right]}^{\mathrm{T}}+K\left(t\right)R{K}^T\left(t\right)$$E8

In (4)–(8), $Q,R$ are the covariance matrix for $w\left(t\right)$ and $v\left(t\right)$, respectively. One sees that $K=K\left(M\right)$—the gain matrix, is a function of $M≔M\left(t+1\right)$—the error covariance matrix (ECM) for the state prediction error (PE) $e\left(t+1/t\right)$ which is defined as $e\left(t+1/t\right)≔x\left(t+1\right)-\widehat{x}\left(t+1\right)$, and $\zeta \left(t+1\right)$ is known as innovation vector. Note from (7) and (8) that the ECM $M\left(t+1\right)$ can be found by solving the matrix nonlinear Algebraic Riccati equation (ARE). Generally speaking, a solution of the ARE is not unique. Conditions must be introduced for ensuring an existence of a unique non-negative definite solution [5]. It is remarkable that the ECMs $P,M$ in (7) and (8) do not depend on observations; therefore, they can be computed in advance, offline, given the system matrices and the noise covariances. The same remark is valid for the gain matrix $K$ in (6). In contrast, the gain in the AF is observation-dependent [6] (see Section 2).

Under the most favorable conditions (perfect knowledge of all system parameters and noise statistics), for a dynamical system with dimension of order ${10}^7$–${10}^8$, it is impossible to solve the ARE (due to computational burden), not to say about storing $M,P$. To overcome these difficulties, the AF is proposed. Mention that the KF is also an MMSE filter in the complete Hilbert space of random variables, $E{\left|\xi \right|}^2<\infty$, with scalar product ${\left({\xi }_1,{\xi }_2\right)}_H=E<{\xi }_1{\xi }_2>$ and the norm $\left|\left|\xi \right|\right|={\left[E\left({\left|\xi \right|}^2\right)\right]}^{\left\{1/2\right\}}$.

For the nonlinear models, there are KF variants, among those are the extended KF (EKF) [7], the unscented KF (UKF, [8]), and the Ensemble Kalman filter (EnKF, [9]). In the EnKF, the ECM is a sampled ECM whose samples are generated using samples of the state variable, and consequently the ECM in the KF becomes a sampled ECM. For an example of application of the EnKF for data assimilation in geophysical data assimilation with high dimensional model, see [10]. Another class of ensemble filtering technique is a class of particle filters (PF, [11]). The basic idea of the PF (also the EnKF) is to use a discrete set of weighted $n$ particles to represent the distribution of $x\left(t\right)$, where the distribution is updated at each time by changing the particle weights according to their likelihoods.

Despite a possible implementation of the KF variants, they might still be seriously biased because the accuracy of the KF update requires linearity of the observation function and Gaussianity of the distribution of system state $x\left(t\right)$. In reality, the KF (4)–(8) may be biased and unstable, even divergent [12]. Today, the PF algorithms are ineffective for HdS data assimilation.

In this chapter, we shall show how the AF can be efficient in dealing with uncertainties existing in the filtering problem (1) and (2). In Section 2, a brief outline of the AF is given. The main features of the AF, which are different to those of the KF, are presented. This concerns the optimal criteria, stabilizing gain structure, optimization algorithms. Section 3 shows in detail how the AF is capable of dealing with uncertainties in the specification of ME covariance. The hypothesis on a subspace of ME is presented in Section 4 from which one sees how one can make order reduction for representing the bias and ME covariance. Simple numerical examples on one- and two-dimensional systems are given in Section 5 to illustrate in details the differences between the AF and the traditional KF. Numerical experiments on low and high dimensional systems are given in Section 6 to demonstrate how the AF algorithm works. The performance comparison between the AF and KF, for both situations of perfect knowledge of ME statistics and that with ME uncertainties, is also presented. Conclusions and perspectives of the AF are summarized in Section 7.

2. Adaptive filter

The AF is originated from [13]. It is constructed for estimating the state of a dynamical system based on partially observed quantities related in some way to the system state. As reported before, for linear systems contaminated by Gaussian noise, the MMSE estimate can be obtained by the KF. Since publication of [1] in 1960, an uncountable number of works are done for solving engineering problems by KF, in all engineering fields, as well as many modifications have been proposed. The reasons for the need in modification of the KF are numerous, but mostly related to nonlinear dynamics, parameter uncertainties in specification of system parameters, bias of ME, unknown statistics of ME, model reduction. With the rapid progress of computer technology (computer power, memory, …), various simplified versions of KF are suggested for solving filtering (or data assimilation) for HdSs, in particular, in meteorology and oceanography.

Direct application of the KF to HdSs is impossible due to the limit in computer power, memory, and computational time. In particular, the KF requires to solve the matrix AREs (7) and (8) for computing ECMs $M\left(t\right),P\left(t\right)$. Storing such matrices is impossible, not to say on computational time.

Different simplified approaches are proposed for overcoming difficulties in the application of the KF. The example of successful tool for solving data assimilation problems in HdSs is the EnKF [9]. In the EnKF, an ensemble of error samples, of small size, is generated on the basis of model states to approximate the ECMs. In practice of data assimilation for HdSs, it is possible to generate only ensembles of moderate sizes (of order $O\left(100\right)$) by model integrations over the assimilation window (time interval between two successive arrivals of observations) since one such integration takes several hours! The other approach like PF is based on sampling from conditional distributions. Theoretically, this approach is more appropriate for nonlinear problems because no linearization is required as in the EKF (Extended KF based on linearization technique). However, even for filtering problems with state dimensions of order $O\left(10\right)$, relatively large ensembles of size $O\left(10000\right)$ would be required in order to produce reasonably good performance.

The AF in [13] is based on the different idea. Here, no linearization is required for nonlinear filtering problems. For the problem (1) and (2), the filter is of the form (4) and (5) but the gain $K=K\left(\theta \right)$is assumed to be of a given stabilizing structure [6]. It means that $K$ is parametrized by some vector of unknown parameters $\theta \in \mathrm{\Theta }$ so that the filter (4) and (5) with the gain $K\left(\theta \right),\forall \theta \in \mathrm{\Theta }$is stable. It is well-known that under mild conditions, the solution of the ARE will tend (quickly) to stationary solution $M_{\infty }$ and so the gain (6), to the stationary gain $K_{\infty }$. Moreover, the innovation $\zeta \left(t+1\right)=z\left(t+1\right)-\widehat{z}\left(t+1/t\right)$, representing the error for the output prediction$\widehat{z}\left(t+1/t\right)≔H\widehat{x}\left(t+1/t\right)=H\Phi \widehat{x}\left(t\right)$, is unbiased and of minimum variance. This fact leads to the idea to seek the optimal vector $\theta$by minimizing

$$J\left(\theta \right)≔E\left[\Psi \left(\zeta \right)\right]\to argmin\theta$$E9

here $E\left[.\right]$ denotes the average in a probabilistic sense. For stationary systems (1) and (2), if we assume the validity of the ergodic hypothesis, the average value in a probabilistic sense, expressed in (7), is almost everywhere equivalent to the time average (for large time of running the dynamical system). The optimal ${\theta }^{\ast }$ can be found by solving the equation

$${\nabla }_{\theta }J\left(\theta \right)={\nabla }_{\theta }E\left[\Psi \left(\zeta \right)\right]=E\left[{\nabla }_{\theta }\Psi \left(\zeta \right)\right]=0$$E10

A stochastic approximation (SA) algorithm for solving (10) can be written out

$$\theta \left(t+1\right)=\theta \left(t\right)-\gamma \left(t\right){\nabla }_{\theta \left(t\right)}\Psi \left[\mathrm{\zeta }\left(\mathrm{t}+1\right)\right]$$E11

Conditions related to the sequence of positive scalar $\gamma \left(t\right)$ for ensuring a convergence $\left\{\theta \left(t\right)\right\}$in the procedure (11) are

$$\gamma \left(t\right)>0,\sum _{t=0}^{\infty }\gamma \left(t\right)=\infty ,\sum _{t=0}^{\infty }{\gamma }^2\left(t\right)<\infty$$E12

One of the most advantages of the SA algorithm (11) is that, instead of computing the gradient of the cost function (9) (which requires knowledge of probability distribution), the algorithm (11) is based on the knowledge of only the gradient of sample cost function $\Psi$ (wrt to $\theta$) which can be easily evaluated numerically.

Comment 2.1. Generally speaking, the convergence rate of the algorithm (11) and (12) is $O\left(1/t\right)$. It is possible to improve the convergence rate of the SA by averaging of the iterates,

$${\theta }_m\left(t\right)=\frac{1}{t}\sum _{t\prime =1}^t\theta \left(t\prime \right)$$E13

For more details, see [14].

Comment 2.2. For high HdSs, even with $\theta$ being of moderate dimension, instead of the algorithm (12) or (13), the SPSA (Simultaneous Perturbation Stochastic Perturbation) algorithms in [15, 16] are of preference. That is due to the fact that integration of HdS over the assimilation window is very expensive. These algorithms generate stochastic perturbation $\mathit{\delta \theta }={\left(\delta {\theta }_1,\dots ,\delta {\theta }_n\right)}^T$ with components as Bernoulli i.i.d. realizations. Each $i^{\mathit{th}}$component of the gradient-like (pseudo-gradient) vector is computed as the divided difference $\delta \Psi /\delta {\theta }_i$ where $\Psi ≔\Psi \left[\mathrm{\theta }+\mathrm{\delta \theta }\right]-$$\Psi \left[\mathrm{\theta }\right]$. This allows to evaluate the gradient-like vector by only two or three time integration of the numerical model.

For details on the SPSA algorithm and its convergence rate, see [15, 16].

3. Covariance uncertainties and AF

3.1 Covariance uncertainties

3.1.1 Adjoint approach

As seen from (11) and (12), implementation of SA algorithms is much simpler for searching optimal gain parameters compared to the other optimization methods. The SA algorithms require only numerical computing derivatives ${\nabla }_{\theta }\Psi \left[\zeta \left(t+1\right)\right]$ of the sample cost function $\Psi \left[\zeta \left(t+1\right)\right]$ wrt $\theta$ evaluated at $\theta \left(t\right)$ and $\gamma \left(t\right)$is a scalar which can be chosen a priori, for example, as $\gamma \left(t\right)=\frac{1}{t}$. That is possible due to introducing the ergodic hypothesis on of the system (1) and (2) from which there exists an asymptotic optimal gain

First, consider the situation when the vector of parameters consists of are all elements of $K,\theta =K.$ Compute the innovation vector,

$$\mathrm{\zeta }\left(\mathrm{t}+1\right)=z\left(t+1\right)-\widehat{z}\left(t+1/t\right)=z\left(t+1\right)-H\Phi \widehat{x}\left(t\right)$$

$$\widehat{x}\left(t\right)=\widehat{x}\left(t-1/t\right)+K\left(\theta \right)\mathrm{\zeta }\left(\mathrm{t}\right)$$

$${\mathrm{\delta }}_{\mathrm{\theta }}\mathrm{\zeta }\left(\mathrm{t}+1\right)=-\mathit{H\Phi }{\mathrm{\delta }}_{\mathrm{\theta }}\widehat{\mathrm{x}}\left(\mathrm{t}\right)=-\mathit{H\Phi }{\mathrm{\delta }}_{\mathrm{\theta }}\mathrm{K}\left(\mathrm{\theta }\right)\mathrm{\zeta }\left(\mathrm{t}\right)=-\mathit{\text{HΦ}}\text{δK}\mathrm{\zeta }\left(\mathrm{t}\right).$$

$${\mathrm{\delta }}_{\mathrm{K}}\Psi \left[\mathrm{\zeta }\left(\mathrm{t}+1\right)\right]={\mathrm{\delta }}_{\mathrm{K}}<\mathrm{\zeta }\left(\mathrm{t}+1\right),\mathrm{\zeta }\left(\mathrm{t}+1\right)>=2<\mathrm{\zeta }\left(\mathrm{t}+1\right),{\mathrm{\delta }}_{\mathrm{K}}\mathrm{\zeta }\left(\mathrm{t}+1\right)>$$

$$=-2<\mathrm{\zeta }\left(\mathrm{t}+1\right),\mathit{\text{HΦ}}\text{δKζ}\left(\mathrm{t}\right)>=-2<{\Phi }^{\mathrm{T}}{\mathit{H}}^{\mathrm{T}}\mathrm{\zeta }\left(\mathrm{t}+1\right),\mathrm{\delta K\zeta }\left(\mathrm{t}\right)>.$$

Let us compute derivatives of the sample cost function $\Psi$ wrt the elements $K_{\mathit{ij}}$ of the gain $K$. To do so, one needs to integrate the adjoint operator ${\Phi }^T$ s.t. the forcing $f≔-H^T\zeta \left(t+1\right)$ which yields $\psi ≔-{\Phi }^T{H}^T\zeta \left(t+1\right)$ and hence

$$\mathit{\delta \Psi }/\delta K_{\mathit{ij}}=-{\psi }_i\zeta \left(t,j\right),$$E14

here ${\psi }_i$ is the$i^{\mathit{th}}$component of $\psi ,\zeta \left(t,j\right)$ the $j^{\mathit{th}}$component of $\zeta \left(t\right)$. The AF now takes the form

$$\widehat{x}\left(t+1\right)=\widehat{x}\left(t+1/t\right)+\mathit{K\zeta }\left(t+1\right),$$E15

$$\widehat{x}\left(t+1/t\right)=\Phi \widehat{x}\left(t\right),\zeta \left(t+1\right)=z\left(t+1\right)-H\widehat{x}\left(t+1/t\right),$$E16

$$K\left(t+1\right)=K\left(t\right)-\gamma \left(t\right){\nabla }_K\Psi \left[\mathrm{\zeta }\left(\mathrm{t}+1\right)\right],$$E17

where ${\nabla }_K\Psi \left(\mathrm{\zeta }\left(\mathrm{t}+1\right)\right)$ is the gradient vector whose components are computed by (14). In the AF (15)–(17), no matrix ARE (see (7) and (8) in the KF) is involved. The AF (15)–(17) is quite realizable for HdSs, since at each assimilation instant we need to integrate only the direct model to produce the forecast (16) and (eventually) an adjoint model over the assimilation window for computing ${\nabla }_K\Psi \left[\mathrm{\zeta }\left(\mathrm{t}+1\right)\right]$ [13].

3.1.2 Simultaneous perturbation stochastic approximation (SPSA) approach

Remark that in the form (14) the adjoint operator ${\Phi }^T$ would be available to implement the AF. It is well-known that construction of numerical code for ${\Phi }^T$ is a very difficult and heavy task, especially for meteorological and oceanic numerical models which are HdSs and nonlinear (linearization is required).

A comparison study of the AF with other assimilation methods is done in [17]. Compared to the AF, the widely used variational method (VM) minimizes the distance between the observations available (for example, the observations of the whole set $Z\left[1:T\right]$) and the outputs of the dynamical system. This optimization problem is carried out in the phase space, hence is very difficult and expensive. Theoretically, a simplification is possible subject to (s.t.) the condition of linearity of the dynamical system: in this case, one can reformulate the VM minimization problem as searching the best estimate for the initial system state $x\left(0\right)$. For HdSs, to ensure a merely high quality estimate for $x\left(0\right)$, it is necessary: (i) to take the observation window as large as possible; (ii) to parameterize the initial state by some parameters (using a slow manifold, for example). Iterative minimization procedures require usually $O\left(10\right)$ iterates involving integrating the direct and adjoint models over the window $\left[1,T\right]$. For an unstable dynamics, integration of direct and adjoint equations over a long period naturally amplify the initial errors during assimilation process. For a more detailed comparison between the AF and VM, see [17].

Thus, if the ergodic conditions hold, there exists an optimal stationary gain and the AF in limit will approach to the optimal one in the given class of stable filters. It is important to emphasize that up to this point, no covariance matrices $Q,R$ are specified. It means that the AF in the form (12) is robust to uncertainties in the specification of the covariances of the ME and ObE.

3.2 Stability of the AF

One of the main features of the AF is related to its stability.

For simplicity of presentation, in the previous section, the AF algorithm is written out under the assumption (13). In practice, application of the AF in the form (13) is not recommended since instability may occur. It is easy to see that the transition matrix of the filter is given by

$$L=\left(I-\mathit{KH}\right)\Phi$$E18

It is evident that if we do not take care on the structure of $K$, varying stochastically all elements of $K$ can lead to instability of $L$ and the filter will be exploded. Moreover, for HdSs, the number of elements of $K$ is very large. It is therefore primordial to choose a parametrized stabilizing structure for $K$ (depending on $\theta$) to ensure a stability of $L$ and reducing a number of tuning parameters. This question is addressed in [6]. One of possible structures for $K$ is of the form

$$K\left(t\right)=P_r{\mathrm{\Theta }K}_e,K_e≔M_e{H}_e^T{\left[H_e{M}_e{H}_e^T+\mathrm{\alpha I}\right]}^{-1},H_e≔H{P}_r$$E19

where $P_r\in R^{\mathit{nxr}}$ is a matrix with dimensions $n\times r$, $r$ is the dimension of the reduced space (equal to the number of unstable eigenvectors (EiVecs)of $\Phi$), the matrix $M_e$ is a strictly positive symmetric definitive playing the role of the ECM in the reduced space $R_e$, $\mathrm{\Theta }$ is a diagonal matrix with diagonal elements ${\theta }_i$ whose values belong to $\left(ϵ,2-ϵ\right)$, i.e. ${\theta }_i$$\in \left(ϵ,2-ϵ\right)$, with $ϵ\in \left(0,1\right)$ whose value depends on the modulus of the first stable eigenvalue (EiV) [6]. We will refer to the filter s.t. (19) with $\Theta =\mathit{Id}$ ($\mathit{Id}$ is the identity matrix of appropriate dimension) as a nonadaptive filter (NAF). In the AF, the parameters ${\theta }_i$ are adjusted each time when a new observation arrives, to minimize the cost function (9). Thus ${\theta }_i$ is a time-varying function. As to the matrix $P_r$, its choice is important to ensure a filter stability. One simple and efficient procedure (called Prediction Error Sampling Procedure—PeSP) to generate $P_r$ is to use the power orthogonal iteration method [18] which allows to compute real leading Schur vectors (SchVecs) of ${\Phi }_i$. The advantage of using the SchVecs compared to the EiVecs, is that they are real and their computation is stable. It is seen that the optimal AF is found in a class of stable filters which is stable even for an unstable numerical model. As to the VM, the optimal trajectory is found on the basis of only the numerical model with the initial state as a control vector. It means that for unstable dynamics, the errors in the forcing or numerical errors arising during computations will be amplified and lead to large estimation error growth. More seriously, the VM requires a large set of observations and large number of iterations (i.e., many forward and backward integrations of the direct and adjoint models) which naturally leads to increase of estimation error too.

3.3 On improving the initial gain

Consider the gain structure (19). Suppose that $M_e$ has been chosen in agreement with the required stability conditions. Before tuning the parameters ${\theta }_i$ to minimize the cost function (9), remark that stability of the filter is still ensured for the following gain:

$$K\left(t\right)=P_r{\mathrm{\Lambda \Theta }K}_e,$$E20

where $\mathrm{\Lambda }=diag\left({\lambda }_1,\dots ,{\lambda }_r\right),{\lambda }_k\in \left(0,1\right)$ since then for $\Gamma =\mathrm{\Lambda \Theta }=diag\left({\gamma }_1,\dots ,{\gamma }_r\right)$ it implies $\gamma \left(t\right)=\lambda \left(t\right)\theta \left(t\right)$$\in \left(ϵ,2-ϵ\right)$, where $ϵ\in \left(0,1\right)$. Writing the equation for the filtered error (FE) $e_f\left(t\right)≔\widehat{x}\left(t\right)-x\left(t\right)$one sees that the matrix $L$ in (18) also represents the transition of the FE $e_f\left(t\right)$. It means that it is possible to choose a more optimal initial gain by solving, for example, the minimization problem

$$J_o\left(\mathrm{\Lambda }\right)={‖L\left(\mathrm{\Lambda }\right)‖}_2^2\to {\mathit{min}}_{\mathrm{\Lambda }}$$E21

The problem (21) is solved without using the observations, hence it is offline. Once the optimal ${\mathrm{\Lambda }}^{\ast }$ has been found, the standard AF is implemented s.t. the filter gain

$$K\left(t\right)=P_r\mathrm{\Theta }{K}_e^{\ast },K_e^{\ast }≔{\mathrm{\Lambda }}^{\ast }{K}_e$$E22

It is seen that using the structure (22) this optimization procedure does not require the information on the ME statistics.

4. Joint estimation of state and model error in AF

The previous section shows how the AF is designed to deal with the difficulty in specification of covariances of the ME and ObE. This is done without exploiting a possibility to determine, more or less correctly, a subspace for the ME. If such a subspace can be determined without major difficulties, it would be beneficial for better estimating the AF gain and improving the filter performance. In [19], the hypothesis of the structure of the ME has been introduced and a number of experiments have been successfully conducted.

There is a long history of joint estimation of state and ME for filtering algorithms, in particular, with the bias and covariance estimation. One of the most original approaches, dealing with the treatment of bias in recursive filtering (known as bias-separated estimation—BSE), is carried out by Friedland in [20]. He has shown that the MMSE state estimator for a linear dynamical system augmented with bias states can be decomposed into three parts: (1) bias-free state estimator; (2) bias estimator; and (3) blender. This BSE approach has the advantage that it requires fewer numerical operations than the traditional augmented-state implementation and avoids numerical ill-conditioning compared to the case of bias-separated estimation by filtering technique.

It is common to treat the bias as part of the system state and then estimate the bias as well as the system state. There are two types of ME—deterministic (DME) and stochastic (SME). Generally speaking, a suitable equation can be introduced for the ME. In the presence of bias, under the assumption on constant $\mathit{b}$, instead of (1) one has

$$x\left(t+1\right)=\Phi x\left(t\right)+b\left(t\right)+w\left(t\right),b\left(t+1\right)=b\left(t\right),t=0,1,2\dots$$E23

To introduce a subspace for the variables $w\left(t\right),b\left(t\right)$ the SME and DME in (23), let

$$w=G_w\rho ,b=G_{b}d$$

$$G_w\in R^{n\times n_w},G_b\in R^{n\times n_d},n\ge n_w,n\ge n_b$$E24

Generally speaking, $G_w,G_b$ are unknown, and finding reasonable hypothesizes for them is desirable but not self-evident. In [19], one hypothesis for $G_w,G_b$ has been introduced (it will be referred to as Hypothesis on model error—HME).

The information on $G_w,G_b$, given in (25), allows to better estimate the DME $b$ and SME $w$ for improving the filter performance, especially for $n_b$ < n, $n_w$ < n in a HdS setting. The difficulty, encountered in practice of operational forecasting systems, is that (practically) nothing is given a priori on the space of the ME values. To overcome this difficulty, one simple hypothesis has been introduced in [19]. This hypothesis is postulated by taking into consideration the fact that for a large number of data assimilation problems in HdSs, the model time step $\mathit{\delta t}$ (chosen for ensuring a stability of numerical scheme and for guaranteeing a high precision of the discrete solution) is much smaller than $\mathrm{\Delta }t$—the assimilation window (time interval between two successive observation arrivals).

Suppose that $\mathrm{\Delta }t=n_a\mathit{\delta t}$ where $n_a$ is a positive integer number.

Hypothesis (on the subspace of ME—HME) [19]. Under the condition that $n_a$ is relatively large, the ME belongs to the subspace spanned by all unstable and neutral EiVecs (or SchVecs) of the system dynamics $\Phi$.

5. Simple numerical examples

5.1 One-dimensional system

To see the difference between the AF and the KF in doing with ME uncertainties, introduce the one-dimensional system

$$x\left(t+1\right)=\mathit{\Phi x}\left(t\right)+w\left(t\right),z\left(t+1\right)=\mathit{hx}\left(t+1\right)+v\left(t+1\right),t=0,1,\dots$$E25

In (25), $\Phi$ is the unique eigenvalue (also the singular value) of the system dynamics.

1. For simplicity, let $\Phi =1,h=1$. This corresponds to the situation when the system is neutrally stable. The filter fundamental matrix (18) now is $L\left(K\right)=\left(1-K\right)$ which is stable if $K\in \left(0,2\right)$. For the KF gain (4)–(8), as $K_{\mathit{kf}}\left(t\right)=\frac{M_{\mathit{kf}}\left(t\right)}{M_{\mathit{kf}}\left(t\right)+R}$ we have $K_{\mathit{kf}}\left(t\right)\in \left(0,1\right)$, $M_{\mathit{kf}}\left(t\right)$ is the solution of (7). That is true for any $M_{\mathit{kf}}\left(t\right)\ge 0$, $R>0$. It means then the KF is stable. Mention that if $Q>0$ always $M_{\mathit{kf}}\left(t\right)>0$. In general, $K_{\mathit{kf}}\left(t\right)=M_{\mathit{kf}}\left(t\right){\left[M_{\mathit{kf}}\left(t\right)+R\right]}^{+}$ where ${\left[A\right]}^{+}$ is the pseudo-inverse of A [21].

   For the AF, we have in this case $P_r=1.$ Consider the gain $K_{\mathit{af}}\left(\theta \right)≔P_r\theta K_e$, where $K_e$ is the gain of the form $K_e=\frac{M_e}{M_e+R}$, $M_e>0$, $R>0$, $M_e$ is constant. We have then for the NAF ($\theta =1)$$0<K_e$ < 1 and $K_{\mathit{naf}}=K_e$.

   For the AF, the transition matrix (18) reads $L_{\mathit{af}}\left(\theta \right)=\left(1-\theta K_e\right)$. For $\theta \in \left(0,2\right)$, $\mid L_{\mathit{kf}}\left(\theta \right)\mid \in \left(0,1\right)$, $K_{\mathit{af}}\left(\theta \right)\in \left(0,2\right)$ and the AF is stable. It is evident that there is a larger margin for varying the gain in the AF than that in the KF since $K_{\mathit{kf}}\left(t\right)\in \left(0,1\right)$. One sees that the stationary KF is a member of the class of stable AFs (19). The performance of AF is optimized by solving the problem (9) using the procedure (11) and (12) or SPSA algorithms (Comment 2.2).

2. Let $\Phi <1$, i.e., the system (1) is stable. The results in (i) are valid for the AF structure. In this situation, the filter is stable even for $K_{\mathit{af}}=0$.

3. Let $⌈\Phi ⌉>1$—the system (1) is unstable. Consider two situations (a) $\Phi >1$; (b) $\Phi <-1$. As before $P_r=1$.

For $\Phi >1$we have

$$\frac{\Phi -1}{K_e\Phi }<\theta <\frac{\Phi +1}{K_e\Phi }$$E26

In particular, when $\Phi \to 1$, approximately $\theta \in \left(0,\frac{2}{K_e}\right)$. When $Q\gg R$ (that is usually in practice), approximately $\theta \in \left(0,2\right)$ as in the situation (i). For large $\Phi \gg 1$, $\frac{\Phi -1}{K_e\Phi }\to \frac{1}{K_e}$ (left-hand limit), $\frac{\Phi +1}{K_e\Phi }\to \frac{1}{K_e}$ (right-hand limit) and there remains no margin for varying $\theta$ (or $Q\gg R$) and $K_{\mathit{af}}\to 1$. It is important to emphasize that as $K_e$ is chosen by designer, we can define the interval for varying $\theta$ if the amplitude of $\Phi$ is more or less known. In practice, one can vary $\theta \in \left[ϵ,2+ϵ\right]$ with small $ϵ>0$ for $\Phi$ close to 1, and with $ϵ$close to 1 for large $\Phi$.

For $\Phi <-1$we have

$$\frac{\Phi +1}{K_o\Phi }<\theta <\frac{\Phi -1}{K_o\Phi }$$E27

It is seen from (27) that when $\Phi \to -1$, approximately $\theta \in \left(0,\frac{2}{K_e}\right)$. As for the situation $\Phi \ll -1$, $\frac{\Phi +1}{K_o\Phi }\to \frac{1}{K_o}$ (left-hand limit), $\frac{\Phi -1}{K_o\Phi }\to \frac{1}{K_o}$ (right-hand limit) when $Q\gg R$ approximately $\theta \to \frac{1}{K_o}$hence $K_{\mathit{af}}\to 1$.

It is important to stress that the KF gain is computed on the basis of $Q$ and $R$ (under the condition that the statistics of the initial state will be forgotten as $t$ becomes large); whereas, the gain of the AF is updated on the basis of samples of the innovation vector. It means that the KF is optimal in the MMSE sense (under the condition of exact knowledge of the required statistics) whereas the AF is optimized during the assimilation process using PE realizations of the system output (innovation vector). The KF gain can be computed in an offline fashion, whereas the AF gain is a function of observation and computed in online.

5.2 Two-dimensional system: specification of covariances

5.2.1 Stable filter

To see the role of the correction subspace $R\left[P_r\right]$ in ensuring a stability of the AF, let us consider the system (1) and (2) s.t.

$$\Phi =\mathit{diag}\left({\Phi }_{11},{\Phi }_{22}\right),H=\mathit{diag}\left(1,1\right)$$E28

Consider the AF with the gain (19) s.t. $P_r=\mathit{Id}$, i.e., two columns of $P_r$ are in fact the EiVecs associated with two EiVs ${\Phi }_{11}$ and ${\Phi }_{22}$. Let us denote the AF gain $K_{\mathit{af}}=P_r\mathrm{\Theta }M{H}^T{\left[\mathit{HM}{H}^T+R\right]}^{-1}=\mathrm{\Theta }{K}_e$ with $K_e≔M_e{H}^T{\left[H{M}_e{H}^T+R\right]}^{-1}$ which is structurally identical to that of the KF. For the nonadaptive filter $K_{\mathit{naf}}=K_e$ and with the choice $M_e=\mathit{diag}\left(M_{11},M_{22}\right)$, taking into account (28) one gets

$$K_{\mathit{naf}}=\mathit{diag}\left(K_{11},K_{22}\right),K_{\mathit{ii}}=\frac{M_{\mathit{ii}}}{M_{\mathit{ii}}+R_i},i=1,2$$E29

$$L_{\mathit{naf}}=\mathit{diag}\left(l_{11},l_{22}\right),l_{\mathit{ii}}=\left(1-\frac{M_{\mathit{ii}}}{M_{\mathit{ii}}+R_i}\right){\Phi }_{\mathit{/ii}},i=1,2$$E30

The filter transition matrix (30) is obtained on the basis of $L_{\mathit{naf}}=\left(I-{\mathit{K}}_{\mathit{naf}}H\right)\Phi$ and the assumption (29). It is easily to see that $L_{\mathit{naf}}$ has two EiVs, ${\lambda }_i=l_{\mathit{ii}},i=1,2$ where $l_{\mathit{ij}}$ est. the $\left(\mathit{ij}\right)$ element of $L_{\mathit{naf}}$.

Stability of the filter depends on the condition $\mid l_{\mathit{ii}}\mid <1$, $i=1,2$. For $M_{\mathit{ii}}>0,R_i>0$, if ${\Phi }_{\mathit{ii}}$ is stable or neutrally stable, i.e. $\mid {\Phi }_{\mathit{ii}}\mid \le 1,i=1,2$, we have $\mid l_{\mathit{ii}}\mid <1$. For unstable${\Phi }_{\mathit{ii}}$ ($\mid {\Phi }_{\mathit{ii}}\mid >1,i=1,2)$, the filter is unstable if ${\Phi }_{\mathit{ii}}>\frac{M_{\mathit{ii}}+R_i}{R_i}$ (situation ${\Phi }_{\mathit{ii}}>1$) or ${\Phi }_{\mathit{ii}}<-\frac{M_{\mathit{ii}}+R_i}{R_i}$ (situation ${\Phi }_{\mathit{ii}}<-1$). These conditions should be taken into account when the EiVs of $\Phi$ are large.

For the AF gain (19) ($P_r=I$),

$$l_{\mathit{ii}}=\left(1-\frac{{\theta }_i{M}_{\mathit{ii}}}{M_{\mathit{ii}}+R_i}\right){\Phi }_{\mathit{ii}},$$E31

From (31) conditions for $\mid l_{\mathit{ii}}\mid <1$ can be obtained as done in Section 5.1 with the one-dimensional system since $l_{\mathit{ii}},i=1,2$ are independent one from another. The length of the interval $I_i$for varying ${\theta }_i$ depends on the value of ${\Phi }_{\mathit{ii}}$ (see (26)).

This example shows that for $P_r=\mathit{Id}$, it is always possible to construct a stable AF whatever are the EiVs of $\Phi$ (stable or unstable). There are some constraints for $M_{\mathit{ii}}$ (they are positive) and for $R_i$ (small positive). Optimality of the AF is obtained by searching recursively (in time) the optimal ${\theta }_i$ during assimilation process. Thus, in the AF, a correct specification of ME and ObE statistics (second order) is not important as happens in the KF.

5.2.2 Unstable filter

Consider the situation when $P_r$ is constructed from only one vector. Let $P_r={\left(1,0\right)}^T$—the EiVec associated with ${\Phi }_{11}$ (the results remain the same if we choose $P_r={\left(0,1\right)}^T$—the EiVec associated with ${\Phi }_{22}$).

$$\Phi =\mathit{diag}\left({\Phi }_{11},{\Phi }_{22}\right),\mid {\Phi }_{11}\mid >1,\left|{\Phi }_{22}\right|>1,H=\mathit{diag}\left(1,1\right)$$E32

We show now that the filter with the gain (19) is unstable. We have (for $\Theta =\mathit{Id}$),

$$H_e=H{P}_r=P_r,K=P_r{K}_e,K_e=M_e{H}_e^T{\left[H_e{M}_e{H}_e^T+R\right]}^{-1}=M_e{P}_r{\left[P_r{M}_e{P}_r^T+\mathit{\alpha I}\right]}^{-1}$$E33

if we put $R=\mathit{\alpha I}$. For $L_{\mathit{naf}}=\left[I-\mathit{KH}\right]\Phi$, taking into account (32), (33) it implies

$$L_{\mathit{naf}}=\left[\begin{array}{cc}\frac{\alpha {\Phi }_{11}}{m_e+\alpha }& 0\\ 0& {\Phi }_{22}\end{array}\right]$$E34

As $\frac{\alpha }{m_e+\alpha }$ can be made as small as desired by choosing small $\alpha >0$, the first EiV $l_{11}=\frac{\alpha {\Phi }_{11}}{m_e+\alpha }$ can be made stable. However, the second EiV in (34)$l_{22}={\Phi }_{22}>1$ is unstable. It implies that the filter with the gain (19) s.t. $P_r={\left(1,0\right)}^T$ is unstable. This happens even for $\Theta \ne \mathit{Id}$. It means that when the projection subspace $R\left[P_r\right]$ does not contain all unstable and neutral EiVecs of the system dynamics, it is impossible to guarantee a stability of the filter.

5.3 Two-dimensional system: estimation of ME

Consider the filtering problem (1) and (2), the dynamical system (1) describes a sequence of system states at time instants $t=0,1,\dots$ when the observations are available. It means that $\Phi$ represents the transition of system state over the (observation) time window $\mathrm{\Delta }t$ separating arrivals of two successive observations. In practice, the interval $\mathrm{\Delta }t$ is much larger than the model time step$\mathit{\delta t}$ which is the step size in approximating the temporal derivative. The choice of $\mathit{\delta t}$ is important for guaranteeing a stability of discretized scheme and having high is important for guaranteeing a stability of discretized scheme and having high precision of the discretized solution (wrt the continuous solution). We have then $\mathrm{\Delta }T=n_a\mathit{\delta t}$, where $n_a$ is a relatively large positive integer. For example, in the HYCOM model at SHOM (French marine) for the Bay of Biscay configuration, the interval $\mathrm{\Delta }t$ between two observation arrivals is 7 days which is equivalent to integrating 1200 model time steps $\mathit{\delta t}$. It means $n_a=1200$. Symbolically we have then the equations for model time step integration

$$x\prime \left(\tau +1\right)={\Phi }^{\prime }x\prime \left(\tau \right)+{\psi }^{\prime }\left(\tau \right),{\psi }^{\prime }\left(\tau \right)≔b^{\prime }\left(\tau \right)+w\prime \left(\tau \right)$$E35

In (35), $\Phi \prime$ represents the integration of numerical model over one model time step $\mathit{\delta t}$. Hence

$$\Phi ={\left[{\Phi }^{\prime }\right]}^{n_a}$$E36

The contribution of ${\psi }^{\prime }\left(\tau \right),$ over the assimilation window $\left[t-1,t\right]$ (for simplicity and without loss of generality, one supposes $t-1≔0,t≔n\_a$) is

$$\psi \left(t\right)=b\left(t\right)+w\left(t\right),b\left(t\right)≔\sum _{\tau =0}^{n_a}{\mathrm{\nu }}_1{\left(\tau \right)}_1,{\mathrm{\nu }}_1\left(\tau \right)≔{\left[{\Phi }^{\prime }\right]}^{n_a-1-\tau }b\prime \left(\tau \right),$$

$$w\left(t\right)≔\sum _{\tau =0}^{n_a}{\mathrm{\nu }}_2\left(\tau \right),{\nu }_2\left(\tau \right)≔{\Phi }^{n_a-1-\tau }{w}^{\prime }\left(\tau \right),{\left[{\Phi }^{\prime }\right]}^{-1}≔0$$E37

The HME in Section 4 says that the SME $w\left(t\right)$ and DME $b\left(t\right)$, as functions of $n_a$, belong to the subspace spanned by leading EiVs (or SchVecs) of $\Phi$ for a relatively large $n_a$. The initial filtering problem now has the form (1) and (2) s.t. (36) and (37) where $t=\tau /n_a$.

To illustrate this HME, continue the two-dimensional system in Section 5.2.2 and suppose that $\mid \Phi {\prime }_{11}\left|>1,|\Phi {\prime }_{22}\right|<1$. Applying HME in this case is equivalent to saying that the values of MEs $b\left(t\right),w\left(t\right)$, approximately, belong to the subspace $R\left[u_1\right]$spanned by the first EiVec $u_1=\mathit{col}\left(1,0\right)$, associated with the EiV $\Phi {\prime }_{11}$. Here $y=\mathit{col}\left(y_1,\dots ,y_n\right)$ denotes the vector-column with components $y_1,\dots ,y_n$. It follows that the covariance matrix of $w\left(t\right)$ is assumed to be of the form $Q={\sigma }_w^2{u}_1{u}_1^T$ and $b\left(t\right)$—of the structure $b\left(t\right)=c{u}_1$, $c$ is a scalar to be estimated. For the algorithm of joint estimation of state and bias (in term of $c$), see [19].

6. Simulation results

6.1 One-dimensional system

In this section, the filtering problem (25) in Section 5.1 is considered s.t.

$$\Phi =0.99,H=1.02,Q=0.09,R=0.01.$$

The true system states and observations are simulated using the initial state $x\left(0\right)=1$ and $w\left(t\right),v\left(t\right)$ are zero mean Gaussian mutually uncorrelated and temporal uncorrelated sequences.

To see the performance of the AF, unknown system states are estimated on the basis of the AF algorithm. To obtain a reference, the standard KF is also implemented for solving this filtering problem. In the filtering algorithms, the estimate of the initial state is $\widehat{x}\left(0\right)=2.$ The gain $K_{\mathit{naf}}$ in the NAF is taken as that of the KF at $t=0$, i.e., $K_{\mathit{naf}}=K_{\mathit{kf}}\left(0\right)$.

Figure 1 shows the temporal evolution of the parameters ${\theta }_m\left(t\right)$ during assimilation process.

Figure 1.

Temporal evolution of the parameter ${\theta }_m\left(t\right)$ in the AF gain.

The gains in the KF and AF during the assimilation process are displayed in Figure 2. Mention that the KF gain is computed s.t. true statistics $Q,R$. In the AF, ${\theta }_m\left(t\right)$ has been used for computation of the AF gain, i.e., $K_{\mathit{af}}={\theta }_m\left(t\right)K$. From Figure 2, one sees that initialized by the same value, the two gains become different during assimilation process. The KF gain has reached a stationary regime very quickly.

Figure 2.

Temporal evolution of gains in KF and AF during data assimilation.

The mean temporal RMS (root mean square) of the innovation is shown in Figure 3. It is interesting to remark that no significant difference is observed between two curves and a slightly better performance is produced by the KF.

Figure 3.

RMS of innovation produced by the KF and AF.

In Figure 4, we show RMS of the state FE produced by the KF and AF under the condition that the variance $Q$ is known exactly. One sees that the KF, as expected, produces the best results.

Figure 4.

RMS of the state FE produced by the KF and AF under the condition that the variance of ME is known exactly. It is seen that when the ME is correctly specified, the KF behaves better than the AF.

Figure 5 shows the RMS of FE as a function of the variance $Q$. Here, the value of $Q$ varies from 0.1 to 1.9. Note that the true value of $Q$ is 0.1. The red curve represents the RMS of FE produced by the KF at the end of the assimilation period (as a function of $Q$). The green curve has the same meaning, but for the FE produced by the AF. It is interesting to note that when $Q$ is correctly specified, the KF behaves better than the AF, but misspecification of $Q$ leads to growing of the error in the KF. The AF is robust wrt the error in the specification of $Q$. This fact says in favor of the AF as an efficient tool for overcoming uncertainties in the ME.

Figure 5.

RMS of FE as a function of $Q$. The true value of $Q$ is equal to 0. It is noted that the KF behaves better than the AF s.t. true Q but is more and more degraded as the ME becomes greater and greater. At the same time, the FE of the AF remains very robust.

6.2 Two dimensional system

6.2.1 Illustration of hypothesis HME

According to the notations in Section 5.3, consider the two-dimensional system (1) with $\Phi \prime =\left[\Phi {\prime }_{\mathit{ij}}\right],\Phi {\prime }_{11}=1.02,\Phi {\prime }_{12}=0.1,\Phi {\prime }_{21}=0,\Phi {\prime }_{22}=0.9,H=\left(1,1\right)$ with the true DME $b\prime \left(\tau \right)=\mathrm{col}\left(\mathrm{0.1,0.1}\right)\prime .$ Thus the first EiV is unstable, the second—stable [19].

Numerically one finds that the first SchVec is equal to $u_1={\left(-1,-7.0E-7\right)}^{\mathrm{T}}$.

Figure 6 [19] shows the simulation results obtained on the basis of (37). One sees that, for $n_a>10$, the second component of $w\left(t\right)$ is close to $0$ whereas the first component becomes bigger and bigger (in absolute value) as $n_a$ increases. Here, $w^{\prime }\left(\tau \right)$ is a sequence of independent two-dimensional Gaussian random vectors of zero mean and variance 1. This means that the values of $w\left(t\right)$ become more and more close to the subspace $R\left[u_1\right]$ spanned by $u_1$, hence the HME is practically valid for $n_a>10$ in this example. Mention that, as a rule, in ocean numerical models, $n_a$ is of order o($100$) ($n_a=800$ or the MICOM model in the experiment in Section 6.3). See also [22].

Figure 6.

Two components of $w\left(t\right)$ as functions of $n_a$.

6.2.2 Assimilation results

First simulation of the sequence of true system states $x\prime \left(\tau \right),\tau =1,\dots ,390$ has been carried out (see (35)). The observations are picked at $\tau =15,30,\dots ,390$. In terms of $x\left(t\right)$, the filtering problem then is of the form (1) and (2) s.t. $t=\frac{\tau }{15},\Phi =\Phi {\prime }^{15}$ (if no bias exists). When there is a bias, instead of $w\left(t\right)$ stands $\psi \left(t\right)≔b\left(t\right)+w\left(t\right).$

In the experiment, the true system states $x\prime \left(t\right)$ are generated by (35) s.t. $b\prime \left(\tau \right)=\mathit{col}\left(\mathrm{0.1,0.1}\right),w\prime \left(t\right)$ is zero mean with the covariance $Q^{\prime }=\mathit{Id}$. The observation error is of zero mean and covariance $R=0.16$. For the state $x\left(t\right)$ in the KF and NAF, the forecast is obtained at each assimilation instant $t$ as $\widehat{x}\left(t+1/t\right)=\Phi \widehat{x}\left(t\right)+\widehat{b}\left(t\right)$. The simulation yields $b\left(t\right)=\left(0.2296,2.0589\mathrm{E}-02\right)$ which results from applying (37) s.t. $b\prime \left(\tau \right)=\mathit{col}\left(\mathrm{0.1,0.1}\right)$.

Figure 7 depicts the time evolution of the KF and AF gains. One sees here as in the experiment with 1D system (Figure 2) that the KF gain is stabilized very quickly compared to that of the AF gain.

Figure 7.

Gain coefficients in KF and AF: The gains in KF and AF are identical at the beginning of the assimilation process.

Figure 8 (from [19]) shows the sample time average RMS of the state FE produced by the three filters NAF, KF, and AF. One sees that the AF outperforms the NAF and KF.

Figure 8.

Sample time average RMS of the state filtered error produced by the NAF, KF, and AF.

6.3 Data assimilation in the high-dimensional ocean model

To illustrate the effectiveness of the AF in dealing with uncertainties in HdSs, this section presents the results on data assimilation in the oceanic numerical model MICOM (Miami Isopycnic Coordinate Ocean Model) [19]. This MICOM describes the oceanic circulation in the North Atlantics. The model has four vertical layers with the state consisting of three variables $x=\left(h,u,v\right)$ where $h$ is layer thickness, $\left(u,v\right)$ are two velocity components. The horizontal grid is $140\times 180$. Totally at each time instant $t$ we have the state $x\left(t\right)$ of dimension $302400\left(140\times 180\times 4\left(\mathit{\text{layers}}\right)\times 3\left(\mathit{\text{variables}}\right)\right)$. For more details on the configuration of this experiment, see [22].

The experiment is carried out on estimating the oceanic circulation using sea surface height (SSH) measurements. The SSH observation is available each 7 days (ds) (hence the observation window $\mathrm{\Delta }T=7\mathit{ds}$). Mention that simulating the circulation over 7 ds requires $800$ model time steps $\left(\mathit{\delta t}\right)$integration.

6.3.1 AF with optimal initial gain

First, in order to examine whether the method of optimal gain initialization, described in Section 3.2, is really useful for improving the filter performance, the optimization problem (21) has been solved. Symbolically, in the gain (20), $P_r={\left[\mathit{Id},{\left(\mathit{QG}\right)}^T\right]}^T$ where $\mathit{Id}$ is the identity operator on the space of layer thickness $h$, $\mathit{QG}$ is the quasi-geostrophy operator computing the correction for velocity using the SSH innovation $\zeta$. The gain $K_e$ computes the correction for $h$ using $\zeta$. The ECM $M=M_v\otimes M_h$—Kronecker product of $M_h$—ECM of horizontal variable, $M_v$—ECM with vertical variable (see below for details). The two parameter matrices $\mathrm{\Theta }$ and $\mathrm{\Lambda }$ are related to parameterization of $M_v$. The problem (21) is solved s.t. $\mathrm{\Theta }=\mathit{Id}$—identity operator.

The optimal parameters ${\lambda }_i,i=1,\dots ,4$ are found by solving the minimization problem (21) using SPSA algorithm. Figure 9 shows the averaged values (see Comment 2.1) of ${\lambda }_i,i=1,\dots ,4$ resulting during the optimization process. All ${\lambda }_i,i=1,\dots ,4$ are initialized as ${\lambda }_i=1,i=1,\dots ,4$.

Figure 9.

Control parameters ${\lambda }_i,i=1,\dots ,4$ during optimization process.

The two NAFs are performed, one (denoted as NAFI) is with the gain (20) s.t. $\mathrm{\Theta }=\mathit{Id}$, $\mathrm{\Lambda }=\mathit{Id}$, and the other (denoted by NAFOI) s.t. $\mathrm{\Theta }=\mathit{Id}$ and ${\lambda }_i,i=1,\dots ,4$ obtained by solving (21) (their values are those displayed at the end of the optimization process in Figure 9). The performances of these two NAFs are shown in Figure 10. One sees here that the NAFOI has improved considerably the quality of estimates of the velocity u-component compared with the NAFI. This result justifies that offline optimization (21) is an interesting strategy for finding the optimal initial gain in the NAF.

Figure 10.

Temporal average RMS of FE for total velocity $u$-component produced by the NAF s.t. the initial gain (red curve) and the NAF s.t. optimal initial gain resulting from solving (22) (green curve).

6.3.2 Estimating the ECM of ME

In practice, for real operational systems, information on the space of ME is not available or very poorly known. Usually, there is a big difference between the model and the real physical process and if the ME statistics are taken more or less properly, in some way, in the filtering algorithm, one can improve the filter performance and reduce the estimation error.

This idea is tested here by applying the HME in Section 4. We carry out the procedure for estimating the ECM of the ME by first constructing the subspace for the ME. For more details on the structure of the ECM $M$ in the AF, see [23]. According to [23], the ECM $M$ is assumed to be of the structure $M=M_v\otimes M_h$—the Kronecker product of $M_h$ with $M_v$ where $M_h$ is the ECM of the horizontal variable, $M_v$—ECM with vertical variable. Figure 11 displays RMS of FE for the $u$ velocity component at the surface resulting from two AFs. The curve $\mathit{AF}0U$ corresponds to the AF whose nonadaptive version has the gain computed on the basis of the ECM $M$ using an ensemble of PE samples (generated by the PeSP in [18]). The curve $\mathit{AF}3U$shows the performance of the AF with the modified ECM (by adding the vertical ME covariance $Q_v$ to the vertical ECM $M_v$). More precisely, $Q_v$ is assumed to belong to the subspace spanned by three leading EiVs of $M_v$. This choice is justified by the fact: the eigenvalue decomposition of $M_v$ has the first three EiVs with the explained variances 67, 17, 15%, respectively. As the fourth EiVec has only the explained variance 0.7E-07%, it is dropped from the subspace constructed for the vertical ME. The better performance of the $\mathit{AF}3U$, in comparison with that of the $\mathit{AF}0U$, is apparently seen in Figure 11.

Figure 11.

Performance of the AF: (i) $\mathit{AF}0U$—no ME ECM has been taken into account; (ii) $\mathit{AF}3U$—with ME ECM computed in accordance with the HME.

The above experiment shows in details how, on the basis of HME, the subspace for the ME can be constructed, and how one estimates the ECM for the model error. The superior performance of the $\mathit{AF}3U$ over that of the $\mathit{AF}0U$ validates the usefulness of the HME which can serve as an important tool for estimating the ME and improving the performance of the AF for solving the data assimilation problems with HdSs.

7. Conclusions

One of the key assumptions to ensure the optimal performance of the KF is that a priori knowledge of the system model is given without any uncertainty. This assumption, however, is never valid in practice for dynamical systems under consideration. The uncertainties exist everywhere in modeling a real process like structural uncertainty, model parameterization, model resolution, model bias or ME statistics. For HdSs, order reduction introduced either in the original numerical model or in the filtering algorithms, inevitably leads to uncertainty in the ME, especially in geophysical numerical models.

Our focus in this chapter is to show how the AF solves efficiently filtering problems for systems operating in an uncertain environment.

As seen from this chapter, the AF has proven to be efficient to deal with uncertainties in the specification of the ME statistics, system bias or model reduction. The reasons of the success of the AF are that (i) it belongs to the class of parametrized stable filters; (ii) it is defined as the best member minimizing mean PE for the system outputs; (iii) The tuning parameters are chosen as elements of stabilizing gain and they are of no physical sense.

It is obvious from this chapter that the performance of the AF is comparable with that of the KF when perfect knowledge of all ME statistics is given, and it outperforms the KF in presence of uncertainties. This happens since the AF acquires knowledge during assimilation process, regardless of uncertainties existing in the filtering problems. From the computational point of view, implementation of the AF consumes much less memory and computational time than the KF or other assimilation methods.

Simple numerical examples and simulation results, presented in Sections 5 and 6, clearly demonstrate the advantages gained through application of the AF in dealing with uncertainties. These positive results encourage a wide application of the AF in different fields of technology and applied sciences like automatic control, finance, aerospace, space exploration, meteorology, and oceanography. A more in-depth and significant research on the capacity of the AF to deal with uncertainties is surely a challenge for the near future.