A Research on the Simultaneous Localization Method in the Process of Autonomous Underwater Vehicle Homing with Unknown Varying Measurement Error

: We propose an acoustic-based framework for automatically homing an Autonomous Underwater Vehicle (AUV) to the ﬁxed docking station (F-DS) and mobile docking station (M-DS). The proposed framework contains a simultaneous localization method of AUV and docking station (DS) and a guidance method based on the position information. The Simultaneous localization and mapping (SLAM) algorithm is not available as the statistical characteristics of the measurement error of the observation system are unknown. To solve this problem, we propose a data pre-processing method. Firstly, the measurement error data of acoustic sensor are collected. Then, We propose a Variational Auto-Encoder (VAE) based Gaussian mixture model (GMM) for estimating the statistical characteristics of measurement error. Finally, we propose a support vector regression (SVR) algorithm to ﬁt the non-linear relationship between the statistical characteristics of measurement error and its corresponding working distance. We adopt a guidance method based on line-of-sight (LOS) and path tracking method for homing an AUV to the ﬁxed docking station (F-DS) and mobile docking station (M-DS). The lake experimental data are used to verify the performance of the localization with the estimated statistical characteristics of measurement error.


Introduction
Autonomous underwater vehicles (AUVs) play important roles in underwater missions. However, the working time is limited by on-board energy and storage capability of AUVs. The autonomous recovery of AUVs is used for recharging, data transfer, and relaunch. Docking stations can either be stationary [1][2][3][4][5][6] or mobile [7][8][9][10][11]. The process of AUV sailing from far-end to the docking station is called homing. The homing process aims at the precise physical contact between AUVs and docking stations . Physical contact can be passive or active. A common kind of passive contact is that an AUV is captured by a robotic manipulator. A common kind of active contact is that an AUV catches a cable or pole [12][13][14] with a latch or hook. The other kind of active contact is using funnel-like [15][16][17] docking stations. In this work, We address the recovery of an AUV by actively docking with a funnel-like docking station.
One of the key technologies in homing is navigation. The positions of AUV and DS should be simultaneously sensed by AUV. Three kinds of sensors are available for the position sensing of DS: (1)electromagnetic sensors, (2) acoustic sensors, and (3) optical sensors. Acoustic sensors have the longest working range, but the lowest precision. Woods Hole Oceanographic Institution (WHOI) [18,19] conducted a research on the AUV homing method based on acoustic guidance. The transceiver installed in the head of the AUV receives an acoustic signal from the beacon installed on the DS, and measures its position relative to the beacon. The absolute position of AUV is calculated indirectly as it knows the absolute position of the beacon. Electromagnetic sensors have a better precision, but much less range, about 30-50 m. Massachusetts Institute of Technology (MIT) [20] proposed the electromagnetic guidance method. In this way, the final docking accuracy can reach 20 cm. However, due to some physical factors, the effective docking range is 20-30 m. Optical sensors have the shortest range, but the highest precision. National Radiology Solutions (NRaD) [21] proposed an optical docking system, which was demonstrated to be accurate and robust for vehicle terminal guidance and provided targeting accuracy on the order of 1 centimeter under real-world conditions, even in turbid bay water. According to the characteristics of the sensors, they usually combined acoustic and optical sensors for AUV recovery. In this work, we focused on the development of acoustic-based AUV recovery algorithms since the vision guidance technology is maturing. More importantly, if the positional error of acoustic sensor is large, it may result in the inability to enter the terminal guidance phase based on vision, resulting in a failed docking.
Acoustic methods can be used at large distance. One such method for localization, is range-only localization. Range-only measurements are extensively used in many underwater homing application [22][23][24][25][26]. However, these measurements present a high-nonlinear and standard Extended Kalman Filters (EKFs) cannot cope with it. The solution to the range-only localization for homing problem using a sum of Gaussian filter is proposed in Reference [27]. The other method for localization, is range-bearing localization. In this work, the range and bearing of beacon are provided by the acoustic sensor. In other words, we can obtain the projections of distance from the transceiver to the beacon in the x and y directions of the observation coordinate system. Some SLAM algorithms for AUV have been applied to sense the position of a robot and map its surroundings simultaneously [28][29][30][31]. These algorithms are not suitable for our recovery system as the statistical characteristics of our measurement error are non-Gaussian and unknown. We adopt the improved SLAM algorithm to simultaneously sense the positions of AUV and DS.
The novel contributions of this paper include: .
We propose a method to collect the measurement error data for the observation system with the unknown statistical characteristics of measurement error. .
We propose a Variational auto-encoder based on a Gaussian mixture model to estimate the overall distribution of all measurement error data. The advantage of this method is that it can be used in both Gaussian and non-Gaussian cases. .
We adopt the support vector regression algorithm to divide the overall distribution into multiple single Gaussian distributions and their corresponding working ranges, and fit the non-linear relationships between them. It is always available when the working range of acoustic sensor exceeds the maximum working range corresponding to the measurement error data.

