Autonomous Underwater Vehicle Localization Using Sound Measurements of Passing Ships

: This paper introduces the localization method of an Autonomous Underwater Vehicle (AUV) in environments (such as harbors or ports) where there can be passing ships near the AUV. It is assumed that the AUV can access the trajectory and approximate source level of a passing ship, while identifying the ship by processing the ship’s sound. This paper considers an AUV which can localize itself by integrating propeller and Inertial Measurement Units (IMU). Suppose that the AUV has been moving in underwater environments for a long time, under the IMU-only localization. To ﬁx long-term drift in the IMU-only localization, we propose that the AUV localization uses sound measurements of passing ships whose trajectories are known a priori. As far as we know, this AUV localization method is novel in using sound measurements of passing ships of which the trajectories are known a priori. The performance of the proposed localization method is veriﬁed utilizing MATLAB simulations. The simulation results show signiﬁcant estimation improvements, compared to IMU-only localization. Moreover, using measurements from multiple ships gives better estimation results, compared to the case where the measurement of a single ship is used.


Introduction
The localization problem of an Autonomous Underwater Vehicle (AUV) is a critical problem for the navigation of the AUV. This paper presents the localization of an AUV in environments (such as harbors or ports) where there can be passing ships near the AUV. In practice, there may be ships with known trajectories. We assume that the AUV can access the trajectory and approximate Source Level (SL) of a passing ship, while identifying the ship by processing the sound of the ship. The proposed localization method uses sound bearing measurements of a passing ship moving along a known trajectory (planned routes).
We use the planned routes of large commercial ships. Large commercial ships make large noise, thus their sound signals travel a long distance in underwater environments. Hence, the sound information of large commercial ships is valuable for the proposed localization. The information of ships' trajectories is accessible, since a ship uses Automatic Identification System (AIS) nowadays. AIS records the position, routes, destinations of large commercial ships worldwide. The AIS of commercial ships is available using the following website: https://www.marinetraffic.com/.
It is feasible that the AUV can access approximate Source Level (SL) of a passing ship, while identifying the ship by processing the sound of the ship. McKenna et al. [1] introduced metrics of ship noise, including spatial radiation patterns, source levels, spectral characteristics, and sound exposure levels for seven merchant ship-types: container ships, vehicle carriers, bulk carriers, open-hatch cargo ships, and chemical, crude oil and product tankers. It is feasible to gather the SL information of large commercial ships using the estimation approach in [2]. Jiang et al. [3] acquired To solve the BOTMA problem, Range Parameterized EKF (RPEKF) was used to reduce the range estimation error [37,38]. The EKF is suitable for real-time processing of nonlinear systems, since it is computationally efficient [34]. The RPEKF is composed of many independent EKFs, such that every EKF has a different initial range estimation. The combined target estimation is derived as a weighted sum of every EKF output.
If the AUV cannot measure the sound of any ship, then it relies on locating itself by integrating propeller and IMU (IMU-only localization). We assume that the AUV has been moving in underwater environments using the IMU-only localization before time step 0. Note that the range between a ship and the AUV is not accessible at time step 0, due to long-term drift of the AUV. Thus, this paper modifies the RPEKF to decrease the range estimation error.
Moreover, there may be a case where the long-term drift of the AUV is too severe. Assume that the AUV can access the approximate SL of a passing ship. Based on the SL, we can utilize an acoustic propagation model, which simulates the sound emanated from the ship. Environmental factors, such as underwater obstacles or topographical maps, can be taken into account in the modeling of an acoustic propagation model. The Range-dependent Acoustic Model (RAM) [39][40][41] has been widely used to simulate how the sound signal of the ship propagates in underwater environments. Using RAM, we can derive the feasible region where the AUV can exist, in order to measure the sound of the ship. Then, the RPEKF is applied to estimate the AUV position, considering the feasible region.
The contributions of this paper are summarized as follows. This paper assumes that the AUV can access the trajectory and source level of a passing ship, while identifying the ship by processing the ship's sound. To fix long-term drift in the IMU-only localization, we introduce the AUV localization method using sound measurements of passing ships that travel along known routes. Under the proposed localization method based on moving ships, we do not have to create additional infrastructure for AUV localization. Moreover, we can localize the AUV in a stealthy and power-efficient manner, since the AUV does not emit a signal during this passive localization. As far as we know, this AUV localization method is novel in using sound measurements of passing ships of which the trajectories are known a priori. This paper is novel in presenting the underwater localization method using RAM and sound measurements of ships moving along known routes. We verify the effectiveness of the proposed localization method utilizing MATLAB simulations. This paper is organized as follows. Section 2 presents the preliminary information of the paper. Section 3 presents the localization method using bearing-only measurements of passing ships. Section 4 presents MATLAB simulations to verify the effectiveness of the proposed localization method. Section 5 provides Conclusions.

