An Adaptive Filtering Method for Cooperative Localization in Leader–Follower AUVs

In the complex and variable marine environment, the navigation and localization of autonomous underwater vehicles (AUVs) are very important and challenging. When the conventional Kalman filter (KF) is applied to the cooperative localization of leader–follower AUVs, the outliers in the sensor observations will have a substantial adverse effect on the localization accuracy of the AUVs. Meanwhile, inaccurate noise covariance matrices may result in significant estimation errors. In this paper, we proposed an improved Sage–Husa adaptive extended Kalman filter (improved SHAEKF) for the cooperative localization of multi-AUVs. Firstly, the measurement anomalies were evaluated by calculating the Chi-square test statistics based on the innovation. The detection threshold was determined according to the confidence level of the Chi-square test, and the Chi-square test statistics exceeding the threshold were regarded as measurement abnormalities. When measurement anomalies occurred, the Sage–Husa adaptive extended Kalman filter algorithm was improved by suboptimal maximum a posterior estimation using weighted exponential fading memory, and the measurement noise covariance matrix was adjusted online. The numerical simulation of leader–follower multi-AUV cooperative localization verified the effectiveness of the improved SHAEKF and demonstrated that the average root mean square and the average standard deviation of the localization errors based on the improved SHAEKF were significantly reduced in the case of the presence of measurement abnormalities.


Introduction
The Global Navigation Satellite System (GNSS) is often unavailable in deep water environments, and inertial/dead reckoning (DR) navigation technology is usually considered. The Doppler velocity log (DVL) can measure the velocity of autonomous underwater vehicle (AUVs) and it can carry out autonomous navigation and localization without installing additional sensors and external reference information when the DVL is combined with the Inertial Navigation System (INS). In a multi-AUV formation, there will be localization error growth over time for the AUVs, which are only equipped with inertial/DR equipment. We can improve the localization accuracy of the whole system by improving the localization accuracy of each AUV, but there is a great sacrifice in cost. Therefore, researchers have proposed leader-follower multi-AUV cooperative localization methods [1][2][3]. AUVs transmit information through sensors, such as underwater acoustic communication equipment. When conducting cooperative localization, various state estimation technologies need to be introduced to estimate the position of the AUV. A Kalman filter (KF) is the best Bayesian estimator for linear systems with Gaussian uncertainty [4]. However, the model of a leader-follower multi-AUV cooperative localization system is often nonlinear. Therefore, an extended Kalman filter (EKF) or unscented Kalman filter (UKF) is usually used for state estimation. In the framework of Kalman filters, these methods require that the statistical characteristics of process noise and measurement noise are known and fixed during cooperative localization. However, in practical applications, due to the low precision of low-cost inertial/DR equipment carried on an AUV, these methods are easily affected by many uncertain factors such as external environment interference, carrier maneuver changes, internal instrument failure, and so on. At present, in order to overcome the uncertainty of the filter in AUV navigation and localization, many researchers have applied adaptive interactive multi-models, fuzzy logic, and other methods [5][6][7][8][9][10]. In Refs. [5,6], the adaptive filtering algorithm was applied to an AUV. When the system noise matrix remained positive or semi positive, the latest measurement data were used to adaptively estimate the changing process noise or measurement noise covariance matrix. Refs. [7,8] introduced interactive multi-model methods in order to solve the problem of state estimation when the statistical characteristics of measurement noise were unknown or easy to change under the complex working environment of an AUV, and the model probability obtained in the process of state estimation was used for decision-making, so as to obtain the expected model closer to the real mode of the system for state estimation. In refs. [9,10], researchers used fuzzy logic technology to adapt to sensor noise changes or communication data loss, enhanced the fault detection and signal recovery algorithm of the navigation system reliability, adjusted the initial noise statistical assumption of the EKF, and maintained the stability and performance of the filter. In the aforementioned research, the improved filter algorithm was mostly used for navigation and localization of a single AUV. This study was concentrated on the improved filter algorithm for the navigation and localization of leader-follower AUVs. In this paper, an improved Sage-Husa adaptive extended Kalman filter (improved SHAEKF) was proposed to improve the adaptability of the filter and the accuracy of cooperative localization.

