Estimation of Shallow Water Flow Based on Kalman Filter FEM

Abstract

In this chapter, we present numerical examples of an estimation of shallow water flow based on Kalman filter finite element method (Kalman filter FEM). Shallow water equations are adopted as the governing equations. The Galerkin method, using triangular elements, is employed to discretize the governing equation in space, and the selective lumping method is used to discretize time. We describe the influence on the numerical results of setting the observation points.

1. Introduction

Conventionally, shallow water flow analysis based on the FEM has been carried out to investigate the marine environment, specifically coastal drift sand, storm surges, tsunamis, and so on [1, 2]. However, the computed results are rarely close to the observed values if the appropriate governing equation is not employed, the appropriate discretization technique is not applied to the governing equation, or the boundary conditions are not defined appropriately. This has prompted a great deal of inverse analysis in recent years, aiming to obtain flow fields that more closely approach the observed values. Boundary control analysis is one of these inverse analytic techniques, in which an unknown boundary condition is numerically determined by iterative computation such that the computed water elevation homes in on the target value or the observed value [3, 4]. Deterministic methods such as the adjoint variable or sensitivity equation methods are commonly employed. On the other hand, the stochastic approach has also been adopted, in which estimation of the flow field is employed using a Kalman filter. Removal of the noise data from the observed values and the flow field estimation are also performed. Conventionally, flow fields are estimated using Kalman filter FEM, and this method is applied to the practical coastal model [5, 6]. However, the effect of boundary conditions on the estimated flow field has not been investigated. This prompted us to research and describe in detail in this chapter the formulation of Kalman filter FEM and numerical experiments for flow estimation when changing the positions of observation points [7].

2. Concept of state estimation

In the Kalman filter, system and observation equations are employed. Eqs. (1) and (2) are sample equations:

${\varphi }^{n+1}=a{\varphi }^n+q^n$E1

$z^{n+1}={\varphi }^{n+1}+r^{n+1}$E2

where ϕ, z, a, q, and r indicate the true value of the state variable, the observation value, the constant coefficient, the system, and the observation noise, respectively. n denotes the number of time steps. Figure 1 shows relationship between the observation and the estimation values.

Figure 1.

Diagram of the relationship between the observation and the estimation values.

Here, $\widehat{\varphi }$ and p indicate the estimation value for the state variable and estimation error. From Figure 1, it is found that the estimation value is expressed as a summation of the true value and the estimation error, and the observation value is expressed by the summation of the true value and the observation error. Based on this concept, our formulation of the Kalman filter is shown below.

3. Discretization of the governing equation

The shallow water equations, shown below, are employed as the governing equations.

${\dot{u}}_i+g\eta {,}_i=0$E3

$\dot{\eta }+h{u}_{i,i}=0$E4

where u_i, η, g, and h indicate flow velocity component for the x and y directions, water elevation, gravitational acceleration, and water depth, respectively. The Galerkin method and the selective lumping method are applied to discretize the governing equation. The discretized governing equation is written as Eq. (5) and is represented by Eq. (6),

$\left\{\begin{array}{l}\left\{u^{n+1}\right\}\\ \left\{v^{n+1}\right\}\\ \left\{{\eta }^{n+1}\right\}\end{array}\right\}={\left[\begin{array}{ccc}\left[\overline{M}\right]& \left[0\right]& \left[0\right]\\ \left[0\right]& \left[\overline{M}\right]& \left[0\right]\\ \left[0\right]& \left[0\right]& \left[\overline{M}\right]\end{array}\right]}^{-1}\left[\begin{array}{ccc}\left[\tilde{M}\right]& \left[0\right]& -g\Delta t\left[S_x\right]\\ \left[0\right]& \left[\tilde{M}\right]& -g\Delta t\left[S_y\right]\\ -\overline{h}\Delta t\left[S_x\right]& -\overline{h}\Delta t\left[S_y\right]& \left[\tilde{M}\right]\end{array}\right]\left\{\begin{array}{c}\left\{u^n\right\}\\ \left\{v^n\right\}\\ \left\{{\eta }^n\right\}\end{array}\right\}$E5