Extended Kalman Filter (EKF)
The range between a ship and the AUV is not accessible at time step 0, due to long-term drift of the AUV. Thus, this paper uses the RPEKF to decrease the range estimation error. The RPEKF is composed of many independent EKFs, such that every EKF has a different initial range estimation.
As the preliminary information of this paper, we introduce the EKF in [34] briefly. We consider discrete-time systems. dt indicates the sampling duration of our discrete-time system. At sample-index k, we can define a nonlinear process model as Here, X k denotes the state vector at sample-index k. Also, w p k is a Gaussian noise with mean 0 and covariance matrix Q k .
At sample-index k, we can define a nonlinear measurement model as Here, w m k is a Gaussian noise with mean 0 and covariance matrix R k . LetX k|k denote the state vector of the EKF, utilizing all measurements up to sample-index k. In addition, P k|k denotes the error covariance ofX k|k . The EKF is composed of the prediction step and the measurement update.
Using Equation (1), we predict the state vector of the EKF as We then predict the error covariance matrix as Here, F k is the partial derivative of f (X) with respect to X evaluated at X = X k|k . Using the measurement z k in Equation (2), the state vector of the EKF is updated as Here, we use where H k is the partial derivative of h(X k ) with respect to X k evaluated at X k = X − k . The error covariance is updated using Then, we update both the state vector and its error covariance usinĝ

Derive the Maximum Sensing Range Associated to a Ship
As the preliminary information of this paper, we introduce the propagation model which simulates a sound field in underwater environments [41]. While the surface target's sound signal propagates in underwater environments, a sound field generated from the target gets complicated due to environmental interference [42]. Numerical propagation models, such as ray theory and parabolic equation, were widely utilized to simulate a sound field [41].
In the parabolic equation, a small-angle approximation approximates the elliptic wave equation obtained from the Helmholtz equation. However, this results in error at steep angle spaces. RAM reduces this error utilizing a wide-angle parabolic equation with Pade approximation [39][40][41]. In this way, RAM can simulate a sound field accurately.
Next, several definitions in [42] are discussed. Source Level (SL) indicates the sound strength generated by a source. We assume that the approximate SL of a ship is known a priori. Transmission Loss (TL) is the reduction of sound while the sound emanates from the source. TL profile is the function of both target range and depth.
TL(C, target) represents TL of the target's signal as the AUV is at a position C and the source is at a position target. TL(C, target) is constructed utilizing the RAM. Let M(C) denote the sound strength as the sound signal from the target at target reaches the AUV in C. M(C) is estimated as Note that the AUV at C fails to measure the sound of the target ship at target in the case where M(C) is less than zero. Using Equation (9), the AUV cannot measure the ship's sound in the case where is met. The AUV's depth can be measured using pressure sensors. Thus, the maximum sensing range can be derived by searching for a target range, which satisfies Equation (10).

