Range-Only Single-Beacon Tracking of Underwater Targets From an Autonomous Vehicle: From Theory to Practice

Underwater localization is one of the main problems that must be addressed in subsea exploration, where no global positioning system (GPS) is available. In addition to the traditional underwater localization systems, such as long base line (LBL), new methods have been developed to increase the navigation performance and flexibility and to reduce the deployment costs. For example, range-only and single-beacon (ROSB) is based on an autonomous vehicle that localizes and tracks different underwater targets using slant range measurements carried out with acoustic modems. This paper presents different strategies to improve ROSB tracking methods. The ROSB target tracking method can be seen as a hidden Markov model (HMM) problem. Using Bayes’ rule, the probability distribution function of the HMM states can be solved by using different filtering methods. Here, we present and compare different methods under different scenarios, both evaluated in simulations and field tests. The main mathematical notation and performance of each algorithm are presented, where best practice has been derived. From a methodological point of view, this paper advanced the understanding of accuracy that can be achieved by using the ROSB target tracking methods with autonomous underwater vehicles.


I. INTRODUCTION
U NDERWATER localisation is one of the main problems which has to be addressed in ocean exploration, where no GPS is available due to the high attenuation that electromagnetic waves suffer in seawater [1]. Therefore, most underwater positioning systems have to be conducted with acoustic communications, despite the channel limitations, such as frequency dependent attenuation, Doppler spread and multipath propagation [2].
One of the first acoustic underwater localisation methods was the LBL, which appears in the 1960s and 1970s [3]. Since then, different alternatives have been developed, such as the Short Base Line (SBL), the Ultra-Short BaseLine (USBL), or the GPS Intelligent Buoys (GIB) systems [4] and [5]. In addition to these methods, new techniques are being developed based on the Range-Only and Single-Beacon (ROSB). For example, the moving long baseline [6] or the cooperative range-only tracking methods [7], which use the improvements in autonomous vehicles' performance, and their capabilities to work in more complex scenarios (e.g. [8]).
In general, the ROSB methods are based on an autonomous vehicle which is used as a tracker (or observer). This vehicle conducts a set of manoeuvres in order to track (or localise) some target(s). In this manoeuvre, the vehicle periodically performs new slant range measurements using the Time Of Flight (TOF) of exchanged messages between the tracker and the target (e.g. [6]), whereas the LBL method uses the Time Difference Of Arrival (TDOA) between different well localized and synchronized transponders deployed previously on the seafloor (e.g. [9] [10] [11]).
The interest in ROSB has increased in recent years as a consequence of the necessity to reduce localisation costs (e.g. transponders' deployment and clocks' synchronization) [12] [13], or to increase the vehicles' autonomy capabilities in complex missions [8] [14].
The ROSB target tracking methods can be studied from different points of view. For example, the tracker's optimal path, where [13] and [15] developed a complete analytical study for an optimal sensor placement in an underwater target localisation scenario. In [6] the optimal path for static target localisation was not only demonstrated analytically and using simulations, but also in real field tests. In both cases a Least Square (LS) and a Maximum Likelihood (ML) estimator where used, and a circumference path was determined as optimal. Later, in [7], the authors derived the same idea but for cooperative range-based underwater moving target localisation. In that case an Extended Kalman Filter (EKF) estimator was used.
On the other hand, if system observability is not taken into account, the algorithms to estimate the target's position have to solve complex multimodal probability density functions, yielding more complex algorithms. For example, in [16] the authors studied and compared the Particle Filter (PF) and the batch-MAP estimators. However, they considered an indoor scenario where a robot with movement constraints tracked a second one.
The ROSB target tracking method can also be seen as a HMM problem. Generally, the HMM is defined as a sequence of states, known as a Markov chain, and a set of observations for each state [17]. Using Bayes' rule the probability distribution function of the HMM states can be derived given a set of observations z ∈ R m , and therefore, the current state x ∈ R 2n can be estimated. Where m indicates the number of observations carried out, and n can be either 2 or 3, which is the space dimension of the problem. And p(x k |z) is the posterior probability distribution, expressed also as p(x k |z :k ); :k subscript denotes all observations up to k. The p(x k−1 ) is the prior probability distribution expressed also as p(x k |z :k−1 ). And finally, p(z) is the total probability of z [18], expressed also as x k p(z|x k )p(x k−1 )dx k , which is used as a normalized factor. However, to compute the predicted state x k , the total probability p(z) can be ignored, which yields in the optimal solution of the following maximization problem In prediction theory and filtering, the posterior distribution can be computed recursively from the prior distribution using a prediction step p(x k |z :k−1 ) and an update step p(x k |z :k ).
In general, the existing filtering methodologies compute either the predictions with respect to the conditional probability distribution p(x k |z :k ), such as PF, or with respect to the probability joint distribution p(x k , z k |z :k−1 ), such as EKF, see [19] and the references therein. One of the differences between these methods is the computational cost. Whereas the computational cost of the first methodology increases exponentially with the state dimension, the second one increases linearly with the state dimensions, therefore in areas with either high state dimensions or computational restrictions, this performance should be taken into consideration.
Hereinafter, the following considerations and parameters will be considered in all filtering methodologies which have been studied. To simplify the notation, a surface or underwater 2D scenario is used, where the tracker conducts manoeuvres on the sea surface to predict the target's position. This is a common procedure in cases where the target's depth can be determined with high accuracy using cheap devices (e.g. used in GPS Intelligent Buoys [20,Chapter 3]), and therefore, a 3D scenario can be projected into a 2D plane. Consequently, the state vector used for both tracker and target is defined as where x and y are the positions in the 2D plane, andẋ andẏ are their associated velocities. Finally, the observation measurement vector is defined as where m denotes the number of observations conducted. In ROSB methods, those are the ranges between the tracker and the target, which will be computed using the slant range measured by acoustic modems and the target's depth provided by a pressure sensor during the exchange message procedure conducted to measure the range between both devices or by the prior knowledge of the target's depth. According to the premises above, the remaining part of this work is structured as follows. In Section II different ROSB methods are presented. In Section III a brief introduction to using system observability to derive the optimal tracker's path is discussed. A set of simulations to characterize and compare different ROSB methods are presented in Section IV. Different field tests are shown in Section V. Finally, discussion and conclusions are addressed in Sections VI and VII, respectively.

