# Dynamic behavior and its typical time of elevators under the random inflow
of passengers

Sakurako Tanida  Graduate School of Science, the University of Tokyo, 7-3-1
Hongo, Bunkyo-ku, Tokyo, Japan  Research Center for Advanced Science and
Technology The University of Tokyo, 4-6-1 Komaba, Meguro-ku, Tokyo, Japan

###### Abstract

Elevators can be regarded as oscillators driven by the calls of passengers who
arrive randomly. We study the dynamic behavior of elevators during the down
peak period by simulation and analytical calculation. We assume that new
passengers arrive according to a Poisson process on each floor and call the
elevators to go down to the ground floor. We numerically examine how the round
trip time depends on the inflow rate of passengers on each floor and reproduce
it by an equation considering the combination of floors where call occurs. We
also examine the degree of the synchronization of two elevators by setting an
order parameter.

##  1 Introduction

Congestion is a notorious problem in our daily lives, but at the same time, it
shows interesting physical phenomena as the out-of-equilibrium system [1, 2,
3]. For example, the traffic jam was empirically known to emerge even in
normal conditions, without disruptions cause by traffic accidents or
bottlenecks. The quantitative measurements revealed the jamming transition
occurs as the vehicle density increases [4, 5]. The traffic flow linearly
increases with the density at low densities, while it varies widely at high
densities. In the intermediate densities, the flow does not depend uniquely on
the density, which implies the existence of the hysteresis effect and meta-
stable states. The traffic jam is also known for showing density waves and
proposed various models to describe it [6, 7, 8, 9, 10]. Such out-of-
equilibrium behavior are often observed in the transportation systems, where
particles moving around by consume energy interact with each other.

Elevators are another example of transportation systems exhibiting nontrivial
out-of-equilibrium behaviors. The multiple elevators move as clusters in the
crowded case, which seem similar to traffic jams [11, 12, 13, 14, 15, 16, 17,
18, 19, 20]. However, the elevators exhibit different spontaneous ordering
from the traffic jam as the out-of-equilibrium oscillators. Elevator and car
systems are also different in the interaction: elevators freely pass other
elevators, but cars cannot go ahead if other cars exist in front of them. Such
difference can be simply interpreted as the difference in the volume exclusion
effect of vehicles. Since elevators interact with each other through the
existence of waiting passengers instead of the volume exclusion effect, it was
found that the clustering of elevators emerges by increasing rather the mean
number of passengers coming in a unit time, or the inflow, than the vehicle
density [11, 12, 13, 14, 15, 16, 17, 18, 19, 20]. The dependence of dynamics
on the inflow is reported in various scenarios, such as during down- and up-
peak periods [11, 12, 13, 14, 15, 16, 17, 18, 19, 20]. Recent studies propose
stochastic models and show the typical time of round trips in steady-state at
each inflow [14, 15, 16, 17, 18, 19]. Those models assume that elevators’
motion is a limit cycle by setting the elevators necessarily stop at both the
highest and ground floors. However, the elevator system is a system whose
motion is essentially caused by the Poisson noise of passengers’ arrivals.
Thus, it is significant to investigate the ordering of elevators spontaneously
emerging from the noise to understand the mechanism of the clustering of
elevators.

In this paper, we investigate elevators’ motion during a down peak period when
passengers go down to the ground floor to exit from the building. Note that we
do not focus on efficiency in this paper. We consider the down peak situation
that new passengers on each floor arrive according to a Poisson process. Since
the inflow rate of passengers on each floor is uniform, and elevators do not
necessarily serve the highest floor. We simulate a stochastic model and show
that the ordering of elevators spontaneously appears from the noise by setting
an order parameter. We also acquire the average time of a round trip at each
inflow rate of passengers. To explain the dependency of the time of a round
trip on the inflow rate, we make a self-consistent equation considering the
combination of calling floors during a round trip.

##  2 Model

Figure 1: Schematic illustration of this model.

The transportation systems of one or two elevators serving $K$ floors of a
building have been investigated. The elevators take 1 timestep to move 1 floor
up or down and $\gamma$ timesteps for passengers to enter or exit. As employed
in previous works [11, 12], we assume that all calls for elevators are from
passengers waiting for the elevators on floor $1\sim K$ in order to move to
the ground floor and exit the building. Once an elevator goes down, it does
not go up again until it arrives at the ground floor. If there are no
passengers in the elevators nor floors, the elevators stay at the floor where
they stopped until the next call. In the case that the next call is single and
there is no other unresolved call and no passenger in elevators, the closest
elevator is set to go to the calling floor. In the case that number of
unresolved calls is more than one and both elevators are free, the higher
elevator moves toward the highest floor. And the lower elevator also moves if
thereare any floors where the lower one is estimated to arrive faster than the
higher one. To estimate the time to go to the target floor $k_{t}$ from the
floor where the eleveter is $k$, we assume following things: An elevator in
going up is assumed to stop at the highest floor among those having unsolved
call $k_{h}$ and then switch to go down and stop at all floors between $k_{t}$
and the highest floor. If an elevator in going down and $k>k_{t}$, it is
assumed to stop at the all floors where unsolved callings occure between $k$
and $k_{t}$. If it in going down and if $k_{t}\leq k$, it is assumed to stop
at the all floors where unsolved callings occure between $k$ and the ground
floor and between $k_{h}$ and $k_{t}$. If an elevator is stopping at the
ground floor, it is assumed to move toward $k_{t}$ without stopping at other
floors. This estimated time is also used to decide whether the lower elevator
goes up or stops at a floore where unsolved callings occurs, in the case that
the higher elevator has passengers and in going down, and lower elevator is in
going up. Note that the actiual time to arrive the target floor is different
to the estimated time since the passengers arrival is updated at every time
step. The number of new passengers on each floor at every time step $n$ is
distributed according to Poisson law of the Poisson parameter $\lambda=\mu/K$
:

$\displaystyle P_{\lambda}(n)=\frac{\lambda^{n}}{n!}e^{-\lambda}\ .$ (1)

The parameter $\mu$ is the inflow rate of passengers of whole the building.
Each elevator could carry no more than $M$ passengers simultaneously. We
assumed the elevators are not smart enough to sense the number of carrying
passengers in it and possibly stop even if they are full.

##  3 Results

First, the relation between the arrival frequency of passengers $\mu$ and the
behaviors of an isolated single elevator in a 10-floor building ($K=10$) is
examined. Figures 2(a)-2(d) show examples of the time evolution of an
elevator’s position in a 10-floor building. In general, when $\mu$ is small,
the elevator stays on the ground floor and solves only one calling in most
round trips. As $\mu$ increases, the elevator makes round trips between the
ground floor and the higher floor continuously. Fig. 2(e) shows a proportion
of the period staying on the ground floor $r$ for various $\mu$. The value of
$r$ asymptotically approaches to 0.09 as $\mu$ increases. The exponent of the
decay toward 0.09 is $\exp(-\mu/0.03)$ in small $\mu$ and switches to
$\exp(-\mu/0.06)$ at around $\mu=0.1$.

Figure 2: [(a)-(d)] Examples of behaviors of a single elevator at $\mu=0.01$
(a), 0.04 (b), 0.06 (c), and 0.2 (d). The dark blue lines represent the
position of the elevator at each timestep. As $\mu$ increases, the elevator
moves periodically. (e) The ratio of the period during the elevator on the
ground floor $r$ vs the inflow rate $\mu$ at $K=10$ and $M=20$. As $\mu$
increases, $r$ approaches to asymptotic value $r_{0}=$0.09 as displayed with a
broken line. Inset: $r-r_{0}$ decays exponentially as $\mu$ increases. The
broken and dot lines represent slopes of $\exp(-\mu/0.03)$ and
$\exp(-\mu/0.06)$, respectively. (f) The mean and standard deviation time
since a departure to the next departure from the ground floor. Black dots
represent the mean . The shaded area represents values within the standard
deviation. The mean decreases in small $\mu$ and turns to increase at about
$\mu=0.04$. The solid cyan and broken orange lines represent the estimated
values driven from Eq.(2) and an exponential distribution, respectively.

Figure 2(e) shows the mean time since a departure to the next departure from
the ground floor $\overline{\tau}$. For small $\mu$, $\overline{\tau}$ decays
as $1/\mu$ since the elevator moves to solve a single call, whose interval
follows an exponential distribution [orange broken line in Fig.2(e)]. At
around 0.05, $\overline{\tau}$ has a minimal value. Above about 0.05,
$\overline{\tau}$ rises and approaches 130, which is time between a departure
and the next departure in the case that the elevator stops at all floors.
Turning into an increase from the $1/\mu$ decay implies that the elevator
motion shifts from the intermittent motion to the periodic motion.

If the elevator moves periodically, $\overline{\tau}$ can be represented by an
equation on some assumptions as follows. The interval time of the $j$-th
roundtrip, $\tau_{j}$, which is defined as the duration between the $j$-th and
$(j+1)$-th departures, can be described as $\tau_{j}=2z_{j}+\gamma(n_{j}+1)$,
where $z_{j}$ is the highest floor to which it goes in the $j$-th round trip,
and $n_{j}$ is the number of floors it stops. Note that adding 1 in the second
term corresponds to a stop at the ground floor. The number of new passengers
who have come during the $j$-th round trip is expressed as $m_{j}$. The number
of new passengers who have come during N($\gg 1$) times round trip can be
approximated as $\sum_{j=0}^{N}m_{j}=\mu\sum_{j=0}^{N}\tau_{j}$. Thus, the
mean number of new passengers on a round trip is:

$\displaystyle\overline{m}=\mu[2\overline{z}+\gamma(\overline{n}+1)]\ ,$ (2)

where $\overline{x}=1/N\sum_{j=0}^{N}x_{j}$.

$\overline{z}$ is estimated by $E_{\overline{m}}(z)$, the expectation of the
highest floor where passengers are waiting in the situation that
$\overline{m}$ passengers are randomly allocated to $K$ floors. Calculating
this expectations belong to the ball-in-box problem (Fig.3). Since the number
of ways that $m$ passengers are allocated to $k$ floors is $k^{\overline{m}}$,
combination satisfying the condition that the highest floor is $k$ can be
represented as $k^{\overline{m}}-(k-1)^{\overline{m}}$. Thus, the expectation
of the highest floor among those stopped during a round trip is

$\displaystyle E_{\overline{m}}(z)$ $\displaystyle=$
$\displaystyle\frac{1}{K^{\overline{m}}}\sum_{k=1}^{K}k[k^{\overline{m}}-(k-1)^{\overline{m}}]$
(3) $\displaystyle=$ $\displaystyle
K-\frac{1}{K^{\overline{m}}}\sum^{K-1}_{k=1}k^{\overline{m}}\ .$

Figure 3: Schematic illustration of the ball-in-box problem corresponding to
our model when $m$ passengers are waiting in a $K$-floors building. The
passengers are labeled according to the order of arrival. Black circles
represent individual passengers. The highest floor where passengers are
waiting is represented as $z$, and the number of floors where passengers are
waiting are represented as and $n$.

$\overline{n}$ can be regarded as the expectation value of the number of
floors where more than zero passengers arrive during a round trip. Assuming
that the typical round trip time is $\tau$ and
$\overline{m}=\mu\overline{\tau}$, the expectation value $E_{\overline{m}}(n)$
is represented as:

$\displaystyle E_{\overline{m}}(n)$ $\displaystyle=$ $\displaystyle
K\left(1-e^{\overline{m}/K}\right)\ .$ (4)

Since both Eqs.(3) and (4) are the functions of $m$, Eq.(2) is the closed-form
of $\overline{m}$. As shown as a solid line in Fig. 2(h), the value of
$2\overline{z}+\gamma(\overline{n}+1)$ derived from Eq.(2) for each $\mu$ can
appropriately reproduce the simulation result above $\mu=0.05$. It indicates
that the elevator motion is in the periodic motion phase above 0.05.
Additionally, the condition that new calls do not occur during a round trip is
$n<1$ and this condition corresponds to $\mu<0.03$. From $\mu=0.03$ to 0.05,
the intermediate motion emerges.

Next, the dynamical behaviors of the two elevators are investigated. Figures
4(a-d) show examples of the time evolution of two elevators for $\mu=0.02$,
0.1, 0.2, and 0.4. As same as the case of a single elevator, the elevators
stay on the ground floor when $\mu$ is small [Fig. 4(a)]. The elevators
exhibit round trips between the ground floor and the higher floor continuously
when $\mu$ is middle [Fig. 4(b)]. When $\mu$ is large, the two elevators
arrive at the ground floor more simultaneously as $\mu$ increases [Figs.
4(c-d)]. To quantify this tendency, we measure the simultaneousness $S$
defined as follows:

$\displaystyle S=\left|\frac{1}{N}\sum^{N}_{j=1}\exp\left(2\pi
i\frac{\tau_{j}}{\overline{\omega}}\right)\right|\ ,$ (5)

where $\tau_{j}$ is time between $j$-th and ($j+1$)-th departure of either
elevator on the ground floor, $\overline{\omega}$ is mean time of a round trip
since a departure to the next arrival of both elevators, and $N$ is a total
number of departures. The dot and cross markers in Fig. 4(e) show $S$ for
parameter sets of $(K,M)=(10,20)$ and (10,10000), respectively. In both cases,
$S$ climbs in small $\mu$ and saturates at $\mu=0.3$. As shown in Fig. 4(f),
for $(K,M)=(10,20)$, the proportion of fully packed round trips is almost 0
below $\mu=0.25$ and climbs slowly above 0.25. On the other hand, it stays 0
for $(K,M)=(10,10000)$. These results indicate that the capacity of the
elevators $M$ does not relate to simultaneous arrival. We also investigate
behavior at $(K,M)=(100,10000)$ and find that $S$ rises in smaller $\mu$ than
at $K=10$ and stabilizes from $\mu=0.2$ [red plus markers in Fig. 4(e)]. It
fluctuates in the region over 0.5, implying that the initial distance between
two elevators does not get smaller because every floor at every timestep has
passengers and the elevators need to stop all floors. The fully-packed
proportion stays 0 for various $\mu$ we examine [Fig. 4(f)].

Figure 4: [(a)-(d)] Examples of behaviors of two elevators at $\mu=0.02$ (a),
0.1 (b), 0.2 (c), and 0.4 (d). The dark blue and bright orange lines represent
the position of each elevator at each timestep. As $\mu$ increases, the
elevators move as a cluster. (e) The simultaneousness $S$ vs $\mu$. The cyan
cross and blue dot markers represent $S$ at $(K,M)=(10,10000)$ and (10,20),
respectively. The capacity $M$ does not change $S$. The red plus markers
represent $S$ at $(K,M)=(100,10000)$. (f) The proportion of fully packed round
trips vs $\mu$. Each marker style is the same as those in (e). Although $S$
shows the same behavior in $(K,M)=(10,20)$ and (100,10000), only the
proportion of $(K,M)=(10,20)$ rises with $\mu$. (g) The simultaneousness $S$
vs the normalized highest floor in a roundtrip $\overline{z}/K$.
$\overline{z}$ is the mean of the highest floors where the elevators stop in a
round trip, and $K$ is the highest floor of the building. [(h)-(i)] $T$ vs
$\mu$ of the simulation and theory at $K=10$ (g) and 100 (h). Each marker
style is the same as those in (e). The broken lines represent the theoretical
values derived from Eq.(6).

To figure out whether the elevators arrive simultaneously even in the
condition that the elevators do not arrive at the highest floor of the
building, we examine the relation between $\overline{z}/K$ and $S$ for various
$\mu$ As shown in Fig. 4(g), $S$ rises in small $z$ regions. The simultaneous
arrival occurs without elevators’ going to the highest floor, which indicates
that the simultaneous arrival also emerges from the noise even without the
given limit cycle.

We define $T$ as the round-trip time since an elevator’s departure from the
ground floor to the elevator’s next arrival on the ground floor. $T$ can be
estimated by the following discussion. As the time since the last elevator has
departed is shorter, the less number of passengers newly come at a floor.
Hence, the number of floors stopped by the elevator coming later becomes
smaller, and the speed of it becomes effectively faster than the leading one.
As a result, the distance between the two elevators becomes small with time.
Since the two elevators are symmetric and the elevator’s system is absent from
volume exclusion, the anteroposterior relation of the two elevators is
frequently switched when the distance is close to 0. Therefore, it can be
assumed that about half of the callings that occurred during the last round
trip are solved by an elevator, and the other half are solved by the other
elevator. Additionally, the probability that the new passengers come on a
floor in the time since the forward elevator has left at that floor is
$1-\exp(-\gamma\mu/K)$, since the typical time since the forward elevator has
left is $\gamma$. Hence, assuming the two move as a cluster, the number of
floors that an elevator stops during a round trip can be represented as
$\frac{n}{2}+\left(z-\frac{n}{2}\right)\left(1-e^{-\frac{\gamma\mu}{K}}\right)$,
where $z$ is the highest calling floor among those that occur during the last
round trip of the cluster, and $n$ is the number of calling floors occurring
in the last round trip of the cluster. Following the same discussion of Eq.(2)
derivation, the mean number of new passengers on a round trip is described as:

$\displaystyle\overline{m}=\mu\left\\{2\overline{z}+\gamma\left[\frac{\overline{n}}{2}+1+\left(\overline{z}-\frac{\overline{n}}{2}\right)\left(1-e^{-\frac{\gamma\mu}{K}}\right)\right]\right\\}\
.$ (6)

Since $\overline{z}$ and $\overline{n}$ are estimated as function of
$\overline{m}$ as shown in Eqs.(3) and (4), Eq.(6) is also described as closed
form of $\overline{m}$. As shown in Fig. 4(h), $T$ calculated from Eq.(6) can
represent the trend of simulated $T$. Since the assumption that the distance
of two elevators is always small cannot be adequate in a small $S$ region, the
difference between theory and simulation is relatively larger below about
$\mu=0.02$

##  4 Summary and Discussion

In summary, we examined the dynamics of elevators in a down peak period in the
case that passengers arrives according to a Poisson process. We numerically
measured the simultaneous arrival of two elevators and showed that,
irrespective of the capasity of elevators, the ordering of elevators
spontaneously emergs from the Poisson noise. We also numerically measured the
round-trip time for each level of inflow rate and reproduced the round-trip
time by a self-consistent equation considering possible combinations of floors
on which the calling occurs during the round-trip time.

The simultaneous arrival of elevators reminds us of the phase synchronization
of oscillators as discussed in previous studies. The elevators travel
periodically and show the limit cycle in the region above a certain inflow
rate, even when they are driven by the random arrival of passengers. The
interaction that the smaller distance between elevators tends to be much
smaller is the same as that of buses in the bus-route model [21, 22, 23].
However, as mentioned in the Introduction, the volume exclusion effect does
not work between elevators, not like the bus-route model. Although high-
density clusters are formed in both cases, this difference of the local effect
affects the cluster speed: The mean speed of multi elevators can be faster
than that in a single elevator system, while that of buses cannot be. We can
regard that the clusters of elevators are results of in-phase synchronization
and those of buses are results of the jamming.

## Acknowledgment

We are grateful to Tetsuya Hiraiwa, Hiroshi Kori, Katsuhiro Nishinari, Claudio
Feliciani, Daichi Yanagisawa and Hisashi Murakami for helpful discussions and
for their kind interest in this work.

## References

  * [1] Takashi Nagatani.  The physics of traffic jams.  Reports on Progress in Physics, 65(9):1331–1386, sep 2002.
  * [2] Dirk Helbing.  Traffic and related self-driven many-particle systems.  Reviews of Modern Physics, 73(4):1067–1141, oct 2001.
  * [3] Debashish Chowdhury, Ludger Santen, and Andreas Schadschneider.  Statistical physics of vehicular traffic and some related systems, may 2000.
  * [4] Fred L. Hall, Brian L. Allen, and Margot A. Gunter.  Empirical analysis of freeway flow-density relationships.  Transportation Research Part A: General, 20(3):197–210, may 1986\.
  * [5] L. Neubert, L. Santen, A. Schadschneider, and M. Schreckenberg.  Single-vehicle data of highway traffic: A statistical analysis.  Physical Review E - Statistical Physics, Plasmas, Fluids, and Related Interdisciplinary Topics, 60(6):6480–6490, dec 1999.
  * [6] Toshimitsu Musha and Hideyo Higuchi.  The 1/f fluctuation of a traffic current on an expressway.  Japanese Journal of Applied Physics, 15(7):1271–1275, jul 1976\.
  * [7] Toshimitsu Musha and Hideyo Higuchi.  Traffic Current Fluctuation and the Burgers Equation.  Japanese Journal of Applied Physics, Part 1: Regular Papers and Short Notes and Review Papers, 17(5):811–816, 1978.
  * [8] B. S. Kerner and P. Konhäuser.  Cluster effect in initially homogeneous traffic flow.  Physical Review E, 48(4):R2335, oct 1993.
  * [9] Teruhisa S. Komatsu and Shin Ichi Sasa.  Kink soliton characterizing traffic congestion.  Physical Review E, 52(5):5574–5582, nov 1995.
  * [10] Douglas A. Kurtze and Daniel C. Hong.  Traffic jams, granular flow, and soliton selection.  Physical Review E, 52(1):218–221, jul 1995.
  * [11] Thorsten Pöschel and Jason A. C. Gallas.  Synchronization effects in the dynamical behavior of elevators.  Physical Review E, 50(4):2654–2659, oct 1994.
  * [12] Takashi Hikihara and Shinichi Ueshima.  Emergent synchronization in multi-elevator system and dispatching control.  IEICE Transactions on Fundamentals of Electronics, Communications and Computer Sciences, E80-A(9):1548–1553, 1997.
  * [13] Takashi Nagatani.  Dynamical behavior in the nonlinear-map model of an elevator.  Physica A: Statistical Mechanics and its Applications, 310(1-2):67–77, jul 2002.
  * [14] Takashi Nagatani.  Complex behavior of elevators in peak traffic.  Physica A: Statistical Mechanics and its Applications, 326(3-4):556–566, aug 2003.
  * [15] Takashi Nagatani.  Dynamical transitions in peak elevator traffic.  Physica A: Statistical Mechanics and its Applications, 333:441–452, feb 2004.
  * [16] Takashi Nagatani.  Complex motion in nonlinear-map model of elevators in energy-saving traffic.  Physics Letters A, 375(20):2047–2050, may 2011.
  * [17] Takashi Nagatani and Kazuhiro Tobita.  Effect of periodic inflow on elevator traffic.  Physica A: Statistical Mechanics and its Applications, 391(18):4397–4405, sep 2012.
  * [18] Takashi Nagatani.  Complex motion induced by elevator choice in peak traffic.  Physica A: Statistical Mechanics and its Applications, 436:159–169, oct 2015.
  * [19] Takashi Nagatani.  Effect of stopover on motion of two competing elevators in peak traffic.  Physica A: Statistical Mechanics and its Applications, 444:613–621, feb 2016.
  * [20] Zhijie Feng and S. Redner.  When Will an Elevator Arrive?  dec 2020.
  * [21] O. J. O’Loan, M. R. Evans, and M. E. Cates.  Jamming transition in a homogeneous one-dimensional system: The bus route model.  Physical Review E, 58(2):1404–1418, aug 1998.
  * [22] Takashi Nagatani.  Interaction between buses and passengers on a bus route.  Physica A: Statistical Mechanics and its Applications, 296(1-2):320–330, jul 2001.
  * [23] Carlos Gershenson and Luis A. Pineda.  Why does public transport not arrive on time? The pervasiveness of equal headway instability.  PLoS ONE, 4(10), 2009.


# Theoretical Advances in Current Estimation and Navigation from a Glider-
Based Acoustic Doppler Current Profiler (ADCP)

###### Abstract

We examine acoustic Doppler current profiler (ADCP) measurements from
underwater gliders to determine glider position, glider velocity, and
subsurface current. ADCPs, however, do not directly observe the quantities of
interest; instead, they measure the relative motion of the vehicle and the
water column. We examine the lineage of mathematical innovations that have
previously been applied to this problem, discovering an unstated but incorrect
assumption of independence. We reframe a recent method to form a joint
probability model of current and vehicle navigation, which allows us to
correct this assumption and extend the classic Kalman smoothing method.
Detailed simulations affirm the efficacy of our approach for computing
estimates and their uncertainty. The joint model developed here sets the stage
for future work to incorporate constraints, range measurements, and robust
statistical modeling.

The trade-offs between cost, power, and mission duration of ocean gliders
allow them to measure underwater currents around the world. We seek to improve
the methods used by seagliders to infer true current profiles and vehicle
position from limited data. Where previous methods applied smoothing
techniques ad hoc, we start with a probability model for current and
navigation from which we derive a more powerful and flexible smoothing method.
In simulation, our method performed well, able to track within 75m a vehicle
submerged and without GPS for three hours. The flexible method allows
experimenting with innovations from recent robust process model literature,
including more difficult range measurements, and potentially real-time
processing by low-powered gliders.

##  1 Introduction

In order to map subsurface currents far away from fixed and mobile
infrastructure, underwater gliders embark small, 1 MHz acoustic Doppler
current profilers (ADCP). These high-resolution profilers, when embarked on a
vehicle such as a seaglider, measure the relative velocity of the water
column. Using these measurements poses a challenge because they reference the
vehicle’s navigation frame and only sample a local slice of the water column.
Nevertheless, one can use these relative velocities together with surface GPS
measurements and a set of useful assumptions, to separate the relative
measurement into its ground-referenced components: the subsurface ocean
current profile and vessel navigation velocity. Doing so eliminates the
reliance on navigation support from acoustic long baseline beacons or surface
tenders, as in Jakuba et al. (2021).

Previous efforts to develop a method for inferring currents and navigational
data from ADCP measurements began with low-resolution sensors that imaged the
entire water column. The sensors typically descended from either a ship or
fixed mooring. Higher-resolution profilers had a limited range, and thus
Firing and Gordon (1990) developed a ”shear” method to stitch together frames:
the change in ADCP returns are averaged across overlapping depths to form the
baroclinic part of the current profile.

Visbeck (2002) developed what is now the most widely used method, a linear
least squares model commonly referred to as the “inversion” method. The
inversion method originally involved one step to measure depth-averaged
current (DAC) and another to solve the matrix inverse problem.

Methods derived from Visbeck (2002) differ according to whether smoothing
occurs through a process-model or through a finite difference regularizer;
whether the model includes other sensors such as bottom-tracking Doppler
velocity log, inertial sensors, or terrain matching; whether the model treats
currents as horizontally and temporally stable; and what measurement errors,
such as sensor misalignment and ADCP bias, the method seeks to control. The
particulars of each method often reflect whether it aims to support seagliders
or propeller-driven AUVs. Medagoda et al. (2016) contains a literature review
of both corpora, while we focus more specifically on the former.111 The more
complex methods generally aim to satisfy the needs of propeller-driven
vehicles. Such vehicles often conduct missions more concerned with precise
navigation than current field modeling. The methods also incorporate sensors
that require more power and space than available on a typical underwater
glider, such as the tactical-grade IMU of Medagoda et al. (2010). Beyond the
review in Medagoda et al. (2016), valuable AUV-focused contributions have
occurred in Arnold and Medagoda (2018) and Medagoda and Kinsey (2016), which
test out different filters, a physics-based model, and Liu et al. (2021))
which suggests a circular calibration dive prior to mission start and a high-
pass filter in order to control IMU biases.

### Contributions and roadmap.

We present an improved inversion method based on Jonker et al. (2019), which
frames the combined smoothing/ADCP inverse problem as a single joint
optimization problem.222 In the above taxonomy, their method, designed for
seagliders, uses a process model smoother, minimal additional sensors, and a
one-dimensional current field.  We start a step further back, with the
stochastic model that gives rise to the physical problem. The more theoretical
starting point allows us to compute the proper smoothing covariance between
the current and the vehicle’s over-the-ground velocity and position. We then
derive the joint optimization problem from the posterior mode of vehicle and
current states. The approach utilizes ADCP data, hydrodynamic model data, and
GPS location fixes (at the start and, optionally, at the end of the dive). We
demonstrate the success of the method in simulations without either a final
GPS point, Doppler Velocity Log (DVL) seabed measurements, or IMU
measurements.

The paper proceeds as follows. We illustrate connections to relevant prior
work in Section 2. Section 3 describes the state-space model for our method,
including the details of covariance calculations between current and vehicle
velocity. It then describes the measurements we include in our model,
assembles the final optimization model variants, and finishes with a
discussion of uncertainty quantification. Section 4 describes our simulations,
numerical experiments, and results. We conclude with a discussion and outline
future work.

##  2 Comparison with Prior Methodologies

We begin this section with an explanation of the inversion method in Visbeck
(2002), focusing on elements that this paper’s antecedents have improved. The
basic method centers on the linear equation:

$\displaystyle\bm{\mathrm{b}}=\bm{\mathsf{A}}\bm{\mathrm{x}}$

where the vector $\bm{\mathrm{b}}$ contains the set of ADCP measurements. The
state vector, $\bm{\mathrm{x}}$, represents the desired information as a
stacked vector:

$\bm{\mathrm{x}}:=\left[\begin{matrix}\bm{\mathrm{x}}_{\upsilon}\\\
\bm{\mathrm{x}}_{c}\end{matrix}\right],$

with $\bm{\mathrm{x}}_{\upsilon}$ as the vehicle velocity at different times
and $\bm{\mathrm{x}}_{c}$ as the true ocean currents at different depths. The
measurement matrix, $\bm{\mathsf{A}}$, selects the appropriate depth of
current and time of vehicle motion for each ADCP measurement of relative
velocity. Binning the state vector makes the system overdetermined, and we can
solve it using least squares:

$\bm{\mathrm{x}}=(\bm{\mathsf{A}}^{T}\bm{\mathsf{A}})^{-1}\bm{\mathsf{A}}^{T}\bm{\mathrm{b}}.$

In the absence of binning, an acceptable condition number of
$\bm{\mathsf{A}}^{T}\bm{\mathsf{A}}$ requires additional measurement or
smoothing terms.333Without binning, the system would be underdetermined.
$\bm{\mathrm{x}}$ would have two entries for each entry of $\bm{\mathrm{b}}$,
because there is both a vehicle velocity and a current to infer at each
measurement point.

Although Visbeck (2002) bins the state vector, the paper still offers a
smoothing term for currents and two additional, mutually-exclusive measurement
terms. The first additional measurement term, that paper’s expression (26),
regularizes the depth-average of current (DAC), or barotropic current, to be
zero. Thus, the regularizer implicitly changes $\bm{\mathrm{x}}_{c}$ to model
just baroclinic currents and shifts $\bm{\mathrm{x}}_{\upsilon}$ by DAC. The
paper then suggests that DAC or surface current could be calculated separately
and added in to correct these values. On the other hand, the alternative
expression (30) in Visbeck (2002) forces the integral of vehicle velocities to
be equal to the displacement of the vehicle, as measured by GPS. Such a term
obviously requires $\bm{\mathrm{x}}_{\upsilon}$ represent true, over-the-
ground velocity and $\bm{\mathrm{x}}_{c}$ represent the true current.

Todd et al. (2017) identifies that the DAC regularization, expression 26 of
Visbeck (2002), can act instead as a measurement term if one calculates DAC
before the inversion method. He calculates the DAC “measurement” from the
difference between GPS measurement and dead reckoning, where a hydrodynamic
model supplies dead reckoning. In addition, to correct the hydrodynamic model,
the work describes a method to pre-calculate ADCP misalignment, pitch error
biases, angle of attack, and sideslip. Thus, while following the method of
expression 26, Todd’s two-step method still models vehicle over-the-ground
velocity. Although the reconciled DAC expression nominally would allow a GPS
displacement measurement term like expression (30) of Visbeck (2002), the
method already expends the GPS-measured displacement in calculating DAC. It
does, however, include a GPS-derived surface velocity measurement. Finally,
both Todd et al. (2011) and Todd et al. (2017) smooth, not just the current
velocity as in Visbeck (2002), but the vehicle velocity as well.

Jonker et al. (2019) reflect that one need not compose all the hydrodynamic
model terms and GPS displacement into a single measurement of depth-averaged
current; these terms may factor more properly as individual measurements in a
combined model. Indeed, Visbeck (2002) advocates the ease with which least
squares models can accommodate additional measurements. However, as we point
out, that ease of adding terms relies on their putative independence as much
as it relies on the linear nature of the problem formulation.444Easy
incorporation of measurements faces additional drawbacks; Jonker et al.
(2019)’s error modeling cannot accommodate all the innovation in error
modeling of Todd et al. (2017) without process terms for the sideslip and
angle-of-attack errors. The paper’s second combined model smooths currents and
vehicle process using a Kalman smoother, replacing the finite difference
smoothing of Visbeck (2002), Todd et al. (2011), and Todd et al. (2017). The
smoother requires modeling both the vehicle horizontal position as well as
velocity, which simplifies the GPS measurement term. Such a process model
approach has prior exploration in Medagoda et al. (2010) and follow-on
papers.555 Medagoda and co-authors test out a variety of filters: a Sparse
Extended Information Filter in Medagoda et al. (2010), Medagoda et al.
(2011a), and Medagoda et al. (2011b), an Extended Kalman Filter in Medagoda et
al. (2015), an Extended Information Filter in Medagoda et al. (2016), and an
Unscented Kalman Filter in Arnold and Medagoda (2018).  These papers use an
inversion approach and focus on propeller-driven AUVs with an IMU in shallower
water. Their process model filter, in distinction to a smoother, calculates
the optimal value from recent, real-time observations, as opposed to the
greater power but greater expense of retaining the entire state history, as in
the smoother of Jonker et al. (2019). As a final note, Jonker et al. (2019)
mentions that vehicle smoothing ought to properly apply to vehicle through-
the-water velocity, not over-the-ground velocity. However, because such
independent smoothing would require modeling through-the-water velocity and
consequent difficulties incorporating GPS measurement, Jonker et al. (2019)
does not do so.

Overall, we contribute several important improvements to the method. We
explicitly model the problem as a maximum likelihood estimator, including
process terms for the current and the vehicle kinematics. This probability
model lets us derive the correct covariances to smooth the vehicle’s through-
the-water velocity while still modeling the vehicle’s over-the-ground
velocity. The process model approach also allows us to avoid binning the of
data as well as to compare the smoothing of different process models. Finally,
viewing the current and vehicle smoothing as a consequence of a stochastic
process illuminates several additional assumptions whose refinement could
improve the model performance.

##  3 Methods

Here, we specify the Gaussian current and navigation random processes and
sensor measurement error. Then, we combine terms in a single optimization
problem for the latent current and vehicle variables and estimate them using
maximum likelihood.

###  3.1 Process Model

A process model for current and vehicle motion smooths available measurements
by applying a prior, or regularizer. We use the vector $\bm{\mathrm{x}}$ to
describe the combined state space of the vehicle and current. The state
vector, $\bm{\mathrm{x}}$, includes the vehicle’s velocity and position on a
local Cartesian grid and the water’s velocity. Moreover, it eschews the need
to bin observations: $\bm{\mathrm{x}}$ has vehicle kinematic entries at every
time observed and current entries at every depth observed. We use a subscript
to denote when selecting some subset of the state; e.g.
$\bm{\mathrm{x}}_{\upsilon}$ represents just the time-indexed vehicle
components of the state (both position and velocity). In implementation, any
subscript manifests as matrix multiplication. A complete list of this notation
follows:

  1. 1.

$\bm{\mathrm{x}}_{\upsilon}$ Vehicle components

  2. 2.

$\bm{\mathrm{x}}_{c}$ Current components

  3. 3.

$\bm{\mathrm{x}}_{q}$ vehicle position components, one part of
$\bm{\mathrm{x}}_{\upsilon}$

  4. 4.

$\bm{\mathrm{x}}_{\dot{q}}$ vehicle velocity components, the other part of
$\bm{\mathrm{x}}_{\upsilon}$

The method does not account for compass measurements or model heading, leaving
such items to preprocessing. Consequently, we can expect an additional 17m of
error for every degree of uncorrected compass error and every kilometer
traveled. Moreover, without an inertial measurement, we do not have a need for
earth curvature terms; a local grid will not supply more than a .1% error over
the distances of a single dive between the 45th parallels. As a final note,
the assumption that preprocessing handles vehicle and geodetic orientation
gives us all measurements in cardinal directions. As such, the northerly and
easterly components of our model comprise two non-interacting problems. When
we specify a model term, we implicitly mean two terms: one for the northerly
component, and one for the easterly component.

In our initial model for $\mathrm{Pr}(\bm{\mathrm{x}})$, we use a standard
Kalman smoother, the maximum likelihood estimator for a particle whose
velocity obeys Brownian motion. This probability model gives a negative log-
likelihood term of:

$\displaystyle-\ell(\bm{\mathrm{x}}_{\upsilon})=\frac{1}{2\sigma_{\upsilon}^{2}}\|\bm{\mathsf{G}}_{\upsilon}\bm{\mathrm{x}}_{\upsilon}\|_{\bm{\mathsf{Q}}_{\upsilon}^{-1}}^{2},$
(1)

where Appendix A contains details of the update matrix
$\bm{\mathsf{G}}_{\upsilon}$ and covariance matrix
$\bm{\mathsf{Q}}_{\upsilon}$. A user-defined covariance scaling,
$\sigma_{\upsilon}^{2}$, completes the term.

The current process adds an analogous term, with three distinctions. Firstly,
current varies by depth $s$ and not time $t$. Second, the basic model only
requires the water velocity and has no need for a current “position.” Finally,
in order to approximate the time/horizontal variance in the water column, we
treat ascent as continuing the dive by reflecting depths on ascent across the
max depth. This innovation rolls variance due to horizontal, vertical, and
time dimensions into a single process. Because of the straight paths that
seagliders tend to dive, this reflection should not degrade results
significantly. The negative log likelihood contribution for the current
process mirrors the one for $\bm{\mathrm{x}}_{\upsilon}$:

$\displaystyle-\ell(\bm{\mathrm{x}}_{c})=\frac{1}{2\sigma_{c}^{2}}\|\bm{\mathsf{G}}_{c}\bm{\mathrm{x}}_{c}\|_{\bm{\mathsf{Q}}_{c}^{-1}}^{2},$
(2)

where Appendix A contains details of the update matrix $\bm{\mathsf{G}}_{c}$
and covariance matrix $\bm{\mathsf{Q}}_{c}$. A user-defined covariance
scaling, $\sigma_{c}^{2}$, completes the term.

We can compare these terms to the smoothers and regularizers employed in
binning methods. If intervals between observations remained constant, our
method would regularize the first forwards difference of velocity.666But more
rigorously, because the Kalman smoother also includes regularization of
vehicle position according to it’s covariance with velocity. Visbeck (2002)
and Todd et al. (2011) instead regularize the 2nd centered difference, which
does not have a probabilistic representation in our model.

###  3.2 Measurement Terms

In addition to the priors responsible for smoothing the current and vehicle
processes, we need measurement terms. We use ADCP measurements of the shear
velocity of the water column at one depth and the vehicle over-the-ground
velocity at the same time. A hydrodynamic model estimates vehicle through-the-
water velocity, mathematically similar to an ADCP measurement, but for which
we allow a different error variance. Finally, a pair of GPS fixes either at
the beginning of the dive or spanning it measure the change in vehicle
position.

We present the measurement error expressions below. $\bm{\mathsf{B}}^{q},$ or
$\bm{\mathsf{B}}^{\dot{q}}$ selects the vehicle kinematics of the state vector
at the appropriate measurement times and $\bm{\mathsf{B}}^{c}$ selects the
currents of the state vector at the appropriate measurement depths. Each term
has a user-specified variance parameter $\sigma^{2}$, representing the
measurement error variance; one can set these parameters based on published
sensor accuracy. Again, we only give the explicit terms for either the
easterly or northerly components; in actuality, we use two of each term.

  1. 1.

$\bm{\mathrm{z}}_{adcp}$: ADCP data measured from the seaglider with error
variance $\sigma^{2}_{adcp}$. This measurement reflects the current at a
certain depth relative to the true velocity of the glider:

$\displaystyle-\ell(\bm{\mathrm{z}}_{adcp}|\bm{\mathrm{x}})=\frac{1}{2\sigma_{adcp}^{2}}\|\bm{\mathrm{z}}_{adcp}-(\bm{\mathsf{B}}^{\dot{q}}_{adcp}\bm{\mathrm{x}}_{\dot{q}}-\bm{\mathsf{B}}^{c}_{adcp}\bm{\mathrm{x}}_{c})\|^{2}.$
(3)

  2. 2.

$\bm{\mathrm{z}}_{ttw}$: Through-the-water estimated velocity of the seaglider
with error variance $\sigma^{2}_{ttw}$. While this quantity nominally relies
on calculations from pitch and depth rate measurements using a vehicle’s
hydrodynamic model, we treat it as direct measurement. It behaves conceptually
the same as ADCP data measured at the exact depth of the vehicle.

$\displaystyle-\ell(\bm{\mathrm{z}}_{ttw}|\bm{\mathrm{x}})=\frac{1}{2\sigma_{ttw}^{2}}\|\bm{\mathrm{z}}_{ttw}-(\bm{\mathsf{B}}^{\dot{q}}_{ttw}\bm{\mathrm{x}}_{\dot{q}}-\bm{\mathsf{B}}^{c}_{ttw}\bm{\mathrm{x}}_{c})\|^{2}.$
(4)

  3. 3.

$\bm{\mathrm{z}}_{gps}$: GPS position measurements at the surface before and
after a dive with error variance $\sigma^{2}_{gps}$.

$-\ell(\bm{\mathrm{z}}_{gps}|\bm{\mathrm{x}})=\frac{1}{2\sigma_{gps}^{2}}\|\bm{\mathrm{z}}_{gps}-(\bm{\mathsf{B}}^{q}_{gps}\bm{\mathrm{x}}_{q})\|^{2}.$

These terms compose the complete likelihood specification for the model:

$\displaystyle-\ell(\bm{\mathrm{x}}|\bm{\mathrm{z}})=-\ell(\bm{\mathrm{x}}_{\upsilon})-\ell(\bm{\mathrm{x}}_{c})-\ell(\bm{\mathrm{z}}_{p}|\bm{\mathrm{x}})-\ell(\bm{\mathrm{z}}_{ttw}|\bm{\mathrm{x}})-\ell(\bm{\mathrm{z}}_{gps}|\bm{\mathrm{x}}).$
(5)

While other authors model more robust measurement terms, we omit terms for
compass error, ADCP bias due to misalignment, or hydrodynamic model bias due
to sideslip and angle-of-attack errors aggravated by biofouling.

###  3.3 Extensions

Here we discuss the extensions and variations of our basic model.

####  3.3.1 Higher-order Smoothing.

In Subsection 3.3.1, we described vehicle over-the-ground velocity and current
as Brownian processes. The vehicle’s over the ground position then becomes the
integral of Brownian motion. Whereas previous authors have described the
concomitant regularizer as a smoothing term, strictly speaking, the assumption
of Brownian motion implies non-smoothness. We can instead expand our state
space model to include higher order terms: vehicle acceleration and the change
in current with respect to depth. These terms then become Brownian. Their
integrals, vehicle velocity and current velocity, become smooth.

This decision changes the process update matrices,
$\bm{\mathsf{G}}_{\upsilon}$ and $\bm{\mathsf{G}}_{c}$, and the covariance
matrices, $\bm{\mathsf{Q}}_{\upsilon}$ and $\bm{\mathsf{Q}}_{c}$, as described
in Appendix B.

Now that our current process has more than one order, we use
$\bm{\mathrm{x}}_{\dot{p}}$ to refer specifically to the velocity components
of the current process. $\bm{\mathrm{x}}_{c}$ then refers more generally to
all the current process terms in the same way that
$\bm{\mathrm{x}}_{\upsilon}$ refers to $\bm{\mathrm{x}}_{q}$,
$\bm{\mathrm{x}}_{\dot{q}}$, and now $\bm{\mathrm{x}}_{\ddot{q}}$

####  3.3.2 Independence.

Unfortunately, our likelihood expression 5 implies an incorrect assumption.
Adding together $-\ell(\bm{\mathrm{x}}_{\upsilon})$ and
$-\ell(\bm{\mathrm{x}}_{c})$ implies that the joint probability
$\mathrm{Pr}(\bm{\mathrm{x}}_{\upsilon},\bm{\mathrm{x}}_{c})=\mathrm{Pr}(\bm{\mathrm{x}}_{\upsilon})\mathrm{Pr}(\bm{\mathrm{x}}_{c})$,
which in turn implies independence of the vehicle’s over-the-ground velocity
from current. Previous authors who modeled over-the-ground vehicle kinematics
and additively smoothed both current and vehicle kinematics made the same
assumption de facto. Instead, the glider’s through-the-water velocity is
independent of current; its over-the-ground velocity obviously depends upon
current, and this dependence becomes more significant in larger currents:

$\displaystyle\bm{\mathrm{x}}_{\dot{q}}=\bm{\mathrm{x}}_{\dot{p}}+\bm{\mathrm{x}}_{\dot{r}},$
(6)

where, $\bm{\mathrm{x}}_{\dot{r}}$ is the relative velocity of the vehicle
through the water.

To correct our likelihood formula, we have two options. (1) We can switch from
modeling $\bm{\mathrm{x}}_{q}$ and $\bm{\mathrm{x}}_{\dot{q}}$ to modeling
$\bm{\mathrm{x}}_{r}$ and $\bm{\mathrm{x}}_{\dot{r}}$, the vehicle kinematics
relative to the water column. (2) Alternatively, we can calculate the correct
covariance between vehicle over-the-ground velocity and current in order to
use the conditional formula,
$\mathrm{Pr}(\bm{\mathrm{x}}_{\upsilon},\bm{\mathrm{x}}_{c})=\mathrm{Pr}(\bm{\mathrm{x}}_{\upsilon}|\bm{\mathrm{x}}_{c})\mathrm{Pr}(\bm{\mathrm{x}}_{c})$.
Both involve similar mathematical computations.

##### Modeling $\bm{\mathrm{x}}_{r}$

In the former variant, we model the vehicle’s kinematics relative to the water
column. However, if $\bm{\mathrm{x}}_{\dot{r}}$ represents the vehicle’s
through-the-water velocity, then $\bm{\mathrm{x}}_{r}$ represents the
vehicle’s position relative to the water. In order to incorporate GPS
measurements, we must calculate the true, over the ground position of the
vehicle by integrating equation 6. To effect that calculation, we initially
let the current velocity field, $x_{\dot{p}}$ vary as a function of depth and
time. We integrate equation 6:

$\displaystyle x_{q}(T)-x_{q}(0)$
$\displaystyle=\int_{0}^{T}x_{\dot{q}}(t)dt=\int_{0}^{T}x_{\dot{r}}(t)+x_{\dot{p}}(s(t),t)dt.$
(7)

We change variables using $\dot{s}=ds/dt$ in the second integral term:

$\displaystyle=x_{r}(T)-x_{r}(0)+\int_{s(0)}^{s(T)}x_{\dot{p}}(s,t(s))\frac{1}{\dot{s}}ds$
(8) $\displaystyle=x_{r}(T)-x_{r}(0)+x_{p}(s(T))-x_{p}(s(0)).$ (9)

Here, we recognize $x_{p}(s)$ as the position of a test particle drifting in
the current but changing depths to experience the same currents as the
vehicle. The vehicle’s over-the-ground position change would then equal the
sum of the test particle’s position and the vehicle’s through-the-water
position. Our state vector $\bm{\mathrm{x}}_{c}$ now comprises both
$\bm{\mathrm{x}}_{\dot{p}}$ and $\bm{\mathrm{x}}_{p}$.

Practically, the additional current process order manifests as a different
$\bm{\mathsf{G}}_{c}$ and $\bm{\mathsf{Q}}_{c}$ as detailed in Appendix C. In
addition to those smoothing terms, the $\bm{\mathsf{B}}_{adcp}$ and
$\bm{\mathsf{B}}_{ttw}$ matrices adjust to select the proper velocities for
their respective measurements. Moreover, the GPS measurement now requires both
a $\bm{\mathsf{B}}_{gps}^{p}$ and a $\bm{\mathsf{B}}_{gps}^{r}$ to add up the
true position change of the vehicle.

##### Covariance between $\mathrm{Pr}(\bm{\mathrm{x}}_{\upsilon})$ and
$\mathrm{Pr}(\bm{\mathrm{x}}_{c})$

If instead we model vehicle kinematics over the ground (as in the basic
method), we must derive our maximum likelihood estimate from:

$\displaystyle\mathrm{Pr}(\bm{\mathrm{x}}_{q},\bm{\mathrm{x}}_{\dot{q}},\bm{\mathrm{x}}_{c})=\mathrm{Pr}(\bm{\mathrm{x}}_{q},\bm{\mathrm{x}}_{\dot{q}}|\bm{\mathrm{x}}_{c})\mathrm{Pr}(\bm{\mathrm{x}}_{c})$
(10)
$\displaystyle\ell(\bm{\mathrm{x}}_{q},\bm{\mathrm{x}}_{\dot{q}},\bm{\mathrm{x}}_{c})=\ell(\bm{\mathrm{x}}_{q},\bm{\mathrm{x}}_{\dot{q}}|\bm{\mathrm{x}}_{c})+\ell(\bm{\mathrm{x}}_{c}).$
(11)

We achieve the conditional likelihood expression by reusing equation 8 to get
the joint probability, then conditioning on $\bm{\mathrm{x}}_{c}$. Appendix D
provides more details.

###  3.4 Optimization and Uncertainty Quantification

The likelihood approach of (5) corresponds to a Bayesian interpretation of the
loss. Our prior belief of the Brownian process terms combined with the
Gaussian measurement likelihoods results in a combined posterior and its
estimator that are Gaussian random variables. Given (5), one can represent the
posterior estimator of $\bm{\mathrm{x}}$ as:

$\displaystyle\hat{\bm{\mathrm{x}}}=\arg\min_{\bm{\mathrm{x}}}$
$\displaystyle\frac{1}{2\sigma^{2}_{\upsilon}}\|\bm{\mathsf{G}}_{\upsilon}\bm{\mathrm{x}}_{\upsilon}\|^{2}_{\bm{\mathsf{Q}}_{\upsilon}^{-1}}+\frac{1}{2\sigma^{2}_{c}}\|\bm{\mathsf{G}}_{c}\bm{\mathrm{x}}_{c}\|^{2}_{{\bm{\mathsf{Q}}_{c}}^{-1}}+$
$\displaystyle\frac{1}{2\sigma^{2}_{adcp}}\left\|\bm{\mathrm{z}}_{adcp}-(\bm{\mathsf{B}}^{\dot{q}}_{adcp}\bm{\mathrm{x}}_{\dot{q}}-\bm{\mathsf{B}}^{\dot{p}}_{adcp}\bm{\mathrm{x}}_{\dot{p}})\right\|^{2}+$
(12)
$\displaystyle\frac{1}{2\sigma^{2}_{ttw}}\left\|\bm{\mathrm{z}}_{ttw}-(\bm{\mathsf{B}}^{\dot{q}}_{ttw}\bm{\mathrm{x}}_{\dot{q}}-\bm{\mathsf{B}}^{\dot{p}}_{ttw}\bm{\mathrm{x}}_{\dot{p}})\right\|^{2}+$
$\displaystyle\frac{1}{2\sigma^{2}_{gps}}\left\|\bm{\mathrm{z}}_{gps}-(\bm{\mathsf{B}}^{q}_{gps}\bm{\mathrm{x}}_{q}-\bm{\mathsf{B}}^{p}_{gps}\bm{\mathrm{x}}_{p})\right\|^{2}.$

One can compute the estimator in closed form, using the normal equations.
Knowing that any subscript on $\bm{\mathrm{x}}$ practically represents matrix
multiplication and that both $\bm{\mathsf{Q}}^{-1}$, as positive definite
matrices, factor to $LL^{T}$, we can rewrite every term as

$\displaystyle\frac{1}{2\sigma^{2}_{i}}\left\|\bm{\mathsf{A}}_{i}x-\bm{\mathrm{b}}_{i}\right\|^{2}.$
(13)

And stacking all $\bm{\mathsf{A}}_{i}/\sigma_{i}$ into matrix
$\bm{\mathsf{A}}$ and all $\bm{\mathrm{b}}_{i}/\sigma_{i}$ into vector
$\bm{\mathrm{b}}$, we can write the optimal value as:

$\displaystyle\hat{\bm{\mathrm{x}}}=\bm{\mathsf{A}}\backslash\bm{\mathrm{b}}=\bm{\mathsf{A}}^{\dagger}\bm{\mathrm{b}}=(\bm{\mathsf{A}}^{T}\bm{\mathsf{A}})^{-1}\bm{\mathsf{A}}^{T}\bm{\mathrm{b}}.$
(14)

Moreover, the variance of the estimator $\bm{\mathrm{\hat{x}}}$ has an
explicit solution:

$\displaystyle\textrm{var}(\hat{\bm{\mathrm{x}}})=\mathbb{E}[(\hat{\bm{\mathrm{x}}}-\bm{\mathrm{x}})(\hat{\bm{\mathrm{x}}}-\bm{\mathrm{x}})^{T}]=(\bm{\mathsf{A}}^{T}\bm{\mathsf{A}})^{-1}.$
(15)

Here, $\mathbb{E}$ is the expectation operator.

##  4 Simulation Results

###  4.1 Simulation setup

We conducted simulations to compare our method and its variants. Each trial
includes a variant with GPS points at both ends of the dive, and an identical
variant trial with two GPS points at the beginning and none at the end. This
second variant behaves similarly to how Todd et al. (2017) includes initial
surface drift velocity, as measured by GPS, as a term in his least-squares
solution. However, we leave the observations as positional in nature which
allows the Kalman machinery to manage the covariance between positions and
velocities. As discussed later in the conclusion, this ability to swap in
different measurements shows the flexibility of the probability model
framework. Table 1 lists simulation parameters.

Parameter | Value
---|---
Duration | 3 hrs
Depth | 750 m
Hydrodynamic model measurements | 500
ADCP measurements | 450
ADCP bins | 4
ADCP range | 12m, upwards-facing
Current | piecewise sinusoid with O(.6) kts
Vehicle TTW velocity | piecewise sinusoid with O(.6) kts
Hydrodynamic model noise | O(1e-2) m/s
ADCP measurement noise | O(1e-2) m/s
GPS measurement noise | O(1) m
Table 1: Simulation Parameters, reflecting likely dive of a Seaglider with
Nortek ADCP

###  4.2 Parameters and grid-search

While published information about physical devices can inform setting the
measurement variance parameters $\sigma^{2}$, the process variance parameters
for smoothing represent prior knowledge of system (vehicle and current)
dynamics. We wish to understand how precisely our model’s success depends upon
knowing the exact process model covariance, a parameter we can only ever
assume imprecisely. Thus, in each simulation, we run a parameter search on the
covariance coefficients for the process terms, $\sigma_{\upsilon}^{2}$ and
$\sigma_{c}^{2}$, and evaluate the effect on accuracy metrics.

###  4.3 Results

We compare four methods based on different process models: (1) the basic one
described in section 3.3.1, (2) the higher-order one described in section
3.3.3.3.3.1, (3) the process with proper vehicle-current covariance described
in section 3.3.3.3.3.2, and (4) a combination of methods 3 and 4: a proper
covariance method that smooths to higher order. While section 3.3.3.3.3.2
argues that modeling through-the-water kinematics is mathematically equivalent
to modeling true kinematics with appropriate vehicle-current covariance on
paper, the higher condition numbers made the former less successful. Thus, we
only display results for the latter method. In comparing methods, we use as
criteria: (a) The horizontal position accuracy during the dive at optimal
process parameters777 When the process parameters that give rise to the most
accurate navigation estimates and those that give rise to the most accurate
current estimates disagree, we choose as optimal the parameters that lead to
the most accurate current estimation. , (b) The current accuracy during the
dive at optimal process parameters, (c) The breadth of process parameters that
achieve satisfactory results, and (d) performance on simulations including
only two closely-placed GPS fixes before the dive and no GPS fix at the end.
We evaluate root mean square error (RMSE) across the entire history, averaged
over twenty independent trials and displayed in Table 2. Fig. 1 shows the
panoply of visualizations produced by each trial, and Appendix E contains the
visualizations for each trial, organized for side-to-side comparisons of
different trials.

Method | Variant | Navigation Error | Current Error | Nav err (no final GPS) | Curr err (no final GPS)
---|---|---|---|---|---
1 | Basic | 102 | 4.13e-2 | 380 | 8.40e-2
2 | Higher-order | 107 | 3.83e-2 | 219 | 6.09e-2
3 | Vehicle-Current covariance | 140 | 5.29e-2 | 377 | 8.42e-2
4 | Covariance and Higher-order | 99.1 | 3.71e-2 | 249 | 6.72e-2
Table 2: Results, with the optimal method for each column boldfaced. Navigation RMSE uses units of meters, current RMSE uses units of meters per second. All methods perform competitively, with a different method taking top spot for each metric. |
---|---
(a) | (b)
|
(c) | (d)
|
(e) | (f)
Figure 1: Visualizations of evaluation criteria for solution method 1 with
both starting and ending GPS measurements. (a) Navigation solution, (b)
Current solution, (c) Navigation error across process covariance parameter
search, with minimizer marked by red dot (d) Current error across process
covariance parameter search, with minimizer marked by red dot(e) Vehicle and
current velocity estimator uncertainty, and (f) Vehicle position estimator
uncertainty. The parameter search plots average twenty simulations. The rest
of the plots represent a single simulation from among those twenty. See
appendix E for all plots from all trials, organized for comparison across
trials.

All solution methods appear to perform adequately, and both higher order
methods appear the best at two particular metrics. Considering first
navigation error, most methods track the vehicle with at most a 75m difference
between true and modeled position. Even with just starting GPS positions, the
fourth method keeps track of the vehicle with an accuracy of 100m after four
hours. For this particularly challenging track, the level of accuracy
displayed in Fig. 2 reflects an impressive success. Looking in more detail at
how the choice of vehicle and current process parameters affect navigation
error in Fig. 3, we see that the higher order methods (2 and 4) allow a much
wider region of parameter space to produce satisfactory results.

|
---|---
(a) | (b)
Figure 2: Navigation error (a) method 1, with start and end GPS points, and (b) method 4, including only starting GPS points. All of this paper’s methods perform well, outperforming DR corrected by DAC and tracking the vehicle to O(25m/hr) before receiving a resurfacing GPS fix. |
---|---
(a) | (b)
Figure 3: Navigation error across process covariance parameter search (a) in
method 1 and (b) method 2. Both cases shown include only starting GPS points,
and the higher-order smoothing methods have a much larger acceptable region of
parameter space.

As for current error, the higher order methods show a much less noisy
solution; Figure 4 illustrates the comparison between methods 3 and 4 for
cases including only starting GPS points. The higher-order methods, methods 2
and 4, in the variant without a final GPS point both outperform the other
methods and appear visually nicer. Likewise, the much wider acceptable
parameter region for higher order methods makes them better candidates for a
production model. We omit a separate plot for current error across the process
covariance parameter search, as it looks similar to Fig. 3.

|
---|---
(a) | (b)
Figure 4: Current error (a) in method 3 and (b) method 4. Both solutions shown
include only starting GPS points. Both the results on prepared metrics and the
much smoother profile of the higher order methods recommend methods 2 and 4.

The estimator variance plots tell several interesting stories. Simulations
including only beginning GPS points, such as in Fig. 5, show position and
velocity uncertainty growing with duration since the initial GPS
observations—as expected. The square-corner covariance plots also appear to
exactly replicate those of Brownian motion, providing more confidence that our
model performs consistent with assumptions.

Figure 5: Method 1 velocity error covariance without a final GPS point appears
appropriately similar to the covariance of Brownian motion. Additionally, the
repeated off-diagonal terms indicate a lack of identifiability. All methods
including only starting GPS points demonstrate these effects.

Simulations that included both start and finish GPS measurements, as in Fig.
6, predictably show low position uncertainty at the endpoints and high
uncertainty at nadir. On the other hand, velocity uncertainty reaches maxima
at the endpoints and a minimum at nadir. Given that GPS measurements affect
the integral of velocity, such terms constrain velocities most at the
midpoint, where the velocity has the most adjacent points for smoothing. This
effect shows that additional measurements reduce uncertainty of the observed
quantity, but how that information passes to different process terms depends
upon the smoother.

|
---|---
(a) | (b)
Figure 6: Method 1 (a) velocity and (b) position estimator error covariance
including both starting and ending GPS measurements. This example shows that
position uncertainty reaches a minimum in vicinity of GPS observations, while
velocity uncertainty reaches a maximum.

In all cases, the off-diagonal blocks show that vehicle velocity terms
correlate highly with their simultaneous currents. Such structure indicates a
rank deficiency in our model, a lack of identifiability in whether an
individual velocity reflects vehicle motion or currents. It suggests that
direct observation of a just a few true vehicle velocities or currents during
the dive would dramatically improve estimates. Indeed, Medagoda et al. (2016)
finds that such a dramatic correction in navigation estimates occurs as soon
as the vehicle gets a DVL lock on the seafloor.

As a final note, while some of our second-order results may appear less smooth
than those of other authors’, much of that roughness merely reflects the
decision not to bin observations.

##  5 Conclusion and Further Research

We generated data and showed that our basic model, essentially a method of
Jonker et al. (2019) without binning, performed reasonably well. Extending the
model to higher order to enforce smoothness of velocities improved the model
substantially. Our other innovation, the proper covariance relationship
between the vehicle process and the current process, yielded only the moral
gain of a more rigorous model. All variants significantly outperformed dead
reckoning. In particular, models succeeded without the final GPS point as long
as two nearby GPS observations before diving helped identify true surface
drift.

Various ancillary corrections may help improve the evidence for our major
contributions. Most significantly, we may prefer to return to binning
observations. Matrix conditioning caused many of the difficulties with higher
order models and even with the proper-covariance models. Our time and distance
scaling might also prevent the proper-covariance models from showing their
merit. Additionally, measurement terms could benefit from reintroducing ADCP
bias, compass error, and biofouling errors. Moreover, we have yet to test our
model on real data.

As a key avenue for improvement, one could replace the direct matrix inversion
of our method with a generic convex solver (which produces equivalent results
in the least squares case). The new solver, in turn, would allow for robust
vehicle and current processes other than Brownian motion, a topic of interest
in Kalman smoothing literature such as Aravkin et al. (2017). Other processes
such as an Ornstein-Uhlenbeck process could also help by introducing the
assumption that current centers around zero and the vehicle velocity centers
around its design drift velocity.

The convex solver would also permit more robust measurement forms than the
linear least-squares case. Thus, we could include nonlinear measurements such
as range in our framework, as from Ultra-Short Baseline (USBL) or Long
Baseline (LBL) systems or even the range-azimuth-elevation of the USBL method
in Jakuba et al. (2021). For measurements close to the $c_{x}$ and $c_{y}$,
the local coordinates of the sensor, range measurement adds a log likelihood
function of:

$-\ell(\bm{\mathrm{z}}_{range}|x)=\frac{1}{2\sigma_{range}^{2}}\|\bm{\mathrm{z}}_{range}^{2}-(\bm{\mathsf{B}}^{q}_{range}\bm{\mathrm{x}}_{q}-\bm{\mathrm{c}}_{x})^{2}-(\bm{\mathsf{B}}^{q}_{range}\bm{\mathrm{x}}_{q}-\bm{\mathrm{c}}_{y})^{2}\|^{2},$

where the vector exponents apply elementwise. These measurements, as
nonlinear, give rise to a nonconvex likelihood expression. Even so, the
expression may still permit finding a good enough local solution, or in cases
a global one.

Finally, automatic tuning of process covariance parameters could benefit
usability of the model. While we obtained stable results in simulation using
grid search for these parameters, a comprehensive approach for simultaneous
state and current estimation with more efficient ways of finding underlying
variance parameters would more completely solve the engineering problem.

The ability to calculate navigation and current state in real time without
external submerged position support allows for more effective use of
seagliders. For example, the vehicle could choose to navigate at a depth with
beneficial currents to extend its range. It could also dispose of the raw
measurements and store the current profile in a compressed form, allowing
longer dives and longer deployments. Such innovations would vastly improve the
ability to understand global subsurface currents.

[A] Kalman Matrices Referring to time index $j$ for a coordinate $x_{q}^{j}$
with velocity $x_{\dot{q}}^{j}$, the Kalman smoother assumes linear dynamics
such that:

$\displaystyle\mathbb{E}\left[x^{j}_{\dot{q}}-x_{\dot{q}}^{j-1}|x_{\dot{q}}^{j-1}\right]$
$\displaystyle=0$
$\displaystyle\mathbb{E}\left[x^{j}_{q}-x_{q}^{j-1}|x_{q}^{j-1},x_{\dot{q}}^{j-1}\right]$
$\displaystyle=\Delta t\cdot x_{\dot{q}}^{j-1}.$

As a matrix equation, this expectation takes the form

$\displaystyle\mathbb{E}[\bm{\mathsf{G}}_{\upsilon}\bm{\mathrm{x}}_{\upsilon}]=0,$

where

$\displaystyle\bm{\mathsf{G}}_{\upsilon}=\left[\begin{matrix}-1&0&1&0&&\\\
-\Delta t_{1}&-1&0&1&&\\\ &&-1&0&1&0\\\ &&-\Delta t_{2}&-1&0&1\\\
&&&&&&\ddots\end{matrix}\right]~{}~{}~{}\text{and}~{}~{}~{}~{}\bm{\mathrm{x}}_{\upsilon}=\left[\begin{matrix}x_{\dot{q}}^{1}\\\
x_{q}^{1}\\\ \vdots\\\ x_{\dot{q}}^{n}\\\ x_{q}^{n}\end{matrix}\right].$ (16)

The Kalman smoother also assumes that observations are normally distributed
and specifies the covariance for each timestep, i.e.
$\bm{\mathrm{x}}^{j}_{\upsilon}|\bm{\mathrm{x}}^{j-1}_{\upsilon}\sim\mathcal{N}((\bm{\mathsf{G}}_{\upsilon}\bm{\mathrm{x}})^{j},\bm{\mathsf{Q}}_{\upsilon}^{j})$,
with successive timesteps independent. The overall covariance matrix
$\bm{\mathsf{Q}}_{\upsilon}$ has block diagonals
$\bm{\mathsf{Q}}_{\upsilon}^{j}$ equal to

$\displaystyle\left[\begin{matrix}\Delta t&\Delta t^{2}/2\\\ \Delta
t^{2}/2&\Delta t^{3}/3\end{matrix}\right].$ (17)

This Kalman update and covariance matrix equate to the covariance matrix of a
Brownian $\bm{\mathrm{x}}_{\dot{q}}$ when multiplied through the likelihood
expression 1.

In $\bm{\mathrm{x}}_{c}$, where the state vector only includes a velocity, one
can simply drop the rows/columns of $\bm{\mathsf{G}}$ and $\bm{\mathsf{Q}}$
corresponding to position. Finally, since current varies by depth $s$, we
replace $\Delta t$ with $\Delta s$.

[B] Higher Order Smoothing In cases when we include an acceleration term, we
can simply integrate the Brownian motion one more time to get
$\bm{\mathsf{G}}_{\upsilon}$ and $\bm{\mathsf{Q}}_{\upsilon}$:

$\displaystyle\bm{\mathsf{G}}_{\upsilon}=\left[\begin{matrix}-1&0&0&1&0&0&&\\\
-\Delta t_{1}&-1&0&0&1&0&&&\\\ -\frac{\Delta{t_{1}}^{2}}{2}&-\Delta
t_{1}&-1&0&0&1&\\\ &&&-1&0&0&1&0&0\\\ &&&-\Delta t_{2}&-1&0&0&1&0\\\
&&&-\frac{\Delta{t_{2}}^{2}}{2}&-\Delta t_{2}&-1&0&0&1\\\
&&&&&&&&&\ddots\end{matrix}\right],\qquad\bm{\mathsf{Q}}_{\upsilon}=\left[\begin{matrix}{\Delta
t}&\frac{\Delta t^{2}}{2}&\frac{\Delta t^{3}}{6}\\\ \frac{\Delta
t^{2}}{2}&\frac{\Delta t^{3}}{3}&\frac{\Delta t^{4}}{8}\\\ \frac{\Delta
t^{3}}{6}&\frac{\Delta t^{4}}{8}&\frac{\Delta t^{5}}{20}\\\
\end{matrix}\right].$ (18)

The current process now has two total orders, and so has similar terms to
vehicle smoothing in Appendix A, but again replacing $t$ with $s$.

[C] Modeling Vehicle Relative Velocity Return to equation 8, where we sought
to understand the implication of smoothing $x_{r}$, rather than $x_{q}$, and
recognized $x_{p}$ terms. If we focus on the second integral,

$\int_{s(t_{1})}^{s(t_{2})}x_{\dot{p}}(s,t(s))\frac{1}{\dot{s}}ds,$

we ought to make several assumptions explicit:

  1. 1.

Current does not change in the time or horizontal strata between observations

  2. 2.

Depth rate remains constant between two nearby observations.

The first assumptions simplifies our $x_{\dot{p}}$ term, and the second allows
the $1/\dot{s}$ becomes a coefficient. Thus, our current process term,
equation 2, now uses the process update matrix:

$\displaystyle\bm{\mathsf{G}}_{c}=\left[\begin{matrix}-1&0&1&0&&\\\ -\Delta
s_{1}/\dot{s}_{1}&-1&0&1&0&\\\ 0&0&-1&0&1&\ddots\\\ 0&0&-\Delta
s_{2}/\dot{s}_{2}&-1&\ddots&\ddots\end{matrix}\right],~{}~{}~{}~{}~{}~{}~{}\bm{\mathrm{x}}_{c}=\left[\begin{matrix}x_{\dot{p}}^{1}\\\
x_{p}^{1}\\\ \vdots\\\ x_{\dot{p}}^{n}\\\ x_{p}^{n}\end{matrix}\right].$ (19)

Likewise, $\bm{\mathsf{Q}}_{c}$ now has blocks:

$\displaystyle\left[\begin{matrix}\Delta s&\Delta s^{2}/(2\dot{s})\\\ \Delta
s^{2}/(2\dot{s})&\Delta
s^{3}/(3\dot{s}^{2})\end{matrix}\right].~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$
(21)

[D] Proper Covariance for Current and Vehicle True Velocity We seek a
specification for
$\mathrm{Pr}(\bm{\mathrm{x}}_{\dot{q}},~{}\bm{\mathrm{x}}_{q}|~{}\bm{\mathrm{x}}_{c})$
in equation 10. We first form the correct covariance matrix for
$\mathrm{Pr}(\bm{\mathrm{x}}_{\dot{q}},~{}\bm{\mathrm{x}}_{q},~{}\bm{\mathrm{x}}_{c})$
and then derive the conditional covariance matrix and expected value. Since
$\bm{\mathsf{G}}_{\upsilon}$ uses the expected value to recenter the process,
we need only concern ourselves with a single timestep $\Delta t$.

Although we do not utilize an $\bm{\mathrm{x}}_{p}$ term, we still return to
equation 6. Writing out all the terms as random variables driven by
independent Brownian motions $W_{1}^{s}$ driving current and $W_{2}^{t}$
driving the vehicle’s speed through the water, we calculate:

$\displaystyle x_{c}^{j}-x_{c}^{j-1}$
$\displaystyle=\int_{s(t_{j-1})}^{s(t_{j})}dW_{1}^{s}$
$\displaystyle\sim\int_{0}^{\Delta s}dW_{1}^{s}$ $\displaystyle
x_{\dot{q}}^{j}-x_{\dot{q}}^{j-1}$
$\displaystyle=\int_{t_{j-1}}^{t_{j}}dW_{2}^{t}+\int_{s(t_{j-1})}^{s(t_{j})}dW_{1}^{s}$
$\displaystyle\sim\int_{0}^{\Delta t}dW_{2}^{t}+\int_{0}^{\Delta s}dW_{1}^{s}$
$\displaystyle x_{q}^{j}-x_{q}^{j-1}$
$\displaystyle=\int_{t_{j-1}}^{t_{j}}dt~{}x_{\dot{q}}(t)$ $\displaystyle
x_{q}^{j}-x_{q}^{j-1}-\Delta t\cdot x_{\dot{q}}^{j-1}$
$\displaystyle=\int_{t_{j-1}}^{t_{j}}dt~{}x_{\dot{q}}(t)-x_{\dot{q}}^{j-1}$
$\displaystyle=\int_{t_{j-1}}^{t_{j}}dt\int_{t_{j-1}}^{t}dW_{2}^{r}+\int_{s(t_{j-1})}^{s(t)}dW_{1}^{r}$
$\displaystyle\sim\int_{0}^{\Delta t}dt\int_{0}^{t}dW_{2}^{r}+\int_{0}^{\Delta
t}dt\int_{0}^{s(t)}dW_{1}^{r}$ $\displaystyle=\int_{0}^{\Delta
t}dt\int_{0}^{t}dW_{2}^{r}+\frac{1}{\dot{s}}\int_{0}^{\Delta
s}ds\int_{0}^{s}dW_{1}^{r}.$

We have used the independence of successive increments and stationarity to
change from $[t_{j-1},t_{j}]$ to simpler variables $[0,\Delta t]$. For further
simplification, we now let $\Delta x_{c}=x_{c}^{j}-x_{c}^{j-1}$, $\Delta
x_{\dot{q}}=x_{\dot{q}}^{j}-x_{\dot{q}}^{j-1}$, and $\Delta
x_{q}=x_{q}^{j}-x_{q}^{j-1}-\Delta t\cdot x_{\dot{q}}^{j-1}$. These mean-zero,
Gaussian random variables have the covariance matrix:

$\displaystyle Q=$ $\displaystyle\left[\begin{matrix}\Delta s&\Delta
s&\frac{\Delta s^{2}}{2\dot{s}}\\\ \Delta s&\Delta t+\Delta s&\frac{\Delta
t^{2}}{2}+\frac{\Delta s^{2}}{2\dot{s}}\\\ \frac{\Delta
s^{2}}{2\dot{s}}&\frac{\Delta t^{2}}{2}+\frac{\Delta
s^{2}}{2\dot{s}}&\frac{\Delta t^{3}}{3}+\frac{\Delta s^{3}}{3\dot{s}^{2}}\\\
\end{matrix}\right].~{}~{}{\color[rgb]{.5,.5,.5}\left.\begin{matrix}\Delta
x_{c}\\\ \Delta x_{\dot{q}}\\\ \Delta x_{q}\end{matrix}\right.}$
$\displaystyle{\color[rgb]{.5,.5,.5}\begin{matrix}~{}~{}\Delta
x_{c}~{}~{}~{}~{}~{}&\Delta x_{\dot{q}}~{}~{}~{}~{}~{}~{}&\Delta
x_{q}\end{matrix}}$

Marginalizing current and then calculating the conditional distribution gives:

$\displaystyle\Delta x_{c}$ $\displaystyle\sim\mathcal{N}(0,\Delta s)$ (22)
$\displaystyle\left.\left[\begin{matrix}\Delta x_{\dot{q}}\\\ \Delta
x_{q}\end{matrix}\right]\right|\Delta x_{c}$
$\displaystyle\sim\mathcal{N}\left(\left[\begin{matrix}\Delta x_{c}\\\
\frac{\Delta s}{2\dot{s}}\Delta
x_{c}\end{matrix}\right],\left[\begin{matrix}\Delta t&\frac{\Delta
t^{2}}{2}\\\ \frac{\Delta t^{2}}{2}&\frac{\Delta t^{3}}{3}+\frac{\Delta
s^{3}}{12\dot{s}^{2}}\end{matrix}\right]\right).$ (23)

When smoothing the vehicle process to higher order and modeling acceleration
in this manner, we must also smooth the current process to the same order.
Assuming the left side of equation 6 to be smooth would preclude allowing the
right side to be non-smooth. Thus, to smooth current to higher order, we re-
notate $x_{c}(s)$ as $x_{\dot{p}}(s)$, the current velocity, and
$x_{\ddot{p}}(s)$, the change in current with depth. The $\Delta$ variables
(e.g. $\Delta x_{q}$, $\Delta x_{\dot{p}}$) then also must adjust to maintain
a mean-zero process. The conditional distribution becomes:

$\displaystyle\left.\left[\begin{matrix}\Delta x_{\ddot{q}}\\\ \Delta
x_{\dot{q}}\\\ \Delta
x_{q}\end{matrix}\right]\right|\left[\begin{matrix}\Delta x_{\ddot{p}}\\\
\Delta
x_{\dot{p}}\end{matrix}\right]\sim\mathcal{N}\left(\left[\begin{matrix}\dot{s}\Delta
x_{\ddot{p}}\\\ \Delta x_{\dot{p}}\\\ \frac{\Delta s}{2\dot{s}}\Delta
x_{\dot{p}}-\frac{\Delta s^{2}\Delta x_{\ddot{p}}}{\dot{s}12}\\\
\end{matrix}\right],\left[\begin{matrix}\Delta t&\frac{\Delta
t^{2}}{2}&\frac{\Delta t^{3}}{6}\\\ \frac{\Delta t^{2}}{2}&\frac{\Delta
t^{3}}{3}&\frac{\Delta t^{4}}{8}\\\ \frac{\Delta t^{3}}{6}&\frac{\Delta
t^{4}}{8}&\frac{\Delta t^{5}}{20}+\frac{\Delta s^{5}}{720\dot{s}^{2}}\\\
\end{matrix}\right]\right).$ (24)

In both expressions 22 and 24, we roll the conditional mean into our
definition of $\bm{\mathsf{G}}_{\upsilon}$ in order to recover a mean-zero
random variable.

[E] Complete Plots for Comparison

|
---|---
(a) | (b)
|
(c) | (d)
|
(e) | (f)
|
(g) | (h)
Figure 7: Representative navigation results from (a, b) method 1, (c, d) method 2, (e, f) method 3, and (g, h) method 4. Trials on the left include start and end GPS points, trials on the right include two starting GPS points only |
---|---
(a) | (b)
|
(c) | (d)
|
(e) | (f)
|
(g) | (h)
Figure 8: Representative current results from (a, b) method 1, (c, d) method 2, (e, f) method 3, and (g, h) method 4. Trials on the left include start and end GPS points, trials on the right include two starting GPS points only |
---|---
(a) | (b)
|
(c) | (d)
|
(e) | (f)
|
(g) | (h)
Figure 9: Navigation parameter search results from (a, b) method 1, (c, d) method 2, (e, f) method 3, and (g, h) method 4. Trials on the left include start and end GPS points, trials on the right include two starting GPS points only. All colors normalized to same values. |
---|---
(a) | (b)
|
(c) | (d)
a  |
(e) | (f)
|
(g) | (h)
Figure 10: Current parameter search results from (a, b) method 1, (c, d) method 2, (e, f) method 3, and (g, h) method 4. Trials on the left include start and end GPS points, trials on the right include two starting GPS points only. All colors normalized to same values. |
---|---
(a) | (b)
|
(c) | (d)
|
(g) | (h)
|
(i) | (j)
Figure 11: Velocity estimator covariance (subset of all states) from (a, b) method 1, (c, d) method 2, (e, f) method 3, and (g, h) method 4. Trials on the left include start and end GPS points, trials on the right include two starting GPS points only. Colors are not normalized to same values. |
---|---
(a) | (b)
|
(c) | (d)
|
(g) | (h)
|
(i) | (j)
Figure 12: Navigation position estimator covariance (subset of all states)
from (a, b) method 1, (c, d) method 2, (e, f) method 3, and (g, h) method 4.
Trials on the left include start and end GPS points, trials on the right
include two starting GPS points only. Colors are not normalized to same
values.

###### Acknowledgements.

This research was supported by the NOAA Office of Exploration and Research
award NA20OAR0110429 and by the Veterans Administration under the GI Bill. One
can generate the figures and data for this paper by downloading our Python
package at https://github.com/UW-AMO/seaglider-navigation and running
”publish_figures.py”. The commit used to generate the exact figures in the
paper is 21452bc.

## References

  * Aravkin et al. (2017) Aravkin, A., J. V. Burke, L. Ljung, A. Lozano, and G. Pillonetto, 2017: Generalized kalman smoothing: Modeling and algorithms. Automatica, 86, 63–86.
  * Arnold and Medagoda (2018) Arnold, S., and L. Medagoda, 2018: Robust model-aided inertial localization for autonomous underwater vehicles. IEEE conference on Robotics and Automation, 4889–4896.
  * Firing and Gordon (1990) Firing, E., and R. Gordon, 1990: Deep ocean acoustic doppler current profiling. Proceedings of the IEEE Fourth Working Conference on Current Measurement, 192–201, 10.1109/CURM.1990.110905.
  * Jakuba et al. (2021) Jakuba, M. V., J. W. Partan, S. E. Webster, D. Giaya, and C. Ramirez, 2021: Performance of a Low-Power One-Way Travel-Time Inverted Ultra-Short Baseline Navigation System.
  * Jonker et al. (2019) Jonker, J., A. Shcherbina, R. Krishfield, L. Van Uffelen, A. Aravkin, and S. E. Webster, 2019: Preliminary Results in Current Profile Estimation and Doppler-aided Navigation for Autonomous Underwater Gliders. 1–8, 10.1109/oceanse.2019.8867108, 1907.02897.
  * Liu et al. (2021) Liu, X., X. Liu, L. Wang, Y. Huang, and Z. Wang, 2021: SINS/DVL integrated system with current and misalignment estimation for midwater navigation. IEEE Access, 9, 51 332–51 342, 10.1109/ACCESS.2021.3069469.
  * Medagoda et al. (2010) Medagoda, L., M. V. Jakuba, O. Pizarro, and S. B. Williams, 2010: Water column current profile aided localisation for autonomous underwater vehicles. OCEANS’10 IEEE Sydney, OCEANSSYD 2010, 1–10, 10.1109/OCEANSSYD.2010.5604016.
  * Medagoda and Kinsey (2016) Medagoda, L., and J. C. Kinsey, 2016: Water-Current and IMU Aided AUV Localization in Deep Mid-Water. IEEE International Conference on Robotics and Automation,  (Dvl), 2–4.
  * Medagoda et al. (2015) Medagoda, L., J. C. Kinsey, and M. Eilders, 2015: Autonomous Underwater Vehicle localization in a spatiotemporally varying water current field. Proceedings - IEEE International Conference on Robotics and Automation, 2015-June (June), 565–572, 10.1109/ICRA.2015.7139235.
  * Medagoda et al. (2011a) Medagoda, L., S. B. Williams, O. Pizarro, and M. V. Jakuba, 2011a: Water column current aided localisation for significant horizontal trajectories with Autonomous Underwater Vehicles. OCEANS’11 - MTS/IEEE Kona, Program Book, 10.23919/oceans.2011.6107240.
  * Medagoda et al. (2011b) Medagoda, L., S. B. Williams, O. Pizarro, and M. V. Jakuba, 2011b: Water column current profile aided localisation combined with view-based SLAM for autonomous underwater vehicle navigation. Proceedings - IEEE International Conference on Robotics and Automation, 3048–3055, 10.1109/ICRA.2011.5980141.
  * Medagoda et al. (2016) Medagoda, L., S. B. Williams, O. Pizarro, J. C. Kinsey, and M. V. Jakuba, 2016: Mid-water current aided localization for autonomous underwater vehicles. Autonomous Robots, 40 (7), 1207–1227, 10.1007/s10514-016-9547-3.
  * Todd et al. (2011) Todd, R. E., D. L. Rudnick, M. R. Mazloff, R. E. Davis, and B. D. Cornuelle, 2011: Poleward flows in the southern California Current System: Glider observations and numerical simulation. Journal of Geophysical Research: Oceans, 116 (2), 1–16, 10.1029/2010JC006536.
  * Todd et al. (2017) Todd, R. E., D. L. Rudnick, J. T. Sherman, W. Brechner Owens, and L. George, 2017: Absolute velocity estimates from autonomous underwater gliders equipped with doppler current profilers. Journal of Atmospheric and Oceanic Technology, 34 (2), 309–333, 10.1175/JTECH-D-16-0156.1.
  * Visbeck (2002) Visbeck, M., 2002: Deep velocity profiling using lowered acoustic Doppler current profilers: Bottom track and inverse solutions. Journal of Atmospheric and Oceanic Technology, 19 (5), 794–807, 10.1175/1520-0426(2002)019¡0794:DVPULA¿2.0.CO;2.


* (25) M. Bonesini, T. Cervi, A. Falcone, U. Kose, R. Mazza, A. Menegolli, C. Montanari, M. Nessi, M. Prata, A. Rappoldi, G. Raselli, M. Rossella, M. Spanu, M. Torti, W. Vollenberg, A. Zani, Journal of Instrumentation 13(12), P12020 (2018).  DOI 10.1088/1748-0221/13/12/p12020.  URL https://doi.org/10.1088%2F1748-0221%2F13%2F12%2Fp12020
  * (26) Polyethylene naphthalate (pen) - film - material information.  http://www.goodfellow.com/E/Polyethylene-naphthalate-Film.html.  Accessed: 2021-10-05
  * (27) D.Mary, et al., J. Phys. D: Appl. Phys. 30, 171 (1997).  DOI 10.1088/0022-3727/30/2/004
  * (28) R. Francini, et al., JINST 8(09), P09006 (2013).  DOI 10.1088/1748-0221/8/09/p09006.  URL https://doi.org/10.1088%2F1748-0221%2F8%2F09%2Fp09006
  * (29) Y. Abraham, et al.  Wavelength-Shifting Performance of Polyethylene Naphthalate Films in a Liquid Argon Environment (2021).  URL arXiv:2103.03232
  * (30) A. Bueno, J. Lozano, A.J. Melgarejo, F.J. Muñoz, J.L. Navarro, S. Navas, A.G. Ruiz, Journal of Instrumentation 3(01), P01006 (2008).  DOI 10.1088/1748-0221/3/01/p01006.  URL https://doi.org/10.1088%2F1748-0221%2F3%2F01%2Fp01006
  * (31) Q. Zhao, M. Guan, P. Zhang, C. Yang, J. Li, Y. Wei, L. Wang, Y. Gan, W. Xiong, Journal of Instrumentation 16(06), T06014 (2021).  DOI 10.1088/1748-0221/16/06/t06014.  URL https://doi.org/10.1088/1748-0221/16/06/t06014
  * (32) C. Benson, et al., Eur. Phys. J. C 78, 329 (2018).  DOI 10.1140/epjc/s10052-018-5807-z
  * (33) J.R. Graybill, C.B. Shahi, M.A. Coplan, A.K. Thompson, R.E. Vest, C.W. Clark, Appl. Opt. 59(4), 1217 (2020).  DOI 10.1364/AO.380185.  URL http://ao.osa.org/abstract.cfm?URI=ao-59-4-1217
  * (34) C.H. Lally, G.J. Davies, W.G. Jones, N.J.T. Smith, (1995)
  * (35) B. Aimard, et al., JINST 16(03), P03007 (2021).  DOI 10.1088/1748-0221/16/03/P03007
  * (36) C.F. Lastoria.  Analysis of the scintillation light production and propagation in the WA105 Dual-Phase demonstrator. Análisis de la producción y propagación de la luz de centelleo en el detector de doble-fase WA105 (2020).  URL http://cds.cern.ch/record/2748990.  Presented 17 Dec 2020
  * (37) E. Segreto, Phys. Rev. C 91(3), 035503 (2015).  DOI 10.1103/PhysRevC.91.035503
  * (38) D. Whittington, et al., JINST 11(05), P05016 (2016).  DOI 10.1088/1748-0221/11/05/P05016
  * (39) E.L. Snider, G. Petrillo, J. Phys. Conf. Ser. 898(4), 042057 (2017).  DOI 10.1088/1742-6596/898/4/042057
  * (40) S. Agostinelli, et al., Nucl. Instrum. Meth. A 506, 250 (2003).  DOI 10.1016/S0168-9002(03)01368-8
  * (41) D. Heck, J. Knapp, J.N. Capdevielle, G. Schatz, T. Thouw, (1998)
  * (42) C. Forti, H. Bilokon, B. d’Ettorre Piazzoli, T.K. Gaisser, L. Satta, T. Stanev, Phys. Rev. D 42, 3668 (1990).  DOI 10.1103/PhysRevD.42.3668.  URL https://link.aps.org/doi/10.1103/PhysRevD.42.3668
  * (43) T.T. Böhlen, F. Cerutti, M.P.W. Chin, A. Fassò, A. Ferrari, P.G. Ortega, A. Mairani, P.R. Sala, G. Smirnov, V. Vlachoudis, Nucl. Data Sheets 120, 211 (2014).  DOI 10.1016/j.nds.2014.07.049
  * (44) B.J.P. Jones, C.S. Chiu, J.M. Conrad, C.M. Ignarra, T. Katori, M. Toups, JINST 8, P07011 (2013).  DOI 10.1088/1748-0221/8/07/P07011.  [Erratum: JINST 8, E09001 (2013)]
  * (45) M. Babicz, et al., JINST 15(09), P09009 (2020).  DOI 10.1088/1748-0221/15/09/P09009
  * (46) E. Grace, A. Butcher, J. Monroe, J.A. Nikkel, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 867, 204 (2017).  DOI https://doi.org/10.1016/j.nima.2017.06.031.  URL https://www.sciencedirect.com/science/article/pii/S0168900217306848
  * (47) M. Antonello, et al., Technical Report ICARUS-TM/06-03 (2006)
  * (48) B. Zatschler, Attenuation of the scintillation light in liquid argon and investigation of the double beta decay of 76Ge into excited states of 76Se in the GERDA experiment.  Ph.D. thesis, Dresden, Tech. U., TU, Dresden (main) (2020)
  * (49) S. Kubota, et al., Phys. Rev. B 20, 3486 (1979).  DOI 10.1103/PhysRevB.20.3486
  * (50) P. Agnes, et al., Phys. Rev. D 97(11), 112005 (2018).  DOI 10.1103/PhysRevD.97.112005
  * (51) E. Segreto, Phys. Rev. D 103(4), 043001 (2021).  DOI 10.1103/PhysRevD.103.043001
  * (52) P. Abratenko, et al., Journal of Instrumentation 16(04), P04004 (2021).  DOI 10.1088/1748-0221/16/04/p04004.  URL https://doi.org/10.1088/1748-0221/16/04/p04004
  * (53) B. Mitrica, R. Margineanu, S. Stoica, M. Petcu, I. Brancus, A. Jipa, I. Lazanu, O. Sima, A. Haungs, H. Rebel, M. Petre, G. Toma, A. Saftoiu, D. Stanca, A. Apostu, C. Gomoiu, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 654(1), 176 (2011).  DOI https://doi.org/10.1016/j.nima.2011.07.006.  URL https://www.sciencedirect.com/science/article/pii/S0168900211014082
  * (54) S. Kubota, M. Hishida, S. Himi, J. Suzuki, J. Ruan, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 327(1), 71 (1993).  DOI https://doi.org/10.1016/0168-9002(93)91413-H.  URL http://www.sciencedirect.com/science/article/pii/016890029391413H
  * (55) A. Neumeier, T. Dandl, T. Heindl, A. Himpsl, L. Oberauer, W. Potzel, S. Roth, S. Schönert, J. Wieser, A. Ulrich, EPL (Europhysics Letters) 109(1), 12001 (2015).  DOI 10.1209/0295-5075/109/12001.  URL https://doi.org/10.1209/0295-5075/109/12001
  * (56) A. Neumeier, T. Dandl, A. Himpsl, L. Oberauer, W. Potzel, S. Schönert, A. Ulrich, EPL (Europhysics Letters) 111(1), 12001 (2015).  DOI 10.1209/0295-5075/111/12001.  URL https://doi.org/10.1209/0295-5075/111/12001
  * (57) R. Acciarri, et al., Nuc. Phys. B - Proc. Supp. 197, 70 (2009).  DOI 10.1016/j.nuclphysbps.2009.10.037