RPEKF
The RPEKF in [37,38] is utilized to decrease the range estimation error. In the RPEKF, we utilize M EKF tracks to decrease the range estimation error. Every EKF track of the RPEKF runs the EKF in Section 3.2 independently, such that every EKF has a different initial range estimate.
We assume that the AUV has been moving in underwater environments using the IMU-only localization before time step 0. Suppose that at time step 0, the IMU-only localization estimates the range between the ship and the AUV as ∆. Also, suppose that the range estimation error at time step 0 is bounded above by δ. Doing underwater experiments, we can derive the relationship between the operation time of the AUV (under the IMU-only localization) and the increase in localization error. Then, the bound for localization error can be set considering the operation time of the AUV in underwater environments. This range estimation error is unavoidable, since the AUV has been moving in underwater environments using the IMU-only localization.
A feasible range interval (r min − r max ) is set considering the localization error of the AUV at time step 0. Here, r min is set as Also, we use Moreover, there may be a case where the long-term drift of the AUV is too severe. Assume that the AUV can access the approximate SL of a passing ship. Based on the SL, we can utilize an acoustic propagation model, which simulates the sound emanated from the ship. Suppose that the AUV measures the sound of the ship. Then, using the RAM in Section 2.2, we can derive the feasible region where the AUV can exist.
For instance, suppose that the maximum sensing range, say mS, is derived by searching for a target range satisfying Equation (10). In this case, we use instead of Equation (11). Also, we use instead of Equation (12).
To maintain a comparable performance of every EKF track, the range interval (r min − r max ) is divided such that the j-th sub-interval has its range interval as (rn(j) − rn(j + 1)) for j ∈ {1, . . . , M}.
k|k denote the state vector of the j-th EKF track, utilizing all measurements up to sample-index k. In addition, P (j) k|k denotes the error covariance ofX k|k−1 denote the state vector of the j-th EKF track, utilizing all measurements up to sample-index k − 1. Let In Equation (17), the likelihood p(z k |j) is derived as Here, The state vector and its error covariance at sample-index k are calculated aŝ At sample-index 0, we set w

System Models
This manuscript considers horizontal AUV motion in two dimensions. In other words, the AUV maneuvers while maintaining a certain depth. The AUV's depth can be measured fairly well by utilizing pressure sensors.
We first introduce some notations used in this paper. I r presents the identity matrix with r rows and r columns. Let diag(a, b, c, . . . ) indicate a diagonal matrix with diagonal elements a, b, c, . . . in this order. Cov[q] is used to present the error covariance of a variable q. Let atan2(y, x) denote the argument (also called phase or angle) of the complex number, say x + iy.
Suppose that the AUV measures the sound of N ships in total. [x t,n k , y t,n k ] (n ∈ {1, 2, . . . , N}) is the vector representing the position of the n-th ship at sample-index k. [ẋ t,n k ,ẏ t,n k ] T is the 2D vector representing the velocity of the n-th ship at sample-index k.
We assume that the AUV can access the trajectory of a passing ship, while identifying the ship by processing the sound of the ship. This implies that the AUV at sample-index k can access [x t,n k , y t,n k ,ẋ t,n k ,ẏ t,n k ] for all n ∈ {1, 2, . . . , N}.
denote the state vector representing the AUV state.
Recall that dt denotes the sampling interval in the discrete-time systems. The discrete-time process model of the AUV is The AUV velocity [ẋ o k ,ẏ o k ] T can be estimated by integrating propeller and IMU. We can estimate the AUV velocity using the integration approach in [26,27]. Considering the process noise in the estimation, Equation (23) can be rewritten as In (24), F is Also, A k is In (24), presents the process noise. Here, [a x , a y ] T presents the AUV's 2D acceleration noise, which is unknown a priori. In w p k , we utilize The process noise w p k is added in (24), considering the error in estimating [ẋ o k ,ẏ o k ] T as we integrate propeller and IMU. In (24), A k + w p k is assumed to be derived by integrating propeller and IMU [26,27]. For instance, unknown sea currents can generate large process noise w p k [26]. In (24), w p k is a Gaussian noise with mean 0 and covariance matrix Q k .
The proposed localization scheme can be considered as the reverse process of the BOTMA. The proposed localization method uses sound bearing measurements of a passing ship moving along a known trajectory (planned routes). Let z n k denote the bearing measurement associated to the n-th ship (n ∈ {1, 2, . . . , N}) at sample-index k. The bearing measurement model associated to the n-th ship (n ∈ {1, 2, . . . , N}) is Here, we utilize h n (X k ) = atan2(y r,n k , x r,n k ).
Here, x r,n k = x t,n k − x o k , and y r,n k = y t,n k − y o k . In (29), w m k is the measurement noise with following distribution: w m k ∼ N(0, R k ). Here, N(0, R k ) represents Gaussian distribution with mean 0 and covariance matrix R k .
Recall that H k is the partial derivative of h(X k ) with respect to X k evaluated at X k = X − k . See Equation (6). We next present how to derive ∂h n ∂X k . We utilize We next present how to derive each term in (31). Using Equation (30), we have Here, r n k = (x r,n k ) 2 + (y r,n k ) 2 . Furthermore, we utilize Using Equation (30), we have

