# Dual-Model Reverse CKF Algorithm in Cooperative Navigation for USV

Authors: Bo Xu; Yong ping Xiao; Wei Gao; Yong gang Zhang; Ya Long Liu; Yang Liu
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

As one of the most promising research directions, cooperative location with high precision and low-cost IMU is becoming an emerging research topic in many positioning fields. Low-cost MEMS/DVL is a preferred solution for dead-reckoning in multi-USV cooperative network. However, large misalignment angles and large gyro drift coexist in low-cost MEMS that leads to the poor observability. Based on cubature Kalman filter (CKF) algorithm that has access to high accuracy and relative small computation, dual-model filtering scheme is proposed. It divides the whole process into two subsections that cut off the coupling relations and improve the observability of MEMS errors: it first estimates large misalignment angle and then estimates the gyro drift. Furthermore, to improve the convergence speed of large misalignment angle estimated in the first subsection, “time reversion” concept is introduced. It uses a short period time to forward and backward several times to improve convergence speed effectively. Finally, simulation analysis and experimental verification is conducted. Simulation and experimental results show that the algorithm can effectively improve the cooperative navigation performance.


---

## Body

## 1. Introduction

In recent years, cooperative navigation technology has been gradually expanded in satellite [1], multirobot [2–4], intelligent underwater vehicles [5], and many other fields. Research team named Stergios I. Roumeliotis from America has conducted a series of theoretical studies for colocation of robots, including observability analysis, estimation error analysis, consistency analysis, and formation optimization [6]. As one of the most promising research directions in navigation and positioning, cooperative navigation is attracting more and more attention in colocation industry and academia.Unmanned surface vehicle (USV) refers to running crafts on surface of water that can be remotely operated by a person or operate independently. It has characteristic of small size, high mobility, and low-cost. Multi-USV formation can take advantage of network using coordination and communication technology. The US Navy issued “the navy unmanned surface craft master plan” [7] that shows the significance of cooperative combat techniques of multiple USV applied in the military. Cooperative navigation of multi-USV can be commonly divided into two types according to configuration: one is master-slave (leader-follower) structure (some with high-precision and others with low precision); the other is parallel structure (navigation accuracy of each USV is equivalent). Compared to master-slave structure, parallel cooperative navigation cannot reconcile the contradiction between high-precision and low-cost fundamentally, so master-slave USV cooperative navigation is the mainstream currently.With the rapid development of computer and driven by demand of high-performance filtering in engineering, nonlinear filtering theory has been greatly developed [8–14]. EKF, UKF, and PF are three commonly used nonlinear filtering methods in practical applications. EKF has high calculating efficiency, but the precision of EKF is limited because of the Taylor expansion that neglects high-order terms. UKF and PF have good performance in some applications; however, covariance calculated in the UKF filtering process might be nonpositive definite that leads to numerical instability and filtering divergence, while PF is easy to fall into the particles degradation problem and heavy computational complexity for high-dimensional applications [12]. CKF can avoid linearization of the nonlinear system by using cubature point sets to approximate the mean and variance. Because of its high accuracy and low calculation load, CKF [14] has become a hotspot in research. In this paper, CKF is introduced in the cooperative navigation with nonlinear model.As a branch of inertial field, inertial systems based on microelectromechanical system (MEMS) gained rapid development in recent years. Because of its low-cost, small size, light weight, and high reliability, there have been more and more applications in low-cost inertial systems. It has become a preferred solution for the inertial navigation device on follower boats in multi-USV cooperative network. However, due to the limitations of accuracy, it is necessary to construct reasonable error model of MEMS based on which estimation and compensation can be preferred [15]. Commonly used error compensation methods include RBF neural network analysis, ARMA model used in the timing model of gyro random error, and piecewise interpolation compensation of scale factor [16]. However, coexistence of larger initial misalignment and larger gyro drift is still unavoidable in practical use that leads to poor observability in cooperative navigation. Based on different character of initial misalignment and gyro drift, estimation with subsection idea is proposed to cut off the coupling relations.In addition, how to obtain accurate estimation performance in short time is also a key technique. Li et al. proposed a new alignment method for AUV with backtracking framework. The fine alignment runs with the recorded data during the process of coarse alignment which effectively improves the convergence speed and improve navigation accuracy [17, 18]. Xixiang Liu also proposed gyrocompass alignment method based on reverse loop to solve compass alignment for SINS on moving base and it greatly reduces the alignment time. To improve the convergence speed of large initial misalignment, “time reversion” concept is introduced in this paper. As will be confirmed by simulation and experimental results, the convergence speed and estimation accuracy are improved greatly.The rest of this paper is organized as follows. Section2 presents cooperative navigation model of single master USV and the process of CKF. In Section 3, based on observability analysis of cooperative navigation, dual-model method and reverse filtering are proposed that improve the convergence speed and improve the accuracy of cooperative location. Simulation of the proposed algorithm is provided in Section 4. Aiming at validating the proposed algorithm further, the process and the results of experiment in Tai Lake are given in Section 5; Section 6 gives data analysis and discussions. In Section 7, a summary is provided, and future research directions are discussed.


## 2. Cooperative Navigation Modeling and Problem Statement

### 2.1. Fundamental Theory of Cooperative Navigation Modeling

Master USV and slave USV are equipped with hydroacoustic communication modem, respectively, at their bottoms to achieve information dissemination and distance measurements between USVs. The masters accomplish navigation using their high-precision device (usually strap-down inertial navigation system combined with GPS) while the slaves use their own equipped Doppler velocity log (DVL) and attitude and heading reference system (AHRS) based on MEMS to proceed dead-reckoning.The principle of multi-USV cooperative navigation system is as follows. Firstly master USV sends out fixed frequency acoustic pulse signal according to preset time interval. With this signal, slave USV can calculate the relative distance according to underwater acoustic propagation delay. Then master USV sends out its position information by hydroacoustic modem in the form of broadcasting. After receiving this information, slave USV can conduct information fusion with dead-reckoning, position of master USV, and the relative distance and then correct its navigation result using filtering algorithms. Moreover, route planning of USV can be carried out to improve the observability and further improve the accuracy of the cooperative navigation system.Positioning errors will accumulate over time (as shown by the purple dotted ellipse shown in Figure1). Therefore, it cannot provide accurate positioning information for a long time. When the information sent by master USV is received, the slave USV can use cooperative navigation algorithm to correct the result of dead-reckoning and effectively reduce the positioning error (as shown by the green dotted ellipse shown in Figure 1). Through this cyclical filtering and correction, slave USV with low-precision navigation device can realize precise location for a long time [19, 20].Figure 1
Principle of cooperative navigation system.


### 2.2. System Model of Single Master USV

When the follower does not receive the information sent by leaders, it can only use its own DVL and low accuracy MEMS to conduct dead-reckoning, and the motion model can be expressed as follows:(1)
x
k
+
1
S
=
x
k
S
+
v
k
·
Δ
t
·
sin
α
k
,
y
k
+
1
S
=
y
k
S
+
v
k
·
Δ
t
·
cos
⁡
α
k
,
where x
k
S, y
k
S represent the position of follower USV at time k, v
k represents the speed, Δ
t means the time interval, and α
k means the heading angle at time k. However, the input model of system in practice is as follows:
(2)
v
k
=
1
+
δ
C
k
v
^
k
+
w
v
k
,
α
k
=
α
^
k
-
δ
α
k
+
w
α
k
,
wherein v
^
k is the velocity measurements and δ
C
k means the scale factor error. α
^
k is the heading angle measured by MEMS, δ
α
k
=
Δ
α
k
+
k
ε
Δ
t is error of α
^
k, Δ
α
k is the initial misalignment angle, and ε is gyro drift. Putting (2) into (1), the actual equation of states is as follows:
(3)
x
k
+
1
S
=
x
k
S
+
1
+
δ
C
k
v
^
k
+
w
v
k
·
Δ
t
·
sin
α
^
k
-
Δ
α
k
-
ε
k
·
k
·
Δ
t
+
w
α
k
,
y
k
+
1
S
=
y
k
S
+
1
+
δ
C
k
v
^
k
+
w
v
k
·
Δ
t
·
cos
⁡
α
^
k
-
Δ
α
k
-
ε
k
·
k
·
Δ
t
+
w
α
k
,
δ
C
k
+
1
=
δ
C
k
,
Δ
α
k
+
1
=
Δ
α
k
,
ε
k
+
1
=
ε
k
.Motion model can be simplified as(4)
X
k
+
1
=
X
k
+
T
u
k
,
w
k
in which X
k
=
x
k
,
y
k
,
Δ
α
k
,
ε
k
,
δ
C
k
T means the state vector of USV at time k, u
k
=
[
v
^
k
α
^
k
]
T, T
(
u
k
,
w
k
) stands for the nonlinear terms, the process noise is w
k
=
[
w
v
k
w
α
k
]
T, w
v
k
~
N
(
0
,
σ
v
k
2
), and w
α
k
~
N
(
0
,
σ
α
k
2
). The noise covariance matrix is
(5)
Q
k
=
E
w
k
w
k
T
=
σ
v
k
2
σ
α
k
2
.Systematic observation equation can be expressed as:(6)
Y
k
=
R
k
+
V
k
=
x
k
M
-
x
k
S
2
+
y
k
M
-
y
k
S
2
+
V
k
,
wherein Y
k stands for the measurement of distance R
k between follower and leader, (
x
k
M
,
y
k
M
) represents the position of leader USV at time k while x
k
S
,
y
k
S represents the position of follower, and V
k
~
N
(
0
,
σ
r
2
) is the measurement noise.


