A New Continuous Rotation IMU Alignment Algorithm Based on Stochastic Modeling for Cost Effective North-Finding Applications

Abstract

Based on stochastic modeling of Coriolis vibration gyros by the Allan variance technique, this paper discusses Angle Random Walk (ARW), Rate Random Walk (RRW) and Markov process gyroscope noises which have significant impacts on the North-finding accuracy. A new continuous rotation alignment algorithm for a Coriolis vibration gyroscope Inertial Measurement Unit (IMU) is proposed in this paper, in which the extended observation equations are used for the Kalman filter to enhance the estimation of gyro drift errors, thus improving the north-finding accuracy. Theoretical and numerical comparisons between the proposed algorithm and the traditional ones are presented. The experimental results show that the new continuous rotation alignment algorithm using the extended observation equations in the Kalman filter is more efficient than the traditional two-position alignment method. Using Coriolis vibration gyros with bias instability of 0.1°/h, a north-finding accuracy of 0.1° (1σ) is achieved by the new continuous rotation alignment algorithm, compared with 0.6° (1σ) north-finding accuracy for the two-position alignment and 1° (1σ) for the fixed-position alignment.

1. Introduction

Cost effective north-finding technology is widely required for many applications. North-finding is sometimes based on Digital Magnetic Compasses (DMCs) [1]. However, DMCs is easily degraded by magnetic interference. Although Dynamically Tuned Gyros and Ring Laser Gyroscopes are suitable for precise north-finding, they are generally bulky and expensive [2,3]. In contrast, Coriolis vibration gyroscopes (e.g., a kind of cost effective medium precision Hemispherical Resonator Gyroscopes (HRGs) [4,5]) are generally compact and low-cost and suitable for a cost effective north-finding system. However, the drift errors of these gyroscopes are big problems, which limit the north-finding accuracy.

To improve the accuracy of the north-finding system using cost effective gyroscopes, several methods have been designed. Lee [6] proposed a multi-position alignment algorithm to increase the azimuth accuracy. For the same purpose, Yu [7] used analytic optimization of Strapdown Inertial Navigation System (SINS) multi-position alignment. Renkoski [8] and Sun [9] improved the accuracy of North-finding system through continuous rotation.

This paper focuses on Inertial Measurement Unit (IMU)-based north-finding systems using a Kalman filter for applications such as dynamic orientation and dead reckoning. Stochastic modeling for a Coriolis vibration gyroscope is obtained using the Allan variance technique. It is shown that the Rate Random Walk (RRW) and Markov noises are the main errors which limit the north-finding accuracy. A new continuous rotation IMU alignment algorithm is therefore proposed using extended observation equations in the Kalman filter to solve this problem. Experimental results as well as theoretical analysis are also presented.

This paper is organized as follows: Section 2 analyses the random error model of a Coriolis vibration gyroscope using the Allan variance technique. The north-finding errors due to the main parts of the gyro drift error are presented. Section 3 presents three different IMU based north-finding algorithms or three different error compensation approaches: two-position alignment, continuous rotation alignment, and a new continuous rotation alignment algorithm with extended observation equations for a Kalman filter. Section 4 presents theoretical and simulation analyses of the performances of the methods mentioned above. Section 5 reports north-finding experimental results and comparisons. The Allan variance analysis results for the equivalent east gyro are presented for the interpretation of effectiveness of the gyro drift error compensation approaches. Section 6 concludes the paper. The appendices show detailed theoretical proofs.

2. Error Model for a Coriolis Vibration Gyroscope

IMU errors can be classified into two types: deterministic errors and random errors. Major deterministic error sources including constant bias, scale factor errors and misalignment can be removed by calibration and compensation [10]. The random constant bias (turn to turn bias) and random noises are the main error sources in the North-finding system. Therefore, we focus on the stochastic modeling for a Coriolis vibration gyroscope.

2.1. Error Model Based on Allan Variance Analysis

Traditionally, random constant bias, ARW (Angle Random Walk), RRW and Markov process are used to develop stochastic error model for gyros. The error model of a gyroscope can be expressed as follows [11,12]:

$$\epsilon ={\epsilon }_b+{\epsilon }_m+w_a+{\epsilon }_r$$

(1)

where $\epsilon$ is the stochastic drift error of the gyroscope measurements, ${\epsilon }_b$ is the random constant bias with the variance of ${\sigma }_b^2$, ${\epsilon }_m$ is the Markov process, $w_a$ is the ARW, ${\epsilon }_r$ is the RRW.

The random bias can be described as an unpredictable random quantity with a constant value, that is:

$${\dot{\epsilon }}_b=0$$

(2)

$$w_a\in N(0,{\sigma }_a^2)$$

(3)

where ${\sigma }_a^2$ is the variance of $w_a$.

The Markov noise is the low-frequency component in the error sources. Usually, the noise is modeled as a First order Gauss-Markov process [11]:

$${\dot{\epsilon }}_m=-\frac{1}{\tau }{\epsilon }_m+w_m,w_m\in N(0,{\sigma }_m^2)$$

(4)

where $\tau$ is the process time constant, $w_m$ is the zero-mean Gaussian white noise, ${\sigma }_m^2$ is the variance of $w_m$:

$${\dot{\epsilon }}_r=w_r,w_r\in N(0,{\sigma }_r^2)$$

(5)

where ${\sigma }_r^2$ is the variance of $w_r$.

In Equation (1), the characteristics of the stochastic errors are usually estimated by an optimal estimation algorithm, such as a Kalman filter [13]. The parameters of the stochastic error model are necessary for a Kalman filter algorithm. Hence, there is a need to determine the parameters of the error model using Allan variance analysis. The sampling data of a HRG in 3 h is present in Figure 1a. The Allan variance results of the HRG are presented in Figure 1b. The sampling frequency is 10 Hz.