$\left\{{\widehat{\varphi }}^{n+1}\right\}=\left[A\right]\left\{{\widehat{\varphi }}^n\right\}$E6

where Δ, $\overline{h}$, [S_i], [M], and $\left[\overline{M}\right]$ denote the time increment, the mean water depth, the matrix for the pressure, the consistent mass matrix, and diagonal mass matrix, respectively. $\left[\tilde{M}\right]$ is represented by Eq. (7) by using the lumping parameter e(0 ≤ e ≤ 1). $\left\{\widehat{\varphi }\right\}$ indicates the estimation value.

$\left[\tilde{M}\right]=\left(1-e\right)\left[\overline{M}\right]+e\left[M\right]$E7

The system equation in the Kalman filter is obtained by adding the vector represented by the multiplication of the driving matrix [Γ] and the system noise vector {g} and is shown as Eq. (8). Here, the vector {ø} indicates the vector of true value. The observation equation shown in Eq. (9) is also introduced. The vectors {z} and {r} and the matrix [H] denote the observation value and the observation noise vectors and the observation matrix.

$\left\{{\varphi }^{n+1}\right\}=\left[A\right]\left\{{\varphi }^n\right\}+\left[\Gamma \right]\left\{q^n\right\}$E8

$\left\{z^{n+1}\right\}=\left[H\right]\left\{{\varphi }^{n+1}\right\}+\left\{r^{n+1}\right\}$E9

4. Derivation of computational equations in the Kalman filter FEM

The process of modification of the estimation vector $\left\{\widehat{\varphi }\right\}$ is referred to as “assimilation,” and the estimation vector before and after assimilation is expressed by “(−)” and “(+)” marks, respectively. First of all, it is assumed that the estimated value after assimilation is expressed by the observation value and the estimated value before assimilation, and that the estimated value after assimilation, i.e., the optimal estimated vector, is expressed by Eq. (10).

$\left\{{\widehat{\varphi }}_{(+)}^{n+1}\right\}=\left[K_1^{n+1}\right]\left\{z^{n+1}\right\}+\left[K_2^{n+1}\right]\left\{{\widehat{\varphi }}_{(-)}^{n+1}\right\}$E10

where $\left[K_1^{n+1}\right]$ and $\left[K_2^{n+1}\right]$ indicate arbitrary matrices. In addition, the optimal estimated vector $\left\{{\widehat{\varphi }}_{(+)}^{n+1}\right\}$ is expressed by true value $\left\{{\varphi }^{n+1}\right\}$ and estimation error $\left\{p_{(+)}^{n+1}\right\}$ as shown in Eq. (11).

$\left\{{\widehat{\varphi }}_{(+)}^{n+1}\right\}=\left\{{\varphi }^{n+1}\right\}+\left\{p_{(+)}^{n+1}\right\}$E11

Consequently, Eq. (12) is obtained from Eqs. (9), (10), and (11).

$\begin{array}{l}\left\{{\varphi }^{n+1}\right\}+\left\{p_{(+)}^{n+1}\right\}=\\ \left[K_1^{n+1}\right]\left(\left[H\right]\left\{{\varphi }^{n+1}\right\}+\left\{r^{n+1}\right\}\right)\\ +\left[K_2^{n+1}\right]\left(\left\{{\varphi }^{n+1}\right\}+\left\{p_{(-)}^{n+1}\right\}\right)\end{array}$E12

Here, we consider the expected value for Eq. (12) in time. In addition, the expected value for vectors $\left\{r^{n+1}\right\}$, $\left\{p_{(+)}^{n+1}\right\}$ and $\left\{p_{(-)}^{n+1}\right\}$ in time is assumed to be zero, and Eq. (12) is represented as Eq. (13).

