Mobile robotic platforms for the acoustic tracking of deep-sea demersal fishery resources

Science Robotics Manuscript Template Page 1 of 47 Mobile robotic platforms for the acoustic tracking of deep-sea 1 demersal fishery resources * 2 3 I. Masmitja1*, J. Navarro2, S. Gomariz1, J. Aguzzi2,3, B. Kieft4, T. O’Reilly4, K. Katija4, P.J. 4 Bouvet5, C. Fannjiang6, M. Vigo2, P. Puig2, A. Alcocer7, G. Vallicrosa8, N. Palomeras8, M. 5 Carreras8, J. Del-Rio1, J.B. Company2 6 7 1 SARTI Research Group, Electronics Department, Universitat Politècnica de Catalunya, 8 Barcelona, Spain. 9 2 Institut de Ciències del Mar (ICM-CSIC), Barcelona, Spain. 10 3 Stazione Zoologica Anton Dohrn, Naples, Italy. 11 4 Research and Development, Monterey Bay Aquarium Research Institute, Moss Landing, U.S.A. 12 5 L@BISEN, ISEN Brest Yncréa Ouest Brest, France. 13 6 Department of Electrical Engineering and Computer Sciences, UC Berkeley, Berkeley, U.S.A. 14 7 Department of mechanical, electronics and chemical engineering, Oslo Metropolitan University, 15 Oslo, Norway. 16 8 Computer Vision and Robotics Institute (VICOROB), Universitat de Girona, Girona, Spain. 17 * Corresponding author. Email: [email protected] 18 19 20 21 "This manuscript has been accepted for publication in Science Robotics. This version has not undergone final editing. Please refer to the complete version of record at robotics.sciencemag.org/content/5/48/eabc3701. The manuscript may not be reproduced or used in any manner that does not fall within the fair use provisions of the Copyright Act without the prior, written permission of AAAS." Science Robotics Manuscript Template Page 2 of 47 Abstract 22 Knowing the displacement capacity and mobility patterns of industrially exploited (i.e. fished) 23 marine resources is pivotal to establish effective conservation management strategies in a 24 progressively anthropized ocean. To establish the sizes and adequate locations of marine 25 protected areas within the framework of large international societal programs (e.g. European 26 Community H2020, as part of the Blue Growth economic strategy), accurate behavioral 27 information of deep-sea fished ecosystems is necessary but currently scarce and poorly accessible 28 to high-frequency and prolonged data collection. A breakthrough in the autonomous capability of 29 mobile platforms to deliver data on animal behavior beyond traditional fixed platform capabilities 30 (e.g. cabled observatories or acoustic long-baseline systems) is overcoming these limitations. 31 Here, we present useful example of that potential in relation to the implementation of autonomous 32 underwater vehicles (AUVs) and remotely operated vehicles (ROVs) as an aid for acoustic long-33 baseline localization systems for autonomous tracking of Norway lobster (Nephrops norvegicus), 34 one of the key living resources exploited in European waters. We reported the outcomes of that 35 monitoring in combination with seafloor moored acoustic receivers to detect and track the 36 movements of 33 tagged individuals at 400 m depth over more than three months. We identified 37 best procedures to localize both the acoustic receivers and the tagged-lobsters, based on cutting-38 edge algorithms designed for off-the-self acoustic tags identification. These procedures represent 39 an important step forward for prolonged, in situ monitoring of deep-sea benthic animal behavior 40 at meter spatial scales. 41 42 Summary 43 Mobile robots with different degrees of platform operability are a key element to improve and 44 extend the traditional acoustic tracking methods to study the spatiotemporal behavior of deep-sea 45 fishery resources. 46 Science Robotics Manuscript Template Page 3 of 47 47 Introduction 48 The marine benthic realm is progressively becoming wired with cabled infrastructures in an 49 attempt to transform strategic or protected areas (i.e. those of commercial or ecological value) 50 into robotized laboratories with permanent monitoring functions (1, 2). At the same time, other 51 relevant oceanic-networks are being established worldwide, seeking to track the large-scale 52 pelagic movements of species over large geographic areas and durations by animal-borne data-53 loggers (3–6). Data like these provide essential behavioral information for applying new cutting-54 edge conservation policies (7). 55 Large marine megafauna (e.g. cetaceans, dolphins, elasmobranchs or sea turtles), which rise to 56 the sea surface habitually, allows the use of the data-loggers with global positioning system (GPS) 57 and remote communication (e.g. Argos satellite network) to determine the duration and 58 trajectories of those movements (8). However, never-surface emerging benthic and pelagic 59 species cannot be tracked using such a methodology since electromagnetic waves suffer the 60 drawbacks of high attenuation in seawater medium (9). 61 For those non-emerging benthic species, acoustic positioning methods from fixed platforms 62 can be used alongside acoustic tags sensors deployed on animals (10), being tracked using long 63 base-line (LBL) triangulation techniques (11). However, deploying benthic anchored receivers 64 may increase operation complexity (e.g. in terms of spatial precision) and economic costs (12), 65 with minimal flexibility (e.g. single location). Moreover, it is necessary to use specific tags in 66 order to maintain synchronization between each receiver into the listening network, which may 67 increase the complexity in data post-processing (13). While this technology has proven useful for 68 behavioral tracking in shallow water scenarios (14), its performance has not yet been fully 69 examined in the deep-sea. Only a few efforts have been conducted to follow populations 70 Science Robotics Manuscript Template Page 4 of 47 movements over kilometer scales (15) with acoustic receivers mounted on moored of curtain or 71 gate typologies (16). 72 A complementary strategy to the use of moored devices is to mount acoustic receivers on 73 autonomous underwater vehicles (AUV), which is used as a virtual LBL, measuring the distance 74 range with the target by acoustic modems (17, 18). Differently from acoustic tags, these modems 75 have bi-directional communications capabilities, and therefore, the time of flight (TOF) and slant 76 range of an acoustic signal can be measured knowing the sound’s velocity. Finally, triangulation 77 localization techniques are applied to estimate the position of tagged-individuals with different 78 algorithms (19). In addition, the bearing information estimated by ultra-short base-line (USBL) 79 systems can also be used, which increases the overall speed response (20, 21). Nevertheless, the 80 investigations are again limited to large animals due to the size of electronic tags (22). Other 81 authors use bearing-only techniques in order to avoid the use of acoustic modems and overcome 82 the size limitations (23, 24), where an AUV borne hydrophones’ array is used to track acoustic 83 tags. Unfortunately, a significant localization uncertainty is produced by the too-close positioning 84 of hydrophones, which often requires larger separation that is not achievable on AUVs (25). 85 However, marine robots have in recent years been used for the tracking of marine species. For 86 example, AUVs equipped with a single hydrophone were used to track fishes with different error 87 ranges and procedures (e.g. SYNAPS and SPLWCA (14, 26, 27)). Some of the studies were 88 conducted in combination with seafloor moored receivers, allowing records of the presence of 89 tagged animals within the area of detection, but with high uncertainty in their position (28–31). 90 Other authors used custom transponders attached to large marine species to increase the efficiency 91 of vehicle tracking capabilities (20, 32), but this approach is impractical for small marine species. 92 Despite this interest, to our best knowledge, no previous study has addressed the tracking 93 methodologies’ performance by using static receivers and underwater vehicles, including testing 94 its accuracy and capabilities in deep waters. 95 Science Robotics Manuscript Template Page 5 of 47 Here, we describe a new procedure for the multi tracking of one of the most important fishery 96 resource in Europe, the Norway lobster (Nephrops norvegicus Linnaeus, 1758) (33), using a set of 97 moored seabed receivers along with a remotely operated vehicle (ROV) and an AUV (Fig. 1). We 98 developed a new area-only target tracking (AOTT) method to achieve active tracking of 99 instrumented individuals, which only uses the detection pings of acoustic tags. Moreover, time 100 difference of arrival (TDOA) algorithms have been adapted and tested to study their accuracy in 101 variable operational scenarios. Specific objectives were: (i) TDOA algorithms performance 102 comparison through Monte Carlo (MC) simulations; (ii) new AOTT algorithm capabilities 103 presentation, where both simulations and field tests have been conducted; and (iii) the results of a 104 3-month campaign using static receivers (i.e. mooring lines with acoustic receivers) and 105 underwater vehicles (an ROV and an AUV), where both TDOA and AOTT algorithms have been 106 used. The tracking potential of this combined mobile and moored technology was tested in a 107 deep-sea, no-take fishing zone under restoration, to show how a long-lasting acoustic-based 108 deployment can provide new behavioral data that can inform the establishment and spatial extent 109 of conservation areas. 110 111 Results 112 113 TDOA tracking algorithms performance 114 The performance of the different TDOA algorithms tested is presented in Fig. 2, where a set of 115 MC simulations has been conducted. We used 4 receivers to localize a target on two-dimensional 116 (2D) scenario, since the depth of the targets under study was known and constant as the species is 117 benthic, with no swimming capability. These simulations are important to demonstrate the 118 capabilities to track benthic tagged animals and to set the appropriate configuration (e.g. number 119 of receivers, receivers’ positions, or acoustic tag transmission period). 120 Science Robotics Manuscript Template Page 6 of 47 The Cramér-Rao bound (CRB) representation is presented in Fig. 2A, where 4 receivers with 121 200 m of baseline distance and a time error of 1 ms have been used (for other array configurations 122 see Fig. S1). The area inside of the receivers’ array showed the lowest expected measurement 123 standard deviation error (< 1 m), whereas the error increased up to 7 m at 250 m off the receivers’ 124 array center. 125 To compare the algorithms’ performance, a predefined target trajectory has been designed 126 (Fig. 2B and Movie S1), according to which the target moves at 1 ms-1 among fourth receivers 127 with a transmission period of 60 s. The root mean square error (RMSE) over the time is shown in 128 Fig. 2C, where all the algorithms were iterated 100 times with a Gaussian error with 1 ms 129 standard deviation (Fig. S2 shows the RMSE evolution with other errors). This result clearly 130 showed that the error is lower inside the receivers’ array, especially for the maximum likelihood 131 (ML) estimation algorithm. That latter registered the greatest error. Due to numerical singularities 132 around the receivers, the ML estimation failed to find the minimum of the cost function, and 133 instead, it reached the local minimum nearby the receiver position. This problem was reduced by 134 choosing a different initial estimation (i.e., closer to the real position), as explored with the 135 weighted least squares ML (WLS-ML) algorithm, where the WLS method is used to initialize the 136 ML estimation algorithm. 137 The algorithms’ RMSE over the 100 MC iterations with different noise added in the time of 138 arrival (TOA) measurements are presented in Fig. 2D. We simulated the algorithms’ performance 139 with noise standard deviation ( σ ) equal to 0.5, 1.0, and 1.5 ms. Moreover, additional tests with a 140 Gaussian TOA error with 1.5 σ = ms plus 5% of outliers were simulated to observe the 141 algorithms’ behavior when facing strong multipath scenarios. These simulations showed that the 142 particle filter (PF) had the best performance under different noise conditions, however, it had 143 more difficulties to handle scenarios with outlier measurements, whereas the WLS excelled. In 144 addition, the use of the WLS-ML combination slightly improved the algorithm’s performance 145 Science Robotics Manuscript Template Page 7 of 47 with different noise configurations. Nonetheless, this benefit was not observed in scenarios with 146 outliers. 147 Finally, the average runtime required to compute one target position is shown in Fig. 2E, were 148 the fastest algorithm was the WLS with 5 ms per iteration. In contrast, the PF required 977 ms, 149 which means an increase of more than two orders of magnitude of the computational resources. 150 151 AOTT algorithm performance 152 A set of simulations were conducted to observe the optimal parameters for the AOTT algorithm 153 (Fig. 3A) and the functions to weight the PF’s particles (Fig. 3B). For example, the results 154 showed a clear relationship between the tracker circumference radius (TCR) and the maximum 155 transmission range (MTR), where the greatest ratio, radius /TCR MTR Γ = , was equal to 0.8 (Fig. 156 3C). This means that the tracker had to conduct circumference maneuvers over the target 157 estimation position with radius less than the MTR but closer to it. Nonetheless, is hard to known a 158 priori the MTR achievable by an acoustic tag, which can be affected by different factors such as 159 the sea state or the acoustic noise. Therefore, different in situ tests should be conducted to 160 estimate its value. In our case, those tests pinpointed a maximum range less than 400 m with only 161 a 20% of successful receptions (Fig. S3). In addition, the MTR is pivotal to spread the PF’s 162 particles, and therefore, different relationships between MTR and the maximum particles range 163 (MPR) were studied, which allowed to identify the relation between the ratio 164 range /MPR MTR Γ = and the AOTT’s performance (see Fig. 3C). Moreover, the behavior of 165 missing some of the tag’s transmissions could also be observed, where the successful reception 166 (SR) over the total transmissions (TT) ratio defined by reception /SR TT Γ = is presented. Finally, 167 random particles were spread around the latest estimated target position (Compound resampling 168 Science Robotics Manuscript Template Page 8 of 47 method), which helped to increase the particles diversity, and emphasized the latest time that the 169 tag was detected, which yielded to an increase in tracking performance (see Fig. 3C). 170 The AOTT’s performance can be observed on Fig. 3D (and Movie S2), where all the 171 recommendations derived from the previous MC simulations presented above were used, which 172 showed an error of 100 m. After these simulations, a field test was conducted on June 27-28, 173 2018 (Fig. 3E) using the Monterey Bay Aquarium Research Institute (MBARI) coastal profile 174 float (CPF) as a target (Fig. 3F) and a Wave Glider as a tracker (Fig. 3G) in Monterey Bay area, 175 (CA, USA) (Fig. 3H). This test lasted more than 15 h, where the CPF conducted 3 immersions at 176 60 m depth. 177 178 Norway lobster tracking 179 The results of the four-step process to adjust the receiver clocks’ drift and offset is shown in Fig. 180 4A–F, where a resolution greater than 2 ms was obtained. Moreover, a small number of outliers 181 were detected during the post-processing (e.g. Fig. 4E), which had a random nature due to the 182 homogeneous bathymetry of the experiment zone (i.e. a quasi-flat slope ground). In addition, the 183 deployment position of the mooring lines using the oceanographic vessel’s GPS, the ROV’s 184 USBL and the positions computed using the acoustic receivers are presented in Fig. 4G. Here, a 185 great difference between the GPS’s and ROV’s positions could be observed, which pinpointed the 186 necessity of the use of underwater vehicles to know final position of the receivers (Table 1). 187 After determining the receiver localizations and calibrating their clock offsets, the tagged 188 Norway lobster positions could be tracked using the TDOA algorithms (see next section). The 189 trajectory showed by each animal can be observed in Fig. 5A (and Movie S3), where the 190 localization of the synchronization acoustic tags attached on the mooring lines, and the acoustic 191 tags attached on the 33 individuals, are shown. After the canister release in the center of the 192 receivers’ array, the individuals show a dynamic dispersion and occupation of the monitored area. 193 Science Robotics Manuscript Template Page 9 of 47 Furthermore, the accumulative distance of each individual was plotted in Fig. 5B, where we could 194 appreciate how some animal went outside of the receivers’ reception range, being therefore not 195 detectable any more, from that moment on. 196 By the use of the two underwater robots (an ROV and an AUV), we could track the presence 197 of some of those area-evading animals. The ROV conducted different lawn pattern movements on 198 the southeast of the area, covering 10 km2, and the AUV conducted a circumference path on the 199 west (see Fig. 5C) with a radius equals to 150 m. During these tests, 4 tags were localized, and 200 moreover, different images could be obtained, which will be used to study the seabed recovery in 201 the protected areas (Fig. 5D and Movie S4). 202 203 Comparison between methods 204 The algorithms studied to track acoustic tags using the TDOA information could be compared 205 together during the entire Norway lobster tracking experiment. Because the “true” position of the 206 tagged lobsters was unknown, the synchronization tags attached on the mooring lines were used. 207 In this case, the tags were not moving but static. Fig. 6A-6E show their estimated position and the 208 error covariance matrix, which are represented as error bars in Fig. 6F (and summarized in Table 209 S1). For example, the accuracy obtained to localize the lobster canister synchronization tag (i.e. 210 base station (BS) D), which was placed in the receiver array center, was similar among all the 211 algorithms (error <1 m). However, the PF had the poorest performance when it came to localize 212 the synchronization tags attached on the mooring lines. This low performance was due to the 213 nature of the PF’s particles distribution near the receivers in a TDOA topology (i.e. eccentricity of 214 the hyperbola close to 1). Moreover, we found that both PF and WLS methods showed higher 215 errors in positioning moorings Vemco acoustic receiver indicated as BS(D) and BS(E), which had 216 different configuration (smaller dead weights and VR2AR-69k receivers). Taking into 217 Science Robotics Manuscript Template Page 16 of 47 other lines and the lobster canister through the synchronization tags and the receivers, their 367 relative positions could be determined by simple trigonometry functions and rotation matrices. 368 369 TDOA algorithms 370 Target localization using TDOA is a well-known problem which has been addressed on both 371 terrestrial and underwater environments during the last decades. The TDOA has been usually 372 used when no synchronization between transmitters and receivers can be enforced, and even 373 more, if transmitters ping-time is irregular (e.g. using Vemco devices). In both cases, the TOF 374 cannot be measured or estimated, and consequently, the TDOA between different pair of receivers 375 is used. 376 In general, TDOA algorithms can be divided in two groups, the ML and least-squares (LS) 377 methods (47). Using 1n+ receivers (where { } 2,3n∈ is the space dimensionality of the problem) a 378 set of hyperbolic equations can be obtained to find the coordinates of the target. The TDOA 379 measurement between two receivers n ∈b and the target at position n ∈q can be written as 380 381 00 1 || 11 () ( ) ( | || ) (|| | || |||| ) ij ji i j t tw cc w c µ =+ − −+ − + = − −− + x ||qb qb qb qb , (1) 382 where , {0, , }ij m∈… and ij≠ , c is the sound velocity in water, and 0 t is the target 383 transmission time. Assuming a zero-mean white Gaussian error noise distribution of the TDOA 384 measurements, i.e. 2 ~ (0, )w σ  with variance 2 σ , the unknown parameter n ∈q can be 385 estimated using the ML estimation method. In this case, the density function for each () ij µ q is 386 given by 387 Science Robotics Manuscript Template Page 17 of 47 () 2 2 2 - ( ) 1 ( ) exp 2 2 ij ij f µµ σ πσ   = −   q q , (2) 388 where µ represents the measured TDOA. Given a vector of observations m ∈μ the function 389 : [0,1] n →⊂ which for any target position n ∈q yields the probability (|) pμq , is 390 referred to as the likelihood function, given by 391 ( ) ( ) ( ) ( ) 2 2 2 1 2 2 21 2 1 2 2 - ( ) 1 : ( | ) exp 2 2 11 exp- - ( ) 2 2 11 exp || - ( || | ) |2 2 k k m k m k mm k k p µµ σ πσ µµ σ πσ π = =   = = −    =   = −   ∏ ∑ R q qq q q R μ  µµ , (3) 392 where 21 || || T− M a aM a  , and R is the covariance matrix, and m I is the identity matrix of dimension 393 mm× . The ML estimator is defined as 394 ˆarg max ( ) n ∈q q = q   . (4) 395 A common practice in ML estimation is to work with the log-likelihood function. Since the 396 logarithm is a strictly increasing function, and ( ) q is strictly positive, maximizing the 397 likelihood and the log-likelihood are equivalent. Neglecting constant terms, the ML estimator can 398 be found by solving the optimization problem 399 ˆarg min ( ) n f ∈q q = q  , (5) 400 where : n f→ is given by the following cost function 401 ( ) ( ) 2 21 11 (): || - ()|| - () - () 22 T f − == R q q qR q µµ µµ µµ . (6) 402 Science Robotics Manuscript Template Page 18 of 47 In general, there is no closed form solution to the previous optimization problem. The cost 403 function is relatively complex, nonlinear and even not differentiable at some points because of the 404 square roots that defines the TDOA measurements. 405 A standard approach for its optimization is to employ Newton-Raphson iterative minimization 406 (48). In order to implement gradient and Newton descent algorithms to minimize the cost function 407 it is necessary to have expressions for its gradient ()f∇q and Hessian 2()f∇q , which are the 408 vector of its first partial derivatives and matrix of its second partial derivatives respectively. This 409 can be done resorting to Matrix Differential Calculus, see (11, 49) and the references therein. 410 Nonlinear estimation problems are also often addressed using linearized estimators, e.g., the 411 extended Kalman filter (EKF) (50). However, linearization-based filtering approach marginalize 412 all but the current state and is hence unable to refine past linearization points. In contrast, a batch 413 maximum a posteriori (MAP) estimator computes the estimates for the states at all-time steps 414 using all available measurements (51). The difference between MAP and ML estimation lies in 415 the assumption of an appropriate prior distribution of the parameters to be estimated (52). The 416 MAP estimator utilizes all available information to estimate the entire target’s trajectory which is 417 represented by stacking all states in the time interval [ ] 0, k as 418 0:k 0 1 k T TT T  = x xx x , (7) 419 where 2 Tn k qk qk qk qk xxyy  = ∈  x  is the target’s position and all the higher order time 420 derivatives (i.e. velocity or acceleration). In addition, a motion model is used, which typically 421 consider that the target moves randomly but assume that a stochastic kinematic model describing 422 its motion (e.g., constant velocity) is known. Thus, the discrete-time state propagation equation is 423 generally given by 424 11 1k kk k−− − =Φ+x xw , (8) 425 Science Robotics Manuscript Template Page 19 of 47 where 1k− w is zero-mean white Gaussian noise with covariance Q , and the state transmission 426 matrix, 1k− Φ , is given by 427 1 1 00 0100 001 0001 k t t − ∆    Φ=  ∆   . (9) 428 Then, the MAP estimator seeks to determine the entire state-space trajectory that maximizes 429 the following posterior probability density function 430 ( ) ( ) ( ) 0|0 ' 1: 0: 2 0 0|0 1 2 0|0 k 2 11 1 1'2 k 2 1 12 2 (|) 11 ˆ exp || - || 2 2|| 11 exp || - || 2 2 || 11 exp || - ( ) || 2 2 || kk n k kk n k kk m k k p π π π −− = = ∝  −    ×−    ×−   ∏ ∏ P Q R μx xx P xx Q q R Φ µµ , (10) 431 where a prior distribution equal to 0 0|0 0|0 ˆ (,) ()p= x xP has been used, and 1:k μ denotes all the 432 measurements in the time interval [ ] 1, k . Using the same procedure as in eq. (3), the cost function 433 is given by 434 0|0 ' 2 0:k 0 0|0 2 k 1 k 1 11 2 1ˆ ( ) : || - || 2 1|| - || 2 1|| - ( ) || 2 k k kk kk k f − = = − = + + ∑ ∑ P Q R x xx xx q Φ µµ . (11) 435 And finally, the solution can also be computed employing Newton-Raphson iterative 436 minimization methods, see (51) and references therein. However, this solution, heavily depends 437 Science Robotics Manuscript Template Page 20 of 47 on the quality of the initial estimate, especially if multi-modal probability density functions are 438 involved (i.e., the solution may lie on local minimum instead of the true target position). 439 To estimate multi-modal distributions, one of the most used methods is the PF (53, 54). The PF 440 solves in a non-parametric way the probability distribution problem using a set of particles, 441 2n ∈x , which are spread on the area in order to represent the true distribution. Each particle 442 represents a hypothesis of the target state. The particles are weighted and normalized based on 443 their measurement likelihood, and resampled accordingly (55, 56). 444 Another method to solve the likelihood function eq. (3) is using a closed-form LS solution. A 445 wide used closed-form method was developed by Chan and Ho (36). They give an alternative 446 solution for hyperbolic position fix by using an approximation of the ML estimation when the 447 TDOA estimation errors are small. The original set of TDOA equations are transformed into 448 another set 1 0 TnT r+  = ∈  x q  , which are linear in source position coordinates q , and adding 449 an extra variable 0 r , which is the range between the target and the reference sensor. Then, the 450 algorithm uses a two-step WLS method to estimate the target position, which is given by 451 ( ) 111TT a aa −− − =xGGGhΨΨ , (12) 452 Where 453 2 22 10 10 1 0 10 2 22 2 20 10 1 0 10 2 22 1 0 10 10 0 || || || || || || || || 1 ,, 2 || || || || T T a aa T L rr rrc r r   −−   −−   =    −−      bbb bbb G h = B RB bb b   Ψ= , (13) 454 and 0 || - || am =B qb I . 455 Further, different authors have improved this technique, for example, in (57) the WLS includes 456 a vertical plane constraint and a cone tangent plane constraint. These two constraints are derived 457 from the initial value and updated again after each iteration. 458 Science Robotics Manuscript Template Page 21 of 47 Finally, in (37) the authors developed the yet another positioning solver (YAPS) method, 459 where they used the TOA instead of the TDOA to estimate the target position. Because of that, 460 they had to also estimate the target transmission time k t . The modelling follows the state space 461 paradigm, which uses the process and observation models as in MAP estimation method. They 462 used a stochastic processes to describe the state propagation as a random walk with different 463 degrees of standard deviation for both transmission time 2 1 12 ,() k k k k bi tt t t σ − −− −− and target 464 position 0.5 1),2( k k xy Dt −∆ qq  , where xy D is the diffusivity. The YAPS method is coded as a 465 C++ file, which evaluates the joint density through the template model builder (TMB) framework. 466 That latter uses the Laplace approximation to find the unobserved random variables (e.g. x , y , 467 and t ), and the parameters (e.g. xy D ) that can be estimated using the ML principle and built-in 468 optimizer in R. Therefore, this model analysis follows a standard ML analysis of non-linear 469 mixed-effects model, suing the TMB as the computational tool to automatize the process with R 470 software. 471 Here, all these algorithms have been compared with the CRB (58), which sets the lowest 472 bound on the performance of unbiased estimators that use observations according to a certain 473 probability density function. This bound is one of the most widely used (59–61), which for a 474 TDOA target localization problem is given by 475 { } 1 ˆ C ( ov ) − Iq q , (14) 476 where I denotes the Fisher Information Matrix (FIM) defined as 477 1 () () () T ff − =∇∇Iq q R q , (15) 478 where ()f∇q is the gradient of the log likelihood function with respect to the unknown 479 parameters, which has been used to compute the target position using the ML estimation. Taking 480 Science Robotics Manuscript Template Page 22 of 47 the trace of I we obtain a new inequality, which sets a fundamental lower bound on the mean-481 square error of any unbiased estimator, given by 482 { } { } ( ) 12 ˆˆ var =E || ( ) | r ()|t µ − −≥q Iq qq (16) 483 TDOA Simulations 484 Different simulations have been conducted in order to characterize the TDOA target localization 485 algorithms explained above under different parameters and scenarios. These simulations have 486 been carried out using the MC simulation method. For all the simulations, the RMSE has been 487 computed using the median, and the 5th and 95th percentile, over 100 iterations, where different 488 TDOA Gaussian noise has been added using 0.5 σ = ms, 1 σ = ms, and 1.5 σ = ms. The 489 parameters of the scenario simulated used were: (i) tag transmission delay = 120 s, (ii) target 490 velocity = 0.2 m/s, and (iii) number of particles (for the PF algorithm) = 6000 particles. 491 Algorithms’ run-time has been obtained using a Processor Intel® Core™ i7-4760HQ CPU @ 492 2.10 GHz with 8 GB of RAM memory. 493 494 Receiver clock drift adjustment and localization 495 Four receivers have been used in this study, where each one has an internal clock which is not 496 synchronized periodically. Consequently, during the campaign they suffered from drift and 497 misalignment. This behavior introduces an error which must be fixed for twofold: (i) to be able to 498 associate independent receptions at separate receivers, corresponding to the same target and 499 emission time, and (ii) to compute the TDOA accurately. The TDOA between two receivers 500 (considering their clocks’ drift) can be modelled as 501 ( ) ( ) || |( ) || | 1 ij k k k j ik jki CC c µ = −−− + −q ||qb qb , (17) 502 Science Robotics Manuscript Template Page 23 of 47 where ik C is the clock’s misalignment of receiver i at time step k . Considering static receivers 503 and a static acoustic tag 0 q (typically localized in the center of the receivers’ array), the 504 measurement 0 () ij k µ q should be constant. However, due to the differences in the clocks’ drift 505 00 ||ij ik jk CC C= − qq this is not true, which would result in target localization errors. Therefore, here 506 we developed a procedure to adjust the drift using a four-step process: (i) using the initial points 507 and a linear regression, (ii) using a polynomial regression with all the points, (iii) using different 508 polynomial regression functions at different segments of data, and (iv) using the distance 509 difference to correct the offset. 510 The first step was used to adjust the main drift, which is necessary to associate independent 511 receptions at separate receivers. If the clock's drift is greater than the acoustic tag transmission 512 interval time, it is not possible to associate the receptions of an acoustic tag transmission at 513 different independent receivers (in long field studies, i.e. more than one month, the drift can reach 514 more than 30 seconds). Thus, only the initial points can be used. Then, the different receptions 515 can be associated and a polynomial fitted curve can be used to eliminate the main clocks' drift for 516 the entire data. In addition, the whole data was segmented into small portions (e.g., by weeks), 517 and a second polynomial fitted curve was used for a fine tune. With this procedure, the drift was 518 adjusted (i.e. the slope of the ij C , aka slope ij C ). Nonetheless, a final step to adjust the clocks’ offset 519 was still necessary. We know that the distance between the receiver pair ij have to be equal to the 520 distance between the receiver pair ji , and therefore, an offset equal to ( ) 2 offset ij ij ji C dd= − can be 521 added. Thus, each clock’s receiver adjustment is given by 522 , ,( ) , 0 N slope r offset i n ij N r i n ij r Clk C Clk C − = = + ∑ , (18) 523 Science Robotics Manuscript Template Page 24 of 47 where ,in Clk is the timestamp value n of receiver i , and N is the polynomial degree of the fitting 524 curve. 525 Once the internal clock drift was adjusted, the position of each receiver i b could be computed. 526 First, the distance among each receiver ij d were calculated using the TOF, which was known due 527 to the fact that each receiver had also a synchronization acoustic tag attached on the mooring line, 528 and therefore, the 0 t was known. Then, the i b positions were computed using these distances and 529 trigonometry. Finally, the positions were adjusted using a rotation matrix and a translation matrix 530 to obtain the final position referenced to the geographic coordinates system, where the mooring 531 anchors’ positions found by the ROV were used. 532 533 AOTT algorithm 534 The AOTT method uses a single moving receiver, and therefore, no TDOA information is 535 present. In its place, the tag’s position is estimated by using the ping detection/no-detection 536 information provided by a receiver. However, the detection of a tag’s transmission is complex due 537 to acoustic noise form platforms’ thrusters, multipath, or distance between the tag and the 538 receiver. Consequently, the AOTT algorithm attributes such as the reception ratio or maximum 539 transmission range have been studied through simulations and field tests before the Norway 540 lobsters’ field survey. 541 Given the acoustic receiver and transmitter tag used for this work, the only information that 542 can be determined is the presence or absence of acoustic tag transmissions in the area of the 543 receiver, without information about the tag’s direction or range of detection. The AOTT method 544 infers the target position by taking the area determined by the maximum reception range as the 545 only filter input (62). Two types of areas can be defined: one where the tag is detected, and one 546 where the tag is not detected. The estimation of the target’s localization can then be computed by 547 Science Robotics Manuscript Template Page 25 of 47 overlapping all of these areas, where the zone with a main coincidence is where the target should 548 be, thereby representing its probability distribution. 549 The AOTT was implemented by using a PF algorithm, where a set of particles 2n ∈x are 550 randomly spread in the area, and then, each particle is moved accordingly to a motion model eq. 551 (8), and each particle’s weight is updated for each new detection (or no-detection) until all of 552 them converge into the target position estimation. Therefore, the probability distribution function 553 can be derived using the Bayes’ rule (63) with the recursion of the prediction step 554 :1 1 1 :1 1Motion model Particles (| ) (| )( | ) k kk kk k k p pp − − −− − =∑ x xz xx x z      , (19) 555 and the update state 556 : :1 Importance weights Particles (|) (|)(| ) kk kk kk p pp − ∝ xz zx xz      , (20) 557 where m ∈z are a set of measurements. 558 The main difference between the range-only (19) and area-only target tracking algorithm based 559 on PF is how the particles are weighted. In a range-only method, the likelihood ratio based on the 560 measurement probability function is described as 561 ( ) 2 2 2 - ( ) 1exp 2 2 n kk n k w w zz W σ πσ   = −   x , (21) 562 in this case, the index {0, , }nN∈ indicates the particle number up to N . 563 Whereas in the area-only method, the measurement probability function is based on the 564 distance that each particle has between each other and the observer, where the particles which are 565 inside the area defined by the maximum range that an acoustic tag can be detected will be more 566 weighted than the particles which are outside of this area. Moreover, if an acoustic tag detection is 567 missed, the particles inside the area will be less weighted than the particles which are outside. 568 Science Robotics Manuscript Template Page 32 of 47 61. C. Radhakrishna Rao, Linear Statistical Inference and its Applications (John Wiley & 812 Sons, Inc., Hoboken, NJ, USA, 1973), Wiley Series in Probability and Statistics. 813 62. I. Masmitja, S. Gomariz, J. Del-Rio, B. Kieft, T. O’Reilly, J. Aguzzi, P.-J. Bouvet, C. 814 Fannjiang, K. Katija, in OCEANS 2019 - Marseille (IEEE, 2019), pp. 1–10. 815 63. M. S. Arulampalam, S. Maskell, N. Gordon, T. Clapp, in Bayesian Bounds for Parameter 816 Estimation and Nonlinear Filtering/Tracking (IEEE, 2009), vol. 50, pp. 723–737. 817 64. J. Zhang, D. Berleant, Envelopes around cumulative distribution functions from interval 818 parameters of standard continuous distributions. Annu. Conf. North Am. Fuzzy Inf. Process. 819 Soc. - NAFIPS. 2003-Janua, 407–412 (2003). 820 65. G. K. Karagiannidis, A. S. Lioumpas, An improved approximation for the Gaussian Q-821 function. IEEE Commun. Lett. 11, 644–646 (2007). 822 66. I. Masmitja, S. Gomariz, J. Del Rio, P. J. Bouvet, J. Aguzzi, Underwater multi-target 823 tracking with particle filters. 2018 Ocean. - MTS/IEEE Kobe Techno-Oceans, Ocean. - 824 Kobe 2018, 1–5 (2018). 825 67. F. Arrichiello, G. Antonelli, A. Pedro Aguiar, A. Pascoal, An observability metric for 826 underwater vehicle localization using range measurements. Sensors (Switzerland). 13, 827 16191–16215 (2013). 828 68. N. Crasta, D. Moreno-Salinas, B. Bayat, A. M. Pascoal, J. Aranda, Range-based 829 underwater target localization using an autonomous surface vehicle: Observability analysis. 830 2018 IEEE/ION Position, Locat. Navig. Symp. PLANS 2018 - Proc., 487–496 (2018). 831 69. J. D. Quenzer, K. A. Morgansen, Observability based control in range-only underwater 832 vehicle localization. Proc. Am. Control Conf., 4702–4707 (2014). 833 70. E. W. Griffith, K. S. P. Kumar, On the observability of nonlinear systems: I. J. Math. Anal. 834 Appl. 35, 135–147 (1971). 835 71. I. Masmitja, S. Gomariz, J. Del Rio, B. Kieft, T. O’Reilly, Range-only underwater target 836 localization: Path characterization. Ocean. 2016 MTS/IEEE Monterey, OCE 2016, 1–7 837 (2016). 838 72. D. Moreno-Salinas, A. Pascoal, J. Aranda, Optimal Sensor Placement for Acoustic 839 Underwater Target Positioning with Range-Only Measurements. IEEE J. Ocean. Eng. 41, 840 620–643 (2016). 841 73. N. Crasta, D. Moreno-Salinas, A. M. Pascoal, J. Aranda, Multiple autonomous surface 842 vehicle motion planning for cooperative range-based underwater target localization. Annu. 843 Rev. Control. 46, 326–342 (2018). 844 74. B. Ha, G. Massion, “Coastal Profiling Float Depth Control” (Moss Landing, CA, 2018), 845 (available at https://www.mbari.org/wp-content/uploads/2018/12/Ha.pdf). 846 75. G. B. Skomal, E. M. Hoyos-Padilla, A. Kukulya, R. Stokey, Subsurface observations of 847 white shark Carcharodon carcharias predatory behaviour using an autonomous underwater 848 vehicle. J. Fish Biol. 87, 1293–1312 (2015). 849 76. D. Haulsee, M. Breece, D. Miller, B. Wetherbee, D. Fox, M. Oliver, Habitat selection of a 850 coastal shark species estimated from an autonomous underwater vehicle. Mar. Ecol. Prog. 851 Ser. 528, 277–288 (2015). 852 77. D. R. Zemeckis, M. J. Dean, A. I. Deangelis, S. M. Van Parijs, W. S. Hoffman, M. F. 853 Baumgartner, L. T. Hatch, S. X. Cadrin, C. H. McGuire, R. O’Driscoll, Identifying the 854 distribution of Atlantic cod spawning using multiple fixed and glider-mounted acoustic 855 technologies. ICES J. Mar. Sci. 76, 1610–1625 (2019). 856 78. J. H. Eiler, T. M. Grothues, J. A. Dobarro, R. Shome, Tracking the movements of juvenile 857 Chinook salmon using an autonomous underwater vehicle under payload control. Appl. Sci. 858 9 (2019), doi:10.3390/app9122516. 859 79. L. Zacher, “Tracking the Alaskan Red King Crab – Post 1-8, NOAA Alaska Fisheries 860 Science Center” (2019). 861 Science Robotics Manuscript Template Page 33 of 47 80. J. S. Horne, E. O. Garton, S. M. Krone, J. S. Lewis, Analyzing animal movements using 862 Brownian bridges. Ecology. 88, 2354–2363 (2007). 863 864 Acknowledgments 865 866 Acknowledgments: We gratefully acknowledge the assistance and support of Larry Bird 867 (MBARI) and the David and Lucile Packard Foundation. We would also like to thank Nixón 868 Bahamón, José Antonio García, and Guiomar Rotllant for their help during the cruises. Ruth 869 Durán helps us with the study area map. This work has been lead and carried out by members of 870 the Tecnoterra associated unit of the Scientific Research Council through the Universitat 871 Politècnica de Catalunya, the Jaume Almera Earth Sciences Institute and the Marine Science 872 Institute (ICM-CSIC). 873 Funding: This work received financial support from different research projects of the Spanish 874 Ministerio de Economía y Competitividad (RESNEP: CTM2017-82991-C2-1-R, RESBIO: 875 TEC2017-87861-R, and SASES: RTI2018-095112-B-I00), of the Generalitat de Catalunya 876 “Sistemas de Adquisición Remota de datos y Tratamiento de la Información en el Medio Marino 877 (SARTI-MAR)” 2017 SGR 371, and NSF-IDBR (#145501 to K. Katija). J. Navarro was funded 878 by the Spanish National Program Ramón y Cajal (RYC-2015-17809). 879 880 Author contributions: IM, JN, SG, JA, JDR and JBC conceived the presented idea, and planned 881 the overall experiments. BK, TOR, CF and KK contributed to the design and implementation of 882 the area-only target tracking research. IM developed the theory and performed the computations 883 and performed the numerical simulations. PJB and AA verified the analytical methods. PP, MV, 884 GV, NP and MC contributed to field tests preparation and the interpretation of the results. IM and 885 JN wrote the manuscript with input from all authors. SG, JBC and KK supervised the project. All 886 authors discussed the results and contributed to the final manuscript. 887 888 Competing interests: No competing financial interests exist. 889 890 Data and materials availability: GitHub github.com/imasmitja/TDOA_algorithms_SR 891 892 893 894 Science Robotics Manuscript Template Page 34 of 47 895 Fig. 1. Tracking deep water benthic marine animals. The strategy designed to track Norway 896 lobsters (Nephrops norvegicus) is represented. Four receivers created an acoustic LBL 897 localization system, where each one was in self-recording mode and was not accessed in real-898 time. The tags transmitted periodically an acoustic ping, which was recorded by the static 899 receivers and the underwater vehicles, both systems were used to track the lobsters’ movements. 900 Moreover, different pictures detailing operations and systems are included: (A) the canister used 901 to release the Norway lobsters, (B) the Super Mohawk II ROV, (C) a tagged Norway lobster 902 showing the Vemco tag glued on its superior portion of the cephalothorax ( manipulation of the 903 lobsters occurred in red light to avoid eye damaging), and (D) the Girona500 AUV. The 904 experiment was conducted in the northwest area of the Mediterranean Sea (E). 905 906 907 Science Robotics Manuscript Template Page 35 of 47 908 Fig. 2. TDOA algorithms performance. The CRB for TDOA target estimation (A), where the 909 red crosses represent static receivers creating an acoustic LBL system. Target trajectory designed 910 to compere the different TDOA algorithms’ performance (B) are presented in relation to the 911 particle filter (PF), the maximum a posteriori (MAP) estimation, the MAP marginalizing the latest 912 measures (MAP-M), the maximum likelihood (ML) estimation, the weighted least squares 913 (WLS), the WLS-ML, and the YAPS. The target estimation RMSE over the time (C), and the 914 RMSE over 100 Monte Carlo iterations for all the algorithms (D), where different TDOA noise 915 has been added ( 0.5 ms, 1 ms, 1.5 ms σ = ), the plots show the median, and 5th and 95th 916 percentiles. Finally, the average run-time required to compute one target position is shown (E). 917 918 Science Robotics Manuscript Template Page 36 of 47 919 Fig. 3. AOTT method applied to the Monterey Bay test-site. AOTT method representation (A). 920 Functions designed to weight the PF’s particles (B). MC simulations to find the optimal value for 921 different parameters (C), such as the circumference radius, the particles range, the reception, and 922 the PF’s resampling ratio, computed for static and moving targets. Simulations conducted to 923 observe the AOTT’s performance under different scenarios (D), where a reception ratio of 100% 924 and 60% were used over 100 MC iterations. Results obtained during a field test (E), where a CPF 925 (F) was used as a target and a Wave Glider (G) as an observer. Map of the study area in Monterey 926 Bay, CA, USA (H). 927 928 Science Robotics Manuscript Template Page 37 of 47 929 Fig. 4. Clock drift results during the Norway lobsters tracking. The synchronization process 930 of the receivers can be observed in the flowchart (A). Then, the four-step process and the results 931 obtained at each step are also presented as: the main drift at the beginning (B), the coarse drift 932 after the first step (C), the fine drift (D), the TDOA error result and its outliers (E), and the final 933 TDOA using a synchronization tag as a reference (F). C12, C13, and C14 denotes the difference 934 between two receiver clock drifts. In addition, the positions of the moorings and the lobster 935 canister in (G), where their initial position using the ship’s GPS, the position obtained using the 936 ROV’s USBL, and the position computed acoustically are also pictured. 937 938 Science Robotics Manuscript Template Page 38 of 47 939 Fig. 5. Norway lobster tracking results. (A) The trajectories conducted by the tagged Norway 940 lobsters during the moored experiment are represented, where the receivers are denoted as BS(X) 941 and each tagged lobster has a different color. (B) The accumulative traveled distance covered by 942 each tagged individual. (C) The different trajectories conducted by the underwater robots, in order 943 to localize and track the Norway lobsters, where the receivers’ localizations are represented by a 944 red X and the detected tags denoted as T0-T3. Finally, an image obtained with the ROV HD 945 camera (D), picturing the slope seabed with serval tunnel entrances and a wandering Norway 946 lobster (15 cm distanced laser green-beams dots can also be observed). 947 948 Science Robotics Manuscript Template Page 39 of 47 949 Fig. 6. Algorithms’ performance during the Norway lobster experiment. Synchronization tag 950 positions computed using different TDOA algorithms. These tags were attached on each mooring 951 line alongside with a Vemco acoustic receiver (BS) (A, B, D and E). A last tag was mounted on 952 the lobster canister, which was deployed on the center of the experiment (C). The error 953 covariance matrix with a confidence interval of 98% is also presented. This information is 954 presented as error bars in (F). The plots show the median, and 5th and 95th percentiles. 955 956 957 Science Robotics Manuscript Template Page 40 of 47 Table 1. Mooring lines positioning error. Position of the moorings obtained using the ROV and 958 the TOA signals (Position 1) compared with the positions computed using the WLS-ML 959 algorithm (Position 2), and the associated error. 960 Moorings Position 1 Position 2 Error (m) SD (m) x (m) y (m) x (m) y (m) BS(A) 543957.71 4651491.78 543958.24 4651491.60 0.74 0.55 BS(B) 544090.23 4651482.54 544089.80 4651482.68 0.68 0.40 BS(C) 544092.14 4651434.09 544092.10 4651433.95 1.29 1.05 BS(D) 543952.24 4651347.94 543952.56 4651348.23 0.51 0.29 BS(E) 544114.21 4651341.26 544114.03 4651341.52 0.48 0.27 961 962 Science Robotics Manuscript Template Page 41 of 47 963 Ref. (26) ( 14, 27 ) (14) (28) (29) ( 21, 75 ) (76) (24) (32) (20) (30) (31) (77) (78) (79) ¶Species names: Acipenser oxyrinchus oxyrinchus, Pseudopleuronectes americanus, Micropogonias undulatus, Anoplopoma fimbria, Paralithodes camtschaticus , Carcharodon carcharias, Carcharias taurus, Triakis semifasciata, Dermochelys coriacea, Epinephelus morio, Lutjanus campechanus, Chionoecetes opilio, Gadus morhu a, Oncorhynchus tshawytscha, and Paralithodes camtschaticus NR = Information not reported Dynamic transponders Vehicle REMUS -100 AUV REMUS -100 AUV REMUS -100 AUV Slocum G2 Glider REMUS 100 AUV REMUS 100 AUV Slocum G2 Glider Iver2 AUV Iver2 AUV REMUS 100 AUV Slocum G1 Glider Wave Glider Slocum G2 Glider REMUS 100 AUV Saildrone ASV Error (m) 25 - NR - - NR - 80 ~10 NR NR NR NR NR NR Algorithm SYNAPS ‡ - SPLWCA - - NR - PF PF NR - DWA ‡‡ BBMM*** SYNAPS ‡ NR Acoustic range (m) NR ~890 NR NR ~500 NR 250 NR NR NR NR ~500 ~1000 ~500 NR Method TDOA Presence SPL Presence Presence USBL Presence Bearing Bearing USBL Presence Presence Presence TDOA NR ††The SmartTag package consist of an IMU, a Lotek MM-M-16-50-PM acoustic tag, a VHF transmitter and a video logger system * Using ALPS: Asynchronous Logging Positioning System software (Lotek Wireless, Inc.) Static transponders Error (m) - ~2 - - NR - - - - - NR NR NR - - Algorithm - TDOA* / SPLWCA - - NR - - - - - - VPS** BBMM*** - - Acoustic range (m) - NR - ~800 NR - ~800 - - - NR NR ~1000 - - Method - LBL - Gate F. § NR - Gate F. § - - - Presence LBL Fisheries F. § - - Depth (m) ~18 < 8 10 ~90 < 585 93 -130 <25 < 100 < 10 0 -20 30 -60 ~116 < 50 30 -100 <100 ‡SYNAPS: synthetic hydrophone array, proprietary software from Lotek Time 2 d 51d 7 h 1095 d 61 d 12 h 12 d 3 d 3 d 36 h 365 d 720 d 720 d 7 d 365 d Tag size (mm) 32x101 ○ ○ 16x54 62x16 76x380 16x54 16x80 200x127 76x380 13x36 13x36 16x54 8.5x43 NR ○= Information not found - = Information not applicable †Designed by WHOI § See M. R. Heupel et al. (16) **VPS: Vemco Positioning System array ***BBMM: Brownian Bridge Movement Model (80) ‡‡DWA = Daily Weighted Average Tag type MAP32 1s MA - TP11 18 MA - TP11 18 V16 69kHz MA -TP16-33 Transp. † V16 69kHz MM -M-16-50 Smart Tag †† Transp. † V13P L V13 and V9 V16 -6H MM -M-8-S0 Vemco # 2 39 1 4 41 6 292 1 3 9 61 164 317 20 150 Species A. O . oxyrinchus P. americanus M. undulatus A. O . oxyrinchus A. fimbria and P.camtschaticus C. carcharias C. taurus T. semifasciata T. semifasciata D. coriacea E. morio and L. campechanus C. opilio G. morhua O. tshawytsch a P. camtschaticus Location Hudson R. Navesink R. Caribbean S. NW Atlantic N Pacific NE Pacific NW Atlantic NE Pacific NE Pacific NE Pacific G . Mexico NW Atlantic NW Atlantic G. Alaska Bering S. Year 2008 2013 2013 2013 2014 2015 2015 2016 2017 2018 2018 2019 2019 2019 2019 964 Table 2. Target tracking experiments using underwater robots. Different campaigns conducted to localize and follow a marine tagged animal using both fixed receivers and/or underwater vehicles (aka dynamic transponders)