### 2.3. Cubature Kalman Filter (CKF)

Cubature Kalman filter (CKF) is proposed based on the spherical-radial cubature criterion [14]. The core of the method is that the mean and variance of probability distribution can be approximated by cubature points. CKF first approximates the mean and variance of probability distribution through a set of 2N (N is the dimension of the input random vector) cubature points with the same weight, propagates the above cubature points through nonlinear functions, and then calculates the mean and variance of the current approximate Gaussian distribution by the propagated cubature points.A set of cubature points are given byξ
i
W
i [14], where ξ
i is the ith cubature point and W
i is the corresponding weight
(7)
ξ
i
=
m
2
l
i
,
W
i
=
1
m
,
wherein i
=
1,2
,
…
,
m, m
=
2
N.Systematic state equation and observation equation are expressed in (4) and (6). Assuming the posterior density at time k
-
1 is known, the steps of CKF are shown as follows [14].Time update is as follows:(8a)
P
k
-
1
∣
k
-
1
=
S
k
-
1
∣
k
-
1
S
k
-
1
∣
k
-
1
T
,
(8b)
X
i
,
k
-
1
∣
k
-
1
=
S
k
-
1
∣
k
-
1
ξ
i
+
X
^
k
-
1
∣
k
-
1
,
(8c)
X
i
,
k
∣
k
-
1
*
=
f
X
i
,
k
-
1
∣
k
-
1
,
(8d)
X
^
k
∣
k
-
1
=
1
2
N
∑
i
=
1
2
N
X
i
,
k
∣
k
-
1
*
,
(8e)
P
k
∣
k
-
1
=
1
2
N
∑
i
=
1
2
N
X
i
,
k
∣
k
-
1
*
X
i
,
k
∣
k
-
1
*
T
-
X
^
k
∣
k
-
1
X
^
k
∣
k
-
1
T
+
Q
k
-
1
.Measurement update is as follows:(9a)
P
k
∣
k
-
1
=
S
k
∣
k
-
1
S
k
∣
k
-
1
T
,
(9b)
X
i
,
k
∣
k
-
1
=
S
k
∣
k
-
1
ξ
i
+
X
^
k
∣
k
-
1
,
(9c)
Z
i
,
k
∣
k
-
1
=
h
X
i
,
k
∣
k
-
1
,
(9d)
Z
^
k
∣
k
-
1
=
1
2
N
∑
i
=
1
2
N
Z
i
,
k
∣
k
-
1
,
(9e)
P
k
∣
k
-
1
z
z
=
1
2
N
∑
i
=
1
2
N
Z
i
,
k
∣
k
-
1
Z
i
,
k
∣
k
-
1
T
-
Z
^
k
∣
k
-
1
Z
^
k
∣
k
-
1
T
+
σ
r
2
,
(9f)
P
k
∣
k
-
1
x
z
=
1
2
N
∑
i
=
1
2
N
X
i
,
k
∣
k
-
1
Z
i
,
k
∣
k
-
1
T
-
X
^
k
∣
k
-
1
Z
^
k
∣
k
-
1
T
.With the new measurement vectorY
k, the estimation of the state vector X
^
k
∣
k and its covariance matrix P
k
∣
k at time k can be achieved by the following equations:(10a)
K
k
=
P
k
∣
k
-
1
x
z
P
k
∣
k
-
1
z
z
-
1
,
(10b)
X
^
k
∣
k
=
X
^
k
∣
k
-
1
+
K
k
Y
k
-
Z
^
k
∣
k
-
1
,
(10c)
P
k
∣
k
=
P
k
∣
k
-
1
-
K
k
P
k
∣
k
-
1
z
z
K
k
T
.CKF uses cubature rule and 2N cubature point sets ξ
i
W
i to compute the mean and variance of probability distribution without any linearization of a nonlinear model. Thus, the model can reach the third-order or higher.


### 2.4. Problem Statement-Poor Observability of Error in MEMS

A global observability analysis method proposed in [20] indicates that whether a state can be observed or not means if there exists unique solution of the equation. The problem that whether the initial misalignment error and MEMS gyro drift has unique solution will be preliminarily discussed in this section. Here in order to facilitate analysis, only the initial misalignment and gyro drift in MEMS are considered in the system temporarily. Equation (2) can be expressed as follows:
(11)
R
k
+
1
=
x
k
S
+
v
k
·
Δ
t
·
sin
α
^
k
-
δ
α
k
+
w
α
k
-
x
k
+
1
M
2
m
m
+
y
k
S
+
v
k
·
Δ
t
·
cos
⁡
α
^
k
-
δ
α
k
+
w
α
k
-
y
k
+
1
M
2
1
/
2
from which we can obtain the expression of δ
α
k:
(12)
δ
α
k
=
arcsin
R
k
+
1
2
-
x
k
S
-
x
k
+
1
M
2
-
y
k
S
-
y
k
+
1
M
2
-
v
k
2
Δ
t
2
2
v
k
Δ
t
x
k
S
-
x
k
+
1
M
2
+
y
k
S
-
y
k
+
1
M
2
-
β
-
α
^
k
,
where β
=
arctan
y
k
S
-
y
k
+
1
M
/
x
k
S
-
x
k
+
1
M. So δ
α
k can be calculated that means it is observable. However, the proportion of initial misalignment and MEMS gyro drift cannot be accurately known. Therefore, the observability is poor.As mentioned above, poor observability of MEMS is the main issue in cooperative navigation. The elaborate observability analysis and corresponding solutions will be provided in Section3 below.


## 2.1. Fundamental Theory of Cooperative Navigation Modeling

Master USV and slave USV are equipped with hydroacoustic communication modem, respectively, at their bottoms to achieve information dissemination and distance measurements between USVs. The masters accomplish navigation using their high-precision device (usually strap-down inertial navigation system combined with GPS) while the slaves use their own equipped Doppler velocity log (DVL) and attitude and heading reference system (AHRS) based on MEMS to proceed dead-reckoning.The principle of multi-USV cooperative navigation system is as follows. Firstly master USV sends out fixed frequency acoustic pulse signal according to preset time interval. With this signal, slave USV can calculate the relative distance according to underwater acoustic propagation delay. Then master USV sends out its position information by hydroacoustic modem in the form of broadcasting. After receiving this information, slave USV can conduct information fusion with dead-reckoning, position of master USV, and the relative distance and then correct its navigation result using filtering algorithms. Moreover, route planning of USV can be carried out to improve the observability and further improve the accuracy of the cooperative navigation system.Positioning errors will accumulate over time (as shown by the purple dotted ellipse shown in Figure1). Therefore, it cannot provide accurate positioning information for a long time. When the information sent by master USV is received, the slave USV can use cooperative navigation algorithm to correct the result of dead-reckoning and effectively reduce the positioning error (as shown by the green dotted ellipse shown in Figure 1). Through this cyclical filtering and correction, slave USV with low-precision navigation device can realize precise location for a long time [19, 20].Figure 1
Principle of cooperative navigation system.


## 2.2. System Model of Single Master USV

When the follower does not receive the information sent by leaders, it can only use its own DVL and low accuracy MEMS to conduct dead-reckoning, and the motion model can be expressed as follows:(1)
x
k
+
1
S
=
x
k
S
+
v
k
·
Δ
t
·
sin
α
k
,
y
k
+
1
S
=
y
k
S
+
v
k
·
Δ
t
·
cos
⁡
α
k
,
where x
k
S, y
k
S represent the position of follower USV at time k, v
k represents the speed, Δ
t means the time interval, and α
k means the heading angle at time k. However, the input model of system in practice is as follows:
(2)
v
k
=
1
+
δ
C
k
v
^
k
+
w
v
k
,
α
k
=
α
^
k
-
δ
α
k
+
w
α
k
,
wherein v
^
k is the velocity measurements and δ
C
k means the scale factor error. α
^
k is the heading angle measured by MEMS, δ
α
k
=
Δ
α
k
+
k
ε
Δ
t is error of α
^
k, Δ
α
k is the initial misalignment angle, and ε is gyro drift. Putting (2) into (1), the actual equation of states is as follows:
(3)
x
k
+
1
S
=
x
k
S
+
1
+
δ
C
k
v
^
k
+
w
v
k
·
Δ
t
·
sin
α
^
k
-
Δ
α
k
-
ε
k
·
k
·
Δ
t
+
w
α
k
,
y
k
+
1
S
=
y
k
S
+
1
+
δ
C
k
v
^
k
+
w
v
k
·
Δ
t
·
cos
⁡
α
^
k
-
Δ
α
k
-
ε
k
·
k
·
Δ
t
+
w
α
k
,
δ
C
k
+
1
=
δ
C
k
,
Δ
α
k
+
1
=
Δ
α
k
,
ε
k
+
1
=
ε
k
.Motion model can be simplified as(4)
X
k
+
1
=
X
k
+
T
u
k
,
w
k
in which X
k
=
x
k
,
y
k
,
Δ
α
k
,
ε
k
,
δ
C
k
T means the state vector of USV at time k, u
k
=
[
v
^
k
α
^
k
]
T, T
(
u
k
,
w
k
) stands for the nonlinear terms, the process noise is w
k
=
[
w
v
k
w
α
k
]
T, w
v
k
~
N
(
0
,
σ
v
k
2
), and w
α
k
~
N
(
0
,
σ
α
k
2
). The noise covariance matrix is
(5)
Q
k
=
E
w
k
w
k
T
=
σ
v
k
2
σ
α
k
2
.Systematic observation equation can be expressed as:(6)
Y
k
=
R
k
+
V
k
=
x
k
M
-
x
k
S
2
+
y
k
M
-
y
k
S
2
+
V
k
,
wherein Y
k stands for the measurement of distance R
k between follower and leader, (
x
k
M
,
y
k
M
) represents the position of leader USV at time k while x
k
S
,
y
k
S represents the position of follower, and V
k
~
N
(
0
,
σ
r
2
) is the measurement noise.