II. RANGE-ONLY AND SINGLE-BEACON METHODS
In [6] we studied the optimal path shape which should be performed by an autonomous vehicle to increase the accuracy of a static target localisation using a ROSB method. The LS and the ML estimators were compared to Cramér-Rao Bound (CRB) and different field tests. These two methods are commonly used when no straightforward solution is possible [21]. For example, for either non-linear, non-smooth, or overdetermined systems (when m > n + 1).
However, when the target to be localised is not static, but moving, and active tracking is desired, the LS and ML estimators are not suitable. These dynamic scenarios are typically modelled in a state-space representation of HMM, where the next state only depends on the current state, and the current measurement depends only on the current state.
In this paper, different filters and methods have been studied and compared, presenting main aspects and some implementation. These methods are: The main aspects of these filters are presented below with a description for their implementation. We have conducted a study of the performance of different ROSB methods through simulations, and validated them through field tests. This extended study has been carried out focusing on performance comparison among different algorithms (EKF, UKF, MAP, and PF), specifically designed for 3 typical underwater scenarios: localising a static target, tracking a dynamic target, and multi-target tracking. For a methodology point of view, this work advanced the understanding of accuracy that can be achievable using ROSB localisation methods and an autonomous vehicle.
The notation employed to develop these algorithms is summarized in Table 1. List of symbols State vector Target state vector estimated at time k q ∈ R n Target position p ∈ R n Tracker position z ∈ R m Vector of ranges F ∈ R n×n State transition matrix Q ∈ R n×n Process noise matrix R ∈ R n×n Range measurement error covariance matrix Dimension of estimation problem m ∈ N Number of measurements

A. EXTENDED KALMAN FILTER
The EKF is the classical inference method for non-linear dynamic systems, which is based on the linearisation of the state and measurement equations along the trajectory [22] and [23]. This deterministic and parametric method estimates the target position based on the probability joint distribution as follows. First of all, the state vector of the target at time-step k is defined by x k = [x T kẋT k y T kẏT k ] T . Then, assuming a constant target velocity, which is a general consideration, the motion model of the target is where F is the state transition matrix, and Q is the process noise, which has variance σ 2 v . Both are related to time-step ∆t, and are described as On the other hand, the measurement model used at timestep k can be described by where q k ∈ R 2 and p k ∈ R 2 are the target and observer positions respectively in a 2D scenario, and w k ∼ N (0, σ 2 wk ) is a zero-mean Gaussian noise, leading to a covariance matrix equal to R = diag[σ 2 wk ]. Finally, the Jacobian matrix of h(x k ) is computed as Algorithm 1 has been designed to track an underwater target using the ROSB method and the EKF, which has been derived using the equations explained above, where the target state estimationx k and its associated covariance P k at each step k are given.

B. UNSCENTED KALMAN FILTER
The UKF was proposed in [24] as a derivative-free alternative to the EKF. Whereas the EKF's linearisation process incorporates inherent flaws (i.e. the expressions are approximated using a firs-order Taylor series), the UKF addresses them by utilizing a deterministic sampling strategy [25], where essentially, a set of points are propagated through the true nonlinearity, without approximation. i.e. the unscented transformation uses a set of appropriately chosen weighted points to parameterise the means and covariances of probability VOLUME 4, 2016 Input: ∆t, z i , New_range Output: Next target state estimationx k if _Init_ then Initialize: F, R, Q, P 0 ,x 0 end Predict step: distributions. These points, called sigma points χ, are propagated through the system using the state transition matrix F presented in (6).
Different methods can be used to choose the sigma points (e.g. [26] and [27]). Here, the method presented in [27] has been used, where χ ∈ R (n x 2n+1) is defined as where γ = √ n + λ, with λ = α 2 (n + k) − n. These constant values are usually set as follows: α is chosen between 1e −4 and 1, which determines the spread of the sigma points; k is set to 0 for state estimation; and β = 2 for Gaussian distributions.
Finally, these sigma points are weighted as follows, The equations (5)-(9) presented in the previous section are also used in the UKF. Following the notation of [24] and [27] the UKF for ROSB tracking is described in Algorithm 2.