System Overview
In this section, we introduce an overview of our recovery system. We illustrate the recovery system in Figure 1. The recovery system includes a mother vessel (called AUV-I) and a sub-AUV (called AUV-I I) for recovery. The DS was rigidly fixed to the underbelly of AUV-I. The entrance of the DS was funnel-like and was 2 m in diameter. When we need an F-DS, AUV-I remains hovering. When we need a M-DS, AUV-I moves at a constant speed in a straight line. AUV-I I was a torpedo-shaped vehicle that was 384 mm in diameter and 5486 mm in length. It was equipped with a Doppler velocity log (DVL), inertial measurement unit (IMU), GPS and an acoustic sensor. The acoustic senor-Evologics 32C R-was used for long-range positioning. The IMU used in our recovery system is manufactured by a Chinese scientific research institution with high precision. We propose a data pre-processing method for SLAM with unknown measurement error. The overview of our system is shown in Figure 2.

Data Pre-Processing
The working medium of an acoustic sensor is sea water. The ocean is a very complex acoustic medium, whose complexity is mainly reflected in different acoustic propagation rules with the change of sea area and season. Due to the viscosity and heat conductivity of seawater, some of the sound energy is absorbed by seawater. Air bubbles, suspended particles, uneven water masses, plankton, and uneven boundaries in seawater can cause acoustic signal scattering. The above phenomenons will cause the energy loss of the signal from the source in the process of transmission. In addition, due to the non-uniform sound velocity in sea water, the sound line will be refracted, so the acoustic source signal may not reach the receiving point at a long distance. On the other hand, due to the existence of a large number of scatterers and uneven interfaces in the ocean, when the acoustic source emits sound waves, it will encounter these scatterers and generate scattered waves, which will interfere with the reception of acoustic signals. In conclusion, due to the working medium and working principle of the acoustic sensor, the measurement results of the acoustic sensor will have certain errors, and the farther the working distance is, the greater the interference will be and the greater the error will be. Table 1 shows the navigation data of our recovery system in a sea trial. (x AUV , y AUV ) and (x DS , y DS ) are the positions of AUV and DS measured by IMU. These positions are true as the IMU used by our AUV has very high precision. (x DS , y DS ) are the true distance between AUV and DS in X and Y direction of observation coordinate system. Z x and Z y are the measurements of acoustic sensor of AUV. We define measurement error as follows: The Figure 3 shows all measurement errors of this trail.  One of the premises of the SLAM algorithm is to know the statistical characteristics of the measurement error. In order to better describe the statistical characteristics of the measurement error, it is assumed that the error consists of two parts, namely, system error and random error, that is, w = ε + . In the formula, w represents the measurement error. ε represents the system error. is random error and follows the zero-mean Gaussian distribution. ∼ N(0, σ 2 ), where σ 2 is the variance of random error. In order to estimate the statistical characteristics of measurement error, firstly, the measurement error data are accumulated, then the statistical characteristics of the error are estimated by Variational auto-encoding (VAE) algorithm, and the relationship between the statistical characteristics of the measurement error and working distance is fitted by support vector regression (SVR).

