A coherent structure approach for parameter estimation in Lagrangian Data Assimilation

TL;DR

This paper presents a novel data assimilation method that estimates model parameters by directly assimilating Lagrangian Coherent Structures using ABC and SMC, outperforming traditional methods in chaotic advection scenarios.

Contribution

The paper introduces a new approach combining Lagrangian Coherent Structures, ABC, and PCA for parameter estimation, avoiding likelihood computation and improving accuracy in chaotic flows.

Findings

Significantly better parameter estimation results than bootstrap particle filter.

Effective identification of coherent patterns using PCA from tracer data.

Robust performance in chaotic advection conditions.

Abstract

We introduce a data assimilation method to estimate model parameters with observations of passive tracers by directly assimilating Lagrangian Coherent Structures. Our approach differs from the usual Lagrangian Data Assimilation approach, where parameters are estimated based on tracer trajectories. We employ the Approximate Bayesian Computation (ABC) framework to avoid computing the likelihood function of the coherent structure, which is usually unavailable. We solve the ABC by a Sequential Monte Carlo (SMC) method, and use Principal Component Analysis (PCA) to identify the coherent patterns from tracer trajectory data. Our new method shows remarkably improved results compared to the bootstrap particle filter when the physical model exhibits chaotic advection.

Equations

$${\mathbf{P}}_{x_1}=\frac{1}{d-1}({\mathbf{A}}_{x_1}{\mathbf{A}}_{x_1}^{\prime }).$$

$$H(f,g)=\frac{1}{\sqrt{2}}\Vert \sqrt{f}-\sqrt{g}{\Vert }_2.$$

$$\pi (\theta |y)\propto f(y|\theta )\pi (\theta ),$$

$$f(y|\theta )=\int \tilde{f}(y|\theta ,{\mathbf{x}}_{0:n})\pi ({\mathbf{x}}_{0:n}|\theta )d{\mathbf{x}}_{0:n},$$

$$f(y|\theta )\approx \tilde{f}(y|{\hat{\mathbf{x}}}_{0:n}(\theta ))\equiv \tilde{f}(y|\theta ,{\hat{\mathbf{x}}}_{0:n}).$$

$$\pi (\theta |y)\propto \tilde{f}(y|{\hat{\mathbf{x}}}_{0:n}(\theta ))\pi (\theta ).$$

$${\pi }_{\epsilon }(\theta ,y|y^o)\propto f(y|\theta )\pi (\theta )I_{\epsilon }(y),$$

$$I_{\epsilon }(y)=\{\begin{array}{l}1\text{for }\rho (\eta (y),\eta (y^o))<\epsilon ,\\ 0\text{otherwise}.\end{array}$$

$${\pi }_{\epsilon }(\theta |y^o)=\int {\pi }_{\epsilon }(\theta ,y|y^o)dy\approx \pi (\theta |y^o)$$

$${\tilde{\pi }}_n(z_{0:n})={\pi }_n(z_n)\prod _{j=0}^{n-1}{L}_j(z_{j+1},z_j),$$

$$w_n^{(i)}\propto w_{n-1}^{(i)}\frac{{\pi }_n(z_n^{(i)})L_{n-1}(z_n^{(i)},z_{n-1}^{(i)})}{{\pi }_{n-1}(z_{n-1}^{(i)})K_n(z_{n-1}^{(i)},z_n^{(i)})}.$$

$$N_{eff}\left(\{w_{n}^{(i)})\}\right):=\biggl{(}\sum_{i=1}^{N}(w_{n}^{(i)})^{2}\biggr{)}^{-1}.$$

$$L_{n-1}(z,z^{\prime })=\frac{{\pi }_n(z^{\prime })K_n(z^{\prime },z)}{{\pi }_n(z)}.$$

$${\tilde{w}}_n^{(i)}\propto \frac{{\pi }_{{ϵ}_n}({\theta }_{n-1}^{(i)},y_{n-1}^{(i)}|y^o)}{{\pi }_{{ϵ}_{n-1}}({\theta }_{n-1}^{(i)},y_{n-1}^{(i)}|y^o)}{w}_{n-1}^{(i)}$$

