# Deep Lagrangian connectivity in the global ocean inferred from Argo floats

Ryan Abernathey  Lamont Doherty Earth Observatory of Columbia University, New
York, USA  Christopher Bladwell  School of Mathematics and Statistics,
University of New South Wales, Sydney NSW 2052, Australia  Gary Froyland
School of Mathematics and Statistics, University of New South Wales, Sydney
NSW 2052, Australia  Konstantinos Sakellariou  School of Mathematics and
Statistics, University of New South Wales, Sydney NSW 2052, Australia  Nodes &
Links Ltd, Leof. Athalassas 176, Strovolos, Nicosia, 2025, Cyprus

###### Abstract

We describe the application of a new technique from nonlinear dynamical
systems to infer the Lagrangian connectivity of the deep global ocean. We
approximate the dynamic Laplacian using Argo trajectories from January 2011 to
January 2017 and extract the eight dominant coherent (or dynamically self-
connected) regions at 1500m depth. Our approach overcomes issues such as
sparsity of observed data, and floats continually leaving and entering the
dataset; only 10% of floats record for the full six years. The identified
coherent regions maximally trap water within them over the six-year time
frame, providing a distinct analysis of the deep global ocean, and relevant
information for planning future float deployment. While our study is concerned
with ocean circulation at a multi-year, global scale, the dynamic Laplacian
approach may be applied at any temporal or spatial scale to identify coherent
structures in ocean flow from positional time series information arising from
observations or models.

##  1 Introduction

The Argo program has transformed oceanography by vastly increasing the
spatiotemporal coverage of ocean observations [RT09, RFR+16, JRZ+17]. The most
widely used aspect of Argo data are the temperature / salinity (T/S) profiles,
which can provide, for example, estimates of ocean heat content e.g. [RCG+15,
RFR+16, CTF+17] and constraints for ocean data assimilation systems e.g.
[FCH+15]. However, Argo floats also provide a secondary type of data: their
drift velocity at 1000 dbar over their 10-day cycle. The aim of this paper is
to use the Argo drift velocities to probe the structure of the mid-depth
lateral circulation using new analytical tools from dynamical systems theory.

Compared to T/S profiles, Argo drift velocities have been relatively under-
utilized. Lebedev  _et al._ [LYMH07] produced the first dataset of global
drift velocities, and Katsumata and Yoshinari [KY10] used this data to provide
an esimate of the mean flow field at 1000 dbar. Gray and Riser [GR14] used
individual Argo drift velocity vectors as a reference level to derive,
together with the Argo geostrophic shear, a climatology of absolute
geostrophic velocities over the upper 2000 dbar. Argo drift velocities have
also been used to study mesoscale processes. Chiswell [Chi13] used these
velocities to make estimates of eddy diffusivity based on the deep Eulerian
velocity decorrelation timescale.

Roach  _et al._ [RBS16] made an important conceptual breakthrough by
considering long-term 1000 dbar Lagrangian trajectories; by analyzing
simulated Argo floats in an ocean model, they concluded that the 10-day
profiling cycle, in which the floats return to the surface to transmit their
data, does not significantly disrupt the Lagrangian trajectory at 1000 dbar.
In other words, simulated floats with profiling cycles had very similar
trajectories and dispersion characteristics to those that did not profile but
just remained at 1000 dbar. This means that Argo floats can be used not only
to study 10-day displacement vectors but also trajectories over much longer
timescales. This property was further exploited by Roach  _et al._ [RBS16],
and subsequent studies by Balwada  _et al._ [BSL+16, BLSF21], to calculate
single and multi-particle dispersion statistics. The notion that Argo
trajectories may represent unbroken water mass trajectories over long time
intervals opens the door to many tools from dynamical systems [Dij05] and
Lagrangian [Wig05, VSGA+18] analysis.

Over the last two decades a prominent subset of Lagrangian methods have aimed
at identifying coherent sets and structures in fluids, with ocean flow as an
important motivating testbed. These objects have influence on for example: the
motion of surface drifters [OBVH+13], the accumulation of plastic on the ocean
surface [FSvS14], and the risk of bycatch [SHJ+18]. One prevalent class of
Lagrangian methods uses linearisations of the flow to identify curves (e.g. on
the surface ocean) that are distinguished in some way: usually as local
extremisers of local spatial expansion, possibly constrained tangentially or
transversally to the curve. These include finite-time Lyapunov exponents
(FTLEs) [Pie94, HY00, WKC12], and variational theory [Hal11, HBV13], to name
just a few of these approaches. The vast majority of these techniques rely on
accurate spatial derivatives of the flow, which are impossible to obtain from
Argo drifter observations because of the latter’s sparsity.

Instead of using spatial derivatives of the nonlinear flow, transfer operator
approaches [FSM10, Fro13] create a global linear description of the flow (the
transfer operator) and extract coherent regions from singular vectors of this
operator. Finite-time coherent sets were introduced in [FSM10] as time-
parameterised families of regions in a fluid that have _persistently small
leakage_ from the family over the finite-time evolution. The transfer operator
method for estimating finite-time coherent sets [FSM10] has been applied in
the ocean setting to identify ocean eddies in three dimensions [FHR+12] and
study their surface decay [FHRVS15]. This approach directly estimates the
coherent regions themselves, e.g. the eddy as a two-dimensional or three-
dimensional object on the surface ocean or in the three-dimensional ocean (or
in two-dimensional surfaces or three-dimensional volumes in the atmosphere
[FSM10]). Because the transfer operator method [FSM10] does not require
spatial derivatives, it is robust to noise, and has less stringent
requirements on the resolution of trajectory sampling.

In a move toward sparser data requirements, Froyland and Padberg-Gehle [FPG15]
showed fuzzy $c$-means clustering could extract coherent sets from sparse data
sets, and provided a 5-year coherent set analysis of the surface ocean from
the Global Drifter Array (GLAD), similar to the analysis we will describe in
this work. A diffusion-maps implementation of the transfer operator method
[Fro13] was described in [BK17], who were able to reproduce the surface ocean
results from [FPG15] using the same trajectory dataset.

In the present work we use the dynamic Laplacian [Fro15] approach, and in
particular the finite-element method (FEM) implementation [FJ18]. The dynamic
Laplacian formally arises as a “zero-noise limit” of the transfer operator
approach [Fro13], and therefore one expects very similar estimates of finite-
time coherent sets. In the next section, we discuss an equivalent
characterisation of finite-time coherent sets as a time-parameterised family
of sets whose _boundary remains persistently small_ (i.e. does not become
filamented, and therefore minimises diffusive leakage) under evolution. Both
the transfer operator and dynamic Laplacian approaches are a form of spectral
clustering, adapted to the Lagrangian dynamics we are investigating.

The FEM implementation [FJ18] handles sparse, scattered data well, including
piecing together trajectory segments formed by short drifter trajectories that
do not last for the full time interval of analysis (adapting an approach from
[FPG15]). The FEM approach also has a number of advantages over other
implementations such as [HKTH16, BK17], for instance (i) being unaffected by
varying spatial densities of trajectories, (ii) providing full coherent set
“fields” on the domain (not only at trajectory points), and (iii) not
requiring any parameter tuning. Because of the above points, we use this
FEM–dynamic Laplacian approach in this work to analyse the ocean at the
$\sim$1000 metre level, using only the sparsely distributed Argo float
position information encoded in their trajectories.

Related analysis of _almost-invariant sets_ , namely regions of the ocean that
remain approximately _fixed in space_ under the ocean flow, have been
undertaken using float trajectories in the surface Gulf of Mexico [MBVO+17]
using a combination of drifters from several drifter arrays, and at 1500–2500m
depth in the Gulf of Mexico [MBVO+19] using RAFOS floats. In these studies,
the underlying ocean currents were assumed to be stationary in a statistical
sense, meaning that the likelihood of drifter movements between similar
regions of the ocean at different times did not change over time. This enabled
a huge multiplicative “reuse” of the data (similar to [VSEF12]), and allowed
the more classical Ulam approach for finding almost-invariant sets, e.g.
[FPET07] to be employed.

As far as we are aware the current study is the first truly Lagrangian study
of coherent sets that directly uses deep ocean float trajectories. The FEM
implementation of the dynamic Laplacian [FJ18] approach we use here is ideally
suited to this type of transport analysis when only sparse, scattered
trajectory data is available. This FEM approach has been previously used in
the oceanographic context to track eddy motion using satellite-derived
velocity fields [FJ18, FRS19], as well as to track the meander of the Gulf
Stream from purely kinematic trajectory positional information [FRS19]; these
experiments were at smaller spatial scales and used richer, synthetic
trajectory information compared to the Argo trajectories used here.

Our goals in this paper are twofold. First, we wish to prototype the dynamic
Laplacian methodology using real-world data and demonstrate that it is
computationally practical and gives reasonable results. We hope that this
demonstration will inspire others to adapt this method to different datasets.
Second, we wish to use the results of the analysis to investigate the dynamic
geography of the mid-depth ocean flow, a topic with multiple applications. The
coherent sets identified by this method have significant implications for the
Argo deployment strategy, since they effectively reveal the regions that are
poorly dynamically connected by Argo floats. We can also explore the
relationships between the Argo coherent sets and large-scale ocean tracers,
although, as described below, these results are somewhat inconclusive.

An outline of the paper is as follows. In section 2 we give a brief background
on finite-time coherent sets and define the dynamic Laplace operator. Section
3.1 describes Argo drifter data and the method by which we create trajectories
from the drifter movements. Section 3.2 details the numerical discretisation
of the dynamic Laplacian and how we compute this discretisation using
trajectories from Argo drifters. In section 3.3 we describe how to isolate
individual coherent sets from the eigenfunction of the dynamic Laplacian using
the Sparse EigenBasis Approximation (SEBA) algorithm. Section 4 opens with a
discussion of the eight most coherent regions we identify in the context of
the ocean circulation at 1000m depth. In section 4.1 we describe how to find
precise boundaries for the eight maximally coherent sets. Section 5 compares
the locations of the eight identified global-scale coherent sets with the
large-scale properties planetary potential vorticity and oxygen concentration.
We conclude in section 6 with a summary of our findings, the advantages of the
dynamic Laplacian analysis approach, and the relevance of coherent set
analysis to elucidating ocean circulation properties and to informing the
placement of drifters in programs such as Argo.

##  2 Finite-time coherent set analysis and the dynamic Laplacian

