Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC2964148

Formats

Article sections

- Abstract
- 1. Introduction
- 2. FFT EnKF
- 3. Morphing EnKF
- 4. Epidemic model
- 5. Computational results
- 6. Conclusion
- References

Authors

Related links

Procedia Comput Sci. Author manuscript; available in PMC 2010 October 26.

Published in final edited form as:

Procedia Comput Sci. 2010 May 1; 1(1): 1221–1229.

doi: 10.1016/j.procs.2010.04.136PMCID: PMC2964148

NIHMSID: NIHMS232482

See other articles in PMC that cite the published article.

The FFT EnKF data assimilation method is proposed and applied to a stochastic cell simulation of an epidemic, based on the S-I-R spread model. The FFT EnKF combines spatial statistics and ensemble filtering methodologies into a localized and computationally inexpensive version of EnKF with a very small ensemble, and it is further combined with the morphing EnKF to assimilate changes in the position of the epidemic.

Starting a model from initial conditions and then waiting for the result is rarely satisfactory. The model is generally incorrect, data is burdened with errors, and new data comes in that needs to be accounted for. This is a well-known problem in weather forecasting, and techniques to incorporate new data by sequential statistical estimation are known as data assimilation [1]. The ensemble Kalman filter (EnKF) [2] is a popular data assimilation method, which is easy to implement without any change in the model. The EnKF evolves an ensemble of simulations, and the model only needs to be capable of exporting its state and restarting from the state modified by the EnKF. However, the ensemble size required can be large (easily in the hundreds), the amount of computations in the EnKF can be significant, special localization techniques need to be employed to suppress spurious long-range correlations in the ensemble covariance matrix, and the EnKF does not work well for problems with sharp coherent features, such as the travelling waves found in epidemics and wildfires.

We propose a variant of EnKF based on the Fast Fourier transform (FFT), which reduces significantly the amount of computations required by the EnKF, as well as the ensemble size. The use of FFT is inspired by spatial statistic: FFT EnKF assumes that the state approximately a stationary random field, that is, the covariance between two points is mainly a function of their distance vector. Then the multiplication of the covariance matrix and a vector is a convolution. In addition, the morphing transform [3] is used here so that changes of the state both in position and in amplitude are possible.

The FFT EnKF with morphing is illustrated here for tracking a simulated epidemic wave. The use of data assimilation techniques can increase the accuracy and reliability of epidemic tracking by using the data as soon as they are available, and some applications of data assimilation in epidemiology already exist [4, 5]. The FFT EnKF with morphing has the potential to reduce complicated simulations and accurate real-time use of data to a laptop or a smartphone in the field.

For FFT EnKF in a wildfire simulation, see [6]. The Fourier domain Kalman filter (FDKF) [7] consists of the Kalman filter used in each Fourier mode separately.

The covariance of a stationary random field can be estimated from a single realization by the covariogram [8], which can be computed efficiently by the FFT [9]. We propose to use the covariogram for an *EnKF with an ensemble of one*, which will be further developed elsewhere.

The EnKF approximates the probability distribution of the model state *u* by an ensemble of simulations *u*_{1}, …, *u _{N}*. Each member is advanced by the simulation in time independently. When new data

$${u}_{k}^{a}={u}_{k}+{C}_{N}{H}^{T}{(H{C}_{N}{H}^{T}+R)}^{-1}(d+{e}_{k}-H{u}_{k}^{f}),\phantom{\rule{1em}{0ex}}k=1,\dots ,N,$$

(1)

to yield the *analysis ensemble* $\left[{u}_{k}^{a}\right]$. Here, *C _{N}* is an approximation of the covariance

When *C _{N}* is the ensemble covariance, the EnKF formulation (1) does not take advantage of any special structure of the model. This allows a simple and efficient implementation [12], but large ensembles, often over 100, are needed [2]. In an application, variables in the state are