$${\tilde{w}}_n^{(i)}=\{\begin{array}{l}{w}_{n-1}^{(i)}{I}_{{ϵ}_n}(y_{n-1}^{(i)})\text{for }{w}_{n-1}^{(i)}>0\\ 0\text{otherwise}\end{array},$$

$$\#\{w_n^{(i)}:w_n^{(i)}>0\}=\alpha \#\{w_{n-1}^{(i)}:w_{n-1}^{(i)}>0\},$$

$$1\wedge \frac{I_{{ϵ}_n}(y^{\ast })}{I_{{ϵ}_{n-1}}(y)}\frac{q_n({\theta }^{\ast },\theta )}{q_n(\theta ,{\theta }^{\ast })}\frac{\pi ({\theta }^{\ast })}{\pi (\theta )}.$$

$$w_n^{(i)}\propto$$

$$p\left({\mathbf{x}}_n^o|{\mathbf{x}}_n^{(i)},{\theta }_n^{(i)}\right)\propto \exp \left(-\frac{1}{2}({\mathbf{x}}_n^o-{\mathbf{x}}_n^{(i)}{)}^T{\mathbf{R}}^{-1}({\mathbf{x}}_n^o-{\mathbf{x}}_n^{(i)})\right),$$

$$d{x}_1$$

$$d{x}_2$$

$${\theta }^{\ast }=\theta +\eta ,\eta \sim \mathcal{N}(0,0.01).$$

Full text

1 Introduction

The problem of assimilating the path travelled by ocean instruments such as drifters to estimate states and/or parameters of dynamical systems is commonly known as Lagrangian data assimilation (LaDA). Every method for LaDA depends on the knowledge of error statistics for the observed position of the drifter paths; typically, it assumes the normal likelihood function and uncorrelated errors of observed position among all drifters. In this paper, we explore a novel idea for LaDA that assimilates a “coherent (or persistent) structure” hidden in the Lagrangian path of the drifters instead of directly using the observed positions. We discuss the main advantage of this idea in the situation where the number of drifters is large, observed position has small variance and the flow is chaotic. We also demonstrate the improved accuracy of the parameter estimates of the new method when comparing with the particle filtering (PF) approach, which is apparent if the drifter paths are dominated by chaotic advection.

Modern satellite-tracked ocean drifters [13, 7] have been recognized as an important source of data for oceanographic and climate research. Drifters travel near the surface and measure physical data such as surface temperature along their paths via drifter sensors [12]. In addition, they offer Lagrangian path data through the Doppler frequency shift on their satellite-based transmission. The accuracy of drifter locations can be less than 150 meters [7], and there is a global ocean drifter set providing measurement data at a mean time interval of 1.2 hours spanning over a decade [13].

Oceanographic flows exhibit patterns at scales from metres to kilometres, and correspondingly drifters are seen to move in eddies or gyres and follow currents, creating small and large-scale patterns which can be used to infer the sea-surface flows [38, 14]. These patterns are associated with fundamental structures of the flow, including stable and unstable directions, centres and saddles. The goal of the current work is to infer the underlying flow state, or the parameter set governing the flow, from the patterns rather than directly from the drifter trajectories.

Several LaDA methods have been presented in the last decade for sequentially assimilating Lagrangian path data into dynamical system models of the ocean [25, 31, 40, 49, 2, 29, 27, 43, 44]. These methods append a Lagrangian model of the drifter locations to the Eulerian flow model in order to make inferences about how drifter locations correlate with model states; this approach also avoids the difficulty in which Lagrangian path data does not correspond to the fixed Eulerian grid locations. The key challenge of Lagrangian Data Assimilation (LaDA) is that the Lagrangian path of the passive tracers or drifters can be strongly nonlinear, which makes it difficult to sequentially estimate model state variables or parameters by computationally efficient methods such as the Ensemble Kalman filter (EnKF) [2]. It has been emphasized in several pieces of work that a nonlinear, Bayesian filtering framework such as the Particle Filter (PF) is more reliable in quantifying the uncertainty of the estimates [2, 49, 29, 1, 43]. However, when the drifter locations are observed at a high accuracy (e.g. the class-three Argos), the observation is considered to be highly informative, corresponding to a likelihood probability with a small variance in the Bayesian framework. Assimilating highly-informative observations in a high-dimensional space is problematic since significant probability is contained in a region of extremely small “volume”. The standard PF tends to perform very poorly in such a regime [34]. With a narrow prior, it is probable that no particles will lie in a region of significant probability after a short amount of time. In this situation the PF weights, which fall off exponentially, assign probability one to a single particle and probability 0 to the rest. Therefore, it would be difficult to apply PF or its variants to high-dimensional models of the ocean, or assimilate data from a large number of drifters, such as the approximately 2200 ocean drifters from the Global Ocean drifter Program (www.aoml.noaa.gov/envids/gld/). The problem of drawing inference from drifter trajectories is compounded by factors which tend to separate the numerically simulated drifters from the true trajectories. One such factor is model error, which is not considered in this paper. A second factor of particular importance to the current work is chaotic advection in the flow, which can readily create a situation where the simulated and observed drifter trajectories will appear unrelated, particularly if the observations are sparse.

A hybrid PF-EnKF method has been proposed to reduce the computational cost in high-dimensional problems [47, 37]; in this approach it is assumed that the Eulerian flow model is high-dimensional and the Lagrangian model for the drifter locations is low-dimensional. A similar idea using a hybrid PF-EnKF to assimilate the Lagrangian path data has also been developed in [43, 44] where EnKF is used to estimate drifter paths in tandem with PF for parameter estimation. However, to the best of our knowledge, the challenge of assimilating a high-dimensional drifter data set is still an open problem.

These issues motivate us to exploit the qualitative structure of complex nonlinear flows - such as coherent structures or persistent patterns, if they exist - to make inferences on the model parameters.

In particular, we will exploit a large-scale “Lagrangian coherent pattern” or just “coherent pattern” hidden in the Lagrangian path data. The coherent pattern is imprecisely defined as a region in state space, for example a coherent vortex or nonlinear jet, that moves along with the flow without dispersing. Coherent patterns must also be robust under small diffusive peturbations, that is, the coherent object should still hold its geometric structure together under some diffusion process. When the flow is known either numerically or in a closed form, these coherent regions can be identified by several approaches. In the probabilistic approach, the almost invariant-set framework was developed for autonomous and periodic flows by using the transfer (or Perron-Frobenius) operator [10, 11, 5, 15] or the infinitesimal generator [17]; the latter reduces the high computational cost of the former approach. For nonautonomous flows, the finite-time coherent set method via a transfer operator was introduced for the first time in [18, 42, 21] and it was applied to detect long-lived vortices such as the stratospheric polar vortex [42, 21] and Agulhas rings [16], both as two and three dimensional coherent objects. This method has a strong connection with spectral clustering as described in Chapter of 4 of [6]. The coherent region as used in these works probabilistically minimizes the mass transport in and out of the region with respect to a reference measure (not necessarily the invariant measure), see [42, 19] for full details.

Unfortunately, all of the above coherent set identification methods require a set of governing equations of the underlying dynamics and cannot be applied to identify the coherent pattern directly from the Lagrangian path data, which is an essential goal of using LaDA in the current work. Recently, new data-driven algorithms have been developed that can extract coherent patterns from possibly sparse Lagrangian path data without having to rely on governing equations. To name a few, these methods include the diffusion-map algorithm [3], fuzzy c-mean algoritm[20], bipartite spectral clustering [23] and dynamic mode decomposition (DMD) [45, 52]. These methods can capture the time-dependent coherent pattern that represents a slowly-decaying mode of a complex flow field.

It is not a primary concern of this paper to select or consider the ideal algorithms to identify the spatio-temporal coherent patterns in a general context, but to advocate for a methodology in which identified coherent patterns are assimilated in a Bayesian framework. For this reason, numerical examples in this paper will be limited to the case of a stationary spatial pattern where an application of standard principal component analysis (PCA) is satisfactory for the coherent pattern identification. We will provide a brief review of PCA in Section 2.

In Bayesian data assimilation, a closed form of the likelihood function for the coherent pattern must be known. However, even if the likelihood of tracer positions is commonly known, it is still a difficult task to derive the likelihood function of any corresponding coherent patterns. We will address this issue by applying the so-called “Approximate” Bayesian computation (ABC) [30], which can be thought of as a “likelihood-free” Monte Carlo method. Originally, the ABC method was proposed to avoid evaluating a likelihood function that is computationally expensive by constructing a “distance function” together with a summary statistic in a rejection-acceptance algorithm. In our context, this summary statistic is the coherent pattern. It is typically impossible to construct the required distance function without knowledge of the likelihood function of the coherent pattern; we will argue that given the likelihood function is unknown, a distance function with minima in the most likely regions of parameter values is a useful surrogate. We will discuss our choice of the distance function in Section 3 and ABC methods in Section 5.

2 Coherent spatial patterns via PCA

As we will make inference about model parameters based on the coherent pattern of the flow, we use the dominant eigenvector (corresponding to the eigenvalue with the largest magnitude) of the covariance matrix as a representative of the coherent pattern, as conventionally performed by the Principal Component Analysis (PCA) [26]. Note that we focus on assimilating only the most “dominant” pattern in this work, though several patterns may exist for the underlying flow. A given trajectory of $d$ tracers, consisting of recordings of the $x_{1}$ and $x_{2}$ coordinates at times $t_{0}<t_{1}<t_{2}<\cdots<t_{n}$, is denoted by $\mathbf{x}_{0:n}:=\{\mathbf{x}_{0:n}^{(i)}\}_{i=1}^{d}$, where $\mathbf{x}_{0:n}^{(i)}=\{(x_{1}^{(i)}(t_{0}),x_{2}^{(i)}(t_{0})),\ldots,(x_{1}^{(i)}(t_{n}),x_{2}^{(i)}(t_{n}))\}$. Let $x_{1,0:n}^{(i)}=[x_{1}^{(i)}(t_{0}),\ldots,x_{1}^{(i)}(t_{n})]$ for $i=1,\ldots,d$ and $\tilde{\mathbf{A}}_{x_{1}}$ be a $d\times n$ matrix whose $i-$th row is $x_{1,0:n}^{(i)}$. We define the anomaly matrix of the $x_{1}$ coordinate by $\mathbf{A}_{x_{1}}=\tilde{\mathbf{A}}_{x_{1}}-\bar{\mathbf{x}}_{1}\otimes\mathbf{1}_{d}$, where the vector $\bar{\mathbf{x}}_{1}$ denotes the average of the tracer position over their trajectory and $\mathbf{1}_{d}$ is a $d$-dimensional vector with all entries 1. The covariance matrix can then be constructed by

[TABLE]

The dominant eigenvector of $\mathbf{P}_{x_{1}}$ gives an approximation of the coherent pattern according to the correlation structure of the $x_{1}$-coordinate of the tracer trajectories. In general, one may have to consider different patterns that may be captured by eigenvectors with smaller eigenvalues. The coherent pattern described by the correlation structure in the $x_{2}-$direction can be similarly extracted from $\mathbf{P}_{x_{2}}$, which can be constructed in a similar manner. For the numerical experiment in this work, it will be adequate to use only $\mathbf{P}_{x_{1}}$ to identify coherent patterns. In more general cases, recent advanced techniques [3, 20, 23] to identify coherent patterns from the data may be required. We remark that the time scale is important to observe the coherent patterns. Thus, $t_{n}$ has to be “long enough”, which will depend on the underlying dynamics.

3 Hellinger distance

In order to implement ABC, we need a distance function. We adopt the Hellinger distance to model the “distance” between the predicted spatial pattern $f\in\mathbb{R}^{d}$, and the observed pattern $g\in\mathbb{R}^{d}$, both of which are identified by performing PCA on the tracer trajectory data. We normalize $f$ and $g$ so that their entries are between 0 and 1 and summed to unity. The Hellinger distance, see [28], between $f$ and $g$ is given by

[TABLE]

Therefore, $H(f,g)$ satisfies the triangle inequality and $0\leq H(f,g)\leq 1$ for all $f$ and $g$ normalized as above. When comparing the patterns, the difference in the observed positions of the tracers and the predicted positions is irrelevant. What is important is that we are able to group together the tracers with similar dynamical fates; for example, tracers in a persistent vortex should be approximately clustered as a group, separating them out from the tracers in a strongly mixing or turbulent region. The preceding remark is especially crucial in assimilating a large number of tracers for noisy dynamical systems or chaotic flows because the model prediction of tracer positions can be inaccurate due to the noisy peturbation in the flow, while the coherent pattern is typically robust under peturbation [24].

4 Likelihood function of coherent structures

Although considerable research has been devoted to the augmented state approach for LaDA (where model parameters are augmented to the state), rather less attention has been paid to the situation where the tracer positions are not of interest as their positions are observed quite accurately; in this situation the tracer positions are “nuisance” random variables in the inference problem.

Suppose that the random variable $\theta$ aggregates all unknown parameters, and a random variable $y$ represents the coherent pattern identified by the principle component. At the heart of our proposed LaDA method for parameter estimation is the inference problem

[TABLE]

where $f(y|\theta)$ is the likelihood function and $\pi(\theta)$ is the prior distribution. Due to the uncertainty of the tracer locations $\mathbf{x}_{0:n}$, the likelihood function of $y$ conditional on $\theta$ is obtained by

[TABLE]

where $\tilde{f}(y|\theta,\mathbf{x}_{0:n})$ is the likelihood function of $y$ conditional on $\theta$ and $\mathbf{x}_{0:n}$. The likelihood function $\tilde{f}(y|\theta,\mathbf{x}_{0:n})$ is usually intractable due to a lack of knowledge of its analytic form, and so is (4.2). One way to simplify the problem is to substitute the random variable $\mathbf{x}_{0:n}$ in (4.2) with some reference value $\hat{\mathbf{x}}_{0:n}(\theta)$ (e.g. the mean, a sample, etc.). By doing so, we have

[TABLE]

As a consequence, the inference problem (4.1) can be approximated by

[TABLE]

We emphasize that the approximated likelihood $f(y|\theta)$ in(4.3) still lacks an analytic form; hence, the evaluation of $f(y|\theta=\theta^{\ast})$ for a given parameter value $\theta^{\ast}$ is still unavailable. We can, however, simulate the coherent pattern $y$ for $\theta^{\ast}$, which allows us to find the Hellinger distance between the predicted pattern (extracted from a sample of $\mathbf{x}_{0:n}$) and the observed pattern $y$ (extracted from the observed tracer trajectories).

Although the statistical approximation properties of (4.4) are not examined herein, we heuristically argue that the approximation (4.3) is reasonable in the current application. As pointed out earlier, the coherent pattern is robust to stochastic peturbation of the velocity field, so most samples of Lagrangian path $\mathbf{x}_{0:n}$ would remain within the same flow regime. For example, even under the stochastic peturbation, most samples of the tracers in the recirculation regions tend to be in the recirculation region and those in the jet stream would remain in the jet stream as long as they are propagated under the same model parameters. It is the model parameters that prominently determine the large-scale coherent patterns of the flow. Thus, the likelihood $\tilde{f}(y|\theta=\theta^{\ast},\mathbf{x}_{0:n})$ for a fixed value $\theta^{\ast}$ would have almost the same value for any sample of $\mathbf{x}_{0:n}$, so the approximation (4.3) should be acceptable. In the next section, we describe a framework to solve the Bayesian computation of (4.4) that only uses a distance function to avoid the difficulty of evaluating the likelihood function.

5 Approximate Bayesian computation (ABC)

Given a prior density $\pi(\theta)$ of the parameter $\theta\in\Theta$, consider a situation where a likelihood function $f(y|\theta)$ of observations $y\in\mathcal{D}$ is not available in a closed form or its evaluation is computationally intractable. It is then difficult to use standard Bayesian computation methods to sample the posterior density $\pi(\theta|y)$. Approximate Bayesian computation (ABC) is a likelihood-free monte carlo approach that only requires a simulation (or sampling) from the likelihood $f(y|\theta)$ but avoids an evaluation of $f(y|\theta)$. The foundation of ABC algorithms was suggested by Rubin as early as 1984 [39]: that if we observe a discrete random variable $y^{o}\sim f(y|\theta)$, then we can sample $\pi(\theta|y)$ by jointly sampling $\theta^{\ast}\sim\pi(\theta)$ and $y^{\ast}\sim f(y|\theta^{\ast})$, accepting $\theta^{\ast}$ and the pseudo observation $y^{\ast}$ only if $y^{\ast}=y^{o}$. The accepted values of $\theta^{\ast}$ will distribute according to the posterior density.

In the ABC framework, the acceptance-rejection step replaces the impractical condition $y^{\ast}=y^{o}$ with $\rho(\eta(y^{\ast}),\eta(y^{o}))<\varepsilon$, where $\rho$ is some measure of discrepancy, $\varepsilon$ is a threshold level, and $\eta$ is some summary statistic. The threshold $\varepsilon$ controls the trade-off between computational cost and accuracy; for a small $\varepsilon$, the sample will give a good representation of the true posterior distribution but require more computation. Thus, all ABC algorithms draw the sample from the approximate posterior density

[TABLE]

where

[TABLE]

Provided that $\varepsilon$ is sufficiently small and $\eta$ is a sufficient statistic, we hope that

[TABLE]

The simplest ABC algorithm draws a sample of $\theta$ by the following acceptance-rejection algorithm [50, 46]:

Algorithm 0

Step 1.

Sample $\theta\sim\pi(\theta)$

Step 2.

Sample $y\sim f(y|\theta)$

Step 3.

Accept $\theta$ if $\rho(\eta(y),\eta(y^{o}))\leq\varepsilon$

In the current application, $\theta$ is the model parameters, $y^{o}$ is the observed dominant spatial pattern, $y$ is the predicted spatial pattern, $\rho(y,y^{o})$ is a distance function designed by the user, which is the Hellinger distance in this work, and $\eta$ is the identity function (i.e. $\eta(y)=y$). In the second step of the algorithm, we use the PCA method for pattern identification. Unfortunately, the above algorithm is ineffective due to a low acceptance rate, especially in the case where the prior and posterior distributions are very different. To improve the algorithm efficiency, several algorithms have been recently proposed, e.g., the adaptive SMC-ABC [4, 51, 33], which incorporates a sequential Monte Carlo (SMC) sampler [32] to approximate the posterior density by a collection of weighted “particles” $\{\theta^{(i)},w^{(i)}\}$ for $i=1,\ldots,N$. We will only outline the SMC-ABC algorithm used in our work and refer the details and discussion of the SMC-ABC in general to [32, 33].

5.1 SMC-ABC

Suppose that we wish to approximate a probability density $\pi(z)$ by random samples, where in the current context $z=(\theta,y)$, the collection of parameter estimates together with the corresponding predicted coherent patterns. The SMC sampler is an approach to approximate a sequence of probability distributions $\tilde{\pi}_{n}(z_{0:n})$ for $n=1,\ldots,T$ and $z_{0:n}\equiv(z_{0},\ldots,z_{n})$ that admits $\pi_{n}(z_{n})$ as a marginal distribution and satisfies $\pi_{T}(z_{T})\approx\pi(z)$. Del Moral et al. [33] proposed the following sequence

[TABLE]

where $L_{n}(z_{n+1},z_{n})$ is the probability density of moving a given $z_{n+1}$ backward to $z_{n}$. Thus, given weighted samples $\{z_{n-1}^{(i)},w_{n-1}^{(i)}\}$ at the time step $n-1$, the SMC sampler draws $z_{n}^{(i)}\sim K_{n}(z_{n-1}^{(i)},\cdot)$ via the (forward) Markov transition kernel $K_{n}$ and updates the weight according to

[TABLE]

The above weights are then normalized so that $\sum_{i=1}^{N}w_{n}^{(i)}=1$. A resampling step may be required to prevent the so-called particle collapse in which most particles have negligible weights. The degree of particle collapse is typically monitored by the effective sample size (ESS) given by

[TABLE]

The resampling step is triggered when $N_{eff}<\delta N$ for some $0<\delta<1$.

Again, we consider $z=(\theta,y)$ and the target distribution is $\pi_{\epsilon}(\theta,y|y^{o})$. We are interested in sampling a sequence $\pi_{\epsilon_{n}}(\theta,y|y^{o})$ for $\epsilon_{0}>\epsilon_{1}>\cdots>\epsilon_{T}=\epsilon$ and the required marginal $\pi_{\epsilon}(\theta|y^{o})$ is obtained by the construction in (5.4). The optimal choice of the backward kernel $L_{n}$ is given in [33] but it is typically difficult to implement. A more convenient choice of $L_{n}$ is an MCMC kernal of invariant distribution $\pi_{n}$ associated with $K_{n}$:

[TABLE]

This leads to a weight update equation

[TABLE]

According to (5.1), the weight update equation is given by

[TABLE]

where $I_{\epsilon}$ is defined in (5.2). The pseudo-weights $\tilde{w}_{n}^{(i)}$ can be normalized, yielding $w_{n}^{(i)}=\tilde{w}_{n}^{(i)}\left(\sum_{i=1}^{N}\tilde{w}_{n}^{(i)}\right)^{-1}.$ It is also possible to use repeated simulations of the pseudo observations for a given parameter, rather than using a single pseudo observation as done here, to reduce the variance of the particle weights. This generalization is discussed in [33].

We also follow [33] to design an adaptive schedule for the tolerance levels, $\epsilon_{0}>\epsilon_{1}>\cdots>\epsilon_{T}=\epsilon$. In particular, $\epsilon_{n}$ is chosen to satisfy the criterion

[TABLE]

for $\alpha\in(0,1)$, and where $\#\{.\}$ denotes the number of elements in the set $\{.\}$. In practice the preceding three equations are solved in reverse order. One sorts the particles and associated trajectories $(\theta_{n-1}^{(i)},y_{n-1}^{(i)})$ according to their (Hellinger) distance $\rho(y_{n-1}^{(i)},y^{o})$. Of the particles with non-zero weight, a proportion $1-\alpha$ with the largest ‘distance’ as measured by $\rho$ has weight set to 0. The resulting weights are the desired weights $w_{n}^{(i)}$ for the $n$-th step; the tolerance level $\epsilon_{n}$ is given by the largest Hellinger distance among the non-zero weight particles, $\max_{\,n}\left\{\rho(y_{n}^{(i)},y^{o})|w_{n}^{(i)}>0\right\}$.

We implement the the MCMC kernel $K_{n}(z,z^{\prime})$ by the Metropolis-Hastings algorithm. In particular, we draw a candidate $z^{\ast}=(\theta^{\ast},y^{\ast})$ from a proposal $q_{n}(\theta,\theta^{\ast})f(y^{\ast}|\theta^{\ast})$, where $q_{n}(\theta,\theta^{\ast})$ is a transition kernel for $\theta$. It is then accepted with probability given by the MH ratio

[TABLE]

Note that the above MH acceptance-rejection process is undertaken only for the particles with non-zero weights since the zero-weighted particles will keep their zero weights until being “resurrected” via the Resampling algorithm.

The weight update equations (5.10) and (5.9) ensure that the maximal distance between the pattern associated with any particle and the observed pattern, as measured by (5.2), decreases. One immediate consequence of this is that one cannot add noise to the particles when resampling. If noise is added, some or indeed most of the refreshed particles will have associated distances larger than the threshold $\epsilon_{n}$, and the monotonic convergence of the algorithm will be violated. Rather than refreshing the particles by resampling, the Metropolis-Hastings algorithm in steps 3-4 serve to continually refresh the available particles. In particular the accept-reject criterion for the algorithm depends on the current threshold $\epsilon_{n}$ and grows more stringent as time passes.

The SMC-ABC algorithm is summarized as follows.

Algorithm 1

Initialization

Set $\epsilon_{0}=\infty$, $w_{0}^{(i)}=1/N$ for $i=1,\ldots,N$ and $n=1$. Draw $\theta_{0}^{(i)}\sim\pi(\cdot)$ and $y^{(i)}\sim f(\cdot|\theta_{0}^{(i)})$.

While

$\epsilon_{n-1}>\epsilon$

step 1

Solve (5.10) and (5.9) jointly for $\epsilon_{n}$ and $w_{n}^{(i)}$.

step 2

If $N_{eff}(\{w_{n}^{(i)})\})<N_{T}$, do Resampling.

