An Observer for an Occluded Reaction-Diffusion System With Spatially Varying Parameters

TL;DR

This paper develops an observer method for a reaction-diffusion system with spatially varying parameters, enabling parameter estimation and synchronization from partial, noisy, and occluded data, inspired by ocean ecology remote sensing.

Contribution

It introduces an autosynchronization approach for estimating parameters and states in a spatially varying reaction-diffusion PDE with occlusions, applicable to ocean ecology modeling.

Findings

Successfully estimates parameters from partial data

Handles noisy and occluded observations effectively

Analyzes stability using a network approach

Abstract

Spatially dependent parameters of a two-component chaotic reaction-diffusion PDE model describing ocean ecology are observed by sampling a single species. We estimate model parameters and the other species in the system by autosynchronization, where quantities of interest are evolved according to misfit between model and observations, to only partially observed data. Our motivating example comes from oceanic ecology as viewed by remote sensing data, but where noisy occluded data are realized in the form of cloud cover. We demonstrate a method to learn a large-scale coupled synchronizing system that represents spatio-temporal dynamics and apply a network approach to analyze manifold stability.

Equations

$$\frac{\partial P}{\partial t}$$

$$\frac{\partial Z}{\partial t}$$

$$k_1(x,y)$$

$$m_1(x,y)$$

$$k_2(x,y)$$

$$m_2(x,y)$$

$${\mathbf{u}}_{\mathbf{t}}(x,y,t)=\mathbf{f}(\mathbf{u}(x,y),\mathbf{p}(x,y)),$$

$${\mathbf{v}}_{\mathbf{t}}(x,y,t)=\mathbf{g}(\mathbf{u}(x,y),\mathbf{v}(x,y),\mathbf{q}(x,y)),$$

$${\mathbf{q}}_{\mathbf{t}}(x,y,t)=\mathbf{s}(\mathbf{u}(x,y),\mathbf{v}(x,y)),$$

$$\frac{\partial \hat{P}}{\partial t}$$

$$\frac{\partial \hat{Z}}{\partial t}$$

$$\frac{\partial \hat{k}}{\partial t}$$

$$\frac{\partial \hat{m}}{\partial t}$$