## 2.3. Cubature Kalman Filter (CKF)

Cubature Kalman filter (CKF) is proposed based on the spherical-radial cubature criterion [14]. The core of the method is that the mean and variance of probability distribution can be approximated by cubature points. CKF first approximates the mean and variance of probability distribution through a set of 2N (N is the dimension of the input random vector) cubature points with the same weight, propagates the above cubature points through nonlinear functions, and then calculates the mean and variance of the current approximate Gaussian distribution by the propagated cubature points.A set of cubature points are given byξ
i
W
i [14], where ξ
i is the ith cubature point and W
i is the corresponding weight
(7)
ξ
i
=
m
2
l
i
,
W
i
=
1
m
,
wherein i
=
1,2
,
…
,
m, m
=
2
N.Systematic state equation and observation equation are expressed in (4) and (6). Assuming the posterior density at time k
-
1 is known, the steps of CKF are shown as follows [14].Time update is as follows:(8a)
P
k
-
1
∣
k
-
1
=
S
k
-
1
∣
k
-
1
S
k
-
1
∣
k
-
1
T
,
(8b)
X
i
,
k
-
1
∣
k
-
1
=
S
k
-
1
∣
k
-
1
ξ
i
+
X
^
k
-
1
∣
k
-
1
,
(8c)
X
i
,
k
∣
k
-
1
*
=
f
X
i
,
k
-
1
∣
k
-
1
,
(8d)
X
^
k
∣
k
-
1
=
1
2
N
∑
i
=
1
2
N
X
i
,
k
∣
k
-
1
*
,
(8e)
P
k
∣
k
-
1
=
1
2
N
∑
i
=
1
2
N
X
i
,
k
∣
k
-
1
*
X
i
,
k
∣
k
-
1
*
T
-
X
^
k
∣
k
-
1
X
^
k
∣
k
-
1
T
+
Q
k
-
1
.Measurement update is as follows:(9a)
P
k
∣
k
-
1
=
S
k
∣
k
-
1
S
k
∣
k
-
1
T
,
(9b)
X
i
,
k
∣
k
-
1
=
S
k
∣
k
-
1
ξ
i
+
X
^
k
∣
k
-
1
,
(9c)
Z
i
,
k
∣
k
-
1
=
h
X
i
,
k
∣
k
-
1
,
(9d)
Z
^
k
∣
k
-
1
=
1
2
N
∑
i
=
1
2
N
Z
i
,
k
∣
k
-
1
,
(9e)
P
k
∣
k
-
1
z
z
=
1
2
N
∑
i
=
1
2
N
Z
i
,
k
∣
k
-
1
Z
i
,
k
∣
k
-
1
T
-
Z
^
k
∣
k
-
1
Z
^
k
∣
k
-
1
T
+
σ
r
2
,
(9f)
P
k
∣
k
-
1
x
z
=
1
2
N
∑
i
=
1
2
N
X
i
,
k
∣
k
-
1
Z
i
,
k
∣
k
-
1
T
-
X
^
k
∣
k
-
1
Z
^
k
∣
k
-
1
T
.With the new measurement vectorY
k, the estimation of the state vector X
^
k
∣
k and its covariance matrix P
k
∣
k at time k can be achieved by the following equations:(10a)
K
k
=
P
k
∣
k
-
1
x
z
P
k
∣
k
-
1
z
z
-
1
,
(10b)
X
^
k
∣
k
=
X
^
k
∣
k
-
1
+
K
k
Y
k
-
Z
^
k
∣
k
-
1
,
(10c)
P
k
∣
k
=
P
k
∣
k
-
1
-
K
k
P
k
∣
k
-
1
z
z
K
k
T
.CKF uses cubature rule and 2N cubature point sets ξ
i
W
i to compute the mean and variance of probability distribution without any linearization of a nonlinear model. Thus, the model can reach the third-order or higher.


## 2.4. Problem Statement-Poor Observability of Error in MEMS

A global observability analysis method proposed in [20] indicates that whether a state can be observed or not means if there exists unique solution of the equation. The problem that whether the initial misalignment error and MEMS gyro drift has unique solution will be preliminarily discussed in this section. Here in order to facilitate analysis, only the initial misalignment and gyro drift in MEMS are considered in the system temporarily. Equation (2) can be expressed as follows:
(11)
R
k
+
1
=
x
k
S
+
v
k
·
Δ
t
·
sin
α
^
k
-
δ
α
k
+
w
α
k
-
x
k
+
1
M
2
m
m
+
y
k
S
+
v
k
·
Δ
t
·
cos
⁡
α
^
k
-
δ
α
k
+
w
α
k
-
y
k
+
1
M
2
1
/
2
from which we can obtain the expression of δ
α
k:
(12)
δ
α
k
=
arcsin
R
k
+
1
2
-
x
k
S
-
x
k
+
1
M
2
-
y
k
S
-
y
k
+
1
M
2
-
v
k
2
Δ
t
2
2
v
k
Δ
t
x
k
S
-
x
k
+
1
M
2
+
y
k
S
-
y
k
+
1
M
2
-
β
-
α
^
k
,
where β
=
arctan
y
k
S
-
y
k
+
1
M
/
x
k
S
-
x
k
+
1
M. So δ
α
k can be calculated that means it is observable. However, the proportion of initial misalignment and MEMS gyro drift cannot be accurately known. Therefore, the observability is poor.As mentioned above, poor observability of MEMS is the main issue in cooperative navigation. The elaborate observability analysis and corresponding solutions will be provided in Section3 below.


## 3. Cooperative Navigation Modeling and Observability Analysis

### 3.1. Observability Analysis

Observability directly determines the colocation accuracy of cooperative navigation system [20]; the observability theory of linear systems will be no longer suitable to multi-AUV cooperative navigation with nonlinear system. The idea here is to linearize the nonlinear model firstly and then to use observability theory of linearized systems for observability analysis. The system in this paper is stochastic system. However, the observability analysis method of stochastic system can be equivalent to the corresponding deterministic system when Q
>
0 and R
>
0 [21, 22], which is satisfied in our paper.Equation (1) can be abbreviated as X
k
+
1
=
f
(
X
k
,
u
k
)
=
X
k
+
T
(
u
k
,
w
k
), Y
k
=
h
[
·
], wherein X
k means the state of USV at time k and T
(
u
k
,
w
k
) stands for the nonlinear terms. Assume X
k
=
[
x
k
S
y
k
S
Δ
α
k
ε
k
δ
C
k
]
T and δ
α
k
=
Δ
α
k
+
ε
·
k
·
Δ
t; then state equations can be written as
(13)
x
k
+
1
S

=

x
k
S
+
1
+
δ
C
k
v
^
k
+
w
v
k
·
Δ
t
·
sin
α
^
k
-
δ
α
k
+
w
α
k
,
y
k
+
1
S

=