for

$i=1,\dots,N$

step 3

$(\theta^{\ast},y^{\ast})\sim q_{n}(\theta_{n}^{(i)},\theta^{\ast})f(y^{\ast}|\theta^{\ast})$.

step 4

Accept $\theta_{n+1}^{(i)}=\theta^{\ast}$ with the probability (5.11); otherwise $\theta_{n+1}^{(i)}=\theta_{n}^{(i)}$.

end for

step 5

Set $n=n+1$.

End while loop

In order to assess the viability of PCA as a tool in data assimilation, we will compare the above SMC-ABC algorithm to several other schemes. In particular, we will consider the standard particle filter (see for instance [8, 53]). We will estimate the likelihood of the observations in the particle filter by assuming that their error statistics are distributed normally.

We outline a particle filter in the next section.

6 Particle Filters

For comparison with SMC-ABC that is used to assimilate coherent patterns, we use the standard (bootstrap) particle filter (PF) [22], where the importance density is the same as the transition density, to assimilate the tracer trajectories in the standard manner done in previous research [49, 1, 29, 43, 44]. It should be emphasized that while SMC-ABC estimates the target distribution corresponding to batch data from which the coherent patterns are extracted, the PF recursively solves the sequential data assimilation problem where the posterior distribution evolves at each measurement time through accumulation of the tracer positions. In particular, PF aims to sequentially approximate $\pi(\theta_{0:n},\mathbf{x}_{0:n}|\mathbf{x}^{o}_{0:n})$, where $\mathbf{x}^{o}$ denotes the observed tracer positions, $\mathbf{x}$ is the predicted positions and $\theta$ is the parameter vector. Since $\pi(\theta_{0:n},\mathbf{x}_{0:n}|\mathbf{x}^{o}_{0:n})$ is approximated by the weighted particles $\{(\theta_{0:n}^{(i)},\mathbf{x}_{0:n}^{(i)}),w_{n}^{(i)}\}$, the marginalization to obtain the particle distribution $\{(\theta_{n}^{(i)},\mathbf{x}_{n}^{(i)}),w_{n}^{(i)}\}$ can be readily done recursively. This particle distribution gives the empirical pdf for $\theta$, $\pi(\theta_{n}|\mathbf{x}^{o}_{0:n})$. As for the SMC-ABC method, $n$ is the time index and $i$ the index for each particle. Note that this is an abuse of notation for the sake of convenience, since the time index in the SMC-ABC algorithm is artificial but the time index in PF represents the actual sequence of the data.