State Model
The navigation motion of AUV formation in a GNSS-denied environment is a complex three-dimensional process. If the formation adopts the leader-follower multi-AUV structure, the system includes multiple leader AUVs and follower AUVs; a leader AUV and a follower AUV exchange and share information through underwater acoustic communication equipment [11]. After the follower AUV obtains the cooperative localization information, the cooperative localization filtering algorithm is used to correct position information. In order to facilitate the study of cooperative localization algorithms with measurement anomalies, the state variables of a leader AUV (the lth leader AUV) and a follower AUV (the jth follower AUV) are analyzed by a simplified level equation under reasonable assumptions, and the non-rotating east, north, and up (ENU) geographic reference system O g − X g Y g Z g is used to analyze the state variables, as shown in Figure 1a.
Assuming that the depth of an AUV can be measured directly and accurately by depth gauge, the discrete-time nonlinear kinematics equations of an AUV are as follows: where x, y respectively represent the coordinates of the longitude and latitude position of the AUV converted into the geographical coordinate system, and ψ, V, ∆T respectively represent the heading angle, forward synthetic motion velocity, and sampling period. Moreover, in this paper, the quantity with superscript " " represents the observed value of the correlation quantity, and the quantity with superscript " " represents the predicted value of the correlation quantity. Figure 1. Schematic diagram of a multi-AUV cooperative localization system. (a) Description of the projection of the leader AUV in the O g − X g Y g Z g coordinate system with respect to the O g − X g Y g horizontal plane in the O g − X g Y g Z g coordinate system; (b) description of the relative position of the leader AUV and the follower AUV in the O g − X g Y g horizontal plane.
The state vector X and the control input vector u of the jth follower AUV at time k are defined as follows: where the right superscript "T " represents the transpose mode of the vector matrix.
The observed values of the control input are as follows: where w j k,ψ , w j k,V represent the observation errors of the heading angle and the forward synthetic motion velocity, respectively. It is assuming that w where δ pq represents the Dirac delta function, δ pq = 1 i f p = q 0 otherwise . The process noise covariance matrix is Q j k and the state model equation of the jth follower AUV can be expressed as follows: where f (·) is a nonlinear function.

Measurement Model
In this study, the cooperative localization information mainly included the relative distance and relative direction angle between AUVs and the position coordinates of the leader AUV. The follower AUV established the measurement equation by obtaining the cooperative localization information shared by the leader AUV. Figure 1b depicts the geometric relationship between the relative distance r lj and the relative direction angle α lj of the lth leader AUV and the jth follower AUV. It was assumed that the forward movement direction of the follower AUV was consistent with the positive direction of the Y b -axis of its own carrier coordinate system. From the positive direction line of the Y b -axis, the relative direction angle was positive when moving clockwise to the relative direction angle line, otherwise, it was negative. Assuming that the position coordinates of the lth leader AUV at time k is [ x L l k y L l k ] T , the geometric relation equations at time k are as follows: where r lj k,x = x L l k − x j k , r lj k,y = y L l k − y j k respectively represent the decomposition of the true relative distance r lj of the AUV in the direction of the X g -axis and the Y g -axis of the geographical coordinate system. The observed values of the relative distance and relative direction angle between the jth follower AUV and the lth leader AUV are as follows: where v lj k,r and v lj k,α respectively represent the observation errors of the relative distance and relative direction angle between the jth follower AUV and the lth leader AUV. It is assumed that v If (6) is substituted into (7) and we make cos v k,α = 0, the following formulas are obtained by simplification: T is the observation noise of the measurement model based on the observation of relative distance and relative direction angle, the covariance matrix is R . The nonlinear measurement model equation of the jth follower AUV is established as follows: where h(·) is a nonlinear function.

Cooperative Localization Algorithm Based on an EKF
From (4) and (9), we obtain the discrete-time nonlinear state equation and measurement equation of the jth follower AUV as follows: The EKF algorithm is used for state estimation. The steps in a single filtering cycle are as follows.
Step 1 (state prior estimation): Step 2 (the innovation sequence update): Step 3 (error covariance updating of state a priori estimation): Step 4 (filter gain update): where H j k+1 is the Jacobian matrix of h with respect to X j k+1|k .
Step 5 (state a posteriori estimation update): Step 6 (error covariance of state a posteriori estimation): where I is the identity matrix.