y
k
S
+
1
+
δ
C
k
v
^
k
+
w
v
k
·
Δ
t
·
cos
⁡
α
^
k
-
δ
α
k
+
w
α
k
,
δ
C
k
+
1
=
δ
C
k
,
Δ
α
k
+
1
=
Δ
α
k
,
ε
k
+
1
=
ε
k
.Therefore, the Jacobi matrix is(14)
F
k
=
∂
f
∂
X
k
=
1
0
-
v
k
Δ
t
cos
⁡
α
^
k
-
δ
α
k
+
w
α
k
-
v
k
k
Δ
t
2
cos
⁡
α
^
k
-
δ
α
k
+
w
α
k
(
v
^
k
+
w
v
k
)
Δ
t
sin
α
^
k
-
δ
α
k
+
w
α
k
0
1
v
k
Δ
t
sin
α
^
k
-
δ
α
k
+
w
α
k
v
k
k
Δ
t
2
sin
α
^
k
-
δ
α
k
+
w
α
k
(
v
^
k
+
w
v
k
)
Δ
t
cos
⁡
α
^
k
-
δ
α
k
+
w
α
k
0
0
1
0
0
0
0
0
1
0
0
0
0
0
1
,
H
k
=
∂
h
[
·
]
∂
X
k
=
x
k
S
-
x
k
M
R
k
y
k
S
-
y
k
M
R
k
0
0
0
=
Δ
x
k
R
k
Δ
y
k
R
k
0
0
0
=
sin
α
^
k
cos
⁡
α
^
k
0
0
0
,where Δ
x
k
=
x
k
S
-
x
k
M, Δ
y
k
=
y
k
S
-
y
k
M.According to rank criterion, observability matrix is(15)
Γ
=
H
k
H
k
-
1
F
k
-
1
H
k
-
2
F
k
-
2
H
k
-
3
F
k
-
3
H
k
-
4
F
k
-
4
=
sin
α
^
k
cos
⁡
α
^
k
0
0
0
sin
α
^
k
-
1
cos
⁡
α
^
k
-
1
-
v
k
-
1
Δ
t
sin
α
k
-
1
′
-
v
k
-
1
k
-
1
Δ
t
2
sin
α
k
-
1
′
v
^
k
-
1
+
w
v
k
-
1
Δ
t
cos
⁡
α
k
-
1
′
sin
α
^
k
-
2
cos
⁡
α
^
k
-
2
-
v
k
-
2
Δ
t
sin
α
k
-
2
′
-
v
k
-
2
k
-
2
Δ
t
2
sin
α
k
-
2
′
v
^
k
-
2
+
w
v
k
-
2
Δ
t
cos
⁡
α
k
-
2
′
sin
α
^
k
-
3
cos
⁡
α
^
k
-
3
-
v
k
-
3
Δ
t
sin
α
k
-
3
′
-
v
k
-
3
k
-
3
Δ
t
2
sin
α
k
-
3
′
v
^
k
-
3
+
w
v
k
-
3
Δ
t
cos
⁡
α
k
-
3
′
sin
α
^
k
-
4
cos
⁡
α
^
k
-
4
-
v
k
-
4
Δ
t
sin
α
k
-
4
′
-
v
k
-
4
k
-
4
Δ
t
2
sin
α
k
-
4
′
v
^
k
-
4
+
w
v
k
-
4
Δ
t
cos
⁡
α
k
-
4
′
,where α
j
′
=
(
δ
α
j
-
ω
α
j
)
(
j
=
k
-
1
,
k
-
2
,
k
-
3
⋯
). The linear relativity of 3rd column and 4th column in Γ would be discussed here. It can be seen that Γ
2,3
·
Δ
t
/
Γ
2,4
=
1
/
k
-
1, Γ
3,3
·
Δ
t
/
Γ
3,4
=
1
/
k
-
2
Γ
4,3
·
Δ
t
/
Γ
4,4
=
1
/
k
-
3, Γ
5,3
·
Δ
t
/
Γ
5,4
=
1
/
k
-
4, and Γ
i
,
j stands for the element in row i, column j of the system observation matrix. Figure 2 shows the variation of corresponding value over time.Figure 2
Relative analysis of the column in observability matrix.As can be seen from Figure2, when time k
>
100, 1
/
k
-
1, 1
/
k
-
2, 1
/
k
-
3, 1
/
k
-
4 is approximately equal. That means the third and fourth columns of the observation matrix Γ become linear-relative and the linear correlation becomes stronger over time. It can be concluded that the observability is poor and it gets worse as time flies.


### 3.2. Resolve Route 1: Dual-Model Filtering Algorithm

Both of initial misalignment and gyro drift can cause the increase of error angle. The error coupled together results in poor observability. Error caused by the initial misalignment angle is a constant value while the error caused by gyro drift gradually increases over time. Their properties are shown in Figure3.Characteristics of initial misalignment angle and gyro drift.
(a)
(b)The proportion of initial misalignment angle and gyro drift in total error angle over time are shown in Table1 and Figure 4 (assuming initial misalignment is 90° and gyro drift is 10°/h).Table 1
The proportion varies in the first 500 s.
Time (s)
Proportion of initial misalignment
Proportion of gyro drift
100
0.9970
0.0030
200
0.9939
0.0061
300
0.9909
0.0091
400
0.9878
0.0122
500
0.9848
0.0152Figure 4
The proportion of initial misalignment angle and gyro drift.From the above analysis, it can be seen that with the passing of time, gyro drift will take increasingly large proportion while the initial misalignment of MEMS occupies absolutely large part of total deviation angle in the beginning. Therefore, we can cut off the entire time into two phases. In the first phase initial misalignment angle is estimated and modified alone while the gyro drift is estimated in the second stage. As long as length of the first phase is short enough, it will not affect the estimation accuracy; that is,(16)
Subsection 1:
X
k
1
=
x
k
S
y
k
S
Δ
α
k
δ
C
k
T
Subsection 2:
X
k
2
=
x
k
S
y
k
S
ε
k
δ
C
k
T
.Here we will take a look at observability of the system after subsection. When the system state vector is selected asX
k
1
=
[
x
k
S
y
k
S
Δ
α
k
δ
C
k
]
T, the Jacobi matrix is(17)
F
k
=
∂
f
∂
X
k
1
=
1
0
-
v
k
Δ
t
cos
⁡
α
^
k
-
δ
α
k
+
w
α
k
v
^
k
+
w
v
k
Δ
t
sin
α
^
k
-
δ
α
k
+
w
α
k
0
1
v
k
Δ
t
sin
α
^
k
-
δ
α
k
+
w
α
k
v
^
k
+
w
v
k
Δ
t
cos
⁡
α
^
k
-
δ
α
k
+
w
α
k
0
0
1
0
0
0
0
1
,
H
k
=
∂
h
[
·
]
∂
X
k
1
=
x
k
S
-
x
k
M
R
k
y
k
S
-
y
k
M
R
k
0
0
=
sin
α
^
k
cos
⁡
α
^
k
0
0
.According to the rank criterion, observability matrix can be written as follows:(18)
Γ
=
H
k
H
k
-
1
F
k
-
1
H
k
-
2
F
k
-
2
H
k
-
3
F
k
-
3
=
sin
α
^
k
cos
⁡
α
^
k
0
0
sin
α
^
k
-
1
cos
⁡
α
^
k
-
1
-
v
k
-
1
Δ
t
sin
δ
α
k
-
1
-
w
α
k
-
1
v
^
k
-
1
+
w
v
k
-
1
Δ
t
cos
⁡
δ
α
k
-
1
-
w
α
k
-
1
sin
α
^
k
-
2
cos
⁡
α
^
k
-
2
-
v
k
-
2
Δ
t
sin
δ
α
k
-
2
-
w
α
k
-
2
v
^
k
-
2
+
w
v
k
-
2
Δ
t
cos
⁡
δ
α
k
-
2
-
w
α
k
-
2
sin
α
^
k
-
3
cos
⁡
α
^
k
-
3
-
v
k
-
3
Δ
t
sin
δ
α
k
-
3
-
w
α
k
-
3
v
^
k
-
3
+
w
v
k
-
3
Δ
t
cos
⁡
δ
α
k
-
3
-
w
α
k
-
3
,whererank
⁡
Γ
=
4; that is to say, the matrix Γ is full rank, so the system is observable. For the second phase of cooperative navigation system, we can get the same conclusion. Therefore, it can be concluded that the system becomes observable whenever in the first phase or the second phase after subsection.


### 3.3. Resolve Route 2: Reverse Algorithm

In Section3.2, it shows that after using the segmented dual-model filtering, state observability can be significantly improved. However, the proportion of gyro drift in total deviation angle gradually increases as time flies. If the time of first stage is too long, the gyro drift will seriously affect the accuracy of model and the precision of initial misalignment estimation. Currently MEMS devices are with low precision that large initial misalignment and large gyro drift coexist. Therefore, how to estimate the large initial misalignment angle in short time is another problem we have to face.Sampling data in the cooperative navigation system is typically a time sequence. The filter works with chronological sequence of real-time processing. The real-time results can be obtained without data storage. However, similar to strap-down inertial navigation system [18, 23, 24], one significant feature of cooperative navigation system is that all the information used in filtering process can be completely stored. That means the data can be copied into many identical copies, and each copy can be processed by different methods without interference between them. This feature is also called “the diverse existence of mathematical platform.”It is easily conceivable that if storage capacity of the navigation computer is large enough and computing power is strong enough, the data can be stored and analyzed many times in very short time. By repeating the analysis in forward and backward, it is possible to improve the precision or shorten the time cost.Based on this idea, we can take a short period of time as the first stage and carry forward filtering at the first time while the data is saved. Then through repeated use of the data for backward and forward filtering, the large initial misalignment angle can be accurately estimated in short time. The dream of improving accuracy of the initial misalignment angle in short time can be achieved. The process is shown in Figure5.Figure 5
Schematic diagram of reverse algorithm.The corresponding forward and reverse state equations of the algorithm are as follows. In order to distinguish between them, we use “→” and “←” to present forward and reverse, respectively.(1) Forward filtering is as follows:(19)
x
→
k
+
1
S

=

x
→
k
S
+
1
+
δ
C
k
v
^
→
k
+
w
v
k
·
Δ
t
·
sin
α
^
k
-
Δ
α
k
+
w
α
k
,
y
→
k
+
1
S

