Online Estimation of Ocean Current from Sparse GPS Data for Underwater Vehicles

TL;DR

This paper introduces a Gaussian process-based EM algorithm that estimates ocean currents from sparse GPS data and dead-reckoned positions, improving underwater navigation by leveraging position drift.

Contribution

The paper presents a novel GP-EM algorithm that exploits ocean current incompressibility to accurately estimate currents from limited GPS data for underwater vehicles.

Findings

Successfully reconstructs ocean current fields in simulation

Validates the approach on real GPS datasets

Demonstrates potential for improved underwater navigation

Abstract

Underwater robots are subject to position drift due to the effect of ocean currents and the lack of accurate localisation while submerged. We are interested in exploiting such position drift to estimate the ocean current in the surrounding area, thereby assisting navigation and planning. We present a Gaussian process~(GP)-based expectation-maximisation~(EM) algorithm that estimates the underlying ocean current using sparse GPS data obtained on the surface and dead-reckoned position estimates. We first develop a specialised GP regression scheme that exploits the incompressibility of ocean currents to counteract the underdetermined nature of the problem. We then use the proposed regression scheme in an EM algorithm that estimates the best-fitting ocean current in between each GPS fix. The proposed algorithm is validated in simulation and on a real dataset, and is shown to be capable of reconstructing the underlying ocean current field. We expect to use this algorithm to close the loop between planning and estimation for underwater navigation in unknown ocean currents.

Equations

$${\dot{\mathbf{x}}}_t={\mathbf{v}}_t+\mathbf{w}({\mathbf{x}}_t),$$

$${\mathbf{x}}_{t+1}={\mathbf{x}}_t+({\mathbf{v}}_t+\mathbf{w}({\mathbf{x}}_t))\Delta t,$$

$${\hat{\mathbf{x}}}_{t+1}={\hat{\mathbf{x}}}_t+{\mathbf{v}}_t\Delta t.$$

$${\mathbf{y}}_{{\tau }_k}$$

$${ϵ}_k$$

$$\Delta {\mathbf{x}}_k={\mathbf{y}}_{{\tau }_k}-{\hat{\mathbf{x}}}_{{\tau }_k}.$$

$${\mathbf{w}}^{\ast }(\mathbf{x})=\underset{\mathbf{w}(\mathbf{x})\in {\mathcal{C}}^{\infty }({\mathbb{R}}^2)}{\mathrm{argmax}}\mathcal{P}(\mathbf{w}(\mathbf{x})\mid {\mathbf{C}}_{1:k}).$$

$$\mathcal{P}(\mathbf{w}(\mathbf{x})\mid {\mathbf{C}}_{1:k})=\int \mathcal{P}(\mathbf{w}(\mathbf{x})\mid {\mathbf{W}}_{1:k})\mathcal{P}({\mathbf{W}}_{1:k}\mid {\mathbf{C}}_{1:k})d{\mathbf{W}}_{1:k},$$

$$\nabla \cdot \mathbf{w}(\mathbf{x})=0,$$

$$\mathbf{w}(\mathbf{x})={\left[\begin{array}{cc}\frac{\partial \varphi (\mathbf{x})}{\delta y}& -\frac{\partial \varphi (\mathbf{x})}{\delta x}\end{array}\right]}^T.$$

$$\varphi (\mathbf{x})\sim GP(0,k(\Vert \mathbf{x}-{\mathbf{x}}^{\prime }\Vert )),$$

$$\mathbf{w}(\mathbf{x})$$

$$\mathbf{K}(\mathbf{x},{\mathbf{x}}^{\prime })$$

$$=\left[\begin{array}{cc}\frac{{\partial }^{2}k}{\partial y^2}& -\frac{{\partial }^{2}k}{\partial x\partial y}\\ -\frac{{\partial }^{2}k}{\partial x\partial y}& \frac{{\partial }^{2}k}{\partial x^2}\end{array}\right].$$

$$\mathcal{P}(\mathbf{W}({\mathbf{X}}^Q)\mid \mathbf{W}({\mathbf{X}}^D))=\mathcal{N}(\mu ({\mathbf{X}}^Q),\mathbf{\Sigma }({\mathbf{X}}^Q)),$$

$$\mu ({\mathbf{X}}^Q)={\mathbf{K}}_{DQ}^T{K}_{DD}^{-1}{\mathbf{W}}_D$$

$$\mathbf{\Sigma }({\mathbf{X}}^Q)={\mathbf{K}}_{QQ}-{\mathbf{K}}_{DQ}^T{K}_{DD}^{-1}{\mathbf{K}}_{DQ},$$

$$\mathcal{P}({\mathbf{W}}_{1:k}\mid {\mathbf{C}}_{1:k})$$

$$=$$

$$\mathcal{P}({\mathbf{W}}_k\mid {\mathbf{C}}_k,{\mathbf{W}}_{1:k-1})$$

$$\propto$$

$$=$$

$${\mathbf{W}}_k^i=\underset{\mathbf{W}}{\mathrm{argmax}}{\mathbb{E}}_{{\mathbf{X}}_k^i\mid {\mathbf{W}}_k^{i-1},{\mathbf{C}}_k}[\mathcal{P}(\mathbf{W}\mid {\mathbf{X}}_k^i,{\mathbf{C}}_k,{\mathbf{W}}_{1:k-1})].$$

$$\mathcal{P}({\mathbf{X}}_k^i\mid {\mathbf{W}}_k^{i-1},{\mathbf{C}}_k)=\delta ({\mathbf{X}}_k^i-({\hat{\mathbf{X}}}_k+B{\mathbf{W}}_k^{i-1})),$$

$${\mathbf{X}}_k^i={\hat{\mathbf{X}}}_k+B{\mathbf{W}}_k^{i-1}.$$

$${\mathbb{E}}_{{\mathbf{X}}_k^i\mid {\mathbf{W}}_k^{i-1},{\mathbf{C}}_k}\left[\mathcal{P}({\mathbf{W}}_k^i\mid {\mathbf{X}}_k^i,{\mathbf{W}}_{1:k-1})\right]=\mathcal{N}(\mu ({\mathbf{X}}_k^i),\mathbf{\Sigma }({\mathbf{X}}_k^i)),$$

$$\Delta {\mathbf{x}}_k=C{\mathbf{W}}_k+{ϵ}_k,$$

$$\begin{gathered} {\mathbb{E}}_{{\mathbf{X}}_k^i\mid {\mathbf{W}}_k^{i-1},{\mathbf{C}}_k}\left[\mathcal{P}({\mathbf{W}}_k^i,\Delta {\mathbf{x}}_k\mid {\mathbf{X}}_k^i,{\mathbf{W}}_{1:k-1})\right] \\ =\mathcal{N}\left(\left[\begin{array}{c}\mu ({\mathbf{X}}_k^i)\\ C\mu ({\mathbf{X}}_k^i)\end{array}\right],\left[\begin{array}{cc}\Sigma ({\mathbf{X}}_k^i)& \Sigma ({\mathbf{X}}_k^i)C^T\\ C\Sigma ({\mathbf{X}}_k^i)& C\Sigma ({\mathbf{X}}_k^i)C^T+{\sigma }_y^2{I}_{2\times 2}\end{array}\right]\right). \end{gathered}$$

$$\begin{array}{r}{\mathbf{W}}_k^i=\mu +\Sigma C^T{\left(C\Sigma C^T+{\sigma }_y^{2}I\right)}^{-1}\left(\Delta {\mathbf{x}}_{T_k}-C\mu \right)\end{array},$$

Full text

I INTRODUCTION

Ocean monitoring offers tremendous economic value with various applications such as oceanographic research [1], military surveillance [2], and oil and gas source localisation [3]. Various autonomous platforms have been used in such ocean monitoring tasks, including autonomous underwater vehicles (AUVs) [4], underwater gliders [5], and even passive platforms without any actuation [6].

A main challenge in underwater robotics is the effect of ambient ocean current. Due to lack of GPS while submerged, ocean currents can cause considerable position drift. This has a strong impact on the utility of the gathered data, as most ocean monitoring tasks concern spatial phenomena [7]. Thus, there has been substantial work on navigation and planning in flow fields [8, 9, 10, 11]. However, most of this work assumes that the flow field is given a priori, e.g. from an external database [12, 13, 14]. Unfortunately, the spatiotemporal resolution or accuracy provided by most databases is insufficient for the purpose of navigation, as noted in [10].

We are interested in estimating the ocean current online, without any prior information. This idea was recently explored as a deterministic, discrete estimation problem based on a prior oceanic model [15] augmented by the drift of the vehicle. Here, we also exploit drift as a source of information, but focus on a continuous, probabilistic form. Drift is measured between the true position, measured by GPS at the surface, and the dead-reckoned trajectory. As GPS is unavailable underwater, the dead-reckoning estimate drifts substantially from the true trajectory. In the case of underwater gliders, this disparity is typically in the order of a few hundred metres, as shown in Figs. 1 and 2.

We propose an expectation-maximisation (EM) algorithm for estimating ocean current given GPS measurements and a dead-reckoned trajectory. The problem is severely underdetermined as the position drift is the sum of current along the trajectory. We present a Gaussian process (GP) regression technique that incorporates the concept of incompressibilty to provide a physically meaningful constraint. Our algorithm is demonstrated both in simulation and using two experimental datasets collected by underwater gliders. The simulation results show that the algorithm is capable of estimating accurately, starting from a uniform prior on current. With the experimental datasets, we could not compare to a baseline due to lack of other sources of data at the site, but the estimated current aligns with observed drift and also aligns with the shoreline, which is an expected pattern. The significance of this result is that we can now attempt to close the loop between estimation and planning, hence enabling underwater navigation in unknown ocean current.

II RELATED WORK

To estimate ocean currents online, one approach is to consider current as a low-frequency disturbance and then apply an extended Kalman filter (EKF) [16] or nonlinear observer [17] in conjunction with acoustic sensors. However, modelling current as a temporal phenomenon clearly overlooks its spatial structure, and acoustic sensors typically require a stationary reference (e.g., the seabed) [18].

An approach that does consider the spatial nature of the problem is presented in [19]. The authors examine the feasibility of ocean current estimation through simply calculating the average current velocity by dividing the position drift by time. Unsurprisingly, the estimate is increasingly unreliable as distance between diving and surfacing locations grows, and no predictive capability is provided.

An improvement to this concept was presented in [15]. The work proposes the ‘motion tomography’ algorithm, which reconstructs the local ocean current from GPS measurements using techniques from the computer tomography (CT) literature. The algorithm represents the ocean as a discrete grid and iteratively solves a linear system based on a prior obtained from an oceanic model, such as [14]. The authors in [20] develop a similar technique with dense GPS measurements. In comparison, the algorithm we propose in this paper does not require prior information other than sparse GPS measurements. Because we use a GP, we can incorporate prior current measurements at any location, if available, and estimate current at any location with associated uncertainty. Further, our proposed method accounts for the incompressibility of the ocean current, which not only enhances the physical fidelity, but also aids with the quality of the estimates as will be shown later. Using a GP also aligns better with planning, as evident in previous work such as [21, 10, 8]. The techniques we use in this paper are in line with the Bayesian system identification literature [22], and are inspired by its rigorous theoretical analysis.

III PROBLEM FORMULATION

Suppose we have a continuous-time dynamic model of an underwater vehicle:

[TABLE]

where $\mathbf{x}_{t}\in\mathbb{R}^{2}$ is the position of the vehicle at time $t$, $\mathbf{v}_{t}\in\mathbb{R}^{2}$ is the velocity through water at time $t$, $\mathbf{w}\in\mathcal{C}^{\infty}(\mathbb{R}^{2})$ is the ocean current modelled as a smooth 2D vector field. For simplicity, we do not take into account the vertical motion of the vehicle. The continuous-time model is discretised as:

[TABLE]

where $\Delta t$ is the sampling time.

The vehicle’s velocity through water $\mathbf{v}_{t}$ is known, whereas the current $\mathbf{w}(\mathbf{x})$ is unknown. The vehicle estimates its own position based on dead-reckoning assuming zero ocean current while submerged. The initial estimate derives from the last GPS measurement prior to dive-in. Namely,

[TABLE]

The vehicle attempts to reach a target point using its dead-reckoned position estimate. When the vehicle’s estimate is within a pre-set tolerance range from the target point, the vehicle climbs up to the surface and updates its position using GPS. We denote the time of surfacing events as $\tau_{k}\in[1,T]$, $\tau_{k}<\tau_{k+1}$. The GPS measurements are assumed to have i.i.d. Gaussian measurement error. Namely,

[TABLE]

For the periods in between each GPS measurement, we use shorthand notation $\mathbf{X}_{k}=\mathbf{X}_{\tau_{k-1}:\tau_{k}}$ for true trajectory, $\mathbf{\hat{X}}_{k}=\mathbf{\hat{X}}_{\tau_{k-1}:\tau_{k}}$ for dead-reckoned trajectory, $\mathbf{W}_{k}=\{\mathbf{w}(\mathbf{x}_{\tau_{k-1}}),\ldots\mathbf{w}(\mathbf{x}_{\tau_{k}})\}$ for current along trajectory, and $\mathbf{y}_{k}=\mathbf{y}_{\tau_{k}}$ for GPS measurements respectively.

As the dead-reckoned estimates do not take ocean current into account, there is a substantial disparity between the dead-reckoned estimate and the GPS measurement. Throughout the rest of the paper, we refer to this disparity as drift, and denote it by $\Delta\mathbf{x}_{k}$. In other words:

[TABLE]

We will use the concept of a cycle to describe the three behaviours: 1) dive-in, 2) manoeuvre and 3) surfacing, as depicted in Fig. 2. Namely, a cycle $\mathbf{c}_{k}=\{\Delta\mathbf{x}_{k},\mathbf{\hat{X}}_{k}\}$ is a tuple containing the dead-reckoned trajectory, $\mathbf{\hat{X}}_{k}$, and the measured drift $\Delta\mathbf{x}_{k}$.

Although the current $\mathbf{w}(\mathbf{x})$ throughout each cycle is unknown, we know that the drift measurements are related to the current. The ultimate aim of this paper is to solve the following maximum a posteriori (MAP) estimation problem.

Problem 1** (Ocean current estimation).**

Suppose we have a sequence of cycles $\mathbf{C}_{1:k}=\mathbf{c}_{1}\mathbf{c}_{2}\cdots\mathbf{c}_{k}$. Find an optimal estimate for ocean current $\mathbf{w}^{*}(\mathbf{x})$ over the space of 2D smooth vector fields $\mathcal{C}^{\infty}(\mathbb{R}^{2})$ that maximises the posterior probability:

[TABLE]

Intuitively, solving the MAP problem implies that we find the ocean current $\mathbf{w}(\mathbf{x})$ that is the best trade-off between fitting 1) the drift measurements, $\Delta\mathbf{x}_{k}$ and 2) a constraint on the general behaviour of ocean current, which we will discuss in Sec. IV. The constraint is necessary because there are infinitely many possibilities of ocean current vectors that sum up to the drift measurement $\Delta\mathbf{x}_{k}$. In other words, the problem is underdetermined.

Finding a direct solution to Problem 1 is difficult because there can be infinitely many relations between drift measurements $\Delta\mathbf{x}_{k}$ and ocean current vectors $\mathbf{w}(\mathbf{x})$ over the trajectory of vehicle $\mathbf{x}_{t}$. In this paper, we make the following assumption about how current vectors at different positions $\mathbf{w}(\mathbf{x})$ and the drift measurements are related.

Assumption 1** (Conditional independence).**

For all $\mathbf{x}\in\mathbb{R}^{2}$ such that $\mathbf{x}\neq\mathbf{x}_{t}$, $\mathbf{w}(\mathbf{x})$ is conditionally independent of $\mathbf{C}_{1:k}$ given the current along trajectory, $\mathbf{W}_{1:k}=\{\mathbf{W}_{1},\cdots,\mathbf{W}_{k}\}$. In other words, $\mathbf{w}(\mathbf{x})$ is indirectly related to $\mathbf{C}_{1:k}$ through $\mathbf{W}_{1:k}$.

With the assumption, the overall problem can be re-written in a form that reveals two sub-problems:

[TABLE]

where the sub-problems are to 1) estimate the current at a remote location given the current along the trajectory (i.e., $\mathcal{P}(\mathbf{w}(\mathbf{x})\mid\mathbf{W}_{1:k})$), and 2) estimate the current experienced along the vehicle’s trajectory given drift measurements (i.e., $\mathcal{P}(\mathbf{W}_{1:k}\mid\mathbf{C}_{1:k})$).

IV GP REGRESSION OF INCOMPRESSIBLE FLOW FIELDS

In this section, we solve for Subproblem 1, where we estimate the oceanic flow at a query position given flow at other locations (i.e., $\mathcal{P}(\mathbf{w}(\mathbf{x})\mid\mathbf{W}_{1:k})$). We model the oceanic flow with a GP, and impose the assumption of incompressibility. We first introduce incompressibility and the concept of a streamfunction, then exploit the properties of the streamfunction to derive an incompressible GP. Incompressibility also serves as a useful constraint for solving Subproblem 2. We demonstrate an example with real ocean dataset to illustrate that the incompressible GP outperforms the standard for modelling ocean currents.

IV-A Incompressibility and Streamfunction

In this work, we model the ocean current as a planar, time-invariant and incompressible flow field. Planarity implies that the ocean current has no $z$-component, which describes the horizontal stratification of oceanic flow well (see, e.g., [23]).

A flow field is incompressible [24] when

[TABLE]

where $\nabla\cdot$ is the divergence operator. Intuitively, incompressibility implies that ‘the amount of water coming into a point is equal to the amount exiting the area’. As we do not expect to see surplus or deficit of water entering an area in the ocean, incompressibility is an appropriate description.

If a planar flow is incompressible, it can be represented by a streamfunction $\phi:\mathbb{R}^{2}\rightarrow\mathbb{R}$. Given a streamfunction $\phi(\mathbf{x})$, one can compute the current $\mathbf{w}(\mathbf{x})$ as:

[TABLE]

IV-B Streamfunction-based GP Representation

In this section, we show how to enforce the incompressibility condition in a GP model using a streamfunction. First, consider a streamfunction modelled as a GP:

[TABLE]

where $k(\|\mathbf{x}-\mathbf{x^{\prime}}\|)$ is a kernel function.

Because the derivative of a streamfunction $\phi$ is flow field $\mathbf{w}$ as shown in (10) and the derivative of a GP is another GP [25], our flow field can be represented by a GP. In infinite-dimensional Bayesian estimation, derivative operators apply to functions as do matrices to vectors [26]. Recall that if $Cov(\mathbf{A})=\Sigma_{\mathbf{A}}$ for a vector-valued random variable $\mathbf{A}$, $Cov(M\mathbf{A})=M\Sigma_{\mathbf{A}}M^{T}$ given a matrix $M$.

The derivative operators can be written as $\mathcal{D}=\begin{bmatrix}\frac{\partial}{\partial y}&-\frac{\partial}{\partial x}\end{bmatrix}^{T}$ and $\mathcal{D}^{\prime}=\begin{bmatrix}\frac{\partial}{\partial y^{\prime}}&-\frac{\partial}{\partial x^{\prime}}\end{bmatrix}$ for the function case, and the flow field is represented by:

[TABLE]

where the kernel function $\mathbf{K}$ is given by:

[TABLE]

From the first line to the second line in (13), we used the stationarity of the kernel. It is important to note that (13) can be computed analytically given a choice of kernel for the streamfunction [25, 27].

Using the GP representation of flow field $\mathbf{w}$ with the kernel function in (13), we can predict a set of current vectors $\mathbf{W}(\mathbf{X}^{Q})=\begin{bmatrix}\mathbf{w}(\mathbf{x}_{1}^{Q})&\ldots&\mathbf{w}(\mathbf{x}_{N}^{Q})\end{bmatrix}$ given previous measurement data $\mathbf{W}(\mathbf{X}^{D})=\begin{bmatrix}\mathbf{w}(\mathbf{x}_{1}^{D})&\ldots&\mathbf{w}(\mathbf{x}_{M}^{D})\end{bmatrix}$. The predictions are given as a set of normal random variables:

[TABLE]

with mean and covariance:

[TABLE]

where the matrices $\mathbf{K}_{DD}^{(i,j)}=\begin{bmatrix}\mathbf{K}(\mathbf{x}_{i}^{D},\mathbf{x}_{j}^{D})\end{bmatrix}$, $\mathbf{K}_{DQ}^{(i,j)}=\begin{bmatrix}\mathbf{K}(\mathbf{x}_{i}^{D},\mathbf{x}_{j}^{Q})\end{bmatrix}$, and $\mathbf{K}_{QQ}^{(i,j)}=\begin{bmatrix}\mathbf{K}(\mathbf{x}_{i}^{Q},\mathbf{x}_{j}^{Q})\end{bmatrix}$ are constructed blockwise.

IV-C Comparison with Standard GP

We discuss how the incompressible GP is a better representation of the oceanic currents with a real dataset in Fig. 3. We selected a representative eddy from the east Australian current data provided by the Australian Bureau of Meteorology. Then, we selected training samples along a line to emulate the current estimated along a trajectory. These training samples were extrapolated with GPs having the proposed incompressible kernel and the standard kernel $\mathbf{K}_{SE}=\mathbf{diag}(k_{SE}(\mathbf{x},\mathbf{x}^{\prime}),k_{SE}(\mathbf{x},\mathbf{x}^{\prime}))$.

The standard kernel only fits a smooth vector field to the training samples. Meanwhile, the proposed incompressible GP extrapolates the ocean current much more accurately even with the limited training samples. An apparent benefit is that we can reconstruct eddy-like patterns [13], which leads to a better extrapolation for the flow along a future trajectory given the estimate along the present trajectory.

V EXPECTATION-MAXIMISATION FOR CURRENT ESTIMATION

In this section, we solve Subproblem 2, which concerns estimating the current along trajectory. For simplicity, we will focus on estimating the flow along trajectory sequentially, given each incoming GPS measurements. In doing so, we are making a Markov assumption, where we fix the estimate of current along previous trajectories, $\mathbf{W}_{1:k-1}$, when estimating the current along trajectory, $\mathbf{W}_{k}$. It substantially reduces the computational complexity of the problem, as the algorithm is incremental. More precisely, we assume:

[TABLE]

which shows the problem reduces to estimating the current along present trajectory, given GPS measurements and previous current estimates (i.e. $\mathcal{P}(\mathbf{W}_{k}\mid\mathbf{C}_{k},\mathbf{W}_{1:k-1})$).

The main challenge in estimating the current along trajectory arises from the strong causality between the trajectory itself and the current along the trajectory. To predict the current along trajectory, we must know the trajectory beforehand, and to predict the trajectory, we must know the current along trajectory. The challenge is solved through an EM algorithm. In each iteration, $i$, the EM algorithm iterates over estimating the trajectory $\mathbf{X}_{k}^{i}$, called the expectation step (E-step) and 6), and estimating the current along trajectory $\mathbf{W}_{k}^{i}$ given $\mathbf{X}_{k}^{i}$, called the maximisation step (M-step). In effect, we iteratively ‘guess’ the true trajectory, estimate the flow using the guess, and refine the guess on trajectory using the estimated flow.

Subproblem 2 is re-written in the EM formulation:

[TABLE]

In the E-step (Sec. V-A), we find $\mathcal{P}(\mathbf{X}_{k}^{i}\mid\mathbf{W}_{k}^{i-1},\mathbf{C}_{k})$ and evaluate the expectation (18). In the M-step (Sec. V-B, we find

[TABLE]

V-A E-step

In the E-step, we need to compute the expectation in (18). In (2), notice that the only source of uncertainty derives from the fact that $\mathbf{w}(\mathbf{x})$ is a GP. Therefore, given current, the trajectory is fully known. Formally, the conditional distribution of the trajectory becomes a Dirac delta distribution:

[TABLE]

where $B^{(i,j)}=\begin{bmatrix}\Delta t\mathbf{1}_{2\times 2}\end{bmatrix}$ if $i\leq j$, and $B^{(i,j)}=\begin{bmatrix}\mathbf{0}_{2\times 2}\end{bmatrix}$ otherwise. $\mathbf{1}_{2\times 2}$ and $\mathbf{0}_{2\times 2}$ denote identity and zero matrices.

As such, the expectation integral collapses to a mere evaluation at:

[TABLE]

V-B M-step

In the M-step, we maximise the expectation taken in Sec. V-A. As shown in Sec. V-A, the expectation (18) is a simple evaluation at (21). As $\mathbf{w}(\mathbf{x})$ is a GP, this yields a prior given by a Gaussian random variable:

[TABLE]

where $\mathbf{\mu}(\mathbf{X}_{k}^{i})$ and $\mathbf{\Sigma}(\mathbf{X}_{k}^{i})$ are calculated using GP prediction equations (15) and (16) given the current estimated with previous drift measurements, $\mathbf{W}_{1:k-1}$.

Notice that we can write $\Delta\mathbf{x}_{k}$ in terms of $\mathbf{W}_{k}$ as:

[TABLE]

where $C=\Delta t\begin{bmatrix}\mathbf{1}_{2\times 2}&\mathbf{1}_{2\times 2}&\ldots&\mathbf{1}_{2\times 2}\end{bmatrix}$. Thus, $\mathbf{W}_{k}^{i}$ and $\Delta\mathbf{x}_{k}$ are joint normal random variables:

[TABLE]

The merit of this formulation is that the maximising solution is now given in closed form, because maximising the posterior (19) is equivalent to finding the conditional mean using (24). This is given by [28]:

[TABLE]

where $\mathbf{\mu}=\mathbf{\mu}(\mathbf{X}_{k}^{i})$ and $\Sigma=\Sigma(\mathbf{X}_{k}^{i})$.

V-C Implementation

The algorithm for solving Subproblem 2 is shown in Alg. 1. We initialise the algorithm with a zero-mean GP without any measurements (Alg. 1 line 1). From lines 4 to 7, we convert each incoming GPS measurement into a pair of true trajectory and the current along trajectory, $\mathbf{X}_{k}$ and $\mathbf{W}_{k}$ using our EM algorithm. Afterwards, the current along trajectory is added to the measurement set of the GP, as ‘pseudo-target’ [22]. As more measurements become available, the GP produces better prior for the iteration. We found that the Markov assumption reduces the computation time substantially with minimal performance sacrifice.

VI RESULTS

In this section, we present case studies of our algorithm with simulated example and real dataset collected from Teledyne Webb G3 Slocum gliders. We used the squared exponential kernel [29] for the streamfunction, having lengthscale of $\ell=35km$, and self-variance of $\sigma_{W}^{2}=0.5m^{2}s^{-2}$. The kernel for the current vectors (13) was computed analytically by taking the partial derivatives [25, 27]. The hyperparameters were found by maximising the data likelihood of the dataset from the Australian Bureau of Meteorology using a standard hyperparameter learning procedure [29].

VI-A Simulation Results

We simulated an underwater vehicle travelling in a flow field, with true and dead-reckoned estimates evolving as (2) and (3). The vehicle is given four waypoints, and surfaces when its dead-reckoned position estimate obtained using (3) is within 100 metres from the current target waypoint. For the purpose of validation, we used a double-gyre model and ran the proposed algorithm. The results are shown in Fig. 4.

From Figs. 4a to 4c, it can be seen that the algorithm actively improves the estimate of current as the mission progresses. In Fig. 4a, it can be seen that the estimate after only one cycle is as good as the average current method in [19]. However, by the fourth cycle, it can be seen that the estimated and the true flow fields are already in good agreement, with minor disparity. By the eighth cycle, the estimated and the true flow fields are almost indistinguishable. The uncertainty of the estimated current also decreases.

In order to examine the convergence of the proposed algorithm further, we performed a Monte Carlo simulation with 100 randomly generated double gyre patterns. We took the error between predicted and true currents on a selected grid, normalised by the magnitude of the true current. The convergence was also compared with a standard GP with kernel function $\mathbf{K}(\mathbf{x},\mathbf{x}_{{}^{\prime}})=\mathbf{diag}(k_{SE}(\mathbf{x},\mathbf{x}_{{}^{\prime}}),k_{SE}(\mathbf{x},\mathbf{x}_{{}^{\prime}}))$, having identical parameters (i.e. no incompressibility assumption or streamfunction). The result is shown in Fig. 5.

Figure 5 shows that the algorithm gradually learns any randomly generated flow field, which is demonstrated by the decrease in normalised error for both standard and the proposed incompressible GP. However, the incompressible GP shows a much faster rate of convergence and a lower steady-state error than the standard. The result clearly indicates that our incompressible GP outperforms over the standard in describing oceanic flows as shown in Fig. 3.

VI-B Field Results

The proposed algorithm was tested in two field trials with a Slocum G3 underwater glider, one near Jervis Bay, Australia, and one in open ocean near Perth, Australia. The result from the Jervis Bay trial is shown in Fig. 6a, and the Perth trial in Fig. 6b. The glider was tasked to visit designated waypoints and communicate the current cycle $\mathbf{c}_{k}$ when on surface. For the purpose of experiment, we disabled the onboard average current correction [19] to ensure the dead-reckoned estimate evolve as (3).

For the Jervis Bay trial in Fig. 6a, it can be seen that the current estimated in Fig. 6a is in good agreement with what is expected near a bay: the ocean flow comes in from the open ocean, and majority of the flow enters the bay in alignment with the bay’s shoreline. An important observation is that there exists only a small perpendicular ocean current to the land which well satisfies our common intuition. The observation clearly indicates that our method accurately models the correlation between ocean currents at different points. This is because the method incorporates incompressibility, a physical attribute of the real ocean, unlike the standard GP-based model.

For the Perth trial in Fig. 6b, the predicted ocean current for Fig. 2 within the operating zone had only one data point, pointing 45$\degree$ southeast. However, the prediction does not seem to agree with the drift observed. Meanwhile, the proposed algorithm can estimate a flow field that best explains the observed drift.

VII Conclusion and Future Work

We have developed an EM algorithm for estimating ocean current from sparse GPS data. The performance of the algorithm is improved by using a GP regression scheme that takes into account the physical intuition of incompressibility. The proposed algorithm was tested with both simulated and experimental data, with positive results. Our future work includes extending the algorithm to consider time- and depth-variation of ocean current, by adding extra dimensions to the kernel function. We would also like to combine the proposed algorithm with motion and task planning algorithms [30, 11], and validate experimentally the performance of the whole framework.