Method for Accumulating Measurement Error Data
The motion model of the AUV is as follows: where v k is the forward velocity of AUV at kth time. G k is the yaw angle of AUV at kth time, which constitutes the input state of AUV, that is, dt is the sampling interval time. The measurement coordinate system is shown in Figure 4: The measurement model corresponding to the coordinate system shown in Figure 4 is as follows: The distance between AUV and DS at kth time is calculated by the Equation (4): Assuming that we know the initial position of AUV, first of all, (x v,k , y v,k ) denotes the position of AUV at kth time. (x f ,k , y f ,k ) denotes the position of F-DS or M-DS at kth time. (z x,k , z y,k ) denotes the measurement of acoustic sensor at kth time. Under the premise of only using an acoustic sensor, the method 1 is used to obtain measurement error data in the process of AUV homing to F-DS: for k = 1, · · · 5. According to the position of AUV at 1 : kth time, (x f , y f ) and Equation (3), the true measurements at 1 : kth time, that are, z x,1 , · · · ,z x,k ,z y,1 , · · · ,z y,k are calculated. 6. The measurement errors at 1 : kth time are the difference between z x,1 , · · · ,z x,k ,z y,1 , · · · ,z y,k and z x,1 , · · · , z x,k , z y,1 , · · · , z y,k . 7. The distances between AUV and F-DS at 1 : kth time, that is, dis 1 , · · · , dis k are calculated.
The method 2 is used to obtain measurement error data in the process of AUV homing to M-DS: for k = 1, · · · 1. According to the initial position, heading and speed of AUV and Equation (2) According to the heading and speed of M-DS, the distances travelled by M-DS in dt time, that is, position, dx f ,k , dy f ,k and Equation (2), the true positions of M-DS at 1 : kth time, that is, (3), the true measurements at 1 : kth time, that is, z x,1 , · · · ,z x,k ,z y,1 , · · · ,z y,k are calculated. 7. The measurement errors at time 1 : kth are the difference between z x,1 , · · · ,z x,k ,z y,1 , · · · ,z y,k and z x,1 , · · · , z x,k , z y,1 , · · · , z y,k . 8. The distances between AUV and M-DS at 1 : kth time, that is, dis 1 , · · · , dis k are calculated.
To verify the effectiveness of method 1 and method 2, we use IMU to sense the positions of AUV, F-DS and M-DS, namely x According to Equation (3), the measurements at 1 : kth time, that is, z x,1 ,z y,1 , · · · , z x,k ,z y,k are calculated. The measurement errors at 1 : kth time are the difference between z x,1 ,z y,1 , · · · , z x,k ,z y,k and z x,1 , · · · , z x,k , z y,1 , · · · , z y,k . We utilize three sets of experimental data of our recovery system. Figures 5 and 6 show the true and estimated measurement error data and their corresponding operating distance in AUV homing to F-DS, respectively. Figure 7 is the comparison between them. The Figures 8 and 9 are the true and estimated measurement error data and their corresponding operating distance in AUV homing to M-DS, respectively. Figure 10 is the comparison between them. The comparison results show that the system error of estimated measurement error is close to that of true measurement error. Although there are some deviations between them, the estimated measurement error can reflect the statistical characteristics of measurement error when only using acoustic sensor. The methods mentioned above are not suitable for the long-distance situation, because the mileage error of AUV for long-distance is large, which results in a large gap between the estimated measurement error and the true measurement error.

Variational Auto-Encoder Based on Gaussian Mixture Model
The measurement error corresponding to each operating distance follows a single Gaussian distribution, and the measurement errors corresponding to multiple working distances follow the Gaussian mixture distribution as a whole. We adopt a variational auto-encoder based on a Gaussian mixture model to estimate the statistical characteristics of the measurement error data.

