Lecture II-4: Filtering, Sequential Estimation, and the Ensemble Kalman Filter

1. Lecture II-4: Filtering, Sequential Estimation, and the Ensemble Kalman Filter • Lecture Outline: • A Typical Filtering Problem • Sequential Estimation • Formal Bayesian Solution to the Nonlinear Filtering Problem • The Linear Case -- Classical Kalman Filtering • Practical Solution of the Nonlinear problem – Ensemble Kalman Filtering • SGP97 Case Study

2. A Filtering Problem -- Soil Moisture Timeti-1 Timeti Timeti+1 Meteorological station State Eq. (Land surface model) Satellite footprint Meas Eq. (Radiative transfer model) Surface soil saturation contour Problem is to characterize time-varying soil saturation at nodes on a discrete grid. Estimates rely on scattered meteorological measurements and on microwave remote sensing observations (brightness temperature). Estimate must be compatible with relevant mass and energy conservation equations. State and meas. equations are nonlinear.

3. Limitations of the Variational Approach for Filtering Applications Variational methods provide an effective way to estimate the mode of f y|z(y|z) for interpolation and smoothing problems. However, they have the following limitations: • The adjoint equations for a complex model may be difficult to derive, especially if the model includes discontinuous functions. • Most variational methods rely on gradient-based search algorithms which may have trouble converging. • Variational data assimilation algorithms tend to be complex, fragile, and difficult to implement, particularly for large time-dependent problems. • Variational algorithms generally provide no information on the likely range of system states around the conditional mode. • Most variational methods are unable to account for time-varying model input errors when problem size is large. Considering these limitations, it seems reasonable examine other alternatives for the filtering problem. One option is sequential estimation.

4. Sequential Estimation For filtering problems, where the emphasis is on characterizing current rather than historical conditions, it is natural to divide the data assimilation process into discrete stages which span measurement times. To illustrate, suppose that the measurements available at all times through the current time ti are assembled in a large array Zi : Zi = [z1, z2, …, zi] = Set of all measurements through time ti Then the conditional PDF of y(ti), given all measurements through ti , is f yi|Zi[ y(ti)|Zi ] and the conditional PDF of the state at some later time (e.g.ti+1 ) given the same measurement informationis f y,i+1|Zi[ y(ti+1)|Zi ]. With these defintions we can define two distinct types of estimation steps: A propagation step -- Describes how the conditional PDF changes between two adjacent measurement times (e.g. between tiandti+1). An update step -- Describes how the conditional PDF changes at a given time (e.g. ti+1) when new measurements are incorporated. These steps are carried out sequentially, starting from the initial time between t0through the final time tI.

5. Propagating and Updating Conditional PDFs Meas. i Meas. 1 Meas. 2 zi z1 z2 Zi= [Zi-1 , zi] Z1= [z1] Z2= [Z1 , z2] t0 t1 t2 ti ti+1 f y1[ y1] f y2| z1[ y2|Z1 ] f yi| zi-1[ yi|Zi-1] f y,i+1| zi[ yi+1|Zi ] Propagation 1 to 2 Propagation i to i+1 Propagation 0 to 1 f y0 [ y0] f y1| z1[ y1|Z1 ] f yi| zi[ yi|Zi ] Update 1 Update i Algorithm initialized with unconditional (prior) PDF at t0 The propagation and update steps of the sequential estimation process are related as follows:

6. Formal Solution to the Nonlinear Bayesian Filtering Problem Propagation fromti to ti+1 : Update at ti+1: Bayes theorem provides the basis for deriving the PDF transformations required during the propgation and update steps of the sequential filtering procedure. For example, when the random input and the measurement error are additive and independent from time to time the transformation from f yi| Zi[ yi|Zi ] to f yi+1| Zi[ yi+1|Zi ] is: The transformation from f yi+1| Zi[ yi+1|Zi ] to f yi+1| Zi+1[ yi+1|Zi+1] is: These expressions constitute a formal sequential solution to the nonlinear filtering problem since everything on the right-hand side of each equation is known from the previous step. Although it is generally not possible to evaluate the required PDFs these Bayesian equations provide the basis for a number of practical approximate solutions to the filtering problem.

7. The Linear Bayesian Filtering Problem random For simplicity, assume a measurement is taken at every discretized time When the state and measurement equations are linear and the input and measurement error PDFs are multivariate normal the Bayesian propagation and update expressions yield a convenient sequential algorithm known as the Kalman filter. Suppose that the state and measurement equations have the following forms: If ui , y0 , and i are independent multivariate normal random variables then the conditional PDFs appearing in the Bayesian propagation and update equations are all multivariate normal. These PDFs are completely defined by their means and covariances, which can be derived from the general expression for the moments of a conditional multivariate normal PDF.