Run the RPEKF Based on Bearing-Only Measurements of Passing Ships
This subsection presents running the RPEKF based on bearing-only measurements of passing ships. We discuss how to update the EKF track in the RPEKF (See Section 2.3).
Recall that the superscript (j) in Section 2.3 denotes the index of each EKF track in the RPEKF (See Section 2.3). The EKF is composed of two steps: prediction step and measurement update step.

Prediction Step
We first present the prediction step. As the process model (1), we utilize the linear model in (24). Since our model is linear, the prediction step is Based on Equation (4), we further predict the error covariance matrix using

Measurement Update Step
We next present the measurement update step in the EKF. Suppose that the AUV measures the bearing of sound of N k ships at sample-index k. Let {n k,1 , n k,2 , . . . , n k,N k } denote the index of the measured ship at sample-index k. Since there are N ships in total, we have {n k,1 , n k,2 , . . . , n k,N k } ⊂ {1, 2, . . . , N}.
Using Equation (29), the measurement associated to the n k,s -th ship (s ∈ {1, 2, . . . , N k }) is In each EKF track, we perform the EKF measurement update associated to each ship in {n k,1 , n k,2 , . . . , n k,N k }.
The EKF measurement update indicates that we apply the EKF Equations (6) and (7) sequentially.
We have N k measurements at sample-index k. Thus, each EKF track keeps running the Equations (6) and (7) N k times iteratively, while using measurements in Equation (38). In other words, each EKF track uses z n k,1 k as the first measurement at sample-index k. The EKF track applies the Equations (6) and (7), while using z n k,1 k as the measurement z k in Equation (5). Next, the track uses z n k,2 k as the second measurement at sample-index k. The EKF track applies the Equations (6) and (7), while using z n k,2 k as the measurement z k in Equation (5). Iterate this until the track uses z n k,N k k as the last measurement at sample-index k. The EKF track applies the Equations (6) and (7), while using z n k,N k k as the measurement z k in Equation (5).
In this way, each EKF track can process the sound of N k ships at sample-index k. Once the measurement update is done using N k measurements at sample-index k, each EKF track updates both the state vector and its error covariance using Equation (8).

Outlier Measurement Handing
In some cases, a bearing measurement z n k,s k may not be accurate enough to update the KF estimation. Considering this aspect, a bearing measurement z n k,s k is discarded in the case where a bearing measurement z n k,s k satisfies the following condition. where , and Σ > 0. Recall the measurement Equation (2). In the case where the measurement z n k,s k has a Gaussian distribution with mean 0 and covariance R k , µ k is also Gaussian with mean 0 and covariance H k P (39) implies that this paper does not utilize a bearing measurement, which is npt reliable considering the KF estimation. If the condition in (39) is satisfied, then we do not perform the measurement update associated to the bearing measurement.

