Model-Aided Localization and Navigation for Underwater Gliders Using Single-Beacon Travel-Time Differences

Abstract

An accurate motion model and reliable measurements are required for autonomous underwater vehicle localization and navigation in underwater environments. However, without a propeller, underwater gliders have limited maneuverability and carrying capacity, which brings difficulties for modeling and measuring. In this paper, an extended Kalman filter (EKF)-based method, combining a modified kinematic model of underwater gliders with the travel-time differences between signals received from a single beacon, is proposed for estimating the glider positions in a predict-update cycle. First, to accurately establish a motion model for underwater gliders moving in the ocean, we introduce two modification parameters, the attack and drift angles, into a kinematic model of underwater gliders, along with depth-averaged current velocities. The attack and drift angles are calculated based on the coefficients of hydrodynamic forces and the sensor-measured angle variation over time. Then, instead of satisfying synchronization requirements, the travel-time differences between signals received from a single beacon, multiplied by the sound speed, are taken as the measurements. To further reduce the EKF estimation error, the Rauch-Tung-Striebel (RTS) smoothing method is merged into the EKF system. The proposed method is tested in a virtual spatiotemporal environment from an ocean model. The experimental results show that the performance of the RTS-EKF estimate is improved when compared with the motion model estimate, especially by 46% at the inflection point, at least in the particular study developed in this article.

1. Introduction

Underwater gliders are a type of autonomous underwater vehicles (AUVs) but without a propeller [1], which are great autonomous platforms fitted in a persistent underwater surveillance system [2,3]. Due to the limited carrying capacity restricted by the design of underwater gliders, unless specifically required, underwater gliders are not equipped with sensors such as an inertial measurement unit (IMU) [4], a Doppler velocity log (DVL) or an acoustic Doppler current profiler (ADCP) sonar, which are usually carried by AUVs [5]. This makes the navigation and positioning of underwater gliders more challenging than for AUVs, but obtaining the glider’s position is very critical for tasks such as target detection [6] and sound field construction [7].

A collection of localization and navigation methods have been proposed for AUVs [8,9,10,11,12]. Any navigation algorithm is based on state estimation [8], which is conducted in a predict-update cycle. The prediction process of state estimation is established mainly based on motion models of AUVs. Based on the installed sensors, AUVs usually adopt inertial navigation by using the IMU-detected acceleration and the DVL-measured relative velocity of the AUV [13]. A carried ADCP provides the relative speed of the local current for AUV navigation [14]. The navigation of an underwater glider depends on the glider’s working mechanism and the installed sensors, including electronic compass and conductivity-temperature-depth (CTD) sensors, but without any velocity sensor. Based on the working mechanism, the dynamic motion of underwater gliders can be modeled [15,16,17]. For practical applications, motion models are usually simplified, for example, applying Lanchester’s phugoid assumptions [18,19] which assuming the angle of attack is fixed [20]. In [21], a localization scheme is proposed for underwater gliders by neglecting the influence of ocean currents and by assuming a zero drift angle. Seaglider and the Sea-Wing underwater glider have the ability to estimate the depth-averaged current velocities for engineering applications [22,23], which could be utilized for a glider localization and navigation scheme. Here, we combine a kinematic model of underwater gliders with the estimated attack and drift angles and add the depth-averaged current velocities [24] to model the prediction process of the glider positions in the ocean.

Measurements are needed for the updating process to correct the prediction errors. For AUVs, in addition to relying on the installed angle and velocity sensors, the acquisition of measurements depends mainly on acoustic beacons and geophysical features [9]. Geophysical techniques use external environmental information as a reference, which can be observed by sonar [25,26,27], optical [28] or magnetic [29] sensors fitted to AUVs. Acoustic techniques are focused on different types of acoustic features, such as the travel time [30,31] and time difference of arrival [32,33], resulting from hydrophone signals. Common measuring techniques include the ultrashort baseline (USBL), short baseline (SBL), long baseline (LBL), single-beacon and acoustic modem techniques [8]. For USBL and SBL navigation, a set of transceivers and transponders is required, one mounted on the AUV and the other mounted on the support vehicle [34,35]. For LBL or single-beacon navigation, one-way travel-time ranging can be conducted on the basis of synchronized clocks [36,37]. For acoustic communication-based underwater navigation, an acoustic communication system is needed for both the transmitter and the AUV receiver [38]. To reduce the requirements for sensors on underwater gliders, a single beacon is selected as the measuring technique. Regarding the glider’s movement over time as a virtual array, we adopt the travel-time difference between signals from the beacon received at two adjacent receiving locations as the measuring acoustic feature, which conventionally refers to the time difference of arrivals at a pair of hydrophones in an acoustic array [39,40]. Herein, it is assumed that the acoustic beacon transmits signals at a regular interval of time instead of satisfying synchronization requirements and that the beacon remains stationary.

Various methods have been used to solve the predict-update cycling state estimation, mainly including the least squares [41], extended Kalman filter (EKF) [42,43], unscented Kalman filter [44], extended information filter [45], and particle filter [46,47] methods. Considering that our proposed algorithm can be performed online, we select an EKF method for the state estimation of underwater gliders, which can approximate nonlinear systems by a linearized first-order Taylor series expansion [48]. The EKF is employed to fuse acoustic measurements and predictions of the process model [49]. Specifically, estimating the glider adjacent positions from the process model enables the algorithm to properly model the range-difference measurements, and the glider position predicted by the process model is corrected within the EKF framework by real range-difference measurements.

Smoothing, which can be viewed as the third part for sequential estimation of glider’s positions, differs from filtering in that further measurements are assimilated [50,51]. Therefore a smoothing process intuitively reduces the EKF estimation error. Many forms of smoothing algorithms have been developed, in which the two-filter smoother and the forward-backward smoother are usually used when numerical approximations are required [52]. Forward-backward smoothers process the measurements first by a forward filter to compute the filtered estimate, and then by a backward smoothing pass to determine the smoothing estimate from the forward filtered estimate [53]. The two-filter smoothing is an alternative approach to the forward-backward smoothing, relying on two independent filter procedures: one runs forward in time and another runs backward [54]. However, the backward filter for the two-filter smoother usually cannot be computed in a closed form, which prohibits the use of numerical techniques to approximate the smoothing [55]. Hence we adopt a forward-backward smoothing method, named the Rauch-Tung-Striebel (RTS) smoothing method [56], for smoothing estimation. The RTS smoother has been widely employed in inertial navigation systems [57] and target tracking problems [58,59].