=

y
→
k
S
+
1
+
δ
C
k
v
^
→
k
+
w
v
k
·
Δ
t
·
cos
⁡
α
^
k
-
Δ
α
k
+
w
α
k
,
δ
C
k
+
1
=
δ
C
k
,
Δ
α
k
+
1
=
Δ
α
k
.(2) Backward filtering is as follows:(20)
x
←
k
+
1
S

=

x
←
k
S
+
1
+
δ
C
k
v
^
←
k
+
w
v
k
·
Δ
t
·
sin
α
^
k
-
Δ
α
k
+
w
α
k
,
y
←
k
+
1
S

=

y
←
k
S
+
1
+
δ
C
k
v
^
←
k
+
w
v
k
·
Δ
t
·
cos
⁡
α
^
k
-
Δ
α
k
+
w
α
k
,
δ
C
k
+
1
=
δ
C
k
,
Δ
α
k
+
1
=
Δ
α
k
.The whole “subsection + reverse” filtering process is summarized as follows (as shown in Figure6): (1) first, initialize the filter parameters; (2) in the first subsection (a short period of time), conduct the initial misalignment estimation ignoring effects of gyro drift and store the data at the same time; (3) use the stored data to estimate the initial misalignment reversely; (4) return to steps (2) and (3) and repeat two or three times until the covariance of the state decreases in demand range (the selection of repeated time is discussed in Section 6); (5) after compensation of initial misalignment using estimation result in the first subsection, estimate the gyro drift and conduct cooperative navigation filtering in the second subsection.Figure 6
“Subsection + reverse” filtering timing diagram.


## 3.1. Observability Analysis

Observability directly determines the colocation accuracy of cooperative navigation system [20]; the observability theory of linear systems will be no longer suitable to multi-AUV cooperative navigation with nonlinear system. The idea here is to linearize the nonlinear model firstly and then to use observability theory of linearized systems for observability analysis. The system in this paper is stochastic system. However, the observability analysis method of stochastic system can be equivalent to the corresponding deterministic system when Q
>
0 and R
>
0 [21, 22], which is satisfied in our paper.Equation (1) can be abbreviated as X
k
+
1
=
f
(
X
k
,
u
k
)
=
X
k
+
T
(
u
k
,
w
k
), Y
k
=
h
[
·
], wherein X
k means the state of USV at time k and T
(
u
k
,
w
k
) stands for the nonlinear terms. Assume X
k
=
[
x
k
S
y
k
S
Δ
α
k
ε
k
δ
C
k
]
T and δ
α
k
=
Δ
α
k
+
ε
·
k
·
Δ
t; then state equations can be written as
(13)
x
k
+
1
S

=

x
k
S
+
1
+
δ
C
k
v
^
k
+
w
v
k
·
Δ
t
·
sin
α
^
k
-
δ
α
k
+
w
α
k
,
y
k
+
1
S

=

y
k
S
+
1
+
δ
C
k
v
^
k
+
w
v
k
·
Δ
t
·
cos
⁡
α
^
k
-
δ
α
k
+
w
α
k
,
δ
C
k
+
1
=
δ
C
k
,
Δ
α
k
+
1
=
Δ
α
k
,
ε
k
+
1
=
ε
k
.Therefore, the Jacobi matrix is(14)
F
k
=
∂
f
∂
X
k
=
1
0
-
v
k
Δ
t
cos
⁡
α
^
k
-
δ
α
k
+
w
α
k
-
v
k
k
Δ
t
2
cos
⁡
α
^
k
-
δ
α
k
+
w
α
k
(
v
^
k
+
w
v
k
)
Δ
t
sin
α
^
k
-
δ
α
k
+
w
α
k
0
1
v
k
Δ
t
sin
α
^
k
-
δ
α
k
+
w
α
k
v
k
k
Δ
t
2
sin
α
^
k
-
δ
α
k
+
w
α
k
(
v
^
k
+
w
v
k
)
Δ
t
cos
⁡
α
^
k
-
δ
α
k
+
w
α
k
0
0
1
0
0
0
0
0
1
0
0
0
0
0
1
,
H
k
=
∂
h
[
·
]
∂
X
k
=
x
k
S
-
x
k
M
R
k
y
k
S
-
y
k
M
R
k
0
0
0
=
Δ
x
k
R
k
Δ
y
k
R
k
0
0
0
=
sin
α
^
k
cos
⁡
α
^
k
0
0
0
,where Δ
x
k
=
x
k
S
-
x
k
M, Δ
y
k
=
y
k
S
-
y
k
M.According to rank criterion, observability matrix is(15)
Γ
=
H
k
H
k
-
1
F
k
-
1
H
k
-
2
F
k
-
2
H
k
-
3
F
k
-
3
H
k
-
4
F
k
-
4
=
sin
α
^
k
cos
⁡
α
^
k
0
0
0
sin
α
^
k
-
1
cos
⁡
α
^
k
-
1
-
v
k
-
1
Δ
t
sin
α
k
-
1
′
-
v
k
-
1
k
-
1
Δ
t
2
sin
α
k
-
1
′
v
^
k
-
1
+
w
v
k
-
1
Δ
t
cos
⁡
α
k
-
1
′
sin
α
^
k
-
2
cos
⁡
α
^
k
-
2
-
v
k
-
2
Δ
t
sin
α
k
-
2
′
-
v
k
-
2
k
-
2
Δ
t
2
sin
α
k
-
2
′
v
^
k
-
2
+
w
v
k
-
2
Δ
t
cos
⁡
α
k
-
2
′
sin
α
^
k
-
3
cos
⁡
α
^
k
-
3
-
v
k
-
3
Δ
t
sin
α
k
-
3
′
-
v
k
-
3
k
-
3
Δ
t
2
sin
α
k
-
3
′
v
^
k
-
3
+
w
v
k
-
3
Δ
t
cos
⁡
α
k
-
3
′
sin
α
^
k
-
4
cos
⁡
α
^
k
-
4
-
v
k
-
4
Δ
t
sin
α
k
-
4
′
-
v
k
-
4
k
-
4
Δ
t
2
sin
α
k
-
4
′
v
^
k
-
4
+
w
v
k
-
4
Δ
t
cos
⁡
α
k
-
4
′
,where α
j
′
=
(
δ
α
j
-
ω
α
j
)
(
j
=
k
-
1
,
k
-
2
,
k
-
3
⋯
). The linear relativity of 3rd column and 4th column in Γ would be discussed here. It can be seen that Γ
2,3
·
Δ
t
/
Γ
2,4
=
1
/
k
-
1, Γ
3,3
·
Δ
t
/
Γ
3,4
=
1
/
k
-
2
Γ
4,3
·
Δ
t
/
Γ
4,4
=
1
/
k
-
3, Γ
5,3
·
Δ
t
/
Γ
5,4
=
1
/
k
-
4, and Γ
i
,
j stands for the element in row i, column j of the system observation matrix. Figure 2 shows the variation of corresponding value over time.Figure 2
Relative analysis of the column in observability matrix.As can be seen from Figure2, when time k
>
100, 1
/
k
-
1, 1
/
k
-
2, 1
/
k
-
3, 1
/
k
-
4 is approximately equal. That means the third and fourth columns of the observation matrix Γ become linear-relative and the linear correlation becomes stronger over time. It can be concluded that the observability is poor and it gets worse as time flies.


## 3.2. Resolve Route 1: Dual-Model Filtering Algorithm