MATLAB Simulations
We use MATLAB (MathWorks, Natick, MA, USA) simulations to demonstrate the performance of the proposed localization approach. The scenario runs for 3000 s. The sampling duration is dt = 1 second.
In Equation (29), the bearing measurement noise is set as σ b = 0.5 degrees. This bearing noise model is commonly used in the literature on target tracking based on bearing measurements [32][33][34][35].
This paper considers the AUV that can locate itself by integrating a propeller and IMU. We set the noise in the simulations as follows. Let σ a represent the standard deviation of a. In Equation (27), we utilize σ a x = σ a y = 2 (m/s 2 ). This process noise is set considering the error in estimating [ẋ o k ,ẏ o k ] T , as the AUV integrates propeller and IMU for its localization [26,27].
Using the IMU-only localization, the AUV has been moving in underwater environments before time step 0. Suppose that at time step 0, the IMU-only localization estimates the range between the 1-st ship and the AUV as ∆ = 1500 m. Also, suppose that the range estimation error at time step 0 is bounded above by δ = 1000 m. Using Equations (11) and (12), we set r min = 500 m, and r max = 2500 m. Also, M = 5 EKF tracks are used in total.
Using Taylor expansion onX j 0|0 (1 : 2) in Equation (41), we obtain where Furthermore, we have Considering the j-th EKF track, P j 0|0 is selected as Here, V max is the maximum speed of an AUV, which is assumed to be known a priori. We utilize V max = 5 m per second. In Equation (45), 0 2 indicates a 2 by 2 zero matrix.
We run M c = 20 Monte-Carlo simulations. E m k where m ∈ {1, 2, . . . , M c } is the AUV position estimation at sample-index k utilizing the m-th Monte-Carlo simulation. E m k − E k is the error of the AUV position estimation at sample-index k.
At the m-th Monte-Carlo simulation (m ∈ {1, 2, . . . , M c }), one considers the estimation error, given by Then, at each sample-index k, the following metric is utilized: Also, at each sample-index k, we define In simulations, we plot meanE k , maxE k , and minE k with respect to sample-index k. Figure 1 presents the scenario using one ship. The speed of the ship is 2 m/s, and that of the AUV is 2 m/s. The AUV starts from the origin. The ship starts from (2000,0) in meters, and the trajectory of the ship at every minute is plotted with red asterisks. The starting point of the ship is marked with a diamond. The trajectory of the AUV at every minute is depicted with blue circles. Consider the case where the proposed localization method is utilized. In Figure 2, we plot meanE k , maxE k , and minE k with respect to sample-index k. The AUV maneuvers at 750 and 2250 s. At these moments, the estimation error decreases sharply. This is similar to the case where the maneuver of an observer makes the BOTMA system observable [36].

One Ship Scenario
To run one Monte-Carlo simulation takes about 1.3 s utilizing MATLAB simulations. See that the computational load of the proposed localization scheme is very low.

Three Ships Scenario
We next consider the scenario with three ships. The speed of every ship is 2 m/s, and that of the AUV is 2 m/s. The AUV starts from the origin. In Figure 3, the starting point of each ship is marked with a diamond. At every minute, the trajectory of every ship is plotted with red asterisks. At every minute, the trajectory of the AUV is depicted with blue circles.
Consider the case where the proposed localization method is utilized, while using three ships. In Figure 4, we plot meanE k , maxE k , and minE k with respect to sample-index k. See that the estimation error decreased considerably, compared to Figure 2. This is due to the fact that the AUV localization uses three bearing lines emitted from three ships, instead of one bearing line emitted from one ship. To run one Monte-Carlo simulation takes about 2.6 s utilizing MATLAB simulations. See that the computational load of the proposed localization scheme is very low. Consider the case where the proposed localization method is utilized. We plot meanE k , maxE k , and minE k with respect to sample-index k. The AUV maneuvers at 750 and 2250 s. At these moments, the estimation error decreases sharply. This is similar to the case where the maneuver of an observer makes the Bearing-Only Target Motion Analysis (BOTMA) system observable [36]. At every minute, the trajectory of every ship is plotted with red asterisks. At every minute, the trajectory of the AUV is depicted with blue circles (three ships scenario). Figure 5 is the enlarged figure of Figure 4. Since we utilize the bearing measurements of three ships, the estimation error is small from the beginning of the scenario.