C. MAXIMUM A POSTERIORI ESTIMATION
The MAP estimation [17] is a well-known method for target tracking problems. Although, the MAP estimator had long been considered to be too computationally intensive for realtime applications [28], it is becoming more commonly used thanks to processor improvements (e.g. [16] and the references therein).
Input: ∆t, z i , New_range Output: Next target state estimationx k if _Init_ then Initialize: F, R, Q, P 0 ,x 0 end Predict step: The EKF addresses the non-linear estimation problems by applying linearization methods, which introduce inherent errors. While UKF has been developed as an alternative strategy to address these errors, it only refines the current state, being unable to refine past linearized points (see the previous Section). In contrast, the MAP estimator computes the estimations of all states at all time steps, by using all available measurements.
The main equations of the MAP estimator are described below (adapted from [17] and [16]). Firstly, as mentioned, all available information is used to estimate the entire target trajectory by stacking all states in the time interval [0, . . . , k], Then, the entire state vector is estimated by maximising the posterior probability density function as follows where p(x 0 ) ∼ N (x 0 , P 0 ) is the prior distribution. By applying Bayes' rule, and assuming a Gaussian and independent noise in both measurement and state functions, plus using the target motion model (5) and the range measurement model (8) explained above, (15) can be rewritten as Using the monotonicity of the negative logarithm, the maximisation of (16) is equivalent to the minimisation of the following cost function Due to the non-linearity of the measurement model (8), there is no straightforward solution. A standard approach for its optimization is to employ iterative algorithms, which can find the solution from an initial estimationx 0 based on the recursion ofx where the parameter δx i 0:k is a correction factor, which indicates the step size and its direction. Employing the Newton-Raphson iterative minimisation method [29], as in previous works (e.g. [20]), and following [16], δx i 0:k can be found solving the linear system where A and b are the Jacobian (∇) and Hessian (∇ 2 ) of the cost function (17) with respect to all stacked states (14), evaluated at the latest state estimation, which can be obtained as and where Π, H, and F are used to adjust the dimension of a single state estimation to the entire stacked state, which have the following structure Therefore, the MAP algorithm can be formulated using these equations as shown in Algorithm 3.
T Refine all the states using Newton-Raphson iterative minimisation algorithm: while i <max or ∇c(x 0:k ) >min do Find b i and A i using (20) and (21)  Finally, a marginalization method can be used to reduce the computational cost of stacking all the states, which at a certain point can be computationally intractable. Different marginalization methods have been developed to discard old states, which are not affected significantly by a new measurement available at the current target position. For example, in [16] the Schur complement is used. However, as will be shown, a simple sliding window can also be applied with good results, where at each time-step k the state vector is updated with a new state, while the oldest one is discarded. As a result, the stacked state vector always has the same size, and therefore, the computational cost does not change. VOLUME 4, 2016

D. PARTICLE FILTER
Despite the benefits of the above algorithms, the EKF, the UKF or the MAP all have difficulties in tracking multimodal probability density functions, which is a usual problem in ROSB tracking methods [30]. Only a few estimators are specifically designed to treat multi-modal distributions. Nowadays, the PF is one of the most commonly used [31] and [32].
The PF solves, in a non-parametric way, the probability distribution problem of the HMM using the Bayes' rule (1) with the recursion of and where a bunch of particles x ∈ R 2n are spread on a 2D area, and are used to represent different possible states. Equation (25) represents the prediction step, which uses the motion model presented in (5) to move each particle with some random noise. In this case, the mean of all these particles represents the prior probability distribution. Then, using (26), each particle is weighted with a likelihood ratio based on the measurement probability function which calculates the probability of the state x n k for one dimensional Gaussian function with mean equal to the distance between the observer and the particle h(x n k ), which is the measurement model described in (8), and variance equal to σ 2 wk . In this case, the index n ∈ {0, . . . , N } indicates the particle number up to N .
Finally, all the particles are resampled according to their weight in order to obtain the posterior probability distribution and to estimate the target's position. Different resampling methods have been developed [33], where the Systematic method offers a good performance in terms of computational complexity and resampling quality. However, [34], demonstrated that other methods, such as the Compound strategy, have better performance under fast target manoeuvre circumstances.
The Compound method consists of a twofold strategy: a standard Systematic resampling method for (N − ) particles; and a Random resampling method for the last ( ) particles, which are dropped randomly inside a circular area around the previous target position that has been estimated asx k−1 . This strategy maintains particles near the target in all directions, improving the PF's time response in front of unexpected target position variations. Moreover, it maintains the particles' spatial variability, which helps to reduce the common degeneracy problem in the PF.
Using all these considerations, Algorithm 4 can be used to track underwater targets using the PF.
The state vector for each particle and its weight associated are also initialised: end Algorithm 4: PF for Range-Only and Single-Beacon target tracking.