The FFT EnKF discussed here uses a very small ensemble, but larger than one. We explain the FFT EnKF in the 1D case; higher-dimensional cases are exactly the same. Consider first the case when the model state consists of one block only. Denote by *u* (*x _{i}*),

$$\mathit{v}\left({x}_{i}\right)=\sum _{j=1}^{n}C({x}_{i},{x}_{j})u\left({x}_{j}\right)=\sum _{j=1}^{n}u\left({x}_{j}\right)c({x}_{i}-{x}_{j}),\phantom{\rule{1em}{0ex}}i=1,\dots ,n.$$

After FFT, convolution becomes entry-by-entry multiplication of vectors, that is, multiplication by a diagonal matrix.

We assume that the random field is approximately stationary, so we neglect the off-diagonal terms of the covariance matrix in the frequency domain, which leads to the the following FFT EnKF method. First apply FFT to each member, ${\widehat{u}}_{k}=F{u}_{k}$. Next, approximate the forecast covariance matrix in the frequency domain by the diagonal matrix with the diagonal entries given by

$$\widehat{{c}_{i}}=\frac{1}{N-1}\sum _{k=1}^{N}{\mid {\widehat{u}}_{ik}-{\overline{\widehat{u}}}_{i}\mid}^{2},\phantom{\rule{1em}{0ex}}\mathrm{where}\phantom{\rule{1em}{0ex}}{\overline{\widehat{u}}}_{i}=\frac{1}{N}\sum _{k=1}^{N}{\widehat{u}}_{ik}.$$

(2)

Then define approximate covariance matrix *C _{N}* by term-by-term multiplication · in the Fourier domain

$$u={C}_{N}\mathit{v}\iff \widehat{u}={F}_{u},\phantom{\rule{1em}{0ex}}\widehat{\mathit{v}}=\widehat{c}\u2022\widehat{u},\phantom{\rule{1em}{0ex}}\mathit{v}={F}^{-1}\widehat{\mathit{v}},\phantom{\rule{1em}{0ex}}{\left(\widehat{c}\u2022\widehat{u}\right)}_{i}=\widehat{{c}_{i}}\widehat{{u}_{i}}.$$

When *H* = *I* and *R* = *rI*, the evaluation of (1) reduces to

$${\widehat{u}}_{k}^{a}=\widehat{{u}_{k}}+\widehat{c}\u2022{(\widehat{c}+r)}^{-1}\u2022(\widehat{d}+\widehat{{e}_{k}}-{\widehat{u}}_{k}^{f}).$$

(3)

In general, the state has more than one variable, and *u*, *C*, and *H* have the block form

$$u=\left[\begin{array}{c}\hfill {u}^{\left(1\right)}\hfill \\ \hfill \vdots \hfill \\ \hfill {u}^{\left(n\right)}\hfill \end{array}\right],\phantom{\rule{1em}{0ex}}C=\left[\begin{array}{ccc}\hfill {C}^{\left(11\right)}\hfill & \hfill \cdots \hfill & \hfill {C}^{\left(1M\right)}\hfill \\ \hfill \vdots \hfill & \hfill \ddots \hfill & \hfill \vdots \hfill \\ \hfill {C}^{\left(M1\right)}\hfill & \hfill \cdots \hfill & \hfill {C}^{\left(MM\right)}\hfill \end{array}\right],\phantom{\rule{1em}{0ex}}H=[{H}^{\left(1\right)}\phantom{\rule{1em}{0ex}}\cdots \phantom{\rule{1em}{0ex}}{H}^{\left(M\right)}].$$

(4)

Here, the first variable is observed, so *H*^{(1)} = *I*, *H*^{(2)} = 0,…, *H*^{(M)} = 0, and (1) becomes

$${u}_{k}^{\left(j\right),a}={u}_{k}^{\left(j\right)}+{C}_{N}^{\left(j1\right)}{({C}_{N}^{\left(11\right)}+R)}^{-1}(d+{e}_{k}-{u}_{k}^{\left(1\right)}),\phantom{\rule{1em}{0ex}}j=1,\dots ,M,$$