Our main innovations are as follows. First, we establish a modified kinematic model of underwater gliders by introducing the attack and drift angles estimated by hydrodynamic coefficients and angles measured by electronic compass or derived from control system records. Second, the range variation from underwater gliders to a static acoustic beacon, calculated from the travel-time difference of two adjacent signals from the beacon, is set as the measurement for the EKF-based localization and navigation system.

The rest of the paper is organized as follows. We describe the localization and navigation issues of underwater gliders in Section 2. In Section 3, we theoretically derive the attack and drift angles to modify a kinematic model of underwater gliders. In Section 4, an EKF-based localization and navigation system is constructed by introducing the depth-averaged currents to the modified kinematic model and calculating the discrete ranging differences as the measurements. The RTS smoothing method is used to further improve the EKF estimate. A simulation experiment is conducted in Section 5. The state data of underwater gliders are based on the records recorded during an experiment conducted in the South China Sea, including angles measured by electronic compass, depth information measured by CTD, and rudder angles derived from control system records. Spatiotemporal ocean currents are generated from an ocean model. The travel times of acoustic signals are simulated via an acoustic simulator based on a spatiotemporal environmental field from an ocean model. In Section 6, we present discussion and conclusions.

2. Problem Statement

The Sea-Wing underwater glider, developed by the Shenyang Institute of Automation, Chinese Academy of Sciences, is shown in Figure 1, and some of the specifications of the Sea-Wing underwater glider are summarized in

Table 1. Some Specifications of the Sea-Wing Underwater Glider.

Size (m)	Hull diameter $0.22$, vehicle length 2, wing span $1.2$
Weight (kg)	65
Inflection Depth (m)	1000
Cruising Speed (m/s)	$0.25$, maximum $0.5$
Range (km)	>1100
Communications	Iridium and radio communication
Navigation	Global positioning system, altimeter and electronic compass
Science Sensor	CTD and hydrophone

. The processes related to posture adjustment on the Sea-Wing underwater glider include buoyancy, pitch and rudder regulations [60]. The gliders move vertically by adjusting their buoyancy and generate a gliding motion through the ocean via a pair of wings. The driving buoyancy is controlled by altering the volume of the underwater glider by feeding hydraulic oil into or bleeding hydraulic oil out of an external oil bladder. The pitch, which determines the vertical gliding angle, is controlled by moving an internal battery package along the longitudinal axis to change the center of gravity of the glider. The rudder, installed on the tail of the glider, alters the course angle by changing the rudder angle to resist the effects of ocean currents. Therefore, to describe the movement of underwater gliders, we need to involve buoyancy, pitch and rudder regulating processes.

The sensors mounted on the Sea-Wing underwater glider provide an effective tool for measuring the status of the glider. The CTD sensor and electronic compass are the two most important sensors for glider navigation in the ocean environment. Through the depth change of the glider over time as measured by the CTD, the vertical speed of the glider in the inertial coordinate system can be estimated, which reflects the buoyancy regulation. From the measurements of electronic compass, the pitch and rudder regulating processes can be incorporated into the glider motion simulation. Furthermore, to simulate the glider motion in the real ocean environment, the external forces and force variations for the glider also need to be considered. The hydrodynamic model is an excellent tool to simulate and analyze the forces of underwater gliders. Based on the glider-mounted sensors and hydrodynamic parameters, we can simulate the glider motion by an improved kinematic model.

A geodetic global positioning system (GPS) receiver is installed on the underwater glider for implementing glider positioning at the sea surface. However, the gliders cannot rely on sensors on board for positioning while underwater. Figure 2 shows an example of the glider’s trajectories. On a vertical section, we can only acquire the glider’s depth over time, instead of over latitude and longitude coordinates. On a horizontal section, we can localize the glider at the sea surface only.

Based on the proposed kinematic model, we can estimate the position of underwater gliders in still water, whereas the flow influence in the ocean is inevitable. With the depth-averaged current velocity estimation module on the glider, the positions of underwater gliders in the ocean can be roughly estimated. Considering the estimate error of depth-averaged currents and the reality of spatiotemporal ocean currents, additional information and strategies need to be adopted. We focus here on using an acoustic beacon to periodically emit acoustic signals and receiving the emitted signal by a hydrophone installed on the glider. With no need to synchronize the hydrophone time and beacon time, the time difference between two adjacent received signals is utilized to calculate the distance variation from the glider to the beacon, and then the calculated distance variation is set as an additional measurement to estimate the glider location.

3. Modified Kinematic Model for the Sea-Wing Underwater Glider

The Sea-Wing underwater glider [61] is driven by an internal buoyancy regulating system and is steered by a movable battery package and a vertical rudder. In this section, we establish a kinematic model for the Sea-Wing underwater glider by considering hydrodynamic parameters to simulate the movement of the glider accurately.

3.1. Kinematic Model

To describe the motion of the Sea-Wing underwater glider, two coordinate frames are established based on a right-handed coordinate system, including the inertial frame and body frame (Figure 3). The body frame $\mathbf{O}:(\mathbf{x},\mathbf{y},\mathbf{z})$ is established at the buoyancy center of the glider. The $\mathbf{x}$-axis points forward along the longitudinal axis of the glider. The $\mathbf{y}$-axis lies in the wing plane, pointing to the right when viewed along the $\mathbf{x}$ direction. The $\mathbf{z}$-axis is set as $\mathbf{x}\times \mathbf{y}$. The inertial frame $\mathbf{E}:(\xi ,\eta ,\zeta )$ is established at a point in space. A rotation matrix ${\mathbf{R}}_{OE}$ that transforms the body coordinate system into the inertial coordinate system can be defined as

$${\mathbf{R}}_{OE}=\left[\begin{array}{ccc}c\theta c\psi & s\varphi s\theta c\psi -c\varphi s\psi & c\varphi s\theta c\psi +s\varphi s\psi \\ c\theta s\psi & s\varphi s\theta s\psi +c\varphi c\psi & c\varphi s\theta s\psi -s\varphi c\psi \\ -s\theta & s\varphi c\theta & c\varphi c\theta \end{array}\right],$$

(1)

in which $\theta$ is the pitch angle, $\psi$ is the heading angle and $\varphi$ is the roll angle. The simplified notations $c·=cos(·)$ and $s·=sin(·)$ are used. By the electronic compass installed on the underwater glider, $\theta$, $\psi$ and $\varphi$ can be measured in real time. When the glider glides downward, $\theta$ is defined as negative. In addition, $\psi$ is the azimuth relative to the north. When viewed along the $\mathbf{x}$ direction, $\varphi$ is negative to the left.

Because there is no velocity sensor, the glider velocity, defined as $\mathbf{V}={[u,v,w]}^T$ in the body frame, needs to be estimated from the dynamic model and available measurements. Based on the attack angle $\alpha$ and the drift angle $\beta$, the relation among the velocity components can be expressed as

$$tan\alpha =\frac{w}{u},$$

(2)

$$tan\beta =\frac{v}{u}.$$

(3)

The available measurement is the vertical speed of the glider $\dot{\zeta }$ in the inertial coordinate system, which can be calculated through the depth change of the glider measured by the CTD. The relation between the rate of change in position $\mathbf{b}={[\xi ,\eta ,\zeta ]}^T$ in the inertial frame and the velocity $\mathbf{V}$ in the body frame can be written as

$$\dot{\mathbf{b}}={\mathbf{R}}_{OE}\mathbf{V}.$$

(4)

Substituting Equations (2) and (3) into Equation (4) gives

$$\dot{\mathbf{b}}={\mathbf{R}}_{OE}{\mathbf{R}}_V\left|\mathbf{V}\right|$$

(5)

where

$${\mathbf{R}}_V=\left[\begin{array}{c}c\beta c\alpha \\ s\beta \\ c\beta s\alpha \end{array}\right]$$

(6)

and $\left|\mathbf{V}\right|=\sqrt{u^2+v^2+w^2}$. From the vertical speed $\dot{\zeta }$, we can derive

$$\left|\mathbf{V}\right|={\left({\mathbf{R}}_{OE}^{\left(3\right)}{\mathbf{R}}_V\right)}^{-1}\dot{\zeta }.$$

(7)

where ${\mathbf{R}}_{OE}^{\left(i\right)}$ denotes the i-th row of ${\mathbf{R}}_{OE}$ and the same denotation is used for the other matrices. Then, the horizontal components of $\dot{\mathbf{b}}$ can be calculated from

$${\dot{\mathbf{b}}}^{(1,2)}={\mathbf{R}}_{OE}^{(1,2)}{\mathbf{R}}_V{\left({\mathbf{R}}_{OE}^{\left(3\right)}{\mathbf{R}}_V\right)}^{-1}\dot{\zeta }.$$

(8)

Therefore, $\alpha$ and $\beta$ are the determining factors for this solution, which can be estimated from the dynamic model of the Sea-Wing underwater glider.

3.2. Solution of Attack Angle

The attack angle, $\alpha$, affects the vertical velocity of the glider. Figure 4 shows the forces of the glider when moving in the vertical plane. Here, $\gamma$ is the gliding angle, $U=\sqrt{u^2+w^2}$ is the gliding speed, L is the lift force, D is the drag force, and $\mathsf{\Delta }B$ is the residual buoyancy difference.

The steady-state dynamic model of the glider in the vertical plane can be expressed as

$$\mathsf{\Delta }Bcos\gamma =-L,$$

(9)

$$\mathsf{\Delta }Bsin\gamma =-D,$$

(10)

in which $\gamma =\theta +\alpha$. In addition, $\alpha$ is defined in the same direction as $\theta$. The hydrodynamic lift and drag forces are modeled in [62] as

$$L=-(K_{L0}+K_L\alpha )U^2,$$

(11)

$$D=-(K_{D0}+K_D{\alpha }^2)U^2,$$

(12)

where $K_{L0}$ and $K_L$ are the lift coefficients and $K_{D0}$ and $K_D$ are the drag coefficients. Then, we can derive the relation between $\theta$ and $\alpha$ as

$$tan(\theta +\alpha )=\frac{D}{L}=\frac{K_{D0}+K_D{\alpha }^2}{K_{L0}+K_L\alpha }.$$

(13)

the lift and drag coefficients can be estimated using computational fluid dynamics (CFD) technology. Specifically, for the Sea-Wing underwater glider, $K_{L0}=0.0105$, $K_L=496.8596$, $K_{D0}=6.9650$, and $K_D=439.7317$.

Because $\theta$ can be measured by the electronic compass in real time during the gliding process, $\alpha$ can be calculated based on Equation (13). An example is shown in Figure 5. Since $\alpha$ is not zero, the actual gliding angle is not equal to $\theta$ which is usually used for dead reckoning. The sharp change of the angle corresponds to the unsteady state of the glider at the beginning of the gliding cycle and near the inflection point due to state transition. In this case, the measured $\theta$ and calculated $\alpha$ are not reliable, so we actually set the angles during the period of unsteady state to interpolation results of adjacent steady-state points.

3.3. Solution of Drift Angle

The drift angle, $\beta$, affects the horizontal velocity of the glider. Figure 6 shows the variables for motion analysis in the horizontal plane. Here, $\vartheta =\psi +\beta$ is the course angle, $C=\sqrt{u^2+v^2}$ is the horizontal speed, and ${\delta }_r$ is the rudder angle. Rudder rotation changes the direction of the horizontal force, which in turn changes the course of the glider. Based on the analysis for a maneuvering submarine [63], the glider motion with respect to the body coordinate system in the horizontal plane can be represented by

$$\begin{array}{cc}\hfill m(\dot{u}-vr)& =\frac{1}{2}\rho l^4{X}_{rr}^{{}^{\prime }}{r}^2+\frac{1}{2}\rho l^3(X_{\dot{u}}^{{}^{\prime }}\dot{u}+X_{vr}^{{}^{\prime }}vr)\hfill \\ \hfill & +\frac{1}{2}\rho l^2(X_{uu}^{{}^{\prime }}{u}^2+X_{vv}^{{}^{\prime }}{v}^2+X_{{\delta }_r{\delta }_r}^{{}^{\prime }}{u}^2{\delta }_r^2),\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill m(\dot{v}+ur)& =\frac{1}{2}\rho l^4(Y_{\dot{r}}^{{}^{\prime }}\dot{r}\hfill \\ \hfill & +Y_{r\left|r\right|}^{{}^{\prime }}r\left|r\right|)+\frac{1}{2}\rho l^3(Y_{\dot{v}}^{{}^{\prime }}\dot{v}+Y_r^{{}^{\prime }}ur+Y_{v\left|r\right|}^{{}^{\prime }}v\left|r\right|)\hfill \\ \hfill & +\frac{1}{2}\rho l^2(Y_0^{{}^{\prime }}{u}^2+Y_v^{{}^{\prime }}uv+Y_{v\left|v\right|}^{{}^{\prime }}v\left|v\right|+Y_{{\delta }_r}^{{}^{\prime }}{u}^2{\delta }_r),\hfill \end{array}$$

(15)

$$\begin{array}{cc}\hfill I_{zz}\dot{r}& =\frac{1}{2}\rho l^5(N_{\dot{r}}^{{}^{\prime }}\dot{r}+N_{r\left|r\right|}^{{}^{\prime }}r\left|r\right|)\hfill \\ \hfill & +\frac{1}{2}\rho l^4(N_{\dot{v}}^{{}^{\prime }}\dot{v}+N_r^{{}^{\prime }}ur+N_{v\left|r\right|}^{{}^{\prime }}v\left|r\right|)\hfill \\ \hfill & +\frac{1}{2}\rho l^3(N_0^{{}^{\prime }}{u}^2+N_v^{{}^{\prime }}uv+N_{v\left|v\right|}^{{}^{\prime }}v\left|v\right|+N_{{\delta }_r}^{{}^{\prime }}{u}^2{\delta }_r),\hfill \end{array}$$

(16)

where m is the glider mass; $\rho$ is the seawater density; l is the glider length; r is the angular velocity about the z-axis; $\dot{r}$ is the angular acceleration about the z-axis; $I_{zz}$ is the moment about the z-axis; $\dot{u}$ and $\dot{v}$ are the accelerations in the x and y directions, respectively; and $X_{rr}^{{}^{\prime }},Y_r^{{}^{\prime }},N_v^{{}^{\prime }},\cdots$ are nondimensional coefficients of hydrodynamic forces. When the glider is in the steady rotation phase, $\dot{u}=\dot{v}=\dot{r}=0$, and Equations (14)–(16) can be expressed as

$$\begin{array}{cc}\hfill 0& =\frac{1}{2}\rho l^4{Y}_{r\left|r\right|}^{{}^{\prime }}r\left|r\right|+\frac{1}{2}\rho l^3(Y_r^{{}^{\prime }}ur+Y_{v\left|r\right|}^{{}^{\prime }}v\left|r\right|)\hfill \\ \hfill & +\frac{1}{2}\rho l^2(Y_v^{{}^{\prime }}uv+Y_{v\left|v\right|}^{{}^{\prime }}v\left|v\right|+Y_{{\delta }_r}^{{}^{\prime }}{u}^2{\delta }_r)-mur\hfill \end{array}$$

(17)

$$\begin{array}{cc}\hfill 0& =\frac{1}{2}\rho l^5{N}_{r\left|r\right|}^{{}^{\prime }}r\left|r\right|+\frac{1}{2}\rho l^4(N_r^{{}^{\prime }}ur+N_{v\left|r\right|}^{{}^{\prime }}v\left|r\right|)\hfill \\ \hfill & +\frac{1}{2}\rho l^3(N_v^{{}^{\prime }}uv+N_{v\left|v\right|}^{{}^{\prime }}v\left|v\right|+N_{{\delta }_r}^{{}^{\prime }}{u}^2{\delta }_r)\hfill \end{array}$$

(18)

According to the assumptions of the linear motion equation, Equations (17) and (18) can be simplified as

$$\left\{\begin{array}{c}\frac{1}{2}\rho l^2{Y}_v^{{}^{\prime }}uv-(mu-\frac{1}{2}\rho l^3{Y}_r^{{}^{\prime }}u)r=-\frac{1}{2}\rho l^2{Y}_{{\delta }_r}^{{}^{\prime }}{u}^2{\delta }_r\\ \frac{1}{2}\rho l^3{N}_v^{{}^{\prime }}uv+\frac{1}{2}\rho l^4{N}_r^{{}^{\prime }}ur=-\frac{1}{2}\rho l^3{N}_{{\delta }_r}^{{}^{\prime }}{u}^2{\delta }_r\end{array}\right.$$

(19)

Assuming that $\beta$ is small, Equation (3) can be approximated as

$$\beta =v/u.$$

(20)

Then, the solution of Equation (19) can be derived as

$$\beta =-\frac{N_{{\delta }_r}^{{}^{\prime }}(m^{{}^{\prime }}-Y_r^{{}^{\prime }})+N_r^{{}^{\prime }}{Y}_{{\delta }_r}^{{}^{\prime }}}{N_v^{{}^{\prime }}(m^{{}^{\prime }}-Y_r^{{}^{\prime }})+N_r^{{}^{\prime }}{Y}_v^{{}^{\prime }}}·{\delta }_r,$$

(21)

where $m^{{}^{\prime }}=m/\left(0.5\rho l^3\right)$ is the nondimensional weight. For the Sea-Wing underwater glider, ${\delta }_r\in [-{30}^{\circ }{30}^{\circ }]$. Here, ${\delta }_r=0^{\circ }$ represents the $\mathbf{x}$ direction, and when viewed along the $\mathbf{x}$ direction, ${\delta }_r$ is positive to the left. ${\delta }_r$ is continuously changed during the gliding process to resist the ocean current, to adjust the course angle of the glider, and to keep the preset trajectory as far as possible.

Using CFD software Fluent [64], we can simulate the solutions of $\beta$ under different ${\delta }_r$. The specific numerical values of the hydrodynamic coefficients used in the simulation are shown in

Table 2. Nondimensional coefficients of hydrodynamic forces used for simulations.

Coefficients	$Y_{r\left|r\right|}^{{}^{\prime }}$	$Y_r^{{}^{\prime }}$	$Y_{v\left|r\right|}^{{}^{\prime }}$	$Y_v^{{}^{\prime }}$	$Y_{v\left|\mathit{v}\right|}^{{}^{\prime }}$	$Y_{{\delta }_r}^{{}^{\prime }}$
Values	0.00611	0.01065	−0.03931	−0.03545	0.01558	−0.00968
Coefficients	$N_{r\left|r\right|}^{{}^{\prime }}$	$N_r^{{}^{\prime }}$	$N_{v\left|r\right|}^{{}^{\prime }}$	$N_v^{{}^{\prime }}$	$N_{v\left|v\right|}^{{}^{\prime }}$	$N_{{\delta }_r}^{{}^{\prime }}$
Values	−0.00303	−0.00523	0.01312	−0.00149	−0.01983	0.00368

. Figure 7a shows both the nonlinear results from Equations (17) and (18) and the linear results from Equation (21). As the absolute value of ${\delta }_r$ decreases, the difference between the nonlinear and linear results gradually decreases. Based on real-time adjustment of the glider, ${\delta }_r$ can be acquired to estimate the corresponding value of $\beta$, and an example is shown in Figure 7b.

3.4. Comparison with Dead Reckoning

For comparison purpose, we compute the dead-reckoning navigation results. The dead reckoning method takes no consideration of $\alpha$ and $\beta$, that is $\alpha =0$ and $\beta =0$. Then from Equation (8), the velocity of the underwater glider, ${\dot{\mathbf{b}}}_{DR}$, could be calculated by

$${\dot{\mathbf{b}}}_{DR}=\left[\begin{array}{c}-\dot{\zeta }cot\theta c\psi \\ -\dot{\zeta }cot\theta s\psi \\ \dot{\zeta }\end{array}\right].$$

(22)

Figure 8 shows an example of dead reckoning result compared with the result of modified kinematic model. The deviation between the preset trajectory and the actual trajectory, represented by the connection between the actual start and end positions, is mainly caused by ocean currents. The estimated trajectory of modified kinematic model is closer to the preset trajectory than the dead reckoning trajectory, indicating that the proposed kinematic model reflects the resistance of the glider to the ocean current through angle adjustment. Then the navigation result of the kinematic model superimposed ocean current can be closer to the real situation.

4. EKF-Based Localization and Navigation System Modeling

Based on the established kinematic model in Section 3, we can estimate the moving trajectories of a Sea-Wing underwater glider in still water. With the estimation module of depth-averaged current velocities on the glider, the position of an underwater glider in the ocean can be approximately estimated. The basic approach for estimating the depth-averaged current velocity is to use the difference between the dead-reckoning and GPS positions of resurfacing when a gliding cycle ends, divided by the gliding time [24]. The calculated depth-averaged current velocity is not accurate due to the dead-reckoning navigation error. Even if we adopt a more accurate model, such as the proposed kinematic model, the estimated positions of underwater gliders using depth-averaged current velocities may still be unreliable because of the reality of spatiotemporal ocean currents. Additional acoustic measurements are required to calibrate the estimated positions.

In this section, an EKF method is imported to recursively estimate the positions of the glider during a gliding cycle. The prediction process is based on the modified kinematic model and depth-averaged current velocities, and the update process is based on the distance variation from the glider to the beacon. The RTS smoothing algorithm is adopted to improve the EFK estimation by introducing subsequent measurements.

4.1. System State Prediction

Because the glider depth, $\zeta$, can be measured by the installed CTD, the components of glider positions that need to be estimated are ${(\xi ,\eta )}^T$, that is, the east and south components in the horizontal plane. We set the system state vector, $\mathbf{s}$, as below,

$$\mathbf{s}={(\xi ,\eta )}^T.$$

Based on Equation (8) and the depth-averaged current velocity, the process model for the system satisfies the following motion equation:

$$\left[\begin{array}{c}\dot{\xi }\\ \dot{\eta }\end{array}\right]=\left[\begin{array}{c}{\mathbf{R}}_{OE}^{\left(1\right)}{\mathbf{R}}_V{\left({\mathbf{R}}_{OE}^{\left(3\right)}{\mathbf{R}}_V\right)}^{-1}\dot{\zeta }+v_x\\ {\mathbf{R}}_{OE}^{\left(2\right)}{\mathbf{R}}_V{\left({\mathbf{R}}_{OE}^{\left(3\right)}{\mathbf{R}}_V\right)}^{-1}\dot{\zeta }+v_y\end{array}\right]+\left[\begin{array}{c}{g}_1\\ g_2\end{array}\right]w_s$$

(23)

that is,

$$\dot{\mathbf{s}}\left(t\right)=\mathbf{f}(\mathbf{s}\left(t\right),\mathbf{u}\left(t\right),\mathbf{v}\left(t\right),t)+{\mathbf{G}}_p{\mathbf{w}}_p\left(t\right)$$

(24)

where

$${\mathbf{G}}_p=\left[\begin{array}{cc}{g}_1& 0\\ 0& g_2\end{array}\right],{\mathbf{w}}_p\left(t\right)=\left[\begin{array}{cc}{w}_s\left(t\right)& 0\\ 0& w_s\left(t\right)\end{array}\right],$$

$\mathbf{v}\left(t\right)={(v_x,v_y)}^T$ and $\mathbf{u}\left(t\right)={(\alpha ,\beta ,\theta ,\psi ,\varphi ,{\delta }_r,\dot{\zeta })}^T$. $\mathbf{u}\left(t\right)$ is the control input for the modified kinematic model. $v_x$ and $v_y$ are the depth-averaged current components along the $\xi$ and $\eta$ directions, respectively. ${\mathbf{w}}_p\left(t\right)$ is the zero-mean Gaussian process noise, representing accelerations that allow the ocean current to deviate from the depth-averaged current. ${\mathbf{G}}_p$ can be considered as a weight matrix, which determines the difference of process noise in different directions.

Then, the propagation of the system state estimate $\widehat{\mathbf{s}}\left(t\right)$ and error covariance $\mathbf{P}\left(t\right)$ using the EKF [48] can be written as

$$\dot{\widehat{\mathbf{s}}}\left(t\right)=\mathbf{f}(\widehat{\mathbf{s}}\left(t\right),\mathbf{u}\left(t\right),\mathbf{v}\left(t\right),t)$$

(25)

$$\dot{\mathbf{P}}\left(t\right)=\mathbf{F}\left(t\right)\mathbf{P}\left(t\right)+\mathbf{P}\left(t\right){\mathbf{F}}^T\left(t\right)+{\mathbf{G}}_p\mathbf{Q}{\mathbf{G}}_p^T,$$

(26)

where $\mathbf{F}\left(t\right)$ is a system dynamic matrix,

$$\mathbf{F}\left(t\right)\equiv \frac{\partial \mathbf{f}}{\partial \mathbf{s}}{|}_{\widehat{\mathbf{s}}\left(t\right)}=\left[\begin{array}{cc}0& 0\\ 0& 0\end{array}\right],$$

(27)

and the process-noise matrix $\mathbf{Q}$ can be written as

$$\mathbf{Q}=\left[\begin{array}{cc}{\mathsf{\Phi }}_s& 0\\ 0& {\mathsf{\Phi }}_s\end{array}\right]$$

(28)

under the consideration that each direction of the ocean current accelerations is uncorrelated. ${\mathsf{\Phi }}_s$ is the spectral density of the white noise $w_s$.

4.2. Measurement Model

A fixed beacon with a known position is needed to transmit signals at regular intervals. The Sea-Wing underwater glider equipped with hydrophones can receive the acoustic signal emitted from the preset acoustic beacon. Assuming that the emitted acoustic signal propagates along a straight line, the distance between the receiver and the beacon can be calculated by the product of the travel-time measurement and the sound speed in water. Time synchronization between the receiver and the beacon is necessary in order to obtain an accurate travel time. To overcome this constraint, we propose here to utilize the time difference between two adjacent receiving locations. Specifically, as long as the beacon transmits a signal according to a preset time interval $\mathsf{\Delta }t$, the receiver can calculate the travel-time difference according to the receiving time of two adjacent signals, thereby calculating the distance difference.

The range from the beacon to the receiver is a nonlinear function of the receiver location and the beacon location. At a time point, $t_k(k=1,2,\dots )$, the location of the receiver, that is, the glider, is ${\mathbf{s}}_k={({\xi }_k,{\eta }_k)}^T$, and the beacon location is ${\mathbf{s}}_{bk}={({\xi }_{bk},{\eta }_{bk})}^T$. The perfect range measurement can be written as

$$R_k=\sqrt{{({\mathbf{s}}_k-{\mathbf{s}}_{bk})}^T({\mathbf{s}}_k-{\mathbf{s}}_{bk})}.$$

(29)

The distance difference corresponding to two adjacent receiving times, $t_k$ and $t_{k-j}(k\ge K+1,j=1,\dots ,K)$, can be modeled as

$$\mathsf{\Delta }{R}_k^{\left(j\right)}=R_k-R_{k-j}.$$

(30)

Then considering the actual noisy range difference measurement, ${\tilde{\mathbf{y}}}_k$, we can set the measurement model as

$${\tilde{\mathbf{y}}}_k=h\left({\mathbf{s}}_k\right)+n_k$$

(31)

where

$${\tilde{\mathbf{y}}}_k^{\left(j\right)}=\mathsf{\Delta }{R}_k^{\left(j\right)}+n_k^{\left(j\right)}=c(t_k-t_{k-j}-j\mathsf{\Delta }t),$$

(32)

$c=1500$ m/s, and K is the virtual array size set by experience. $n_k\sim \mathcal{N}(0,N_r)$ is the measurement noise, which is in units of distance and represents the imprecision in timing multiplied by c. The variance $N_r$ is assumed to be independent of time. The Jacobian of the range difference measurement, Equation (31), at $t_k$ with respect to $\mathbf{s}$ is ${\mathbf{H}}_k\equiv \frac{\partial h}{\partial \mathbf{s}}{|}_{{\widehat{\mathbf{s}}}_k}$ with

$${\mathbf{H}}_k^{\left(j\right)}={\left[\begin{array}{c}\frac{{\xi }_k-{\xi }_{bk}}{R_k}-\frac{{\xi }_{k-j}-{\xi }_{b(k-j)}}{R_{k-j}}\\ \frac{{\eta }_k-{\eta }_{bk}}{R_k}-\frac{{\eta }_{k-j}-{\eta }_{b(k-j)}}{R_{k-j}}\end{array}\right]}^T.$$

(33)

4.3. Recursive Estimation of the System State

EKF methods can operate recursively on noisy input data to produce a statistically suboptimal estimate of the system state [65]. When the range difference measurements are available, the system states are updated. If no measurement is acquired, the system states are predicted by the propagation model, Equation (25), which performs over the sampling step $T_s$.

Given the state vector at time point $t_{i-1}=(i-1)T_s,(i=2,3,\dots )$, a priori state estimation ${\widehat{\mathbf{s}}}_i^{-}$ by integrating Equation (25) is performed. If at $t_i$, the k-th $(k=K+1,K+2,\dots )$ beacon emitted signal is received, then the range-difference measurements ${\tilde{\mathbf{y}}}_k$ are available. From the a priori state estimation ${\widehat{\mathbf{s}}}_i^{-}$ and the estimated state at previous time point $t_{k-j}$, we model the range-difference measurements by Equations (29)–(31). By weighting the deviation between the model-calculated range differences $h\left({\widehat{\mathbf{s}}}_i^{-}\right)$ and the true range-difference measurements ${\tilde{\mathbf{y}}}_k$, we can correct the a priori state estimate ${\widehat{\mathbf{s}}}_i^{-}$ to acquire an a posteriori state estimate ${\widehat{\mathbf{s}}}_i^{+}$ using

$${\widehat{\mathbf{s}}}_i^{+}={\widehat{\mathbf{s}}}_i^{-}+{\mathbf{K}}_k[{\tilde{\mathbf{y}}}_k-h\left({\widehat{\mathbf{s}}}_i^{-}\right)],$$

(34)

where

$${\mathbf{K}}_k={\mathbf{P}}_i^{-}{\mathbf{H}}_k^T\left({\widehat{\mathbf{s}}}_i^{-}\right){[{\mathbf{H}}_k\left({\widehat{\mathbf{s}}}_i^{-}\right){\mathbf{P}}_i^{-}{\mathbf{H}}_k^T\left({\widehat{\mathbf{s}}}_i^{-}\right)+N_r]}^{-1}$$

(35)

is the Kalman gain corresponding to the k-th measurement and ${\mathbf{P}}_i^{-}$ is a priori estimate of the error covariance matrix calculated from Equation (26). The Kalman gain ${\mathbf{K}}_k$ is regarded as a weighting function to make the EKF algorithm recursive. Predicted system states with smaller uncertainty, that is, ${\widehat{\mathbf{s}}}_i^{-}$ with smaller ${\mathbf{P}}_i^{-}$, are given more weight. The weighted result ${\widehat{\mathbf{s}}}_i^{+}$ and its covariance, which is calculated by inform the prediction used in the following time step.

$${\mathbf{P}}_i^{+}=[\mathbf{I}-{\mathbf{K}}_k{\mathbf{H}}_k\left({\widehat{\mathbf{s}}}_i^{-}\right)]{\mathbf{P}}_i^{-},$$

(36)

This online recursive estimation process for a diving cycle of underwater gliders starts from the initial position at the sea surface and continues until the glider rises to the surface again. The purpose of this process is to make the glider know its position and then to better detect the target [6]. The limitation of this EKF system is that the state estimation error introduced by the inaccurate depth-averaged current velocity, which is estimated from the previous gliding cycle, cannot be fully corrected and the cumulative error will increase over time.

4.4. Estimation Improvement by RTS Smoothing

To reduce the EKF estimation error, we introduce the RTS smoothing method to the EKF system. Measurements acquired at all available future and past time points are used to estimate the system state. Because additional measurements are utilized in the estimation algorithm, the RTS smoothing intuitively gives better estimates [49,66] to support the spatiotemporal field construction of target radiation signal [7]. Smoothing is the postprocessing process. We can calculate the depth-averaged current velocity for the current gliding cycle from the GPS-localized resurfacing position of the glider, which further reduces the state estimation error.

Let ${\widehat{\mathbf{s}}}_s$ denote the smoothed estimation. Based on the EKF estimated results ${\widehat{\mathbf{s}}}_i^{+}$, the RTS smoothed result ${\widehat{\mathbf{s}}}_{si}^{+}$ can be estimated by the dynamic model, ${\mathbf{F}}_i$, with a correction term, ${\mathbf{K}}_{si}$, as below,

$${\widehat{\mathbf{s}}}_{si}^{+}={\widehat{\mathbf{s}}}_{si}^{-}-({\mathbf{F}}_i+{\mathbf{K}}_{si})({\widehat{\mathbf{s}}}_{si}^{-}-{\widehat{\mathbf{s}}}_i^{+}),$$

(37)

where

$${\widehat{\mathbf{s}}}_{si}^{-}={\widehat{\mathbf{s}}}_{s(i+1)}^{+}-\mathbf{f}({\widehat{\mathbf{s}}}_s\left(t\right),\mathbf{u}\left(t\right),\mathbf{v}\left(t\right),t){|}_{t_{i+1}}(t_{i+1}-t_i),$$

(38)

and

$${\mathbf{K}}_{si}\equiv {\mathbf{G}}_p\mathbf{Q}{\mathbf{G}}_p^T{\left({\mathbf{P}}_i\right)}^{-1}$$

(39)

is the smoothed gain. When ${\tilde{\mathbf{y}}}_k$ is available, ${\mathbf{P}}_i={\mathbf{P}}_i^{+}$; otherwise, ${\mathbf{P}}_i={\mathbf{P}}_i^{-}$.

Although for the RTS smoothing process, we use the depth-averaged current velocity in the current gliding cycle, there is still a gap with the spatiotemporal flow field that actually affects the position of the glider. To alleviate the problems caused by accumulation errors, we divide a gliding cycle into two parts by taking the inflection point as the boundary. For the second half of the gliding cycle, that is, after the inflection point, the RTS smoothing method is used to perform the estimation on the EKF results from the GPS-localized resurfacing position. For the first half, we still use the EKF results from the initial GPS-localized position with the depth-averaged current velocity in the current gliding cycle. Therefore, it is inevitable that a larger error within the gliding cycle appears near the inflection point. We call this combined estimate the RTS-EKF Estimate.

5. Simulation Based on Experimental and Model Data

Simulations show the experimental scene of a Sea-Wing underwater glider in a sea region with a submersible mooring beacon. The simulation site is shown in Figure 9. The number of gliding cycles for the entire operation is 16. The real trajectories of the underwater glider were generated by our proposed kinematic model and simulated ocean currents. The angle and depth data for the kinematic model in the first gliding cycle were based on records recorded during an experiment conducted in the South China Sea because we had no measurement model for the electronic compass and CTD sensor installed on the underwater glider, while the data used in subsequent gliding cycles were based on the conversion of the adopted data in the first cycle. Ocean currents in the South China Sea experimental site had no influence on the accuracy of these adopted angle and depth data. The acoustic data for the EKF-based localization and navigation algorithm were generated by an acoustic simulator.

5.1. Simulated Ocean Currents

To simulate the motion of an underwater glider in a flow field, we generated a 3D spatiotemporal flow field based on the hydrological basic data acquired from an ocean model, named POM [67]. The generated flow field was combined with the modified kinematic model of the glider for real position setting. Figure 10 shows the flow along the glider trajectory, which reflects that the simulated spatiotemporal flow field varies drastically during the glider operation.

5.2. Acoustic Travel-Time Simulation

A virtual acoustic beacon was deployed within the operation region. To simulate the transmission of acoustic signals emitted by the acoustic beacon, we set the environmental parameters. The bathymetry data in the region are obtained from geospatial data at the NOAA website [68]. The sound speed distribution in this region is calculated from the basic hydrological data acquired from a POM South China Sea $1/{15}^{\circ }$ analysis provided by the South China Sea Institute of Oceanology, Guangzhou, P.R. China [67]. The acoustic travel time of the emitted signals is generated using the toolbox Bellhop 3D, which is published by the Ocean Acoustics Library [69].

The transmission time interval of the acoustic beacon, $\mathsf{\Delta }t$, is set to 10 s. However, considering the simulation situation, we need to obtain the acoustic travel time based on the glider position. Therefore, we set the receiving interval for the underwater glider to 10 s instead of the beacon emitting interval. The travel time of the direct-path wave can be acquired for each receiving location. Then, the travel-time differences between adjacent receiving locations can be calculated, which are taken as the system measurements by multiplying the sound velocity. So instead of calculating the noisy range-difference measurement ${\tilde{\mathbf{y}}}_k$ by Equation (32), in simulations, ${\tilde{\mathbf{y}}}_k$ is acquired by

$${\tilde{\mathbf{y}}}_k^{\left(j\right)}=c({\tau }_k-{\tau }_{k-j}),$$

(40)

where ${\tau }_k(k\ge K+1,j=1,\dots ,K)$ is the simulated one-way travel time corresponding to the receiving time $t_k$.

5.3. EKF Estimation

The EKF needs to be initialized on startup. The initial value of the system state ${\widehat{\mathbf{s}}}_0$ is set to the real GPS location. The depth-averaged current velocity for the first gliding cycle is set to $v_x{{|}_0=v_y|}_0=0$. For the m-th $(m=2,3,\dots )$ gliding cycle, the adopted depth-averaged current velocity is calculated from the $(m-1)$-th gliding cycle, that is, $(v_x{|}_{m-1},v_y{{|}_{m-1})}^T$. Considering the actual situation, we introduce zero-mean Gaussian noise $\mu \sim \mathcal{N}(0,N_v)$ into the real depth-averaged current velocity. Here, $N_v={\left(0.05\mathrm{m}/\mathrm{s}\right)}^2$.

The combined effect of ${\mathbf{G}}_p$ and $\mathbf{Q}$ explains the error induced by the depth-averaged current velocity. We set $\mathbf{Q}$ to be a constant matrix with ${\mathsf{\Phi }}_s=2\times {\left(0.2\mathrm{m}/\mathrm{s}\right)}^2$ for all gliding cycles. ${\mathbf{G}}_p$, regarded as a weighting matrix, is assigned varying values for different gliding cycles based on the acceleration of the depth-averaged current velocity. For the first cycle, without any prior knowledge of the current, we set $g_1=g_2=0.5$. For the m-th $(m=2,3,\dots )$ gliding cycle, the weights $g_1$ and $g_2$ are set to

$$g_1=\frac{a_x{|}_{m-1}}{{a|}_{m-1}},\mathrm{and}{g}_2=\frac{a_y{|}_{m-1}}{{a|}_{m-1}},$$

(41)

where

$$a_x{|}_{m-1}=\frac{\left|v_x{{|}_{m-1}-v_x|}_{m-2}\right|}{T_{m-1}},$$

(42)

$$a_y{|}_{m-1}=\frac{\left|v_y{{|}_{m-1}-v_y|}_{m-2}\right|}{T_{m-1}},$$

(43)

$${a|}_{m-1}=a_x{{|}_{m-1}+a_y|}_{m-1},$$

(44)

and $T_{m-1}$ is the gliding period for the $(m-1)$-th gliding cycle.

The measurement noise variance $N_r$ also needs to be given. We set $N_r={\left(50\mathrm{m}\right)}^2$ to account for the uncertainty during ranging based on the time difference multiplied by c.

Because the calculation of depth-averaged current velocity does not rely on the EKF process, we use the modified kinematic model with the superimposed depth-averaged current velocity for comparison, labeled as “Motion Model Estimate” in Figure 11. The trajectory estimated by the EKF method is closer to the real trajectory than the motion model estimate. At the resurfacing position, the estimated error is reduced from 800 m to 630 m. However the estimated trajectory is still far from the actual trajectory.

5.4. RTS-EKF Estimation

The RTS backward smoothing algorithm is applied to improve the EKF estimate. For postprocessing, the resurfacing position of the glider is known. Then, the depth-averaged current velocity for the current gliding cycle can be calculated. Therefore, for the m-th $(m=1,2,\dots )$ gliding cycle, the adopted depth-averaged current velocity is $(v_x{|}_m,v_y{{|}_m)}^T$. The weights $g_1$ and $g_2$ are assigned to

$$g_1=\frac{a_x{|}_m}{{a|}_m},\mathrm{and}{g}_2=\frac{a_y{|}_m}{{a|}_m},$$

(45)

where

$$a_x{|}_m=\frac{\left|v_x{{|}_m-v_x|}_{m-1}\right|}{T_m},$$

(46)

$$a_y{|}_m=\frac{\left|v_y{{|}_m-v_y|}_{m-1}\right|}{T_m},$$

(47)

$${a|}_m=a_x{{|}_m+a_y|}_m,$$

(48)

and $T_m$ is the gliding period for the m-th gliding cycle.

Considering the accumulation of errors, we use the RTS smoothing from Equations (37)–(39) for the second half of the gliding cycle, that is, after the inflection point, and for the first half of the gliding cycle, we keep the EKF estimate with $(v_x{|}_m,v_y{{|}_m)}^T$. The performance of the estimate state $\widehat{\mathbf{s}}$ is evaluated by the root-mean-square error (RMSE) as follows:

$$RMSE=\sqrt{\frac{1}{N}{(\widehat{\mathbf{s}}-\mathbf{s})}^T(\widehat{\mathbf{s}}-\mathbf{s})},$$

(49)

where N is the number of trials. Figure 12 shows the estimate results with $N=1$. For comparison, the motion model estimate and the EKF estimate for the second half of the gliding cycle is calculated by integrating the dynamic model in reverse time from the resurfacing position.

Because we divide the gliding cycle into two parts, a larger error within the gliding cycle appears near the inflection point. At the inflection point, the RMSE of the RTS-EKF estimate is reduced from $211.4$ m to $114.1$ m compared with the EKF result, and when compared with the motion model estimate, the performance of the RTS-EKF estimate is improved by $46\%$. The reason for RTS performance improvement is that the cumulative estimation error of the EKF from the starting point is corrected, and new measurements up to the resurfacing position are introduced.

5.5. Discussion

The main error induced to the model-aided estimation process is due to the depth-averaged current velocity. Therefore, except the spatiotemporal characteristics of the ocean current, the error $\mu$ of the estimated depth-averaged current velocity impacts the estimate performance. We discussed one trial for the case of $N_v={\left(0.05\mathrm{m}/\mathrm{s}\right)}^2$ in Section 5.4. Here, multiple trials are conducted for the case of $N_v={\left(0.05\mathrm{m}/\mathrm{s}\right)}^2$ and for other cases as $N_v=(0,0.1,0.15,0.2){\left(\mathrm{m}/\mathrm{s}\right)}^2$ to explore the statistical performance of the RTS-EKF estimate.

Figure 13 shows the RMSE results over different $N_v$. The specific percentages are given in

Table 3. RMSE comparison over different errors of the depth-averaged current velocity.

$\sqrt{{\mathit{N}}_{\mathit{v}}}$ m/s	RTS-EKF Estimate			Motion Model Estimate
	≤100 m	100 m$<·<$500 m	≥500 m	≤100 m	100 m$<·<$500 m	≥500 m
0	100%	0	0	100%	0	0
0.05	62%	38%	<1%	51%	49%	<1%
0.1	39%	52%	9%	27%	64%	9%
0.15	30%	50%	20%	21%	55%	24%
0.2	22%	49%	30%	15%	47%	38%

. For each $N_v$, we calculate the percentages of the RMSE values within different error ranges for all gliding cycles. The motion model estimate is shown for comparison. As $N_v$ increases, the percentage for $RMSE\le$ 100 m gradually decreases, while the percentage for $RMSE\ge$ 500 m gradually increases. Under the same case of $N_v$, the percentage for $RMSE\le$ 100 m corresponding to the RTS-EKF estimate is always larger than that of the motion model estimate, while the percentage for $RMSE\ge$ 500 m corresponding to the RTS-EKF estimate is always not greater than that of the motion model estimate.

6. Conclusions

This paper has formulated a dynamic localization and navigation method for underwater gliders based on a modified kinematic model combined with acoustic measurements from a single beacon. The adopted acoustic measurement is based on the travel-time differences between adjacent signals from the beacon, multiplied by the sound speed in water. Therefore, there is no need to synchronize the hydrophone time and beacon time. The depth-averaged current velocity is introduced into the modified kinematic model to generate a state estimate of underwater gliders in a flow field. An EKF method is preformed based on the estimated state and the calculated range differences. To improve the EKF estimate, the RTS smoothing algorithm is adopted for each gliding cycle after the inflection point by introducing subsequent measurements.

From simulation results, the EKF method can estimate the glider positions better than the motion model. However, the estimated glider resurfacing point is still 630 m from the actual resurfacing position. The proposed RTS-EKF method can further reduce the gap between the estimated trajectory and the true trajectory from $211.4$ m to $114.1$ m at the inflection point compared with the EKF result. When compared with the motion model estimate, the performance of the RTS-EKF estimate is improved by $46\%$ at the inflection point. The influence of the error of the adopted depth-averaged current velocity is discussed, which further illustrates the superior performance of the proposed RTS-EKF method.