IMU-Only Localization
For comparison, we consider the case where the AUV localizes itself by integrating propeller and IMU. Note that no ships are used in this IMU-only localization. The movement of the AUV is identical to that in Figure 1. Similarly to Equation (24), the AUV can locate itself using    Figure 4. We plot meanE k , maxE k , and minE k with respect to sample-index k. Since we utilize the bearing measurements of three ships, the estimation error is small from the beginning of the scenario (three ships scenario).
Here, A k is in Equation (26), and w p k is in Equation (27). A k + w p k in (50) is assumed to be derived by integrating a propeller and IMU [26,27]. The error in the IMU integration is modeled by adding Gaussian noise w p k in Equation (50). As we increase the variance in w p k , the error in the IMU-only localization increases. Recall that in Equation (27), we utilize σ a x = σ a y = 2 (m/s 2 ).
Recall that at time step 0, the IMU-only localization estimates the range between the 1-st ship and the AUV as ∆ = 1500 m. Hence, similarly to Equation (41), X 0 in Equation (50) is initialized using Consider the case where the IMU-only localization is used. In Figure 6, we plot meanE k , maxE k , and minE k with respect to sample-index k. Since no ship is used to fix long-term drift, the estimation error does not decrease as time goes on. Compared to Figures 2 and 4, the estimation error increases. To run one Monte-Carlo simulation takes only 0.6 seconds utilizing MATLAB simulations. We plot meanE k , maxE k , and minE k with respect to sample-index k. Since no ship is used to fix long-term drift, the estimation error does not decrease as time goes on. Compared to Figures 2 and 4, the estimation error increases.

Conclusions
This paper presents the localization method of an AUV in environments, such as harbors or ports, where passing ships can exist. Assume that source levels and trajectories of the passing ships are known a priori. This paper introduces the AUV localization method based on bearing-only measurements of passing ships moving along known trajectories (planned routes). The RPEKF is used to reduce the range estimation error in bearing-only measurements.
As far as we know, this paper is novel in presenting the underwater localization method using RAM and sound measurements of ships moving along known routes. If we do not use the proposed localization, the AUV needs to depend on IMU-only localization. It is inevitable that the IMU-only localization leads to the integration of localization error as time goes on. The simulation results in Sections 4.1 and 4.3 show that using the measurements of a single ship can improve the localization accuracy, compared to IMU-only localization.
Moreover, the simulation results in Sections 4.1 and 4.2 show that using measurements from multiple ships gives better localization accuracy, compared to the case where the measurement of a single ship is used. As our future works, we will verify the effectiveness of the proposed localization method using experiments with a real AUV.
There may be a case where there exist many ships of the identical type (for example, container ships). In this case, it is not trivial to identify a ship by processing the sound of the ship. In the case where a ship is not correctly identified, we may generate false measurements, similar to clutters in BOTMA problems [43]. In the future, we will tackle the case where there exist many ships of the identical type.
The proposed localization scheme can be further applied for localization of various autonomous vehicles, such as an Unmanned Aerial Vehicle (UAV). The UAV can measure the directions of signals emitted from airplanes (or sites) of which trajectories (locations) are known a priori. Then, bearing direction measurements can be used for localization of the UAV. This localization scheme can be used when Global Positioning Systems are not available for the UAV.