(5)

and in the frequency domain

$${\widehat{u}}_{k}^{\left(j\right),a}={\widehat{u}}_{k}^{\left(j\right)}+{\widehat{c}}^{\left(j1\right)}\u2022{({\widehat{c}}^{\left(11\right)}+r)}^{-1}\u2022(\widehat{d}+\widehat{{e}_{k}}-\widehat{{u}_{k}}).$$

(6)

The cross-covariance between field *j* and field 1 is approximated by neglecting the off-diagonal terms of the sample covariance in the frequency domain as well,

$${\widehat{c}}_{i}^{\left(j1\right)}=\frac{1}{N-1}\sum _{k=1}^{N}\left({\widehat{u}}_{ik}^{\left(j\right)}-{\overline{\widehat{u}}}_{i}^{\left(j\right)}\right)\left({\widehat{u}}_{ik}^{\left(1\right)}-{\overline{\widehat{u}}}_{i}^{\left(1\right)}\right),\phantom{\rule{thickmathspace}{0ex}}\text{where}\phantom{\rule{thickmathspace}{0ex}}{\overline{\widehat{u}}}_{i}^{\left(\ell \right)}=\frac{1}{N}\sum _{k=1}^{N}{\widehat{u}}_{jk}^{\left(\ell \right)},\ell =1,j.$$

(7)

In the computations reported here, we have used the real sine transform, so all numbers in (7) are real. Also, the use of the sine transform naturally imposes no change of the state on the boundary.

Given an initial state *u*, the initial ensemble in the morphing EnKF [3, 12] is given by

$${u}_{k}^{\left(i\right)}=({u}_{N+1}^{\left(i\right)}+{r}_{k}^{\left(i\right)})\phantom{\rule{thinmathspace}{0ex}}o\phantom{\rule{thinmathspace}{0ex}}(I+{T}_{k}),\phantom{\rule{1em}{0ex}}k=1,\cdots ,N,$$

(8)

with an additional member *u*_{N+1} = *u*, called the *reference member*. In (8), ${r}_{k}^{\left(i\right)}$ are random smooth functions on Ω, *T _{k}* are random smooth mappings

The data *d* is an observation of *u*^{(1)}. The first blocks of all members *u*_{1},…, *u _{N}* and

$${u}_{k}^{\left(1\right)}\approx {u}_{N+1}^{\left(1\right)}\phantom{\rule{thinmathspace}{0ex}}o\phantom{\rule{thinmathspace}{0ex}}(I+{T}_{k}),\phantom{\rule{1em}{0ex}}{T}_{k}\approx 0,\phantom{\rule{1em}{0ex}}\nabla {T}_{k}\approx 0,\phantom{\rule{1em}{0ex}}k=0,\dots ,N,$$

${u}_{0}^{\left(1\right)}=d$ and *T _{k}* : Ω → Ω

$${u}_{k}\mapsto {\stackrel{~}{u}}_{k}={M}_{{u}_{N+1}}\phantom{\rule{thinmathspace}{0ex}}\left({u}_{k}\right)=({T}_{k},{r}_{k}^{\left(1\right)},\dots ,{r}_{k}^{\left(M\right)}),$$

(9)

