Problem 2: Lagrangian Observation Mapping#
Problem Statement#
The CryoSat-2 and ICESat-2 satellites are our primary sources for remotely acquired ice thickness data. While the former has wide coverage, its accuracy is low. Conversely, ICESat-2 has high accuracy but a low footprint. The goal of data fusion is to achieve the best of both worlds: We want high coverage and high accuracy. However, in order to combine these data sets, we must first solve a preliminary problem.
Satellite observations are taken at different times, and the ice topology (i.e., its spatial configuration on the surface of the ocean) will have changed during the interim through a combination of advection as well as various mechanical deformations, such as splitting and ridging. If we ignore the effects of ridging, then, over small time scales, thermodynamic effects will be negligible, and we can assume that the thickness at any particular point on the ice floe or ice sheet in question remains constant.
Considering the movement of particular points of ice gives us a Lagrangian, or particle-based, coordinate system, as opposed to an Eulerian coordinate system, where we measure changes in ice position relative to a fixed background grid. Under the Lagrangian perspective, we can summarize our assumption from above by saying: Thickness remains constant along flow lines. Thus, if we can register the data—that is, if we can undo the effects of the various transformations to obtain a static view of the ice topology—we will have a much larger data set to work with.
More precisely, our problem is this: Let \(\Omega\subseteq\mathbb{R}^2\) be a domain, and let \(\Gamma_t\subseteq\Omega\) be an open subset representing the ice topology at time \(t\). Further, let \(H:\mathbb{R}_{+}\times \Gamma_t\to\mathbb{R}\) be a (random) field denoting ice thickness at time \(t\) and location \((x, y)\). Our goal is to register data of the form \(\{H(t_i,x_i)\}_{i=1}^N\) to obtain a non-time-varying ice thickness field \(\tilde{H}:\Gamma\to\mathbb{R}\) defined on a common reference configuration \(\Gamma\subseteq \Omega\) of the ice.
For simplicity, normalize the time interval to \([0, 1]\), and assume that thickness observations are taken only at initial and final times. That is, assume we have \(N\) thickness measurements \(\{H(0, x_i)\}_{i=1}^N\) at time \(t=0\), and \(M\) thickness measurements \(\{H(1,y_j)\}_{j=1}^M\) at time \(t=1\). In order to chain together sequences of observations, we will typically take \(\Gamma_1\) to be the reference configuration. Thus, we seek a diffeomorphism, or merely a measure-preserving map, \(T:\Gamma_0\to\Gamma_1\). Given such \(T\), we set
Using this data, we can then apply further inference, approximation, or model-based methods (see problems 1 and 2) to define a registered ice thickness field \(\tilde{H}:\Gamma\to\mathbb{R}\) over the entire reference configuration \(\Gamma=\Gamma_1\).
Methodology#
Registration via Optimal Transport (OT)#
We now describe our method for obtaining the map \(T\) used to register the data. As shown by Parno et al. [Parno2019], optimal transport methods are well-suited for solving sea ice motion estimation and registration problems, and it is this methodology that we follow here. There are, of course, other methods one can consider, such as image registration in the large deformation diffeomorphic metric matching (LDDMM) framework (see, e.g., [Beg2005]), or optical flow algorithms (e.g., [Petrou2016]), to name only two. Such methods, however, can struggle with the large displacements and topological changes (namely splitting) that are common in sea ice scenarios.
We refer to [Parno2019] for complete details of the method. Briefly, the idea is this: Suppose we have scalar-valued source and target images \(p\) and \(q\), both of dimension \(N=n\times m\). After vectorizing and normalizing the images, we can view \(p\) and \(q\) as probability mass functions on \(\{1, \dots, N\}\). A transport plan, or coupling, from \(p\) to \(q\) is a joint distribution with marginals given by \(p\) and \(q\), and in the discrete setting, a coupling is given by an \(N\times N\) matrix \(\gamma\) whose rows and columns sum to \(p\) and \(q\), respectively. That is,
where \(p_i\) denotes the \(i\)-th entry of \(p\), and similarly for \(q\). To select the plan that rearranges the source mass into the target mass most efficiently, we define a ground cost \(c_{ij}\) that gives the cost of moving a unit of mass from pixel \(i\) to pixel \(j\). Typically, \(c_{ij}\) is taken as the squared Euclidean distance between the pixel centers. Then, the cost of a transport plan is given by \(\sum_{ij}c_{ij}\gamma_{ij}\), and the optimal transport problem is to find the transport plan \(\gamma\) with minimum cost.
Extensive effort has gone into solving optimal transport problems efficiently; see, e.g., [Peyre2019]. For numerical experiments, we use the solver developed by Zhou and Parno [Zhou2022], which can be used for more general multimarginal optimal transport problems without regularization.
Given the optimal coupling \(\sigma\), we obtain a transport map \(T\) via the barycentric projection map, as defined in [Peyre2019]; see [Parno2019] for a quick summary of details. We can see the results of this method in the synthetic example shown below:
The point in the lower right-hand corner corresponds to a location where we have thickness data available. As can be seen, the transport map accurately tracks the motion of this point as the ice floe undergoes various changes. As such, we can map this and other measurements from the source image into the target image to obtain a registered data set suitable for further processing.
Registration via Multimarginal Optimal Transport (MMOT)#
When we are given a sequential of images:
we are able to use those finer information from the intermediate images to tone the motion, such that images at \(t=0.25\), \(t=0.5\), \(t=0.75\) match with \(\mu_1\), \(\mu_2\), \(\mu_3\) respectively.
The top plot is done via the OT interpolation between \(\mu_0\) and \(\mu_4\). The bottom plot is done via the MMOT interpolation that satisfies all marginals \(\{\mu_i\}_{i=0}^4\). As the OT-based interpolation follows constant-speed geodesics, the interpolation at \(t=0.125\) for example, may not capture that the motion is slow from marginal \(\mu_0\) at \(t=0.0\) to marginal \(\mu_1\) at \(t=0.25\). On the opposite, this feature is captured via MMOT interpolation.
References#
- Beg2005
Beg, M. F., Miller, M. I., Trouvé, A., & Younes, L. (2005). Computing large deformation metric mappings via geodesic flows of diffeomorphisms. International journal of computer vision, 61(2), 139-157.
- Parno2019(1,2,3)
Parno, M. D., West, B. A., Song, A. J., Hodgdon, T. S., & O’Connor, D. T. (2019). Remote measurement of sea ice dynamics with regularized optimal transport. Geophysical Research Letters, 46(10), 5341-5350.
- Petrou2016
Petrou, Z. I., & Tian, Y. (2016). High-resolution sea ice motion estimation with optical flow using satellite spectroradiometer data. IEEE Transactions on Geoscience and Remote Sensing, 55(3), 1339-1350.
- Peyre2019(1,2)
Peyré, G., & Cuturi, M. (2019). Computational optimal transport: With applications to data science. Foundations and Trends® in Machine Learning, 11(5-6), 355-607.
- Zhou2022
Zhou, B., & Parno, M. (2022). Efficient and Exact Multimarginal Optimal Transport with Pairwise Costs. arXiv preprint arXiv:2208.03025.