In a nutshell, the (bootstrap) PF updates the weight according to

[TABLE]

As done for the SMC-ABC algorithm, we monitor the Effective Sample Size and resample if $N_{eff}<\delta N$ for some $0<\delta<1$. We employ an adaptive resampling scheme from [35], which merges particles together in the resampling step. In particular, instead of resampling $N$ samples from the ensemble, the merging scheme resamples multiple samples of size $N$ and generates a new sample of size $N$ using a weighted average of these multiple samples. This scheme provides a way to gain sample diversity and the weights are designed to preserve the sample and variance of the original sample. In all of our experiments, we use 3 samples of size $N$ to generate a new sample, merged with the weights found in [35].

For implementation purposes, we will model the likelihood of the observations by

[TABLE]

where $\mathbf{R}$ is the observational error covariance and, on the right of the equation, $\mathbf{x}_{n}^{(i)}$ depends on $\theta_{n}^{(i)}$. As described in the Introduction, we are interested in the highly-informative situation where the drifter locations are measured accurately. As such we take $\mathbf{R}=10^{-4}\,\mathbf{I}$ in the following experiments.

7 Case Study: Meandering Jet

7.1 Coherent pattern

Our test model follows the example in [41], Chapter 5. We consider a kinematic traveling wave model in the co-moving frame $(x_{1},x_{2})$ that is deterministically perturbed by an oscillatory disturbance and stochastically perturbed only in the $x_{1}$-direction:

[TABLE]

where the additive noise is described by a Wiener process $dW$ with variance $\sigma^{2}$. Here $A>0$ is the amplitude, $c$ is the phase speed of the propagation of the primary wave and $K$ is the wavenumber in the $x_{1}-$direction. The paramter $c_{1}$ is the phase speed of the peturbation term in the co-moving frame, $k_{1}$, $l_{1}$ are the the $x_{1}$- and $x_{2}$- wavenumbers, respectively. We assume the periodic boundary $[0,2\pi)$ in the $x$-direction. For $\sigma=0$, $\varepsilon=0$ and $c$ chosen appropriately, the flow is integrable and consists of two closed circulation regions adjacent to the rigid boundaries $x_{2}=0$ and $x_{2}=\pi$ and a jet flowing regime in the $x_{1}$-direction passing between the two circulation regions. The fluid exchanges under the lobe dynamic due to the oscillatory disturbance of (7.1)–(7.2) are given in [41].

In the following we use $\sigma^{2}=0.01$, $\varepsilon=0.3$, $A=1$, $K=1$, $c=0.5$, $k_{1}=1$, $l_{1}=2$, $c_{1}=\pi$, as given in [41].

For these parameters, the flow around the two circulation regions is dominated on a long time scale by chaotic advection [36]. A line of drifters placed in one of the gyres and initially spaced $\mathcal{O}(10^{-2})$ apart, that is spaced according to the observation error covariance $\mathbf{R}$, will elongate into a filament that persists to $t\approx 20$. By $t\approx 30$ the filament breaks and the drifters appear to be distributed homogenously within the gyre. The flow in the jet between the gyres is not dominated by chaotic advection on these time scales.