Both of initial misalignment and gyro drift can cause the increase of error angle. The error coupled together results in poor observability. Error caused by the initial misalignment angle is a constant value while the error caused by gyro drift gradually increases over time. Their properties are shown in Figure3.Characteristics of initial misalignment angle and gyro drift.
(a)
(b)The proportion of initial misalignment angle and gyro drift in total error angle over time are shown in Table1 and Figure 4 (assuming initial misalignment is 90° and gyro drift is 10°/h).Table 1
The proportion varies in the first 500 s.
Time (s)
Proportion of initial misalignment
Proportion of gyro drift
100
0.9970
0.0030
200
0.9939
0.0061
300
0.9909
0.0091
400
0.9878
0.0122
500
0.9848
0.0152Figure 4
The proportion of initial misalignment angle and gyro drift.From the above analysis, it can be seen that with the passing of time, gyro drift will take increasingly large proportion while the initial misalignment of MEMS occupies absolutely large part of total deviation angle in the beginning. Therefore, we can cut off the entire time into two phases. In the first phase initial misalignment angle is estimated and modified alone while the gyro drift is estimated in the second stage. As long as length of the first phase is short enough, it will not affect the estimation accuracy; that is,(16)
Subsection 1:
X
k
1
=
x
k
S
y
k
S
Δ
α
k
δ
C
k
T
Subsection 2:
X
k
2
=
x
k
S
y
k
S
ε
k
δ
C
k
T
.Here we will take a look at observability of the system after subsection. When the system state vector is selected asX
k
1
=
[
x
k
S
y
k
S
Δ
α
k
δ
C
k
]
T, the Jacobi matrix is(17)
F
k
=
∂
f
∂
X
k
1
=
1
0
-
v
k
Δ
t
cos
⁡
α
^
k
-
δ
α
k
+
w
α
k
v
^
k
+
w
v
k
Δ
t
sin
α
^
k
-
δ
α
k
+
w
α
k
0
1
v
k
Δ
t
sin
α
^
k
-
δ
α
k
+
w
α
k
v
^
k
+
w
v
k
Δ
t
cos
⁡
α
^
k
-
δ
α
k
+
w
α
k
0
0
1
0
0
0
0
1
,
H
k
=
∂
h
[
·
]
∂
X
k
1
=
x
k
S
-
x
k
M
R
k
y
k
S
-
y
k
M
R
k
0
0
=
sin
α
^
k
cos
⁡
α
^
k
0
0
.According to the rank criterion, observability matrix can be written as follows:(18)
Γ
=
H
k
H
k
-
1
F
k
-
1
H
k
-
2
F
k
-
2
H
k
-
3
F
k
-
3
=
sin
α
^
k
cos
⁡
α
^
k
0
0
sin
α
^
k
-
1
cos
⁡
α
^
k
-
1
-
v
k
-
1
Δ
t
sin
δ
α
k
-
1
-
w
α
k
-
1
v
^
k
-
1
+
w
v
k
-
1
Δ
t
cos
⁡
δ
α
k
-
1
-
w
α
k
-
1
sin
α
^
k
-
2
cos
⁡
α
^
k
-
2
-
v
k
-
2
Δ
t
sin
δ
α
k
-
2
-
w
α
k
-
2
v
^
k
-
2
+
w
v
k
-
2
Δ
t
cos
⁡
δ
α
k
-
2
-
w
α
k
-
2
sin
α
^
k
-
3
cos
⁡
α
^
k
-
3
-
v
k
-
3
Δ
t
sin
δ
α
k
-
3
-
w
α
k
-
3
v
^
k
-
3
+
w
v
k
-
3
Δ
t
cos
⁡
δ
α
k
-
3
-
w
α
k
-
3
,whererank
⁡
Γ
=
4; that is to say, the matrix Γ is full rank, so the system is observable. For the second phase of cooperative navigation system, we can get the same conclusion. Therefore, it can be concluded that the system becomes observable whenever in the first phase or the second phase after subsection.


## 3.3. Resolve Route 2: Reverse Algorithm

In Section3.2, it shows that after using the segmented dual-model filtering, state observability can be significantly improved. However, the proportion of gyro drift in total deviation angle gradually increases as time flies. If the time of first stage is too long, the gyro drift will seriously affect the accuracy of model and the precision of initial misalignment estimation. Currently MEMS devices are with low precision that large initial misalignment and large gyro drift coexist. Therefore, how to estimate the large initial misalignment angle in short time is another problem we have to face.Sampling data in the cooperative navigation system is typically a time sequence. The filter works with chronological sequence of real-time processing. The real-time results can be obtained without data storage. However, similar to strap-down inertial navigation system [18, 23, 24], one significant feature of cooperative navigation system is that all the information used in filtering process can be completely stored. That means the data can be copied into many identical copies, and each copy can be processed by different methods without interference between them. This feature is also called “the diverse existence of mathematical platform.”It is easily conceivable that if storage capacity of the navigation computer is large enough and computing power is strong enough, the data can be stored and analyzed many times in very short time. By repeating the analysis in forward and backward, it is possible to improve the precision or shorten the time cost.Based on this idea, we can take a short period of time as the first stage and carry forward filtering at the first time while the data is saved. Then through repeated use of the data for backward and forward filtering, the large initial misalignment angle can be accurately estimated in short time. The dream of improving accuracy of the initial misalignment angle in short time can be achieved. The process is shown in Figure5.Figure 5
Schematic diagram of reverse algorithm.The corresponding forward and reverse state equations of the algorithm are as follows. In order to distinguish between them, we use “→” and “←” to present forward and reverse, respectively.(1) Forward filtering is as follows:(19)
x
→
k
+
1
S

=

x
→
k
S
+
1
+
δ
C
k
v
^
→
k
+
w
v
k
·
Δ
t
·
sin
α
^
k
-
Δ
α
k
+
w
α
k
,
y
→
k
+
1
S

=

y
→
k
S
+
1
+
δ
C
k
v
^
→
k
+
w
v
k
·
Δ
t
·
cos
⁡
α
^
k
-
Δ
α
k
+
w
α
k
,
δ
C
k
+
1
=
δ
C
k
,
Δ
α
k
+
1
=
Δ
α
k
.(2) Backward filtering is as follows:(20)
x
←
k
+
1
S

=

x
←
k
S
+
1
+
δ
C
k
v
^
←
k
+
w
v
k
·
Δ
t
·
sin
α
^
k
-
Δ
α
k
+
w
α
k
,
y
←
k
+
1
S

=

y
←
k
S
+
1
+
δ
C
k
v
^
←
k
+
w
v
k
·
Δ
t
·
cos
⁡
α
^
k
-
Δ
α
k
+
w
α
k
,
δ
C
k
+
1
=
δ
C
k
,
Δ
α
k
+
1
=
Δ
α
k
.The whole “subsection + reverse” filtering process is summarized as follows (as shown in Figure6): (1) first, initialize the filter parameters; (2) in the first subsection (a short period of time), conduct the initial misalignment estimation ignoring effects of gyro drift and store the data at the same time; (3) use the stored data to estimate the initial misalignment reversely; (4) return to steps (2) and (3) and repeat two or three times until the covariance of the state decreases in demand range (the selection of repeated time is discussed in Section 6); (5) after compensation of initial misalignment using estimation result in the first subsection, estimate the gyro drift and conduct cooperative navigation filtering in the second subsection.Figure 6
“Subsection + reverse” filtering timing diagram.


## 4. Simulation Analysis

To verify effectiveness of the proposed algorithms and validity of observability analysis above, simulation is conducted in this section. Simulation environment settings are as follows.(1)
A and C are leaders and B is follower. Three USVs are all equipped with underwater acoustic equipment for communication and ranging. A and C send their position alternately; time intervalT is set to 5 s.
(2)
The initial positions of A, B, and C are (0,200), (0,0), and (200,0). A and C are equipped with GPS receiver that can accomplish self-positioning, and the positioning error does not accumulate. The mean of positioning error is 0 and variance is set to 2 m.
(3)
USV B can only conduct dead-reckoning navigation when the information of A or C is not received. The speed and heading angle are given by DVL and MEMS, respectively; the initial misalignment is set to 90°; gyro drift is 10°/h.
(4)
δ
C
k of speed measured by DVL in USV B is 0.005, and the noise is set to zero mean and variance is 0.5 m/s.
(5)
Noise of ranging measurement is set with zero mean and variance is set to 2 m.
(6)
A, B, and C are in linear motion with speedup and speed-down. The speed changes between 8 m/s and 15 m/s by ±0.05 m/s2 without considering the dynamics and current disturbance.The filter is initialized as follows:P
0
=
diag
⁡
10
m
2
10
m
2
0.01
2
8
0
°
2
2
0
°
/
h
2. The trajectories of three USVs are shown in Figure 7.Figure 7
Trajectories of USVs.Results comparison of corresponding methods is given in Section4.1. It includes comparison between methods with segmentation and without segmentation and also comparison of method without reverse filtering and with reverse filtering (reverse once and reverse twice).


### 4.1. Comparison of Filtering Performance with and without Subsection

The whole simulation time is set as 6000 seconds and we take 1500 seconds for the first stage. According to “dual-model algorithm with subsection” proposed in Section4, the curves of estimation are shown in Figures 8 and 9.Figure 8
Estimation comparison with and without subsection.Comparison of filtering performance with and without subsection.
(a)
Initial alignment estimation (enlarge view)
(b)
Gyro drift estimation in long time without subsection
(c)
Colocation error comparisonFrom the comparison in Figures8 and 9, the following conclusions can be drawn.(1)
As can be seen from Figure8, the initial misalignment angle without segmentation begins to converge within a certain range but diverges later (as previously analyzed, observability grows gradually worse over time). The error estimation accuracy is better than that without segment, and the modification at the end of first stage effectively suppresses initial misalignment angle divergence.
(2)
It can be seen in Figures8 and 9(b) that the estimation performance of gyro drift without segmentation is poor. It gradually converges in a very long time while it converges quickly after subsection.
(3)
As can be seen in Figure9(a), the estimation of the initial misalignment angle in subsection 1 is better than that without subsection, but there is still a tendency of divergence.
(4)
Figure9(c) tells us that due to better estimation of the initial misalignment angle and gyro drift after segmented, the performance of colocation is better than nonsegments. However, such divergent trend still exists in colocation estimation error due to certain bias of initial heading angle estimation.In summary, the performance with segment is better than that without segmentation. But there is still space for further improvement especially in how to improve estimation accuracy of the initial misalignment angle of the first phase.


### 4.2. “Subsection + Reverse” Algorithm