III. OPTIMAL PATH
One of the first problems to solve in underwater target tracking is to determine the path that should be followed by the observer in order to increase the accuracy of the target estimations. The ROSB methods suffers from the multimodal distribution estimation, which is difficult to solve using standard algorithms such as the EKF. Different solutions have been found to solve this problem. For example, EKF or MAP filters have been used in parallel, where each filter tracks one possible trajectory, and a cost function is derived to find the most probable target path. One example of this is the RP-EKF [35], where each EKF uses a different initial range estimation to track a target using bearings-only, and the bank-MAP in [16], where the authors used different MAP estimators for Range-Only target tracking.
Other authors (e.g. [13], [6], [36], and [37]) have solved the multimodal problem from the system's observability perspective: i.e. by driving specific paths to maximise the amount of information or quality of the measurements conducted, affecting the accuracy of the estimated target position. Using an optimum path, the multimodal problem can be avoided, and consequently, the tracking algorithms can compute the correct target position. These studies, which are based on the Fisher Information Matrix [38], determined that a circular trajectory centred over the target maximises the system's observability. In [6], the optimal circumference for target tracking using a surface vehicle was derived not only analytically, but with field tests. Whereas the ideal circumference is one with a radius as large as possible, in real scenarios the maximum radius is typically a few hundred metres due to time constraints. For example, for a typical velocity of 1 m/s, an boat will need more than 50 min to conducts one trajectory with 500 m of radius, and moreover, the power consumption and battery limitations should also be taken into consideration. This circumference can be written as follows where p * k is the optimal position of the tracker, q k is the target's position, where m is the number of range measurements per circle, and d * k is the optimal circumference's radius.
On the other hand, in [7] the optimal radius was estimated when the tracker's velocity was taken into account according to the larger magnitude of where α k is the tracker's displacement at time k, defined by the time elapsed ∆T and the tracker's velocityv as α k = ∆Tv k , and θ k is the angle between the vectors (q k − p * k−1 ) and g(β k ).
The maximum radius in (29) has two main boundaries: the tracker's maximum velocityv; and the number of range measurements m to be conducted for each circumference, both implicitly defined in α k and θ k . As the tracker's velocity is constrained by the specifications of the vehicle being used, the only parameter that can be adjusted is the number of measurements. If a small number of measurements is used, a small circular radius will be obtained, otherwise the number of measurements must be increased. However, another option is to introduce additional time steps M ∈ {1, ..., M } between measurements to increase the circumference's radius while respecting the vehicle's maximum velocity. In this case the circumference is defined by β k = 2(M m) −1 πk.
Therefore, considering all these factors, a circular tracker path with a constant radius of 100 m (∼10 min per circle) and range measurements every 40 s have been used in our study, in order to evaluate the performance of each ROSB method.

IV. SIMULATIONS
A set of simulations have been conducted for the following scenarios: (a) Static Target Tracking; (b) Mobile Target Tracking; and (c) Multi-Target Tracking. These simulations have been performed in order to characterize the performance of each method described above, and tested them under different noise levels.
In general, Gaussian noise with zero mean and unit variance during the process of measurement is typically assumed [19] in analytical developments, whereas it has been observed that in some cases this is not accurate [6]. A Gaussian noise with non-zero mean and non-unit variance errors, w ∼ N ( w , σ 2 w ), introduces a systematic error and a random uncertainty respectively in field measurements. Moreover, some measurements can differ substantially from the true range (i.e. outliers), with a potentially strong influence on the estimations. During simulations, different outlier measurements were randomly introduced by multiplying the real range by four. The total number of outlier measurements did not exceed the 1% of the total number of measurements. Therefore, each scenario has been studied using Monte Carlo Simulation (MCS) methods [39], where different noise levels have been added to each range measurement to evaluate the tracking method's robustness in the face of them, and to obtain a more realistic simulation.
To characterise the tracking filters' performance, the step response criteria has been used. The step response concept is generally used in control system analysis to characterise the time evolution of a dynamic system [40,. It is known that the system's response has two components: transient response and steady state response. The transient response is present in the short period of time immediately after the system is turned on or a change is conducted on the input control. If the system is asymptotically stable, then the transient disappears and the system is determined by its steady state component only. Under this assumption, the ROSB methods presented in this paper can be characterised by • Settling Time (T S ): Time required to reach and stay below a threshold error • Recovery Time (T R ): Time required to reach and stay below a threshold error after a step response • Steady State Error (ε SS ): The error between the real target position and its estimation in the limit as time goes to infinity • Root Mean Square Error (RMSE): The X-Y error between the true target position and its estimation com- Finally, with all these considerations, the main parameters used to conduct the simulations are described below: • The target localisation algorithms explained in the previous section were first tested in a static scenario. This scenario is used as an initial test to evaluate the performance of the ROSB methods. Moreover, we were able to compare them against standard target localisation methods such as the LS [6], which is a good estimator to localise targets in static scenarios. Fig. 1 shows the simulated RMSE average value and its Standard Deviation (STD) after 100 MCS iterations, showing the filters' time response. This test was conducted with a range noise equal to w = 1% and σ w = 4 m. The fastest algorithm to reach an RMSE lower than 15 m is the LS (2.0 min), which is also the algorithm with the lowest ε SS (0.8 m). On the other hand, the EKF provided the worst performance with an ε SS of 9.5 m, and a T S of 30 min, i.e. it did not estimate the target's position with very high accuracy.  Table 2 shows the simulated filters' performance under other range noise parameters. The LS method is the best one in many scenarios. However, the most robust filter in the face of outliers is the PF, which has an average T S = 11 min and ε SS = 8.8 m. In contrast, the inclusion of a systematic error was not relevant to the filters' performance. It is known that concentric circumferences around a target constitute the most robust path in the face of systematic errors, as explained in [6].