Variational Auto-Encoder
Variational Auto-Encoder (VAE) [32,33] is a deep generation model, which aims to build a distribution model p model to approximate unknown data distribution p data . The data set contains N independent and identically distributed discrete variables. Assuming that the variable x is generated by the latent variable z, the generating process is as follows: (1) The variable z is generated, which is denoted as p θ (z); (2) The variable x is generated, which is denoted as p θ (x | z). In addition, q φ (z | x) is used to approximate a posteriori probability p θ (z | x) which is difficult to calculate. The whole process is shown in Figure 11.
The edge likelihood of X is obtained by summing the edge likelihood of each independent data. Namely: Since q φ (z | x) approximates p θ (z | x), then: Figure 11. The generation process of variable x.
The first term on the right side of the equation is the Kullback-Leibler (KL) divergence of approximate from the true posterior. Since the KL divergence is non-negative, the second item on the right side of the equation is called the (variational) lower bound on the marginal likelihood of data-point i, and can be written as: Let the latent variable z be the normal Gaussian distribution and p θ (x | z) be the Gaussian distribution. So formula (3) can be written as follows: The reparameterization is z (i,l) = µ (i) + σ (i) (l) . Where (l) is an auxiliary noise variable, and (l) ∼ N(0, I)). The objective of optimization is to make logp θ x (i) as large as possible. We use a neural network for the probabilistic encoder q φ (z | x) and where the parameters φ and θ are optimized jointly with the auto-encoding variational Bayes algorithm. The variational approximation distribution can be regarded as an encoder, which maps observable variables to latent variables. The generated model can be regarded as a decoder, which maps latent variables into observable variables. The structure of VAE is shown in Figure 12. q x z p The traditional VAE algorithm assumes that the posterior distribution p θ (z | x) of the latent variable z satisfies a single Gaussian distribution, which is easy to cause the low-dimensional representation is too simple and cannot fit the space of the latent variable well, leading to the low accuracy of the generated model. Based on this, the Gaussian mixture model is used to fit the latent variable space.

Gaussian Mixture Model
The Gaussian model is a commonly used variable distribution model. The probability density function of one-dimensional Gaussian distribution is as follows: In the formula, x is a variable, µ and σ 2 are the mean and variance of Gaussian distribution, respectively.The essence of the Gaussian mixture model is to fuse several single Gaussian models to make the model more complex. The one-dimensional Gaussian mixture distribution is expressed as follow: In the formula, K is the number of single Gaussian distribution in the Gaussian mixture model. ω k is the weight of the kth Gaussian distribution, and 0 ≤ ω k ≤ 1, ∑ K k=1 ω k = 1. Theoretically, if there are enough single Gaussian models fused by the Gaussian mixture model and the weight design between them is reasonable enough, the Gaussian mixture model can fit arbitrarily distributed samples [34].
Since we cannot directly observe which distribution the data x comes from, we introduce a latent variable z k represents the probability that data x comes from the kth Gaussian distribution. Assuming that z k is independent and identically distributed, then The probability form of conditional probability based on latent variable z is as follows: Thus, the form of p (x) is as follows:

Improved Variational Auto-Encoder
In this paper, the measurement error data set is used as the input of the VAE algorithm. The measurement error corresponding to each working range follows a single Gaussian distribution. Therefore, the Gaussian mixture model is used to describe the overall distribution of all measurement error data. The output of Decoder are the parameters of the Gaussian mixture model of measurement error data.
The latent variable z is assumed to have a normal Gaussian distribution in the traditional VAE algorithm. The Variational auto-encoder based on Gaussian mixtrue model (VAE-GMM) algorithm assumes that the hidden variable z is the Gaussian mixture distribution. Namely: Since the coefficients of the GMM cannot be trained and updated as trainable parameters by back propagation algorithm, the coefficients of the GMM are super-parameters, that is, α i = π i = 1 N . Then: Due to the true distribution of measurement error data p data(w) is unknown, we use distribution p model(w) to approximate it.
where ε i and σ 2 wi are the mean and variance of the i th Gaussian distribution of the Gaussian mixture model optimized by Decoder respectively, that is, w i = ε i + i , and i ∼ N 0, σ 2 wi where ε i is the system error, σ 2 wi is the variance of random error. The coefficient of i th Gaussian distribution is P(p(dis i ) = 1) = π i . It can be understood that when the working range of acoustic sensor is dis i , the measurement error follows a single Gaussian distribution, that is N w | ε i , σ 2 wi , namely the system error is ε i , the variance of random error is σ 2 wi . The structure of VAE-GMM is shown in Figure 13. The Support Vector Regression (SVR) algorithm is used to divide the overall distribution into multiple Gaussian distributions and its corresponding working range, and fit the relationship between them.