Using trajectories of the Argo float array, we wish to identify coherent
regions in the deep global ocean. The trajectories of the Argo floats at 1000m
depth induce a nonautonomous dynamical system on a two-dimensional phase space
at this depth. Finite-time coherent sets were introduced in [FSM10, Fro13] as
subsets of phase space that _minimally mix with (or leak into) the rest of the
phase space_. Coherent features such as mesoscale eddies [FHR+12, FHRVS15] and
the Antarctic polar vortex [FSM10] have been identified using transfer
operator methods using this least-mixing / least-leaking criterion.

An alternate, but strongly related criterion for finite-time coherence is
based on geometric considerations [Fro15]. In the presence of small diffusive
processes, the extent to which a region will mix with the surrounding phase
space is proportional to average boundary size over the time interval of
interest. Let $X$ denote the 1000m depth level (with land removed); this is a
subset of $\mathbb{R}^{2}$. By $[0,t_{f}]$ we denote the time interval over
which we analyse the float dynamics, where we nominally begin at an initial
time of 0. Let $\phi:[0,t_{f}]\times X\to X$ denote the nonlinear time-
dependent flow map that evolves a location $x\in X$ at time $0$ to its future
location $\phi_{t}(x)$ at time $t$. We note that in our application in this
paper, the domain $X$ is invariant under the flow, i.e. $\phi_{t}(X)=X$, but
this is not necessary for the methodology to be applied. Let $A\subset X$ and
we denote by $\partial A$ and $\partial X$ the boundaries of $A$ and $X$
respectively. We will be concerned with regions $A$ at the 1000m depth level
so that $A\cap\partial X=\emptyset$. Let $\ell(\partial A)$ denote the length
of the boundary of $A$ and $a(A)$ denote the area111If $X$ were three-
dimensional, $\ell(\partial A)$ and $a(\partial A)$ would be replaced with the
area and volume, respectively, of the boundary of $A$. of $A$. Finite-time
coherent sets are sets $A\subset X$ that produce close to minimal values
of222The quantity $h(A)$ is the dynamic Cheeger value of $A$ [Fro15], written
here in the “Dirichlet boundary condition” version; see [FJ18] for details.

$h(A):=\frac{1}{t_{f}}\int_{0}^{t_{f}}\frac{\ell(\partial(\phi_{t}(A)))}{a(\phi_{t}(A))}\
dt.$ (1)

Thus, we seek $A\subset X$ whose time-averaged (under forward evolution of $A$
by $\phi$) ratio of boundary length to enclosed area is minimised. Such a set
will have minimal mixing in the presence of small diffusion.

Figure 1: Illustration of a region of the ocean (light gray), containing a set
$A$ (dark gray), whose boundary does not grow significantly under forward
evolution.

To find these regions, we will use the _dynamic Laplace operator_ [Fro15],
which we now describe. We denote by $\Delta$ the usual Laplace operator acting
on functions $f:X\to\mathbb{R}$ defined on our 1000m depth level $X$. It is
well known that the spectrum and eigenfunctions of $\Delta$ encode geometric
features of the domain $X$ [GN13]. This is related to the question posed in
the 1966 paper of Mark Kac: “Can one hear the shape of a drum?” [Kac66], which
asks if a surface can be classified by the spectrum of the Laplace operator on
that surface. Further, the spectrum and eigenfunctions of $\Delta$ appear in
_isoperimetric theory_ (see for example [Cha84]), which concerns sets with
minimal boundary size relative to volume.

The dynamic Laplace operator extends these ideas to dynamical systems to
answer _dynamic_ isoperimetric questions like the minimal values of (1). Let
$f:X\to\mathbb{R}$ denote a function that we think of as being “anchored” at
time $0$. There is a natural way to push forward this function under the
dynamics, namely by precomposing with the inverse flow map $\phi_{-t}$. That
is, the function $f\circ\phi^{-t}:X\to\mathbb{R}$ is anchored at a time $t$ in
the future and strictly speaking is a function defined on $\phi^{t}(X)$, but
in our deep global ocean application, we have $\phi^{t}(X)=X$ for all
$t\in[0,t_{f}]$. To access the geometry of $\phi^{t}(X)$ (at time $t$), which
has been deformed by the nonlinear $\phi^{t}$, we apply the usual Laplace
operator $\Delta$ to the function $f\circ\phi^{-t}$, which is anchored at time
$t$. Having formed $\Delta(f\circ\phi^{-t})$ (which is just another function
anchored at time $t$), we need to pull it back to the initial time 0. We do
this by precomposing with $\phi^{t}$, leading finally to
$(\Delta(f\circ\phi^{-t}))\circ\phi^{t}$ as the contribution to the dynamic
Laplacian from time $t$. We average these contributions across the interval
$[0,t_{f}]$ to obtain

$\Delta^{D}_{[0,t_{f}]}:=\frac{1}{t_{f}}\int_{0}^{t_{f}}(\Delta(f\circ\phi^{-t}))\circ\phi^{t}\
dt=\frac{1}{t_{f}}\int_{0}^{t_{f}}(\phi_{t})^{*}\circ\Delta\circ(\phi_{t})_{*}\
dt,$ (2)

where the second equality uses the standard push-forward and pull-back
notation $(\phi_{t})_{*}f:=f\circ(\phi_{t})^{-1}$ and
$(\phi_{t})^{*}f:=f\circ\phi_{t}$. Just as the eigenfunctions of the standard
Laplace operator $\Delta$ contain information about regions with minimal
boundary length relative to enclosed area, the leading eigenfunctions of the
dynamic Laplace operator $\Delta^{D}_{[0,t_{f}]}$ encode regions with _minimal
average boundary length under evolution with $\phi_{t}$_, namely _finite-time
coherent sets_. Using Argo float trajectories, we will infer $\phi_{t}$,
numerically construct $\Delta^{D}_{[0,t_{f}]}$, and compute its leading
eigenfunctions.

##  3 Approximating dynamic Laplacian eigenfunctions from Argo data

###  3.1 Description of the Argo-based trajectory data

The Argo program is an international collaboration which began in the late
1990s with the goal of providing freely available subsurface temperature and
salinity measurements to understand ocean climate variability. By deploying a
fleet of autonomous devices called floats, the Argo program has been able to
capture measurements in the upper 2000m of the water column over most of the
global ocean. Currently, Argo is the sole global subsurface data set. The
floats are approximately one metre in length and are able to alter their
density by extruding mineral oil contained inside them into an external
bladder. This mechanism allows them to descend to a predetermined pressure,
called the parking pressure, which allows them to drift at a depth of
approximately 1000m below sea level. Figure 2 shows the number of Argo floats
that were active in the 72 months from January 2011 to January 2017, and the
distribution of lifetimes of the Argo floats.

Figure 2: Left: Active floats vs. month. Right: Histogram of float lifetimes
over six years.