B. MOBILE TARGET TRACKING
In this section we discuss the results of simulations, where a mobile target with constant velocity of 0.2 m/s plus a 90°r ight turn after the 100 th step (i.e. at 35 min from starting) was used as the second testing scenario. Besides the T S and the ε SS , the T R was also computed as the filter's time response after the right turn, when the accuracy of the target's trajectory is lost, see Fig. 2  The result of these indicators (T S , T R , and ε SS ) after 100 MCS and for different configurations of range noise are shown in Table 3. In general, the PF algorithm out-performed all the other methods, followed by the MAP algorithm. On the other hand, the LS was unable to track the mobile target as expected. The same performance is observed in all the noise cases, and for all the indicators studied.
We now discuss the performance of the PF using different configurations, e.g. the contribution of different resampling methods, the number of particles used, and the computational time required.

1) PF: Resampling Method
The Compound (C) method has the advantage of sudden response in case of fast changes in the target's direction (see Section II-D), by combining Systematic and Random resampling algorithms. Therefore, it can be adjusted by modifying the ratio between the number of particles used for the Systematic method (N sys ) and for the Random method (N rand ) respectively.
We conducted a set of simulations using different resampling methods as: (1) Multimodal, (2) Systematic, and (3) Compound. The results obtained in the time domain are presented in Table 4. After different iterations, we observed The range error introduced at each range measurement was as follows:   The range error introduced at each range measurement was as follows: that the best ratio (r(%) = N rand /N sys * 100) for the Compound resampling method was ∼6.7 (i.e. only 6.7% of particles are resampled using the Random method).

2) PF: Number of Particles
The number of particles also has an important impact on the PF's performance. The more particles used to represent the target's position, the more accurate its estimation will be. In Table 5 the results for 1000, 3000, 6000, and 10000 particles are shown. In all these simulations, the Compound method with ratio 6.7 was used.

3) Processing Time Required
Finally, the processing time can be an important constraint and a decisive factor to choose one or another method. While processing time may not be a limiting factor in some underwater scenarios due to slow dynamic processes involved, in some cases this may not be true. For example, in centralised multi-target tracking situations, the total time required to compute all the targets' positions could increase significantly with the number of targets, and therefore, the processing time must be taken into consideration.
The algorithms' runtime performance is shown in Table 6, where clearly the PF is the most expensive method from the VOLUME 4, 2016  2) a see [33] and [34] for a detailed description of these methods b (σ 2 w = 1m) computational time point of view, whereas the MAP algorithm appears as good compromise between performance and computation time. Nonetheless, both methods are suitable for this application due to the slow dynamics in most underwater scenarios.

C. MULTI-TARGET TRACKING
Multi-target Monte Carlo simulations using both static and dynamic target were used to characterise the filters' performance, in order to determine: (1) the filter's response when the tracker is not conducting its manoeuvres directly over the target, but with some offset; and (2) the feasibility of tracking multiple targets simultaneously.
Firstly, 49 targets were spread on a grid of 7x7, each one separated 100 m from its immediate neighbour. Then, a tracker conducted circular manoeuvres over the centre of the grid with a radius of 100 m and a velocity of 1 m/s. Every 40 seconds a new range measurement was computed between the tracker and each of the 49 targets, updating an individual filter for each target in order to estimate its position. This procedure was repeated 5 times, one for each tracking algorithm.