Support Vector Regression
When support vector machine (SVM) is used for classification, its basic idea is to find an optimal classification surface to separate the two classes of samples. When support vector machine is used for regression fitting analysis, its basic idea is to find an optimal classification surface so that all training samples have the least error from the optimal classification surface [35]. Support vector regression (SVR) is also divided into linear regression and non-linear regression. This paper adopts non-linear regression. Let us consider a training set S = (x 1 , y 1 ) , · · · , (x l , y l ) |x i ∈ R n , y j ∈ R, (i, j = 1 · · · l) consists of l training samples. For the sample set S which cannot be linearly separated in the original space R n , using a non-linear mapping Φ to map data S into a high-dimensional feature space H, to make Φ (S) has good linear regression characteristics in feature space H. Then linear regression is carried out in feature space H, and finally returns to the original space R n . The specific implementation steps of non-linear regression are as follows: 1. Looking for a kernel function K (s, t), makes K x i , x j = Φ (x i ) , Φ x j . The commonly used kernels are radial basis function, polynomial function, spline curve function, sigmoid function and so on. 2. Finding the solution of the optimization problem.
In the formula, ε specifies the error requirement of the regression function. The smaller ε is, the smaller the error of the regression function is. C is a penalty factor. The larger C is, the greater the penalty for training error greater than ε. α i , α * i are the optimal solution of the equation.
4. The regression function is: The structure of SVR is similar to that of a neural network, as shown in the Figure 14. The output is a linear combination of intermediate nodes, each of which corresponds to a support vector. The training samples set in this paper are as follows: Where d i represents the working distance of an acoustic sensor. ε i represents the system error of the acoustic sensor at this working distance. σ 2 wi represents the variance of random error at this working distance. In the regression algorithm of fitting working range and system error, the radial basis function is used as the kernel function. In the regression algorithm of fitting working range and the variance of random error, the B-spline curve function is used as the kernel function. With 1000 training samples, the fitting results are shown in Figures 15 and 16.

Simultaneous Localization of AUV and Fixed Docking Station
The states of AUV and DS are estimated simultaneously by the FastSLAM2.0 algorithm in the process of homing an AUV to the fixed DS. FastSLAM algorithm assumes that the position and attitude of AUV follow the known hypothesis distribution, and uses particle filter to estimate the posterior probability distribution of the position and attitude of AUV, where each particle represents a possible position and attitude. An EKF filter is used to estimate the posterior probability distribution of features in the environment. The joint status of SLAM is expressed as follows: where s k = s 1 , s 2 , · · · s k represents the state sequence of AUV, and s k = [x k , y k , Ψ k ] T are the position and heading angle information of AUV, respectively. Θ m = θ 1 , θ 2 , · · · , θ m are the features in the environment. z k = z 1 , z 2 , · · · , z k represents the measurement sequence. u k = u 1 , u 2 , · · · , u k represents the input sequence. FastSLAM algorithm includes the following four basic iteration steps: The formulas corresponding to the above algorithm are sorted out in Table 2: Table 2. FastSLAM2.0.
Step Formula The motion model of AUV and measurement model are shown in Equations (2) and (3).

Simultaneous Localization of AUV and Mobile Docking Station
Assuming that M-DS does a uniform linear motion and has good communication with AUV. The SLAM with mobile object (SLAMMO) problem based on sampling mechanism is to calculate the joint posterior probability distribution of AUV and M-DS.
where x k−1 and x k are the states of AUV at 1 : k − 1 time and 1 : k time, respectively. Θ k−1 and Θ k are the states of M-DS at 1 : k − 1 time and 1 : k time, respectively. z k represents the measurement at kth time. u x,k represent the control input of AUV at 1 : k time. u θ,k represent the control input of M-DS at 1 : k time. Firstly, the posterior probability distribution of AUV is calculated, and the proposed distribution is expanded as follows by Bayesian principle: According to the motion model of AUV, the above formula is simplied to: Using the theory of total probability and convolution theorem, the approximate form of the proposed distribution is obtained.
The distribution represented by the formula (23) is Gaussian, and its mean and variance are given by the minimum value and curve of the formula (25). By calculating the first order and second order of differential of y k to x k , the fllowing formulas are obtained.
where µ x,k is the estimated state of AUV at kth time. Σ x,k is the estimated variance of AUV at kth time. Based on the sampling mechanism, the posterior probability of M-DS is calculated by the fllowing formula: Using the same method as calculating the sate of AUV, the estimated state and variance of M-DS are as follows: where µ Θ,k is the estimated state of M-DS. Σ Θ,k is the estimated variance of M-DS. Due to the deviation between the proposed distribution and the expected distribution, this difference is corrected by importance sampling. In the state probability distribution of AUV, the importance weight of each particle is calculated as follows: In the state probability distribution of M-DS, the importance weight of each particle is calculated as follows: After calculating the weight of each particle, the random resampling is used to resample. The purpose of resampling is to retain weighted particles in order to reduce the degree of particle degradation.