Adaptive Cooperative Localization Algorithm
Ref. [13] analyzes the position error covariance of a group of agents equipped with an inertial measurement unit (IMU) and relative distance and azimuth sensors. Researchers regard all the agents in the group as a unified system. Agents exchange their information, including measured and estimated location data. Therefore, the exchange of each external perceptual measurement results in an improvement in the overall position error estimation. The analysis showed that the system position error covariance increased with the cube of time. Accordingly, as the number of agents increased, the growth rate of covariance decreased. However, due to the different construction methods of the system state equation and the measurement equation, the conclusions obtained from the analysis were also different. For example, in the agent equipped with an IMU, the covariance growth of position error did not depend on the path. However, the position error covariance of an agent equipped with an encoder depended on the path and changed with the path direction.
Refs. [12,14] studied the transfer equation of the overall localization uncertainty with respect to the relative position measurement error. By solving the algebraic Riccati equation of the evolution of the system localization errors with time, the analytical expression of the upper bound of the expected positioning uncertainty was determined. The analysis of ref. [12] shows that the upper bound of the variance of the localization errors in the steady state depended on the measurement accuracy of the AUV's speed, heading angle, and relative position, and was independent of the initial filter variance of the system.
All the above studies were under the steady state of the system. We know that the noise error of an AUV sensor is related to the process noise error and measurement noise error in the cooperative localization algorithm based on an EKF, but there are differences between process noise error and measurement noise error. Through error uncertainty propagation analysis [12][13][14], it was found that the covariance matrix of system process noise and measurement noise were time-varying.The process noise covariance matrix Q j k is related to the observation value of heading angle, which is mainly determined by the internal mechanism of a system and is relatively stable. The measurement noise covariance matrix R lj k+1 is not only related to the error covariance matrix R lj k,Z lj of the measurement model, but also to the observation value of heading angle and the noise of the external information, which is easy to change, and there is a certain unpredictability. In the model of the leader-follower multi-AUV cooperative localization system studied in this paper, we introduced and improved the SHAEKF to help solve the adverse impact of the follower AUV measurement anomalies on localization. A measurement anomaly might occur in the prior characteristics of measurement noise based on relative distance and relative heading angle, or when the leader AUV transmits position data to the follower AUV. This problem has received less attention in previous research on cooperative localization algorithms.

Sage-Husa Adaptive Kalman Filter Algorithm
The Sage-Husa adaptive filtering algorithm is an improvement based on the classical Kalman filtering algorithm. The Sage-Husa filter adaptively estimates system process noise and measurement noise online, including the mean vector q k , the covariance matrix Q k of the process noise, the mean vector r k , and the covariance matrix R k of the measurement noise. The estimators of the Sage-Husa adaptive filtering algorithm can be expressed as follows [6,15]: Since the innovation vector affects the calculation of Q k and R k at the same time, this easily led to filter divergence. Therefore, the noise estimator could not estimate the statistical properties of the process noise and measurement noise at the same time. Moreover, there was a negative sign in the above formula and the positive definite of Q k and R k could not be fully guaranteed.
We focused on finding a method to detect an anomaly and re-estimate the statistical characteristics of the measurement noise when the measurement noise was abnormal, which is discussed next in this paper.

Measurement Anomaly Detection
This innovation can reflect the relationship between the observed value and the estimated value of a measurement model, so it is often used to measure the filtering performance of a Kalman filter [16,17]. According to Kalman filter theory, if the discrete-time nonlinear state space model of process and measurement noise is established under the assumption of Gaussian distribution and the measurement is not subject to abnormal inter-ference, the innovation sequence obeys the standard Gaussian distribution. The expectation and covariance of the innovation have the following characteristics: The probability density function of the innovation can be written in the following form: where m is the degree of freedom, and P s lj k+1 is the determinant of P s lj k+1 .
If there are some observed outliers in the measurement model, or the measurement noise is affected by other noises and no longer conforms to the Gaussian distribution, (18) and (19) will no longer be tenable. It can be considered that there are some violations of assumptions or some modeling errors. Specifically, hypothesis testing can be performed to detect measurement anomalies. The purpose of hypothesis testing is to check whether the actual measurement is compatible with the hypothetical model or, in other words, the zero hypothesis. According to the orthogonality principle of the innovation in Kalman filtering, the square of Mahalanobis distance based on innovation obeys the Chi-square distribution. We can took (19) as the relevant zero distribution that did hold under the hypothetical model, and then constructed the test statistics based on the following innovation: (20) ξ k+1 follows the Chi-square distribution with degree of freedom m if the assumption holds. The significance level is set as γ, which indicates the probability threshold that the null hypothesis below this threshold will be rejected. The critical value of the corresponding Chi-square test is χ 2 γ (m). When the innovation sequence is calculated with the actual observation value, the actual judgment index ξ k+1 is obtained. When ξ k+1 > χ 2 γ (m) rejects the original hypothesis, it is considered that the measurement is abnormal. For example, from the Chi-square distribution table, we obtained that the probability of {χ 2 γ (m = 3) > 11.345} was only 1%, i.e., γ = 0.01. The problem of measurement anomaly detection is expressed as:

Adjustment of Measurement Noise Covariance Matrix
It is generally believed that we should pay attention to the position of recent measurements in the current filtering in order to prevent filtering divergence; hence, we need to pay special attention to the innovation sequence [18][19][20][21].
In the filtering process of the improved SHAEKF, we used (21) to judge the filtering state. If (21) was not tenable, this indicated that there was measurement abnormality and it was necessary to re-estimate the measurement noise parameters. If (21) held, there was no need to re estimate the measured noise parameters. In the conventional Sage-Husa adaptive filtering algorithm, the measurement noise parameters were adaptively estimated by the equal-weighted time average: There is a negative sign in the above formula. In order to prevent the measurement noise estimation from losing its positive definiteness, the measurement noise covariance matrix R lj k+1 is estimated based on the maximum a posteriori (MAP) estimation criterion: Replace X j i|k+1 with X j i|i to obtain the following approximation: If (24) is substituted into (23), we get the following equation: It is easy to verify that the mean value of R lj k+1 is as follows: The suboptimal MAP estimates of the measurement noise covariance matrix are obtained from (25) and (26): Considering that K j k+1 and P j k+1 must be obtained after R lj k+1 and become stable during the filtering process, K j k and P j k+1|k can be used here to approximate K j k+1 and P j k+1 respectively.
From (27), the estimation of R lj k+1 could be improved by using the exponential fading memory weighted average recurrence method in the following equation: where, b is the fading memory factor, 0.95 < b < 0.99. It could be seen that when b was larger, R lj k accounted for a larger proportion of the estimated value of R lj k+1 . After the improvement, if the estimated value of the last filter gain was less than 1, the estimated value of R lj k+1 could be guaranteed to be positive and definite, which suppressed the possibility of filter divergence.
In addition, considering that the state posterior estimation P j k+1 = (I − K j k+1 H j k+1 )P j k+1|k , in general the EKF algorithm could only be established when the optimal filter gain was used. When the filter was not always stable or when a non-optimal filter gain was used, an un-simplified state posterior estimation error covariance equation was required: Figure 2 describes the flow chart of the improved SHAEKF, and evaluates the measurement anomaly using the calculated the Chi-square test statistic based on innovation. The detection threshold was determined according to the confidence of the Chi-square test, and the Chi-square test statistics exceeding the threshold were regarded as measurement anomalies. When an observation anomaly occurred, the suboptimal MAP estimation weighted by the exponential decay memory was used to adjust the observation noise covariance matrix online. When no observation anomaly occurred, the measurement noise covariance at the previous time was maintained. It could then adjust the measurement noise covariance adaptively and ensure the stability of the filter.

Monte Carlo Simulation
In order to verify the feasibility of the improved SHAEKF algorithm, a leader AUV and a follower AUV were simulated in this study. Monte Carlo method was used in the simulation, and the simulation parameters were set as shown in Table 1. The initial error covariance was set to Q j 0 = diag( 0.025 2 0.5 2 ), R j 0 = diag( 10 2 10 2 ), P j 0 = 100 · I 4 , where diag (·) represents a diagonal matrix, and I 4 is the identity matrix of order 4. Table 1. Simulation parameter setting.