We now establish that for this flow there is a coherent pattern, that the coherent pattern can be detected using relatively few tracers, and that the Hellinger distance between simulated and observed coherent patterns is robust to chaos in a sense that the metric employed by the particle filter is not.

In Figure 1, we show the PCA-based coherent pattern for various values of $\varepsilon$ and for $\sigma=0$. We initialize tracer positions using $50\times 25$ uniform grid points on $(0,2\pi)\times(0,\pi)$ and integrate them to the time $t=10T$ for $T=2$ (the period of the flow under the above parameter values) using the standard Euler-Maruyama method with step size $\delta t=0.01$. The integration time is chosen to be long enough to capture distinct dynamics of the tracers and to inspect the robustness of the coherent pattern. The covariance matrix $\mathbf{P}_{x_{1}}$ from (2.1) is then constructed from the coordinates of the tracer trajectories taken as observations every $10\delta t$ time steps. As shown in Figure 1, although changing its shape with $\varepsilon$, the coherent pattern is very persistent in the sense that those initial tracers labelled in similar colours according to the PCA dominant mode are still tightly grouped together and retain their collective geometrical structure after a long period of time; this gives an intuition of what we mean here by the coherent pattern. Note also that, dynamically speaking, under the small (deterministic) disturbance, the flow is near-integrable and the coherent pattern is determined by the persisting geometry of the circulation regime and jet flowing regime.