Analysis of Experimental Data
Our recovery system was tested on a lake at the xinanjiang experimental field of the institute of acoustics, Chines Academy of Sciences in December 2018. AUV sensed the position of DS by using the ultra short baseline (USBL) system and camera. The intersection of the estension line of AUV initial heading and the center of DS is the initial target point. AUV reached the target point by line-of-sight (LOS) guidance, and tracked the center line of DS. More detailed descriptions are given in Reference [36]. The experimental results are shown in Figure 17.  We use these experimental data to verify the effectiveness of the proposed algorithm. Firstly, the method 1 and method 2 are applied to collect the measurement error data of acoustic sensor. Then the VAE-GMM algorithm is applied to estimate the overall distribution of measurement error data. Finally, the SVR algorithm is used to fit the relationship between the ststistical characteristics of measurement error data and its coressponding working distance. After data pre-processing, SLAM algorithm is applied to estimate the positions of AUV and DS. The actual measurement error in the experiment and the corresponding to estimated statistical characteristics are partly shown in Table 3. The data analysis results are shown in Figures 18-21. Figures 18 and 20 show the deviation between the estimated position and true position of AUV sensed by IMU. We can know that the position of AUV estimated by SLAM algorithm is close to the position sensed by IMU. It is applicable if without IMU. In Figures 19 and 21, the blue line shows the deviation between the position of DS estimated by SLAM and true position of DS sensed by IMU. The red line shows the deviation between the positon of DS directly estimated by the measurement of sensor and true position. We can know that the localization precision of the SLAM algorithm is higher than that of the original data.
When the initial distance between AUV and DS is larger than the maximum operating distance corresponding to measurement error data, the above method can still be used. This is because the SVR method has the ability of prediction after training. In the experiment of AUV homing to the M-DS mentioned above, the initial distance between AUV and M-DS was about 207 m. Assuming that we only use measurement error data with the maximum working distance less than 150 m, the experimental data analysis results are as follows:   In Figure 22, the blue line and red line respectively show the deviation between the estimated position and the true position of AUV sensed by IMU with the the maximum working range corresponding to measurement error data less than 150 m and less than 250 m. In the Figure 23, the blue line and red line respectively show the the deviation between the estimated position and true position of M-DS sensed by IMU with the maximum working range corresponding to the measurement error data less than 150 m and less than 250 m. The black line shows the position of M-DS directly estimated by the measurement of the sensor. We can know that the positioning accuracy of the SLAM algorithm with the maximum working range corresponding to measurement error data less than 150 m is lower than that with the maximum working range corresponding to measurement error data less than 250 m, but higher than that of the original algorithm. So even if the initial distance between AUV and DS exceeds the maximum working distance corresponding to the measurement error data, the algorithm proposed in this paper is still effective.

Conclusions
We propose a data pre-processing method for SLAM with unknown measurement error. The measurement error data collected by method 1 and method 2 can reflect the statistical characeristics of the measurement error of an acoustic sensor. The measurement error corresponding to each working range follows a single Gaussian distribution, and the measurement errors corresponding to multiple working distances follow the Gaussian mixture distribution as a whole. The overall distribution of measurement error data is non-Gaussian, and its statistical characeristics can be estimated by a VAE-GMM algorithm. The SVR algorithm is applied to divide the overall distribution into multiple single Gaussian distributions and their corresponding working range, and fit the relationships between them. After data pre-processing, the SLAM algorithm is applied to estimate the positions of AUV and F-DS and M-DS. The analysis of experimental data shows that the position estimated by the SLAM algorithm is closer to the true position than that directly estimated by the sensor measurement.
Author Contributions: L.D. conceived the method and wrote the manuscript draft; H.X. and X.F. helped to modify it; X.H. and C.Y. provided experimental data.
Funding: This research was funded by Joint fund for equipment pre-research of the Chinese academy of sciences with Grant No. 6141A01060101.

Conflicts of Interest:
The authors declare no conclict of interest.