class: center, middle .title[Progress in Studying Ocean Turbulence from Lagrangian Trajectories]           .author[Jonathan Lilly⁰] <!--.center[[www.jmlilly.net](http://www.jmlilly.net)]--> .coauthor[Jeffrey Early¹, Adam Sykulski², Sofia Olhede³] .institution[⁰Theiss Research, ¹NorthWest Research Associates, ²Lancaster University, ³University College London]   .date[October 25, 2018]   .center[[{www.jmlilly.net}](http://www.jmlilly.net)] .footnote[Created with [{Remark.js}](http://remarkjs.com/) using [{Markdown}](https://daringfireball.net/projects/markdown/) + [{MathJax}](https://www.mathjax.org/)] --- class: center ## The Global Surface Drifter Dataset <video preload="auto" width="100%" height="auto" data-setup="{}" autoplay loop controls><source src="../videos/driftermovie.mp4" type="video/mp4" /></video> Note: non-uniform sampling distribution, temporal variation, superposition of spatial scales --- class: center ## Mean Surface Current Speed <img style="width:100%;margin-top:-0.7em;margin-bottom:-0.8em" src="../figures/globalinertialspeed.png"> Formed by binning in latitude and longitude, then averaging. Easy to compute maps of low-order statistics: mean, variance, etc. Does not capture full richness of dataset. What else can be done? -- Strategy: create mathematical/statistical methods to isolate particular processes: (i) turbulent background (ii) eddies. --- class: left ## Perspective Lagrangian trajectories offer a greater potential for studying transport in ocean turbulence than is generally appreciated. This is significant because there is a lot of such data—both surface drifters and RAFOS-type deep floats. This talk will give an overview of two promising methods for rigorous and objective analysis of massive Lagrangian datasets: 1. Stochastic modeling of dispersive Lagrangian trajectories 1. Time/frequency extraction of coherent eddy properties The development of these methodologies involved adapting and advancing modern ideas from both signal processing and statistics. These methods complement more traditional tools such as pseudo-Eulerian maps, dispersion estimates, FTLE's, etc. There are outstanding physical questions that would make promising directions for further research. --- class: center, middle #Part I: Stochastic Modeling --- class: center ##Ocean Turbulence <video preload="auto" width="60%" height="auto" data-setup="{}" autoplay loop controls><source src="../videos/turbulencemovie.mp4" type="video/mp4" /></video> An idealized model of oceanic turbulence by Jeffrey Early. Shading is speed with particle trajectories shown in color. --- class: center ##Modeling Lagrangian Trajectories <video style="margin-top:-0.5em;" preload="auto" width="100%" height="auto" data-setup="{}" autoplay loop controls><source src="../videos/dispersionmovie.mp4" type="video/mp4" /></video> <!-- One of these is 2D turbulence, one is a 3-parameter stochastic model.--> A simple stochastic model for trajectories in oceanic turbulence is the *Matérn process* or <i>damped</i> Fractional Brownian motion. This three-parameter model lets us match (i) the kinetic energy, (ii) the dispersion, and (iii) the degree of small-scale roughness. --- class: left ##The Matérn Process A relatively little-known random process called the Matérn process .cite[(Matérn, 1960)] has a spectrum given by `\[S_{zz}(\omega) = \frac{A^2}{(\omega^2+\lambda^2)^\alpha}\]` where `$A$` sets the energy level, `$\alpha$` sets the high-frequency slope, and `$\lambda$` determines a low-frequency transition to a constant spectral value. -- Recently, we investigated this process in detail, and showed that it is an excellent model for trajectories in ocean turbulence. We speculate it will prove suitable for many other processes as well. We showed the Matérn process is a generalization of *fractional Brownian motion* .cite[(Mandelbrot and Van Ness, 1968)] to include a damping term, `$\lambda$`, that is the essential ingredient for modeling dispersion. .cite[Lilly, Sykulski, Early, and Olhede (2017). Fractional Brownian motion, the Matérn process, and stochastic modeling of turbulent dispersion. <a href="./papers/lilly17-npg.pdf"><nobr>[3.2 Mb pdf]</nobr> </a>] --- class: left ## An Hierarchy of Stochastic Models Let `\(z(t)=u(t)+iv(t)\)` represent a complex-valued <i>velocity</i>. Then we have the following correspondences between stochastic differential equations (SDEs) and velocity spectra `$S_{zz}(\omega)$`. Brownian motion: `\[\frac{dz}{d t} = A\,\frac{d w}{d t} \quad\longrightarrow\quad S_{zz}(\omega) = \frac{A^2}{\omega^2}\]` Fractional Brownian motion (fBm): `\[\frac{d^\alpha z}{d t^\alpha} = A\,\frac{d w}{d t} \quad\longrightarrow\quad S_{zz}(\omega) = \frac{A^2}{\omega^{2\alpha}}\]` Damped Fractional Brownian motion a.k.a. the Matérn process: `\[\left[\frac{d}{d t} + \lambda \right]^\alpha z = A\,\frac{d w}{d t} \quad\longrightarrow\quad S_{zz}(\omega) = \frac{A^2}{(\omega^2+\lambda^2)^\alpha}\]` (These should be written as the corresponding *integral* equations.) --- class: left ## Another Hierarchy of Stochastic Models Damped Brownian motion <i>a.k.a.</i> Ornstein-Uhlenbeck (OU) process: `\[\frac{dz}{d t} + \lambda z = A\,\frac{d w}{d t} \label{maternsde}\quad\longrightarrow\quad S_{zz}(\omega) = \frac{A^2}{\omega^2+\lambda^2}\]` Damped Brownian motion + spin <i>a.k.a.</i> complex OU (cOU) process: `\[\frac{dz}{d t} - i\Omega z+ \lambda z = A\,\frac{d w}{d t} \quad\longrightarrow\quad S_{zz}(\omega) = \frac{A^2}{(\omega-\Omega)^2+\lambda^2}\]` Damped fBm + spin <i>a.k.a.</i> Matérn + spin <i>a.k.a.</i> fractional cOU: `\[\left[\frac{d}{d t} - i\Omega+ \lambda \right]^\alpha z = A\,\frac{d w}{d t} \quad\longrightarrow\quad S_{zz}(\omega) = \frac{A^2}{[(\omega-\Omega)^2+\lambda^2]^\alpha}\]` <!--Note spin is a deterministic effect!--> <!--`\[z(t) = e^{i\Omega t} \tilde z(t)\longrightarrow R_{z z}= e^{i\Omega \tau} R_{\tilde z\tilde z}(\tau) \longrightarrow S_{zz}(\omega)=S_{\tilde z\tilde z}(\omega-\Omega) \]`--> Matérn + spin generalizes cOU to fractional orders, thus unifying OU+cOU+fBm into a single larger family. --- class: left ##Slope and Roughness <img style="width:100%;margin-top:-0.2em;margin-bottom:-0.2em" src="../figures/fbmroughness.png"> The `$\alpha$` parameter sets the slope in fractional Brownian motion, or at high frequencies in the Matérn process, as `$S(\omega)\sim \omega^{-2\alpha}$`. It also sets the degree of small-scale roughness in the trajectories. In the above, we go from slopes of `$\omega^{-3}$` at left to `$\omega^{-1}$` at right. Steeper spectral slope = more smooth. Spectral slope can be formally linked to a mathematical measure of roughness called the *fractal dimension*, though in my opinion, fractal dimension seems to have little practical value. --- class: center ##Slope and Self-Similarity <img style="width:90%;margin-top:-0.2em;margin-bottom:-0.2em" src="../figures/fbmselfsimilarity.png"> The `$\alpha$` parameter also sets *aspect ratio* of statistical self-similarity. --- class: left ##Diffusivity is the Spectrum's Value at Zero The velocity autocovariance and spectrum are defined as `\[R_{zz} (\tau)\equiv\left\langle z(t+\tau)\,z^*(t)\right\rangle = \frac{1}{2\pi}\int_{-\infty}^\infty e^{i \omega \tau} S_{zz}(\omega) \, d\omega\]` The corresponding diffusivity is given by `\[\kappa \equiv \lim_{t\longrightarrow\infty} \frac{1}{4} \frac{d}{d t} \left\langle x^2(t) + y^2(t)\right\rangle =\frac{1}{4} \int_{-\infty}^\infty R_{zz}(\tau)\, d \tau =\frac{1}{4}S_{zz}(0)\]` with the last equality following from `\(\int_{-\infty}^\infty e^{-i \omega \tau} R_{zz}(\tau) \, d\tau= S_{zz}(\omega)\)`. Matérn: `\( \,\,\,\,S_{zz}(\omega) = \frac{A^2}{(\omega^2+\lambda^2)^\alpha} \quad\longrightarrow\quad \kappa = \frac{1}{4}\frac{A^2}{\lambda^{2\alpha}}\)` fBm: `\( \quad\quad S_{zz}(\omega)= \frac{A^2}{\omega^{2\alpha}} \quad\quad\,\,\longrightarrow\quad \kappa = \infty\)` Infinity means the diffusivity increases without bound. Thus the Matérn, unlike fBm, can model diffusive processes. --- class: center ##Comparison of Models <img style="width:65%;margin-top:-0.2em;margin-bottom:-0.2em" src="../figures/randomwalk_dispersion.png"> Top left: numerical simulations. Top right: Matérn process. Bottom left: white noise. Bottom right: red noise (Brownian). --- class: left ##Argument for Generality We conjecture that the Matérn may be of general usefulness: “Systems are often characterized by a pressure to grow—represented by a forcing—together with some drag or resistance on that growth, represented by a damping. After a sufficiently long time, the forcing and the damping equilibrate and one reaches a bounded state. This leads to the speculation that many time series that are well described as fBm over relatively short timescales may be better matched by the Matérn process over longer timescales. More generally, the Matérn process adds a third parameter (damping) to the two parameters (amplitude together with spectral slope or the Hurst parameter) of fBm, thus permitting a wider range of spectral forms to be accommodated. It is therefore reasonable to think that the Matérn process could be of broad interest in many areas in which fBm has already proven itself useful.” .cite[Lilly, Sykulski, Early, and Olhede (2017)] --- class: center, ## The Global Surface Drifter Dataset <video preload="auto" width="100%" height="auto" data-setup="{}" autoplay loop controls><source src="../videos/driftermovie.mp4" type="video/mp4" /></video> Now we will apply the Matérn process to the global drifter dataset. --- class: center, ## Spectra of Drifter Velocities <img style="width:45%;margin-top:-0.7em;margin-bottom:-0.8em" src="../figures/highreslatitudespectra.png"> We see a low-frequency maximium and tidal maxima—all vertical lines—and a maximum along a curving line, the inertial frequency. <!--Inertial frequency becomes the vertical line `\(\omega/f_o=1\)`. Tides become curved. Note roughly universal appearance in new coordinates.--> --- class: center, ##A Conceptual Model for the Spectrum <img style="width:65%;margin-top:-0.7em;margin-bottom:-0.8em" src="../figures/schematictransparency2.png"> \begin{equation} S(\omega) = \overset{\mathrm{background}}{\overbrace{\frac{A_b^2}{\left[\omega^2+\lambda_b^2\right]^\alpha}}} +\overset{\mathrm{inertial\,\, oscillations}}{\overbrace{\frac{A_o^2}{\left(\omega-f_o\right)^2+\lambda_o^2}}}+ \overset{\mathrm{semidiurnal\,\,tide}}{\overbrace{\frac{A_s^2}{\left(\omega-f_s\right)^2+\lambda_s^2}}} \end{equation} This is just several versions of the Matérn spectrum added together. --- class: left ## Details of the Stochastic Modeling We have developed a method for stochastic modeling of trajectories in ocean turbulence, and inferring parameters from large datasets. This is an example of *parametric spectral analysis*: rather than trying to estimate the spectrum directly, one begins with a *conceptual model*, controlled by some parameters, and then finds the choice of parameters which provide the best match to the data. The full method consists of several ingredients: * Create an appropriate stochastic model for particle trajectories. .cite[Sykulski, Olhede, Lilly, and Danioux (2015)] * Employ the Matérn process as a building block. .cite[Lilly, Sykulski, Early, and Olhede (2017)] * Parameter estimation is best done in the frequency domain. .cite[Whittle (1953)] * Parameter estimation must be adjusted to handle *complex-valued* time series. .cite[Sykulski, Olhede, Lilly, and Early (2017)] * Parameter estimation must be corrected for subtle bias due to small sample size effects. .cite[Sykulski, Olhede, Guillaumin, Lilly, and Early (2018), in review.] --- class: center ## Spectra in Latitude Bands <img style="width:75%;margin-top:-0.7em;margin-bottom:-0.8em" src="../figures/drifterspectralinestransparency.png"> Averages in `\(\pm 5 ^\circ \)` bins. Frequency is nondimensionalized as `\(\omega/f_o\)`. Inertial peak broadens equatorward. --- class: center ## Comparison with Spectral Fit <img style="width:75%;margin-top:-0.7em;margin-bottom:-0.8em" src="../figures/drifterspectralinestransparencywithfullfit.png"> In general, spectral model compares very well with estimated spectra. Inertial and semidiurnal peaks are well captured. --- class: center ## Fit to Background vs. Cyclonic Side <img style="width:75%;margin-top:-0.7em;margin-bottom:-0.8em" src="../figures/drifterspectralinestransparencycyclonicwithfit.png"> “Background” fit to anticyclonic (inertial side), matches well the spectrum on the cyclonic (anti-inertial) side, apart from tides. --- class: center ##Application to the Global Drifter Dataset <img style="width:100%;margin-top:-0.2em;margin-bottom:-0.2em" src="../figures/globalantiinertiallambda.png"> The nondimensionalized damping parameter from a parametric fit to the Matérn, only on the anti-inertial side. Shorter ‘memory’ in more energetic regions. We are not aware of any theory for this. --- class: center ##Application to the Global Drifter Dataset <img style="width:100%;margin-top:-0.2em;margin-bottom:-0.2em" src="../figures/globalantiinertialalpha.png"> The first global map of the slope parameter `\(\alpha\)`, with `\(|\omega|^{-2\alpha}\)`. Slopes vary from `\(|\omega|^{-1}\)` to `\(|\omega|^{-3}\)`, with `\(|\omega|^{-2}\)` over major currents. Does not fit the conventional wisdom of only `\(|\omega|^{-2n}\)` slopes, as in e.g. .cite[Berloff and McWilliams (2002b)] <!--Lagrangian integral timescale: `\[T \equiv \frac{1}{2}\frac{\int_{-\infty}^\infty R_{zz}(\tau) \,d \tau}{R_{zz}(0)} = \frac{1}{2}\frac{S_{zz}(0)}{\frac{1}{2\pi} \int_{-\infty}^\infty S_{zz}(\omega) \,d \omega} =\frac{2\kappa}{\sigma^2}\]`--> <!--Spectral slope / fractal dimension: `\[R_{zz}(\tau) \sim |\tau|^{2\alpha-1}, \quad\tau \approx 0 \quad \longleftrightarrow\quad S_{zz}(\omega) \sim |\omega|^{-2\alpha},\quad \omega\gg 0\]`--> <!--which holds for $1/2<\alpha<3/2$; see \citet{goff88-jgr},--> <!-- $\sigma^2 =c_\alpha A^2/\lambda^{2\alpha-1}$ and $\kappa = \frac{1}{4}A^2/(\Omega^2+\lambda^2)^\alpha$, where $c_\alpha\equiv\Gamma\left(\alpha-\frac{1}{2}\right)/\Gamma(\alpha)/(2\sqrt{\pi})$ with $\Gamma(\cdot)$ being the gamma function, see \citet{lilly16-itit}, while the Lagrangian integral timescale is just~$T=2\kappa/\sigma^2$. --> --- class: left ##Summary of Part I The Matérn process is the simplest random process that can simultaneously reproduce (i) the velocity variance (ii) the diffusivity and (iii) the spectral slope or degree of small-scale roughness. -- This is a promising stochastic model for Lagrangian trajectories, where it is useful for parameteric spectral modeling, as a model for synthetic trajectories, as a null hypothesis in eddy studies. -- The Matérn process is conjectured to be of broad applicability. -- The paper describing the Matérn process is intended to serve as an accessible introduction to stochastic processes in general. .cite[Lilly, Sykulski, Early, and Olhede (2017). <a href="https://doi.org/10.5194/npg-24-481-2017">{Fractional Brownian motion, the Matérn process, and stochastic modeling of turbulent dispersion}</a>, <i>Nonlinear Processes in Geophysics</i>.] -- Code can be found in [{jLab}](http://www.jmlilly.net/jmlsoft.html). --- class: center, middle #Part II: Vortex Extraction --- class: center ##Lagrangian Analysis of a Nonlinear Eddy <video preload="auto" width="95%" height="auto" data-setup="{}" autoplay loop controls><source src="../videos/nonlineareddyparticlemoviewithellipses.mp4" type="video/mp4" /></video> An eddy at 24°N tracked by 512×256=131,072 particles, by J. Early. Color is estimated enclosed vorticity in inferred ellipses. --- class: center ##Estimated vs. Actual Enclosed Vorticity <img style="width:70%;margin-top:-0.2em;margin-bottom:-0.2em" src="../figures/nonlineareddygoodvorticity.png"> Log10 histogram of enclosed vorticity estimated from “good” ridges. --- class: center ##Estimated vs. Actual Normal Strain <img style="width:70%;margin-top:-0.2em;margin-bottom:-0.2em" src="../figures/nonlineareddygoodnormalstrain.png"> Log10 histogram of estimated enclosed normal strain. --- class: left ##Elements of Lagrangian Vortex Extraction This analyis method consists of three steps: 1. *Extracting* oscillatory velocity features using a technique called *wavelet ridge analysis*. .cite[Lilly, Scott, and Olhede (2011)] 2. *Thresholding* with a suitable error quantity to remove false positives. .cite[Lilly and Olhede (2012a)] 3. *Linking* ellipse properties to spatially-integrated properties using an extended version of Stokes' theorem. .cite[Lilly (2018)] --- class: left ## A Model for a Trajectory with Eddies First, we conceptualize the orbit of a particle trapped in an eddy, `$\mathbf{x}_o=[x_o(t)\,\, y_o(t)]^T$` or `$z_o(t)=x_o(t) + \mathrm{i} y_o(t)$`, as a <i>time-varying ellipse.</i> `\[z_o(t)=x_o(t)+\mathrm{i}y_o(t)=\mathrm{e}^{\mathrm{i}\theta(t)}\left[a(t)\cos\phi(t) + \mathrm{i} b(t)\sin \phi(t)\right]\]` <img style="width:100%;margin-top:-0.2em;margin-bottom:-0.2em" src="../figures/ellipsestability.png"> We construct a type of “filter” that can isolate such signals. This lets us approximately split a time series `\(z(t)=z_o(t)+z_\epsilon(t)\)` into an oscillatory `\(z_o(t)\)` and stochastic `\(z_\epsilon(t)\)` portion. --- class: left ##Recovering Modulated Oscillations The wavelet transform is a tool for <i>recovering</i> or <i> estimating</i> the properties of modulated oscillation immersed in background noise. The idea is to project the time series onto an oscillatory test signal, or <i>wavelet</i>, to find a “best fit” frequencies at each moment. The fits are then chained together into a continuous curve called a <i>ridge.</i> <!--This uses <i>analytic wavelets</i>, complex-valued wavelets that combine a bandpass operation together with the analytic transform. --> For a vector-valued signal `\(\mathbf{x}_o(t)=\begin{bmatrix}x_1(t) & x_2(t) & \cdots & x_N(t)\end{bmatrix}^T\)` in noise `\(\mathbf{x}_\epsilon(t)\)`, `\(\mathbf{x}(t)=\mathbf{x}_o(t)+\mathbf{x}_\epsilon(t)\)`, define the wavelet transform `\[\mathbf{w}_\mathbf{x}(t ,s) \equiv \int_{-\infty}^{\infty} \frac{1}{s} \psi^*\left(\frac{\tau-t}{s}\right)\,\mathbf{x}(\tau)\,\mathrm{d} \tau\]` (also a vector) and then find the <i>wavelet ridges</i> `\(s(t)\)` from `\[\frac{\partial}{\partial s}\, \left\|\mathbf{w}_\mathbf{x}(t ,s)\right\| = 0,\quad\quad \frac{\partial^2}{\partial s^ 2}\, \left\|\mathbf{w}_\mathbf{x}(t ,s)\right\| < 0.\]` The oscillation is estimated simply by `\(\widehat{\mathbf{x}_o}(t)\equiv\Re\left\{\mathbf{w}_\mathbf{x}(t,s(t))\right\}.\)` .cite[See Delprat et al. (1992), Lilly and Olhede (2009, 2010, 2012a). ] --- class: center ##Example of Vortex Extraction <img style="width:100%" src="../figures/bandwidth_example.png"> The <i>ridge</i> is the curve made by tracing out the maximum modulus of the wavelet transform. --- class: center ##Example of Vortex Extraction <img style="width:100%" src="../figures/bandwidth_looper.png"> The transform value along this curve is an estimate of the oscillatory signal component, visualized as ellipses. Ridge analysis separates modulated elliptical signal from low-frequency meandering and higher-frequency variability. Physically the ellipses quantify the properties of oceanic vortices. --- class: center ##Extraction of Oscillatory Signals <img style="width:100%;margin-top:-0.2em;margin-bottom:-0.2em" src="../figures/nonlineareddyallridges.png"> All trajectory segments that contain a quasi-oscillatory signal. `\[z(t)=x(t) + \mathrm{i} y(t) =z_\epsilon(t) + \boxed{z_o(t)}\]` `$z_\epsilon(t) =$` background, `$z_o(t) =$` oscillatory (e.g. eddy) At each moment, the oscillatory signal is represented by an ellipse. This ellipse is *assigned to* the particle by the analysis method. It is a material ring that the particle is inferred to belong to. --- class: left ##Analysis of Modulated Oscillations It often happens that one wishes to analyze a signal of the form `\[x(t)=a(t)\cos \phi(t) \]` which is known as a *modulated oscillation* or *AM/FM* or *quasi-periodic* signal. An extremely powerful method, the *analytic signal* lets an *instantaneous* amplitude, frequency, and bandwidth be uniquely associated with a given modulated oscillation. These ‘instantaneous moments’ decompose the frequency-domain moments of the signal's Fourier spectrum across time. For example, the power-weighted time average of the instantaneous frequency is the mean frequency of the spectrum. This theory is primarily due to .cite[Gabor (1946)], with important contributions from .cite[Bedrosian (1963)], .cite[Vakman (1996)], .cite[Picinbono (1997)], and others. --- class: left ##The Analytic Signal Method The instantenous moments are defined by using the <i>analytic signal</i> `\[x_+(t)\equiv x(t) +\mathrm{i}\mathcal{H} x(t) \equiv x(t)+\mathrm{i}\frac{1}{\pi}\,\int_{-\infty}^\infty\frac{x(u)}{t-u}\,\mathrm{d} u\]` with `\(\mathcal{H}\)` being the Hilbert transform. The amplitude is uniquely defined as `$a(t)=|x_+(t)|$` and frequency as `$\omega(t) = \frac{d}{dt}\arg x_+(t)$`. **Uniqueness:** If you wish to assign the instantaneous moments using a <i>linear filter</i>, and if you wish your method to recover the correct results for a strictly periodic signal `$x(t)=a_o\cos(\omega_o t)$`—a property called *harmonic correspondence*— there is no alternative to using the analytic signal .cite[(Vakman, 1996)]. **Recoverability:** With `$x(t)$` generated by `$x(t)=a(t)\cos \phi(t)$`, the generating amplitude `$a(t)$` and phase `$\phi(t)$` are exactly recoverable from `$x(t)$` provided `$a(t)$` is supported on strictly lower frequencies than is `$e^{i\phi(t)}$` .cite[(Bedrosian, 1963)]. --- class: left ##Generalization to 2D and 3D Instantaneous moments generalize to multivariate oscillations. A bivariate modulated oscillation `$\mathbf{x}=[x(t)\,\, y(t)]^T$` is the same as a time-varying ellipse .cite[(Lilly and Olhede, 2010a)]. <img style="width:100%;margin-top:-0.6em;margin-bottom:-0.2em" src="../figures/ellipsestability.png"> A trivariate modulated oscillation `$\mathbf{x}=[x(t)\,\, y(t)\,\, z(t)]^T$` is equivalent to a modulated ellipse within a rotating, tilting plane .cite[(Lilly, 2011)]. The time-varying ellipse properties may be recovered from the vector-valued analytic signal `$\mathbf{x}_+(t)\equiv\mathbf{x}(t)+\mathrm{i}\mathcal{H}\mathbf{x}(t)$`. --- class: left ##Recoverability Condition in 2D Q. Under what conditions are the time-varying properties of an ellipse exactly recoverable from the trajectory of one particle? A. When the ellipse geometry evolves more slowly than circulation of particles around the periphery. <img style="width:100%;margin-top:-0.7em;margin-bottom:-0.5em" src="../figures/bedrosian-ellipses-all.png"> Ellipse parameters for top half are exactly recoverable. --- class: left ##Modulated Oscillations in Noise With noise present, we write `$\mathbf{x}(t)=\mathbf{x}_o(t)+\mathbf{x}_\epsilon(t)$`, where `$\mathbf{x}_o(t)$` is an oscillation and `$\mathbf{x}_\epsilon(t)$` is noise. In this case, applying the analytic signal method directly fails. The wavelet transform can be used as basis for <i>recovering</i> or <i> estimating</i> the properties of modulated oscillation immersed in background noise. Note that the analytic wavelet transform implicitly combines filtering, and creation of the analytic signal, into one step. *Wavelet ridge analysis*, due to .cite[Delprat (1992)] and investigated further by .cite[Mallat (1999)], allows the recovery of the oscillation. The idea is to project the time series onto an oscillatory test signal, or <i>wavelet</i>, to find a “best fit” frequencies at each moment. The fits are then chained together into a continuous curve called a <i>ridge.</i> This method lets us strip off the noise, resulting in an estimate of the analytic part of the oscillatory portion of the signal, `$\mathbf{x}_o(t)$`. <!-- *Optimal* in that separation is done as well as is possible subject to the constraints of the uncertainty principle (time vs. frequency).--> --- class: left ##A Bias Estimate The wavelet ridge method has two sources of error, *variance* arising form the noise, and *bias* arising from the degree of modulation of the oscillatory signal. An internal estimate of the bias, formed for an unknown modulated oscillation from the wavelet transform itself, was derived by .cite[Lilly and Olhede (2010b)]. The first-order departure from a pure sinusoid is called the *bandwidth*, which has terms like `$a'(t)/a(t)$` (amplitude modulation). The *second-order* departure has terms like `$a''(t)/a(t)$`, which we call the *curvature*. This turns out to control the bias of the estimate. Wavelet ridge analysis was generalized to handle modulated multivariate oscillations by .cite[Lilly and Olhede (2010a)]. A bias estimate for the multivariate case was derived by .cite[Lilly and Olhede (2012)], and is a multivariate analogue of the curvature. Thresholding on the bias is critical for removing false positives! --- class: center ##An Extended Version of Stokes' Theorem We can connect the time-varying properities of an ellipse, to the integral properties within the ellipse, using Stokes' theorem. `\[\oint_C \mathbf{u}^T\mathrm{d}\mathbf{x} = \iint_A \mathbf{k} \cdot \nabla \times \mathbf{u}\,\mathrm{d} A \]` -- `\[\oint_C \mathbf{u}\,\mathrm{d}\mathbf{x}^T = \iint_A \nabla \mathbf{u}\,\mathrm{d} A \begin{bmatrix} 0 & 1 \\ -1 & 0\end{bmatrix} \]` Contains four relations between contour integrals of velocity (left) and spatial integrals of velocity gradient (right): (i) Kelvin-Stokes theorem, (ii) divergence theorem, and similar results for (iii) normal strain and (iv) shear strain. An alternate presentation: `\[\oint_C \mathbf{u}\,\mathrm{d}\mathbf{n}^T = \iint_A \nabla \mathbf{u}\,\mathrm{d} A \]` Lilly (2018), <a href="http://www.mdpi.com/2311-5521/3/1/16">{Kinematics of a fluid ellipse in a linear flow}</a>, <i>Fluids</i>. --- class: left ##Ellipse / Flow Equivalence The principle of *ellipse / flow equivalence* states that, for any ellipse of fluid particles marked out within a linear flow: * (i) The linear flow determines ellipse rates of change * (ii) The ellipse rates of change imply a unique linear flow The ellipse rates of change include three geometry parameters (size, shape, orientation), plus the flow of particles around the periphery. This builds on the well known result that ‘a linear flow maps an ellipse into another ellipse’, e.g. .cite[(Kida, 1981)]. Ellipse/flow equivalence can be shown to be a consequence of the extended Stokes' theorem. Thus, if we know the rates of change of a fluid ellipse, we know the spatially averaged vorticity, divergence, and strain terms within it. This implies a deep connection between Lagrangian and Eulerian perspectives, suggesting untapped potential from Lagrangian data. .cite[Lilly (2018), <a href="http://www.mdpi.com/2311-5521/3/1/16">{Kinematics of a fluid ellipse in a linear flow}</a>, <i>Fluids</i>.] --- class: center ##Lagrangian Analysis of a Nonlinear Eddy <video preload="auto" width="95%" height="auto" data-setup="{}" autoplay loop controls><source src="../videos/nonlineareddyparticlemoviewithellipses.mp4" type="video/mp4" /></video> An eddy at 24°N tracked by 512×256=131,072 particles, by J. Early. Color is estimated enclosed vorticity in inferred ellipses. --- class: center ##Estimated vs. Actual Enclosed Vorticity <img style="width:70%;margin-top:-0.2em;margin-bottom:-0.2em" src="../figures/nonlineareddyallvorticity.png"> Computed for all inferred ellipses, prior to thresholding. --- class: center ##Estimated vs. Actual Normal Strain <img style="width:70%;margin-top:-0.2em;margin-bottom:-0.2em" src="../figures/nonlineareddyallnormalstrain.png"> Computed for all inferred ellipses, prior to thresholding. --- class: center ##Distribution on Radius / Velocity Plane <img style="width:72.5%;margin-top:-0.2em;margin-bottom:-0.2em" src="../figures/nonlineareddyallrvdist.png"> Distribution of radius and velocity of all inferred ellipses. --- class: center ##Estimated Bias on Radius / Velocity Plane <img style="width:80%;margin-top:-0.2em;margin-bottom:-0.2em" src="../figures/nonlineareddyallrverr.png"> Bias is estimated within each recovered oscillation using a curvature condition, equation (62) of .cite[Lilly and Olhede (2012)]. --- class: center ##Vorticity Error on Radius / Velocity Plane <img style="width:80%;margin-top:-0.2em;margin-bottom:-0.2em" src="../figures/nonlineareddyallrvzetaerr.png"> Actual error between estimated enclosed vorticity, and actual enlosed vorticity, normalizd by the actual enclosed vorticity. --- class: center ##Distribution After Threshholding <img style="width:72.5%;margin-top:-0.2em;margin-bottom:-0.2em" src="../figures/nonlineareddygoodrvdist.png"> Including only those portions of inferred signals that fall below the value 1/15`$\approx 6.6\%$` for the estimated bias. --- class: center ##Effect of Thresholding on Vorticity Error <img style="width:47.5%;margin-top:-0.2em;margin-bottom:-0.2em" src="../figures/nonlineareddyallvorticity.png"> <img style="width:47.5%;margin-top:-0.2em;margin-bottom:-0.2em" src="../figures/nonlineareddygoodvorticity.png"> Including only those portions of inferred signals that fall below the value 1/15`$\approx 6.6\%$` for the estimated bias. --- class: center ##Effect of Thresholding on Strain Error <img style="width:47.5%;margin-top:-0.2em;margin-bottom:-0.2em" src="../figures/nonlineareddyallnormalstrain.png"> <img style="width:47.5%;margin-top:-0.2em;margin-bottom:-0.2em" src="../figures/nonlineareddygoodnormalstrain.png"> Including only those portions of inferred signals that fall below the value 1/15`$\approx 6.6\%$` for the estimated bias. --- class: center ##Evolution on the Radius/Velocity Plane <video preload="auto" width="70%" height="auto" data-setup="{}" autoplay loop controls><source src="../videos/nonlineareddyrvmovie.mp4" type="video/mp4" /></video> Color is vorticity error. Instantaneous radial profile is shown. --- class: center ##Interpretation of Low-Frequency Band <img style="width:47.5%;margin-top:-0.2em;margin-bottom:-0.2em" src="../figures/nonlineareddyallrverr.png"> <img style="width:47.5%;margin-top:-0.2em;margin-bottom:-0.2em" src="../figures/nonlineareddyallrvduration.png"> The right-hand side shows duration in units of the period. Short-duration signals are difficult to estimate due to edge effects. The low-frequency band is associated with signals that are so brief that edge effects dominate their entire duration. This suggests that large bias is due to edge effects associated with entrainment and detrainment. --- class: center ##Physical Meaning of Estimated Bias <img style="width:47.5%;margin-top:-0.2em;margin-bottom:-0.2em" src="../figures/nonlineareddyallrverr.png"> <img style="width:47.5%;margin-top:-0.2em;margin-bottom:-0.2em" src="../figures/nonlineareddyallrvspeedstd.png"> The estimated bias is a property of an inferred or estimated oscillatory signal. What does this correspond to physically? <!--There are a bunch of things is does *not* correspond to.--> Very good agreement is found with the Lagrangian average speed deviation, that is, `$\sqrt{\langle V-\langle V\rangle\rangle}$`, normalized by `$\langle V\rangle$`. Here `$\langle \cdot\rangle$` is an average over the actual particle trajectory during `$\pm$` one half-period of the inferred signal. --- class: center ##Particle Trapping in a Nonlinear Eddy <video preload="auto" width="70%" height="auto" data-setup="{}" autoplay loop controls><source src="../videos/nonlineareddytrappingmovie.mp4" type="video/mp4" /></video> Color is vorticity gradient strength. Black = stable, pink = unstable. --- class: left ## Three Types of Particles Combining vortex signal extraction with classical analysis, one finds particles may be classified into one of three types at each moment. 1. *Trapped and Stable:* These particles are moving with the vortex, and can be used to accurately infer vortex properties. 1. *Trapped and Unstable:* These particles are also moving with the vortex, but are too variable to accurately infer vortex properties. 1. *Untrapped and Unstable:* These particles do not move with the vortex, and vortex properties cannot be inferred from them—but oscillatory signals linked to the vortex are still detected. All trapped particles are either stable or unstable ridge points, i.e. there are no non-oscillatory trapped particles. Untrapped and *stable* particles occur very rarely and are short-lived. The distinction between trapping and recoverability is important. It means not all particles trapped in vortices contain useful information about the vortex (when considered in isolation). An implication of this is that not all vortices are equally observable. Highly unstable or variable vortices (e.g. rapidly precessing ellipses) are more difficult to recover information about. --- class: center ##Composite Eddy in First and Second Half <img style="width:47.5%;margin-top:-0.2em;margin-bottom:-0.2em" src="../figures/nonlineareddyspeedofcompositefirsthalf.png"> <img style="width:47.5%;margin-top:-0.2em;margin-bottom:-0.2em" src="../figures/nonlineareddyspeedofcompositesecondhalf.png"> The eddy evolves to a smaller size from the first half to the second. --- class: center ##Vorticity Error in First and Second Half <img style="width:47.5%;margin-top:-0.2em;margin-bottom:-0.2em" src="../figures/nonlineareddyzetaerrfirsthalf.png"> <img style="width:47.5%;margin-top:-0.2em;margin-bottom:-0.2em" src="../figures/nonlineareddyzetaerrsecondhalf.png"> Error is controlled by the zero vorticity line in the second half only. Particles in the trapping region, but outside the zero vorticity line, are able to accurately estimate vorticity during the first half only. This represents a qualitative change in the boundary between stable and unstable trapping, arising from the onset of vortex instability. Distinction between an entrainment region and a stirring region? --- class: center ## Application to the Global Drifter Dataset <img style="width:100%" src="../figures/highresspaghetti.png"> Apply multivariate ridge detection algorithm to 20K time series. --- class: center ### Ridge Detection in the Global Drifter Dataset <img style="width:90%" src="../figures/alleddies.png"> Red = cyclonic rotation, blue = anticyclonic. --- class: center ### Removing False Positives <img style="width:90%" src="../figures/significanteddies.png"> After comparison with a null hypothesis of red noise, employing the Matérn process of .cite[Lilly, Sykulski, Early, and Olhede (2017)]. --- class: left ## Summary of Part II These results represent progress in bridging the gap between Lagrangian data analysis and vortex dynamics. 3. An extended version of Stokes' theorem allows the estimated rate of change of an ellipse to be converted to *spatial averages*, for all four components of the velocity gradient tensor. 1. The application of a bias, or error, metric allows enclosed vorticity to be estimated with very high accuracy, and enclosed strain to be estimated fairly well also. 2. In this approach, we asses the analysis method with respect to elliptical curves *that the particles believe they are observing.* 1. Application to a model of a nonlinear eddy shows very good agreement with classical particle trapping analysis. Further work: application to other simple systems; idealized vortex solutions; and global datasets. A very interesting direction would be to connect the single-particle Lagrangian signal theory to the Lagrangian Coherent Structures. Nonlinear mesoscale eddies are not well understood at large times! --- class: center, middle # Part III: Element Analysis --- class: left ## Motivation The ridge analysis is intended for signals consisting of modulated oscillations in noise. Another type of signal that may occur is short-duration, isolated ‘bursts’ or ‘impulses’—that is, events which may be represented as a wavelet or the temporal integral of a wavelet. This suggests a model for a real, univariate signal of the form `\begin{equation}\label{signalmodel} x(t) = \sum_{n=1}^N \Re\left\{c_n \psi\left(\frac{t-t_n}{\rho_n}\right)\right\} +x_\epsilon(t) \end{equation}` consisting of superpositions of amplitude scaled, stretched, and phase-shifted versions of some basis signal `$\psi(t)$`, called the *element*. It is assumed that the different realizations of the element are sufficiently separated in time and frequency such that they do not interfere, in a way that will be made precise later. --- class: left ## Ingredients Element analysis consists of several steps: 1. Choice of element 1. Event detection 2. Rejection of statistically insignificant events 3. Rejection of non-isolated events Developed in .cite[Lilly, J. M. (2017). Element analysis: a wavelet-based method for analyzing time-localized events in noisy time series. Proceedings of the Royal Society of London, Series A, 473 (2200): 20160776, 1–28. [[10.7 Mb Pdf]](http://rspa.royalsocietypublishing.org/content/473/2200/20160776)] --- class: center ## Localized Wavelets <img style="width:100%;margin-top:-0.7em;margin-bottom:-0.5em" src="../figures/timedomainwavelets.jpg"> The method is suitable for events that themselves resemble wavelets, or the temporal integral of a wavelet (the `$\beta=0$` case). --- class: left ## Signal Model Our signal model, using the generalized Morse wavelets, is `\begin{equation}\label{morseelementmodel} x(t) = \sum_{n=1}^N \Re \left\{ c_n \psi_{\mu,\gamma}\left(\frac{t-t_n}{\rho_n}\right)\right\} +x_\epsilon(t) \end{equation}` where `$\mu$` is the most suitable value of the order or `$\beta$` parameter. The `$n$`th event is then characterized by coefficient `$c_n$`, time `$t_n$`, and scale `$\rho_n$`. We wish to *estimate* these unknown quantities. We then take the wavelet transform with an order `$\beta$` wavelet in the same `$\gamma$` family, leading to `\begin{equation}\label{transformofelementmodel} w_{\beta,\gamma}(\tau,s)=\frac{1}{2}\sum_{n=1}^N c_n\int_{-\infty}^{\infty} \frac{1}{s} \psi_{\beta,\gamma}^*\left(\frac{t-\tau}{s}\right)\psi_{\mu,\gamma}\left(\frac{t-t_n}{\rho_n}\right)\,d t+\varepsilon_{\beta,\gamma}(\tau,s). \end{equation}` Owing to the properties of the generalized Morse wavelets, the integral in the above equation has a closed-form expression in terms of a rescaled `$\psi_{\beta+\mu,\gamma}(t)$` wavelet. --- class: left ## Event Detection We then find *isolated maxima* of the transform, that is, time-scale `$(\tau,s)$` points at which `\begin{multline} \quad\quad\quad\quad\frac{\partial}{\partial \tau}\left|w_{\beta,\gamma}(\tau,s)\right| =\frac{\partial}{\partial s}\left|w_{\beta,\gamma}(\tau,s)\right|=0, \\\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\frac{\partial^2}{\partial \tau^2}\left|w_{\beta,\gamma}(\tau,s)\right|<0,\quad\frac{\partial^2}{\partial s^2}\left|w_{\beta,\gamma}(\tau,s)\right|<0.\label{maxconditions} \end{multline}` Because we have a simple expression for `$w_{\beta,\gamma}(\tau,s)$`, we can relate the event properties `$c_n$` `$t_n$`, and `$\rho_n$` to the value of `$w_{\beta,\gamma}(\tau,s)$` at a maxima. This lets us work backwards from the maxima points to estimates of the event properties. --- class: center ##An Example of Events in White Noise <img style="width:100%;margin-top:-0.7em;margin-bottom:-0.5em" src="../figures/impulses-morsetrain-whitenoisy.jpg"> The events are based on the `$\psi_{2,2}(t)$` wavelet in this case. Grey dots are all maxima, black are significant and isolated. --- class: left ##Assessing Statistical Significance This is the key step for large-scale applications. The obvious approach is to simulate many noise realizations suitable for the time series of interest, and then apply the method. This works great for one time series but not for 20,000. Fortunately, there is a workaround. It turns out that all we need to know is the statistics of a 5-vector, containing this point, the next point, the previous point, the point above, and the point below. Actually, given the assumption of power-law Gaussian noise, having a spectrum of the form `$A^2/\omega^{2\alpha}$`, the covariance statistics of this 5-vector are known analytically. This means that all one needs to do is simulation realizations of this 5-vector in order to determine the probability that a point at a given scale will be a maximum of a certain magnitude, thus establishing statistical significance. --- class: center ##Event Distributions for the Example <img style="width:100%;margin-top:-0.7em;margin-bottom:-0.5em" src="../figures/impulses-examplenoisedist.jpg"> Color = 5-vector realizations, gray dots = all detected events Gray squares = values of true events Black circles = estimated values of significant events. Lines are significance curves for event detection rates. <!--<img style="width:100%;margin-top:-0.7em;margin-bottom:-0.5em" src="../figures/impulses-noisedist-22.jpg">--> --- class: left ##Assessing Isolation This is done through the identification of *regions of influence* around isolated maxima, related to but distinct from concentration regions. Analytic expressions for these regions of influence may be found. An event is rejected if it is within the region of influence (subject to a threshold) of another, larger-amplitude event. This turns out to be mostly important for rejecting artifacts due to closely-spaced events. --- class: center ## A Real-World Application <img style="width:100%;margin-top:-0.7em;margin-bottom:-0.5em" src="../figures/tpjaos_tracks.jpg"> Measurement tracks for a sea surface height satellite. One of the main windows into the ocean circulation. --- class: center ## A Real-World Application <img style="width:100%;margin-top:-0.7em;margin-bottom:-0.5em" src="../figures/impulses-labseathreepanel.jpg"> From an short altimeter track segment crossing the Labrador Sea. --- class: center, middle # Thank you! .center[This talk is available at].center[[http://www.jmlilly.net/talks/lilly18-awi.html](http://www.jmlilly.net/talks/lilly18-awi.html)] .center[.footnote[P.S. Like the way this presentation looks? Check out [{Liminal}](https://github.com/jonathanlilly/liminal).]] --- class: center ##Comparison with the Naive Estimator <img style="width:45%;margin-top:-0.2em;margin-bottom:-0.2em" src="../figures/nonlineareddygoodrvdist.png"> <img style="width:45%;margin-top:-0.2em;margin-bottom:-0.2em" src="../figures/nonlineareddyallpointsrvdist.png"> Right-hand side is estimated from *all* Lagrangian data points by locating vortex center, then computing azimuthal velocity. --- class: center ##Comparison with the Naive Estimator <img style="width:45%;margin-top:-0.2em;margin-bottom:-0.2em" src="../figures/nonlineareddygoodvorticity.png"> <img style="width:45%;margin-top:-0.2em;margin-bottom:-0.2em" src="../figures/nonlineareddyallpointsvorticity.png"> As on the previous slide and assuming solid-body rotation. Right-hand side gives estimated vorticity within a circle.