We now confirm that the coherent pattern can be detected using relatively few tracers, and that the coherent pattern is in this case robust to the noise and chaotic advection in the dynamical system. We investigate this by limiting the number of tracers to 25, initially randomly placed. We label each tracer and use three realizations of the noise to propagate the tracers to the time $t=10T$. The coherent pattern is again extracted by using the dominant mode of PCA, and the tracers in each of the three realizations are coloured accordingly. We consider PCA to be robust to the chaotic, stochastic dynamics if tracers with the same label have similar colours across each of the 3 realizations. We show the initial conditions and 3 typical realizations in Figure 2.

We now establish the usefulness of the Hellinger distance as a metric to infer model parameters from ‘observed’ and simulated patterns, compared to the direct comparison of the simulated and observed drifter locations used by the particle filter in (6.2). Figure 3 compares these two metrics, using a single observation taken at either $t=20$ or $t=30$. We observe that at $t=20$, i.e. before chaotic advection has separated simulated drifters from observed, both the Hellinger distance and the distance between simulated and observed drifters has a clear minimum at the true value $\varepsilon=0.3$. If the observation is instead taken at $t=30$, after chaotic advection has separated the simulated and observed tracer trajectories, then the distance between simulated and observed drifters no longer has a minimum near the correct value of $\varepsilon$, but the Hellinger distance does. We used 50 drifters deployed in a circle of radius 0.1 around the centres $\left(\frac{\pi}{2},1\right)$, $\left(\frac{3\pi}{2},\pi-1\right)$ of each gyre.

7.2 Data assimilation by SMC-ABC

In this section we explore the robustness of the SMC-ABC algorithm, measured against the standard particle filter (PF). Following the discussion in Section 1 and the results of Figure 3, we define two ‘failure states’ for the PF. We expect the PF to become degenerate when there are many accurately observed drifters, meaning that one particle will have weight 1 and all other particles will have weight 0. However, we are estimating static parameters; it is quite possible for the (degenerate) particle filter to accurately estimate these parameters with a single particle. We would like to observe a second failure on the part of the PF, that if the drifters are located in regions of chaotic advection in the model, and if the time between observations is sufficient for the initial observation error to completely separate the simulated drifter trajectories from the true trajectories, then the PF may be expected to gain no information over the prior guess for the parameter(s).

We focus on estimating the parameters $\varepsilon,\,k_{1},\,l_{1},\,c_{1}$ (particularly $\varepsilon,\,k_{1}$) that characterise the time-dependent peturbation. In order to observe the above failure state in the PF, we will typically initialise the drifters in two rings around the centre of each gyre, and take observations on intervals at least $t=30$. By taking observations on time intervals less than 30, or as we shall see by initialising the drifters uniformly in phase space, the degenerate PF can estimate parameters more accurately than SMC-ABC. It is specifically the robustness of the coherent pattern to chaotic advection that will provide an advantage to SMC-ABC.

7.2.1 Experimental setup

In the following experiments, we use a final time of $t=30$. The PF and SMC-ABC schemes use a model time step of $0.1$, and a single observation is taken at the final time. The observations are synthetically generated from a model integration using a time step of $0.01$.

Our observations are of the trajectories of $50$ tracers; the initial positions of these tracers are typically circles of radius 0.1 around the centres of each gyre. We construct the “truth” by taking a realization produced by the above parameter values.

We are interested in comparing the performance of PF and SMC-ABC at the same computational cost. Accordingly, we modify SMC-ABC as follows: we choose the number of particles $N$ in SMC-ABC to be one tenth the number of particles used in the PF. Furthermore, we restrict the SMC-ABC algorithm to at most ten steps in the artificial time $n$; that is if $n=10$ in Algorithm 1 the While loop is halted and we return the final particle distribution $\{\theta^{(i)}_{10}\}$. With this modification the computational cost of SMC-ABC is controlled to be roughly equivalent to, or less than, the cost of PF.

In the following experiments, the transition kernel $q_{n}(\theta,\theta^{\ast})$ is given by the random walk

[TABLE]

Therefore, the MH ratio (5.11) is either 0 or 1; the weight of the $i-$th particle at the $n-$th step is 1 if it has survived all $n$ rejection steps, or if it has been resampled and then survived the intervening rejection steps, and it is 0 otherwise.

7.2.2 One parameter experiments