$\left(\left[K_2^{n+1}\right]+\left[K_1^{n+1}\right]\left[H\right]-\left[I\right]\right)⟨\left\{{\varphi }^{n+1}\right\}⟩=\left\{0\right\}$E13

where [I] denotes the unit matrix, and $⟨\left\{{\varphi }^{n+1}\right\}⟩$ indicates the expected value of $\left\{{\varphi }^{n+1}\right\}$ in time. Because of $⟨\left\{{\varphi }^{n+1}\right\}⟩\ne \left\{0\right\}$, the following relation equation is obtained (see Eq. (14)).

$\left[K_2^{n+1}\right]+\left[K_1^{n+1}\right]\left[H\right]-\left[I\right]=\left[0\right]$E14

Substituting Eq. (14) for Eq. (10), Eq. (15) is obtained.

$\left\{{\widehat{\varphi }}_{(+)}^{n+1}\right\}=\left\{{\widehat{\varphi }}_{(-)}^{n+1}\right\}+\left[K_1^{n+1}\right]\left(\left\{z^{n+1}\right\}-\left[H\right]\left\{{\widehat{\varphi }}_{(-)}^{n+1}\right\}\right)$E15

where $\left[K_1^{n+1}\right]$ is referred to as the Kalman gain matrix.

Next, let us consider the Kalman gain matrix $\left[K_1^{n+1}\right]$ so as to minimize the trace norm of the error covariance matrix after assimilation. The estimation vector is represented by Eq. (16).

$\left\{p_{(+)}^{n+1}\right\}=\left\{{\widehat{\varphi }}_{(+)}^{n+1}\right\}-\left\{{\varphi }^{n+1}\right\}$E16

Eq. (17) is obtained by substituting Eq. (15) with Eq. (16), and substituting Eq. (9) with the obtained equation.

$\begin{array}{l}\left\{p_{(+)}^{n+1}\right\}=\\ \left(\left[I\right]-\left[K_1^{n+1}\right]\left[H\right]\right)\left\{p_{(-)}^{n+1}\right\}+\left[K_1^{n+1}\right]\left\{r^{n+1}\right\}\end{array}$E17

Here, the error covariance matrix after assimilation $\left[P_{(+)}^{n+1}\right]$, i.e., the predicted error covariance matrix, is expressed by the expected value of the error covariance matrix as shown in Eq. (18).

$\left[P_{(+)}^{n+1}\right]=⟨\left\{p_{(+)}^{n+1}\right\}{\left\{p_{(+)}^{n+1}\right\}}^T⟩$E18

Also, in calculating Eq. (18), it is assumed that the expected value of the covariance matrix of the observation error and the estimation error before assimilation is a zero vector (see Eq. (19)).

$⟨\left\{r^{n+1}\right\}{\left\{p_{(-)}^{n+1}\right\}}^T⟩=\left[0\right]$E19

Consequently, the predicted error covariance matrix can be expressed, such as in Eq. (20),

$\begin{array}{l}\left[P_{(+)}^{n+1}\right]=\\ \left(\left[I\right]-\left[K_1^{n+1}\right]\left[H\right]\right)\left[P_{(-)}^{n+1}\right]{\left(\left[I\right]-\left[K_1^{n+1}\right]\left[H\right]\right)}^T\\ +\left[K_1^{n+1}\right]\left[R^{n+1}\right]{\left[K_1^{n+1}\right]}^T\end{array}$E20

where $\left[R^{n+1}\right]$ is the observation error covariance matrix, and is represented by $⟨\left\{r^{n+1}\right\}{\left\{r^{n+1}\right\}}^T⟩$. Here, let us consider the Kalman gain matrix $\left[K_1^{n+1}\right]$ so as to minimize the trace norm of the predicted error covariance matrix such as $\frac{\partial }{\partial \left[K_1^{n+1}\right]}\text{tr}\left[P_{(+)}^{n+1}\right]=0$. Consequently, we can obtain the equation for calculating the Kalman gain matrix, as shown in Eq. (21).

$\begin{array}{l}\left[K_1^{n+1}\right]=\\ \left[P_{(-)}^{n+1}\right]{\left[H\right]}^T{\left(\left[H\right]\left[P_{(-)}^{n+1}\right]{\left[H\right]}^T+{\left[R\right]}^{n+1}\right)}^{-1}\end{array}$E21

In addition, Eq. (20) is represented by Eq. (22) taking Eq. (21) into account.

$\left[P_{(+)}^{n+1}\right]=\left[P_{(-)}^{n+1}\right]-\left[K_1^{n+1}\right]\left[H\right]\left[P_{(-)}^{n+1}\right]$E22

Next, let us consider the computation of the estimated error covariance matrix before assimilation $\left[P_{(-)}^{n+1}\right]$, i.e., the estimated error covariance matrix. First of all, Eq. (6) is represented as Eq. (23) based on the consideration expressed in Eq. (11).

$\left\{{\varphi }^{n+1}\right\}+\left\{p_{(-)}^{n+1}\right\}=\left[A\right]\left(\left\{{\varphi }^n\right\}+\left\{p_{(-)}^n\right\}\right)$E23

Subtracting Eq. (8) from Eq. (23) gives us Eq. (24).

$\left\{p_{(-)}^{n+1}\right\}=\left[A\right]\left\{p_{(+)}^n\right\}-\left[\Gamma \right]\left\{q^n\right\}$E24

Here, the estimated error covariance matrix $\left[P_{(-)}^{n+1}\right]$ is also represented by the expected value of the error covariance matrix, as shown in Eq. (25).

$\left[P_{(-)}^{n+1}\right]=⟨\left\{p_{(-)}^{n+1}\right\}{\left\{p_{(-)}^{n+1}\right\}}^T⟩$E25

Considering the relation equations $⟨\left\{p_{(+)}^n\right\}{\left\{q^n\right\}}^T⟩=\left[0\right]$ and $⟨\left\{q^n\right\}{\left\{p_{(+)}^n\right\}}^T⟩=\left[0\right]$ in the calculation in Eq. (25) allows Eq. (26) to be derived.

$\left[P_{(-)}^{n+1}\right]=\left[A\right]\left[P_{(+)}^{n+1}\right]{\left[A\right]}^T+\left[\Gamma \right]\left[Q^n\right]{\left[\Gamma \right]}^T$E26

where $\left[Q^n\right]$ is the system error covariance matrix and is represented by $⟨\left\{q^{n+1}\right\}{\left\{q^{n+1}\right\}}^T⟩$.

5. Computational algorithm based on the Kalman filter FEM

If the matrices [A], [Q], and [R] are steady coefficient matrices, the computational algorithm can be divided into two parts: the computations of the Kalman gain matrix and the estimated value in time. The computational algorithm for the Kalman filter FEM is shown below.

5.1. Computation of Kalman gain matrix

1. Set input data: [A], $\left[P_{(+)}^0\right]$, $\left\{{\widehat{\varphi }}_{(+)}^0\right\}$, [Γ], [Q], [R], ε, (n = 0 ∼ imax)

2. Calculation of estimated error covariance matrix: $\left[P_{(-)}\right]=\left[A\right]\left[P_{(+)}\right]{\left[A\right]}^T+\left[\Gamma \right]\left[Q\right]{\left[\Gamma \right]}^T$

3. Calculation of Kalman gain matrix: $\left[K_1\right]=\left[P_{(-)}\right]{\left[H\right]}^T{\left(\left[H\right]\left[P_{(+)}\right]{\left[H\right]}^T+\left[R\right]\right)}^{-1}$

4. Calculation of predicted error covariance matrix: $\left[P_{(+)}\right]=\left[P_{(-)}\right]-\left[K_1\right]\left[H\right]\left[P_{(-)}\right]$

5. Check of convergence: if $\sqrt{tr\left(\left(\left[P_{(+)}^{k+1}\right]-\left[P_{(+)}^k\right]\right){\left(\left[P_{(+)}^{k+1}\right]-\left[P_{(+)}^k\right]\right)}^T\right)}\le \epsilon$ then go to Step 6, else go to Step 2.

5.2. Computation of estimated value in time

6. Calculation of estimated value: $\left\{{\widehat{\varphi }}_{(-)}^{n+1}\right\}=\left[A\right]\left\{{\widehat{\varphi }}_{(+)}^n\right\}$

7. Calculation of optimal estimated value: $\left\{{\widehat{\varphi }}_{(+)}^{n+1}\right\}=\left\{{\widehat{\varphi }}_{(-)}^{n+1}\right\}+\left[K_1\right]\left(\left\{z^{n+1}\right\}-\left[H\right]\left\{{\widehat{\varphi }}_{(-)}^{n+1}\right\}\right)$

6. Numerical example 1

Figures 2 and 3, respectively, show the computational model and the position of observation points in the computation using the Kalman filter FEM. The estimation of shallow water flow was carried out based on the Kalman filter FEM. Table 1 shows the computational conditions. In this study, we simulated the observed data from the shallow water flow analysis based on the FEM. The boundary condition for the water elevation is given by η(t) = sin (t) on the left-hand side of the open channel, and the artificial observed data is obtained by summation of the numerical result and observation noise, i.e., the white noise. In the computation of the Kalman filter FEM shown in this section, the inflow boundary condition is ignored, and the artificial observation data flow velocities u_x, u_y, and water elevation η, i.e., the computational result by the FEM, are given at the observation points.

Figure 2.

Computational model.

Figure 3.

Numerical examples.

Time increment $\Delta t$, s	0.001
Time steps	2000
Number of nodes	153
Number of elements	200
Gravitational acceleration g, m/s^2	9.8
Lumping parameter e	0.8
Initial of estimated error covariance $P_{(+)}^0$	1.0
Initial of estimated value ${\widehat{\phi }}_{(+)}^0$	0
Convergence determination constant ε	0.01

Table 1.

Computational conditions.

The computational results are shown below. Figures 4 and 5 show the convergence criterion expressed by the Frobenius norm. The equation for the Frobenius norm is shown in Step 5 of the flowchart in Section 4. It is seen that the convergence criterion monotonically decreases and converges. Figures 6 and 7 show the time history of water elevation at (x, y) = (5.0 m, 0.4 m) in Case A and Case B, respectively. It is found that the observation noise can be removed in both cases. Figures 8 and 9 show the distribution of water elevation on line y = 0.4 m in Case A and Case B, respectively. From Figures 8 and 9, it can be seen that if the observation point is not set at the upstream side in the channel, the distribution of the water elevation cannot be appropriately obtained. Conversely, from the result of Case A, it can be said that the distribution of water elevation can be correctly obtained if the position of the observation point is appropriately given in the computational model. Figures 10 and 11 show the distribution of flow velocity u and v and water elevation η at T = 2.0 s. It can be seen that water elevation on the upstream side in Case B is not appropriately obtained in comparison with the result of Case A. Therefore, from this numerical example, it can be said that it is necessary to set the observation points uniformly in the computational domain.

Figure 4.

Convergence criterion in Case A.

Figure 5.

Convergence criterion in Case B.

Figure 6.

Time history of water elevation at (x, y) = (5.0 m, 0.4 m) in Case A.

Figure 7.

Time history of water elevation at (x, y) = (5.0 m, 0.4 m) in Case B.

Figure 8.

Distribution of water elevation on line y = 0.4 m in Case A.

Figure 9.

Distribution of water elevation on line y = 0.4 m in Case B.

Figure 10.

Distribution of flow velocity u and v and water elevation η at T = 2.0 s in Case A.

Figure 11.

Distribution of flow velocity u and v and water elevation η at T = 2.0 s in Case B.

7. Numerical example 2

In this section, the numerical experiment is carried out with the inflow boundary conditions given in Case B, shown in the previous section. An image diagram of the numerical example is shown in Figure 12. The inflow boundary condition is given by η(t) = sin (t). The computational conditions in this section are the same as in Table 1. As the computational result, Figure 13 shows the convergence criterion expressed by the Frobenius norm. It can be seen that the convergence criterion decreases and converges monotonically. Figure 14 shows the time history of water elevation at (x, y) = (5.0 m, 0.4 m) in Case C. It is found that the observation noise can be removed. Figure 15 shows the distribution of water elevation on line y = 0.4 m in Case C. From Figure 15, it is found that the distribution of the water elevation can be appropriately obtained in comparison with the result shown in Figure 9. In addition, Figure 16 shows the distribution of flow velocity u and v and water elevation η at T = 2.0 s. It is seen that water elevation on the upstream side can be more accurately obtained than with the result for Case B. Therefore, from this numerical experiment, it is necessary to give the boundary conditions if the computation of the flow field estimation is carried out based on the Kalman filter FEM, and the observation point cannot be uniformly given.

Figure 12.

Numerical example.

Figure 13.

Convergence criterion in Case C.

Figure 14.

Time history of water elevation at (x, y) = (5.0 m, 0.4 m) in Case C.

Figure 15.

Distribution of water elevation on line y = 0.4 m in Case C.

Figure 16.

Distribution of flow velocity u and v and water elevation η at T = 2.0 s in Case C.

8. Comparison of estimation accuracy in Cases A, B, and C

Based on the estimation result at (x, y) = (5.0 m, 0.4 m) shown in Figures 6, 7, and 14, the numerical accuracy is investigated by comparison of data between the true and the estimated values in each case. The L₂ norm with respect to the difference between the true and the estimated values in each case is shown in Table 2.

Case	Number of observation points	Inflow boundary condition	L^2 norm
Case A	10	Not given	1.899
Case B	5 (only downstream side)	Not given	4.161
Case C	5 (only downstream side)	Given	2.463

Table 2.

Comparison of L² norm.

A comparison of the results in Cases A and B shows that it is necessary to set the observation points uniformly in the computational domain to obtain a highly accurate estimation result. In addition, when comparing the results in Case B and Case C, it is found that the estimation accuracy increases if the inflow boundary condition is taken into account.

9. Conclusions

In this study, flow field estimation analysis in an open channel was carried out based on the Kalman filter FEM. The shallow water equations were employed as the governing equations, and the Galerkin and the selective lumping methods were used to discretize the governing equation in space and time, respectively. In the numerical experiments, the observed value is created by adding white noise to the simulation result in the open channel model. The following conclusions were obtained from this study.

• The observation noise can be removed from the observed value by using Kalman filter FEM.

• When the observation points are set only on the downstream side of an open channel, the estimated water elevation on the upstream side is almost constant.

• If the inflow boundary condition is given in the case where the observation points are positioned only on the downstream side of an open channel, the distribution of the water elevation can be appropriately obtained, including the upstream side in the open channel.

Acknowledgments

This work was supported by Grants-in-Aid for Scientific Research (C) Grant Number 15K05786, and the contents of this work, i.e., formulation by FEM and the theory of the Kalman filter FEM, is based on seminars by Emeritus Professor Mutsuto Kawahara at Chuo University's Department of Civil Engineering. The computations were mainly carried out using the Fujitsu PRIMERGY CX400 computer facilities at Kyushu University's Research Institute for Information Technology. We wish to thank Emeritus Professor Mutsuto Kawahara and the staff at Kyushu University's Research Institute for Information Technology.