$$\displaystyle\hat{P}=\Bigg{\{}\begin{matrix}\hat{P}&:&0<\hat{P}<2\\ 0&:&\hat{P}\leq 0\\ 2&:&\hat{P}\geq 2\end{matrix}\textrm{ \hskip 11.38092ptand}$$

$$\displaystyle\hat{Z}=\Bigg{\{}\begin{matrix}\hat{Z}&:&0<\hat{Z}<2\\ 0&:&\hat{Z}\leq 0\\ 2&:&\hat{Z}\geq 2\end{matrix}$$

$$\frac{\partial \omega }{\partial t}+\nu \frac{\partial \omega }{\partial x}=0,$$

$$\frac{\partial \hat{P}}{\partial t}$$

$$\frac{\partial \hat{Z}}{\partial t}$$

$$H[P]=\{\begin{array}{ccc}\hat{P}& ,& P=I\\ P& ,& P\ne I.\end{array}$$

$$\frac{\partial \hat{P}}{\partial t}$$

$$\frac{\partial \hat{Z}}{\partial t}$$

$$\frac{\partial \hat{k}}{\partial t}$$

$$\frac{\partial \hat{m}}{\partial t}$$

$$\frac{\partial \hat{P}}{\partial t}$$

$$\frac{\partial \hat{Z}}{\partial t}$$

$$\frac{\partial \hat{k}}{\partial t}$$

$$\frac{\partial \hat{m}}{\partial t}$$

$$G_n(t)=\frac{1}{(dx)(dy)}\sum _{x,y\in S_n}(P(x,y,t)-\hat{P}(x,y,t)),$$

$$a_{ij}=\{\begin{array}{l}{w}_{ij}\mathrm{with}\mathrm{probability}{p}_{ij}\\ 0\mathrm{with}\mathrm{probability}1-p_{ij},\end{array}$$

$$\displaystyle L_{2}(t)=\left[\begin{array}[]{cc}\textbf{0}&\textbf{0}\\ I_{N}&-I_{N}\end{array}\right],$$

$${\dot{\mathbf{u}}}_q(t)$$

$$\dot{\delta u}(t)=(F+{\sigma }_1{L}^1\otimes B_1+{\sigma }_2{L}^2(t)\otimes B_2)\delta u(t),$$

$$F$$

$$J_i$$

$$\delta u(t)=(W\otimes I_2)\zeta (t)+e\otimes \delta u_s(t),$$

$$\zeta (t)=(W\otimes I_2{)}^T\delta u(t),$$

Full text

Remote sensing, synchronization, autosynchronization, parameter identification, assimilation, moving neighborhood network

††preprint: AIP/123-QED

Research in large-scale oceanic phenomena is made possible by remote sensing instruments mounted on ocean-observing satellites. These instruments provide datasets that can be filtered to study sizable ecological events, including harmful algal blooms. The fact that datasets are often patchy when clouds hide regions in the spatial domain is a substantial difficulty when attempting to parameterize a dynamical system. To attack this problem we extend a recently developed autosynchronization method. Model parameters and states are evolved in a drive-response pattern, on a-priori known model equations, to learn model states and parameters even while data are considerably spatially occluded. It has been shown that, assuming the model structure to be known, a synchronization system can be designed to effectively act as an observer to identify system parameters, even in a large scaled network system. While a discretized PDE can be interpreted as a particular lattice network, the realistic problem of cloud occlusions will cause times where the observer network is essentially disconnected. Our prior work has shown that synchronization can exist even in a large scale network that is not fully connected but rather has a so-called fast blinking structure. The method is analyzed by interpreting the discretized PDE as a large-scale coupled moving neighborhood network.

I Introduction

Algae form the basis of the food chain in the oceans and are ultimately responsible for providing nourishment for other marine life further up the food chain Fasham (1993). Seasonal environmental heterogeneities such as nutrient replenishment, predation, and temperature provide favorable conditions for recurring algal blooms, often called spring blooms. More localized bloom events are commonly observed in estuaries and coastal regions Truscott and Brindley (1994). Certain bloom events, especially harmful algal blooms, elicit widespread repercussions on regional communities including human sickness, shellfish poisoning, and fish kills Backer and McGillicuddy (2006). These harmful algal blooms are detrimental to regional ecology and economies through fishing losses and tourism depletion. Models for near-shore algal blooms would be extremely useful for forecasting during such events and might help inform short-term management decisions.

Parameter and state identification based on observed data remains an important topic in both dynamical systems and control theory. Several powerful methods for parameter estimation of spatio-temporal systems include Kalman filter methods Schiff and Sauer (2008); Annan et al. (2005); Wan and Van Der Merwe (2000), multiple shooting methods Müller and Timmer (2004, 2002), and synchronization methods Pecora and Carroll (1990); Parlitz (1996); Parlitz, Junge, and Kocarev ; Stojanovski, Kocarev, and Parlitz (1996); Yu et al. ; Yu and Parlitz (2008); Sorrentino and Ott (2009); Quinn et al. (2009); Schumann-Bischoff et al. ; Berg, Luther, and Parlitz (2011); Schumann-Bischoff and Parlitz (2011). Autosynchronization is a special variation of synchronization methods based on an approach to force a response model to adapt to observed data by developing additional equations for the parameters that depend on the synchronization error Parlitz (1996); Sorrentino and Ott (2009). Our implementation of the method assumes prior knowledge of the model structure. Recently, it has been shown that it is possible to estimate spatially dependent parameters for a PDE system by autosynchronization using a combination of diffusive and complete replacement coupling of observed data (drive model) to force the response model and parameters to synchronize with observables Kramer and Bollt (2013).

Our interest here is to exploit these ideas toward modeling ocean ecology as informed by hyperspectral remote sensing data captured by ocean observing satellites. Many well accepted ocean ecology models include predator-prey dynamics between at least two components: zooplankton, the predator, and phytoplankton, the prey Edwards and Bees (2001); Upadhyay, Kumari, and Rai (2009); Malchow et al. (2000); Medvinski et al. (2002); Freund et al. (2006); Malchow et al. (2005); Truscott et al. (1994); Steffen, Malchow, and Medvinsky (1997); Scheffer et al. (1997); Edwards and Brindley (1996); Truscott et al. (1994); Matthews and Brindley (1997); Edwards and Brindley (1999); Fasham, Ducklow, and McKelvie (1990); Scheffer (1991). Data observations are often noisy or patchy, particularly when observing spatio-temporal systems. The usual hurdle to fitting and subsequently solving a predator-prey reaction-diffusion system as informed by remote sensing data is the inability to observe zooplankton. As of now, there exists no method to estimate zooplankton densities based on hyperspectral inferences. Here, we adapt the method of autosynchronization of PDEs to be used with less available information, where noisy data are occluded by clouds.

At the heart of the problem is the observability of the dynamical system based on available sampling data, in this case phytoplankton. The problem of observability on nonlinear systems has been a topic of research over the past decade and is now much better understood Letellier and Maquet (2005); Letellier and Aguirre (2005); Bianco-Martinez, Baptista, and Letellier (2015). We therefore demonstrate that the system we study is observable from the variable provided by remote sensing data. We note that one might first check if the corresponding ODE system is observable by investigating the invertibility of the Jacobian of the differential embedding map of observed samples Letellier and Maquet (2005). Such a result would provide hope that a search for an autosynchronization scheme is worthwhile.

We begin by introducing the reaction-diffusion equations used to create synthetic observed data. Next, we assign the response system and discuss an autosynchronization configuration. We show the method can work with significant proportions of the data unobservable, e.g. data occluded by cloud cover. Finally, we consider the large-scale coupled synchronization system as a moving neighborhood network and apply a theorem for synchronization based on the rate of switching between network topologies to prove that our system can synchronize. It is shown that as long as the average network corresponding to the graph Laplacian supports synchronization and the switching epoch between new samples of network topologies is small enough, synchronization is achieved. Therefore, it is feasible to realize model fitting and data assimilation for multi-component ecological systems with realistic remote sensing data.

II Model Dataset

Satellite data of plankton blooms often reveal complex mesoscale structures such as ocean gyres and eddies for which there are several theories. As a synthetic dataset, the spatiotemporal model for plankton ecology should have the capability to render mesoscale structures. Medvinksi, et al, Medvinski et al. (2002) describe a two-component predator-prey model, including phytoplankton and zooplankton, over a rectangular two-dimensional region. Given perturbed initial conditions, the model exhibits spiral patterns on a spatial scale comparable to that which is observed in nature. By sampling snapshots from the solution of this model, we emulate a satellite image dataset. The dataset is complicated by including spatially varying parameters. This is a valid consideration when modeling mesoscale ocean ecology. Consider the system of two PDEs as given in Medvinski et al. (2002),

[TABLE]

where $P(x,y,t)$ represents phytoplankton density, $Z(x,y,t)$ represents zooplankton density, and both are observed on a compact connected two-dimensional domain, $\Omega$, with zero-flux boundary conditions.

These equations represent a dimensionless reaction-diffusion model for phytoplankton-zooplankton ecology, invoking predator prey dynamics in the reaction term. The ecology is considered over a horizontal layer with homogeneous vertical distributions in the water column. Our simulations are computed over a grid of size $\Omega=864\times 288$. The model assumes that phytoplankton obey a logistic growth and are grazed upon by zooplankton following a Holling-type @slowromancapii@ functional response. The Holling-type @slowromancapii@ functional response Holling (1959) assumes a decelerating growth rate wherein the predator is limited by its ability to efficiently process food. Zooplankton grow at a rate, $k$, proportional to phytoplankton mortality and die according to a natural mortality rate $m$. For scalar parameters, $k=2$, $h=0.4$, and $m=0.6$, and nonuniform initial conditions, this system gives rise to transient spiral pattern behavior, and progresses into spatially irregular patchy patterns Medvinski et al. (2002). We perform numerical simulations with a basic forward-time and central-space discretization using the perturbed initial conditions found in Medvinski et al. (2002).

The system Eq (1) is modified as found in Medvinski et al. (2002) by allowing the parameters to be nonnegative $C^{0}(\Omega)$ functions. Generally, we may allow $\Omega\subset\mathbf{R}^{2}$ to be a compact domain such as a rectangle for simplicity or a realistic domain representing a coastal region obtained from a satellite. Two examples are found in Figure 1, where high a concentration of phytoplankton appears as a greenish coloring of the water. Imaging sensors mounted on satellites measure light in discrete bandwidths, including several bandwidths outside of the visible range. These bandwidths are subsequently combined to build certain products of interest. To reconstruct an image as the eye would see it, bandwidths in the visible spectrum are combined to build what is called a “quasi-true” image. The quasi-true color image at the top of Figure 1 was taken on July 8, 2010 from the HICO (Hyperspectral Imager for the Coastal Ocean) instrument mounted on the Japanese Experiment Module Exposed Facility on the International Space Station. It is the first such imaging spectrometer specifically designed to sample the coastal ocean HIC . The image captures the Columbia River mouth bordering Oregon and Washington. The domain is large enough to render mesoscale and small scale patterns, which may result from complex intra-species and fluid dynamics. The image at the bottom of Figure 1 was taken by the MERIS (Medium Resolution Imaging Spectrometer) instrument on board the Envisat satellite. Here again high phytoplankton concentrations appear as a greenish coloring in the water. This image highlights a presently unavoidable issue with hyperspectral satellite data: the presence of cloud coverage.

In many systems, it is quite reasonable to expect that model parameters need not be spatially homogeneous. And therefore, taking our problem of interest, spatial inhomogeneity in parameter values may be an important theoretical assumption when constructing models for coastal algal blooms, since the plankton growth rate is affected by near-shore nutrient runoff and upwelling Thomann (1975); Medvinski et al. (2002); Mattern, Fennel, and Dowd (2012). More to that point, ocean fronts and eddies cause flow-induced long-term inhomogeneities in the ocean which results in a formidable spatial structure for density profiles in the ocean Medvinski et al. (2002). Whether inhomogeneities be the result of the flow dynamics or of boundary conditions from nutrient runoff, they are an important consideration for modelling ecology over large coastal domains. Thus, depending on the scale and resolution, it may be prudent to include spatially dependent parameters.

Therefore, we develop synthetic datasets with spatially varying parameters to challenge our methods. To push our methods we add random noise to each parameter as displayed in Figure 2. Spatially dependent parameters are chosen to be in the range given in Medvinski et al. (2002) for spatially irregular behavior. Three different functional forms for the parameters are tested for variety. First, we define a Gaussian parameter function,

[TABLE]

where $a=2,c=0.6,m=300,n=900$, and $\sigma=400$. Appropriate parameters are chosen to maintain $m(x,y)$ and $k(x,y)$ in the target range. Figure 2 shows the three parameter forms discussed above, where only $k(x,y)$ is plotted since the parameters differ by a scalar multiple. For example, Eq (2) is displayed in Figure 2a. Next, we define,

[TABLE]

where $a=0.2,b=\pi/(m/2),c=0.6,d=\pi/2,s=0.5$, and $t=1.5$, to test the quality of the autosynchronization method to resolve fine spatial structures in model parameters. The surfaces produced by Eq (3) are displayed in Figure 2b.

Finally, we build a swirly parameter function in order to simulate spiral-like behavior in parameter values as might be expected in turbulent coastal regions. A time instance is sampled from a simulation of the original PDE, Eq (1), is scaled appropriately, and is treated as a parameter function. These spiral parameters, $k_{3}(x,y)$, are shown in Figure 2c.

We discretize the modified system, Eq (1), with explicit finite differences, using a five-point center difference stencil for spatial derivatives and forward Euler time stepping. The spatial and temporal step sizes are chosen as $dx=2$ and $dt=0.2$. The model output $P(x,y,t)$ is treated as an image sequence given by a particular (known) model form but with parameters $k(x,y)$ and $m(x,y)$ and component function $Z(x,y,t)$ to be determined.

In order to properly mimic our target application of remote sensing oceanographic data of hyperspectral images filtered to reveal plankton blooms, we add random noise and “moving cloud cover” to the dataset by occluding large proportions of the image from direct observation. Clouds are a natural occurrence when studying a large terrestrial area over several days, and luckily the clouds tend to move.

III Autosynchronization

Two model systems are required in order to estimate unknown model states and parameters by autosynchronization, a drive system and response system. One-way direct replacement and diffusive coupling are combined so that observables are coupled directly into the response model as it evolves. Samples are taken from the drive system,

[TABLE]

with parameters $\mathbf{p}(x,y)\in C^{0}(\Omega)$ and $u\in H^{2}(\Omega)$. A response system is formed,

[TABLE]

with $\mathbf{q}(x,y)\in C^{0}(\Omega)$. We formulate an associated system of PDEs for the parameters of Eq (5),

[TABLE]

with the goal that $\mathbf{(v,q)}\rightarrow\mathbf{(u,p)}$ as $t\rightarrow\infty$. If successful, the method is called autosynchronization Kocarev and Parlitz (1996) since the parameters are evolved deterministically along with the response model.

Generally, some model variables from the drive system need not be sampled. For a two-species system, we write $\displaystyle\textbf{u}(x,y,t)=(u1(x,y,t),u2(x,y,t))^{T}$ and we do not require that $\displaystyle u2(x,y,t)$ is sampled. An associated response system $\displaystyle\textbf{v}(x,y,t)=(v1(x,y,t),v2(x,y,t))^{T}$ is built wherein both equations are fed samples from $\displaystyle u1(x,y,t)$. A schematic diagram for this type of simulation might be helpful and is found in Figure 3, where dots denote time derivatives.

For our synthetic dataset given by Eq (1) with parameters Eq (2) or Eq (3), we form a response system to be synchronized to the observations as,

[TABLE]

where we assume $\hat{P}(x,y,0)\neq P(x,y,0)$, $\hat{Z}(x,y,0)\neq Z(x,y,0)$, $\hat{k}(x,y,0)\neq k(x,y)$, and $\hat{m}(x,y,0)\neq m(x,y)$.

Parameters are updated as diffusively coupled PDEs during the synchronization process as,

[TABLE]

where $s=30$ and $\kappa=2.4$ are chosen for specificity and for which we observe good convergence results. For these experiments, we sample the drive system at every time step, but note that a larger sampling time will work Kramer and Bollt (2013). The parameter equations are evolved simultaneously by Eq (III) with a forward Euler discretization and the same time step. As we vary $s$ and $\kappa$, autosynchronization may fail as commonly observed with diffusively-coupled systems.

To begin the simulation, parameters are initialized as the constant function, e.g. $\hat{k}(x,y,0)=5$ and $\hat{m}(x,y,0)=5$. We evolve Eq (1) forward and count the model output as observed data. Initial conditions for the response system are $\hat{P}(x,y,0)=2$ and $\hat{Z}(x,y,0)=2$. Furthermore, to avoid values outside the normal range of Eq (1), we enforce that

{parcolumns}

[colwidths=1=1 in,nofirstindent]2 \colchunk[1]

[TABLE]

\colchunk

[2]

[TABLE]

$\forall x,y\in\Omega$ during the simulation. As noted above, autosynchronization is observed for the test set of parameters in Figure 2 and the spatial inhomogeneities in each case are effectively resolved. We emphasize that zooplankton are not observed in Eq (III)-(III).

IV Hidden Data

Ocean-observing satellite imagery often includes significant amounts of cloud cover Martin (2014). In other words, a large fraction of that data may be occluded. Furthermore, we have found that level 2 mapped and processed images may include striping or other defects from projecting a sphere onto a uniform grid. The lack of data presents a challenge to data assimilation and model filtering by synchronization methods. Suppose $\omega\subset\Omega$ is the set of unobservable data. We allow for $\omega=\omega(x,y,t)$ so that the set of unobservables varies with space and time like a cloud. We consider a simple case where the dynamics of $\omega(x,y,t)$ are governed by the advection equation

[TABLE]

with periodic boundary conditions, so that clouds move in the x-direction with speed $\nu$. We couple the systems only on the complement of $\omega$. That is, we turn the driving signal off when the image is unobservable, allowing the two systems to oscillate independently, and switch it on after the clouds have passed. We do this only in the subregion $\omega\subset\Omega$ that is unobservable in order that data contained in the complement of $\omega$ may continue to be driven by observables toward the synchronization manifold.

Here we build on the method described in section III, where zooplankton densities and model parameters are estimated by observing solely the phytoplankton. Now we observe phytoplankton and clouds. However, if we couple at every spatial grid point, the synchronization manifold is de-stabilized by incident cloud coverage. With a large enough amount of cloud coverage over $\Omega$, the systems fail to synchronize. We say large enough in deference to the case where the occluded region is small enough such that diffusion allows information to pass into any hidden regions.

As a remedy we allow the drive and response models to oscillate independently, or uncoupled, while the drive model is hidden by clouds. The pixels representing cloud cover in remote sensing data are typically set to some large fixed integer, $I$. We represent this formally,

[TABLE]

where $H[P]$ represents a switching function given by

{parcolumns}

[colwidths=1=3 in,nofirstindent]2 \colchunk[1]

$\displaystyle H[P]=\left\{\begin{matrix}\hat{P}&,&P=I\\ P&,&P\neq I.\end{matrix}\right.$

The form of response model switches off the coupling when a cloud mask is detected in the image and allows the systems to oscillate independently in the corresponding pixels, while being driven over pixels that are observed. Eq (IV) is slightly different from temporal subsampling of data, where models are not coupled for a given number of time steps. Here the models are always coupled somewhere in $\Omega$, which is determined by time-varying clouds.

We first demonstrate model state synchronization given occluded data before addressing the estimation of model parameters. Let Eq (1) be the drive model and Eq (IV) be the response model. Figure 4 represents a partially observed dataset from Eq (1), with a field of 30 randomly placed synthetic clouds evolving from left to right with periodic boundary conditions resulting in $65.8\%$ of $\Omega$ occluded at all times. The clouds repeatedly scroll from left to right and parts of the image are always occluded, but every element in the domain is eventually driven, causing the drive and response to systems to synchronize. The response system is initialized by $\hat{P}(x,y,0)=2$ and $\hat{Z}(x,y,0)=2$, and we choose $\kappa=2.6$.

Once synchronized, even hidden phytoplankton are revealed for initializing short-term forecasts, demonstrating the utility of this result. Figure 5 demonstrates that despite $65.8\%$ of the drive system hidden, the two PDE systems eventually evolve toward identical synchronization. In Figures 4e and 4f nearly all evidence of clouds is “synchronized away” from the response system and the globally averaged error between the two has been driven to be less than $2.6\times 10^{-12}$. We remark that the choice of coupling strength, $\kappa$, varies with the amount of data occluded.

Given a model form, we advance the method to sample a single species toward parameter estimation and nonlinear data assimilation for a two-species PDE model, regardless of clouds. That is, by stating the response system

[TABLE]

where $k(x,y)$, $m(x,y)$, $Z(x,y,t)$, and $P(x,y,t)|_{\omega}$ are to be estimated by $\hat{k}(x,y)$, $\hat{m}(x,y)$, $\hat{Z}(x,y,t)$, and $\hat{P}(x,y,t)$ by sampling only $P(x,y,t)|_{\omega^{C}}$. As before, the coupling is turned off completely for the pixels on which clouds are detected. The drive system is still Eq (1), but we allow for spatially dependent model parameters $k(x,y),m(x,y)$ in the form of Figure 2. For robustness, we consider that the model parameters need not have the same functional form. We choose $m(x,y)$ defined by Eq (3) and $k(x,y)$ with form shown in Figure 2c and also add random noise to both parameters.

Figures 6 and 7 demonstrate a comparison between drive and response models. In the top of Figure 6a, we see the observed system $P(x,y,t)|_{\omega^{C}}$ wherein $25.5\%$ of the data on $\Omega$ is not observable. Figures 6 and 7 demonstrate that phytoplankton, zooplankton, and both spatially dependent parameters $\hat{k}(x,y,t)$ and $\hat{m}(x,y,t)$ are estimated to high precision.

Figure 8 describes the globally-averaged relative error between the true system and the response system. The rate of convergence to the synchronization manifold is slower than with cloudless data as a result of allowing the systems to oscillate independently while not driven on $\omega$. For the simulation in Figures 6 and 7, we choose $\kappa=0.625$, $s_{1}=0.2$, and $s_{2}=0.6$ for good autosynchronization results. Summarizing, we have demonstrated that it is possible to fill in missing data when hidden by clouds and, as an added bonus, estimate noisy spatially-dependent model parameters with different functional forms. Similar results are obtained by testing other combinations of the two forms of model parameters in Figure 2.

Simulations are run for varying percentages of hidden data. Figure 9 shows the synchronization errors for simulations after a fixed time epoch of $t=2400$ for all simulations. Specifically, the globally-averaged relative error between drive and response systems is plotted against the percentage of hidden data. It is clear that the speed of convergence of assimilation slows with respect to the degree of occlusion. A counter-intuitive side note is that the assimilation quality actually improves by hiding data through about $13\%$ before worsening as a larger percentage of data is hidden. For these simulations, the same initial conditions are used throughout for consistency.

We acknowledge inherent noise in remote sensing data and next demonstrate the method can work with fairly noisy observations. We add random noise to $P(x,y,t)$ for the length of the simulation, occlude $P(x,y,t)$ with $25.5\%$ cloud coverage, and consider the same noisy, mixed functional form model parameters as above. As expected, as more noise is added to observations, the rate of synchronization slows and the error after a fixed epoch increases. Figure 10 includes numerical experiments in which increasing noise is added to $P(x,y,t)$ and synchronization errors are compared over a fixed time epoch of $t=2400$.

To be considered practical, the method should provide decent results if data are available on a courser grid than that on which the model is evolved.

V Samples on a Coarse Grid

Model simulations require data on an appropriately resolved grid, which may not line up with the grid on which samples are available. Often is the case that sampled data exist on a coarsened grid relative to a required simulation grid. We demonstrate that a simple modification of the technique will produce results similar to those above. We imagine that the domain is sampled in discrete patches, denoted by $S_{n}$, and on the patches we have available only local averages of true data as depicted in Figure 11. In this way, we sample a coarser subset of the domain and take local averages to be the new driving signal.

To adapt the problem to such data, we modify the response system, Eq (IV)

[TABLE]

where $\tilde{P}$ represents locally averaged observations from the drive system and

[TABLE]

where $S_{n}$ is the rectangular “sensor” on which local averaging occurs.

Note, for best results, diffusion is added to the parameter equations in Eq (V) in order that data from the driven regions, $S_{n}$ will diffuse into the occluded regions. To further mimic our target problem of remote sensing data, we add white Gaussian noise to the observations. We now demonstrate the method reconstructs parameters and unknown states on incomplete, noisy, patchy experimental data.

Figure 12 includes the results of simulations with coupling given by Eq (V), with a sensor size of 2 pixels by 2 pixels, and a spacing of 1 pixel between sensors. Here we show the end of the simulation for brevity. It is obvious the method suffers from noise, local averaging, and missing data between sensors. In Figure 12a, it is clear a considerable amount of data is occluded, the data that are available include noise, and still unknown states and parameters are reconstructed fairly well.

In Figure 13 we provide simulations to demonstrate the dependency of synchronization quality on the sensor spacing. The synchronization errors for zooplankton and the parameter $k(x,y)$ are plotted against time for simulations admitting zero, one, and two pixels between sampling sensors on the domain. If sensor spacing is too sparse, the method struggles to fill in data between patches on which local averages are provided, and ultimately fails. We note that the parameter $m(x,y)$ and phytoplankton exhibit similar behavior.

These results demonstrate potential for future application to real data. In fact the autosynchronization method is capable of data assimilation by revealing hidden phytoplankton and zooplankton densities and estimating model parameters on noisy, coarse, cloudy data. Furthermore, our simulations indicate that state estimation by synchronization is more robust to coarse data than autosynchronization, which may be of use if parameter estimates are not required. This is merely a demonstration of autosynchronization applied to a remote sensing problem. These results require analytical reinforcement, including a discussion of the basin of attraction for the synchronization manifold and allowable coupling strengths to observe synchronization and parameter estimation. In the following, we provide analysis of manifold stability.

VI Analysis

To better understand this tendency to synchronize despite hidden data, we are inspired by a method from network theory. We represent the system as a moving neighborhood network and define each pixel in the image domain $\Omega$, to be an individual oscillator, $\mathbf{u_{j}}$. We include the drive and response images in the network so that an $m\times n$ image provides $2N=2mn$ oscillators.

Therefore, each drive oscillator $\mathbf{u_{j}}$, for some $j=1:N$, feeds a response oscillator corresponding to the same spatial pixel, $\mathbf{\hat{u}_{j}}$, for some $j=N+1:2N$. The drive system is hidden for a time epoch $\Delta$ to represent intermittent cloud cover. Pixels over which there are clouds are uncoupled while covered. In time the network topology shifts thereby shifting the coupling between drive and response oscillators over $\Omega$.

The adjacency matrix of our directed weighted random graph has zeros on the main diagonal and is defined as:

[TABLE]

for $i\neq j$. We subtract the adjacency matrix from the matrix $D$, with nonzero elements on the main diagonal $d_{i}=\sum_{j=1}^{n}a_{ij}$, to obtain the graph Laplacian $L$. Our network requires edges between nearby neighbors for diffusion, and between images for both direct replacement and diffusive drive-response coupling. We choose to order the 2$N$ vertices representing image pixels from left to right, taking a row at a time, and placing the drive image first, followed by the response image. That is we stack subsequent rows of the drive image followed by the response image to build the vector of 2$N$ components.

Here we analyze the occluded synchronization system (IV). We define two Laplacians to represent diffusion and drive-response coupling. $\displaystyle L_{1}\in\mathbb{R}^{2N\times 2N}$ is a sparse matrix with weighted entries corresponding to a chosen diffusion stencil and boundary conditions. $\displaystyle L_{2}(t)\in\mathbb{R}^{2N\times 2N}$, with elements $l^{2}_{ij}(t)$ is a time-dependent sparse matrix that represents diffusive coupling with a switching network topology and is fixed in the time interval $T_{k}=[t_{k},t_{k+1})$

[TABLE]

where 0 is an $N\times N$ zero matrix. Next, we define two matrices to model coupling between species, $B_{1}=I_{2\times 2}$ and $B_{2}=[1\ 0;0\ 0]$. The drive dynamics for the first $N$ oscillators and the response dynamics for the remaining $N$ oscillators are

[TABLE]

respectively, where $\mathbf{u}_{q}\in{\mathbb{R}}^{2}$ is the state of the $q$th oscillator, $\displaystyle\sigma_{1}=\frac{1}{h^{2}}$, $\displaystyle\sigma_{2}=\kappa$, and $\mathbf{f}:\mathbb{R}^{2}\to\mathbb{R}^{2}$ describes the individual dynamics of each oscillator.

We linearize the system about the synchronization manifold as

[TABLE]

where

[TABLE]

[TABLE]

We decompose (17) into a component that evolves along the synchronization manifold and a transverse component with a matrix $\displaystyle W\in\mathbb{R}^{2N\times(2N-1)}$. The state vector $\displaystyle\delta u(t)$ is decomposed, with $e$ denoting the standard basis vector in $\displaystyle\mathbb{R}^{2}$, as

[TABLE]

where

[TABLE]

and

[TABLE]

The linearized dynamics (17) are partitioned as

[TABLE]

where almost sure asymptotic synchronization of (16) is observed if $\displaystyle\zeta(t)$ almost surely converges to zero Kushner (1967).

It has been proven Porfiri, Stilwell, and Bollt (2008) that if the un-occluded system is uniformly asymptotically stable, then so is the occluded system, provided that the system is observable often enough. For completeness we restate Theorem 1 found in Porfiri, Stilwell, and Bollt (2008):

Theorem 1**.**

Consider the deterministic dynamic system:

[TABLE]

representing the linearized transverse dynamics of

[TABLE]

where $E[L]$ is the time-averaged graph Laplacian. Assume that F(t) and B(t) are bounded and continuous for all $t\geq 0$. If (21) is uniformly asymptotically stable, there is a time-scale $\Delta^{*}>0$ such that for any shorter time-scale $\Delta<\Delta^{*}$, the stochastic system (16) locally asymptotically synchronizes almost surely.

Now we check that the hypothesis of this theorem holds in the scenario of interest to us in this paper, which is that the deterministic graph Laplacian supports synchronization and the switching period $\Delta$ between network topologies is sufficiently small.

Corollary 1**.**

Consider the stochastic system (16), where $L_{2}(t)$ represents the switching network topology induced by moving clouds with speed $\nu$. Suppose $\displaystyle\omega$ is the set of unobservable data over the domain $\displaystyle\Omega$. If (21) is uniformly asymptotically stable, there exists a $C>0$ and $\nu^{*}(C)$, such that if $||\omega||<C$ and $\nu>\nu^{*}(C)$, then the time-scale between switching network topologies is sufficiently small, $\Delta<\Delta^{*}$, and the stochastic system locally asymptotically synchronizes almost surely.

Proof.

The proof of Theorem 1 found in Porfiri, Stilwell, and Bollt (2008) need only be altered slightly for our modification. We modify the definition of $\displaystyle M(t)$ so that

[TABLE]

We note that both $B_{1}$ and $B_{2}$ are bounded and continuous for $t>0$, and there are positive constants $\displaystyle\beta_{1},\beta_{2},$ and $\displaystyle\beta$, such that for any $\displaystyle t\geq 0$, $\displaystyle||B_{1}||\leq\beta_{1},\textrm{ and }||B_{2}(t)||\leq\beta_{2}$. Then setting $\displaystyle\beta=\mathrm{max}(\beta_{1},\beta_{2})$, the remainder of the proof follows. Therefore, we are guaranteed a time-scale, $\triangle^{*}>0$, below which the system will asymptotically synchronize. ∎

The referenced time scale is the time epoch under which the drive system is occluded. That is, if the drive system is occluded for too long, corresponding to clouds that are too large or plentiful, the systems will not exhibit synchronization.

To demonstrate Corollary 1 with respect to our problem, we simulate with varied rates of cloud movement. As cloud movement slows, $\Delta$ increases, and the graph Laplacian switches less frequently. If cloud movement slows too much, then the systems will not exhibit synchronization. Figure 14 includes the results of simulations over a fixed time epoch of $t=2400$. Parameter and state errors increase as the rate of cloud movement decreases.

VII Summary

We have extended the method of autosynchronization for PDEs with spatially dependent parameters to partially observable noisy PDEs with noisy parameters. We have shown that two PDEs can synchronize together even if the drive system is largely unobservable. Furthermore, all model states and parameters are estimated by sampling only one partially available noisy state. The work is a step toward modeling ocean ecology by remote sensing data over coastal regions or regions with recurring algal blooms. Remote sensing data in the form of hyperspectral satellite imagery often suffer from cloud coverage occluding parts of the domain over which we observe. Parameters and model states are estimated by treating the drive and response systems as independent oscillators on the unobservable set and periodically driving the dynamics toward the synchronization manifold.

Future work in this area includes adapting the method to work with reaction-diffusion-advection models, and building tools to optimize parameter and state estimation given extremely sparse data. To find a Lyapunov function for the system would provide strong theoretical backing. Furthermore, one might study the observability of the system, or similar systems, in order to understand which state variables are necessary to be sampled for successful estimation Letellier and Maquet (2005); Letellier and Aguirre (2005)

Algal blooms, especially harmful algal blooms, can have widespread negative consequences on local fisheries and tourism. In effect, models could inform management decisions and provide forecasts for local communities. Coupled with optical flow techniques Luttman et al. , advection could be added to reaction-diffusion models for additional accuracy, particularly over regions for which there are no vector field data for ocean currents. With advection data, techniques such as Finite-time/size-Lyapunov exponents Haller (2002); Bollt et al. (2012) or coherent sets based on transfer operator theory Haller and Poje (1998); Haller (2000); Froyland and Padberg (2009); Bollt et al. (2012) could be used to analyze the coastal dynamics to uncover Lagrangian coherent structures that might inhibit transport between regions of the ocean. Several such techniques could be used in concert to build bloom forecasts for coastal communities, informed by remote sensing data.

Acknowledgements.

This work has been supported by the Office of Naval Research under N00014-15-1-2093, the Army Research Office under N68164-EG, and a faculty development grant from Norwich University. We would like to thank Dr. Nicholas Tufillaro at Oregon State University’s College of Earth, Ocean, and Atmospheric Sciences for providing the satellite image products in the paper.