Floats in the Argo program follow cycles that involve an extended period at a
depth followed by a brief period at the surface. A typical approximately
10-day cycle begins as the device descends from the surface to its parking
pressure at around 1000m, where it spends on average 9 days. The device then
descends to a maximum depth of 2000m, after which it returns to the surface.
The average displacement of each of the floats during the subsurface component
of each cycle is 37.64km and each float completes on average 140 cycles over
its lifespan. During the surface component of each cycle, lasting an average
of 8 hours, the device transmits its location on average 9 times, with an
accuracy of approximately 100m. We will use the initial location of the float
upon surfacing. The raw netCDF data is freely available from the GDAC websites
(http://www.coriolis.eu.org/ or http://www.usgodae.org/Argo/). Full
specifications can be found in [Arg19].

The Argo floats are deployed individually, and their surfacing times are not
coordinated. We create a set of Lagrangian trajectories, approximately
recorded monthly as follows. There are $I=8140$ unique floats that are active
at _some_ time in the 6-year period under consideration; we index these floats
by $i\in\\{1,2,\ldots,I\\}$. For each month, indexed by
$t\in\\{1,2,\ldots,T\\}$, $T=72$, we record the surfacing locations of floats
that surface between the 1st and 12th days of that month. We then create
trajectories of monthly float locations $\\{x_{i,t}\\}_{1\leq i\leq I,1\leq
t\leq T}$ by

$x_{i,t}=\left\\{\begin{array}[]{ll}\mbox{first surfacing location},&\hbox{if
float $i$ surfaces between the 1st and 12th of month $t$;}\\\
\emptyset,&\hbox{otherwise.}\end{array}\right.$ (3)

For later use we define the set of reporting floats $R_{t}$ at time $t$ by
$R_{t}=\\{i:x_{i,t}\neq\emptyset\\}$, $t=1,\ldots,T$. We will additionally use
a discrete set of $C=2712$ points $\\{y_{j}\\}_{1\leq j\leq C}$ sampled on the
coastlines of all continents and large islands. These $C$ points were obtained
from the Matlab Mapping Toolbox’s `worldmap` dataset, subsampling every 5th
point and manually ensuring large islands were well represented. We emphasise
that we use each space/time-stamp of each float only once in our computations,
in contrast to e.g. [VSEF12, MBVO+17, MBVO+19] where a float position at one
time is reused many times over.

###  3.2 Numerical approximation of dynamic Laplacian and its eigenproblem

Techniques of numerically approximating the dynamic Laplacian eigenproblem
include [Fro15, FJ15, BK17]; related constructions include [HKTH16] and some
aspects of [PGS17]. To handle the Argo data we use the “adaptive TO” (adaptive
transfer operator) approach from [FJ18] based on the finite-element method. In
addition to having a very well-developed theoretical [Fro15] and numerical
[FJ18] foundation, this approach has several practical advantages:

  * •

the trajectories can be scattered in space and sparse (the Argo floats are
sparsely distributed),

  * •

the domain can have irregular boundaries (the ocean is bounded by irregular
continents),

  * •

complete trajectories are not required (Argo floats only report over some
fraction of the timespan $[0,t_{f}]$),

  * •

no spatial derivatives are required (this would be impossible with sparse Argo
data),

  * •

eigenfunctions are reconstructed on the full phase space, not only at data
points,

  * •

there are no parameters to select.

We consider the eigenproblem of the dynamic Laplacian with Dirichlet boundary
conditions, namely $\Delta^{D}_{[0,t_{f}]}f=\lambda f$ on the interior of $X$
and $f\equiv 0$ on $\partial X$. The dynamic Laplacian is a symmetric,
elliptic operator [Fro15, FJ18] and its eigenspectrum is real and nonpositive.
We are interested in the eigenfunctions $f$ corresponding to the large values
of $\lambda$ (i.e. those negative eigenvalues closest to 0).

Using the central term of (2), this eigenproblem may be written as

$\frac{1}{t_{f}}\int_{0}^{t_{f}}(\Delta(f\circ\phi_{-t}))\circ\phi_{t}\
dt=\lambda f.$

Multiplying through by a function $g:X\to\mathbb{R}$, integrating over $X$,
and applying integration by parts to move one spatial derivative from $f$ to
$g$, we have the weak form of this eigenproblem:

$-\frac{1}{t_{f}}\int_{0}^{t_{f}}\int_{X}\nabla(f\circ\phi_{-t})(x)\cdot\nabla(g\circ\phi_{-t})(x)\
dx\ dt=\lambda\int_{X}f(x)g(x)\ dx.$ (4)

Our numerical approximation will be based upon the float locations
$\\{x_{i,t}\\}$ in (3) for a discrete set of $t$ in the interval $[0,t_{f}]$.
We follow Sections 3, 3.2, 3.2.3, and 3.2.4 in [FJ18]. The main preparatory
steps are:

  1. 1.

We discretise the interval $[t_{0},t_{f}]$ into times indexed by
$t=1,\ldots,T$. In the computations presented here, the time index $t$
indicates an integer number of months, and $T=72$.

  2. 2.

At each $t=1,\ldots,T$ we mesh the positions of the reporting floats and
coastline points $\\{x_{i,t}:i\in R_{t}\\}\cup\\{y_{j}:1\leq j\leq C\\}$; see
Figure 3. This mesh covers a region that is very close to $X$ and we will
henceforth not distinguish between the domain covered by the mesh and $X$.

  3. 3.

At each $t=1,\ldots,T$ using the mesh at time $t$, we build a collection (a
linear basis) of piecewise linear hat functions
$\\{\varphi_{1}^{t},\ldots,\varphi_{n_{t}}^{t}\\}$, where $n_{t}=|R_{t}|+C$.
We will approximate eigenfunctions of the dynamic Laplacian on $X$ by linear
combinations of the $\\{\varphi_{i}^{1}\\}_{i=1}^{n_{t}}$. A sketch of one of
these hat functions is illustrated in Figure 4.

Figure 3: Mesh at $t=36$ months. Vertices of the mesh are either the positions
of Argo floats at 36 months or points on coastlines. Figure 4: Example graph
of a single piecewise linear hat function $\varphi^{t}_{i}$ on a two-
dimensional mesh. The function $\varphi_{i}^{t}$ takes the value 1 at the
single vertex $x_{i,t}$ and is zero at neighboring vertices (and all other
vertices). Reproduced from [AFW10].

By linearity of (4) in $f$ and $g$ we may consider basis pairs
$f=\varphi^{t}_{i}$ and $g=\varphi_{j}^{t}$, and for each $t=1,\ldots,T$
compute the $n_{t}\times n_{t}$ “stiffness” matrices

$D^{t}_{ij}=\int_{X}\nabla\varphi_{i}^{t}(x)\cdot\nabla\varphi_{j}^{t}(x)\
dx,$ (5)

for $i,j\in R_{t}$, setting $D^{t}_{ij}=0$ otherwise. Similarly, for each
$t=1,\ldots,T$ we compute “mass” matrices

$M^{t}_{ij}=\int_{X}\varphi_{i}^{t}(x)\varphi_{j}^{t}(x)\ dx,$ (6)

for $i,j\in R_{t}$, setting $M^{t}_{ij}=0$ otherwise. We note that because the
$\varphi_{i}^{t}$ are piecewise linear functions of space, the integrands in
(5) and (6) are piecewise constant and quadratic, respectively and can be
efficiently evaluated numerically [FJ18].

We now wish to sum the matrices $D^{t}$ and $M^{t}$ over $t$ and so for each
$t$ we inflate each $D^{t}$ and $M^{t}$ from size
$(|R^{t}|+C)\times(|R^{t}|+C)$ to size $(I+C)\times(I+C)$ by adding zero rows
and columns at indices in $I\setminus R_{t}$. For $u\in\mathbb{R}^{I+C}$, the
discrete form of (4) is the eigenproblem

$-\left(\frac{1}{T}\sum_{t=1}^{T}D^{t}\right)u=\lambda\left(\frac{1}{T}\sum_{t=1}^{T}M^{t}\right)u,\qquad
u_{i}=0\mbox{ for }i=I+1,\ldots,I+C,$ (7)

where we fix $u_{i}$ to 0 for all indices corresponding to coastline points in
order to enforce Dirichlet boundary conditions (namely the last $C$ points
according to our indexing). In this study we have selected Dirichlet boundary
conditions because we are seeking slowly mixing structures that do not
intersect with coastlines.

###  3.3 Eigenvector and finite-time coherent set computations

We use (7) to find eigenvalue, eigenvector pairs $(\lambda,u)$. From the
coefficient vector $u$, we approximate the corresponding eigenfunction $f$ by
$f(x)=\sum_{i=1}^{I+C}u_{i}\varphi_{i}^{1}(x)$, which is anchored at time 0.
We compute the leading eight eigenvectors, denoted $u^{(k)}$, $k=1,\ldots,8$,
and from these we construct the leading eight eigenfunctions
$f^{(k)}=\sum_{i=1}^{I+C}u^{(k)}_{i}\varphi_{i}^{1}$; see Figure 5.

Figure 5: Leading eight eigenfunctions $f^{(k)}$, $k=1,\ldots,8$, of the
dynamic Laplacian computed from Argo trajectory data over six years. These
eigenfunctions are represented as scalar fields at January 2011 (month # 1 for
our trajectories). Deep red and deep blue regions correspond to finite-time
coherent sets.

We may estimate the forward evolved eigenfunction $f\circ\phi_{-t}$ at month
$t$ by $(f\circ\phi_{-t})(x)\approx\sum_{i=1}^{I+C}u_{i}\varphi_{i}^{t}(x)$.

The finite-time coherent sets are encoded in Figure 5 as regions in the ocean
that have extreme positive or negative values in one or more of the
eigenfunctions $f^{(k)}$, $k=1,\ldots,8$. Typically the spatial scale of the
features highlighted in the eigenfunctions decreases as one proceeds down the
eigenspectrum. Some of these extreme value regions occur in more than one
eigenfunction and so to produce a clearer set of images, where there is
_exactly one highlighted feature per function_ , we apply the SEBA algorithm.
SEBA (Sparse EigenBasis Approximation) [FRS19]333Matlab code is listed in
[FRS19], and Matlab and julia code may be downloaded at
http://www.maths.unsw.edu.au/~froyland/software.html. produces a sparse basis
that approximately spans a given eigenbasis; in the present case the span of
$f^{(1)},\ldots,f^{(8)}$. We input the eigenvectors $u^{(1)},\ldots,u^{(8)}$
to SEBA, and the algorithm outputs sparse vectors $s^{(1)},\ldots,s^{(8)}$ of
the same size. The sparse vectors $s^{(k)}$ mostly take values between 0 and 1
and the number $s^{(k)}_{i}$ can be interpreted as the likelihood that the
trajectory $x_{i,t},t=1,\ldots,T$ belongs to a coherent set. From these sparse
coefficient vectors, we may construct sparse functions
$\hat{f}^{(k)},k=1,\ldots,8$ as
$\hat{f}^{(k)}=\sum_{i=1}^{I+C}s^{(k)}_{i}\varphi_{i}^{1}$, shown in Figure 6.

Figure 6: Result of the application of the SEBA algorithm [FRS19] to the
leading 8 eigenvectors to separate the eight main coherent sets. These sparse
basis functions $\hat{f}^{(k)}$, $k=1,\ldots,8$, are represented as scalar
fields at January 2011 (month # 1 for our trajectories). Highlighted regions
are the locations of the finite-time coherent sets at January 2011.

##  4 Discussion of the results of the finite-time coherent set analysis

In Figure 7 we combine the eight coherent features displayed in Figure 6 into
a single image by maximising across the sparse vectors $s^{(k)}$.

Figure 7: Maximum function $f^{\rm max}$ formed from maximising the sparse
basis functions $\hat{f}^{(1)},\ldots,\hat{f}^{(8)}$ shown in Figure 6,
represented as a scalar field at January 2011 (month # 1 for our
trajectories). Highlighted regions are the locations of the finite-time
coherent sets at January 2011.

Set $s^{\rm max}:=\max\\{s^{(1)},\ldots,s^{(8)}\\}$, where entry-wise
maximisation is meant; i.e. $s^{\rm
max}_{i}=\max\\{s^{(1)}_{i},\ldots,s^{(8)}_{i}\\}$ for each $i$. We generate
the corresponding maximum function as $\hat{f}^{\rm
max}=\sum_{i=1}^{I+C}s^{\rm max}_{i}\varphi_{i}^{1}$.

In Figure 7 we see as colour highlights the eight most coherent deep ocean
objects over a six-year time frame; these correspond to several known
oceanographic features. In the Atlantic, the 6-year coherent sets highlight
the North and South Atlantic Subtropical Gyres in red, while the equatorial
countercurrent and North Atlantic Current are excluded in dark blue. In the
North Pacific, the Subpolar Gyre is separated from the Subtropical Gyre, while
in the South Pacific, the South Pacific Gyre is separated into eastern and
western components. In the Indian Ocean, the northern extent of the Indian
Ocean Gyre is limited by stronger mixing due to the equatorial countercurrent.
In each of the Atlantic, Pacific, and Indian Oceans, the western extents of
the highlighted features is limited by the stronger mixing of the western
boundary currents. The Weddell Gyre off the Antarctic coast is visible for the
first 20 months of the 72-month flow, but after that time all floats in our
dataset disappear from the Weddell Gyre and therefore its signature can no
longer be detected. In Figure 8 we show average velocity of the floats across
the six-year period, computed directly from monthly spatial increments of the
floats, and interpolated using the basis functions $\varphi^{t}_{i}$.

Figure 8: Plot of RMS speed in km/month, computed from Argo floats monthly
across all 72 months, illustrated with basis functions $\varphi_{i}^{36}$ at
month 36.

There is a partial negative correlation between average float speed according
to Figure 8 and the likelihood of float membership in a coherent set shown in
Figure 7, but fast currents need not imply the absence of coherent dynamics,
nor need the converse implication be true.

###  4.1 Estimating optimal boundaries of the six-year coherent sets

The likelihood representation of Figure 7 is useful to gauge the relative
coherence of regions, with red indicating greater coherence. We can determine
a boundary of each of these coherent sets highlighted in Figure 7 by using
(1). Recall that each $\hat{f}^{(k)}$ takes positive values on a coherent set.
We will use superlevel sets $A^{k}_{c}:=\\{x\in X:\hat{f}^{(k)}\geq c\\}$ as
candidates for $A$ in (1). Recall that we are interested in sets $A$ that
minimise $h(A)$ because such regions have minimal average evolved boundary
length relative to enclosed area, and therefore such a set will have minimal
mixing in the presence of small diffusion. By using level sets we reduce the
minimisation problem to considering a one-parameter family of sets $A_{c}$
indexed by the scalar $0<c<1$ and minimise over $c$. In practice, for a fixed
$k$ we (i) discretise the range of $c$ into steps of 1/100, (ii) evaluate the
right-hand-side of (1) for each discrete value of $c$, and (iii) select a
local or global minimiser $c$, denoted $c_{\min}$. In step (ii), for a given
$c$ we must estimate the boundary length and area of the evolved $A^{k}_{c}$
at each discrete time $t=1,\ldots,T$, namely
$\ell(\partial(\phi_{t}(A^{k}_{c})))$, and $a(\phi_{t}(A^{k}_{c}))$.

Estimating the length of the boundary of $\phi_{t}(A^{k}_{c})$: To estimate
$\ell(\partial(\phi_{t}(A^{k}_{c})))$ for each $t=1,\ldots,T$ we interpolate
the values $s^{(k)}_{i}$ at the locations $x_{i,t}$ onto a fine, uniform
spatial longitude/latitude grid (1 degree by 1 degree) and then use a
contouring routine. For example, with Matlab, we used `contourcs` to obtain
`contourc` output as a struct
array444https://www.mathworks.com/matlabcentral/fileexchange/28447-contourcs-
to-obtain-contourc-output-as-a-struct-array. This contouring routine outputs a
sequence of points so that the polygonal arcs connecting those points
approximates the level set $\\{x\in
X:\sum_{i=1}^{I+C}s_{i}^{(k)}\varphi_{i}^{t}=c\\}$ (recall the
$\varphi_{i}^{t}$ are the piecewise affine hat function basis constructed from
the mesh at time $t$), which is the boundary of $A^{k}_{c}(t)$. The length of
the boundary of the resulting polygon (or union of polygons) is then easily
computed by the Euclidean lengths of its edges, scaled appropriately according
to latitude.

Estimating the area of $\phi_{t}(A^{k}_{c})$: To estimate the area of
$A^{k}_{c}(t)$ we first apportion an area to each grid point in our 1 degree
by 1 degree uniform grid, taking latitude into account. We then simply sum
those areas corresponding to grid points whose interpolated $s_{i}^{(k)}$
value is above $c$.

Putting the above boundary length and area calculations together, we have for
a given threshold $c$, estimates of $\ell(\partial(\phi_{t}(A^{k}_{c})))$, and
$a(\phi_{t}(A^{k}_{c}))$ for each $t=1,\ldots,72$. These values are input to
(1), where the integral in (1) becomes a discrete sum across $t=1,\ldots,72$.
The result is a value $h(A^{k}_{c})$. We create these numbers for a discrete
collection of $c$ between 0 and 1 and select a local or global minimiser
$c_{\min}$. Once such a $c_{\min}$ has been found, this defines our coherent
set $A_{c_{\min}}$ at our initial time. Having fixed $c=c_{\min}$, to
construct the evolution of this coherent set over the 72-month time frame, we
create a sequence of sets using the same contouring approach for each
$t=1,\ldots,T$, using the _same_ threshold value $c_{\min}$. These 72
boundaries for the coherent set associated with the eastern component of the
South Pacific Subtropical Gyre, the North and South Atlantic, and the Indian
oceans are plotted in red in Figure 9 (upper).

Figure 9: Boundaries of finite-time coherent sets (SEBA function contours) in
the eastern South Pacific, North and South Atlantic, and Indian Oceans (red)
together with climatalogical absolute value of planetary potential vorticity
$fN^{2}$ calculated from the World Ocean Atlas [LMB+18, ZSB+19]. The upper
panel shows the potential vorticity at 1000 db (the Argo float parking depth),
with contours of potential density $\sigma_{1}$. The middle panel shows the
same field on the $\sigma_{1}=31.9$ potential density surface, with contours
of depth. The lower panel shows a cross section through the middle of the SEBA
coherent region in the eastern South Pacific, along the section indicated by
the dashed line. Black contours are $\sigma_{1}$ isopycnals.

We emphasise that the concept of coherent sets is specifically designed to
allow _spatial motion_ of the sets over time. This is distinct from an earlier
concept of almost-invariant sets [DJ99], where approximately _fixed-in-space_
objects, such as the Antarctic gyres, are mapped [FPET07, DFH+09]. In the
present experiment, we have relatively sparse Argo float information and our
dynamic Laplacian approach allows us to reliably identify large basin-scale
and sub-basin-scale coherent features. While the concept of coherent sets
allows these features to move in time, because of their large spatial scale in
this application, they naturally tend to remain approximately in place, as
illustrated by the red curves in Figure 9 (upper).

##  5 Comparison with large-scale ocean tracers

In this section, we look for relationships between the coherent sets detected
by the dynamic Laplacian method and large-scale oceanographic tracers. Our
working hypothesis is that the coherent sets represent regions within which
water parcels can exchange freely but are relatively isolated from exchange
outside of the region. A complicating factor, however, is that the Argo
trajectories are isobaric and thus approximately constrained to the 1000 m
depth plane, while water parcel trajectories are in general three dimensional.
Specifically, in the quasi-adiabatic ocean interior, we expect water parcel
trajectories to lie on the neutral plane.

A central feature of the mid-depth ocean circulation is the existence of
“shadow zones” on the Eastern boundaries. Shadow zones are predicted by
classic ventilated thermocline theory [LPS83, Ped86] as regions for which the
geostrophic streamfunction forms closed interior contours and does not connect
with the surface, thereby isolating the shadow zones from surface ventilation.
One prediction of this theory is that, within such isolated regions, potential
vorticity becomes homogenized due to eddy mixing [RY82]. If the coherent set
detected by Argo floats corresponds with the shadow zone, we might then expect
potential vorticity to be homogenized therein.

To test this, we examine planetary potential vorticity (PV), defined as
$fN^{2}$, where $f$ is the Coriolis parameter and $N$ is the Brunt-Väisälä
frequency (proportional to stratification). The frequency $N$ was calculated
from climatalogical temperature and salinity fields provided by the World
Ocean Atlas [LMB+18, ZSB+19]. In Figure 9 (top), we plot the absolute value of
PV at 1000 m depth, together with the boundaries of the coherent sets in four
major ocean basins (Indian, eastern South Pacific, North Atlantic, and South
Atlantic). The South Atlantic set is anomalously small, and therefore does
represent a region of relatively homogeneous PV. In the three largest sets
(Indian, eastern South Pacific, and North Atlantic) rather than being
homogenized within the coherent set, PV shows a strong meridional gradient
across the region. To test whether these gradients are an artifact of
projection onto a depth surface, we also examine PV on an isopycnal surface
that intersects the 1000 m surface in mid-latitudes ($\sigma_{1}=31.9$ kg
m-3); the same large scale gradients are present within the coherent sets.
This indicates that the dynamically isolated regions identified by the dynamic
Laplacian analysis do not correspond closely with the classical shadow zones.
Moreover, textbook shadow zones—with homogenized PV along the eastern boundary
of the subtropical gyre—are not evident in the PV field at these depths.

Figure 10: As for Fig. 9 except color shading now shows World Ocean Atlas
Oxygen saturation [GWP+19] in blue.

Another prominent tracer feature at these depths are Oxygen Minimum Zones
(OMZs). OMZs are mid-depth layers of low O2 concentration found in most ocean
basins [PRP09] and references therein. OMZs are generally anchored to the
Eastern boundary. Their vertical extent begins near the surface but can reach
as deep as 4000m in some ocean basins. Most OMZs thus intersect with the 1000m
depth surface sampled by Argo trajectories. OMZs are formed through the
combined effect of biology, circulation, and mixing; biological consumption is
ultimately responsible for depleting oxygen, but the effects of circulation
help determine the shape and position of the OMZs [Wyr62]. In particular,
regions of weak deep circulation and low mixing are associated with OMZs. It
is reasonable to ask, therefore, whether the dynamically isolated regions
revealed by the dynamic Laplacian analysis correspond with OMZs.

To test this hypothesis, in Figure 10 we examine the boundaries of the
coherent sets together with the oxygen saturation level from the gridded World
Ocean Atlas dataset [GWP+19]. Similarly to PV, there is no obvious overlap
between low and / or homogenous oxygen saturation levels and the coherent
sets. Instead, the South Pacific and South Indian sets cross directly over the
strong subtropical oxygen gradient, while the North Atlantic set occupies the
entire basin. The coherent sets do not correspond to the zonally elongated,
Eastern intensified regions characteristic of OMZs [PRP09].

The overall lack of correspondence between large-scale ocean tracers (PV and
oxygen saturation) and Argo coherent sets indicates that processes besides
lateral transport at 1000 dbar must be playing a significant role in
determining these tracer distributions. Most notably the Argo floats
experience no vertical transport, their ballasting confining them to the 1000
dbar surface. Real water parcels, in contrast, move predominantly along
neutral surfaces (isopycnals) at these depths. While the local difference
between the neutral angle and the 1000 dbar surface may be small at any
particular point, over long trajectories, parcels originating at 1000 dbar may
drift significantly away from this depth as they follow isopycnals over the
basin scales representative of the coherent sets. The potential mismatch
between isobaric and isopycnal trajectories is therefore a major caveat around
the coherent sets we have inferred from Argo trajectories. Furthermore,
diapycnal advection and mixing likely also play a strong role in tracer
budgets at 1000 db, causing the observed tracer distributions to diverge from
the patterns revealed by the Argo-based coherent sets.

##  6 Conclusions

The dynamic Laplacian approach [Fro15] to identifying coherent sets and its
FEM implementation [FJ18] has previously been deployed on eddy detection in
the Agulhas current [FJ18] and in the North Atlantic [FRS19], where in both
cases an altimetry-derived velocity field was used. The eddy tracking
experiment in [FJ18] used trajectory sampling that was rich enough to compute
reliable spatial derivatives. This was relaxed in [FRS19] below the level
where reliable spatial derivatives of the flow map could be calculated,
without impacting on the performance or resolution of the eddy tracking. In
the present paper we have dramatically reduced the input trajectory resolution
to demonstrate that the dynamic Laplacian approach is well suited to the
extremely spatially sparse trajectory data that is typically recorded from
drifter datasets such as Argo. Moreover, the approach easily accommodates the
irregular temporal sampling of this dataset.

Using Argo trajectories over six years we identified the eight most coherent
subsurface ocean features on this timescale (Figure 6). Because of the global
nature of the trajectory data, many of these dominant coherent regions are
aligned with the major ocean boundaries, but with a bias toward the east of
the basin due to the stronger mixing western boundary currents that destroy
coherent water motion. In particular, the vicinities of the Mozambique,
Brazil, Kuroshio, and East Australian currents are conspicuously not
highlighted in Figure 6. There is a separation between the North and South
Atlantic associated with the equatorial countercurrent and for the same reason
the main coherent region in the Indian Ocean is pushed southwards. In the
North Pacific, the separation between the subpolar and subtropical gyres is
evident, and the South Pacific gyre splits into eastern and western
components. To our knowledge, the dynamical separation between these two parts
of the subtropical gyre has not previously been identified and merits further
investigation. In contrast to ocean eddies, which move considerably in space
over time, the global-scale coherent regions we have identified move only
slightly throughout the six-year duration, reflecting the relatively
stationary nature of the ocean dynamics at these large spatial and temporal
scales.

We hope that this analysis can help inform future deployment of Argo floats.
The coherent sets we have identified represent geographical regions from which
floats are unlikely to escape. Floats deployed near the center of the coherent
sets are nearly certain to remain confined to the regions shown in Figure 6.
In contrast, floats deployed outside the sets, or near set boundaries, may be
free to wander into different regions. This information could help target
float deployment to sample specific regions.

Another potential application of the method introduced in this paper would be
to study the dynamic geography using other Lagrangian float ensembles. The
most obvious candidate would be the surface global drifters. It would also be
interesting to explore the regional dynamic geography of ocean basins with a
high density of Lagrangian floats. For example, the Consortium for Advanced
Research on the Transport of Hydrocarbon in the Environment deployed over 300
drifters in the Gulf of Mexico in 2012 [MRH+16], the so-called Grand
Lagrangian Deployment, GLAD. The high density of drifters in the Nordic Seas
[KLA+11] makes this another attractive region to explore.

We hypothesized that the Argo coherent sets could be somehow aligned with
large-scale ocean tracers and tested this by visualizing the relationship
between set boundaries and isosurfaces of potential vorticity and oxygen
saturation. In general, both tracers and coherent sets show zonal asymmetry,
with more coherence / tracer homogenization towards the Eastern boundary.
However, there was no clear smoking-gun relationship found in the tracers we
examined. Instead, the coherent sets mostly tended to occupy large regions, up
to the entire basin scale, while the tracers contain notable features (e.g.
shadow zones, oxygen minimum zones) on smaller scales. The coherent sets
tended to overlie the strong subtropical gradient present in both tracers,
suggesting that the sets do not represent regions of clear tracer
homogenization. The mismatch between the tracer isosurfaces and the coherent
sets is likely due to three-dimensional effects (e.g. flow along isopycnals
rather than isobars, overturning circulation) and non-conservative processes
(e.g. mixing, oxygen consumption). Future work examining coherent sets in
three dimensions could help resolve these questions. Such work would have to
rely on numerical trajectories from an ocean model, as all existing real
floats and drifters are constrained to 2D motion.

##  7 Acknowledgements

RPA acknowledges support from NASA award NNX 80NSSC19K1252. GF thanks the
Banff International Research Station for supporting a five-day workshop in
January 2017 where initial discussions for this work took place. The research
of GF has been partially supported by two ARC Discovery Projects over the
course of this research. KS was supported by an ARC Discovery Project.

## References

  * [AFW10] Douglas Arnold, Richard Falk, and Ragnar Winther.  Finite element exterior calculus: from Hodge theory to numerical stability.  Bulletin of the American Mathematical Society, 47(2):281–354, 2010\.
  * [Arg19] Data Management Team Argo.  Argo user’s manual v3.3.  Report, 2019.
  * [BK17] Ralf Banisch and Péter Koltai.  Understanding the geometry of transport: Diffusion maps for Lagrangian trajectory data unravel coherent sets.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 27(3):035804, 2017.
  * [BLSF21] Dhruv Balwada, Joseph H. LaCasce, Kevin G. Speer, and Raffaele Ferrari.  Relative dispersion in the Antarctic Circumpolar Current.  Journal of Physical Oceanography, 51(2):553–574, Feb 2021.
  * [BSL+16] Dhruv Balwada, Kevin G. Speer, Joseph H. LaCasce, W. Brechner Owens, John Marshall, and Raffaele Ferrari.  Circulation and Stirring in the Southeast Pacific Ocean and the Scotia Sea Sectors of the Antarctic Circumpolar Current.  Journal of Physical Oceanography, 46(7):2005–2027, Jul 2016.
  * [Cha84] Isaac Chavel.  Eigenvalues in Riemannian geometry.  Academic press, 1984.
  * [Chi13] Stephen M. Chiswell.  Lagrangian time scales and eddy diffusivity at 1000 m compared to the surface in the South Pacific and Indian Oceans.  Journal of Physical Oceanography, 43(12):2718–2732, Dec 2013.
  * [CTF+17] Lijing Cheng, Kevin E. Trenberth, John Fasullo, Tim Boyer, John Abraham, and Jiang Zhu.  Improved estimates of ocean heat content from 1960 to 2015.  Science Advances, 3(3):e1601545, Mar 2017.
  * [DFH+09] Michael Dellnitz, Gary Froyland, Christian Horenkamp, Kathrin Padberg-Gehle, and A Sen Gupta.  Seasonal variability of the subpolar gyres in the southern ocean: a numerical investigation based on transfer operators.  Nonlinear Processes in Geophysics, 16(6), 2009.
  * [Dij05] Henk A Dijkstra.  Nonlinear physical oceanography: a dynamical systems approach to the large scale ocean circulation and El Nino, volume 28.  Springer Science & Business Media, 2005.
  * [DJ99] Michael Dellnitz and Oliver Junge.  On the approximation of complicated dynamical behavior.  SIAM Journal on Numerical Analysis, 36(2):491–515, 1999.
  * [FCH+15] G. Forget, J.-M. Campin, P. Heimbach, C. N. Hill, R. M. Ponte, and C. Wunsch.  Ecco version 4: an integrated framework for non-linear inverse modeling and global ocean state estimation.  Geoscientific Model Development, 8(10):3071–3104, Oct 2015.
  * [FHR+12] Gary Froyland, Christian Horenkamp, Vincent Rossi, Naratip Santitissadeekorn, and Alex Sen Gupta.  Three-dimensional characterization and tracking of an Agulhas Ring.  Ocean Modelling, 52:69–75, 2012.
  * [FHRVS15] Gary Froyland, Christian Horenkamp, Vincent Rossi, and Erik Van Sebille.  Studying an Agulhas ring’s long-term pathway and decay with finite-time coherent sets.  Chaos, 25(8):083119, 2015.
  * [FJ15] Gary Froyland and Oliver Junge.  On fast computation of finite-time coherent sets using radial basis functions.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 25(8):087409, 2015.
  * [FJ18] Gary Froyland and Oliver Junge.  Robust FEM-based extraction of finite-time coherent sets using scattered, sparse, and incomplete trajectories.  SIAM Journal on Applied Dynamical Systems, 17(2):1891–1924, 2018\.
  * [FPET07] Gary Froyland, Kathrin Padberg, Matthew H England, and Anne Marie Treguier.  Detection of coherent oceanic structures via transfer operators.  Physical Review Letters, 98(22):224503, 2007.
  * [FPG15] Gary Froyland and Kathrin Padberg-Gehle.  A rough-and-ready cluster-based approach for extracting finite-time coherent sets from sparse and incomplete trajectory data.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 25(8):087406, 2015.
  * [Fro13] Gary Froyland.  An analytic framework for identifying finite-time coherent sets in time-dependent dynamical systems.  Physica D: Nonlinear Phenomena, 250:1–19, 2013.
  * [Fro15] Gary Froyland.  Dynamic isoperimetry and the geometry of Lagrangian coherent structures.  Nonlinearity, 28(10):3587, 2015.
  * [FRS19] Gary Froyland, Christopher P Rock, and Konstantinos Sakellariou.  Sparse eigenbasis approximation: Multiple feature extraction across spatiotemporal scales with application to coherent set identification.  Communications in Nonlinear Science and Numerical Simulation, 77:81–107, 2019.
  * [FSM10] Gary Froyland, Naratip Santitissadeekorn, and Adam Monahan.  Transport in time-dependent dynamical systems: Finite-time coherent sets.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 20(4):043116, 2010.
  * [FSvS14] Gary Froyland, Robyn M Stuart, and Erik van Sebille.  How well-connected is the surface of the global ocean?  Chaos: An Interdisciplinary Journal of Nonlinear Science, 24(3):033126, 2014.
  * [GN13] Denis S Grebenkov and B-T Nguyen.  Geometrical structure of Laplacian eigenfunctions.  SIAM Review, 55(4):601–667, 2013.
  * [GR14] Alison R. Gray and Stephen C. Riser.  A global analysis of Sverdrup balance using absolute geostrophic velocities from Argo.  Journal of Physical Oceanography, 44(4):1213–1229, Apr 2014.
  * [GWP+19] HE Garcia, KW Weathers, CR Paver, I Smolyar, TP Boyer, Mm Locarnini, MM Zweng, AV Mishonov, OK Baranova, D Seidov, et al.  World ocean atlas 2018, volume 3: Dissolved oxygen, apparent oxygen utilization, and dissolved oxygen saturation.  2019\.
  * [Hal11] George Haller.  A variational theory of hyperbolic Lagrangian coherent structures.  Physica D: Nonlinear Phenomena, 240(7):574–598, 2011.
  * [HBV13] George Haller and Francisco J Beron-Vera.  Coherent Lagrangian vortices: The black holes of turbulence.  Journal of Fluid Mechanics, 731, 2013.
  * [HKTH16] Alireza Hadjighasem, Daniel Karrasch, Hiroshi Teramoto, and George Haller.  Spectral-clustering approach to Lagrangian vortex detection.  Physical Review E, 93(6):063107, 2016.
  * [HY00] George Haller and Guocheng Yuan.  Lagrangian coherent structures and mixing in two-dimensional turbulence.  Physica D: Nonlinear Phenomena, 147(3-4):352–370, 2000.
  * [JRZ+17] Steven R Jayne, Dean Roemmich, Nathalie Zilberman, Stephen C Riser, Kenneth S Johnson, Gregory C Johnson, and Stephen R Piotrowicz.  The Argo program: present and future.  Oceanography, 30(2):18–28, 2017.
  * [Kac66] Mark Kac.  Can one hear the shape of a drum?  The American Mathematical Monthly, 73(4P2):1–23, 1966.
  * [KLA+11] I. Koszalka, J.H. LaCasce, M. Andersson, K.A. Orvik, and C. Mauritzen.  Surface circulation in the Nordic seas from clustered drifters.  Deep Sea Research Part I: Oceanographic Research Papers, 58(4):468–485, Apr 2011.
  * [KY10] Katsuro Katsumata and Hiroshi Yoshinari.  Uncertainties in global mapping of argo drift data at the parking level.  Journal of Oceanography, 66(4):553–569, Jul 2010.
  * [LMB+18] Mm Locarnini, AV Mishonov, OK Baranova, TP Boyer, MM Zweng, HE Garcia, D Seidov, Kw Weathers, Cr Paver, I Smolyar, et al.  World ocean atlas 2018, volume 1: Temperature.  2018\.
  * [LPS83] J. R. Luyten, J. Pedlosky, and H. Stommel.  The ventilated thermocline.  Journal of Physical Oceanography, 13(2):292–309, Feb 1983.
  * [LYMH07] Konstantin V Lebedev, Hiroshi Yoshinari, Nikolai A Maximenko, and Peter W Hacker.  Velocity data assessed from trajectories of Argo floats at parking level and at the sea surface.  IPRC Technical Note, 4(2), 2007.
  * [MBVO+17] P Miron, Francisco J Beron-Vera, Maria J Olascoaga, J Sheinbaum, P Pérez-Brunius, and G Froyland.  Lagrangian dynamical geography of the Gulf of Mexico.  Scientific reports, 7(1):1–12, 2017.
  * [MBVO+19] P Miron, Francisco J Beron-Vera, Maria J Olascoaga, G Froyland, P Pérez-Brunius, and J Sheinbaum.  Lagrangian geography of the deep Gulf of Mexico.  Journal of Physical Oceanography, 49(1):269–290, 2019.
  * [MRH+16] A. J. Mariano, E. H. Ryan, H. S. Huntley, L.C. Laurindo, E. Coelho, A. Griffa, T. M. Özgökmen, M. Berta, D. Bogucki, S. S. Chen, and et al.  Statistical properties of the surface velocity field in the northern gulf of mexico sampled by glad drifters.  Journal of Geophysical Research: Oceans, 121(7):5193–5216, Jul 2016.
  * [OBVH+13] Maria J Olascoaga, Francisco J Beron-Vera, George Haller, J Trinanes, Mohamed Iskandarani, EF Coelho, Brian K Haus, HS Huntley, G Jacobs, AD Kirwan, et al.  Drifter motion in the Gulf of Mexico constrained by altimetric Lagrangian coherent structures.  Geophysical Research Letters, 40(23):6171–6175, 2013.
  * [Ped86] Joseph Pedlosky.  The buoyancy and wind-driven ventilated thermocline.  Journal of Physical Oceanography, 16(6):1077–1087, Jun 1986.
  * [PGS17] Kathrin Padberg-Gehle and Christiane Schneide.  Network-based study of lagrangian transport and mixing.  Nonlinear Processes in Geophysics, 24(4):661–671, 2017.
  * [Pie94] RT Pierrehumbert.  Tracer microstructure in the large-eddy dominated regime.  Chaos, Solitons & Fractals, 4(6):1091–1110, 1994.
  * [PRP09] A. Paulmier and D. Ruiz-Pino.  Oxygen minimum zones (omzs) in the modern ocean.  Progress in Oceanography, 80(3–4):113–128, Mar 2009.
  * [RBS16] Christopher J. Roach, Dhruv Balwada, and Kevin Speer.  Horizontal mixing in the Southern Ocean from Argo float trajectories.  Journal of Geophysical Research: Oceans, 121(8):5570–5586, Aug 2016.
  * [RCG+15] Dean Roemmich, John Church, John Gilson, Didier Monselesan, Philip Sutton, and Susan Wijffels.  Unabated planetary warming and its ocean structure since 2006.  Nature Climate Change, 5(3):240–245, Feb 2015.
  * [RFR+16] Stephen C. Riser, Howard J. Freeland, Dean Roemmich, Susan Wijffels, Ariel Troisi, Mathieu Belbéoch, Denis Gilbert, Jianping Xu, Sylvie Pouliquen, Ann Thresher, and et al.  Fifteen years of ocean observations with the global Argo array.  Nature Climate Change, 6(2):145–153, Jan 2016.
  * [RT09] Dean Roemmich and Argo Steering Team.  Argo: the challenge of continuing 10 years of progress.  Oceanography, 22(3):46–55, 2009.
  * [RY82] Peter B. Rhines and William R. Young.  Homogenization of potential vorticity in planetary gyres.  Journal of Fluid Mechanics, 122:347–367, 1982.
  * [SHJ+18] Kylie L Scales, Elliott L Hazen, Michael G Jacox, Frederic Castruccio, Sara M Maxwell, Rebecca L Lewison, and Steven J Bograd.  Fisheries bycatch risk to marine megafauna is intensified in Lagrangian coherent structures.  Proceedings of the National Academy of Sciences, 115(28):7362–7367, 2018.
  * [VSEF12] Erik Van Sebille, Matthew H England, and Gary Froyland.  Origin, dynamics and evolution of ocean garbage patches from observed surface drifters.  Environmental Research Letters, 7(4):044040, 2012.
  * [VSGA+18] Erik Van Sebille, Stephen M Griffies, Ryan Abernathey, Thomas P Adams, Pavel Berloff, Arne Biastoch, Bruno Blanke, Eric P Chassignet, Yu Cheng, Colin J Cotter, et al.  Lagrangian ocean analysis: Fundamentals and practices.  Ocean Modelling, 121:49–75, 2018.
  * [Wig05] Stephen Wiggins.  The dynamical systems approach to lagrangian transport in oceanic flows.  Annu. Rev. Fluid Mech., 37:295–328, 2005.
  * [WKC12] Darryn W Waugh, Shane R Keating, and Mei-Lin Chen.  Diagnosing ocean stirring: Comparison of relative dispersion and finite-time Lyapunov exponents.  Journal of physical oceanography, 42(7):1173–1185, 2012.
  * [Wyr62] Klaus Wyrtki.  The oxygen minima in relation to ocean circulation.  Deep Sea Research and Oceanographic Abstracts, 9(1–2):11–23, Jan 1962.
  * [ZSB+19] MM Zweng, D Seidov, Tp Boyer, Mm Locarnini, He Garcia, Av Mishonov, OK Baranova, Kw Weathers, CR Paver, I Smolyar, et al.  World ocean atlas 2018, volume 2: Salinity.  2019\.


# Reflection Removal Using
Recurrent Polarization-to-Polarization Network

###### Abstract

This paper addresses reflection removal, which is the task of separating
reflection components from a captured image and deriving the image with only
transmission components. Considering that the existence of the reflection
changes the polarization state of a scene, some existing methods have
exploited polarized images for reflection removal. While these methods apply
polarized images as the inputs, they predict the reflection and the
transmission directly as non-polarized intensity images. In contrast, we
propose a polarization-to-polarization approach that applies polarized images
as the inputs and predicts “polarized” reflection and transmission images
using two sequential networks to facilitate the separation task by utilizing
the interrelated polarization information between the reflection and the
transmission. We further adopt a recurrent framework, where the predicted
reflection and transmission images are used to iteratively refine each other.
Experimental results on a public dataset demonstrate that our method
outperforms other state-of-the-art methods.

Index Terms—  Reflection Removal, Polarization Imaging, Recurrent Neural
Network

##  1 Introduction

Reflections caused by semi-reflectors such as glass are commonly seen in daily
life. When light passes through semi-reflectors, a camera inevitably captures
the reflection and the transmission components at the same time. Nevertheless,
most computer vision applications such as object detection, segmentation, and
depth estimation assume that each pixel value is derived only from the scene
corresponding to the transmission. Therefore, reflection removal is a crucial
task to improve the robustness of real-world applications.

Fig. 1: Different input-output models for the reflection removal. (a) Both the
input and the output of standard single-image methods are intensity images.
(b) Existing polarization-based methods apply polarized images only to the
input. (c) Our proposed polarization-to-polarization approach predicts the
output reflection and transmission as polarized images as well.

Most existing reflection removal methods are based on a single grayscale or
color image, where both the input and the outputs (reflection and
transmission) are an intensity domain, as illustrated by the intensity-to-
intensity model of Fig. 1(a). While recent deep-learning-based methods have
shown great progress [1, 2, 3, 4, 5, 6], the separation of the reflection and
the transmission is still challenging due to an ill-posed problem that an
infinite number of the transmission and the reflection image combinations is
possible to reproduce the same mixed image. Other approaches attempt to solve
this problem by using multi-view color images [7, 8, 9]. However, these
methods typically necessitate image alignment as a pre-processing step, which
imposes constraints on their practical application.

Meanwhile, as the price of one-shot polarization cameras has decreased, one-
shot acquisition of polarized images has become much easier in recent years
[10, 11]. Considering that the existence of reflection components changes the
polarization state of a scene, some non-learning-based [12, 13, 14, 15] or
learning-based [16, 17, 18] methods solve the reflection removal by using a
set of polarized images with different polarizer orientations (typically, four
orientations of $0^{\circ}$, ${45}^{\circ}$, ${90}^{\circ}$, and
${135}^{\circ}$). While these polarization-based methods apply polarized
images as the inputs, they predict the reflection and the transmission images
directly as non-polarized intensity images, as illustrated by the
polarization-to-intensity model of Fig. 1(b).

Fig. 2: The overall structure of our proposed RP2PN.

In this paper, we propose a polarization-to-polarization approach for deep-
learning-based reflection removal, as illustrated in Fig. 1(c). To effectively
learn the polarimetric relationships among the input image and the separated
reflection and transmission images, our approach takes polarized images as the
inputs and predicts the reflection and the transmission images also as the
polarized images. Then, the final reflection and transmission outputs are
derived as the intensity images by averaging the polarized images.

Regarding the network structure, inspired by [5, 6, 18], we propose a two-
stage sequential approach within our polarization-to-polarization framework,
which uses one recurrent network to predict the reflection and then feeds the
reflection result to another network to predict the transmission, as better
transmission estimation is also beneficial for reflection estimation and vice
versa. We also utilize the difference images between different polarizer
angles as the network inputs, because they exhibit an informative feature for
the reflection removal.

Experimental results on a public dataset [18] demonstrate that our method
outperforms other state-of-the-art intensity-based and polarization-based
methods. Additionally, we highlight the significance of our polarization-to-
polarization framework and the effectiveness of the integrated recurrent unit
from the ablation study.

##  2 Proposed Method

###  2.1 Network Structure

Figure 2 shows the overall structure of our proposed recurrent polarization-
to-polarization network (RP2PN), which sequentially and iteratively predicts
the polarized reflection and the polarized transmission images. For network
training, we use Lei et al. real-world dataset [18], which was obtained using
Lucid PHX050S-P one-shot monochrome polarization camera equipped with Sony
IMX250MZR sensor [10]. This dataset provides the triplets of aligned polarized
images $\\{I_{\phi},R_{\phi},T_{\phi}\\}$, where $\phi\in\\{0^{\circ}$,
${45}^{\circ}$, ${90}^{\circ}$, ${135}^{\circ}\\}$ is the polarizer angle, $I$
represents the input image with mixed reflection and transmission, and $R$ and
$T$ represent the corresponding ground-truth reflection and transmission
images, respectively.

Our RP2PN consists of two sequential networks, namely R-LSTM-Net for the
reflection estimation and T-Net for the transmission. As for the network
inputs, from four input polarized images ($I_{0}$, $I_{45}$, $I_{90}$,
$I_{135}$), the intensity image $I$, the degree-of-polarization image ($DoP$)
are calculated by a standard polarimetric calculation. In addition, we
introduce a polarized difference image $I_{diff}$, which serves as informative
cues for the reflection removal.

For a mixed polarized image $I_{\phi}=T_{\phi}+R_{\phi}$ captured under a
certain polarizer orientation $\phi$, the polarized components of $T_{\phi}$
and $R_{\phi}$, denoted as $T^{p}_{\phi}$ and $R^{p}_{\phi}$, will change with
the variation of $\phi$, while the unpolarized components of $T_{\phi}$ and
$R_{\phi}$ remain invariant. Thus, for two mixed images with the polarizer
orientations $\phi_{1}$ and $\phi_{2}$, their difference is formed only by the
polarized components as

$I_{\phi_{1}}-I_{\phi_{2}}=T_{\phi_{1}}^{p}-T_{\phi_{2}}^{p}+(R_{\phi_{1}}^{p}-R_{\phi_{2}}^{p}).$
(1)

We observed in Lei’s real-world dataset [18] that there is a tendency for the
strength of polarization (i.e., $DoP$) of the transmission image to be weaker
than that of the reflection image. For an example depicted in the first row of
Fig. 3, the average $DoP$ values for the ground-truth transmission and
reflection are approximately 0.1 and 0.5, respectively. This indicates that
the polarized $T$ component is considerably weaker than the polarized $R$
component, resulting in the dominance of the $R$ component in the polarized
difference image $I_{diff}$, as shown in the second row of Fig. 3. Based on
this observation, we employ all possible combinations of the four polarizer
angles to compute the difference images, yielding a total of six polarized
difference images ($I_{0,45}$, $I_{0,90}$, $I_{0,135}$, $I_{45,90}$,
$I_{45,135}$, $I_{90,135}$), where
$I_{\phi_{1},\phi_{2}}=|I_{\phi_{1}}-I_{\phi_{2}}|_{1}$.

Fig. 3: An example of $I_{diff}$ image. $I_{0,45}$ demonstrates closer
features to the features of the ground-truth $R$ than either $I_{0}$ or
$I_{45}$. This is due to the polarized $T$ component being considerably weaker
than the polarized $R$ component. The brightness of $I_{0,45}$ is adjusted
solely for the visualization purpose.

An over-exposure binary mask ($M$) is also derived at each pixel as

$M=\left\\{\begin{aligned} 0,\ \ &\text{if\ \
}max(I_{0},I_{45},I_{90},I_{135})>\tau,\\\ \ 1,\ \
&\text{otherwise}.\end{aligned}\right.$ (2)

where $\tau$ is a threshold and set to 0.98 for the pixel value range of
[0,1]. In addition, the features of $I_{0}$, $I_{45}$, $I_{90}$, $I_{135}$,
and $I$ are respectively extracted from a pre-trained VGG-19 network [19]. All
of the above and the original four polarized images are concatenated and fed
to the networks as the inputs.

As for our recurrent structure, we first build R-LSTM-Net, which consists of a
U-Net architecture [20] using a 10-block convolutional encoder and an 8-block
decoder with a long short-term memory (LSTM) unit added in the bottleneck [5]
to predict four polarized reflection images $\hat{R}_{0}$, $\hat{R}_{45}$,
$\hat{R}_{90}$, and $\hat{R}_{135}$. Then, they are concatenated as a part of
the inputs for T-Net with similar U-Net architecture as R-LSTM-Net except for
the LSTM unit to predict four polarized transmission images $\hat{T}_{0}$,
$\hat{T}_{45}$, $\hat{T}_{90}$, and $\hat{T}_{135}$. At last, the predicted
polarized transmission images are used as the inputs to further refine the
reflection result on the next iteration.

###  2.2 Loss Functions

We apply three loss functions to the last iteration’s result of our RP2PN. The
total loss $L_{total}$ is defined as

$L_{total}=\lambda_{1}L_{pixel}+\lambda_{2}L_{percep}+\lambda_{3}L_{pncc},$
(3)

where $\lambda_{1}$, $\lambda_{2}$ and $\lambda_{3}$ are the weighting
parameters.

$L_{pixel}$ is the pixel-wise $L_{1}$ loss between the predicted
($\hat{R},\hat{T}$) and the ground-truth ($R,T$) images to ensure pixel-level
similarity. Different from the existing polarization-based methods [17, 16,
18], which only consider intensity-domain losses, we evaluate the losses for
four polarized images of the reflection and the transmission as

$L_{pixel}=\sum_{\phi\in A}|R^{M}_{\phi}-\hat{R}^{M}_{\phi}|_{1}+\sum_{\phi\in
A}|T^{M}_{\phi}-\hat{T}^{M}_{\phi}|_{1},$ (4)

where $A=\\{0,45,90,135\\}$. The superscript $M$ represents a masked image,
e.g., $R_{\phi}^{M}=R_{\phi}\circ M$, where $\circ$ is the pixel-wise
production.

$L_{percep}$ is the perceptual loss [21] to help the networks to learn high-
level contextual features. Similar to $L_{pixel}$, we here calculate the
losses in the polarized-domain as

$\displaystyle L_{percep}=\sum_{\phi\in
A}\sum_{j}^{N}\gamma_{j}|w_{V}^{j}(R_{\phi}^{M})-w_{V}^{j}(\hat{R}_{\phi}^{M})|_{1}$
(5) $\displaystyle\ +\sum_{\phi\in
A}\sum_{j}^{N}\gamma_{j}|w_{V}^{j}(T_{\phi}^{M})-w_{V}^{j}(\hat{T}_{\phi}^{M})|_{1},$

where $w^{j}_{V}$ expresses the $j$-th layer’s feature map from the pre-
trained VGG-19 network and $\gamma_{j}$ is the weighting parameter of the
$j$-th layer.

$L_{pncc}$ is the perceptual normalized cross-correlation loss [18], which is
applied to minimize the correlation between the predicted reflection and
transmission images, assuming their independency. This loss is applied to the
final intensity output domain as

$L_{pncc}=\sum_{j}^{N}f_{ncc}(w_{V}^{j}(\hat{R}^{M}),w_{V}^{j}(\hat{T}^{M})),$
(6)

where $\hat{R}$ and $\hat{T}$ are the intensity images, which are calculated
by the average of four polarized images, and $f_{ncc}$ is the operator to
calculate the normalized cross-correlation.

##  3 Experimental Results

###  3.1 Implementation Details of Our RP2PN

We used Lei et al. dataset [18], which contains 600, 184, and 107 real-scene
polarized image triplets $\\{I_{\phi},R_{\phi},T_{\phi}\\}$ for training,
validation, and testing, respectively. The weighting parameters in Eq. (3)
were experimentally set as
$\\{\lambda_{1},\lambda_{2},\lambda_{3}\\}=\\{0.1,0.1,6.0\\}$. For the VGG-19
features in Eqs. (5) and (6), we adopted the same six layers ($N=6$) and
weights for each layer as [18]. The number of iterations for RP2PN was
experimentally set to three. To train RP2PN, the learning rate was set to
$1e^{-4}$ at the first 300 epochs with batch size 1. Then, it was reduced to
$1e^{-5}$ for additional 50 epochs. The training took 40 hours using one
Nvidia Geforce RTX 3080 Ti GPU.

Table 1: Quantitative comparisons on Lei et al. dataset [18]. * Non-learning-
based methods (Implementation from [16]). † Learning-based methods using pre-
trained models.

Methods |  With Polar |  Train Data | Transmission | Reflection
---|---|---|---|---
PSNR | SSIM | PSNR | SSIM
Farid* [12] | Yes | - | 25.56 | 0.828 | 24.79 | 0.742
Schechner* [13] | Yes | - | 24.62 | 0.827 | 23.94 | 0.621
BDN†[3] | No | [3] | 24.09 | 0.756 | 23.62 | 0.692
Dong†[6] | No | [6] | 28.30 | 0.864 | 28.79 | 0.659
ReflectNet†[16] | Yes | [16] | 24.76 | 0.821 | 25.03 | 0.715
Lyu†[17] | Yes | [17] | 24.82 | 0.820 | 25.06 | 0.737
Zhang [2] | No | [18] | 32.15 | 0.919 | 32.20 | 0.883
IBCLN [5] | No | [18] | 32.84 | 0.928 | 32.80 | 0.897
Lei [18] | Yes | [18] | 35.00 | 0.950 | 34.58 | 0.921
RP2PN (Ours) | Yes | [18] | 35.87 | 0.954 | 35.63 | 0.933

###  3.2 Comparison with Other Methods

Table 1 summarizes the quantitative results for the real-world Lei et al.
dataset [18]. We categorize the compared methods into three groups: (i) Non-
learning-based methods [13, 12], (ii) learning-based methods using pre-trained
models because of the lack of training codes [3, 6, 17, 16], and (iii)
learning-based methods re-trained using Lei et al. dataset and provided
training codes [2, 5, 18]. For the non-polarization-based methods of [3, 6, 2,
5], we used a single-channel intensity image (the average of four polarized
images) as the input.

Although the direct outputs of our RP2PN are polarized reflection and
transmission images, we evaluated the results in the averaged intensity domain
to compare RP2PN with other existing methods. Because there are scale
differences in the result images from different methods, we also re-scaled the
result images of all the methods as
$\hat{T}^{\prime}_{i}=\alpha_{i}\hat{T}_{i}$ and
$\hat{R}^{\prime}_{i}=\alpha_{i}\hat{R}_{i}$, where $i$ is the scene index in
the testing dataset. The scaling factor $\alpha_{i}$ was determined for each
scene and each method as
$\alpha_{i}=\overline{I_{i}}/\overline{I^{\prime}_{i}}$, where
$\overline{I_{i}}$ is the mean pixel value of the input image and
$\overline{I^{\prime}_{i}}$ is the mean pixel value of the derived mixed image
of $I^{\prime}_{i}=\hat{T}_{i}+\hat{R}_{i}$. With this re-scaling based on the
same input image’s scale, all methods are more fairly compared.

Fig. 4: Qualitative comparison with existing methods. Fig. 5: Example of
qualitative results on polarization outputs.

The PSNR and SSIM results in Table 1 show that non-learning polarization-based
methods of [13, 12] exhibit low performance, due to their idealized physical
assumptions which are often broken in real-world scenarios. Learning-based
polarization methods of [17, 16] do not achieve the expected performance
because the provided pre-trained models were trained on synthetic datasets and
showed limited generalizability to Lei et al. dataset. The results of Lei et
al. method [18] and our RP2PN demonstrate higher performance than the non-
polarization-based methods of [2, 5] using the same training data, which
validates the effectiveness of using the polarization. Furthermore, our RP2PN
achieves the best PSNR and SSIM results and significant improvement,
especially for the reflection. Figure 4 shows the qualitative results (for
selected competivie methods due to limited space), where details of each
result are shown in green rectangles. Compared with other methods, our
transmission result can recover building walls better, while the reflection
result preserves the clear edges of the stairs. The results for other scenes
can be seen in the supplementary material111Link:
https://github.com/wjbian/RP2PN/blob/main/supp.pdf.

Since our RP2PN provides the polarized outputs, we show one example of these
outputs in Fig. 5. From the results, we can confirm that the polarized
reflection and transmission images, as well as the calculated intensity, AoP,
and DoP images are reasonably close to the ground truths, which demonstrates
that our RP2PN can successfully learn the polarization information of the
reflection and the transmission.

###  3.3 Ablation Study

Table 2 summarizes the ablation study results. In models 1 and 2, we replaced
the inputs of four polarized images with the standard intensity image to
investigate the influence of the polarization input. In models 1 to 3, we
replaced the network outputs from four-channel polarized images to one-channel
intensity image to investigate the effect of the polarization output. In
models 1, 3, and 4, we removed the iteration to investigate the impact of the
recurrent framework. Comparing models 1 and 3, utilizing polarized images as
the inputs significantly improves the performance. Comparing models 3 and 4,
incorporating the polarization output also enhances the separation. Comparing
model 4 and ours, it becomes evident that the incorporation of LSTM iterations
offers a substantial improvement in predictions, particularly for the
reflection. From all of these results, we can confirm that both our
polarization-to-polarization approach and recurrent framework are effective.

Table 2: Ablation study.

Model |  Polar Input |  Polar Output |  With Iteration | Transmission | Reflection
---|---|---|---|---|---
PSNR | SSIM | PSNR | SSIM
1 | No | No | No | 32.93 | 0.928 | 32.52 | 0.892
2 | No | No | Yes | 32.97 | 0.928 | 32.87 | 0.897
3 | Yes | No | No | 34.95 | 0.949 | 34.56 | 0.921
4 | Yes | Yes | No | 35.14 | 0.950 | 34.78 | 0.923
Ours | Yes | Yes | Yes | 35.87 | 0.954 | 35.63 | 0.933

##  4 Conclusion

In this paper, we have proposed a novel recurrent polarization-to-polarization
network, named RP2PN, for reflection removal. Compared with existing
polarization-to-intensity approaches, our RP2PN can better utilize the mutual
polarimetric relationship between the reflection and the transmission by
learning the polarized outputs and incorporating a recurrent framework. The
quantitative and qualitative results have validated that our RP2PN is superior
to other state-of-the-art methods.

## References

  * [1] Qingnan Fan, Jiaolong Yang, Gang Hua, Baoquan Chen, and David Wipf,  “A generic deep architecture for single image reflection removal and image smoothing,”  in Proceedings of IEEE International Conference on Computer Vision (ICCV), 2017, pp. 3238–3247.
  * [2] Xuaner Zhang, Ren Ng, and Qifeng Chen,  “Single image reflection separation with perceptual losses,”  in Proceedings of IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2018, pp. 4786–4794.
  * [3] Jie Yang, Dong Gong, Lingqiao Liu, and Qinfeng Shi,  “Seeing deeply and bidirectionally: A deep learning approach for single image reflection removal,”  in Proceedings of European Conference on Computer Vision (ECCV), 2018, pp. 654–669.
  * [4] Kaixuan Wei, Jiaolong Yang, Ying Fu, Wipf David, and Hua Huang,  “Single image reflection removal exploiting misaligned training data and network enhancements,”  in Proceedings of IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2019, pp. 8170–8179.
  * [5] Chao Li, Yixiao Yang, Kun He, Stephen Lin, and John E Hopcroft,  “Single image reflection removal through cascaded refinement,”  in Proceedings of IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2020, pp. 3565–3574.
  * [6] Zheng Dong, Ke Xu, Yin Yang, Hujun Bao, Weiwei Xu, and Rynson WH Lau,  “Location-aware single image reflection removal,”  in Proceedings of IEEE International Conference on Computer Vision (ICCV), 2021, pp. 5017–5026.
  * [7] Xiaojie Guo, Xiaochun Cao, and Yi Ma,  “Robust separation of reflection from multiple images,”  in Proceedings of IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2014, pp. 2187–2194.
  * [8] Byeong-Ju Han and Jae-Young Sim,  “Reflection removal using low-rank matrix completion,”  in Proceedings of IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2017, pp. 5438–5446.
  * [9] Chengxuan Zhu, Renjie Wan, and Boxin Shi,  “Neural transmitted radiance fields,”  in Proceedings of Advances in Neural Information Processing Systems (NeurIPS), 2022, vol. 35, pp. 38994–39006.
  * [10] Yasushi Maruyama, Takashi Terada, Tomohiro Yamazaki, Yusuke Uesaka, Motoaki Nakamura, Yoshihisa Matoba, Kenta Komori, Yoshiyuki Ohba, Shinichi Arakawa, Yasutaka Hirasawa, Yuhi Kondo, Jun Murayama, Kentaro Akiyama, Yusuke Oike, Shuzo Sato, and Takayuki Ezaki,  “3.2-MP back-illuminated polarization image sensor with four-directional air-gap wire grid and 2.5-$\mu$m pixels,”  IEEE Transactions on Electron Devices, vol. 65, no. 6, pp. 2544–2551, 2018.
  * [11] Miki Morimatsu, Yusuke Monno, Masayuki Tanaka, and Masatoshi Okutomi,  “Monochrome and color polarization demosaicking based on intensity-guided residual interpolation,”  IEEE Sensors Journal, vol. 21, no. 23, pp. 26985–26996, 2021.
  * [12] Hany Farid and Edward H Adelson,  “Separating reflections and lighting using independent components analysis,”  in Proceedings of IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 1999, vol. 1, pp. 262–267.
  * [13] Yoav Y Schechner, Joseph Shamir, and Nahum Kiryati,  “Polarization and statistical analysis of scenes containing a semireflector,”  Journal of the Optical Society of America. A, vol. 17, no. 2, pp. 276–284, 2000.
  * [14] Naejin Kong, Yu-Wing Tai, and Joseph S Shin,  “A physically-based approach to reflection separation: From physical modeling to constrained optimization,”  IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 36, no. 2, pp. 209–221, 2013.
  * [15] Takuma Aizu and Ryo Matsuoka,  “Reflection removal using multiple polarized images with different exposure times,”  in Proceedings of European Signal Processing Conference (EUSIPCO), 2022, pp. 498–502.
  * [16] Patrick Wieschollek, Orazio Gallo, Jinwei Gu, and Jan Kautz,  “Separating reflection and transmission images in the wild,”  in Proceedings of the European Conference on Computer Vision (ECCV), 2018, pp. 90–105.
  * [17] Youwei Lyu, Zhaopeng Cui, Si Li, Marc Pollefeys, and Boxin Shi,  “Reflection separation using a pair of unpolarized and polarized images,”  in Proceedings of Advances in Neural Information Processing Systems (NeurIPS), 2019, pp. 14559–14569.
  * [18] Chenyang Lei, Xuhua Huang, Mengdi Zhang, Qiong Yan, Wenxiu Sun, and Qifeng Chen,  “Polarized reflection removal with perfect alignment in the wild,”  in Proceedings of IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2020, pp. 1747–1755.
  * [19] Karen Simonyan and Andrew Zisserman,  “Very deep convolutional networks for large-scale image recognition,”  arXiv preprint:1409.1556, 2014.
  * [20] Olaf Ronneberger, Philipp Fischer, and Thomas Brox,  “U-net: Convolutional networks for biomedical image segmentation,”  in Proceedings of Medical Image Computing and Computer-Assisted Intervention (MICCAI), 2015, pp. 234–241.
  * [21] Justin Johnson, Alexandre Alahi, and Li Fei-Fei,  “Perceptual losses for real-time style transfer and super-resolution,”  in Proceedings of the European Conference on Computer Vision (ECCV), 2016, pp. 694–711.