8. The Classical Kalman Filter forecast correction Propagation fromti to ti+1 : • Note that the expression for the updated conditional mean (the estimate ) consists of two parts: • A forecastE(yi+1|Zi) based on data collected through ti • A correction term which depends on the difference between the forecast and the new measurement zi taken at ti • The Kalman gainKi+1 determines the weight given to the correction term. The conditional mean and covariance given by the Bayesian filter for the linear multivariate case are: Kalman gain Update at ti+1:

9. Kalman Filter Example - 1 random Propagation: Kalman Gain Update: A simple example illustrates the forecast - correction tradeoff at the heart of the Kalman filter. Consider the following scalar state and measurement equations: In this special case the only input is a random initial condition with a zero mean and a specified covariance (actually a scalar variance) Cyy,0. The state remains fixed at this initial value but is observed with a set of noisy measurements. The filter equations are: Note that the filter tends to ignore the measurements (Ki+10 ) when the measurement error variance is large and it tends to ignore the forecast (Ki+11) when this variance is small.

10. Kalman Filter Example - 2 2.5 True value 2 Updated estimate 1.5 1 0.5 0 20 40 60 80 100 0 10 Kalman gain 1/t (Sample mean) -1 10 Time -2 40 60 80 100 20 0 10 In this example the covariance and Kalman gain both decrease over time and the updated estimate converges to the true value of the unknown initial condition. The Kalman gain is smaller and the convergence is slower when the measurement error variance is larger. Note that the Kalman gain approaches 1/t for large time. This implies that the estimate of the unknown initial condition is nearly equal to the sample mean of the measurements for the conditions specified here (variance of both initial condition and measurement error = 1.0).

11. Practical Solution of the Nonlinear Bayesian Filtering Problem Propagation fromti to ti+1 : Obtain L random replicates of the state yifrom the updateat ti. Generate L random replicates of the input vector ui during the interval [ti , ti+1). Solve the state equation for each replicate to construct a population of L random replicates for the state yi+1 at the end of the interval. Generate L random replicates of the measurement error i+1 at the new measurement time ti+1. Then use Bayes theorem(or a suitable approximation)to update each replicate of the state yi+1 . These replicates serve as the initial conditions for the next propagation step. Update at ti+1: The general Bayesian nonlinear filtering algorithm is not practical for large problems. The Kalman filter is intended for linear problems but can be applied to weakly nonlinear problems if the state and measurement equations are linearized about the most recent updated estimate. Unfortunately, the resulting extended Kalman filter tends to be unstable. In order to solve the large nonlinear problems encountered in data assimilation we need to take a different approach. One option is to use ensemble (or Monte Carlo) methods to approximate the probabilistic information conveyed by the conditional PDFs f yi| Zi[ yi|Zi ] and f yi+1| Zi[ yi+1|Zi ] . The basic concept is to generate populations of random inputs, measurement errors, and states which have the desired PDFs. This is easiest to illustrate within the two step sequential structure adopted earlier:

12. Ensemble filtering p(yi+1|Zi) Propagate forward in time Update with new measurement (zi+1) p(yi+1|Zi+1) p(yi|Zi) ti+1 Time Update with new measurement Propagate forward in time Time ti+1 ti ti Propagation of conditional probability density (formal Bayesian filtering): Evolution of random replicates in ensemble (Ensemble filtering): y l(ti+1| Zi+1) y l(ti| Zi) y l(ti+1| Zi) Ensemble filtering propagates only replicates (no PDFs). But how should update be performed? It is not practical to construct complete multivariate PDF and update with Bayes theorem.

13. The Ensemble Kalman Filter The updating problem simplifies greatly if we assumep[y(ti+1)| Zi+1] is Gaussian. Then update for replicate l is: Ki+1 = Kalman gain derived from propagated ensemble sample covariance Cov[y(ti+1)| Zi] and li is a random measurement error generated in the filter. After each replicate is updated it is propagated to next measurement time. There is no need to update the sample covariance. This is the ensemble Kalman filter (EnKF). • Potential Pitfalls • Appropriateness of the Kalman update for non-Gaussian density functions? • Need to construct, store, and manipulate large covariance matrices (as in classical Kalman filter)

14. Ensemble Kalman Filtering Example - 1 • Plot below shows ensemble Kalman filter trajectories for a simple two state time series model with perfect initial conditions, additive input error, and noisy measurements at every fifth time step: • 20 replicates are shown in green • Measurements are shown as magenta dots • The true state (used to generate the measurements) is shown in black • The mean of the replicates is shown in red • Note how ensemble spreads over time until update, when replicates converge toward measurement. 2.5 2 • 1.5 1 • • • 0.5 State 1 • • • 0 • -0.5 • • -1 -1.5 -2 0 5 10 15 20 25 30 35 40 45 50 Time