The reverse filtering is introduced on the base of subsection in this section. It gives comparison of methods without reverse, reverse once, and reverse twice. The velocity curve of USV varies between 8 m/s and 15 m/s without reverse. It moves in a straight line with acceleration and deceleration. Figure10 shows velocity with reverse once and with reverse twice. The corresponding performance comparison is shown in Figure 11.Velocity of follower USV.
(a)
Velocity with reverse once
(b)
Velocity with reverse twiceComparison with no reverse, reverse once, and reverse twice.
(a)
Estimation comparison in whole process
(b)
The partial enlarged view of initial misalignment
(c)
An enlarged view of the gyro drift estimation
(d)
Cooperative navigation error estimationTo make the description more clearly, several axes with different colors are used in the following figures. In Figures10, 15, and 16, the blue time-axis stands for the algorithm with reverse once and with reverse twice. And in Figures 11, 17, and 18, the time-axes with blue color, green color, and red color represent the time axis of the algorithm without reverse, with reverse once, and with reverse twice, respectively.


### 4.3. Analysis of Simulation Results

From Figure11, it can be seen that performance of “subsection + reverse” filtering algorithm is significantly better than that without reverse. Specific analyses are as follows.(1)
From Figures11(a) and 11(b), it can be seen that when reverse filtering is not used, the initial error converges to a certain angle and diverges later with time. At the end of the first subsection (1500 seconds), estimation reaches 92°. When using reverse filtering once, the length of sampling data segment is shortened to 500 seconds that reduces the impact of gyro drift. The precision of initial misalignment angle achieves significant improvement. It reaches 90.5° at the end of the first phase. When reverse filtering is used twice, length of sample data is shortened to 300 seconds that reduce the impact of gyro drift further. Estimation reaches 90.05° in 300 seconds (which is the end of the first stage). Therefore, estimation accuracy is significantly improved; it effectively suppresses the divergence of initial misalignment angle.
(2)
As is shown in Figure11(c), without the introduction of reverse filtering, the estimated effect of the gyro drift is poor. It costs about 6000 seconds to converge. With introduction of reverse filtering once, gyro drift estimation converges in 3500 seconds while the error estimation of gyro drift converges in 1500 seconds after using reverse algorithm twice. The convergence speed becomes significantly faster.
(3)
It can be seen from Figure11(c) that, at the end of the first stage, the initial misalignment is corrected. Because of the introduction of reverse filtering, the model becomes more accurate and the estimation accuracy becomes much higher. Therefore, the precision of cooperative navigation is always within 20 m as shown in Figure 11(d).


## 4.1. Comparison of Filtering Performance with and without Subsection

The whole simulation time is set as 6000 seconds and we take 1500 seconds for the first stage. According to “dual-model algorithm with subsection” proposed in Section4, the curves of estimation are shown in Figures 8 and 9.Figure 8
Estimation comparison with and without subsection.Comparison of filtering performance with and without subsection.
(a)
Initial alignment estimation (enlarge view)
(b)
Gyro drift estimation in long time without subsection
(c)
Colocation error comparisonFrom the comparison in Figures8 and 9, the following conclusions can be drawn.(1)
As can be seen from Figure8, the initial misalignment angle without segmentation begins to converge within a certain range but diverges later (as previously analyzed, observability grows gradually worse over time). The error estimation accuracy is better than that without segment, and the modification at the end of first stage effectively suppresses initial misalignment angle divergence.
(2)
It can be seen in Figures8 and 9(b) that the estimation performance of gyro drift without segmentation is poor. It gradually converges in a very long time while it converges quickly after subsection.
(3)
As can be seen in Figure9(a), the estimation of the initial misalignment angle in subsection 1 is better than that without subsection, but there is still a tendency of divergence.
(4)
Figure9(c) tells us that due to better estimation of the initial misalignment angle and gyro drift after segmented, the performance of colocation is better than nonsegments. However, such divergent trend still exists in colocation estimation error due to certain bias of initial heading angle estimation.In summary, the performance with segment is better than that without segmentation. But there is still space for further improvement especially in how to improve estimation accuracy of the initial misalignment angle of the first phase.


## 4.2. “Subsection + Reverse” Algorithm

The reverse filtering is introduced on the base of subsection in this section. It gives comparison of methods without reverse, reverse once, and reverse twice. The velocity curve of USV varies between 8 m/s and 15 m/s without reverse. It moves in a straight line with acceleration and deceleration. Figure10 shows velocity with reverse once and with reverse twice. The corresponding performance comparison is shown in Figure 11.Velocity of follower USV.
(a)
Velocity with reverse once
(b)
Velocity with reverse twiceComparison with no reverse, reverse once, and reverse twice.
(a)
Estimation comparison in whole process
(b)
The partial enlarged view of initial misalignment
(c)
An enlarged view of the gyro drift estimation
(d)
Cooperative navigation error estimationTo make the description more clearly, several axes with different colors are used in the following figures. In Figures10, 15, and 16, the blue time-axis stands for the algorithm with reverse once and with reverse twice. And in Figures 11, 17, and 18, the time-axes with blue color, green color, and red color represent the time axis of the algorithm without reverse, with reverse once, and with reverse twice, respectively.


## 4.3. Analysis of Simulation Results

From Figure11, it can be seen that performance of “subsection + reverse” filtering algorithm is significantly better than that without reverse. Specific analyses are as follows.(1)
From Figures11(a) and 11(b), it can be seen that when reverse filtering is not used, the initial error converges to a certain angle and diverges later with time. At the end of the first subsection (1500 seconds), estimation reaches 92°. When using reverse filtering once, the length of sampling data segment is shortened to 500 seconds that reduces the impact of gyro drift. The precision of initial misalignment angle achieves significant improvement. It reaches 90.5° at the end of the first phase. When reverse filtering is used twice, length of sample data is shortened to 300 seconds that reduce the impact of gyro drift further. Estimation reaches 90.05° in 300 seconds (which is the end of the first stage). Therefore, estimation accuracy is significantly improved; it effectively suppresses the divergence of initial misalignment angle.
(2)
As is shown in Figure11(c), without the introduction of reverse filtering, the estimated effect of the gyro drift is poor. It costs about 6000 seconds to converge. With introduction of reverse filtering once, gyro drift estimation converges in 3500 seconds while the error estimation of gyro drift converges in 1500 seconds after using reverse algorithm twice. The convergence speed becomes significantly faster.
(3)
It can be seen from Figure11(c) that, at the end of the first stage, the initial misalignment is corrected. Because of the introduction of reverse filtering, the model becomes more accurate and the estimation accuracy becomes much higher. Therefore, the precision of cooperative navigation is always within 20 m as shown in Figure 11(d).


## 5. Water Test Verification

To verify the feasibility of proposed cooperative navigation algorithms further, water experiment is conducted at Tai Lake in Wuxi city of Jiangsu Province. For ease of operation, one self-made USV is used for followers and two patrol crafts are chosen to play the role of leader USVs instead.


### 5.1. Overview of the Test

The equipment installation is as follow. As shown in Figure12, each leader USV is equipped with GPS, acoustic communication devices, and data collection system. GPS is used for providing accurate positioning information, and acoustic communications equipment (Teledyne Benthos’s ATM-885 sonar is selected that can achieve 360° sonar signal transmission and reception) is used to complete information transmission and measurement of distance between the USVs. Follower USV is equipped with PHINS, GPS, DVL, MEMS, and also ATM-885 sonar. PHINS works in combination mode with GPS. It is benchmark that provides high-accurate heading, attitude, and position and speed information. Dead-reckoning is conducted by speed (provided by DVL) and heading (provided by MEMS).Unmanned craft as follower (a) and two patrol crafts as leaders (b).
(a)
(b)The performance of each equipment is as shown in Table2.Table 2
The performance of equipment.
MEMS
Gyroscope bias stability
<30°/h
Gyroscope scale factor error
0.1%
Gyroscope measuring range
±300°/s
Accelerometers bias stability
0.2 mg
Accelerometers threshold
±
5
×
10
-
4 g
Accelerometers measuring range
±10 g
PHINS
Position accuracy (CEP50%)
0.5–3 m
Heading accuracy (1σ value)
0.01° secant latitude
Attitude accuracy (1σ value)
Less than 0.01°
GPS
Time accuracy
1 us
Velocity accuracy
0.1 m/s
Position accuracy
<2.5 m
DVL
Working range
−150 m/s–200 m/s
Measurement accuracy
0.1%
Acoustic device
Working range
Up to 8000 m
Error rate (with correction algorithm)
Less than 10−7Equipment installation and signal flow are as shown in Figure13.Figure 13
The equipment installation in follower (USV C).The test procedure is as follows: (1) two water patrol boats are selected as pilot USV (labeled as A and B) and a USV as follower (labeled as C); the computer clock is synchronized precisely. (2) Dead-reckoning is conducted using speed provided by DVL and heading angle provided by MEMS to obtain the real-time location. (3) The two leaders send their location information to C alternately using acoustic device; the interval is selected as 10 s. (4) After receiving aided information and distance, information fusion is made by using the dead-reckoning information, location information from leaders, and the distance measurement.In the experiment, complete test data are stored to make offline analysis and comparison of different filtering methods. The position information provided by GPS installed on USV C is used as benchmark to make comparison of different methods. Similarly, PHINS on C is used as the benchmark for heading. Shot by Google earth map photo, three USVs trajectories are shown in Figure14 (about 20 kilometers of the whole voyage).Figure 14
Trajectory of three vehicles.The trajectory and speed of follower (reverse once).
(a)
The speed of follower (reverse once)
(b)
Trajectory of follower USV (reverse once)The trajectory and speed of follower (reverse two time).
(a)
Velocity of follower USV (reverse one time)
(b)
Trajectory of follower USV (reverse one time)Estimation comparison.
(a)
Estimation comparison in subsection 1
(b)
Estimation comparison in subsection 2Comparison of the positioning errors.
(a)
Colocating error in subsection 1
(b)
Colocating error of the whole process


