Sequential Monte Carlo Methods for Data Assimilation In Strongly Nonlinear Dynamics Zhiyu Wang B.Sc, Nanjing University, (China) 2005 Thesis Submitted In Partial Fulfillment Of The Requirements For The Degree Of Master of Natural Resources and Environmental Studies The University of Northern British Columbia April 2009 © Zhiyu Wang, 2009. 1*1 Library and Archives Canada Bibliotheque et Archives Canada Published Heritage Branch Direction du Patrimoine de I'edition 395 Wellington Street Ottawa ON K1A0N4 Canada 395, rue Wellington Ottawa ON K1A0N4 Canada Your file Votre reference ISBN: 978-0-494-48741-9 Our file Notre reference ISBN: 978-0-494-48741-9 NOTICE: The author has granted a nonexclusive license allowing Library and Archives Canada to reproduce, publish, archive, preserve, conserve, communicate to the public by telecommunication or on the Internet, loan, distribute and sell theses worldwide, for commercial or noncommercial purposes, in microform, paper, electronic and/or any other formats. AVIS: L'auteur a accorde une licence non exclusive permettant a la Bibliotheque et Archives Canada de reproduire, publier, archiver, sauvegarder, conserver, transmettre au public par telecommunication ou par Plntemet, prefer, distribuer et vendre des theses partout dans le monde, a des fins commerciales ou autres, sur support microforme, papier, electronique et/ou autres formats. The author retains copyright ownership and moral rights in this thesis. Neither the thesis nor substantial extracts from it may be printed or otherwise reproduced without the author's permission. L'auteur conserve la propriete du droit d'auteur et des droits moraux qui protege cette these. Ni la these ni des extraits substantiels de celle-ci ne doivent etre imprimes ou autrement reproduits sans son autorisation. In compliance with the Canadian Privacy Act some supporting forms may have been removed from this thesis. Conformement a la loi canadienne sur la protection de la vie privee, quelques formulaires secondaires ont ete enleves de cette these. While these forms may be included in the document page count, their removal does not represent any loss of content from the thesis. Bien que ces formulaires aient inclus dans la pagination, il n'y aura aucun contenu manquant. Canada Abstract General Bayesian estimation theory is reviewed in this study. The Bayesian estimation provides a general approach to handle nonlinear, non-Gaussian, as well as linear, Gaussian state estimation problems. The Sequential Monte Carlo (SMC) methods are presented to solve the nonlinear, non-Gaussian estimation problems. We compare the SMC methods with the Ensemble Kalman Filter (EnKF) method by performing data assimilation in the nonlinear, non-Gaussian dynamics. The Lorenz 1963 and 1996 models serve as test beds for examining the properties of these two estimation methods. Although EnKF computes only mean and variance based on the assumption of Gaussian dynamics, the SMC methods do not outperform EnKF in practical applications of the nonlinear non-Gaussian cases as we expect in theoretical insights. The reasons behind the experimental results that the SMC methods perform as well as EnKF in data assimilation and the future applications for high dimensional realistic atmospheric and oceanic models are discussed. To my beloved parents Acknowledgement s First and foremost, it is difficult to overstate my gratitude to my supervisor, Dr. Youmin Tang. Without his enthusiasm, his inspiration, and his great efforts, I would have been lost. I would also like to thank my graduate supervisory committee members Dr. Stephen J. Dery and Dr. Liang Chen for all the support they have given to me during the past years. I am indebted to my colleagues Dr. Xiaobing Zhou, Dr. Ziwang Deng, Mr. Yanjie Cheng, Mr. Jaison T. Ambadan, and Ms. Xiaoqin Yan for providing a stimulating and fun environment in which to learn and grow. In particular I would like to acknowledge the help of Dr. Youqin (Jean) Wang for her support at HPC Lab, UNBC. I am grateful to my professors, colleagues, and friends, Dr. Chris Johnson, Dr. Andrew Clifton, Dr. Jinjun Tong, and Mr. Joshua Laurin, from the University of Northern British Columbia and the City of Prince George for their continued support in these years. Finally, I am forever thankful to my parents for their understanding, endless patience and encouragement when it was most required. Contents 1 Introduction 1.1 Data Assimilation: A Brief Review 1 1 1.2 State-Space Form 1 1.3 1.4 1.5 Existing Methods General Case: Nonlinear, Non-Gaussian Research Objective 2 3 4 1.6 Outline of Thesis 4 2 Sequential Data Assimilation Methods 2.1 2.2 2.3 Nonlinear State Estimation Kalman Filter Framework 2.2.1 Extended Kalman Filter (EKF) 2.2.2 Ensemble Kalman Filter (EnKF) Sequential Monte Carlo (SMC) Methods 6 8 8 9 10 2.3.1 Perfect Monte Carlo Sampling 11 2.3.2 2.3.3 Sequential Importance Sampling (SIS) Sequential Monte Carlo Methods/Particle 12 12 2.3.3.1 Sequential Monte Carlo Algorithm 3 Assimilation Experiment I: Lorenz 1963 Model 3.1 3.2 3.3 3.4 6 Observation Interval 5t0bs = 0.50 Observation Interval 1) i) 6 2.1 N o n l i n e a r S t a t e E s t i m a t i o n p(xt\Yt-i) = / p(x t |a;t_i)p(a; t _i|r t _i)dx t _i (2.2) / p(j/ t |x t )p(x t |y t _i)dx t (2.3) where p(j/t|y t _i) = From (2.1), to obtain the conditional probability density function p(xt\Yt), need the observational noise probability density function p(yt\xt), forecast probability density function p(xt\Yt-i) we one step ahead which is the prior knowledge of the state variable, and marginal observational probability density functionp{y t \Y t -i). Since the probability density function p(yt\xt) can be calculated from the measurement model, the one step ahead forecast probability density function can be calculated from the dynamic model, and p(yt\Yt-i) p(xt\Yt-i) can be calculated according to (2.3). The target probability density function p(xt\Yt) can be esti- mated. After that, with the new observation coming in and the dynamic model forward evolution, this estimation algorithm becomes recursive (Doucet ei al, 2001; Schon, 2008). However, in general, there is no analytical solution to the nonlinear recursive estimation problem. This implies that we are forced to make approximations to approach this problem. The approximations suggested in the literature so far, can roughly be divided into two different classes, local approach and global approach (Schon, 2008). It is a matter of either approximating the nonlinear model and using the linear, Gaussian model estimator such as Extended Kalman Filter (EKF) or using the original nonlinear model and approximating the optimal solution such as Sequential Monte Carlo (SMC) methods. Despite the fact that there is a lot of different nonlinear estimators available, the local approximation approach is still the most commonly used nonlinear estimator when it comes to practical applications (Smith el aL, 1962; Schmidt, 1966; Evenseii & van Leeuwen, 1996; Houtekamer i=i J=I *w ^ ° > V i ( 2 - 17 ) where t denotes time, <$(•) is the Dirac delta function and qt^' denotes the weights associated with the particles x\ . The Dirac delta function S(-) can be defined as a function on the real line which is zero everywhere except at the origin, where it is infinite, and which is constrained to satisfy the identity +oo / 5[x)dx = 1 (2.19) -co The Dirac delta function <$(•) has the fundamental property that +oo / 2.3.1 f(x)5(x - a)dx = f(a) (2.20) -oo Perfect Monte Carlo Sampling In the perfect Monte Carlo sampling, all the random samples, also known as particles {x\ : i = 1 , . . . , M} are independent and identically distributed (i.i.d.) from the PDF p(xt\Yt), and every sample has equal weight, which is 1/M. The probability density function can be estimated by these samples according to (2.17). However, it is usually impossible to get i.i.d. samples from the PDF p(xt\Yt) at any time t. Nevertheless, the perfect Monte Carlo method shows the key idea in Sequential Monte Carlo (SMC) methods. 11 2.3 Sequential Monte Carlo (SMC) Methods 2.3.2 Sequential Importance Sampling (SIS) Unlike the perfect Monte Carlo sampling, all the i.i.d. samples are equally weighted, in Sequential Importance Sampling (SIS), all the i.i.d. samples are weighted according to importance weights qt^K In Sequential Importance Sampling, the importance weights qt^> contain the information on how probable it is that the corresponding sample was generated from the target PDF p(xt\Yt). Sequential Importance Sampling is a more general Monte Carlo method than the perfect Monte Carlo method. In the SIS implementation, as t increases, the importance weights become more and more skewed and tend to degenerate, which is called the sample impoverishment problem or the weight degeneracy problem. To avoid this weight degeneracy problem, one needs to introduce an additional selection step. In the selection step, the importance weights can be used as the acceptance probabilities, which allows us to generate approximately independent samples {x^}^ from the target density function to be estimated. This implies that the process of generating the samples from the target density function is limited to these samples. More specifically this is realized by resampling among the samples according to Pr(x{i) = x (i) ) = g(x (i) ), i = l,...,M (2.21) where q(x^) is the weights associated with the particles, and Pr(-) is the probability evaluation. Resampling step is first introduced in SIS by Rubin (1988) and the modified SIS is renamed after Sampling Importance Resampling (SIR). The SIR algorithm is closely related to the bootstrap procedure, introduced by Efron (1979). This relation is discussed in Smith & Cel.fa.nd (1992). 2.3.3 Sequential Monte Carlo Methods/Particle filter In SMC methods, predicted particles { i S J ^ are generated from the underlying dynamic model and the filtered particles from the previous time {xf^u^^. Conceptually, the predicted particles are obtained simply by passing the filtered particles through the system dynamics. Since the weight function reveals how 12 2.3 Sequential M o n t e Carlo ( S M C ) M e t h o d s probable the obtained measurement is given the present state, the more a certain particle explains the received measurement, the more probability that the particle was in fact drawn from the true density. Furthermore, a new set of particles {x^fjfij approximating p(xt \ Yt) is generated by resampling with replacement among the predicted p a r t i c l e s j : ! ; ^ ^ } ^ , belonging to the sampling density Pr{x% = 4 } _ J = q(x^), where q(x^') i = l,...,M (2.22) is the weights associated with the particles, and Pr(-) is the prob- ability evaluation. This procedure can be repeated over time, which forms the algorithm of SMC methods. This algorithm was first successfully implemented in practice by Gordon et al. (1993). Later it was independently rediscovered by Kitagawa (1996) and Isard & Blake (1998). Further references see Doucet et al. (2000), Kitagawa (1996), Liu & Chen (1998), Aruktrnpalam et al (2002). 2.3.3.1 Sequential M o n t e Carlo A l g o r i t h m This algorithm is used in Gordon et al (1.993) and Schon (2006). S t e p 1. Initialize the particles, {X-QVJJ^-, ~ pxo{xo) an d set t := 0 The particle filter is initialized by drawing samples from the prior density function p xo ( x 'o)S t e p 2. Measurement update: calculate the importance weights {q\ }fix ac- cording to qf^p^tlxf^), i = l,...,M (2.23) and normalize q\1' = q\• / Ylj=i Qt • In the measurement update, the new measurement is used to assign the probability, represented by the normalized importance weight q±', to each particle. This probability is calculated using the likelihood functionp(yt\%t\t-i)> wrncri describes how likely it was to obtain the measurement given the information available in the particle. 13 2.3 Sequential Monte Carlo (SMC) Methods Step 3. Calculate target probability density function p(xt\Yt), according to M M (0 p(xt\Yt) « J2 9* *(^ - 4 ° ) , E * W = X' i=l ®W ^ 0, Vi (2.24) i=l where t denotes time, (•) is the Dirac delta function and qt^' denotes the weights (i) associated with particles xt . The normalized importance weights and the corresponding particles constitute an approximation of the filtering probability density function p(xt\Yt). Step 4. Resampling: draw M particles, with replacement, according to Pr(x§ = x^_x) = q(x{j)), i = l,...,M (2.25) The resampling step will then return particles which are equally probable. Step 5. Time update: predict new particles according to 4li\t~P(xt+i\t\x§), i = l,...,M (2.26) The time update is just a matter of predicting new particles according to the underlying dynamic model and the filtered particles from the previous time {^t-iit-ili^i- Conceptually, the predicted particles are obtained simply by passing the filtered particles through the system dynamics. Step 6. Set t := t + 1 and iterate from step 2. Together with the new observations, these predicted particles form the starting point for another iteration of the assimilation algorithm. 14 Chapter 3 Assimilation Experiment I: Lorenz 1963 Model Both Sequential Monte Carlo (SMC) Methods (Gordon et aL, 1993) and Ensemble Kalman Filter (EnKF) (Evenseii, 1994) are sequential data assimilation methods and of stochastic nature, and both of them rely on Monte Carlo integration of the statistical behavior of the dynamic and measurement model system. Therefore, they have some similarities, and we can make some comparison to investigate the properties of these methods, especially in nonlinear and non-Gaussian dynamics. It is qtiite common that the atmospheric and oceanic dynamic systems are nonlinear and non-Gaussian. In this study, we choose the Lorenz model as a test bed (Lorenz. 1.983). It describes to some extent the nonlinear and chaotic nature of the atmosphere and ocean. The renowned Lorenz 1963 model was introduced by Edward Lorenz in 1963, who derived it from the simplified equations of convection rolls arising in the equations of the atmosphere. The Lorenz 1963 model consists of a system of three coupled and nonlinear ordinary differential equations (Lorenz, 1983), ~ = a(y-x) (3.1) dy_ dt px — y — xz (3.2) xy - Qz (3.3) dz_ lit 15 where, x(t), y(t), and z(t) are the dependent variables, and we have chosen the following commonly used values for the parameters in the equation: a = 10, p = 28, and (3 = 8/3. Our goal is to compare the properties and capabilities of Sequential Monte Carlo (SMC) methods and Ensemble Kalman Filter (EnKF) in strongly nonlinear non-Gaussian dynamics. What is the non-Gaussian dynamics? Let us define the Gaussian dynamics first. In the mathematical theory of probability, a Gaussian process {xt,t 6 T} is a stochastic process, of which the probability density function p is normally distributed. p(xi,x2,X3,xt,...,xt-a,xt-2,xt-i,xt) ~ N(ij,,a) (3.4) where p is the joint probability density function of the dynamic process, x is the random variable, and t refers to time. If the process is a not Gaussian process, it is a non-Gaussian process. The Gaussian process is distributed as the normal distribution, also called the Gaussian distribution, which is defined by two parameters, location and scale: the mean and variance (or standard deviation) respectively. The arithmetic mean is the average, often simply called the mean. The variance is one measure of statistical dispersion, averaging the squared distance of its possible values from the expected value (mean). Whereas the mean is a way to describe the location of a distribution, the variance is a way to capture its scale or degree of being spread out. The unit of variance is the square of the unit of the original variable. The positive square root of the variance, called the standard deviation, has the same units as the original variable. Mean \i is defined as: Standard deviation a is defined as: a "i i£pq-,i)2 (3.6) n . 8=1 16 where [i is the mean of the dynamic process. For a Gaussian process, the mean and the standard deviation fully characterize the probability density function of the process; however, for non-Gaussian process, the mean and the standard deviation only are insufficient, and higher order moments of the process are needed. Usually the coefficient of skewness and the coefficient of kurtosis are employed. In probability theory and statistics, the coefficient of skewness is a measure of the asymmetry of the probability distribution. A negative coefficient of skewness means the left tail is longer; the mass of the distribution is concentrated on the right of the figure. The distribution is said to be left-skewed. A positive coefficient of skewness means the right tail is longer; the mass of the distribution is concentrated on the left of the figure. The distribution is said to be rightskewed. The Coefficient of kurtosis is a measure of the peakedness of the probability distribution. The higher coefficient of kurtosis means more of the variance is due to infrequent extreme deviations, as opposed to the frequent, modestly-sized deviations. The coefficient of skewness 71 is defined as: 7l = nElLl(*i »f (3 7) where \i is the mean of the dynamic process. The coefficient of kurtosis 72 is defined as: i2 (3 8) -^m^^f? - where Li is the mean of the dynamic process. Through this thesis, we will use the mean [x, the standard deviation a, the coefficient of skewness 71, and the coefficient of kurtosis 72 as criteria to compare the assimilation results from EnKF method and SMC methods. Meanwhile, the first quartile, the second quartile (median), the third quartile, and the range of data are employed to check the assimilation results. 17 3.1 Observation Interval St0bs = 0.50 We design four different scenarios for the Lorenz 1963 model data assimilation. They axe assimilations with observation interval of 0.50, with observation interval of 0.25, with the initial error probability density function of Beta Distribution, and with the initial error probability density function of Gamma Distribution, see Table 3.1. Table 3.1: Experiment Design for Lorenz 1963 Assimilation Method Scenario 1 Scenario 2 Scenario 3 Scenario 4 3.1 SMC(250 particles) and EnKF(250 ensembles) Observation Interval 5t0bs = 0.50 Observation Interval 5t0bs = 0.25 Non-Gaussian Initial Error: Beta Non-Gaussian Initial Error: Gamma O b s e r v a t i o n I n t e r v a l St^,s = 0.50 The parameters of the Lorenz 1963 model in this case study are a — 10, p = 28, and P = 8/3. In this experiment we choose the same initial conditions as in Miller et ol (1994). The initial condition (x, y, z) is given by (1.508870, -1.531271, 25.46091), and the integration duration of the experiment is 50 dimensionless time units, with an integration time step of 0.01. The true value (reference resolution) is created by integrating the model with the above configurations. The distance between two nearest measurements is 5t0i>s = 0.50 and observations are made on the x, y and z coordinates. In this case study, the system initial error is Gaussian iV(0.0,2.0), and the observational error is also Gaussian iV(0.0, 2.0). The observations are simulated by adding normally distributed noise with zero mean and variance equal to 2.0 to the true value (reference solution). Initial conditions are also simulated by adding normally distributed noise with zero mean and variance equal to 2.0 to the true value (reference solution). This system of equations is integrated by Numerical Algorithms Group (NAG) Numerical Libraries with the fourth-order Runge-Kutta method. The assimilation 18 3.1 Observation Interval 5t0t,s = 0.50 experiments are run on an SGI Altix 3000 (64 Intel Itanium - 2 1500 MHz CPUs) global shared memory supercomputer. The filter performance will be evaluated by three factors: 1) root mean square error (RMSE); 2) CPU computation time; 3) statistics of the probability density function (PDF), which is estimated from the true resolution and assimilated estimates. The root mean square error (RMSE) is calculated between the reference solution and the filtering estimate (analysis) averaged over the whole assimilation period. We performed both the SMC methods and the EnKF method data assimilation in the Lorenz 1963 model, with different numbers of SMC particles and EnKF ensemble members. The number of SMC particles is 250. The number of EnKF ensemble members is also 250. The assimilation for x, y, and z is performed. However, the assimilation method is independent of state variables, thus the assimilation results for three variables x, y, and z are quite similar. Therefore, only the assimilation result for x is showed in this thesis. Table 3.2: Computation time and RMSE for Lorenz 1963 (Case: 8t0i,s = 0.50) Assimilation Method Time (In Seconds) RMSE (X) RMSE (Y) RMSE (Z) SMC (250 particles) 4.830 1.8520 2.9850 2.7383 19 EnKF (250 ensembles) 11.246 2.0208 3.2572 2.9262 Statistics Mean Standard deviation Coefficient of skewness Coefficient of kurtosis First Quartile Second Quartile Third Quartile Range True PDF SMC (250 particles) EnKF(250 ensembles) True - SMC True - EnKF 0.6415 0.7211 0.8688 -0.0795 -0.2272 7.8752 8.0119 -0.0645 7.9397 -0.1366 -0.1544 -0.2073 0.0222 -0.1766 0.0528 0.0252 -0.6077 -0.6330 -0.7009 0.0931 -4.3494 -0.0363 0.0809 -4.3857 -4.4666 -0.4978 1.0567 0.8517 1.5545 0.2050 -0.6212 -0.3047 5.9689 6.2736 6.5900 36.0129 -0.4030 -0.5929 35.4200 35.8230 Table 3.3: Statistics of PDF of x-component of Lorenz 1963 (Case: 5t0bs = 0.50) ! o 3 i-! gCD i—i 0 3 I CD i-i Ol O cr CO 3.1 Observation Interval Stobs = 0.50 State Estimate Error Variance 1 0.008 • 0.006 • 0.004 0.002 1j - Jj i L1u ^ I u* Jl 5 • 10 15 J. 11 if 1 „ ft. f L m MiUbi J i k u4«,y 20 25 30 35 :Aj 40 ,.*J\*M iUvk 45 50 time t State Estimate Error Variance 1 0.008 - 0.006 - r1 jrc , - 0.004 0.002 - .AMJI I I J 1 J I ,JUl JL* 1Hh LiIII . U.JidJl ...lii JIII J J JL ILLJJ 10 15 20 25 time t 30 35 40 45 J 50 Figure 3.1: State estimate and error variance of ^-component of Lorenz 1963 model for EnKF and SMC methods with filter size of 250, 5i0t,s = 0.50 21 3.1 Observation Interval 5tobs = 0.50 Density 0.01 10 15 20 25 10 15 20 25 Density -25 -20 -15 -10 -5 0 5 Figure 3.2: Probability density function of as-component of Lorenz 1963 model for EriKF and SMC methods with filter size of 250, 5tobs = 0.50 22 3.1 Observation Interval Stobs = 0.50 Figure 3.1 shows the variable x estimate and the error variance estimate with time step from both SMC methods and EnKF method with 250 particles and ensemble members. The error variance which defined, as in Evetisen (.1997) is scaled by N where N is a scalar quantity of time steps of the total assimilation period Error Variance = _L(Xf stemate - X f i e ) 2 (3.9) Both the SMC methods and the E n K F method do reasonably good jobs in tracking the phase transitions and also in reproducing the correct amplitudes of the reference solution. There are some locations where the filter estimates start to diverge from the reference solution. To compare the assimilation results, we divide the total assimilation period into two, the first half and the second half, so that we want to examine whether the assimilation results become better with more observations coming into the data assimilation system. Meanwhile, since the error variance is already rescaled, the error variance in all figures is just used to compare its relatively range within two methods, E n K F and SMC methods. To compare the RAISE variation with time, we choose 0.006 as a standard. For the E n K F method, at time 3, 11, 17, 23, 27, 30, 35, 37, 39, 42, 44, and 49, the error variance is greater t h a n 0.006, which means the filter estimates deviate from the true solution. Among these locations, 8 out of 12 are in the second half of the assimilation period. For SMC methods, at time 5, 9, 14, 17, 23, 30, 35, 42, 43, and 47, the error variance is greater than 0.006. Among these locations, 5 out of 10 are in the second half of the assimilation period. Despite these divergences, both methods recover quickly and track the reference solution again. Table 3.2 indicates the CPU computation time and the RMSE for both meth- ods in this case. The CPU computation time is 4.830 s for SMC methods, while 11.246 s for EnKF method. The E n K F method takes almost twice longer than the SMC methods. The RMSE is 1.8520, 2.9850, and 2.7383 for x, y, and z for SMC methods, while it is 2.0208, 3.2572, and 2.9262 for x, y, and z for E n K F method. SMC methods are slightly better than the EnKF method in this case. 23 3.2 Observation Interval 5tobs = 0.25 Figure 3.2 shows the probability density function of x component of the Lorenz 1963 dynamic system for both EnKF and SMC methods. The probability density function is calculated by the kernel density estimation method (Parzen, 1982 and Silverman, 1986). In Matlab, the kernel density estimation is implemented through the ksdensity function. In this thesis, all probability density functions are estimated by Matlab ksdensity function. From Figure 3.2, we can see clearly that the probability density functions are non-Gaussian. Both methods can assimilate it quite well, but, they both have some difficulties to reach the exact peaks of the probability density function. Table 3.3 shows the statistics of the probability density function of x component of the Lorenz 1963 model. The mean, the standard deviation, the coefficient of skewness, and the coefficient of kurtosis for the SMC methods are closer to those of the true state than those of the EnKF method. The quartiles and range from the SMC methods are closer to those of the true resolution than those of the EnKF method. Therefore, the SMC methods estimate the probability density function slightly better than the EnKF in this case. 3.2 Observation Interval 8t0bs = 0.25 In this second experiment, the experimental setup is the same as the previous one, except that the observation interval between two measurements decreases from St0bs = 0.50 to St0bs ~ 0.25 for this case. That means we have more observations in the assimilation process in case 2 than in case 1. Table 3.4: Computation time and RMSE for Lorenz 1963 (Case: 5t0bs = 0.25) Assimilation Method Time (In Seconds) RMSE (X) RMSE (Y) RMSE (Z) SMC (250 particles) 5.309 1.4003 2.1914 1.9663 24 EnKF(250 ensembles) 18.169 1.0156 1.6157 1.5755 Statistics Mean Standard deviation Coefficient of skewness Coefficient of kurtosis First Quartile Second Quartile Third Quartile Range True PDF SMC(250 particles) EnKF(250 ensembles) True - SMC True - EnKF 0.6415 0.7080 0.7562 -0.0664 -0.1146 7.8370 -0.1137 7.9889 7.8752 0.0382 -0.1544 -0.0340 -0.1200 0.0240 -0.1785 -0.6077 0.0775 -0.6853 0.0229 -0.6307 0.0054 -4.3857 -4.3910 0.4691 -4.8548 -0.1929 1.0567 1.3272 1.2496 -0.2705 -0.0937 5.9689 6.0626 -0.1371 6.1060 0.9314 34.4886 35.4200 35.3180 0.1020 Table 3.5: Statistics of PDF of x-component of Lorenz 1963 (Case: St0bs — 0.25) en o a O 13 I a CO cr O 3.2 Observation Interval St0hs = 0.25 State Estimate Error Variance 0.008 0.006 0.004 0.002 0 LJii 5 10 ikhkiu 15 il 20 L iKk.l.Ai I J-J-Xai 25 time t 30 35 d.-JlLL ,,., i h 40 45 50 State Estimate Error Variance , -J SMC | " 0.008 0.006 - - 1 - 0.004 ! 0.002 'I l WfhJ k*iU •1/1 1 . 5 10 1 30 L .J35,1 i>b ,Ak.kM Urn20ijj L1U1I 25 40 45 MSJ. 15 1 lift.] *,S .h. 50 time t Figure 3.3: State estimate and error variance of x-component of Lorenz 1963 model for EnKF and SMC methods, Filter size = 250, 5tohs = 0.25 26 3.2 Observation Interval 5t0i,s — 0.25 Density 0.06 1 ,„„ 1 • ^ 0.05 // \ - 0.04 0.03 - v^ / s~~~y 0.02 - if 0.01 -25 - \ -20 -15 -10 -5 0 5 10 15 20 25 Density 0.06 1 1 1 1 --' - .T. «„ c. 1| 1 /y 0.05 / I 0.04 1/ 0.03 i I / \ \ \ " \ " // // N^>~\ 0.02 / \ If 0.01 -25 - \ -20 -15 -10 -5 0 5 10 15 20 25 Figure 3.4: Probability density function of as-component of Lorenz 1963 model for EnKF and SMC methods, Filter size = 250, Stobs = 0.25 27 3.2 Observation Interval 5t0bs — 0.25 Figure 3.3 shows the variable x estimate and the error variance estimate with time step from both the SMC methods and the EnKF method with the filter size of 250. In tracking the phase transitions, there are some locations where the filter estimates diverge from the reference solution. For EnKF method, at time 5, 11, and 30, the error variance is greater than 0.006, which means filter estimates deviate from true solution. Among these locations, 1 out of 3 is in the second half of the assimilation period. For SMC methods, at time 4, 5, 11, 20, 24, and 34, the error variance is greater than 0.006. Among these locations, 1 out of 6 is in the second half of the assimilation period. With more observations, the EnKF method outperforms the SMC methods. For both methods, the transition in the second half becomes smoother than that in the first half. Despite these divergences, both methods recover quickly and track the reference solution again. From the error variance with time, we can see the error variance decreases with time in this case. Table 3.4 indicates the CPU computation time and the RMSE for both methods in the case. The CPU computation time is 5.309 s for SMC methods, while it is 18.169 s for EnKF method. The EnKF method takes almost 3 times longer than the SMC methods. The RMSE is 1.4003, 2.1914, and 1.9663 for x, y, and z for SMC methods, while it is 1.0156, 1.6157, and 1.5755 for x, y, and z for EnKF method. The EnKF method is significantly better than the SMC methods in this case with more observations available. Figure 3.4 shows the probability density function of x component of the Lorenz 1963 dynamic system for both the EnKF method and the SMC methods. From Figure 3.4, we can see clearly that the probability density function is non-Gaussian. Both methods can assimilate this nonlinear dynamic process quite well; however, the EnKF method almost reaches the exact peak of the probability density function, which is better than the SMC methods. Table 3.5 shows the statistics of the probability density function of x component of Lorenz 1963 model. The mean for the SMC methods are closer to the true resolution than that of the EnKF method, while the standard deviation, the coefficient of skewness, and the coefficient of kurtosis for the EnKF method are closer to the true state than the SMC methods, as well as the quartiles and range of the data. 28 3.3 Non-Gaussian Initial Error: Beta In Table 3.2 and Table 3.4, EnKF takes more than twice the time than the SMC methods do. That means both methods can achieve reasonably good results, but the SMC methods is more efficient than EnKF in this case. The EnKF algorithm used in this thesis is explained in Evensen (2003). We need to perform an analysis algorithm to each individual member, which is why EnKF takes much more time than SMC methods. In addition to the reason above, the EnKF analysis algorithm requires the calculation of the inverse of matrix, which is quite time consuming. One way to reduce the computational time for the EnKF is to reduce the ensemble size. For this, one possible option is to replace random perturbation in Kalman Filter by a deterministic perturbation, which turns out to be Unscented (Sigma-Point) Kalman Filter (Julier & Uhlmami, 1998). 3.3 Non-Gaussian Initial Error: Beta For the third case, we keep the experimental setup the same as the first one; expect that the system initial error is non-Gaussian, Beta Distribution. In Case Beta, Beta (2.0, 5.0) is used as the initial probability density function for model integration. In the first two cases, the model starts with a Gaussian error, and the probability density function may become non-Gaussian after the model iteration starts. In this case, in the beginning, the model starts with a non-Gaussian probability density function, and it will remain non-Gaussian after the model iteration. The EnKF method always uses Gaussian marginal probability density function to represent non-Gaussian marginal probability density function during the assimilation process; theoretically it is not sufficient, because only lower-order moments (mean, variance) are considered. While SMC methods directly estimate non-Gaussian marginal density function, theoretically it is much better than EnKF. Figure 3.5 shows the variable x estimate and the error variance estimate with time from both SMC methods and EnKF method with 250 particles and ensemble members. In tracking the phase transitions there are some locations where the filter estimates diverge from the reference solution. In spite of these divergences, both methods recover quickly and track the reference solution again. For the 29 3.3 Non-Gaussian Initial Error: Beta EnKF method, at time 2, 5, 6, 11, 16, 23, 30, 37, 43, 44 and 45, the error variance is greater than 0.006, which shows filter estimates deviate from true solution. Among these locations, 5 out of 11 are in the second half of the assimilation period. For the SMC methods, at time 2, 3, 6, 9, 11, 16, 17, 23, 27, 30, 37, and 44, the error variance is greater than 0.006. Among these locations, 4 out of 12 are in the second half of the assimilation period, which indicates that the second half assimilation transition for the SMC method is smoother than that for the EnKF method. Table 3.6 indicates the CPU computation time and the RMSE for both methods in the case. The CPU computation time is 4.877s for SMC methods, while it is 11.297s for EnKF method. The RMSE is 2.1640, 3.3007, and 3.8650 for x, y, and z for SMC methods, while it is 2.3394, 3.5670, and 2.9842 for x, y, and z for EnKF method. For variables x and y, the SMC methods is better, while for variables z, the EnKF method is better. Figure 3.6 shows the probability density function of the x component of the Lorenz 1963 dynamic system for both EnKF and SMC methods. From Figure 3.6, we can see clearly that the probability density function is non-Gaussian, which has one major peak and two weak peaks. Both methods can assimilate it reasonably well, but both of them have trouble to track the exact peaks of the probability density function. Table 3.7 shows the statistics of the probability density function of x of Lorenz 1963 model. The EnKF assimilated result is closer to the true resolution than SMC methods in the mean, the coefficient of skewness, the median, the third quartile and the range; while the SMC methods are better in the standard deviation, the coefficient of kurtosis, and the first quartile than the EnKF method. 30 3.3 Non-Gaussian Initial Error: Beta State Estimate Error Variance 0.01 ' ' ' I EnKF | 0.008 - " 0.006 - " 0.002 - | " 0.004 / 0 d I ..J An 5 f L m,, M«JN J JLJ km.k i Mi J l', «.AJ U 10 15 20 25 time t 30 35 30 35 40 uAJA.~. . , ^ 1 45 50 State Estimate Error Variance 0.008 0.006 0.004 0.002 UJ 5 | JAAL 10 15 iM»J\ 20 25 time t wML,40 Jh .W ki 45 50 Figure 3.5: State estimate and error variance of rc-component of Lorenz 1963 model for EnKF and SMC methods, Filter size = 250, Initial Beta distribution 31 3.3 Non-Gaussian Initial Error: B e t a Table 3.6: Computation time and RMSE for Lorenz 1963 (Case: Beta) Assimilation Method SMC(250 particles) Time (In Seconds) RMSE (X) RMSE (Y) RMSE (Z) 4.877 2.1640 3.3007 3.8650 32 EnKF(250 ensembles) 11.297 2.3394 3.5670 2.9842 Statistics Mean Standard deviation Coefficient of skewness Coefficient of kurtosis First Quartile Second Quartile Third Quartile Range True PDF SMC(250 particles) EnKF(250 ensembles) True - SMC True - EnKF 0.6415 0.9465 0.7676 -0.3049 -0.1260 7.8752 8.0014 -0.1120 -0.1262 7.9873 -0.1544 -0.1659 0.0115 0.0747 -0.2291 -0.6077 -0.6902 0.0824 0.1348 -0.7426 0.0270 0.3221 -4.3857 -4.4127 -4.7078 1.3963 -0.6278 1.0567 1.6845 -0.3396 -0.7714 5.9689 6.7403 6.5498 -0.5809 35.2279 0.6623 0.1921 35.4200 34.7577 Table 3.7: Statistics of PDF of ^-component of Lorenz 1963 (Case: Beta) 3.3 Non-Gaussian Initial Error: Beta Density 0.06 ' ' ' 0.05 1 1 EnKF 1 - //"\ L-j 0.04 1 x 0.03 1 1 - - \\ \\ ~S-"\ 1 1 0.02 - - / It II II II 0.01 ' n -25 -20 -15 -10 -5 0 5 10 15 20 25 5 10 15 20 25 Density -25 -20 -15 -10 -5 0 Figure 3.6: Probability density function of x-component of Lorenz 1963 model for EnKF and SMC methods, Filter size — 250, Initial Beta distribution 34 3.4 Non-Gaussian Initial Error: Gamma 3.4 Non-Gaussian Initial Error: G a m m a For the fourth case, we keep the experimental setup the same as the third one; expect that the system initial error is non-Gaussian, Gamma Distribution. In Case Gamma, Gamma (2.0, 2.0) is used to integrate the model forward as the initial error probability density function. Figure 3.7 shows the variable x state estimate and the error variance estimate variation with time forward for both the SMC methods with 250 particles and the EnKF method with 250 ensemble members. In tracking the phase transitions there are some locations where the filter estimates diverge from the reference solution. For example, for the EnKF method, at time 4, 5, 18, 24, 30, 36, 37, 44, and 49, the error variance is greater than 0.006, which indicates filter estimates deviate from true solution. Among these locations, 5 out of 9 are in the second half of the assimilation period, or rather, the estimate deviation occurs all through the assimilation process. For the SMC methods, at time 2, 4, 5, 9, 17, 20, 21, 24, 37, 39, 44, and 49, the error variance is greater than 0.006. Among these locations, 4 out of 12 are in the second half of the assimilation period, which shows the second half assimilation is better than the first half. Despite these divergences, both methods recover quickly and track the reference solution again. In general, both the SMC methods and the EnKF method do the case experiment well. Table 3.8 indicates the CPU computation time and the RMSE for both methods in the case. The CPU computation time is 4.925 s for SMC methods, while it is 11.295 s for EnKF method. The EnKF method takes almost as three times long as the SMC methods. The RMSE is 2.5674, 4.0898, and 3.7482 for x, y, and z for SMC methods, while it is 2.3808, 3.5195, and 3.4267 for x, y, and z for EnKF method. The EnKF method is slightly better than that of the SMC methods in this case with Gamma Initial Distribution. Figure 3.8 shows the probability density function of the x component of the Lorenz 1963 dynamic system for both EnKF and SMC methods. From Figure 3.8, we can see clearly that the probability density function is non-Gaussian, since it has one major peak and another two small peaks. Both methods can assimilate it reasonably well, but neither of them can reach the exact peaks of the probability 35 3.4 Non-Gaussian Initial Error: Gamma density function. Based on the 4 cases above, the process probability density functions are quite similar to one another in the 4 experiments, which indicate that the process PDF is mainly determined by dynamics itself. Still we could try to find the optimal initial probability density function for the specific data assimilation sytem. Table 3.9 shows the statistics of the probability density function of x of the Lorenz 1963 model. Both the mean for the SMC methods and the EnKF method are slightly larger than the true resolution; the standard deviation for the SMC methods is smaller than that of the true resolution, while the standard deviation for the EnKF method is slightly larger than that of the true resolution, the coefficients of skewness are -0.1544, -0.1819, and -0.1897; the difference of coefficients of kurtosis between true resolution and assimilated estimate for the EnKF method is 10 times larger than that of the SMC methods, both absolute values are quite small though. Moreover, all the quartiles in the EnKF method are better than that of the SMC methods, except the range in the SMC methods are better. The SMC methods have the theoretical advantage, why do the experimental results show similar estimates? Why do not the SMC methods outperform EnKF methods? p(xi,x2, • • • ,xn) is the probability density function of the dynamic process, which is non-Gaussian. The marginal probability density function p(xn) is estimated through the data assimilation process is also non-Gaussian. In EnKF, we assume all the marginal probability density function is Gaussian, which is not true for non-Gaussian dynamics, no matter whether it is linear or nonlinear dynamics, and the mean and covariance of marginal probability density is used to fully characterize the dynamics. However, in SMC methods, the marginal probability density function is estimated directly from sample particles. The RMSE is calculated from true resolution and assimilation resolution, it is only mean value or rather the first order of the moment of marginal probability density function. EnKF employs Gaussian marginal probability density function to represent non-Gaussian marginal probability density functions, which is not sufficient theoretically. However, it may be sufficient to represent the mean value of marginal probability density function. That is why RMSE for both EnKF and SMC are quite similar. Since we do not have a true non-Gaussian marginal 36 3.4 Non-Gaussian Initial Error: Gamma probability density function as a reference, it is difficult to verify whether the marginal density function in SMC methods can fully represent true dynamic marginal probability density function or not. Table 3.8: Computation time and RMSE for Lorenz 1963 (Case: Gamma) Assimilation Method Time (In Seconds) RMSE (X) RMSE (Y) RMSE (Z) SMC(250 particles) 4.925 2.5674 4.0898 3.7482 37 EnKF(250 ensembles) 11.295 2.3808 3.5195 3.4267 Statistics Mean Standard deviation Coefficient of skewness Coefficient of kurtosis First Quartile Second Quartile Third Quartile Range True PDF SMC(250 particles) EnKF(250 ensembles) True - SMC True - EnKF 0.6415 0.7903 0.7485 -0.1487 -0.1069 0.0995 -0.1047 7.8752 7.7757 7.9799 -0.1544 -0.1897 0.0275 0.0353 -0.1819 -0.6077 -0.6958 -0.0097 0.0880 -0.5980 -0.1002 -4.5641 -4.3857 0.1785 -4.2855 1.0532 -0.2378 0.0035 1.0567 1.2945 6.2274 6.4560 -0.4871 -0.2585 5.9689 -0.2634 35.4200 35.6834 35.0820 0.3380 Table 3.9: Statistics of PDF of re-component of Lorenz 1963 (Case: Gamma) 3.4 Non-Gaussian Initial Error: Gamma State Estimate x 0 Error Variance 0.008 0.006 0.004 0.002 1I fJMB^jii, 10 15 Mx)J 20 Jk± 25 time t M.J 30 35 ii 40 At,JIM 45 50 State Estimate Error Variance 0.01 r 0.008 0.006 0.004 0.002 0 J l i l i y j L , JLll Jji .Jlli i ILdL J.K Jffi i. .1,,1 Mi, ihi 10 15 20 25 time t 30 35 40 45 50 Figure 3.7: State estimate and error variance of x-component of Lorenz 1963 model for EnKF and SMC methods. Filter size = 250, Initial Gamma distribution 39 3.4 Non-Gaussian Initial Error: Gamma Density 1 1 1 ' j EnKF 1 - \ \\ i w. — \ / " - /^2^ ~^z^y - -25 0.06 - -20 -15 -10 1 -5 0 5 1 Density 1 1 f~\ f 10 15 20 ' I s•ol T 1 t -1 \ l/\ 0.05 25 \ /I ' 0.04 - - \ // / 0.03 - - V. \ - 0.02 1/ >' 0.01 -25 -20 -15 -10 -5 0 5 10 15 V, 20 25 Figure 3.8: Probability density function of ^-component of Lorenz 1963 model for EnKF and SMC methods. Filter size = 250. Initial Gamma distribution 40 Chapter 4 Assimilation Experiment II: Lorenz 1996 Model The second experimental design employs the Lorenz 1996 model as the test bed. The Lorenz 1996 model (Lorenz, 1996) represents an atmospheric variable X at J equally spaced points around a circle of the constant latitude. The j t h component is propagated forward in time following the differential equation f]X - ^ = (Xj+1 - XJ-JXJ-L - Xj + F (4.1) where j = 0 , . . . , J — 1 represents the spatial coordinates (longitude). F is a constant external forcing term, which indicates the dynamics is weakly chaotic when F — 5 or F = 6, it is highly chaotic when F — 8, and it is fully turbulent when F = 16. Note that this model is not a simplification of any atmospheric system, however, it is designed to satisfy three basic properties: it has linear dissipation (the — Xj term) that decreases the totally energy , an external forcing term F that can increase or decrease the total energy and a quadratic advection term that conserves the total energy just like many atmospheric models. In its configuration, J — 40 variables and boundary conditions are cyclic, i.e. X-\ = Xj-i, Xo = X^ and Xj+i = X\. which means the distance between two adjacent grid points roughly represents the midlatitude Rossby radius (about 800 km), assuming the circumference of the midlatitude belt is about 30000 km (Majda k Harlim, 2008). 41 In this experimental design, we define three different forcing term F scenarios. The first category is F — 5, which indicates that the model is weakly chaotic, the second category is F — 8, which indicates that the model is highly chaotic, and the last category is F — 16, which indicates that the model is fully turbulent, see Table 4.1. This dynamic model is integrated by the Numerical Algorithms Group (NAG) Numerical Libraries with the fourth-order Runge-Kutta method, and the integration time step of 0.05, corresponding to 6 hours in the realistic atmospheric physics. The initial condition is given after a spin up integration for 10 years. The duration of the experiment setup is 40 dimensionless time units. The observation interval between two measurements is 5t0f,s — 0.50 and observations are available for all 40 variables. In this case study, the system initial error is Gaussian A(0.0,2.0), and observational error is also Gaussian N(0.0,2.0). The observations are simulated by adding normally distributed noise with zero mean and variance equal to 2.0 to the reference solution. Initial conditions are simulated by adding normally distributed noise with zero mean and variance equal to 2.0 to reference solution. The filter performance will also be evaluated by the root mean square error (RMSE) between the true value (reference solution) and the filtering estimate averaged over the whole assimilation period, the CPU computational time, and the statistics of the probability density functions. The assimilation experiments ran on an SGI Altix 3000 (64 Intel Itanium - 2 1500 MHz CPUs) global shared memory supercomputer. We performed both the SMC methods and the EnKF method data assimilation in the Lorenz 1996 model, with different numbers of SMC particles and EnKF ensemble members equal 250. For the Lorenz 1996 model, 40 variables are functionally equal. Also the assimilation method is independent of model state variables. Therefore, only the assimilation result for state variable XI, A20, and X30 are shown in this thesis. 42 4.1 Weakly Chaotic F = 5 Table 4.1: Experiment Design for Lorenz 1996 Assimilation Method Scenario 1 Scenario 2 Scenario 3 4.1 SMC(250 particles) and EnKF(250 ensembles) Weakly Chaotic F = 5 Highly Chaotic F = 8 Fully Turbulent F = 16 Weakly Chaotic F = 5 In Fig 4.1, both the SMC methods and the EnKF method perform reasonably well in tracking the phase transitions and also in reproducing the correct amplitudes of the reference solution. We divide the whole assimilation period into two. Both the SMC methods and the EnKF method take almost half of the whole period to start to track the reference solution accurately. When the dynamics are weakly chaotic, both methods can assimilate the dynamic process well and quickly, since the second half is obviously better than the first half. However, there are some locations where the filter estimates start to diverge from the reference solution. For the EnKF method, at time 3, 12, and 13, the error variance is greater than 0.006, which means filter estimates deviate from true solution. Among these locations, none of 3 is in the second half of the assimilation period. For the SMC methods, at time 3, 4, 16, 22, and 33, the error variance is greater than 0.006. Among these locations, 2 out of 10 are in the second half of the assimilation period. Despite these divergences, both methods recover quickly and track the reference solution again. From the error variance variation with time, we can see the strong error growth at those phase transition locations. Furthermore, the noisy level of RMSE for the EnKF method is lower than that for the SMC methods. Table 4.2 indicates the CPU computation time and the RMSE for both methods in the case. The CPU computation time is 15.809 s for the SMC methods, while it is 142.960 s for the EnKF method. The EnKF method takes almost 9 times longer than the SMC methods. The RMSE is 1.2777, 1.1791, and 1.1207 for XI, X20, and X30 for the SMC methods, while it is 0.9404, 0.9297, and 43 4.1 Weakly Chaotic F = 5 0.9310 for XI, X20, and X30 for the EnKF method. The SMC methods are slightly worse than the EnKF method in this case. From Figure 4.2, we can see clearly that the probability density function of XI is non-Gaussian, not strongly non-Gaussian though. Both methods can assimilate it quite well, but, they both have some difficulties to reach the exact peaks of the probability density function. Table 4.3 shows the statistics of the probability density function of XI of the Lorenz 1996 model. The mean and the standard deviation for the SMC methods assimilated result are closer to the true state than that of the EnKF method. However, for the higher order moments of the probability density function, the coefficient of skewness for the SMC methods is closer to the true resolution than the EnKF method, while the coefficient of kurtosis for the EnKF method is closer to the true resolution than the SMC methods. All three quartiles in the SMC methods are better than those in the EnKF method except the range. Table 4.2: Computation Time and RMSE of Lorenz 1996 {F = 5) Assimilation Method Time (In Seconds) RMSE (XI) RMSE (X20) RMSE (X30) SMC(250 particles) 15.809 1.2777 1.1791 1.1207 44 EnKF(250 ensembles) 142.960 0.9404 0.9297 0.9310 Statistics Mean Standard deviation Coefficient of skewness Coefficient of kurtosis First Quartile Second Quartile Third Quartile Range True PDF 1.6262 2.3589 0.0602 -0.6937 -0.0816 1.5014 3.3845 11.0311 SMC(250 particles) EnKF(250 ensembles) True - SMC True - EnKF 1.6123 1.6853 0.0138 -0.0591 -0.0918 0.1198 2.4508 2.2390 -0.0397 0.0164 0.1000 0.0438 -0.4271 -0.6882 -0.2666 -0.0055 0.0664 0.0074 -0.1480 -0.0890 -0.0302 -0.1108 1.5316 1.6122 -0.0491 -0.1336 3.4336 3.5181 11.5739 -2.8826 -0.5428 13.9137 Table 4.3: Statistics of PDF of XI of Lorenz 1996 (F = 5) 4.1 Weakly Chaotic F = 5 State Estimate Error Variance 0.01 l' E„KF| 0.008 - 0.006 - 0.004 i 0.002 1" M KMh.Lk., MLMA LAid! 10 15 20 time t 25 30 35 State Estimate 35 40 Error Varianc e i i i i |' 0.006 0.004 0.002 - j SMC | - 0.008 - " 1 II f 1 Jkw wJ, Ji 10 \ » JI ' fljl LiJ 4\i 15 20 time t 25 NULLA^ 30 35 40 Figure 4.1: State estimate and error variance of XI of Lorenz 1996 model (F = 5) for EnKF and SMC methods with filter size of 250 46 4.1 Weakly Chaotic F = 5 Probability Density Function Probability Density Function 0.18 i 1 i 1 1 1 1— 0.16 0.14 - 0.12 - 0.1 - / 0.08 \ " / \\ 0.04 n - 11 0.06 0.02 \\ / " j " " \ / Tme | " ,J 1 _J ^W'-—,,„ 1 8 10 12 Figure 4.2: Probability density function of XI of Lorenz 1996 model (F — 5) for EnKF and SMC methods with, filter size of 250 47 4.2 Highly Chaotic F = 8 4.2 Highly Chaotic F = 8 When F = 8, the dynamics become highly chaotic. In Fig 4.3, both the SMC methods and the EnKF method can track the phase transitions and also in reproducing the correct amplitudes of the reference solution reasonably well, not as good as in the case with F = 5 though. If we still divide the whole period into two. There is not much difference between the two halves in performance. When degree of chaos increases, the difficulty to track the true resolution increases. From the error variance variation with time, we also can see the strong error growth at those phase transition locations. For the EnKF method, at more than 1/3 of the whole time period, the error variance is greater than 0.006, which means filter estimates deviate from true solution quite frequently. These locations are distributed in the whole assimilation period. For SMC methods, at more than 1/3 of the whole assimilation period, the error variance is greater than 0.006. Those locations also exist in the whole assimilation period. Compared to Fig 4.1, the data assimilation performance for both methods is worse than that in case F — 5. The noisy level of RMSE is also greater than that in case F = 5. Despite these divergences, both methods recover quickly and track the reference solution again in general. Table 4.4 indicates the CPU computation time and the RMSE for both methods in the case. The CPU computation time is 21.968 s for SMC methods, while it is 148.003 s for EnKF method. The CPU computation time for the SMC methods increase significantly from 15.809 s to 21.968 s with the degree of chaos from F = 5 to F = 8. The EnKF method takes almost 9 times longer than the SMC methods. The RMSE is 1.8116, 1.8128, and 1.8967 for X I , X20, and X30 for the SMC methods, while it is 1.6783, 1.6406, and 1.8820 for X I , X20, and X30 for the EnKF method. The EnKF method outperforms the SMC methods. From Figure 4.4, the probability density function of XI is still weakly nonGaussian. Both methods can assimilate the probability density functions quite well in general, however, the EnKF method cannot reach the exact peak of the probability density function, while the SMC methods produced two false peaks of the probability density function. 48 4.2 Highly Chaotic F = 8 Table 4.5 shows the statistics of the probability density function of XI of Lorenz 1996 model. The mean, the standard deviation, and the coefficient of skewness for the SMC methods are closer to the true state than those of the EnKF method. However, the coefficient of kurtosis is -0.598373, -0.502227, and -0.661721, which shows the EnKF method assimilates better. Besides, the first and third quartiles and ranges estimate in the EnKF method are better than those from the SMC methods except the second quartile (median). Table 4.4: Computation Time and RMSE of Lorenz 1996 (F = 8) Assimilation Method Time (In Seconds) RMSE (XI) RMSE (X20) RMSE (X30) SMC (250 particles) 21.968 EnKF (250 ensembles) 148.003 1.8116 1.8128 1.8967 1.6783 1.6406 1.8820 49 Statistics Mean Standard deviation Coefficient of skewness Coefficient of kurtosis First Quartile Second Quartile Third Quartile Range True PDF SMC(250 particles) EnKF(250 ensembles) True - SMC True - EnKF 1.6262 1.6123 1.6853 0.0138 -0.0591 2.3589 2.4508 2.2390 -0.0918 0.1198 0.0602 -0.0397 0.0164 0.0438 0.1000 -0.6937 -0.4271 -0.6882 -0.2666 -0.0055 -0.7972 -0.7196 -0.0965 -0.0189 -0.8161 1.8654 1.9207 -0.0086 1.8568 -0.0639 4.5321 4.5672 0.0561 4.4760 -0.0351 19.8032 19.3982 -0.6533 20.4565 0.4050 Table 4.5: Statistics of PDF of XI of Lorenz 1996 (F = 8) 4.2 Highly Chaotic F = 8 State Estimate 15 10 5 0 -5 -10 0 5 J I I L 10 15 20 timet 25 30 35 40 25 30 35 40 25 30 35 40 Error Variance 0.01 0.008 0.006 0.004 0.002 0 5 10 15 -j | | 5 10 15 20 time t State Estimate 15 r 10 5 0 -5 -10 0 20 time t Error Variance time t Figure 4.3: State estimate and error variance of XI of Lorenz 1996 model (F = 8) for EnKF and SMC methods with filter size of 250 51 4.2 Highly Chaotic F = 8 Probability Density Function Probability Density Function Figure 4.4: Probability density function of A'l of Lorenz 1996 model (F = 8) for EnKF and SMC methods with filter size of 250 52 4.3 F u l l y T u r b u l e n t F = 16 4.3 Fully Turbulent F = 16 If we continue to increase the external forcing F, when F = 16, the dynamics become fully turbulent. In Fig 4.5, both the SMC methods and the E n K F method can track the phase transitions and also reproduce the correct amplitudes of the reference solution. Obvio\isly it is worse t h a n the first two cases with F = 5 and F = 8. Since the degree of chaos increases, the whole dynamic process becomes highly unstable and fully turbulent, which increases the difficulties for the data assimilation system. The locations where the filter estimates start to diverge from the reference solution become more frequent t h a n in the first two cases. These locations appear also in the whole assimilation period. From the error variance plot with time, we can see the strong error growth at those phase transition locations. The noisy level of RMSE increases significantly. We cannot use 0.006 as a standard any more. In this case, we choose 0.01 as the standard. For the E n K F method, at more than 1/2 of the whole time period, the error variance is greater than 0.01, which means filter estimates deviate from the true solution significantly and frequently. These locations are in the whole assimilation period. For SMC methods, at more than 1/2 of the assimilation time period, the error variance is greater than 0.01. Table 4.6 indicates the CPU computation time and the RMSE for both meth- ods in this case. The CPU computation time is 33.493 s for SMC methods, while it is 161.101 s for E n K F method. The E n K F method takes almost 4 times longer than the SMC methods. One interesting feature is that the CPU computation time used by the SMC methods increased significantly with the degree of chaotic nature, while the E n K F method does not. This may indicate that the SMC methods depend on model chaotic nature to some extent. The RMSE 3.9858, and 3.5495 for XI, is 3.9114, X20, and X30 for SMC methods, while it is 3.3883, 3.7520, and 3.7619 for XI, X20, and X30 for E n K F method. Both RMSE from two methods are quite similar. From Figure 4.6, we can see clearly that the probability density function is non-Gaussian, the same as the previous two cases. Both methods can assimilate it quite well, but, they both have some difficulties to reach the exact peaks of the probability density function. 53 4.3 Fully Turbulent F = 16 Table 4.7 shows the statistics of the probability density function of XI of Lorenz 1996 model. The difference of the mean between true state and estimate for the SMC methods are much larger than that of the EnKF method as well as the coefficients of skewness and kurtosis. While the standard deviation of the estimate for the SMC methods are closer to that of the true state than the EnKF. Besides, the first and third quartiles and ranges estimate in the EnKF method are better than those from the SMC methods except the median. The range estimate in the SMC method is almost 10 times larger than that in the EnKF estimate, which means the SMC methods generate some extreme values in the assimilated process, more unstable than the EnKF. Table 4.6: Computation Time and RMSE of Lorenz 1996 (F = 16) Assimilation Method Time (In Seconds) RMSE (XI) RMSE (X20) RMSE (X30) SMC(250 particles) 33.493 3.9114 3.9858 3.5495 54 EnKF(250 ensembles) 161.101 3.3883 3.7520 3.7619 Statistics Mean Standard deviation Coefficient of skewness Coefficient of kurtosis First Quartile Second Quartile Third Quartile Range True PDF SMC (250 particles) EnKF(25Q ensembles) True - SMC True - EnKF 3.5728 3.2766 3.5611 0.2961 0.0116 6.2902 6.2494 6.0941 0.0408 0.1961 0.0651 0.0018 0.0669 -0.0039 0.0709 -0.3906 -0.4688 -0.4474 0.0782 0.0568 -0.7733 -1.1512 -0.9289 0.1556 0.3779 3.0752 3.5134 3.5129 0.4377 -0.0005 8.0694 7.9025 7.7888 0.1137 -0.1669 36.4173 40.1501 35.8220 -3.7328 0.5953 Table 4.7: Statistics of PDF of XI of Lorenz 1996 {F = 16) 4.3 Fully Turbulent F = 16 State Estimate Error Variance 0.05 0.04 0.03 0.02 0.01 10 mm 15 20 time t 111! ill I m II j 25 30 35 30 35 State Estimate Error Variance 0.05 0.04 0.03 0.02 0.01 LI Mil I 15 A A iM yi 20 time t AA 25 40 Figure 4.5: State estimate and error variance of XI of Lorenz 1996 model (F 16) for EnKF and SMC methods with filter size of 250 56 4.3 Fully Turbulent F = 16 Probability Density Function -20 -15 -10 -5 0 5 10 15 20 25 30 Probability Density Function -30 -20 Figure 4.6: Probability density function of XI of Lorenz 1996 model (F — 16) for EnKF and SMC methods with filter size of 250 57 4.4 Filter Size Comparison 4.4 Filter Size Comparison In the EnKF method, the error statistics such as mean and variance (covariance) are represented using the model state ensemble. In the SMC methods, the error statistics of the probability density function are estimated using a set of sample particles. In both methods, the larger the ensemble size, the better estimate of the error statistics. As the ensemble size approaches infinity, the estimate of error statistics reaches the optimal estimate. In practice, it is impossible to employ an ensemble of infinite model members or sample particles to perform data assimilation. In realistic applications, especially in the atmospheric and oceanic fields, the ensemble size varies from 10 to 1000, because the computational cost limits the ensemble size. Since the ensemble size is limited, the estimate based on limited ensembles is not optimal, it is suboptimal. In this section, we compare the different ensemble size of two data assimilation methods. The filter size varies from 5, 10, 25, 50, 75, 100, 250, 500, 750 and 1000. Figure 4.7 shows the effect of ensemble size on data assimilation in the Lorenz 1963 model. Figures 4.8 and 4.9 show the effect of ensemble size on data assimilation in the Loren 1996 model. It is clear that TV = 100 is the critical point of the ensemble size in both the EnKF method and the SMC methods. The RMSE is quite large and oscillates when the ensemble size N is smaller than 100, while the RMSE is quite stable when the ensemble size N is larger than 100. It indicates that in Lorenz models, the ensemble size of 100 is sufficient to achieve reasonably good estimates. However, in Fig 4.7, the ensemble size of the SMC methods decreases more quickly than that of the EnKF method when the ensemble size is smaller than 100, and it still decreases slowly when the ensemble size is larger than 100, while EnKF does not. In Fig 4.8 and 4.9, the ensemble size of the EnKF method decreases more quickly than that of the SMC methods when the size is smaller than 100. After N greater than 100, both of them are stable. However, the number of particles needed to perform a successful Sequential Monte Carlo methods increases exponentially with the size of the assimilated dynamic system (Snyder et «•!., 2008). It is not practical to implement the SMC 58 4.4 Filter Size Comparison methods directly to atmospheric and oceanic models at present, since the model has 107 degrees of freedom. 59 as o 10 -©• ensemble size 10' 10 Figure 4.7: Filter size of Lorenz 1963 model for EnKF and SMC methods 10 10 ensemble size Lorenz 63, Variable z ensemble size Lorenz 63, Variable y Lorenz 63, Variable x P 2 5" o O o N CO W CD >-S r? • >*= V- : I :T ' 10 ensemble size ^ 10' 10 ensemble size Lorenz 96, F = 8.0, Variable X30 10 : i ~l • - * — , ....:...:.:..:..:.:. f 10 10° ..—@.~ — ^ s ^ ™ - - - ^ . - ^ ' — 1 „—£ij.„__ E n K F 1 Figure 4.8: Filter size of Lorenz 1996 {F = 8) model for EnKF and SMC methods 10 a: 2 CO ° 10 a. 2 CO ° 10 ensemble size Lorenz 96, F = 8.0, Variable X20 10' Lorenz 96, F = 8.0, Variable X1 P O En" *"S_ 3 O o CD w 5" CD *-i to *-~ ©-- •Sjt' 10 ; i 10 TT*>, ~T" ensemble size 10' ensemble size Lorenz 96, F = 16.0, Variable X30 10 ~a^ ; ~ EnKF I ; ;; -'•^i-.'i 10 10 m -t- Figure 4.9: Filter size of Lorenz 1996 (F = 16) model for EnKF and SMC methods 10 10 10 10 10 ensemble size Lorenz 96, F = 16.0, Variable X20 10 Lorenz 96, F = 16.0, Variable X1 .; 5* o P O o 5' a CO a>s Chapter 5 Discussions and Conclusions Although significant progress has been made in the data assimilation field, it is still difficult to deal with the nonlinear and non-Gaussian state estimation problems, which cannot be resolved with the traditional data assimilation methods. The Sequential Monte Carlo (SMC) methods are among the latest innovations that attempt to bridge the existing gap between the Gaussian dynamics estimation and the non-Gaussian dynamics estimation in the data assimilation process. It has been shown that they perform quite well in the complex practical scenarios. In the Kalman Filter framework, nonlinearity and non-Gaussianity state estimation problems cannot be solved theoretically; they can only employ the Gaussian probability density function to characterize the non-Gaussian probability density function. Therefore, to tackle this estimation problem of the nonlinear, non-Gaussian dynamic system, SMC methods directly approximate the probability density function (PDF) associated with the dynamic system with finite samples, which is a powerful tool to characterize the uncertainty of the dynamic system instead of only mean and covariance of the system. Different order statistical moments such as mean and variance can be calculated directly from the probability density function. In this thesis, the SMC methods were used to perform data assimilation in strongly nonlinear dynamic systems, the Lorenz 1963 and 1996 models. Comparison in the same scenarios is made to the EnKF data assimilation method. In the experiment design with the three variables Lorenz 1963 model, we choose 63 4 scenarios: observation intervals with 0.50 and 0.25, initial non-Gaussian error probability density functions of Beta distribution and Gamma, distribution. In the experiment design with the 40 variables Lorenz 1996 model, we use 3 scenarios: the weakly chaotic dynamics with external forcing F — 5, the highly chaotic dynamics with external forcing F — 8, and the fully turbulent dynamics with external forcing F — 16. Both Lorenz models are relatively low dimensional simplified dynamic models, but they are highly nonlinear dynamics of stochastic nature. Different model parameters represent different situations in the realistic world. However, in the experiments with Lorenz 1963 and 1996 models as test beds, the SMC methods perform almost as well as the EnKF method, which does not outperform the EnKF method as it is in the theoretical aspect. The reasons may be: 1) the non-Gaussianity of all the probability density functions either in the Lorenz 1963 dynamics or the Lorenz 1996 dynamics are not significantly strong. Those probability density functions are not obviously bimodal or multimodal, though the coefficients of skewness and kurtosis exist. Since the Gaussian probability density function is completely characterized by the first two moments, it is clear that one can obtain the probability density function of a Gaussian process from its mean and covariance information. That explains why the EnKF method generates good results as well as the SMC methods. 2) The criteria RMSE and assimilated probability density function are calculated only from mean values of the true resolution and the assimilated resolution. There are no criteria considering the higher order moments of the error statistics, which is the advantage of the SMC methods. Since all the data assimilation work are based on Bayesian statistics, the sensitivity analysis of different prior probability density function could be done in the future work to generate the optimal prior probability density function for the data assimilation system. Although Sequential Monte Carlo methods are suitable for the most general case: nonlinear, non-Gaussian dynamics, for linear, Gaussian dynamics, Kalman Filter will still be the first approach for its easy implementation; for linear or nonlinear, weakly non-Gaussian dynamics as in our experiments, Extended Kalman Filter and Ensemble Kalman Filter could be employed to achieve reasonably 64 good estimate; for highly nonlinear, strongly non-Gaussian dynamics, Sequential Monte Carlo methods will the best choice to fully capture the non-Gaussianity of the dynamics. In the Meantime, the choice of the assimilation method is limited by the computation power. One interesting feature found in our experiments is that the computational cost of the SMC methods is much cheaper than the EnKF method. This seems inconsistent with the perception that the SMC methods are very expensive computationally. This is mainly because of relatively low dimensionality of the dynamic systems used in this study. As the dimensionality increases, the ensemble size required for SMC methods increases exponentially (Snyder el ai, 2008). Within the theoretical insight, Beugtssoii et al. (2008) point out that Sequential Monte Carlo (SMC) methods may fail in large scale dynamic systems. Their simulations suggest that the convergence to unity occurs unless the ensemble grows super-exponentially in the system dimension. At present there is no SMC application in realistic atmospheric and oceanic models because of the high dimension of dynamic models. This is an obstacle to high-dimensional SMC data assimilation. According to Snyder et al. (2008), Gaussian errors, simulations indicate that the required ensemble size scales exponentially with the state dimension. In his example, the particle filter requires at least 1011 members when applied to a 200-dimensional state for the posterior mean from the particle filter to have an expected error smaller than either the prior or the observations. However, in some cases, the system model has some substructure which can be tractable and analytically marginalized out. The advantage of this strategy is that it can drastically reduce the size of the space over which we need to sample and reduce the filter size. Marginalizing out some of the variables is a process which is called Rao-Blackwellisation, because it is related to the Rao-Blackwell formula: see (Casella & Robert, 1996) for a general discussion. In this thesis, the properties and capabilities of the SMC methods is investigated and compared to the EnKF method using the low dimensional, highly nonlinear dynamic systems, the Lorenz 1963 and 1996 models. Despite the interesting fact that the EnKF method performs as well as the SMC methods in the highly nonlinear dynamics, the SMC methods have theoretical advantages and 65 potential practical significance, which is helpful when we design data assimilation systems for nonlinear. non-Gaussian realistic models. 66 Appendix A F O R T R A N Code for the E n K F Method The EnKF algorithm implemented in this thesis is explained in Evensen (2003). For the detailed description, please refer to Evensen (2003). The following FORTRAN code provides a detailed implementation of the EnKF analysis scheme. It assumes access to the Numerical Algorithms Group (NAG) Numerical Libraries, which specializes in the provision of software for the solution of mathematical, statistical and data mining problems. 1 subroutine EnKF(R, OBS, PCL, x h a t ) 3 i m p l i c i t none r. 0 7 integer : : ndim / dimension integer :: nrens ! number of integer :: nrobs ! number of integer : : l a t , Ion of model state ensemble members observations 9 ! number 67 of Latitude , Longitude integer : : l a t . o b s , lon_obs / number integer :: ndim.T / grid l o n = 3 , l a t _ o b s = l, of Latitude of the , Longitude whole parameter (lat=l, parameter ( n d i m = l a t *lon , n r e n s = 2 5 0 , n r o b s = l) real PCL(n r e n s , real domain lon_obs=3) / input ensemble OBS(ndim) ! input obs real Xhat (ndim) ! output analysis real A(ndim , n r e n s ) / Ensemble matrix real X_A(ndim,l) ! A n a I y s i s m a trix real X_F(ndim,l) ! forecast real Y(ndim,l) ! observations real XJVl(ndim,l) ! Ensemble mean :: X_A2(ndim , n r e n s ) ! updated all :: Y_P2(ndim , n r e n s ) ! perturbed real :: K(ndim , ndim) / Kalman real :: D(ndim ,1) ! Innovation real :: P_A(ndim , ndim) / analysis error covariance : : P_F(ndim ,ndim) ! forecast error covariance :: P_R(ndim ,ndim) ! observations real nd i m) matrix matrix matrix matrix matrix member analysis matrix real observations matrix Gain matrix matrix matrix real matrix real 68 error covariance matrix real :: H ( n d i m , ndim) / observation real :: / identity ! Local real, I ( n d i m , ndim) operator matrix variables a l l o c a t a b l e , dimension (:, :) :: XI, X2, X3, X4, X5, X6, & X7, X8, X9, integer :: parameter (NMAX=lat *lon , LDA=NMAX, LWDRK=64*NMAX) allocatable, dimension (:) :: WORK integer, allocatable, dimension (:) :: IPIV external :: F07ADF, F07AJF, F06YAF integer :: t, real :: a l p h a = 1.0, j , kappa, u beta=0.0, a l l o c a t a b l e , dimension character :: d i m e n s i o n (8) real s t a r t , finish , R integer parameter (: , :) t h e t a = 1 . 0 / ( n r e n s —1) :: X_F3 d a t e , t i m e , zone integer, :: Xll NMAX, LDA, LWORK, INFO, N real, real, matrix :: values :: N02, N03 (N02 = n d i m * n r e n s , N03 = ndim + n r e n s ) integer :: I F A I L , IGEN real :: X02(N02) , X03(N03) integer :: ISEED(4) external G05KCF, G05LAF Assimilation Cycle Starts, do t = 1, ndim do j = 1, A(t, nrens j) = PCL(j, t) e n d do e n d do Perturbed Initialize Observation the seed to a un—repeat able sequence I S E E D ( l ) = 1762543 1SEED(2) = 9324783 ISEED(3) = 42344 ISEED(4) = 742355 IGEN identifies the stream, IGEN = 1 c a l l G05KCF(IGEN, ISEED) IFAIL = 0 c a l l G 0 5 L A F ( 0 . 0 e 0 , R, N02, X02, IGEN, ISEED, 70 IFAIL) 97 Y.P2 = d o u = 1, r e s h a p e (X02, (/ndim, nrens/)) nrens 101 / *#*#*#*****#*#*****#*******#****#*#**#*******#** 103 d o t = 1, ndim X_F(t, 1) = A ( t , 105 end do 107 d o t = 1, ndim Y(t, 109 u) 1) = O B S ( t ) + Y J P 2 ( t , u) e n d do 111 / generate the identical matrix ! II(ndim , ndim) 113 d o t = l,ndim 115 do j=l,ndim if ( t = j ) then 117 H(t,j)=1.0 else 119 H(t,j)=0.0 end 121 if end do e n d do 123 125 / Compute Background error covariance 71 P-F / ^s^^*^^^^^:^^*;)!:^*^:^^***************** / A(ndim , nrens) 129 / XJA(ndim, 1) 131 d o t = l,ndim ************************* 127 ensemble mean X_M(t , l ) = s u m ( A ( t , : ) ) / n r e n s 133 end d o 135 / background error ! X-F(ndim, ndim) covariance 137 / 139 error do t = l , 141 nrens allocate (X_F3(ndim, 1)) X.F3 ( 1 : ndim , 1) = A ( l : n d i m , t ) - X_M(1: ndim , 1) 143 / error covariance 145 / X.F3(ndim,l) ! X.F3(ndim,l) 147 / F-F(ndim, 149 allocate ndim) ( X l l (ndim , n d i m ) ) call F06YAF('n', 151 't', ndim, ndim, a l p h a , X_F3 , ndirn , & X.F3, ndim, & 153 155 beta, Xll, ndim) P.F=X11 + P_F 72 1, & 157 deallocate (X_F3) deallocate (Xll) 159 end do 161 P_F=P_F*theta 163 165 / step 1. compute innovation d 167 / 169 / 1) compute H(ndim, X1=H"X„F ,ndim) ! X.F(ndim,l) 171 / Xl(ndim,l) 173 allocate (Xl(ndim , 1) ) call F06YAF('n', 175 'n', ndim, 1, n d i m , & a l p h a , H, ndim, & X_F, ndim, & 177 beta, 179 / ! 2) compute XI, ndim) d=yO-II'X.F=yO~Xl Y(ndim,l) 181 / XI (ndim, 1) ! D(ndim, 1) 183 D=Y-X1 185 deallocate (XI) 73 1 0 J / ?Jc; 5^ ^H ^K ^ ^ ^ ^ ^ ^ H1^ ^ H< ^ >K ^f; 3^; ^ ^; 3^ s^c ^ ^ H^ ^ ^ H^ ^ ^ H ^ ^ ^ ^ K ^ ^ H ^ ^ ^ H ^ ^ ^ H^ H^ ^ ^ ^ ^ ^ ^ ^ ^ K ^ s f c / step 189 . 2. compute gain matrix K ^^^^^^^^^^^^^^^^^^^^^^^^^HS^^^^^^^^^^^^H^ 191 / 1) compute ! P.F(ndim, Pf'H' ndim) 193 / H(ndim, ndim) ! X2(ndim, ndim) 195 allocate 197 (X2(ndim , n d i m ) ) call F06YAF('n\ 't\ ndim, ndim, ndim, & a l p h a , P_F , ndim , & 199 H, ndim., & b e t a , X2, ndim) 201 / 2) compute 203 / H(ndim, H'P.F ndim) ! P-F(ndim, ndim) 205 / X3(ndim, ndim) 207 allocate (X3(ndim , ndim) call F06YAF('n', 209 'n', ndim, ndim, ndim, & a l p h a , H, ndim , & P_F , ndim , & 211 b e t a , X3, ndim ) 213 .' 3) compute XS'H' ! X3(ndim , ndim) 215 / H(ndim, ndim) / (about 74 20 min) 5jc ^c ?f; s|= =|c ^f; sf; >(< ^< ^c s}c ^c / X4 (ndirn, ndim) 217 allocate 219 (X4(ndim , n d i m ) ) call F06YAF('n', 't', ndim, ndim, ndim, & a l p h a , X3, ndim , & 221 II, ndim, & b e t a , X4, ndim ) 223 deallocate 225 / 4) ! P.R compute is (X3) R+X4 the diagonal matrix 227 do t = l,ndim 229 do j = l , n d i m if (t=j) 231 then P_R(t,j)=R else 233 P.R(t,j)=0.0 end 235 if end do end do 237 / X5 (n dim , n dim ) 239 allocate 241 (X5(ndim , n d i m ) ) X5=PJl+X4 deallocate (X4) 243 / 5) compute inverse of X5 245 75 allocate (WORK(LWORK)) 247 allocate (IPIV(NMAX)) 249 / Factorize X5 N = ridim 251 call if 253 / F07ADF(N,N,X5,LDA,1PIV,INFO) (INFO.EQ.O) then inverse X5 compute of c a l l F07AJF(N,X5,LDA, IPIV ,WORK,LWORK,INFO) 255 257 endif d e a l l o c a t e (WORK) deallocate (IPIV) 259 / 6) compute 261 / X2(ndim ! X5-inverse 263 / Kfndim, 265 gain K= ,ndim) =(ndim, ndim) ndim) call F06YAF('n', 'n', n d i m , ndim, ndim, & a l p h a , X2, ndim , & 267 X5, n d i m , & b e t a , K, 269 deallocate (X2) deallocate (X5) ndim) 271 273 / Step 3. compute X-a, 275 76 / 1) compute 277 / K(ndim, K"d ndim) ! D(ndim, 1) 279 / X7(ndim,l) 281 allocate (X7(ndim,l)) call F06YAF('n', 283 'n', ndim, 1, ndim, & a l p h a , K, ndim, fe D, ndim, 285 beta, 287 / 2) compute & X7,ndim) X„a ! X.F(ndim,l) 289 / X.A(ndim,l) 291 X_A=X_F+X7 deallocate (X7) 293 do t = 1., ndim 295 X_A2(t , u) = X_A(t , 1) end d o 297 299 / Step 4. compute P„A 301 / 1) compute K'H 303 / K(ndim , ndim,) ! II'(ndim, ndim) 305 / X8(ndim, ndim) 77 307 allocate (X8(ndim , n d i m ) ) call F06YAF('n', 309 'n', ndim, ndim, ndim, & a l p h a , K, ndim, & H, ndim , 311 & b e t a , X8, ndim) 313 / 2) compute I—K"H 315 / I (ndim, ndim) ! X9(ndim , ndim) 317 allocate 319 (X9(ndim , n d i m ) ) X9=I-X8 deallocate (X8) 321 / 3) compute P„a 323 .' P. A (ndim, ndim) ! P-F(ndim,, ndim) 325 / X9(ndim, 327 ndim) call F06YAF('n', 'n', ndim, ndim, n d i m , & a l p h a , X9, ndim, & 329 P_F, ndim, & b e t a , P_A, 331 deallocate ndim) (X9) end do 78 do t = 1, ndim Xhat(t) = sum(X_A2(t , :))/nrens end do do t = 1, ndim do j = 1, nrens PCL(j, t) = X_A2(t, j) end do end do P.F = 0.0 P_A = 0.0 return end subroutine EnKF 79 Appendix B F O R T R A N Code for t h e SMC Methods The SMC algorithm implemented in this thesis is explained in Gordon et al. (1993). For the detailed description, please refer to Gordon et al. (1993). The following FORTRAN code provides a detailed implementation of the SMC analysis scheme. subroutine Particle .Filter integer i, :: (R, y , x , xhat) j , M parameter (M = 250) r e a l , d i m e n s i o n (M, 3) :: x , e , q_w, qn, : : y , PROB00, x h a t , q.sum real , dimension (3) r e a l , parameter : : pi = 10 real 12 / step 1. ind , tempOl 3.1415926 :: R 80 tempi 1 , prob / Calculate difference between Observations - x(i, :) and model 14 do i = 1 , M 16 e(i, :) = y ( : ) end do 18 / step 2. 20 / Calculate weights from Gaussian Distribution do i = 1, M 22 call Std_Normal(e(i , :), q.w(i, : ) , R) end do 24 / step 3. 26 / summation of all the weights do i = 1. , 3 28 q.sum ( i ) = sum (q_w (: , i ) ) end do 30 / step 4- 32 / Normalize importance weights do j = 1, 3 34 do i = 1, M q_w (i , j ) = q.w (i , j ) / q.sum (j ) 36 end do end do 38 / step 5. 40 / Get particle * weight do j = 1, 3 42 do i = 1 , M 81 forecasts qn ( i . J) = q- w (i . J ) 44 * x (i > J ) end do end do 46 / step 6. ! Get analysis do i = 1, 3 50 x h a t (i ) = sum (qn (: , i ) ) end do 52 .' step 54 7 ** Resampling ** call (q_w, resampling ind) do j = 1, 3 do i = 1, M 58 x(i , j ) = x ( i n d ( i , j ) , j) end do 60 end do 62 return end s u b r o u t i n e Particle.Filter 64 66 / Resampling Process 68 subroutine resampling implicit none integer :: M (q_w, ind) 70 72 parameter (M = 250) 82 real, d i m e n s i o n (M, 3) :: q_w, ind , qc , x x x , y y y , u , real, d i m e n s i o n (3) : : ssum, integer : : i , j , k, n ssum = 0.0 ranOl do j = 1, 3 do i = 1, M ssum (j ) = ssum (j ) + q_w (i , j ) qc ( i , j ) = ssum ( j ) end do end do call random_seed () call random_number (ranOl) do j = 1 , 3 do i = 1, M u(i , j ) = (i - 1 + raii01(j))/M end do end do do n = 1 , 3 do j = 1, M k = 1 do w h i l e (qc(k, n) < u ( j , n ) ) k = k + 1 end do tempOl (j , n) = k end do 83 tempOl e n d do 104 ind = tempOl 106 return end s u b r o u t i n e resampling 108 110 / END 84 References ANDERSON, B.D.O. & M O O R E , J.B. (1979). Optimal Filtering. Information and System Science Series, Prentice Hall, Englewood Cliffs, NJ, USA. 8 M.S., MASKELL, S., GORDON, N. & C L A P P , T. (2002). A tutorial on particle filters for online nonlinear/non-Gaussian Bayesian tracking. IEEE Transactions on Signal Processing, 50(2), 174-188. 8, 13 ARULAMPALAM, BENGTSSON, T., BICKEL, P. & Li, B. (2008). Curse-of-dimensionality revisited: Collapse of the particle filter in very large scale systems. Probability and Statistics: Essays in Honor of David A. Freedman, IMS Collections, 2, 316334. 65 BURGERS, G., VAN LEEUWEN, P.J. & EVENSEN, G. (1998). Analysis scheme in the Ensemble Kalman Filter. Monthly Weather Review, 126, 1719 1724. 3, 7,9 G. & ROBERT, C.P. (1996). Rao-Blackwellisation of sampling schemes. Biometrika, 83(1), 81-94. 65 CASELLA, COHN, S.E. (1997). An introduction to estimation theory. Journal of the Meteorological Society of Japan, 75, 147-178. 6 COURTIER, VIC, D., P., ANDERSSON, HAMRUD, M., E., HECKLEY, W., HOLLINGSWORTH, A., PAILLEUX, J., VASILJE- RABIER, F. & FISHER, M. (1998). The ECMWF implementation of three-dimensional variational assimilation (3D-Var). I: Formulation. Quarterly Journal of the Royal Meteorological Society, 124(550), 1783 1808. 2, 3 85 REFERENCES DIMET, F.X.L. & TALAGRAND, O. (1986). Variational algorithms for analysis and assimilation of meteorological observations: Theoretical aspects. Tellus, 38A, 97-100. 2 DOUCET, A., GODSILL, S.J. & ANDRIEU, C. (2000). On sequential Monte Carlo sampling methods for Bayesian filtering. Statistics and Computing, 10(3), 197208. 8, 13 D O U C E T , A., DE FREITAS, N. & GORDON, N. (2001). Sequential Monte Carlo Methods in Practice. Springer-Verlag, New York, USA. 7, 8 EFRON, B. (1979). Bootstrap methods: Another look at the jackknife. Annals of Statistics, 7(1), 1-26. 12 G. (1994). Sequential data assimilation with a nonlinear quasigeostrophic model using Monte Carlo methods to forecast error statistics. Journal of Geophysical Research, 99(C5), 10143-10162. 3, 9, 15 EVENSEN, EVENSEN, G. (1997). Advanced data assimilation for strongly nonlinear dynam- ics. Monthly Weather Review, 125, 1342-1354. 23 EVENSEN, G. (2003). The Ensemble Kalrnan Filter: Theoretical formulation and practical implementation. Ocean Dynamics, 53, 343 367. 3, 7, 29, 67 G. & VAN LEEUWEN, P.J. (1996). Assimilation of geosat altimeter data for the Agulhas current using the Ensemble Kalman Filter with a quasigeostrophic model. Monthly Weather Review, VIA:, 85-96. 3, 7 EVENSEN, N.J., SALMOND, D.J. & SMITH, A.F.M. (1993). Novel approach to nonlinear/non-Gaussian Bayesian state estimation. IEE proceedings. F, Radar and signal processing, 140, 107-113. 4, 8, 13, 15, 80 GORDON, H.OUTEKAMER, P.L. & MITCHELL, H.L. (1998). Data assimilation using an Ensemble Kalman Filter technique. Monthly Weather Review, 126, 796-811. 3, 7 HOUTEKAMER, P.L. & MITCHELL, H.L. (2005). Ensemble Kalman filtering. Quarterly Journal of the Royal Meteorological Society, 131, 3269-3289. 10 86 REFERENCES HOUTEKAMER, P . L . , MITCHELL. H.L., FELLER™, G., BUEHNER, M., CHARRON, M., SPACER, L. & HANSEN, B. (2005). Atmospheric data assimila- tion with an Ensemble Kalman Filter: Results with real observations. Monthly Weather Review, 133, 604-620. 2 ISARD, M. & BLAKE, A. (1998). A smoothing filter for condensation. Proceed- ings of 5th European Conference on Computer Vision, 1, 767-781. 13 A.H. (1970). Stochastic Processes and Filtering Theory. Academic Press, New York, USA. 3, 8 JAZWINSKI, S.J. & UHLMANN, J.K. (1996). A General Method for Approximating Nonlinear Transformations of Probability Distributions. Technical Report, RRG, Dept. of Engineering Science, University of Oxford. 29 JULIER, KALMAN, R.E. (1960). A new approach to linear filtering and prediction. Trans- actions of the ASME. Series D, Journal of Basic Engineering, 82, 35-45. 2, 3, 8 G. (1996). Monte Carlo filter and smoother for non-Gaussian nonlinear state space models. Journal of Computational and Graphical Statistics, 5(1), 1-25. 8, 13 KITAGAWA, E., RABIER, F., KELLY, G. & MAHFOUF, J.F. (2000). The ECMWF operational implementation of four-dimensional variational assimilation - Part III: Experimental results and diagnostics with operational configuration. Quarterly Journal of the Royal Meteorological Society, 126, 1191 1215. KLINKER, 2 Liu, J.S. & CHEN, E. (1998). Sequential Monte Carlo methods for dynamics systems. Journal of the American Statistical Association, 93, 1032-1044. 8, 13 LORENC, A.C. (2003). The potential of the Ensemble Kalman Filter for NWP: A comparison with 4D-Var. Quarterly Journal of the Royal Meteorological Society, 129, 3183-3203. 3, 7 87 REFERENCES LORENZ, E.N. (1963). Deterministic non-periodic flow. Journal of the Atmo- spheric Sciences, 20, 130-141. 15 LORENZ, E.N. (1996). Predictability: A problem partially solved. Proceedings of the Seminar on Predictability, ECMWF, 1, 1-18. 41 J.F. k, RABIER, F. (2000). The ECMWF operational implementation of four-dimensional variational assimilation - Part II: Experimental results with improved physics. Quarterly Journal of the Royal Meteorological Society, 126,1171-1190. 2 MAHFOUF, M A J D A , A.J. & HARLIM, J. (2008). Filtering nonlinear dynamical systems with linear stochastic models. Nonlinearity, 21, 1281 1306. 41 MILLER, R.N., GHIL, M. & GAUTHIEZ, F. (1994). Advanced data assimilation in strongly nonlinear dynamical systems. Journal of the Atmospheric Sciences, 51, 1037-1056. 18 E. (1962). On estimation, of a probability density function and mode. Annals of Mathematical Statistics, 33, 1065-1076. 24 PARZEN, RABIER, F., JARVINEN, H., KLINKER, E., MAHFOUF, J . F . & SIMMONS, A. (2000). The ECMWF operational implementation of four-dimensional variational assimilation - Part I: Experimental results with simplified physics. Quarterly Journal of the Royal Meteorological Society, 126, 1143-1170. 2 RUBIN, D.B. (1988). Using the SIR algorithm to simulate posterior distributions. Bayesian Statistics, Oxford University Press, 3, 395 402. 12 SCHMIDT, S.F. (1966). Application of state-space methods to navigation prob- lems. Advances in Control Systems, 3, 293-340. 3, 7, 8 T.B. (2006). Estimation of Nonlinear Dynamic Systems - Theory and Applications. Linkoping Studies in Science and Technology, Ph.D. Dissertation. SCHON, 7, 8, 11, 13 88 REFERENCES SILVERMAN, B.W. (1986). Density Estimation for Statistics and Data Analysis. Monographs on Statistics and Applied Probability, Chapman and Hall, London. 24 SMITH, A.F.M. & GELPAND, A.E. (1992). Bayesian statistics without tears: A sampling resampling perspective. American Statistician, 46(2), 84-88. 12 SMITH, G.L., SCHMIDT, S.F. & M C G E E , L.A. (1962). Application of statistical filter theory to the optimal estimation of position and velocity on board a circumlunar vehicle. Technical Report TR, NASA, R-135. 3, 7, 8 C., BENGTSSON, T., BICKEL, P. & ANDERSON, J. (2008). Obstacles to high-dimensional particle filtering. Monthly Weather Review, 136, 4629 4640. 58, 65 SNYDER, TALAGRAND, O. (1997). Assimilation of observations, an introduction. Journal of the Meteorological Society of Japan, 75, 191-209. 1 TIPPETT, M.K., WHITAKER, ANDERSON, J.L., BISHOP, C.H., HAMILL, T.M. & J.S. (2003). Ensemble square root filters. Monthly Weather Re- view, 131, 1485-1490. 3, 7 89