We perform experiments to estimate the time-dependent peturbation parameter $\varepsilon$, assuming that the prior density $\pi(\varepsilon)$ is a uniform distribution on the interval $[0,1]$.

We investigate how the absolute error $\left|\sum_{i=1}^{N}\theta_{n}^{(i)}w_{n}^{(i)}-\varepsilon\right|$ scales with the number of particles $N$, when all other parameters are kept fixed. We present results for PF with $N=\{250,500,1000,2000\}$ particles, and SMC-ABC with $N=\{25,50,100,200\}$ particles. As discussed in Section 7.2.1, under this set-up SMC will have lesser or equal computational cost to PF, and we overlay the results accordingly. We suppose that the true parameter is $\varepsilon=0.3$, and present results for two situations, one in which the drifters are initially distributed uniformly and a second in which the drifters are initially distributed in a circle of radius 0.1 around the centres $\left(\frac{\pi}{2},1\right)$, $\left(\frac{3\pi}{2},\pi-1\right)$ of each gyre. These drifters will all be contained within the gyres (see Figure 1). As discussed previously, in the first situation some drifters will be in nonchaotic or weakly chaotic regions of the flow, and the PF should gain information; in the second situation all trajectories will be dominated by chaotic advection and the PF should be totally uninformative. Figure 4 clearly shows the expected behaviour from PF, highlighting that SMC-ABC is unaffected by the chaotic flow. We mention that each run of the PF is degenerate (having one particle with weight 1).

7.2.3 Multiple parameter experiments

In this section we use SMC-ABC and PF to simultaneously estimate the perturbation parameter $\varepsilon$ and the amplitude and frequency parameter $k_{1}$ in (7.1)–(7.2). We label the set of parameters to be estimated by ${\Theta}=(\varepsilon,k_{1})$, with true values $\epsilon=0.3$ and $k_{1}=1$. We take the prior density $\pi({\Theta})$ to be a uniform distribution on the interval $[0,1]\times[0.2,1.2]$.

We investigate how the absolute error $\left|\sum_{i=1}^{N}\theta_{n}^{(i)}w_{n}^{(i)}-{\Theta}\right|$ scales with the number of particles $N$, when all other parameters are kept fixed. The drifters are again initialised inside the boundaries of each gyre, as done previously. We present results for 20 runs each of a PF with $N=\{250,500,1000,2000\}$ particles, and SMC-ABC with $N=\{25,50,100,200\}$ particles. Again we overlay the results in view of their roughly equivalent computational cost. Figure 5 confirms that SMC-ABC can accurately estimate multiple parameters, while again the PF fails to improve over the prior guess. We remark that the broad spread of the results for SMC-ABC compared to PF is a consequence of the adjustment to the same computational cost; SMC-ABC is (also) a sequential monte carlo method employing one tenth the number of particles as PF, and consequently outputs particles with relatively high variance.

8 Discussion

This paper presents, to our best knowledge, the first attempt to develop an approach to assimilate the Lagrangian coherent structure for parameter estimation. In contrast to previous research for LaDA, the new approach does not directly use the error statistics of observed positions of the tracers in order to make a parameter inference but assimilates the coherent structures as observations instead. Therefore, our approach assumes the existence of coherent structures and requires a tool that can extract that structure from the trajectory data. The structure is extracted via a machine learning method, which is PCA in the current work, and the assimilation algorithm is based on the Approximation Bayesian computation (ABC), which allows a bayesian inference to be performed without computing the typically unavailable likelihood function of the coherent structure. However, a distance function must be appropriately designed to measure the discrepancy between the observed coherent structure and the coherent pattern predicted by the model simulation based on various parameter values. We use the Hellinger distance and empirically demonstrate its robustness in the situation where the underlying flow exhibits chaotic advection and is stochastically perturbed.

Our numerical experiments demonstrate that this new approach is remarkably superior to the trajectory-based Lagrangian DA (employing the particle filter) in the situation where the number of tracers is large and the drifter trajectories are dominated by chaotic advection. In order to achieve this situation, we initialised drifters in the chaotic regions of a peturbed jet flow. The comparison uses the average mean-square error to measure the accuracy of the estimate. The improved accuracy of the ABC approach is attributed in large part to the following facts: (1) the Hellinger distances for the coherent structure is usually less sensitive to the stochastic perturbation than the likelihood function of tracer trajectories, which is assumed to be highly informative (i.e. to have small variance). (2) Even in the chaotic flow, the large-scale coherent structure is still geometrically “regular” while the tracer trajectories may experience high dispersiveness due to the stretching and folding process.

We now briefly discuss convergence results for the SMC-ABC algorithm. On the SMC side, [9] discusses how to obtain convergence results for the stochastically resampled SMC scheme by connecting it to deterministic resampling schemes (convergence results for which are obtained via the Central Limit Theorem and are given in, for instance, [48]). It may be more difficult to establish convergence for the ABC component, which would typically entail confirming that the coherent pattern (as approximated by some numerical method) is a summary statistic for the tracer statistics, in the sense that the likelihood function for the coherent pattern is the same as the likelihood function for the tracer trajectories; see for instance the Discussion in [30]. We comment that for the flow considered in this paper, the coherent pattern is more than just a summary statistic. In other words, the likelihood function given by the tracer trajectories is different to the likelihood function given by the coherent pattern.

In future work we will examine our approach in the case of transitory coherent sets, which exist for a finite time before losing coherency due to a change of parameters. Detecting the parameter change is of interest in many applications and it will require developing a new algorithm that exploits the coherent structure to track the change.

9 Acknowledgements

JM is grateful for the support of the MURI grant A100752, ONR grant N00014-15-1-2112, and NSF grant DMS-1312906. NS acknowledges the support of the Department of Mathematics at Univeristy of Surrey.