Simulation Parameters Parameter Values
Simulation time , s 1200 Sampling frequency, Hz ±0.025 Random error of heading angle, (rad · s −1 ) We designed the AUVs to move along a straight line or a curve at a uniform constant speed, the actual trajectory and the DR calculation results of each AUV in a simulation test of the whole simulation are shown in Figure 3. As shown in Figure 4, we study the filtering effects of the EKF and the improved SHAEKF with or without measurement anomaly sequence under the above conditions. As shown in Figure 5a,b, in the steady state, we assumed that the observed values of the measurement model were continuously disturbed by a normally distributed random sequence with mean value of 0 and standard deviation of 10. Under the above simulation conditions, we studied the filtering effects of the EKF and the improved SHAEKF with and without a measurement anomaly sequence. As shown in Figure 5c,d, the normally distributed random anomaly sequence with a mean value of 0 and standard deviation of 150 was added to the observation values of the measurement models in the time period of 200-1200 s. Then, the EKF and the improved SHAEKF algorithms were used for the simulation by calculating the Chi-square test statistics based on innovation as the judgment index. Considering that the measurement model in this study involved relative distance, relative direction angle, and heading angle, the degree of freedom was taken to be m = 3 and the confidence level was set at γ = 0.01. The results of evaluating the measurement anomalies are shown in Figure 6. As shown in Figures 7 and 8, the localization results of DR, EKF, and improved SHAEKF could be obtained respectively, including the root mean square and the standard deviation of the localization errors. The DR results of the follower AUV were omitted in Figures 7b,d,f,h  and 8b,d,f,h in order to more intuitively compare the localization errors of the EKF and the improved SHAEKF.  Figure 3 shows that the localization errors of the follower AUV became larger and larger with the passage of time only based on DR, and the localization errors of the follower AUV were much larger than that of the leader AUV.

Comparison and Analysis of Simulation
As can be seen from Figure 4, after obtaining the relative distance, relative direction angle, and high-precision DR localization data of the leader AUV, the follower AUV could improve its localization accuracy by using the EKF or the improved SHAEKF cooperative localization algorithm. It can be seen from Figure 4b,d, that when there were measurement anomalies, the filtering stability of the EKF was not as stable as that of the improved SHAEKF.  It can be seen from Figures 5 and 6 that after adding the measurement anomaly sequence, the measurement model vector components of the follower AUV were disturbed and the measurement anomaly could be evaluated by calculating the Chi-square test statistics based on innovation. In order to evaluate the performance of the improved algorithm more effectively, we introduced the average root mean square of localization error (ARMSE) and the average standard deviation of localization error (ASDE). ARMSE and ASDE can be defined as follows: where n T and n MC represent the sampling times and the Monte Carlo simulation times, respectively. ARMSE and ASDE respectively reflect the localization accuracy and the discreteness of the localization error dataset (i.e., the stability of the localization algorithm).  As can be seen from Figure 7, both the conventional EKF and improved SHAEKF work well if no measurement abnormality occurred. In this case, compared with the conventional EKF, the performance of the improved SHAEKF was not reduced, which meant that the contribution of these good measurements to the cooperative localization effect was not reduced, indicating the feasibility of the improved SHAEKF. It can be seen from Figure 8, when there were outliers in the measurement, the root mean square value and standard deviation of the localization errors obtained by EKF increased, and the root mean square value and standard deviation of the localization errors obtained by the improved SHAEKF were reduced to a certain extent. The improved SHAEKF could effectively resist the influence of outliers and had better performance than the EKF in localization accuracy, which indicated that the improved SHAEKF was more stable than the EKF. The above results could also be obtained by analyzing the data in Tables 2 and 3.

Conclusions
The conventional EKF was only the best estimator when some preconditions were true, for example, when the Gaussian distribution measurement noise had a completely known mean and covariance. However, in the practical application of multi-AUV cooperative localization, the actual observation might be vulnerable to outliers. Therefore, a leader-follower multi-AUV adaptive cooperative localization filtering method based on innovation was proposed in this paper. Simulations showed that the improved SHAEKF could effectively resist the influence of measurement anomalies. Although the improved SHAEKF could ensure the filtering fault tolerance and localization accuracy of a leader AUV and a follower AUV in the case of abnormal measurement, we did not further study the influence of abnormal process noise on cooperative localization. In the future, we will further study the cooperative navigation and localization of multiple leader AUVs and multiple follower AUVs and expand from a two-dimensional plane to a three-dimensional space. All AUVs in the group were regarded as a unified system to improve the navigation and localization system model. AUVs exchanged their information, including measured and estimated position data, we then proposed an interactive multi-model fusion estimation method for future work to correctly estimate some parameters that could not be accurately obtained in the AUV dynamic model, including multi-observation information and multi-motion model information, so as to increase the robustness of collaborative navigation and localization.