15. Ensemble Kalman Filtering Example - 2 It is difficult and not very useful to try to estimate the complete multivariate conditional PDF of y from the ensemble. However, it is possible to estimate means, median, modes, covariances, marginal PDFs and other useful properties of the state. Plot below shows histogram (approximate marginal PDF) for state 1 of the time series example. This histogram is derived from a 200 replicate ensemble at t = 18. 30 25 20 Number of replicates 15 10 5 0 -2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 State 1

16. Case Study Area Aircraft microwave measurements SGP97 Experiment - Soil Moisture Campaign

17. Mean land-atmosphere boundary fluxes Random input error Soil properties and land use Land surface model Radiative transfer model “True” soil, canopy moisture and temperature “True” radiobrightness Mean initial conditions Random meas. error Estimation error Random initial condition error Ensemble Kalman Filter Soil properties and land use, mean fluxes and initial conditions, error covariances Estimated radiobrightness and soil moisture Ensemble Kalman Filter Test – SGP97 Soil Moisture Problem Observing System Simulation Experiment (OSSE) “Measured” radiobrightness

18. Synthetic Experiment (OSSE) based on SGP97 Field Campaign Synthetic experiment uses real soil, landcover, and precipitation data from SGP97 (Oklahoma). Radiobrightness measurements are generated from our land surface and radiative transfer models, with space/time correlated model error (process noise) and measurement error added. SGP97 study area, showing principal inputs to data assimilation algorithm:

19. Area-average top node saturation estimation error Normalized error for open-loop prediction (no microwave meas.) = 1.0 Compare jumps in EnKF estimates at measurement time to variational benchmark (smoothing solution). EnKF error generally increases between measurements. Increasing ensemble size

20. Error vs. Ensemble Size top node saturation: mean-square-error (mse) difference at last update [-] mse(EnKF) - mse(smoother) -2 10 regression -3 10 1 2 3 4 10 10 10 10 ensemble size (Ne) Error decreases faster up to ~500 replicates but then levels off. Does this reflect impact of non-Gaussian density (good covariance estimate is not sufficient)?

21. Actual and Expected Error Ensemble Kalman filter consistently underestimates rms error, even when input statistics are specified perfectly. Non-Gaussian behavior? top node saturation rms error, Ne = 500 [-] 0.1 actual error (rms) expected forecast error (ensemble-derived) expected analysis error (ensemble-derived) 0.08 0.06 0.04 0.02 170 172 174 176 178 180 182 day of year

22. Ensemble Error Distribution Errors appear to be more Gaussian at intermediate moisture values and more skewed at high or low values. Uncertainty is small just after a storm, grows with drydown, and decreases again when soil is very dry. pixel 283: ensemble of top node saturation pixel 283: ensemble of top node saturation 0 0 0.7 0.7 5 5 forecast Ne = 500 forecast Ne = 500 10 10 0.6 0.6 analysis Ne = 500 analysis Ne = 500 15 15 0.5 0.5 0.4 0.4 top node saturation [-] precipitation [mm/h] top node saturation [-] precipitation [mm/h] 0.3 0.3 0.2 0.2 0.1 0.1 0 0 1 1 2 2 3 3 4 4 5 5 6 6 7 7 8 8 9 9 10 10 11 11 12 12 13 13 14 14 observation time observation time

23. Model Error Estimates pixel 283: model error in moisture flux upper b.c. [mm/d] 5 open loop true 4 EnKF Ne = 500 variational benchmark 3 2 1 0 -1 -2 170 172 174 176 178 180 182 day of year Ensemble Kalman filter provides discontinuous but generally reasonable estimates of model error but sample problem. Compare to smoothed error estimate from variational benchmark. Specified error statistics are perfect.

24. Summary • Variational methods can work well for interpolation and smoothing problems but have conceptual and computational deficiencies that limit their applicability to filtering problems • Sequential Bayesian filtering provides a formal solution to the nonlinear filtering problem but is not feasible to use for large problems. • Classical Kalman filtering provides a good solution to linear multivariate normal problems of moderate size but it is not suitable for nonlinear problems. • Ensemble Kalman filtering provides an efficient option for the solving large nonlinear filtering problems encountered in data assimilation applications. • Ensemble propagation characterizes distribution of system states (e.g. soil moisture) while making relatively few assumptions. Approach accommodates very general descriptions of model error. • Most ensemble filter updates are based on Gaussian assumption. Validity of this assumption is problem-dependent.