The parameters of the error models for the Coriolis vibration gyroscopes in an IMU based north-finding system are given in

Table 1. Parameters of the error models for the Coriolis vibration gyroscopes (a kind of cost effective HRG in this paper).

Bias Instability ${\sigma }_b$	$0.1°/h$
ARW ${\sigma }_a$	$0.01°/\sqrt{h}$
RRW ${\sigma }_r$	$0.3°/h^{3/2}$
Markov time constant $\tau$	$60s$
Markov process driving noise ${\sigma }_m$	$0.02°/h/\sqrt{s}$

.

Consider the error models in Figure 1, the major parts of the gyroscope errors are ARW, Markov process, bias instability and RRW, which indicates that the error model in Equation (1) is sufficient to characterize the gyroscope. The parameters of the models show that the primary error source for the gyroscope are Markov noise and RRW.

2.2. Propagation of Gyroscope Errors in a North-Finding System

The drift error of the equivalent east gyroscope ${\epsilon }_E$ in an IMU based north-finding system propagates to the azimuth misalignment ${\varphi }_D$, which can be expressed as follows [14]:

$${\varphi }_D=\frac{{\epsilon }_E}{\mathrm{\Omega }\cos L}$$

(6)

where $\mathrm{\Omega }$ is the earth rotation rate, $L$ is the local latitude.

Similar to Equation (1), ${\epsilon }_E$ can be expressed as follows:

$${\epsilon }_E={\epsilon }_{bE}+w_{aE}+{\epsilon }_{rE}+{\epsilon }_{mE}$$

(7)

where the random constant bias ${\epsilon }_{bE}$, the ARW $w_{aE}$, the RRW ${\epsilon }_{rE}$ and the Markov process ${\epsilon }_{mE}$ correspond to ${\epsilon }_b$, $w_a$, ${\epsilon }_r$ and ${\epsilon }_m$ in Equation (1).

The RMS (Root Mean Square) of azimuth misalignment ${\sigma }_{{\varphi }_{Db}}$, ${\sigma }_{{\varphi }_{Da}}$, ${\sigma }_{{\varphi }_{Dr}}$ and ${\sigma }_{{\varphi }_{Dm}}$ due to ${\epsilon }_{bE}$, $w_{aE}$, ${\epsilon }_{rE}$ and ${\epsilon }_{mE}$ can be expressed as [15,16]:

$${\sigma }_{{\varphi }_{Db}}=\frac{{\sigma }_{bE}}{\mathrm{\Omega }\cos L}$$

(8)

$${\sigma }_{{\varphi }_{Da}}=\frac{{\sigma }_{aE}}{\sqrt{t_n}\mathrm{\Omega }\cos L}$$

(9)

$${\sigma }_{{\varphi }_{Dr}}=\frac{{\sigma }_{rE}\sqrt{t_n}}{\sqrt{3}\mathrm{\Omega }\cos L}$$

(10)

$${\sigma }_{{\varphi }_{Dm}}=\frac{\sqrt{\frac{1}{2}{\tau }^2{\sigma }_{mE}^2(2t-\tau e^{-\frac{2t}{\tau }}+4\tau e^{-\frac{t}{\tau }}-3\tau )+P_{{\epsilon }_{mE}}(0)\tau (1-2e^{-\frac{1}{\tau }t}+e^{-\frac{2}{\tau }t})}}{t_n\mathrm{\Omega }\cos L}$$

(11)

where $t_n$ is the alignment time. ${\sigma }_{bE}^2$, ${\sigma }_{aE}^2$, ${\sigma }_{rE}^2$ and ${\sigma }_{mE}^2$ are the variances corresponding to ${\epsilon }_{bE}$, $w_{aE}$, $w_{rE}$ and $w_{mE}$ in the error model of the equivalent east gyroscope. $P_{{\epsilon }_{mE}}(0)$ is the variance of the initial value of Markov process. The proofs of Equations (8)–(11) are shown in Appendix A.

It should be explained that the initial value of RRW noise can be regarded as part of a constant bias. Thus the RRW starts from zero.

Assuming the alignment time $t_n$ is 10 min, the local latitude is 28.22° N, the RMS values of the azimuth misalignment can be obtained from Equations (8)–(11). The azimuth misalignment due to the equivalent east gyroscope errors are shown in

Table 2. The azimuth misalignment due to the equivalent east gyroscope errors in 10 min at 28.22° N.

Gyroscope Errors		RMS of Azimuth Misalignment
Bias Instability ${\sigma }_{bE}$	$0.1°/h$	${\sigma }_{{\varphi }_{Db}}=0.43°$
ARW ${\sigma }_{aE}$	$0.01°/\sqrt{h}$	${\sigma }_{{\varphi }_{Da}}=0.10°$
RRW ${\sigma }_{rE}$	$0.3°/h^{3/2}$	${\sigma }_{{\varphi }_{Dr}}=0.31°$
Markov process $\tau$	$60s$	${\sigma }_{{\varphi }_{Dm}}=0.20°$
Markov process ${\sigma }_{mE}$	$0.02°/h/\sqrt{s}$

.

Although the azimuth misalignment are most affected by the bias instability, the random constant bias can be easily eliminated through north-finding algorithms (such as two-position alignment [6] and continuous rotation alignment [9]). And compared with RRW and Markov noise, the azimuth misalignment due to ARW is slim. RRW and Markov process are the main error source in a north-finding system.

3. Error Compensation Approach for IMU Based North-Finding System

3.1. System Error Model for IMU Based North-Finding

A local level NED (North-East-Down) frame is used as the navigation frame. The common SINS error equations in the navigation frame can be expressed as follows [14]:

$${\dot{\mathbf{\varphi }}}^n=-({\omega }_{en}^n+{\omega }_{ie}^n)\times {\mathbf{\varphi }}^n-C_b^n\delta {\omega }_{ib}^b$$

(12)

$$\delta {\dot{v}}^n=f^n\times {\mathbf{\varphi }}^n-({\omega }_{en}^n+2{\omega }_{ie}^n)\times \delta v^n+C_b^n\delta f^b$$

(13)

where ${\mathbf{\varphi }}^n$ is the attitude error, ${\mathbf{\varphi }}^n={\left[\begin{array}{ccc}{\varphi }_N& {\varphi }_E& {\varphi }_D\end{array}\right]}^T$, $N$, $E$ and $D$ represent north, east and down in navigation frame respectively; $\delta v^n$ is the velocity error, $\delta v^n={\left[\begin{array}{ccc}\delta v_N& \delta v_E& \delta v_D\end{array}\right]}^T$. ${\mathbf{\varphi }}^n$ can be estimated by the observation of $\delta v^n$ in an alignment process. $f^n$ is the measurement of specific force in frame n, $C_b^n$ is the coordinate transformation matrix from the IMU frame b to the navigation frame n, ${\omega }_{en}^n$ is the turn rate of the navigation frame to the earth frame in the frame n, ${\omega }_{ie}^n$ is the turn rate of the earth frame to the inertial frame in the frame n, $\delta {\omega }_{ib}^b$ is the error of the gyroscope measurements, $\delta f^b$ is the error of the specific force measurements.

In the North finding scenario discussed here, since the IMU is stationary on the Earth:

$${\omega }_{en}^n=0$$

(14)

The SINS error model for alignment or IMU based north-finding can be written as:

$$\dot{x}(t)=\left[\begin{array}{cc}{F}_1& F_2\\ 0_{8\times 5}& \mathrm{\Gamma }\end{array}\right]\dot{x}(t)+\left[\begin{array}{cc}{G}_1& 0_{5\times 3}\\ 0_{8\times 5}& \begin{array}{c}{0}_{5\times 3}\\ I_{3\times 3}\end{array}\end{array}\right]w(t)$$

(15)

where:

$$\begin{array}{l}x(t)=[\begin{array}{ccccc}\delta v_N& \delta v_E& {\varphi }_N& {\varphi }_E& \begin{array}{cc}{\varphi }_D& \begin{array}{cc}{\nabla }_{accx}& {\nabla }_{accy}\end{array}\end{array}\end{array}\\ {\begin{array}{cccccc}({\epsilon }_{bx}+{\epsilon }_{rx})& ({\epsilon }_{by}+{\epsilon }_{ry})& ({\epsilon }_{bz}+{\epsilon }_{rz})& {\epsilon }_{mx}& {\epsilon }_{my}& {\epsilon }_{mz}\end{array}]}^T\end{array}$$

(16)

${\nabla }_{accx}$ and ${\nabla }_{accy}$ are the bias error states of the accelerometers, ${\epsilon }_{bx}$, ${\epsilon }_{by}$ and ${\epsilon }_{bz}$ are the random constant bias error states of the gyroscopes, ${\epsilon }_{rx}$, ${\epsilon }_{ry}$ and ${\epsilon }_{rz}$ are the rate random walk of the gyroscopes, ${\epsilon }_{mx}$, ${\epsilon }_{my}$ and ${\epsilon }_{mz}$ are the error states for the Markov process of the gyroscopes.

For the filter noise vector:

$$w(t)={\left[\begin{array}{ccccc}{w}_{accx}& w_{accy}& w_{ax}& w_{ay}& \begin{array}{ccccccc}{w}_{az}& w_{rx}& w_{ry}& w_{rz}& w_{mx}& w_{my}& w_{mz}\end{array}\end{array}\right]}^T$$

(17)

where $w_{accx}$ and $w_{accy}$ are the white noise of the accelerometer x and the accelerometer y, respectively. That is:

$$w_{accx},w_{accy}\in N(0,{\sigma }_{acc}^2)$$

(18)

where ${\sigma }_{acc}^2$ is the variance of the white noise $w_{accx}$ and $w_{accy}$.

$w_{ax}$, $w_{ay}$ and $w_{az}$ are the angular random walk of the gyroscope x, the gyroscope y and the gyroscope z, $w_{mx}$, $w_{my}$ and $w_{mz}$ are the driving noise in the Markov process of the gyroscope x, the gyroscope y and the gyroscope z.

$$F_1=\left[\begin{array}{ccccc}0& -2\mathrm{\Omega }\sin L& 0& g& 0\\ 2\mathrm{\Omega }\sin L& 0& -g& 0& 0\\ 0& 0& 0& -\mathrm{\Omega }\sin L& 0\\ 0& 0& \mathrm{\Omega }\sin L& 0& \mathrm{\Omega }\cos L\\ 0& 0& 0& -\mathrm{\Omega }\cos L& 0\end{array}\right]$$

(19)

where $g$ is the local gravity.

The matrix $F_2$ is defined as follows:

$$F_2=\left[\begin{array}{ccc}-{C_b^n}_{2\times 2}& 0_{2\times 3}& 0_{2\times 3}\\ 0_{3\times 2}& -C_b^n& -C_b^n\end{array}\right]$$

(20)

where ${C_b^n}_{2\times 2}$ is defined as:

$${C_b^n}_{2\times 2}=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\end{array}\right]C_b^n{\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\end{array}\right]}^T$$

(21)

The matrices $G_1$ and $\mathrm{\Gamma }$ are defined as follows:

$$G_1=\left[\begin{array}{cc}{C_b^n}_{2\times 2}& 0_{2\times 3}\\ 0_{3\times 2}& -C_b^n\end{array}\right],\mathrm{\Gamma }=\left[\begin{array}{cc}{0}_{5\times 5}& 0_{5\times 3}\\ 0_{3\times 5}& -\frac{1}{\tau }{I}_{3\times 3}\end{array}\right]$$

(22)

where $\tau$ is the Markov time constant of the gyroscope.

As shown in the analysis above, based on the condition that the system is stationary on the earth, the horizontal velocity errors are used as observation states. Thus, the observation model can be written as:

$$z(t)=\left[\begin{array}{c}{\tilde{v}}_N\\ {\tilde{v}}_E\end{array}\right]=Hx(t)+\mathbf{\upsilon }(t)H=\left[\begin{array}{cc}{I}_{2\times 2}& 0_{2\times 11}\end{array}\right]$$

(23)

where $\mathbf{\upsilon }(t)={\left[\begin{array}{cc}{\upsilon }_{v_N}& {\upsilon }_{v_E}\end{array}\right]}^T$ is the observation noise vector. ${\tilde{v}}_N$ and ${\tilde{v}}_E$ represent north and east components of the estimated velocity, respectively.

3.2. Traditional Two-Position Gyrocompassing

Two-position alignment is demonstrated in Figure 2 [6].

As shown in Figure 2, the axis $x_b$ and $y_b$ of the IMU frame lie on the turntable plane, the axis $z_b$ coincides with the rotation axis. We define the $b_0$ frame when $x_b$ coincides with the turntable null indicator:

$$C_n^b=C_{b_0}^b{C}_n^{b_0}$$

(24)

where $C_n^{b_0}$ is the coordinate transformation matrix from the frame $n$ to the frame $b_0$.

$C_{b_0}^b$ can be written as:

$$C_{b_0}^b(t)=\{\begin{array}{c}\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right]\begin{array}{c}t<t_1\end{array}\\ \left[\begin{array}{ccc}-1& 0& 0\\ 0& -1& 0\\ 0& 0& 1\end{array}\right]\begin{array}{c}t>t_2\end{array}\end{array}$$

(25)

where $\left[t_1,t_2\right]$ is the short time period when the IMU changes the angular position through the turntable rotation.

3.3. Continuous Rotation Gyrocompassing

As an alternative to the two-position alignment, continuous rotation is another efficient method to reduce the alignment errors.

In contrast to the two-position alignment, the coordinate transformation matrix $C_b^n$ is varying as $C_b^{b_0}$ changes by continuous rotation, that is:

$$C_b^{b_0}=\left[\begin{array}{ccc}\cos ({\omega }_{0}t)& \sin ({\omega }_{0}t)& 0\\ -\sin ({\omega }_{0}t)& \cos ({\omega }_{0}t)& 0\\ 0& 0& 1\end{array}\right]{\omega }_{0}T=2\pi$$

(26)

where ${\omega }_0$ is the rotation rate of the turntable. T is the rotation cycle.

Except for the coordinate transformation matrix $C_b^n$ and $C_b^{b_0}$, the error model and the observation equation between the continuous rotation gyrocompassing are the same as that of the two-position gyrocompassing.

3.4. A New Continuous Rotation North-Finding Method Based on an Extended Observation Model

Although the constant random biases of gyroscopes are mostly eliminated by the above compensation approaches, the noise of the gyroscopes will also still affect the efficiency of the Kalman filter. For Coriolis vibration gyroscopes, the noise level is high. It is difficult to estimate the drift errors of the gyroscopes exactly. The accuracy of the North-finding system is limited. To solve the problem, we present an extended observation model for the continuous rotation alignment.

After each $360°$ turn, the integration of the measurements of the gyroscopes can be written as:

$${\displaystyle {\int }_t^{t+T}{\tilde{\mathbf{\omega }}}_{ib}^{b}dt}={\displaystyle {\int }_t^{t+T}{\mathbf{\omega }}_{ib}^{b}dt+}{\displaystyle {\int }_t^{t+T}\left[\begin{array}{c}({\epsilon }_{bx}+{\epsilon }_{rx})+{\epsilon }_{mx}\\ ({\epsilon }_{by}+{\epsilon }_{ry})+{\epsilon }_{my}\\ ({\epsilon }_{bz}+{\epsilon }_{rz})+{\epsilon }_{mz}\end{array}\right]dt}\approx {\displaystyle {\int }_t^{t+T}\left({\mathbf{\omega }}_{eb}^b+C_{b_0}^b{C}_n^{b_0}{\mathbf{\omega }}_{ie}^n\right)dt}+T\left[\begin{array}{c}({\epsilon }_{bx}+{\epsilon }_{rx})+{\epsilon }_{mx}\\ ({\epsilon }_{by}+{\epsilon }_{ry})+{\epsilon }_{my}\\ ({\epsilon }_{bz}+{\epsilon }_{rz})+{\epsilon }_{mz}\end{array}\right]$$

(27)

While the integration of the estimated measurements of the gyroscopes can be written as:

$$\begin{array}{l}{\displaystyle {\int }_t^{t+T}{\widehat{\mathbf{\omega }}}_{ib}^{b}dt}={\displaystyle {\int }_t^{t+T}\left({\mathbf{\omega }}_{eb}^b+{\widehat{\mathbf{\omega }}}_{ie}^b\right)dt}\\ ={\displaystyle {\int }_t^{t+T}{\mathbf{\omega }}_{eb}^{b}dt}+{\displaystyle {\int }_t^{t+T}{\tilde{C}}_n^b{\widehat{\mathbf{\omega }}}_{ie}^{n}dt}\\ ={\displaystyle {\int }_t^{t+T}{\mathbf{\omega }}_{eb}^{b}dt}+{\displaystyle {\int }_t^{t+T}\{C_{b_0}^b{C}_n^{b_0}[I+{\mathbf{\varphi }}^n\times ]{\mathbf{\omega }}_{ie}^n+C_{b_0}^b{C}_n^{b_0}\delta {\mathbf{\omega }}_{ie}^n\}dt}\\ ={\displaystyle {\int }_t^{t+T}\left({\mathbf{\omega }}_{eb}^b+C_{b_0}^b{C}_n^{b_0}{\mathbf{\omega }}_{ie}^n\right)dt}+{\displaystyle {\int }_t^{t+T}\{C_{b_0}^b{C}_n^{b_0}[-{\mathbf{\omega }}_{ie}^n\times ]{\mathbf{\varphi }}^n+C_{b_0}^b{C}_n^{b_0}\delta {\mathbf{\omega }}_{ie}^n\}dt}\end{array}$$

(28)

$$\begin{array}{c}{\mathbf{\omega }}_{ie}^n=\left[\begin{array}{c}\mathrm{\Omega }\cos L\\ 0\\ -\mathrm{\Omega }\sin L\end{array}\right]\delta {\mathbf{\omega }}_{ie}^n=\left[\begin{array}{c}-\delta L\mathrm{\Omega }\sin L\\ 0\\ -\delta L\mathrm{\Omega }\cos L\end{array}\right]\\ C_n^{b_0}=\left[\begin{array}{ccc}\cos {\theta }_0\cos {\phi }_0& \cos {\theta }_0\sin {\phi }_0& -\sin {\theta }_0\\ -\cos {\gamma }_0\sin {\phi }_0+\sin {\gamma }_0\sin {\theta }_0\cos {\phi }_0& \cos {\theta }_0\cos {\phi }_0+\sin {\gamma }_0\sin {\theta }_0\sin {\phi }_0& \sin {\gamma }_0\cos {\theta }_0\\ \cos {\gamma }_0\sin {\theta }_0\cos {\phi }_0+\sin {\gamma }_0\sin {\phi }_0& -\sin {\gamma }_0\cos {\phi }_0+\cos {\gamma }_0\sin {\theta }_0\sin {\phi }_0& \cos {\gamma }_0\cos {\theta }_0\end{array}\right]\end{array}$$

(29)

where ${\displaystyle {\int }_t^{t+T}{\tilde{\mathbf{\omega }}}_{ib}^{b}dt}$ represents the integration of the gyroscope measurements in a rotation cycle of the turntable, ${\widehat{\mathbf{\omega }}}_{ib}^b$ represents the estimated measurements of the gyroscopes in the b-frame, ${\mathbf{\omega }}_{ie}^n$ represents the earth rotation rate in the n-frame. ${\phi }_0$, ${\theta }_0$ and ${\gamma }_0$ are the Euler angles of the $b_0$-frame relative to the n-frame. ${\tilde{C}}_n^b$ is the estimated coordinate transformation matrix with attitude errors.

Considering ${\varphi }_N$ and ${\varphi }_E$ are very small after coarse alignment:

$${\mathbf{\varphi }}^n=\left[\begin{array}{c}{\varphi }_N\\ {\varphi }_E\\ {\varphi }_D\end{array}\right]\approx \left[\begin{array}{c}0\\ 0\\ {\varphi }_D\end{array}\right]$$

(30)

From Equations (29) and (30)

$${\displaystyle {\int }_t^{t+T}{C}_{b_0}^b{C}_n^{b_0}{\mathbf{\omega }}_{ie}^{n}dt}=\mathrm{\Omega }T\left[\begin{array}{c}0\\ 0\\ \cos L(\cos {\gamma }_0\sin {\theta }_0\cos {\phi }_0+\sin {\gamma }_0\sin {\phi }_0)-\sin L\cos {\gamma }_0\cos {\theta }_0\end{array}\right]$$

(31)

$$\begin{array}{l}{\displaystyle {\int }_t^{t+T}\delta {\widehat{\mathbf{\omega }}}_{ib}^{b}dt}\\ ={\displaystyle {\int }_t^{t+T}\{C_{b_0}^b{C}_n^{b_0}[-{\mathbf{\omega }}_{ie}^n\times ]{\mathbf{\varphi }}^n+C_{b_0}^b{C}_n^{b_0}\delta {\mathbf{\omega }}_{ie}^n\}dt}\\ ={\varphi }_D\mathrm{\Omega }T\cos L\left[\begin{array}{c}0\\ 0\\ \left(\begin{array}{l}-\sin {\gamma }_0\cos {\phi }_0+\\ \cos {\gamma }_0\sin {\theta }_0\sin {\phi }_0\end{array}\right)\end{array}\right]+\delta L\mathrm{\Omega }T\left[\begin{array}{c}0\\ 0\\ \left(\begin{array}{l}-\sin L(\cos {\gamma }_0\sin {\theta }_0\cos {\phi }_0+\sin {\gamma }_0\sin {\phi }_0)-\\ \cos L\cos {\gamma }_0\cos {\theta }_0\end{array}\right)\end{array}\right]\end{array}$$

(32)

Under static conditions, we have:

$${\displaystyle {\int }_t^{t+T}{\mathbf{\omega }}_{eb}^{b}dt}=\left[\begin{array}{c}0\\ 0\\ 2\pi \end{array}\right]$$

(33)

Substituting Equations (28), (31)–(33) into Equation (27) gives:

$$\begin{array}{l}T\left[\begin{array}{c}({\epsilon }_{bx}+{\epsilon }_{rx})+{\epsilon }_{mx}\\ ({\epsilon }_{by}+{\epsilon }_{ry})+{\epsilon }_{my}\\ ({\epsilon }_{bz}+{\epsilon }_{rz})+{\epsilon }_{mz}\end{array}\right]-{\varphi }_D\mathrm{\Omega }T\cos L\left[\begin{array}{c}0\\ 0\\ \left(\begin{array}{l}-\sin {\gamma }_0\cos {\phi }_0+\\ \cos {\gamma }_0\sin {\theta }_0\sin {\phi }_0\end{array}\right)\end{array}\right]-\delta L\mathrm{\Omega }T\left[\begin{array}{c}0\\ 0\\ \left(\begin{array}{l}-\sin L\cos {\gamma }_0\sin {\theta }_0\cos {\phi }_0\\ -\sin L\sin {\gamma }_0\sin {\phi }_0\\ -\cos L\cos {\gamma }_0\cos {\theta }_0\end{array}\right)\end{array}\right]\\ ={\displaystyle {\int }_t^{t+T}{\tilde{\mathbf{\omega }}}_{ib}^{b}dt}-{\displaystyle {\int }_t^{t+T}{\widehat{\mathbf{\omega }}}_{ib}^{b}dt}\\ ={\displaystyle {\int }_t^{t+T}{\tilde{\mathbf{\omega }}}_{ib}^{b}dt}-\mathrm{\Omega }T\left[\begin{array}{c}0\\ 0\\ (\cos {\gamma }_0\sin {\theta }_0\cos {\widehat{\phi }}_0+\sin {\gamma }_0\sin {\widehat{\phi }}_0)\cos \widehat{L}-\cos {\gamma }_0\cos {\theta }_0\sin \widehat{L}\end{array}\right]-\left[\begin{array}{c}0\\ 0\\ 2\pi \end{array}\right]\end{array}$$

(34)

When there is latitude error and heading error, the estimated measurements of the gyroscopes are inaccurate. After each $360°$ turn of the turntable, the equivalent east gyroscope error caused by these errors can be calculated as follows:

$$\delta {\omega }_{ibE}^n=\frac{1}{T}\left[\begin{array}{ccc}0& 1& 0\end{array}\right]{\displaystyle {\int }_t^{t+T}{C}_{b_0}^n{C}_b^{b_0}\delta {\widehat{\mathbf{\omega }}}_{ib}^{b}dt}$$

(35)

The equivalent east gyroscope error caused by heading error and latitude error is shown in Equations (36) and (37) respectively:

$$\delta {\omega }_{ibE,{\varphi }_D}^n={(-\sin {\gamma }_0\cos {\phi }_0+\cos {\gamma }_0\sin {\theta }_0\sin {\phi }_0)}^2{\varphi }_D\mathrm{\Omega }\cos L$$

(36)

$$\begin{array}{cc}\hfill \delta {\omega }_{ibE,\delta L}^n& =(\sin {\gamma }_0\cos {\phi }_0-\cos {\gamma }_0\sin {\theta }_0\sin {\phi }_0)[\cos {\gamma }_0\sin {\theta }_0\cos {\phi }_0\tan L\hfill \\ & +\sin {\gamma }_0\sin {\phi }_0\tan L+\cos {\gamma }_0\cos {\theta }_0]\delta L\mathrm{\Omega }\cos L\hfill \end{array}$$

(37)

where $\delta {\omega }_{ibE,{\varphi }_D}^n$ is the equivalent east gyroscope error caused by heading error, $\delta {\omega }_{ibE,\delta L}^n$ is the equivalent east gyroscope error caused by latitude error $\delta L$.

Assuming that:

$${\gamma }_0={\theta }_0=5°$$

(38)

Equations (36) and (37) can be written as:

$$\delta {\omega }_{ibE,{\varphi }_D}^n\approx 0.01{\varphi }_D\mathrm{\Omega }\cos L,\delta {\omega }_{ibE,\delta L}^n\approx 0.1\delta L\mathrm{\Omega }\cos L$$

(39)

In general, the initial heading error is less than 5° (${\varphi }_D(0)<5°$) and the latitude error is less than 0.1° ($\delta L<0.1°$). Considering Equation (39), the equivalent azimuth error caused by initial heading error and latitude error can be ignored when tilt is smaller than 5°.

The additional observation can be obtained using the integration measurements of the gyroscopes in each $360°$ turn of the turntable.

The observation model can be written as:

$$z(t)=Hx(t)+\mathbf{\upsilon }(t)$$

(40)

$$\begin{array}{l}z(t)=\{\begin{array}{c}\left[\begin{array}{c}{\tilde{v}}_N\\ {\tilde{v}}_E\end{array}\right]\begin{array}{c}when2(k-1)\pi <{\omega }_{0}t<2k\pi \end{array}\\ \left[\begin{array}{c}{\tilde{v}}_N\\ {\tilde{v}}_E\\ {\displaystyle {\int }_t^{t+T}{\tilde{\mathbf{\omega }}}_{ib}^{b}dt}-{\displaystyle {\int }_t^{t+T}{\widehat{\mathbf{\omega }}}_{ib}^{b}dt}\end{array}\right]\begin{array}{c}when{\omega }_{0}t=2k\pi \end{array}\end{array}\\ k=1,2,3,\dots \end{array}$$

(41)

$$\begin{array}{l}H=\{\begin{array}{c}\left[\begin{array}{ccc}{I}_{2\times 2}& 0_{2\times 3}& 0_{2\times 8}\end{array}\right]\begin{array}{c}when2(k-1)\pi <{\omega }_{0}t<2k\pi \end{array}\\ \left[\begin{array}{ccccc}{I}_{2\times 2}& 0_{2\times 3}& 0_{2\times 2}& 0_{2\times 3}& 0_{2\times 3}\\ 0_{3\times 2}& D& 0_{3\times 2}& T{I}_{3\times 3}& T{I}_{3\times 3}\end{array}\right]\begin{array}{c}when{\omega }_{0}t=2k\pi \end{array}\end{array}\\ k=1,2,3,\dots \\ D=\left[\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& 0& (-\sin {\gamma }_0\cos {\phi }_0+\cos {\gamma }_0\sin {\theta }_0\sin {\phi }_0)\mathrm{\Omega }T\cos L\end{array}\right]\end{array}$$

(42)

$$\begin{array}{l}\mathbf{\upsilon }(t)=\{\begin{array}{c}{\left[\begin{array}{cc}{\upsilon }_{v_N}& {\upsilon }_{v_E}\end{array}\right]}^T\begin{array}{c}when2(k-1)\pi <{\omega }_{0}t<2k\pi \end{array}\\ \begin{array}{cc}{\left[\begin{array}{ccccc}{\upsilon }_{v_N}& {\upsilon }_{v_E}& {\upsilon }_{{\omega }_x}& {\upsilon }_{{\omega }_y}& {\upsilon }_{{\omega }_z}\end{array}\right]}^T& when{\omega }_{0}t=2k\pi \end{array}\end{array}\\ k=1,2,3,\dots \end{array}$$

(43)

where ${\upsilon }_{{\omega }_x}$, ${\upsilon }_{{\omega }_y}$ and ${\upsilon }_{{\omega }_z}$ are the observation noise corresponding to ${\displaystyle {\int }_t^{t+T}{\tilde{\mathbf{\omega }}}_{ib}^{b}dt}-{\displaystyle {\int }_t^{t+T}{\widehat{\mathbf{\omega }}}_{ib}^{b}dt}$.

4. Comparisons of the Kalman Filter Convergence Rapidity and North-Finding Accuracy

Comparisons of the Kalman filter convergence rapidity and the north-finding accuracy between the proposed algorithms and the traditional alignment methods can be made with the covariance matrix for the estimated states in the Kalman filter.

For the piecewise constant time varying system the covariance matrix of the estimated states $P$ can be obtained by calculating the discrete Riccati matrix equation [7]:

$$\begin{array}{l}{P}^{-1}(k)={[{\mathrm{\Phi }}^T(k,k-1)P(k-1)\mathrm{\Phi }(k,k-1)+G^{T}QG]}^{-1}\\ +H^T{R}^{-1}Hk=1,2,3\dots ,n\end{array}$$

(44)

which is based on the continuous system error model and observation equations (Equations (15)–(43)) as follows:

$$\begin{array}{l}\mathrm{\Phi }(k,k-1)\approx e^{A(t_{k-1})T_s}\approx I+A(t_{k-1})T_s\\ G(k,k-1)\approx B(t_{k-1})\\ Q=q{T}_s\\ q=E\{w^T(t)w(t)\}\\ R=E\{{\mathbf{\upsilon }}^T(t)\mathbf{\upsilon }(t)\}\end{array}$$

(45)

where $T_s=0.04s$ is the sampling time.

In this study, an initial covariance matrix $P_1(0)$, spectral density matrix $Q$ of system noise and measurement noise covariance matrix $R$ are given as follows:

$$\begin{array}{l}{P}_1(0)=diag\{{(0.1 \mathrm{m}/\mathrm{s})}^2,{(0.1 \mathrm{m}/\mathrm{s})}^2,{(1°)}^2,{(1°)}^2,{(1°)}^2,\\ {({10}^{-4}g)}^2,{({10}^{-4}g)}^2,{(0.2°/h)}^2,{(0.2°/h)}^2,{(0.2°/h)}^2,{(0.2°/h)}^2,{(0.2°/h)}^2,{(0.2°/h)}^2\}\\ Q=diag\{{(2\times {10}^{-5} \mathrm{m}/\mathrm{s})}^2,{(2\times {10}^{-5} \mathrm{m}/\mathrm{s})}^2,{({0.12}^{″})}^2,{({0.12}^{″})}^2,{({0.12}^{″})}^2,\\ {(0.0001°/h)}^2,{(0.0001°/h)}^2,{(0.0001°/h)}^2,{(0.004°/h)}^2,{(0.004°/h)}^2,{(0.004°/h)}^2\}\\ R=diag\{{(0.01 \mathrm{m}/\mathrm{s})}^2,{(0.01 \mathrm{m}/\mathrm{s})}^2\}\end{array}$$

(46)

When using the continuous rotation method based on the extended observation model, measurement noise covariance matrix $R$ is expressed as follows:

$$\begin{array}{l}R=\{\begin{array}{c}\begin{array}{c}diag\{{(0.01 \mathrm{m}/\mathrm{s})}^2,{(0.01 \mathrm{m}/\mathrm{s})}^2\}\begin{array}{c}when2(k-1)\pi <{\omega }_{0}t<2k\pi \end{array}\end{array}\\ \begin{array}{cc}diag\{{(0.01 \mathrm{m}/\mathrm{s})}^2,{(0.01 \mathrm{m}/\mathrm{s})}^2,{(4^{″})}^2,{(4^{″})}^2,{(4^{″})}^2\}& when{\omega }_{0}t=2k\pi \end{array}\end{array}\\ k=1,2,3,\dots \end{array}$$

(47)

The rotation rate of the turntable is ${\omega }_0=10°/s$. The number of iterations performed for calculating $P$ using Equation (44) is 15,000 which is equivalent to 600 s. For two-position alignment, the IMU changes position at 300 s. Since the heading error ${\varphi }_D$ is the most crucial error state in the north-finding system, we focus on the RMS value of ${\varphi }_D$.

Figure 3 shows the RMS values of the heading error in the north-finding process. Obviously, the new continuous rotation alignment with the extended observation is more efficient than the existing north-finding algorithms.

In order to analyze the gyroscope error compensation effect of the new continuous rotation alignment approach, we use Allan variance technique to compare the compensated data with the uncompensated data of the equivalent east gyroscope, which determines the north-finding accuracy in a north-finding system.

The uncompensated equivalent east gyroscope data, denoted as ${\tilde{\omega }}_{ibE}^n$ is the measurement of the equivalent east gyroscope in the n frame, when the turntable is not rotating, that is:

$${\tilde{\omega }}_{ibE}^n=[\begin{array}{ccc}0& 1& 0\end{array}]C_{b_0}^n{C}_b^{b_0}{\tilde{\mathbf{\omega }}}_{ib}^b\begin{array}{c}{C}_b^{b_0}=I\end{array}$$

(48)

The compensated equivalent east gyroscope data, denoted as ${\widehat{\omega }}_{ibE}^n$ is the measurement of the equivalent east gyroscope in the n frame, when the turntable is rotating. The compensated data is obtained after the Kalman filter has converged. The drift error of the gyroscope has been estimated and compensated by the Kalman filter. That is:

$${\widehat{\omega }}_{ibE}^n=[\begin{array}{ccc}0& 1& 0\end{array}]C_{b_0}^n{C}_b^{b_0}\left({\tilde{\mathbf{\omega }}}_{ib}^b-\left[\begin{array}{c}{\widehat{\epsilon }}_{bx}+{\widehat{\epsilon }}_{rx}\\ {\widehat{\epsilon }}_{by}+{\widehat{\epsilon }}_{ry}\\ {\widehat{\epsilon }}_{bz}+{\widehat{\epsilon }}_{rz}\end{array}\right]-\left[\begin{array}{c}{\widehat{\epsilon }}_{mx}\\ {\widehat{\epsilon }}_{my}\\ {\widehat{\epsilon }}_{mz}\end{array}\right]\right)$$

(49)

The sampling data are collected over 3 h as shown in Figure 2, and the sampling frequency is 10 Hz. As shown in Figure 4, after compensation, the bias instability of the equivalent east gyroscope is almost eliminated, but the ARW remains as before. It should be noticed that RRW is almost eliminated through the continuous rotation modulation.

The experiment demonstrated that the RRW and Markov noise could be compensated by continuous rotation alignment, but ARW remained unchanged. The theoretical proofs are shown in Appendix B.

5. Experimental Results

The experimental platform is shown in Figure 5.

Considering the installation error, it is difficult to determine the absolute north. The previous north-finding experimental result was used as a reference to evaluate the performance of the approaches. The assumed azimuth was the mean value of 15 experimental results in two weeks north-finding tests. In this study, the experimental north-finding system stayed on a fixed azimuth. For each north-finding algorithm, the north-finding process was repeated five times.

Since the errors of the gyroscopes and accelerometers are unobservable in the fixed-position alignment, which may cause the divergence of the Kalman filter in the practice. We used 5-state Kalman filter for the fixed-position alignment. From Equations (15)–(23), the model can be expressed as follows:

$$\begin{array}{l}x(t)={\left[\delta v_N,\delta v_E,{\varphi }_N,{\varphi }_E,{\varphi }_D\right]}^T\\ \dot{x}(t)=F_1(t)x(t)+G(t)w(t)\\ w(t)={\left[\begin{array}{ccccc}{w}_{accx}& w_{accy}& w_{ax}& w_{ay}& w_{az}\end{array}\right]}^T\\ z(t)=H_{1}x(t)+\mathbf{\upsilon }(t)H_1=\left[\begin{array}{cc}{I}_{2\times 2}& 0_{2\times 3}\end{array}\right]\end{array}$$

(50)

The coarse alignment method using the gravity in the initial frame as a reference was employed in the experiments [17].

As shown in Figure 6a–d, the azimuth errors converged with time, the experimental results are coincident with the simulation analysis as shown in Figure 3 in which the new continuous rotation alignment with extended observation is the most efficient algorithm for a Coriolis vibration gyroscope based north-finding system.

In order to further compare the performances of the north-finding methods, we changed the azimuth of the north-finding system to 6 different directions as shown in Equation (51):

$${\phi }_1=-20.337°,{\phi }_k={\phi }_1+(k-1)60°k=2,\dots 6$$

(51)

For each azimuth, the north-finding process was repeated for 5 times with the 4 different north-finding algorithms. Then, the RMS of heading errors for each of these algorithms was calculated. As shown in Figure 7, the new approach (continuous rotation alignment with the extended observation model) is the best one, the north-finding accuracy is 0.1° (1σ).

6. Conclusions

As analyzed in this paper, it is the gyroscope random drift errors that make it a challenge for a cost effective gyroscope based north-finding systems to be achieved. Since it is the equivalent east gyroscope that determines the north-finding accuracy, Allan variance analysis of the equivalent east gyroscope before and after error compensation provides an efficient technique for the evaluation of the gyroscope error estimation.

Comparisons of the Kalman filter convergence rapidity and north-finding accuracy have been made to evaluate the north-finding algorithms. Compared with the other traditional approaches, the new continuous rotation alignment approach based on the extended observation model can improve the north-finding accuracy and convergence rapidity effectively. The experiments have shown that a heading accuracy of 0.1° (1σ) can be achieved in 10 min at 28.22° north latitude using a HRG IMU with gyro bias instability of $0.1°/h$, compared with 0.6° (1σ) north-finding accuracy for the two-position alignment and 1° (1σ) for the fixed-position alignment.

In fact, ARW, RRW and Markov noise are the main error source of many gyroscopes (e.g., fiber optic gyroscopes [18]). The new continuous rotation IMU alignment algorithm is not only applicable to the Coriolis vibration gyros (a kind of cost effective HRGs in this paper), but is also suitable for many other gyroscopes with similar stochastic error models.