1) Multiple Static Targets
The filter's performance in a multi-static target localisation scenario is represented through coloured maps (Fig.  3). Those maps indicate the RMSE between the true target position (black triangles) and its estimation. This test was conducted 50 times using MCS iterations with a range noise equal to w = 1% and σ w = 4 m. The average value among all RMSE is presented in Fig. 3. This style of X-Y representation is commonly used [20], which indicates the target's estimation error obtained as a function of its position with respect to the tracker's path centre. For example, the best performance achievable using the LS algorithm for static target localisation was presented in [6, Fig. 4], where the CRB was used, and then verified through real field tests. Here, not only the LS but also the MAP and PF algorithms followed a similar performance.
The simulations show that the LS, MAP, and PF have a superior performance, where the targets close to the centre were better estimated. This behaviour is due to the observability of the system (see Section III). On the other hand, both the EKF and UKF have the poorest performance, with more accurate estimated positions close to the centre of the tracker's path. This behaviour can be explained by the state's initialization, where the first tracker position is used as an initial target estimation. Therefore, the targets which are close to the tracker have the best initial estimation.
Finally, the average values of the 49 target estimations for the T S and the ε SS are shown in Fig. 3 (top left) and in Table  7. The LS algorithm exhibited the best performance with a T S = 4.2 min and a ε SS = 6 m. This result is similar to the result obtained in the previous section. However, here we also showed what could be expected in the case of tracking only one target which is not directly below the tracker's path centre.

2) Multiple Mobile Targets
The performance of each tracking method using mobile targets is presented in Fig. 4. The EKF and the UKF had difficulties to track the targets when these were not directly below the centre of the tracker's path. Obviously, the LS algorithm cannot accurately estimate their position either.
In Table 8    the PF and the MAP algorithms were able to track all the targets with an acceptable accuracy, where the PF had the best performance, also when they conduct a right turn (∼30 min after the simulation's beginning).

V. FIELD TESTS
After the study conducted using the MCS methods, different field tests were conducted in order to validate the results and conclusions derived. These tests were divided into two groups: (A) Tests carried out in the OBSEA underwater observatory (www.obsea.es) of the Universitat Politècnica de Catalunya (UPC) 1) Static test 2) Dynamic test (B) Dynamic tests conducted in Monterey Bay with the support of the Monterey Bay Aquarium Research Institute (MBARI) (www.mbari.org) During those tests, the main parameters were the same ones used in Section IV to be able to perform the appropriate comparisons (e.g. number of particles (PF) = 3000).

A. OBSEA TESTS
These tests consisted in localizing three static targets and tracking a mobile one located at the coastal cabled observatory OBSEA. In both cases, S2C-18/34 acoustic modems from the EvoLogics company were used to measure ranges between observer and targets. These modems use the Sweep-Spread Carrier (S2C) technology, delivering an excellent performance, and working at 18 -34 kHz. For these tests, a small boat was employed as observer, which computes the target position using the different methods studied. Fig. 5 shows one of the modems deployed in the seafloor used in the static test (left picture), and the drifter buoy used in the dynamic test (right picture).

1) Static test at OBSEA
The first experiment carried out at OBSEA was designed to localise three acoustic modems previously deployed. One was attached to the observatory's buoy (M3) at 5 m depth, and two other modems deployed on the seafloor (M1 and M2), near the observatory's junction box at 20 m depth. Moreover, one of the seafloor modems had USBL capabilities (M1).
The slant range between the boat and each modem is represented in Fig. 6, where only one outlier (out of more than 300 measurements) between the boat and the USBL was obtained. However, this outlier must be taken into consideration, which has an important implication in the performance of the localisation algorithms, especially for the LS method as observed in Table 2. Therefore, it has to be removed in order to obtain an accurate estimation, whereas the PF is more robust even with the outlier included. Fig. 7 shows a geographic Cartesian coordinate map, where the target position estimations using the LS, MAP and PF, and the boat path conducted are represented. The PF algorithm has been executed 10 times in order to observe the prediction's variability due to the inherent random processes involved. The estimation of the second modem (M2) using the MAP algorithm has an important error, which can be caused by the lack of observability in the system. The M2 was the farther modem from the centre of the boat path, and as observed in Section IV-C, not only the steady state error is worse, but also the settling time, which is important in such situations. Under these circumstances the MAP algorithm needs more time to obtain an accurate estimation of the target's position. Finally, the results obtained during this experiment can be compared against the results obtained through simulations as depicted in Table 9. Unfortunately, the true deployment position was not available, and therefore, the RMSE could not be computed. In this case, the error between the slant range measured with the modems and the theoretical slant range computed using the targets' prediction was used, denoted as Real ε SS (Mean). On the other hand, the STD variation of the latest 10 estimations of each filter were used to obtain an indicator of its variability, denoted as Real ε SS (STD). Therefore, the values are not equal, however, the error's proportion is the same. On the other hand, both the EKF and the UKF were not taken into consideration during the field tests because of their lower performance, especially for   Table 2) b (σ 2 w = 4 m, w = 1 %) from static multi-target scenario (Table 7) targets not centred below the tracker's path (Table 7).

2) Dynamic test at OBSEA
The second experiment carried out in the OBSEA was designed to track a dynamic target, which was a drifting buoy with an acoustic modem. The results obtained are presented in Fig. 8, where the boat path (blue dotted line), the range measurements (blue dots), the real target position (black dotted line), the PF estimation (red dots), and the MAP estimation (green dots) are represented. On the other hand, the inset graphic (in the bottom-right corner) shows the RMSE between the estimated target position and its real position. Whereas the communication with the drifter was lost around 10:10 h UTC, the boat was able to track the drifter as soon as the communication was available again. However, it has to be taken into consideration that the real target position was interpolated using its initial deployment and recovery positions, and the sea currents present during the test, where the GPS position was not available during this experiment. Therefore, the mean error computed during this test should not be taken strictly into consideration, but as an indicator of the filter's performance. For example, as demonstrated in simulations, the MAP algorithm has a Recovery Time greater than the PF algorithm, and such performance is also observed in this field test, see Fig. 9.

B. MBARI TESTS
Finally, a last test in Monterey Bay California in collaboration with the Monterey Aquarium Research Institute (MBARI) was conducted in order to validate the algorithms and observe their performance under real conditions. This test was performed using a Wave Glider (WG) from the Liquid Robotics company, which tracked a Coastal Profiling Float (CPF) [41] for more than 15 hours, Fig. 10.
The CPF is a device which spends the majority of its time static, resting on the seabed. It periodically goes to the surface to fix a GPS position, then it drifts with the sea currents a few metres until it conducts another immersion to return to the seabed. During the test, the CPF conducted three immersions to ∼60 m depth, as depicted in Fig. 11 (red line).
In order to know the CPF's position, the WG conducted circular paths around the area while periodically measuring the slant range to the CPF. An acoustic modem (ATM-900) from the Teledyne Benthos company was installed in the CPF and a Benthos DAT modem installed in the WG to measure the ranges. Both devices work at 16-21 kHz, and use a phase shift keying modulation technique. The Benthos DAT modem is a standard acoustic modem which also has USBL capabilities. Because of that, a comparative study between the target's position obtained with the USBL and the position estimated using the ROSB tracking algorithms explained above could be performed. Fig. 12 shows the path conducted by the Wave Glider and the CPF, both obtained using their own GPS. Moreover, the CPF estimated path computed using the PF (red dots) and the MAP (green dots) algorithms are also shown. In addition, Fig. 13 shows the RMSE between the real CPF position and its estimations provided by the USBL system, and the PF and the MAP algorithms. The real CPF position has been computed by using the GPS positions while on the surface. However, it should be taken into consideration that no "true" CPF position while on the seabed was available. Therefore, this can cause an increase in the average error. For example, the CPF's displacement produced by sea currents during the immersion has not been taken into consideration to compute the real CPF position.
On the other hand, the error obtained from USBL is much greater than the error obtained from the PF or MAP algorithms. In general, a USBL system has to be calibrated in advance, especially to eliminate the attitude misalignment between the acoustic transducer and the Inertial Measurement    (Table 3) b (σ 2 w = 4 m, w = 1 %) from dynamic multi-target scenario (Table 8) Unit (IMU), and also to adjust their internal clocks. However, despite the calibrations, large errors can be expected due to the sea state when a USBL instrument is installed on small platforms, such as a Wave Glider. Nonetheless, the error measured during this test is something unexpected, and therefore, indicates a poor calibration or some undetected misalignment. Finally, the inset of Fig. 13 shows a zoom of PF' and MAP' RMSE results, where a mean of 21.9 m (PF) and 22.4 m (MAP) with an STD of 0.8 (PF) and 2.3 (MAP) have been obtained, which can be compared to the results obtained through the simulations as shown in Table 10.

VI. DISCUSSION
The aim of this paper is the study and development of different algorithms to track targets with autonomous marine vehicles moving along an horizontal plane by using rangeonly methods, such approach reduces the cost, power requirements, and complexity over other methods (e.g. using a USBL system which also requires an IMU). The methods presented in this paper may improve autonomous target tracking as a key factor for maritime and industries activities. For example, in the framework of fishery management (i.e. producing ancillary data for fishery management of relevant commercial items, as the Norway lobster (Nephrops norvegicus), as well as snow crabs [42]), where multiple platforms intercommunication protocols and autonomous navigation capabilities should be developed (i.e. the acoustic tracking of emitters placed on freely moving animals, spreading out from repopulating marine reserves).
Our data may also be useful in developing autonomous networks to monitor and quantify human impacts, as described by the Marine Strategy Framework Directive of the European Commission [43]. The spatial scaling of data gathered at fixed observatories, could be complemented by the use of flexible and adaptive networks of monitors and autonomous underwater vehicles. Our data could help toward the implementation of multi-parametric coordinated monitoring.
In this work, an extended study has been carried out focusing on performance comparisons among different ROSB algorithms (LS, EKF, UKF, MAP, and PF), specifically designed for 3 typical underwater scenarios: localising a static target, tracking a dynamic target, and multi-target tracking. The MCS method provides a close comparison between the simulations and real field tests conducted. Simulations are powerful tools which allow a close study of more complex and noisy/real scenarios, compared to strictly deterministic analytical studies.
For example, Table 9 summarises the field test results conducted in the cabled observatory station OBSEA (www. obsea.es) to localise multiple targets, which have been compared with the Steady State Error (ε SS ) among 100 MCS for single target localization (with 4 m of Gaussian error and 1% of systematic error), and 100 MCS for multi-target localisation (with 4 m of Gaussian error and 1% of systematic error). The results show that the target which is in the centre of the tracker path is estimated better. Finally, the real results VOLUME 4, 2016 show a greater error than the simulations, however, the error's proportion is the same, where the LS is the best algorithms whereas the MAP is the worst method. One can observe that neither the EKF nor UKF algorithms have been considered for real field tests. Moreover, the variability of these filters (STD) also follows the same trend presented in the simulations. It should also be noted that the real position of each modem was not available, and the mean error presented has been obtained using the slant range measured with the acoustic modems, and the slant range computed using their estimated positions. Furthermore, the greatest error presented in the field tests may be due to the lower number of ranges used, since the more measurements used, the greater the accuracy can be achieved.
The dynamic target tracking test conducted in the OBSEA, where a drifter buoy was used as a target, and presented in Fig. 8, shows the performance of PF and MAP algorithms in real field tests. The inset shows the evolution of the RMSE over time, where one can clearly observe that the PF has a settling time faster than that of the MAP. This behaviour has been observed previously using MCS (e.g. Fig. 4).
Finally, the test conducted in Monterey Bay with the PF and the MAP algorithms can be compared as before with simulations, see Table 10. The field test performance is shown to be quite similar to that of the simulations. However, in this case one should take into consideration that the real position of the CPF was only obtained when it was on the surface, using a GPS. Therefore, its displacement while it conducted the immersion and before it settled down on the seafloor could not be computed. This means that the RMSE presented in Fig. 13 (inset) may have an inherent error, although the general performance was demonstrated.

VII. CONCLUSION
This work shows the performance of different algorithms under different scenarios with the objective of tracking underwater target by using autonomous vehicles. The main mathematical notation of each algorithm, and their performance under simulations and field tests have been conducted, and the best practice has been derived. From a methodological point of view, this work advances the understanding of accuracy that can be achieved by using ROSB target tracking methods with autonomous vehicles.
The algorithms considered in this study are LS, EKF, UKF, MAP, and PF. All these algorithms have been compared with each other. Simulations and experimental results suggest that that an accuracy of a few metres can be achieved using the PF, which we demonstrated to be the fastest and the most accurate algorithm with respect to other studied approaches to estimate an underwater target position especially when this target is moving. For example, in a simulated dynamic scenario with a quasi ideal noise measurement of 1m, the PF achieves a settling time equal to 1.7 min, a recovery time equal to 5.8 min, and a steady state error of 1 m, but it also has more accurate values than the other algorithms in noisier cases.
IVAN MASMITJA MS in electronic engineering at the UPC. Currently, he is conducting his PhD in underwater target tracking using acoustic-rangeonly methods and moving platforms at the UPC. His field of expertise is underwater technology, involved in both underwater cabled observatories and autonomous vehicles. He has been working at SARTI-MAR research group for the UPC since 2011, participating in different national and international projects, and collaborating with different institutions: the MBARI, the Institut Supérieur de l'Électronique et du Numérique (ISEN), and the Institute of Marine Sciences (ICM-CSIC).
SPARTACUS GOMARIZ received the bachelor's, M.Sc., and Ph.D. degrees in telecommunication engineering from the UPC, Barcelona, Spain, in 1990, 1995, and 2003, respectively. He is currently an Associate Professor for the Department of Electronic Engineering, the UPC, and a member of the research group remote acquisition systems and data processing (SARTI He has worked on various platforms, including mooring controllers, benthic instruments, ASVs, and several AUVs and their associated payloads. Brian currently works on development of the Tethys AUV-a long-range, upperwater-column AUV designed primarily for biological sensing. Apart from development, Brian also takes part in mission planning and payload integration for ongoing collaborative field programs and engineering tests. VOLUME 4, 2016 TOM O'REILLY is a software engineer at MBARI. His interests include software and applications for autonomous underwater vehicles (AUVs) and autonomous surface vehicles, including navigation and tracking, adaptive sampling, and development of small low-power cytometers and microscopes for use by microbial oceanographers. He is also interested in standard protocols for marine sensor interfaces and networks, and helped to standardize the OGC PUCK interface originally developed at MBARI. Tom has been working with science colleagues to use Android devices for "crowd-sourcing" and data acquisition applications. JACOPO AGUZZI is a marine ecologists working with fisheries that is to date broadly interested in the monitoring of marine environments and the status of the services they provide by autonomous robotic technologies. He dedicated his research to the understanding of how biodiversity and ecosystem functioning are affected by behavioural rhythms of megafauna, which bias our sampling outcomes due to massive populational movements into the benthopelagic 3D scenario. Since the past decade, he is conducting an extensive monitoring activity on species behavioural rhythms and resulting community dynamism by cabled seafloor video-observatories and their docked mobile platforms, from the coastal areas to the deep-sea. As a result of that, he progressively got more interested in robotic technological solutions allowing the establishment of autonomous monitoring networks of cabled observatories and mobile platforms as well as in automated video-imaging for biological data extra ction (i.e. species tracking and classification) and storage in multidisciplinary data banks.