### 5.2. Data Processing

We take 1250 seconds for the first stage in which initial misalignment is estimated, and gyro drift estimation is conducted 1250 seconds later. Corresponding curves of three methods (including algorithm without reverse, reverse once, and reverse twice) are given in Figures15 to 18.


### 5.3. Analysis

From the simulation results in Figures17 and 18, discussion is made as follows.(1)
Figure17(a) shows that when backward filtering is not used, the initial misalignment converges to certain angle but becomes divergence later. Estimation reaches 95° at the 1250th second (which is the end of the first phase). When reverse filtering is used once, sampling data length is shortened to 417 seconds and the impact of gyro drift is suppressed sharply. Accuracy of initial misalignment estimation is improved (that is 92° at 417th second). To reduce the impact of gyro drift further, reverse filtering is used twice. Sampling data length is shortened to 250 seconds while estimation becomes 89.9° in 250 seconds. Accuracy of estimation is significantly improved, and the error is corrected in the end of the first subsection.
(2)
From Figure17(b), it can be seen that estimation performance of gyro drift is poor without the introduction of reverse filtering. It gradually converges in 3500 seconds. By using the reverse filtering once, the convergence of gyro drift costs about 1500 seconds while it converges in 500 seconds by using the reverse filtering twice. The convergence speed become much faster.
(3)
Figure18(b) shows the divergent trend of cooperative navigation error in 6000 seconds without inverse algorithm and the maximum error is larger than 500 meters. After the introduction of reverse filtering, cooperative navigation and positioning accuracy is significantly improved as initial misalignment angle error and gyro drift are estimated more accurately. The maximum positioning error is less than 400 meters with reverse once while it is within 300 meters with reverse twice.


## 5.1. Overview of the Test

The equipment installation is as follow. As shown in Figure12, each leader USV is equipped with GPS, acoustic communication devices, and data collection system. GPS is used for providing accurate positioning information, and acoustic communications equipment (Teledyne Benthos’s ATM-885 sonar is selected that can achieve 360° sonar signal transmission and reception) is used to complete information transmission and measurement of distance between the USVs. Follower USV is equipped with PHINS, GPS, DVL, MEMS, and also ATM-885 sonar. PHINS works in combination mode with GPS. It is benchmark that provides high-accurate heading, attitude, and position and speed information. Dead-reckoning is conducted by speed (provided by DVL) and heading (provided by MEMS).Unmanned craft as follower (a) and two patrol crafts as leaders (b).
(a)
(b)The performance of each equipment is as shown in Table2.Table 2
The performance of equipment.
MEMS
Gyroscope bias stability
<30°/h
Gyroscope scale factor error
0.1%
Gyroscope measuring range
±300°/s
Accelerometers bias stability
0.2 mg
Accelerometers threshold
±
5
×
10
-
4 g
Accelerometers measuring range
±10 g
PHINS
Position accuracy (CEP50%)
0.5–3 m
Heading accuracy (1σ value)
0.01° secant latitude
Attitude accuracy (1σ value)
Less than 0.01°
GPS
Time accuracy
1 us
Velocity accuracy
0.1 m/s
Position accuracy
<2.5 m
DVL
Working range
−150 m/s–200 m/s
Measurement accuracy
0.1%
Acoustic device
Working range
Up to 8000 m
Error rate (with correction algorithm)
Less than 10−7Equipment installation and signal flow are as shown in Figure13.Figure 13
The equipment installation in follower (USV C).The test procedure is as follows: (1) two water patrol boats are selected as pilot USV (labeled as A and B) and a USV as follower (labeled as C); the computer clock is synchronized precisely. (2) Dead-reckoning is conducted using speed provided by DVL and heading angle provided by MEMS to obtain the real-time location. (3) The two leaders send their location information to C alternately using acoustic device; the interval is selected as 10 s. (4) After receiving aided information and distance, information fusion is made by using the dead-reckoning information, location information from leaders, and the distance measurement.In the experiment, complete test data are stored to make offline analysis and comparison of different filtering methods. The position information provided by GPS installed on USV C is used as benchmark to make comparison of different methods. Similarly, PHINS on C is used as the benchmark for heading. Shot by Google earth map photo, three USVs trajectories are shown in Figure14 (about 20 kilometers of the whole voyage).Figure 14
Trajectory of three vehicles.The trajectory and speed of follower (reverse once).
(a)
The speed of follower (reverse once)
(b)
Trajectory of follower USV (reverse once)The trajectory and speed of follower (reverse two time).
(a)
Velocity of follower USV (reverse one time)
(b)
Trajectory of follower USV (reverse one time)Estimation comparison.
(a)
Estimation comparison in subsection 1
(b)
Estimation comparison in subsection 2Comparison of the positioning errors.
(a)
Colocating error in subsection 1
(b)
Colocating error of the whole process


## 5.2. Data Processing

We take 1250 seconds for the first stage in which initial misalignment is estimated, and gyro drift estimation is conducted 1250 seconds later. Corresponding curves of three methods (including algorithm without reverse, reverse once, and reverse twice) are given in Figures15 to 18.


## 5.3. Analysis

From the simulation results in Figures17 and 18, discussion is made as follows.(1)
Figure17(a) shows that when backward filtering is not used, the initial misalignment converges to certain angle but becomes divergence later. Estimation reaches 95° at the 1250th second (which is the end of the first phase). When reverse filtering is used once, sampling data length is shortened to 417 seconds and the impact of gyro drift is suppressed sharply. Accuracy of initial misalignment estimation is improved (that is 92° at 417th second). To reduce the impact of gyro drift further, reverse filtering is used twice. Sampling data length is shortened to 250 seconds while estimation becomes 89.9° in 250 seconds. Accuracy of estimation is significantly improved, and the error is corrected in the end of the first subsection.
(2)
From Figure17(b), it can be seen that estimation performance of gyro drift is poor without the introduction of reverse filtering. It gradually converges in 3500 seconds. By using the reverse filtering once, the convergence of gyro drift costs about 1500 seconds while it converges in 500 seconds by using the reverse filtering twice. The convergence speed become much faster.
(3)
Figure18(b) shows the divergent trend of cooperative navigation error in 6000 seconds without inverse algorithm and the maximum error is larger than 500 meters. After the introduction of reverse filtering, cooperative navigation and positioning accuracy is significantly improved as initial misalignment angle error and gyro drift are estimated more accurately. The maximum positioning error is less than 400 meters with reverse once while it is within 300 meters with reverse twice.


## 6. Discussion

Through the above analysis and comparison of experimental data, it can be concluded as follows.(1)
Because dead-reckoning cannot provide precise positioning information in long period, the cooperative navigation algorithm with CKF is introduced and “subsection + reverse” solution is also proposed to improve accuracy effectively while time cost is reduced.
(2)
There are still some problems to be solved in practical applications. For example, ranging accuracy and the communication quality greatly affect the performance of colocation. In speedup, a lot of bubbles appear around underwater acoustic equipment and information cannot be successfully transferred. More stable prediction algorithm or effective compensation is required and the data packet is lost. In addition, there are still some other problems such as uncertain noise and the influence of unknown current.
(3)
Two new questions appear in front of us. The first one is whether the continuously shortened sampling data could result in continuously better performance. The second is whether length of the first stage can be shortened unlimitedly. For these two questions we do two more simulations: the length of the first phase is taken as 150 seconds and 30 seconds separately. The reversion is set as 40 times and 200 times correspondingly. Initial misalignment estimation and colocation error convergence curve are shown in Figure19.Performance of reverse algorithm with shorter period.
(a)
Initial alignment estimation curve
(b)
Colocation error curveFrom Figure19 we can see that the estimation accuracy can be improved by shortening the length of sampling data period (such the case of 150 seconds), but the improvement is not that significant. When length of sampling data becomes less than certain value (e.g., 30 seconds) the estimation accuracy of the initial misalignment cannot be improved but be depressed. Then in the second subsection, divergent tendency appears in the estimation of gyro drift and position error. It is mainly because when time period is too short, effect of noise would be amplified rapidly. Therefore, it is very important to make further study to select the data length more reasonably.


## 7. Summary

In this paper, we research on cooperative navigation based on dead-reckoning with MEMS/DVL. The nonlinear system model of cooperative navigation is constructed and CKF is proposed to make information fusion. Aiming at the coexistence of large misalignment and large gyro drift, the observability analysis is carried out. On the basis, dual-model estimation idea is introduced to cut off the coupling relations and improve the observability. Furthermore, reverse filtering estimation method is proposed. The result of simulation and experiments show that it can effectively improve the navigation algorithm performance. The problems in the experiment are also stated and future research directions are discussed.



---

*Source: 186785-2014-11-06.xml*