compute analysis
ua = ub + K @ (w-H@ub)
P = (np.eye(n) - K@H) @ B
return ua, P
```
Different forms for evaluating the Kalman gain and the covariance matrices are presented in literature. Some of them are favored for computational cost aspects, while others maintain desirable properties (e.g., symmetry and positive definiteness) for numerically stable implementation [27]. Since we are more interested in nonlinear dynamical models, we shall discuss extensions for standard Kalman filters to account for nonlinearity in the following sections.

✝ ✆

#### **7. Extended Kalman Filter**

Instead of dealing with linear stochastic dynamics, we look at the nonlinear case with general (nonlinear) observation operator written as

$$\mathbf{u}\_{l}(t\_{k+1}) = M(\mathbf{u}\_{l}(t\_{k}); \boldsymbol{\theta}) + \mathfrak{J}\_{p}(t\_{k+1}),\tag{44}$$

$$\mathbf{w}(t\_k) = h(\mathbf{u}\_t(t\_k)) + \mathfrak{J}\_m(t\_k). \tag{45}$$

The first challenge of applying Kalman filter for this system is the propagation of the background covariance matrix in the forecast step. The main clue behind the extended Kalman filter (EKF) to address this issue is to locally linearize *<sup>M</sup>*(**u**(*tk*)) by expanding it around the estimate **<sup>u</sup>***<sup>b</sup>*(*tk*) at *tk* using the first-order Taylor series as follows,

$$M(\mathbf{u}\_t(t\_k); \boldsymbol{\theta}) \approx M(\hat{\mathbf{u}}\_b(t\_k); \boldsymbol{\theta}) + \mathbf{D}\_M(\hat{\mathbf{u}}\_b(t\_k)) \mathbf{f}\_b(t\_k),\tag{46}$$

where **<sup>D</sup>***M*(**<sup>u</sup>***<sup>b</sup>*(*tk*)) is the Jacobian (also known as the tangent linear operator) of the forward model *<sup>M</sup>*(·; ·) evaluated at **<sup>u</sup>***<sup>b</sup>*(*tk*) and *<sup>ξ</sup>b*(*tk*)=(**u***t*(*tk*) <sup>−</sup> **<sup>u</sup>***<sup>b</sup>*(*tk*)) defining the error estimate at *tk*, with zero mean and covariance matrix of **B***<sup>k</sup>*. Thus, the difference between the background forecast and true state at *tk*<sup>+</sup><sup>1</sup> can be written as follows,

$$\begin{split} \mathfrak{J}\_{b}(t\_{k+1}) &= \mathbf{u}\_{l}(t\_{k+1}) - \mathbf{u}\_{b}(t\_{k+1}) \\ &= M(\mathbf{u}\_{l}(t\_{k}); \boldsymbol{\theta}) + \mathfrak{J}\_{p}(t\_{k+1}) - M(\hat{\mathbf{u}}\_{b}(t\_{k}); \boldsymbol{\theta}) \\ &\approx M(\hat{\mathbf{u}}\_{b}(t\_{k}); \boldsymbol{\theta}) + \mathbf{D}\_{M}(\hat{\mathbf{u}}\_{b}(t\_{k})) \hat{\mathfrak{J}}\_{b}(t\_{k}) + \mathfrak{J}\_{p}(t\_{k+1}) - M(\hat{\mathbf{u}}\_{b}(t\_{k}); \boldsymbol{\theta}) \\ &\approx \mathbf{D}\_{M}(\hat{\mathbf{u}}\_{b}(t\_{k})) \hat{\mathfrak{J}}\_{b}(t\_{k}) + \mathfrak{J}\_{p}(t\_{k+1}). \end{split}$$

Similar to the derivation in the linear case, with the assumption of uncorrelation between *ξb*(*tk*) and *ξ <sup>p</sup>*(*tk*<sup>+</sup>1), the background covariance matrix at *tk*<sup>+</sup><sup>1</sup> can be computed as follows,

$$\mathbf{B}\_{k+1} = \mathbf{D}\_M(\hat{\mathbf{u}}\_b(t\_k)) \hat{\mathbf{B}}\_k \mathbf{M}\_k^T + \mathbf{Q}\_{k+1}.\tag{47}$$

The next challenge regarding the analysis step is the computation of the Kalman gain in case of nonlinear observation operator. Again, *h*(**u***t*(*tk*<sup>+</sup>1)) is linearized using Talylor series expansion around **u***b*(*tk*<sup>+</sup>1) (i.e., the background forecast) as follows,

$$h\left(\mathbf{u}\_{l}(t\_{k+1})\right) \approx h(\mathbf{u}\_{b}(t\_{k+1})) + \mathbf{D}\_{h}(\mathbf{u}\_{b}(t\_{k+1})) \mathfrak{F}\_{b}(t\_{k+1}),\tag{48}$$

where **D***h*(**u***b*(*tk*<sup>+</sup>1)) is the Jacobian of the observation operator *h*, computed with the forecast **u***b*(*tk*<sup>+</sup>1). The Kalman gain is thus computed using this first-order approximation of *h* as follows,

$$\mathbf{K}\_{k+1} = \mathbf{B}\_{k+1} \mathbf{D}\_{h} (\mathbf{u}\_{b}(t\_{k+1}))^{T} \left( \mathbf{D}\_{h} (\mathbf{u}\_{b}(t\_{k+1})) \mathbf{B}\_{k+1} \mathbf{D}\_{h} (\mathbf{u}\_{b}(t\_{k+1}))^{T} + \mathbf{R}\_{k+1} \right)^{-1},\tag{49}$$

with an analysis estimate and analysis covariance matrix defined as

$$\mathbf{u}\_{\mathbf{d}}(t\_{k+1}) = \mathbf{u}\_{\mathbf{b}}(t\_{k+1}) + \mathbf{K}\_{k+1} \left( \mathbf{w}(t\_{k+1}) - h(\mathbf{u}\_{\mathbf{b}}(t\_{k+1})) \right), \tag{50}$$

$$\mathbf{P}\_{k+1} = \left(\mathbf{I}\_n - \mathbf{K}\_{k+1} \mathbf{D}\_h(\mathbf{u}\_k(t\_{k+1}))\right) \mathbf{B}\_{k+1}.\tag{51}$$

A summary of the EKF algorithm is described as follows,

$$\begin{aligned} \text{Inputs:} \qquad \hat{\mathbf{u}}\_b(t\_k), \mathbf{B}\_k, M(\cdot; \cdot), \mathbf{Q}\_{k+1}, \mathbf{w}(t\_{k+1}), \mathbf{R}\_{k+1}, h(\cdot) \\ \text{Forcesat:} \qquad \mathbf{u}\_b(t\_{k+1}) &= M(\hat{\mathbf{u}}\_b(t\_k); \theta) \\ \mathbf{B}\_{k+1} = \mathbf{D}\_M(\hat{\mathbf{u}}\_b(t\_k)) \hat{\mathbf{B}}\_k \mathbf{D}\_M(\hat{\mathbf{u}}\_b(t\_k))^T + \mathbf{Q}\_{k+1} \end{aligned}$$

$$\begin{array}{ll}\textbf{Kalman gain:} & \textbf{K}\_{k+1} = \textbf{B}\_{k+1}\textbf{D}\_{h}(\textbf{u}\_{b}(t\_{k+1}))^{T}\left(\textbf{D}\_{h}(\textbf{u}\_{b}(t\_{k+1}))\_{k+1}\textbf{B}\_{k+1}\textbf{D}\_{h}(\textbf{u}\_{b}(t\_{k+1}))^{T} + \textbf{R}\_{k+1}\right)^{-1} \\\\ \textbf{Analys:} & \textbf{u}\_{d}(t\_{k+1}) = \textbf{u}\_{b}(t\_{k+1}) + \textbf{K}\_{k+1}\left(\textbf{w}(t\_{k+1}) - h(\textbf{u}\_{b}(t\_{k+1}))\right) \\\\ & \textbf{P}\_{k+1} = (\textbf{I}\_{n} - \textbf{K}\_{k+1}\textbf{D}\_{h}(\textbf{u}\_{b}(t\_{k+1})))\textbf{B}\_{k+1}. \end{array}$$

and a Python implementation of the data assimilation step is presented in Listing 15.

**Listing 15.** Implementation of the (first-order) EKF with nonlinear dynamics and nonlinear observation operator. ✞ ☎

```
import numpy as np
def EKF(ub,w,ObsOp,JObsOp,R,B):