whare ${r}_{k}^{\left(j\right)}={u}_{k}^{\left(j\right)}o{(I+{T}_{k})}^{-1}-{u}_{N+1}^{\left(j\right)}$, *k* = 0,…, *N*, are *registration residuals*. Likewise, the extended data vector is defined by $d\mapsto \stackrel{~}{d}=({T}_{0},{r}_{0}^{\left(1\right)})$ and the the observation operator is (*T, r*^{(1)},…, *r ^{(M)}*) (

$${u}_{k}^{a,\left(i\right)}={M}_{{u}_{N+1}}^{-1}\phantom{\rule{thinmathspace}{0ex}}\left({\stackrel{~}{u}}_{k}^{a}\right)=({u}_{N+1}^{\left(i\right)}+{r}_{k}^{a,\left(i\right)})\phantom{\rule{thinmathspace}{0ex}}o\phantom{\rule{thinmathspace}{0ex}}(I+{T}_{k}^{a}),\phantom{\rule{1em}{0ex}}k=1,\dots ,N+1,$$

(10)

where the new transformed reference member is given by

$${\stackrel{~}{u}}_{N+1}^{a}=\frac{1}{N}\sum _{k=1}^{N}{\stackrel{~}{u}}_{k}^{a}.$$

(11)

The epidemic model that we used for this study is a spatial version of the common S-I-R dynamic epidemic model. A person is *susceptible* or *infectious* in this context if he or she can contract or transmit the disease, respectively. The *removed* state includes those who have either died, have been quarantined, or have recovered from the disease and become immune. The state variables are the susceptible (*S*), the infectious (*I*), and the removed (*R*) population densities. The core ideas for this model date back to the 1957 spatial formulation by Bailey [14], but the specific version that we have employed here is due to Hoppenstaedt [15, p. 64].

The population is considered to be dispersed over a planar domain $\Omega \subset {\mathbb{R}}^{2}$, and it is labelled according to its position with respect to the spatial coordinates *x* and *y*. The (deterministic) evolution of the state (*S* (*t*), *I* (*t*), *R* (*t*)) is given by

$$\phantom{\{}\begin{array}{c}\frac{\partial S(x,y,t)}{\partial t}=-S(x,y,t)\underset{\Omega}{\iint}w(x,y,u,v)I(u,v,t)dudv,\hfill \\ \frac{\partial I(x,y,t)}{\partial t}=S(x,y,t)\underset{\Omega}{\iint}w(x,y,u,v)I(u,v,t)dudv-q(x,y,t)I(x,y,t),\hfill \\ \frac{\partial R(x,y,t)}{\partial t}={q}_{i}(x,y,t)I(x,y,t).\hfill \end{array}\}$$

(12)

The function *q* (*x*, *y*, *t*) gives the rate of removal of infectives due to death, quarantine, or recovery. The weight function *w* (*x*, *y*, *u*, *v*) measures the influence of infectives at spatial position (*u*, *v*) on the exposure of susceptibles at position (*x*, *y*); in this simulation we used the function *w* (*x*, *y*, *u*, *v*) = α exp [−((*x* − *u*)^{2} + (*y* − *v*)^{2})^{1/2}/λ], which expresses the idea that the influence of nearby infectives decays as an exponential function of Euclidean distance, with constant, λ characteristic of the distance at which the disease spreads. More mobile societies will have larger values of λ. The parameter α measures the infectiousness of the disease.

A stochastic cell model is created by treating the quantities on the right-hand-side of (12) as the intensities of a Poisson process and by piecewise constant integration over the cells. The domain Ω is decomposed into nonover-lapping cells Ω_{i} with centers (*x _{i}*,

$${S}_{i}\phantom{\rule{thinmathspace}{0ex}}(t+\Delta t)={S}_{i}\phantom{\rule{thinmathspace}{0ex}}\left(t\right)-\Delta {S}_{i},\phantom{\rule{1em}{0ex}}{I}_{i}\phantom{\rule{thinmathspace}{0ex}}(t+\Delta t)={I}_{i}\phantom{\rule{thinmathspace}{0ex}}\left(t\right)+\Delta {S}_{i}-\Delta {R}_{i},\phantom{\rule{1em}{0ex}}{R}_{i}\phantom{\rule{thinmathspace}{0ex}}(t+\Delta t)={R}_{i}\phantom{\rule{thinmathspace}{0ex}}\left(t\right)+\Delta {R}_{i},$$

where the random increments Δ*S _{i}* and Δ

$$\begin{array}{cc}\hfill & \Delta {S}_{i}~\mathrm{Pois}\left({S}_{i}\phantom{\rule{thinmathspace}{0ex}}\left(t\right){\Sigma}_{j=1}^{K}w({x}_{i},{y}_{i},{x}_{j},{y}_{j}){I}_{j}\phantom{\rule{thinmathspace}{0ex}}\left(t\right)\phantom{\rule{thinmathspace}{0ex}}A\phantom{\rule{thinmathspace}{0ex}}\left({\Omega}_{i}\right)\phantom{\rule{thinmathspace}{0ex}}\Delta t\right),\hfill \\ \hfill & \Delta {R}_{i}~\mathrm{Pois}\left({q}_{i}\phantom{\rule{thinmathspace}{0ex}}\left(t\right)\phantom{\rule{thinmathspace}{0ex}}{I}_{i}\phantom{\rule{thinmathspace}{0ex}}\left(t\right)\phantom{\rule{thinmathspace}{0ex}}A\phantom{\rule{thinmathspace}{0ex}}\left({\Omega}_{i}\right)\phantom{\rule{thinmathspace}{0ex}}\Delta t\right),\hfill \end{array}$$

(13)

and *q _{i}* (

We have chosen to model an epidemic disease that first emerges in Congo. The computational domain is a square portion of central Africa. In Figure 1 (a), we see the epidemic wave 120 model time steps after the emergence of the disease. The behavior of the model is such that any spurious infection will tend to grow into a secondary infection wave. This is problematic for data assimilation because the occurrence of spurious features is virtually guaranteed. We attempt to reduce the occurrence and magnitude of these features using the morphing transformation and FFT EnKF; however, some amount of residual artifacts will remain. We have found that by processing the model state in the following manner, we can further reduce these artifacts. We begin by scaling the absolute quantities contained in the model variables to a percentage of the local population before performing the data assimilation. After data assimilation, we truncate the variables to the range [0, 1], and we apply a threshold so that any infection rate below 1% is set to 0. Finally, we rescale the output in absolute units ensuring that the number of people at each grid cell is preserved. We have applied the FFT EnKF to the epidemic model described in Section 4 with an ensemble of size 5. Each ensemble simulation was started with the same initial conditions, but with different random seeds, and advanced in time by 100 model time units, then perturbed randomly to obtain the initial ensemble. The analysis ensemble and data were advanced in time an additional 20 model time steps for further assimilation cycles. In total, 3 assimilation cycles were performed in this manner.

(a) The number of people per kilometer squared infected, susceptible, and removed after 120 time steps in a simulation of an epidemic disease spreading through central Africa. These images correspond to variables *I*, *S*, and *R* in Equation (12). (b) Number **...**

We have perturbed each member of the initial ensemble randomly in space by applying (10) to the each variable of the morphing representation of the model. The mappings *T _{k}* for this perturbation were generated from a space of smooth functions that are zero at the boundary. While the residuals

The output of the observation function used in this example consists of the *Infected* field of the model. In this case, the data is a spatial “image” of the number of infected persons in each grid cell. The data were generated synthetically from a model simulation, which was initialized in the same manner as the ensemble.

Four variants of the EnKF were then applied: the standard EnKF and FFT EnKF and morphing EnKF and morphing FFT EnKF. The same initial ensemble and the same data were used for each method. The deviation of the initial ensemble and the model error were chosen so that the analysis should be about half way between the forecast and the data. In the morphing variants, the data deviation in the amplitude was taken very large, so that the filter updates essentially only the position. Ensemble of size 5 was used. The result in the first assimilation cycle for each method is shown in Figures 2 and and3.3. The first image in each column is the forecast mean. In the morphing variants, the mean is taken over all ensemble members in all fields of the morphing representation (9) and it plays the role of the comparison state for registration. Thus, in the morphing variants, both the amplitude and the position of the infection wave in the ensemble members are averaged. The second image in each column is the data, which is a model trajectory started from the same initial state for each method. Because the model is itself stochastic, the data images are slightly different. The third image in each column is the analysis mean, which is taken in the morphing representation (11) for two morphing filters, so that both the amplitude and the location are averaged.

The number of people infected per kilometer squared in analysis cycle 1 using the standard EnKF and FFT EnKF, each with ensemble size of 5. Both approaches are unable to move the location of the infection in the simulation state.

The number of people infected per kilometer squared in analysis cycle 1 using the morphing EnKF and morphing FFT EnKF, each with ensemble size of 5. Both approaches are able to move the state spatially and perform similarly. However, EnKF suffers from **...**

We see that both standard EnKF and FFT EnKF filters cannot move the state towards the data; a much larger ensemble would be needed. The morphing EnKF does move the state towards the data, but there are strong artifacts due to the poor approximation of the covariance by the covariance of the small ensemble. Finally, the morphing FFT-EnKF is capable of moving the state towards the data better.

We have introduced morphing FFT EnKF and presented preliminary results on data assimilation for an epidemic simulation. Morphing was essential to move the state towards the data, but it resulted in artifacts for the small ensemble size used, yet small ensemble size is important to perform simulations with data assimilation on general computing devices instead of supercomputers. We have observed that the estimation of the covariance matrix in the frequency domain results in better forecast covariance in the algorithm, which has the potential to reduce the artifacts due to small ensemble size.

This work was partially supported by NIH grant 1 RC1 LM01641-01 and NSF grants CNS-0719641 and ATM-0835579.

[1] Kalnay E. Atmospheric Modeling, Data Assimilation and Predictability. Cambridge University Press; 2003.

[2] Evensen G. Data Assimilation: The Ensemble Kalman Filter. 2nd Edition Springer Verlag; 2009.

[3] Beezley JD, Mandel J. Morphing ensemble Kalman filters. Tellus. 2008;60A:131–140.

[4] Bettencourt L, Ribeiro R, Chowell G, Lant T, Castillo-Chavez C. Intelligence and Security Informatics: Biosurveillance, Vol. 4506 of Lecture Notes in Computer Science. Springer; 2007. Towards real time epidemiology: data assimilation, modeling and anomaly detection of health surveillance data streams; pp. 79–90.

[5] Rhodes C, Hollingsworth T. Variational data assimilation with epidemic models. Journal of Theoretical Biology. 2009;258(4):591–602. [PubMed]

[6] Mandel J, Beezley JD, Kondratenko VY. Fast Fourier transform ensemble Kalman filter with application to a coupled atmosphere-wildland fire model, arXiv:1001.1588. International Conference on Modeling and Simulation (MS'2010); 2010. accepted.

[7] Castronovo E, Harlim J, Majda AJ. Mathematical test criteria for filtering complex systems: plentiful observations. J. Comput. Phys. 2008;227(7):3678–3714.

[8] Cressie NAC. Statistics for Spatial Data. John Wiley & Sons Inc.; New York: 1993.

[9] Marcotte D. Fast variogram computation with FFT. Computers & Geosciences. 1996;22(10):1175–1186.

[10] Burgers G, van Leeuwen PJ, Evensen G. Analysis scheme in the ensemble Kalman filter. Monthly Weather Review. 1998;126:1719–1724.

[11] Mandel J, Cobb L, Beezley JD. On the convergence of the ensemble Kalman filter, arXiv:0901.2951. Applications of Mathematics. 2009 January; to appear.

[12] Mandel J, Beezley JD, Coen JL, Kim M. Data assimilation for wildland fires: Ensemble Kalman filters in coupled atmosphere-surface models. IEEE Control Systems Magazine. 2009;29:47–65.

[13] Furrer R, Bengtsson T. Estimation of high-dimensional prior and posterior covariance matrices in Kalman filter variants. J. Multivariate Anal. 2007;98(2):227–255.

[14] Bailey N. Mathematical Theory of Epidemics. Griffin; 1957.

[15] Hoppenstaedt F. Mathematical Theories of Populations, Demographics, and Epidemics, CBMS-NSF Regional Conference Series in Applied Mathematics. SIAM. 1975

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |