An Improved Online Fast Self-Calibration Method for Dual-Axis RINS Based on Backtracking Scheme

Abstract

In the field of high accuracy dual-axis rotational inertial navigation system (RINS), the calibration accuracy of the gyroscopes and accelerometers is of great importance. Although rotation modulation can suppress the navigation error caused by scale factor error and bias error in a static condition, it cannot suppress the scale factor errors thoroughly during the maneuvering process of the vehicle due to the two degrees of rotation freedom. The self-calibration method has been studied by many researchers. However, traditional calibration methods need several hours to converge, which is unable to meet the demand for quick response to positioning and orientation. To solve the above problems, we do the following work in this study: (1) we propose a 39-dimensional online calibration Kalman filtering (KF) model to estimate all calibration parameters; (2) Error relationship between calibration parameters error and navigation error are derived; (3) A backtracking filtering scheme is proposed to shorten the calibration process. Experimental results indicate that the proposed method can shorten the calibration process and improve the calibration accuracy simultaneously.

1. Introduction

Rotational inertial navigation systems (RINS) have reached high accuracy of navigation in recent years [1,2]. A RINS consists of three gyroscopes, three accelerometers, and one rotational table. the inertial measurement unit (IMU) is installed inside the two-axis rotation table [3]. The accuracy of RINS highly depends on the accuracy of the accelerometers and gyroscopes, hence, the calibration parameters of the IMU are very important [4]. The RINS makes the online self-calibration of IMU a real, self-calibration is an effective way to improve the navigation performance of RINS [5].

Many researchers have studied the INS calibration method [6,7,8,9]. Ren et al. proposed a multi-position self-calibration method for dual-axis RINS [10]. In [11], Zhang et al. analyzed the multi-position calibration method for IMU. In [12], an eight-position self-calibration method for a dual-axis RINS. Syed et al. proposed a multi-position calibration method for MEMS inertial navigation systems [13]. Jiang et al. in [14] have proposed a 36-dimensional KF for systematic calibration. Wen et al. have added gyro-accelerometer asynchronous time as a state variable in self-calibration KF [15]. Jing et al. in [1] have analyzed the data of the gyroscope to improve the accuracy of the calibration. In [16], Song et al. have added inner lever-arm parameters for RINS, which denotes the distance between the sensitive center of the accelerometer and the center of rotation of the IMU. Based on previous research, the authors have derived a complete IMU calibration error model. However, a higher filtering dimension means slower computation and slower convergence. Nowadays, the demand for the rapid response of navigation equipment is higher, and it will be very inappropriate to self-calibration for a few hours before utilization. It should be noted that in high accuracy INS applications (the accuracy of the gyro is better than 0.01 ${}^{\circ }$/h, generally ring laser gyro or fiber optic gyro), the systematic calibration process can last for 3 h.

Hence, many researchers have utilized the backtracking navigation method to shorten the filtering time, this method has been applied in initial alignment widely to shorten the alignment time. Yan et al. proposed a reverse navigation algorithm for INS gyrocompass in-movement alignment [17]. Li has applied the backtracking navigation method in underwater vehicle applications [18]. In [19,20], the authors have proposed a backtracking integration for fast attitude determination-based initial alignment. Li et al. utilized a backtracking navigation scheme for autonomous underwater vehicles [20]. In the application field of RINS, Song et al. proposed a rapid initial alignment scheme for dual-axis RINS [21]. Based on the previous studies, the backtracking navigation scheme has been proved that can effectively shorten the alignment time in the INS initial alignment process. However, there are few studies have applied the backtracking scheme in RINS calibration method.

The backtracking scheme can be applied to the RINS calibration method in theory. However, we need to propose a rotation scheme that can excite all the calibration parameter errors during the rotation process. Hence, a detailed observability analysis is needed. In [15], the authors have analyzed the observability of the 19-position rotation scheme. Tang et al. has utilized piecewise constant systems (PWCS) and singular value decomposition (SVD) method to provide the observability analysis of KF [22]. Cai et al. provided an observability analysis for IMU in a three-axis rotation table [23]. In [24], the authors have studied the observability of strapdown INS alignment. However, the decoupling relationship between the calibration parameters and navigation velocity error has not been studied in previous studies. Therefore, there is a lack of a specific theoretical guide for the rotation path design to excite all the calibration parameter errors.

To solve the above problems, an improved online fast calibration method for dual-axis RINS based on a backtracking scheme is proposed in this paper. The self-calibration method has been studied by many researchers. However, traditional calibration methods need several hours to converge, which is unable to meet the demand for quick response to positioning and orientation. Kalman filtering has both accuracy and real-time performance, but the convergence is slow, so our method can ensure short-term convergence and the real-time performance of the algorithm. The processer we utilize in dual-axis RINS is DSP6748, it has lower computing power because it is a low-power chip. So KF is used in the system. To solve the above problems, we do the following work in this study: (1) we propose a 39-dimensional online calibration Kalman filtering (KF) model to estimate all calibration parameters; (2) Error relationship between calibration parameter errors and navigation error is derived, which can be a theoretical guide for the design of calibration rotation path; (3) A backtracking filtering scheme is proposed to shorten the calibration process. Experimental results indicates that the proposed method can shorten the calibration process to 1 h and improve the calibration accuracy (within 0.5 ppm) simultaneously.

The remainder of this paper is organized as follows: Section 2 gives reference definitions of the proposed calibration method. In Section 3, we derive a 39-dimensional filtering model and analyze the error relationship between calibration parameters error and navigation errors. In Section 4, the online self-calibration method based on the backtracking scheme, is provided. A rotation test is carried out in Section 5 to verify the effectiveness of the proposed method. In Section 6, we give the conclusions.

2. Reference Definitions

In the filtering process of the calibration method, the reference definitions are shown as follows:

The earth-centered inertial frame (i-frame): o locates on the center of the earth, $o{x}_i$ points to the vernal equinox, $o{z}_i$ is the earth’s axis of self-rotation, the inerital information are measured in this reference.

Earth-centred frame (e-frame): o locates on center of the earth, $o{x}_e$ points to the central meridian, $o{z}_e$ is the earth’s axis of self rotation.

IMU frame (I-frame): o locates on the rotation center of the IMU, $o{x}_I$, $o{y}_I$ and $o{z}_I$ points to the IMU’s right, forward, and upward, respectively.

vehicle frame (b-frame): o locates on the rotation center of the vehicle, $o{x}_I$, $o{y}_I$ and $o{z}_I$ points to the vehicle’s right, forward, and upward, respectively.

navigation frame (n-frame): The navigation coordinate system is the selected coordinate system when calculating the navigation parameters. In this paper, the east-north-up (E-N-U) geographic coordinate system is used as the navigation coordinate system.

gyro sensitive frame (g frame): Non-orthogonal reference frame aligned with gyro-sensitive axes.

accelerometer-sensitive frame (a frame): Non-orthogonal reference frame aligned with accelerometer-sensitive axes.

3. Sefl-Calibration Filtering Method Design

3.1. IMU Error Model

The main measurement components of the IMU are the three-axis gyro component and the three-axis accelerometer component. The input and output models of the gyro and accelerometer can be expressed as [15]:

$$\left\{\begin{array}{c}{\mathit{K}}_G^{-1}{\mathit{N}}^g={\mathit{T}}_g^{g^{\prime }}{\mathit{\omega }}_{ig}^g+{\mathit{b}}^g+{\mathit{v}}^g\hfill \\ {\mathit{K}}_A^{-1}{\mathit{N}}^a={\mathit{T}}_a^{a^{\prime }}{\mathit{f}}_{sf}^a+{\mathit{f}}_{A2}^I+{\mathit{b}}^a+{\mathit{v}}^a\hfill \end{array}\right.$$

(1)

where, ${\mathit{K}}_G$ is the scale factor matrix (consists of scale factors of three gyros) of gyroscope, ${\mathit{K}}_A$ is the scale factor matrix of accelerometers, ${\mathit{v}}^g$ is the noise vector of gyros, ${\mathit{v}}^a$ is the noise vector of accelerometers, ${\mathit{\omega }}_{ig}^g$ is the angular velocity sensed by gyros. g-frame and a-frame denotes the gyroscope and accelerometer’s sensitive axes, g-frame and a-frame are non-orthogonal references. ${\mathit{N}}^g={\left[\begin{array}{ccc}{N}_x^g& N_y^g& N_z^g\end{array}\right]}^{\mathrm{T}}$ is the primitive output of the gyroscopes, ${\mathit{N}}^a={\left[\begin{array}{ccc}{N}_x^a& N_y^a& N_z^a\end{array}\right]}^{\mathrm{T}}$ is the primitive output of the accelerometers. ${\mathit{\omega }}_{ig}^g$ and ${\mathit{f}}_{sf}^a$ are angular velocity and specific force vector. ${\mathit{K}}_G=diag\left(K_{gx},K_{gy},K_{gz}\right)$ and ${\mathit{K}}_A=diag\left(K_{ax},K_{ay},K_{az}\right)$ are scale factor matrix of the gyroscopes and accelerometers respectivley. ${\mathit{b}}^g={\left[\begin{array}{ccc}{K}_{G0x}& K_{G0y}& K_{G0z}\end{array}\right]}^{\mathrm{T}}$ and ${\mathit{b}}^a={\left[\begin{array}{ccc}{K}_{A0x}& K_{A0y}& K_{A0z}\end{array}\right]}^{\mathrm{T}}$ are the bias vectors of gyroscopes and accelerometers. Accelerometer second-order nonlinear coefficient specific force sensitive term (which can be estimated in high accuracy accelerometers like quartz accelerometer) is described as: ${\mathit{f}}_{A2}^I={\mathit{M}}_f{\mathit{K}}_{A2}$, where ${\mathit{M}}_f=diag\left({\left(f_{sfx}^a\right)}^2,{\left(f_{sfy}^a\right)}^2,{\left(f_{sfz}^a\right)}^2\right)$ is the matrix composed of the square term of the theoretical specific force sensitive to the triaxial accelerometer. ${\mathit{K}}_{A2}={\left[\begin{array}{ccc}{K}_{A2x}& K_{A2y}& K_{A2z}\end{array}\right]}^{\mathrm{T}}$ is the second-order nonlinear coefficient vector of the accelerometer, it is a part of scale fator.

Write the orthogonal gyro and accelerometer coordinate systems after coordinate frame $g^{\prime }$ and $a^{\prime }$ are compensated by the axis misalignment angle matrix ${\mathit{T}}_g^{g^{\prime }}$ and ${\mathit{T}}_a^{a^{\prime }}$, respectively. The Equation (1) can be rewritten as:

$$\left\{\begin{array}{c}{\mathit{\omega }}_{ig}^g={\left({\mathit{T}}_g^{g^{\prime }}\right)}^{-1}\left({\mathit{K}}_G^{-1}{\mathit{N}}^g-{\mathit{b}}^g-{\mathit{v}}^g\right)\hfill \\ {\mathit{f}}_{sf}^a={\left({\mathit{T}}_a^{a^{\prime }}\right)}^{-1}\left({\mathit{K}}_A^{-1}{\mathit{N}}^a-{\mathit{f}}_{A2}^I-{\mathit{b}}^a-{\mathit{v}}^a\right)\hfill \end{array}\right.$$

(2)

As shown in Figure 1, since the measurement information of the gyroscope and accelerometer needs to be projected into the I system in the inertia calculation, this paper defines that the measurement center of the gyroscope coincides with the sensitive center of the IMU coordinate system, and the $o{x}_I$ axis of the IMU coordinate system coincides with the sensitive axis $o{x}_g$ of the gyroscope.

Hence, the installation matrix can be written as:

$${\mathit{T}}_g^I\approx \left[\begin{array}{ccc}1& 0& 0\\ -{\gamma }_{yz}^g& 1& 0\\ {\gamma }_{zy}^g& -{\gamma }_{zx}^g& 1\end{array}\right]$$

(3)

$${\mathit{T}}_a^I={\mathit{C}}_p^I{\mathit{T}}_a^p\approx \left[\begin{array}{ccc}1& -{\eta }_z& {\eta }_y\\ {\eta }_z-{\vartheta }_{yz}& 1& -{\eta }_x\\ -{\eta }_y+{\vartheta }_{zy}& {\eta }_x-{\vartheta }_{zx}& 1\end{array}\right]=\left[\begin{array}{ccc}1& {\gamma }_{xz}^a& -{\gamma }_{xy}^a\\ -{\gamma }_{yz}^a& 1& {\gamma }_{yx}^a\\ {\gamma }_{zy}^a& -{\gamma }_{zx}^a& 1\end{array}\right]$$

(4)

where $\gamma$ denotes the installation angles, ${\mathit{T}}_a^p\approx \left[\begin{array}{ccc}1& 0& 0\\ -{\vartheta }_{yz}& 1& 0\\ {\vartheta }_{zy}& -{\vartheta }_{zx}& 1\end{array}\right]$, ${\vartheta }_{ij}$ is the sensitive non-orthogonal angles from p-frame to a-frame, p-frame is the othogonal reference (The p-frame is a transition orthogonal coordinate system), ${\mathit{C}}_p^I\approx \left[\begin{array}{ccc}1& -{\eta }_z& {\eta }_y\\ {\eta }_z& 1& -{\eta }_x\\ -{\eta }_y& {\eta }_x& 1\end{array}\right]$, ${\eta }_i$ is the non-orthogonal angles from I-frame to p-frame.

Furthermore, the effect of the inner lever arm ${\mathit{r}}^I$ is considered. The inner lever arm refers to the distance vector from the IMU sensitive center (IMU rotation center) to the three-axis accelerometer. When the carrier moves angularly, the inner lever arm will cause the tangential acceleration and centripetal acceleration that the accelerometer is sensitive to. Since the installation error angle of the accelerometer is a small angle, only the inner lever arm in the direction of the sensitive axis of the accelerometer is considered, and the rotational angular acceleration is not considered, the specific force-sensitive term ${\mathit{f}}_r^I$ of the inner lever arm of the accelerometer can be expressed as [15]:

$${\mathit{f}}_r^I={\mathit{M}}_{\omega }{\mathit{r}}^I$$

(5)

where, ${\mathit{M}}_{\omega }=diag\left[{\left({\omega }_{iIy}^I\right)}^2+{\left({\omega }_{iIz}^I\right)}^2,{\left({\omega }_{iIx}^I\right)}^2+{\left({\omega }_{iIz}^I\right)}^2,{\left({\omega }_{iIx}^I\right)}^2+{\left({\omega }_{iIy}^I\right)}^2\right]$, ${\mathit{r}}^I={\left[\begin{array}{ccc}{r}_x^I& r_y^I& r_z^I\end{array}\right]}^{\mathrm{T}}$.

Ignoring the noise term in Equation (1), the calibration models of the gyro and accelerometer components in the IMU coordinate system can be expressed as:

$$\left\{\begin{array}{c}{\mathit{\omega }}_{iI}^I={\mathit{K}}_G^{-1}{\mathit{N}}^g-{\mathit{\omega }}_0^I\hfill \\ {\mathit{f}}_{sf}^I={\mathit{K}}_A^{-1}{\mathit{N}}^a-{\mathit{f}}_{A2}^I-{\mathit{f}}_0^I\hfill \end{array}\right.$$

(6)

The error equation is the differential equation of the above equations, based on Equations (2) and (6). ${\mathit{\omega }}_0^I$ is the gyro bias vector, ${\mathit{f}}_0^I$ is the accelerometer bias vector.

3.2. 39-Dimensional Kalman Filtering Model

The error parameters of the system-level filter self-calibration include two categories. The first category is the various error parameters involved in the IMU measurement error model; the second category is the observation error parameters represented by the outer lever arm error. The outer lever arm error refers to the sensitivity of the IMU. The distance error vector from the center to the center of the turntable will stimulate the observation velocity error and the observation position error when the IMU rotates [15]. The outter lever arm are written as $\delta {\mathit{l}}^I$, The relationship between the corresponding true value and the measured value can be expressed as: ${\mathit{l}}^I={\tilde{\mathit{l}}}^I-\delta {\mathit{l}}^I$. Bsed on the previous analysis, combined with the inertial navigation error equation, attitude misalignment angle error ${\mathit{\varphi }}^n={\left[\begin{array}{ccc}{\varphi }_E^n& {\varphi }_N^n& {\varphi }_U^n\end{array}\right]}^{\mathrm{T}}$, the velocity error $\delta {\mathit{v}}_{}^n={\left[\begin{array}{ccc}\delta v_E^n& \delta v_N^n& \delta v_U^n\end{array}\right]}^{\mathrm{T}}$ and the position error $\delta \mathit{p}={\left[\begin{array}{ccc}\delta L& \delta \lambda & \delta h\end{array}\right]}^{\mathrm{T}}$ into the state vector into consideration, the state variable of the KF can be written as:

$${\mathit{X}}_C={\left[\begin{array}{cccccc}{\left({\varphi }^n\right)}^{\mathrm{T}}& {\left(\delta {\mathit{v}}_{}^n\right)}^{\mathrm{T}}& {\left(\delta \mathit{p}\right)}^{\mathrm{T}}& {\left({\mathit{X}}_g\right)}^{\mathrm{T}}& {\left({\mathit{X}}_a\right)}^{\mathrm{T}}& {\left(\delta {\mathit{l}}^I\right)}^{\mathrm{T}}\end{array}\right]}^{\mathrm{T}}$$

(7)

where, ${\mathit{X}}_g={\left[\begin{array}{cccccc}\delta k_{11}^g& \delta k_{21}^g& \delta k_{31}^g& \delta k_{22}^g& \delta k_{32}^g& \delta k_{33}^g\end{array}\begin{array}{c}\end{array}\begin{array}{ccc}{\epsilon }_x& {\epsilon }_y& {\epsilon }_z\end{array}\right]}^T$.

$\begin{array}{c}{\mathit{X}}_a=[\begin{array}{cccccc}\delta k_{11}^a& \delta k_{21}^a& \delta k_{31}^a& \delta k_{12}^a& \delta k_{22}^a& \delta k_{32}^a\end{array}\hfill \\ \begin{array}{ccccccccc}\delta k_{13}^a& \delta k_{23}^a& \delta k_{33}^a& {\nabla }_x& {\nabla }_y& {\nabla }_z& \delta K{T}_x^a& \delta K{T}_y^a& \delta K{T}_z^a\end{array}{]}^T\hfill \end{array}$. $\delta k_{ij}^g$ is the element of ${\mathit{K}}_G^{-1}$, ${\epsilon }_i$ (${\mathit{\omega }}_0^I$) is the gyro bias, $\delta k_{ij}^a$ is the element of ${\mathit{K}}_A^{-1}$, and ${\nabla }_x$ (${\mathit{f}}_0^I$) is the accelerometer bias.

The state transform funciton can be written as:

$${\dot{\mathit{X}}}_C={\mathit{F}}_C{\mathit{X}}_C+{\mathit{G}}_C{\mathit{W}}_C$$

(8)

where, ${\mathit{F}}_C$ represents the state transition matrix, ${\mathit{G}}_C$ represents the system noise driving matrix, ${\mathit{W}}_C$ denotes the system nosie matrix, ${\mathbf{\epsilon }}_w^I={\left[\begin{array}{ccc}{\epsilon }_{wx}^I& {\epsilon }_{wy}^I& {\epsilon }_{wz}^I\end{array}\right]}^{\mathrm{T}}$ and ${\nabla }_w^I={\left[\begin{array}{ccc}{\nabla }_{wx}^I& {\nabla }_{wy}^I& {\nabla }_{wz}^I\end{array}\right]}^{\mathrm{T}}$ are the random noise matrix of gyroscopes and accelerometers. ${\mathit{F}}_C$ amd ${\mathit{G}}_C$ can be written as:

$${\mathit{F}}_C=\left[\begin{array}{c}\begin{array}{cccccc}-\left({\mathit{\omega }}_{in}^n\times \right)& {\mathit{F}}_{12}& {\mathit{F}}_{13}& {\mathit{F}}_{14}& 0_{3\times 18}& 0_{3\times 3}\\ \left({\mathit{f}}_{sf}^n\times \right)& {\mathit{F}}_{22}& {\mathit{F}}_{23}& 0_{3\times 9}& {\mathit{F}}_{25}& 0_{3\times 3}\\ 0_{3\times 3}& {\mathit{F}}_{32}& {\mathit{F}}_{33}& 0_{3\times 9}& 0_{3\times 18}& 0_{3\times 3}\end{array}\\ 0_{30\times 39}\end{array}\right]$$

(9)

$${\mathit{G}}_C=\left[\begin{array}{c}\begin{array}{cc}-{\mathit{C}}_I^n& 0_{3\times 3}\\ 0_{3\times 3}& {\mathit{C}}_I^n\end{array}\\ 0_{33\times 6}\end{array}\right]$$

(10)

where, ${\mathit{C}}_I^n$ is the transform matrix from I-frame to n-frame. ${\mathit{\omega }}_{ie}^{}$ represents the earth self-rotatoin angular velocity. ${\mathit{\omega }}_{in}^n$ is the angular velocity of n-frame to i-frame. L is the latitude. $R_M$ denotes the radius of curvature of the meridian circle, $R_N$ represents the the radius of curvature of the unitary circle, h is the height of the IMU. The elements of the ${\mathit{F}}_C$ are shown as: ${\mathit{F}}_{12}={\mathit{M}}_2$, ${\mathit{F}}_{13}={\mathit{M}}_1+{\mathit{M}}_3$, ${\mathit{F}}_{14}=-{\mathit{C}}_I^n{\mathit{M}}_g$, ${\mathit{F}}_{22}=-\left[\left({\mathit{v}}^n\times \right){\mathit{F}}_{12}+\left(2{\mathit{\omega }}_{ie}^n+{\mathit{\omega }}_{en}^n\right)\times \right]$, ${\mathit{M}}_1=\left[\begin{array}{ccc}0& 0& 0\\ -{\mathit{\omega }}_{ie}sinL& 0& 0\\ {\mathit{\omega }}_{ie}cosL& 0& 0\end{array}\right]$,

$${\mathit{M}}_2=\left[\begin{array}{ccc}0& -\frac{1}{R_M+h}& 0\\ \frac{1}{R_N+h}& 0& 0\\ \frac{tanL}{R_N+h}& 0& 0\end{array}\right],{\mathit{M}}_3=\left[\begin{array}{ccc}0& 0& \frac{v_N}{{\left(R_N+h\right)}^2}\\ 0& 0& -\frac{v_E}{{\left(R_N+h\right)}^2}\\ \frac{v_E{sec}^{2}L}{R_N+h}& 0& -\frac{v_{E}tanL}{{\left(R_N+h\right)}^2}\end{array}\right],$$

$${\mathit{F}}_{23}=\left({\mathit{v}}^n\times \right)\left(2{\mathit{M}}_1+{\mathit{M}}_3\right),$$

$${\mathit{F}}_{25}={\mathit{C}}_I^n{\mathit{M}}_a,{\mathit{F}}_{32}=\left[\begin{array}{ccc}0& \frac{1}{R_M+h}& 0\\ \frac{1}{(R_N+h)cosL}& 0& 0\\ 0& 0& 1\end{array}\right],{\mathit{F}}_{33}=\left[\begin{array}{ccc}0& 0& -\frac{v_N}{{\left(R_N+h\right)}^2}\\ \frac{v_{E}secLtanL}{R_N+h}& 0& -\frac{v_{E}secL}{{\left(R_N+h\right)}^2}\\ 0& 0& 0\end{array}\right].$$

$${\mathit{M}}_g=\left[\begin{array}{cccc}{\tilde{\omega }}_{ibx}^I{\mathit{I}}_{3\times 3}& \left[\begin{array}{c}{0}_{1\times 2}\\ {\tilde{\omega }}_{iIy}^I{\mathit{I}}_{2\times 2}\end{array}\right]& \left[\begin{array}{c}{0}_{2\times 1}\\ {\tilde{\omega }}_{iIz}^I\end{array}\right]& {\mathit{I}}_{3\times 3}\end{array}\right],$$

$${\mathit{M}}_a=\left[\begin{array}{cccccc}{\tilde{f}}_{sfx}^I{\mathit{I}}_{3\times 3}& {\tilde{f}}_{sfy}^I{\mathit{I}}_{3\times 3}& {\tilde{f}}_{sfz}^I{\mathit{I}}_{3\times 3}& {\tilde{\mathit{M}}}_{f2}& {\tilde{\mathit{M}}}_{\omega }& {\mathit{I}}_{3\times 3}\end{array}\right],\mathrm{where},$$

$${\tilde{\mathit{M}}}_f=diag\left[{\left({\tilde{f}}_{sfx}^I\right)}^2,{\left({\tilde{f}}_{sfy}^I\right)}^2,{\left({\tilde{f}}_{sfz}^I\right)}^2\right],$$

$${\tilde{\mathit{M}}}_{\omega }=diag\left[{\left({\tilde{\omega }}_{iIy}^I\right)}^2+{\left({\tilde{\omega }}_{iIz}^I\right)}^2,{\left({\tilde{\omega }}_{iIx}^I\right)}^2+{\left({\tilde{\omega }}_{iIz}^I\right)}^2,{\left({\tilde{\omega }}_{iIx}^I\right)}^2+{\left({\tilde{\omega }}_{iIy}^I\right)}^2\right].$$

To establish the measurement equation of the dual-axis RINS system-level self-calibration filtering method, the measurement information should be determined first. The velocity and position measurement equation is expressed as follows:

$$\left\{\begin{array}{c}{\tilde{\mathit{v}}}_{mea}^n={\tilde{\mathit{v}}}_{}^n+{\tilde{\mathit{C}}}_I^n\left[({\tilde{\mathit{\omega }}}_{iI}^I-{\tilde{\mathit{C}}}_n^I{\tilde{\mathit{\omega }}}_{ie}^n)\times {\tilde{\mathit{l}}}^I\right]\hfill \\ {\tilde{\mathit{p}}}_{mea}^n=\tilde{\mathit{p}}+{\tilde{\mathit{F}}}_{32}{\tilde{\mathit{C}}}_I^n{\tilde{\mathit{l}}}^I\hfill \end{array}\right.$$

(11)

where, ${\tilde{\mathit{C}}}_I^n=\left[{\mathit{I}}_{3\times 3}-\left({\varphi }^n\times \right)\right]{\mathit{C}}_I^n$ represents the calculated attitude matrix with attitude misalignment.

The measurement transform equation can be written as:

$$\begin{array}{c}{\mathit{Z}}_C^n=\left[\begin{array}{c}{\mathit{Z}}_{Cv}^n\hfill \\ {\mathit{Z}}_{Cp}^n\hfill \end{array}\right]=\left[\begin{array}{c}{\tilde{\mathit{v}}}_{mea}^n-{\mathit{v}}_{obv}^n\\ {\tilde{\mathit{p}}}_{mea}-{\mathit{p}}_{obv}\end{array}\right]\\ ={\mathit{H}}_C{\mathit{X}}_C+{\mathit{V}}_C\end{array}$$

(12)

where ${\mathit{Z}}_C^n$ denotes the measurement vector, ${\mathit{H}}_C$ is the measurement matrix,

$${\mathit{H}}_C=\left[\begin{array}{cccccc}{\mathit{H}}_{11}& {\mathit{I}}_{3\times 3}& {\mathit{H}}_{13}& {\mathit{H}}_{14}& 0_{3\times 18}& {\mathit{H}}_{16}\\ {\mathit{H}}_{21}& 0_{3\times 3}& {\mathit{H}}_{23}& 0_{3\times 9}& 0_{3\times 18}& {\mathit{H}}_{26}\end{array}\right]$$

(13)

where ${\mathit{H}}_{21}={\mathit{F}}_{32}\left[\left({\mathit{C}}_I^n{\mathit{l}}^I\right)\times \right]$, ${\mathit{H}}_{23}=\left[\begin{array}{ccc}1& 0& -\frac{{\left({\mathit{C}}_I^n{\mathit{l}}^I\right)}_y}{{\left(R_M+h\right)}^2}\\ \frac{sinL{\left({\mathit{C}}_I^n{\mathit{l}}^I\right)}_x}{(R_N+h){cos}^{2}L}& 1& -\frac{{\left({\mathit{C}}_I^n{\mathit{l}}^I\right)}_x}{{(R_N+h)}^{2}cosL}\\ 0& 0& 1\end{array}\right]$, ${\mathit{H}}_{26}={\mathit{F}}_{32}{\mathit{C}}_I^n$.

So far, the derivation of the 39-dimensional filtering equation is completed.

3.3. Self-Calibration Parameter Error Excitation and Coupling Analysis

To carry out the coupling anaysis of each calibration parameters, the error propagation law of each self-calibration parameters are derived first. Assuming that the initial state b-frame and n-frame are coincide (${\mathit{C}}_b^n\left(t\right)={\mathit{I}}_3$), considering the case where the y-axis of the IMU continues to rotate in the north direction at time period $\left[0,T\right]$ with an angular velocity ${\omega }_r({\omega }_r\gg {\omega }_{ie})$. The measured angular velocity can be expressed as: ${\tilde{\mathit{\omega }}}_{iI}^I\approx {\left[\begin{array}{ccc}0& {\omega }_r& 0\end{array}\right]}^{\mathrm{T}}+\delta {\mathit{\omega }}_{iI}^I$. The measured specific force can be expressed as ${\tilde{\mathit{f}}}_{sf}^I={\left[\begin{array}{ccc}-gsin\left({\omega }_{r}T\right)& 0& gcos\left({\omega }_{r}T\right)\end{array}\right]}^{\mathrm{T}}+\delta {\mathit{f}}_{sf}^I$, the attitude transform matrix can be described as ${\mathit{C}}_I^n\left(t\right)={\mathit{C}}_b^n\left(t\right){\mathit{C}}_I^b\left(t\right)$. Ignoring higher-order episilon (more than two-order), The measurement error due to the $\delta {\mathit{K}}_G$ and bias ${\epsilon }^I$ can be described as:

$$\begin{array}{c}\delta {\dot{\mathit{Z}}}_{Cv}^n\left(\delta {\mathit{K}}_G,{\epsilon }^I\right)\approx \left({\mathit{g}}^n\times \right){\int }_0^T{\mathit{C}}_b^n\left(t\right){\mathit{C}}_I^b\left(t\right){\mathit{M}}_g\left(t\right)dt·{\mathit{X}}_g\hfill \\ \approx \left({\mathit{g}}^n\times \right){\int }_0^T\left(\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{ccc}cos\left({\omega }_{r}T\right)& 0& sin\left({\omega }_{r}T\right)\\ 0& 1& 0\\ -sin\left({\omega }_{r}T\right)& 0& cos\left({\omega }_{r}T\right)\end{array}\right]\right)dt·\left[\begin{array}{c}{\epsilon }_x^I\\ {\omega }_r\delta k_{yy}^g+{\epsilon }_y^I\\ {\omega }_r\delta k_{zy}^g+{\epsilon }_z^I\end{array}\right]\hfill \\ =\left[\begin{array}{ccc}0& g& 0\\ -g& 0& 0\\ 0& 0& 0\end{array}\right]\left[\begin{array}{ccc}\frac{sin\left({\omega }_{r}T\right)}{{\omega }_r}& 0& -\frac{cos\left({\omega }_{r}T\right)-1}{{\omega }_r}\\ 0& T& 0\\ \frac{cos\left({\omega }_{r}T\right)-1}{{\omega }_r}& 0& \frac{sin\left({\omega }_{r}T\right)}{{\omega }_r}\end{array}\right]\left[\begin{array}{c}{\epsilon }_x^I\\ {\omega }_r\delta k_{yy}^g+{\epsilon }_y^I\\ {\omega }_r\delta k_{zy}^g+{\epsilon }_z^I\end{array}\right]\hfill \\ =\left[\begin{array}{c}g{\omega }_{r}T\delta k_{yy}^g+gT{\epsilon }_y^I\\ gcos\left({\omega }_{r}T\right)\left(\delta k_{zy}^g+\frac{{\epsilon }_z^I}{{\omega }_r}\right)-\frac{gsin\left({\omega }_{r}T\right){\epsilon }_x^I}{{\omega }_r}-g\delta k_{zy}^g-\frac{g{\epsilon }_z^I}{{\omega }_r}\\ 0\end{array}\right]\hfill \end{array}$$

(14)

Ignoring higher-order epsilon, The measurement error due to the $\delta {\mathit{K}}_A$ and ${\nabla }^I$ can be described as:

$$\begin{array}{c}\delta {\dot{\mathit{Z}}}_{Cv}^n\left(\delta {\mathit{K}}_A,{\nabla }^I\right)={\mathit{C}}_I^n{\mathit{M}}_a(:,1:12){\mathit{X}}_a(:,1:12)\hfill \\ \approx \left[\begin{array}{ccc}cos\left({\omega }_{r}T\right)& 0& sin\left({\omega }_{r}T\right)\\ 0& 1& 0\\ -sin\left({\omega }_{r}T\right)& 0& cos\left({\omega }_{r}T\right)\end{array}\right]\left[\begin{array}{c}-gsin\left({\omega }_{r}T\right)\delta k_{xx}^a+gcos\left({\omega }_{r}T\right)\delta k_{xz}^a+{\nabla }_x^b\\ -gsin\left({\omega }_{r}T\right)\delta k_{yx}^a+gcos\left({\omega }_{r}T\right)\delta k_{yz}^a+{\nabla }_y^b\\ -gsin\left({\omega }_{r}T\right)\delta k_{zx}^a+gcos\left({\omega }_{r}T\right)\delta k_{zz}^a+{\nabla }_z^b\end{array}\right]\hfill \\ =\left[\begin{array}{c}\begin{array}{c}\frac{gsin\left(2{\omega }_{r}T\right)}{2}\left(\delta k_{zz}^a-\delta k_{xx}^a\right)+\frac{g}{2}\left(\delta k_{xz}^a-\delta k_{zx}^a\right)+\frac{gcos\left(2{\omega }_{r}T\right)}{2}\left(\delta k_{xz}^a+\delta k_{zx}^a\right)+\hfill \\ {\nabla }_x^{b}cos\left({\omega }_{r}T\right)+{\nabla }_z^{b}sin\left({\omega }_{r}T\right)\hfill \end{array}\\ -gsin\left({\omega }_{r}T\right)\delta k_{yx}^a+gcos\left({\omega }_{r}T\right)\delta k_{yz}^a+{\nabla }_y^b\\ \begin{array}{c}\frac{g}{2}\left(\delta k_{xx}^a+\delta k_{zz}^a\right)+\frac{gcos\left(2{\omega }_{r}T\right)}{2}\left(\delta k_{zz}^a-\delta k_{xx}^a\right)-\frac{gsin\left(2{\omega }_{r}T\right)}{2}\left(\delta k_{xz}^a+\delta k_{zx}^a\right)-\hfill \\ {\nabla }_x^{b}sin\left({\omega }_{r}T\right)+{\nabla }_z^{b}cos\left({\omega }_{r}T\right)\hfill \end{array}\end{array}\right]\hfill \end{array}$$

(15)

Similarly, the measurement error due to the $\delta {\mathit{K}}_{A2}$ can be described as:

$$\begin{array}{c}\delta {\dot{\mathit{Z}}}_{Cv}^n\left(\delta {\mathit{K}}_{A2}\right)={\mathit{C}}_I^n{\tilde{\mathit{M}}}_{f2}{\mathit{X}}_a(:,13:15)\hfill \\ \approx \left[\begin{array}{ccc}cos\left({\omega }_{r}T\right)& 0& sin\left({\omega }_{r}T\right)\\ 0& 1& 0\\ -sin\left({\omega }_{r}T\right)& 0& cos\left({\omega }_{r}T\right)\end{array}\right]\left[\begin{array}{c}{g}^2{sin}^2\left({\omega }_{r}T\right)\delta K_{Ax2}\\ 0\\ g^2{cos}^2\left({\omega }_{r}T\right)\delta K_{Az2}\end{array}\right]\hfill \\ =\left[\begin{array}{c}{g}^{2}sin\left(2{\omega }_{r}t\right)\left(sin\left({\omega }_{r}t\right)\delta K_{Ax2}+cos\left({\omega }_{r}t\right)\delta K_{Az2}\right)/2\\ 0\\ -g^2{sin}^3\left({\omega }_{r}t\right)\delta K_{Ax2}+g^2{cos}^3\left({\omega }_{r}t\right)\delta K_{Az2}\end{array}\right]\hfill \end{array}$$

(16)

The measurement error due to the $\delta {\mathit{r}}^I$ can be written as:

$$\begin{array}{c}\delta {\dot{\mathit{Z}}}_v^n\left(\delta {\mathit{r}}^I\right)={\mathit{C}}_I^n{\tilde{\mathit{M}}}_{\omega }{\mathit{X}}_a(:,16:18)\hfill \\ \approx \left[\begin{array}{ccc}cos\left({\omega }_{r}T\right)& 0& sin\left({\omega }_{r}T\right)\\ 0& 1& 0\\ -sin\left({\omega }_{r}T\right)& 0& cos\left({\omega }_{r}T\right)\end{array}\right]\left[\begin{array}{c}{\left({\omega }_r\right)}^2\delta r_x^I\\ 0\\ {\left({\omega }_r\right)}^2\delta r_z^I\end{array}\right]\hfill \\ =\left[\begin{array}{c}{\left({\omega }_r\right)}^2\left(cos\left({\omega }_{r}T\right)\delta r_x^I+sin\left({\omega }_{r}T\right)\delta r_z^I\right)\\ 0\\ {\left({\omega }_r\right)}^2\left(-sin\left({\omega }_{r}T\right)\delta r_x^I+cos\left({\omega }_{r}T\right)\delta r_z^I\right)\end{array}\right]\hfill \end{array}$$

(17)

Based on Equations (14)–(17), The error propagation form of self-calibration parameters under transposition excitation are summarized in

Table 1. Error propagation form of self-calibration parameters under transposition excitation.

Parameters	$\mathit{\delta }{\dot{\mathit{Z}}}_{\mathit{C}{\mathit{v}}_{\mathit{E}}}^{\mathit{n}}$	$\mathit{\delta }{\dot{\mathit{Z}}}_{\mathit{C}{\mathit{v}}_{\mathit{N}}}^{\mathit{n}}$	$\mathit{\delta }{\dot{\mathit{Z}}}_{\mathit{C}{\mathit{v}}_{\mathit{U}}}^{\mathit{n}}$
$\delta k_{zy}^g$	–	$gcos\left({\omega }_{r}T\right)$	–
$\delta k_{yz}^a$	–	$gcos\left({\omega }_{r}T\right)$	–
$\delta k_{yx}^a$	–	$-gsin\left({\omega }_{r}T\right)$	–
$\delta k_{xz}^a$	$\frac{gcos\left(2{\omega }_{r}T\right)}{2}$	–	$-\frac{gsin\left(2{\omega }_{r}T\right)}{2}$
$\delta k_{zx}^a$	$\frac{gcos\left(2{\omega }_{r}T\right)}{2}$	–	$-\frac{gsin\left(2{\omega }_{r}T\right)}{2}$
$\delta k_{yy}^g$	$g{\omega }_{r}T$	–	–
$\delta k_{xx}^a$	$-\frac{gsin\left(2{\omega }_{r}T\right)}{2}$	–	$-\frac{gcos\left(2{\omega }_{r}T\right)}{2}$
${\nabla }_x^I$	$cos\left({\omega }_{r}T\right)$	–	$-sin\left({\omega }_{r}T\right)$
${\nabla }_z^I$	$sin\left({\omega }_{r}T\right)$	–	$cos\left({\omega }_{r}T\right)$
${\epsilon }_x^I$	–	$-\frac{gsin\left({\omega }_{r}T\right)}{{\omega }_r}$	–
${\epsilon }_z^I$	–	$\frac{gcos\left({\omega }_{r}T\right)-g}{{\omega }_r}$	–
$\delta K_{Ax2}^{}$	$g^{2}sin\left(2{\omega }_{r}T\right)sin\left({\omega }_{r}T\right)/2$	–	$-g^2{sin}^3\left({\omega }_{r}T\right)$
$\delta K_{Az2}^{}$	$g^{2}sin\left(2{\omega }_{r}T\right)cos\left({\omega }_{r}T\right)/2$	–	$g^2{cos}^3\left({\omega }_{r}T\right)$
$\delta r_x^I$	${\left({\omega }_r\right)}^{2}cos\left({\omega }_{r}T\right)$	–	${\left({\omega }_r\right)}^{2}sin\left({\omega }_{r}T\right)$
$\delta r_z^I$	$-{\left({\omega }_r\right)}^{2}sin\left({\omega }_{r}T\right)$	–	${\left({\omega }_r\right)}^{2}cos\left({\omega }_{r}T\right)$

:

Based on Table 1, the error propagation forms of the $\delta k_{zy}^g$ and $\delta k_{yz}^g$, $\delta k_{xz}^a$ and $\delta k_{zx}^a$ are the same respectively. When the IMU rotates along the y-axis, the above two sets of installation error parameters will always have a coupling relationship. Similarly, the IMU rotates along the x-axis, The error propagation forms of the three groups of parameters $\delta k_{xy}^a$ and $\delta k_{yx}^g$, $\delta k_{xz}^a$ and $\delta k_{zx}^g$, $\delta k_{zy}^a$ and $\delta k_{yz}^a$ are respectively the same. When the IMU rotates along the z-axis, the error propagation forms of the $\delta k_{yx}^a$ and $\delta k_{xy}^a$ are the same.

Although there are 1–3 sets of coupling relationships in the installation errors of the gyro and the accelerometer during the horizontal rotation of the three sensitive axes of the IMU, such coupling relationships can be decoupled in turn through a specific indexing method.

Except for the coupling term in the installation error, there is no error term in the same propagation form during the horizontal rotation of the three sensitive axes of the IMU. According to the error propagation form of each system-level self-calibration parameter, the longest error propagation period is $2\pi /{\omega }_r$ (derived by Equation (17)).

Based on the previous analysis, the rotation path of the self-calibration can be designed as

Table 2. Rotation path of self-calibration process.

Time	Rotation Axis (Inner (I) (z-Axis of IMU)/Outer (O) (x-Axis of IMU))	Rotation Angle along I/O Axis	Attitude after Rotation (XYZ) (East-North-Upward)
0 s	-	-	ENU
180 s	O	+90${}^{\circ }$	EUS
270 s	O	+180${}^{\circ }$	EDN
360 s	O	+180${}^{\circ }$	EUS
450 s	I	+90${}^{\circ }$	UWS
540 s	I	+180${}^{\circ }$	DES
630 s	I	+180${}^{\circ }$	UWS
720 s	O	+90${}^{\circ }$	SWD
810 s	O	+180${}^{\circ }$	NWU
900 s	O	+180${}^{\circ }$	SWD
990 s	O	+90${}^{\circ }$	DWN
1080 s	O	+90${}^{\circ }$	NWU
1170 s	O	+90${}^{\circ }$	UWS
1260 s	I	+90${}^{\circ }$	WDS
1350 s	I	+90${}^{\circ }$	DES
1440 s	I	+90${}^{\circ }$	EUS
1530 s	O	+90${}^{\circ }$	ESD
1620 s	O	+90${}^{\circ }$	EDN
1710 s	O	+90${}^{\circ }$	ENU

[15]:

The angular rotation rate is set as 20 ${}^{\circ }$/s (the angular rate is not unique, usually 5 ${}^{\circ }$/s to 20 ${}^{\circ }$/s are commonly used), after finishing the last rotation, stay static until 1800 s to end the calibration process.

In [15] the error parameters errors of IMU scale factors and biases are proved observable. Therefore, we need to analyze the observability of inner lever arms and accelerometer second-order factors in this study. Using the PWCS and SVD methods, the observability using the designed excitation path is shown as follows in

Table 3. The observability degree of each state.

State Variable	Singular Value	State Variable	Singular Value
$\delta K_{Ax2}^{}$	232.6811	$\delta K_{Ay2}^{}$	197.5225
$\delta K_{Az2}^{}$	51.2133	$\delta r_x^I$	43.6781
$\delta r_x^I$	30.7752	$\delta r_x^I$	22.1765

.

4. Self-Calibration Process Based on Backtracking Scheme

The reverse navigation algorithm is a time-reversed algorithm based on the forward navigation algorithm. It is usually used to prolong the filtering time so that the filter can complete the convergence in a relatively short time. In the process of reverse navigation, the speed of the solution needs to be reversed. The reverse navigation algorithm can be rewritten as the following:

$$\begin{array}{c}{\mathit{C}}_{bk-1}^n={\mathit{C}}_{bk}^n({\mathit{I}}_3+T_s{\tilde{\mathbf{\Omega }}}_{nbk}^b)\hfill \\ -{\mathit{v}}_{k-1}^n=-{\mathit{v}}_k^n+T_s{\mathit{a}}_{nbk-1}^n\hfill \end{array}$$

(18)

where, ${\mathbf{\Omega }}_{nbk}^b=({\mathit{\omega }}_{nbk}^b\times )$.

The detailed derivation of reverse navigation method has been fully explained in reference [17], hence, reverse navigation process can be summarized as:

$$\begin{array}{c}{\overleftarrow{\mathit{C}}}_{bp}^n={\overleftarrow{\mathit{C}}}_{bp-1}^n(\mathit{I}+T_s{\overleftarrow{\mathbf{\Omega }}}_{nbp}^b)\hfill \\ {\overleftarrow{\mathit{v}}}_p^n={\overleftarrow{\mathit{v}}}_{p-1}^n+T_s{\overleftarrow{\mathit{a}}}_{p-1,p}\hfill \\ {\overleftarrow{L}}_p={\overleftarrow{L}}_{p-1}+T_s{\overleftarrow{v}}_{Np-1}^n/(R_M+{\overleftarrow{h}}_{p-1})\hfill \\ {\overleftarrow{\lambda }}_p={\overleftarrow{\lambda }}_{p-1}+T_s{v}_{Ep-1}^{n}sec{\overleftarrow{L}}_{p-1}/(R_N+{\overleftarrow{h}}_{p-1})\hfill \\ {\overleftarrow{h}}_p={\overleftarrow{h}}_{p-1}+T_s{\overleftarrow{v}}_{Up-1}^n\hfill \end{array}$$

(19)

where, $p=m-k+1$.

As shown in Equation (19), the algorithm form of the reverse navigation algorithm is similar to that of the forward navigation. It only needs to invert the speed, the angular rate of the earth’s rotation, and the stored gyro output.

In the process of reverse navigation, due to the minus output of the gyroscope, the output of the gyroscope can be written as:

$${\mathit{F}}_{14}=-{\mathit{C}}_b^n\left[\begin{array}{cccc}{N}_x^g{\mathit{I}}_{3\times 3}& \left[\begin{array}{c}{0}_{1\times 2}\\ N_y^g{\mathit{I}}_{2\times 2}\end{array}\right]& \left[\begin{array}{c}{0}_{2\times 1}\\ N_z^g\end{array}\right]& -{\mathit{I}}_{3\times 3}\end{array}\right]$$

(20)

Other filter variables are consistent with the forward self-calibration filter model.

To complete the self-calibration process of the dual-axis RINS in a short time, it is necessary to propose a method to enable the self-calibration filter to be able to convergence is accomplished rapidly. Therefore, this paper proposes a fast self-calibration algorithm architecture based on a backtracking scheme, the architecture is shown in Figure 2 as follows:

The 39D KF method is to decouple the error through the specified indexing sequence, to decouple the relationship between the navigation error and the IMU error, and use Kalman filtering to estimate. To ensure the linearity of the KF, we usually set initial calibration parameters at the beginning of the filtering process. The self-calibration filtering process-based backtracking scheme proposed in this paper makes full use of IMU information. Under the condition that the initial calibration parameters are set (inaccurate) (the coarse alignment is necessary, providing a coarse attitude can guarantee the linearity of the Kaman filter), the RINS can obtain an initial attitude after coarse alignment. The raw data of the IMU is stored at the same time. Then, utilize the attitude value after coarse alignment as the initial attitude value, utilize the stored gyroscope and accelerometer data to perform the reverse self-calibration process, and continue to store the IMU data. After the reverse filtering process end, the forward filtering process is performed from the initial moment, because the calculation rate of the self-calibration filtering process is faster than the rate of storing the data, and the final forward filtering process will end at time $t_{calib}$, which can catch up with the data stored in real-time and finally enter into the navigation process.

Given the above explanations, the self-calibration based on the backtracking scheme can be summarized as follows:

• The information from the accelerometer and gyroscope are stored in the memory in real-time during the time period from the start (time 0) to the end (time $t_{coarse}$) of the coarse alignment stage of the dual-axis RINS, and the coarse alignment process ends at time $t_{coarse}$. We obtain attitude information with an acceptable error.
• The process of reverse self-calibration filtering starts at the end of the coarse alignment of the dual-axis RINS (time $t_{coarse}$), and the process of reverse self-calibration filtering towards the start of the coarse alignment stage (time 0), using the stored gyro data, the angular velocity of the earth’s rotation, and the reverse self-calibration The velocity at the initial moment of filtering needs to be negated. At the same time, in the process of reverse self-calibration filtering, the information of the gyroscope and accelerometer are still stored in the memory of the navigation computer.
• After the reverse self-calibration filtering process ends, the program executes to the start time of the coarse alignment stage (time 0), and from this time onwards, the forward self-calibration filtering process is performed without modifying all the parameters of the Kalman filter used in the previous stage, the velocity at time 0 calculated using the reverse self-calibration filter needs to be reversed. At the same time, the stored gyro information and the Earth’s rotation angular rate also need to be returned to the normal state from the previous inversion state (no need to invert). Since the calculation rate of the navigation computer is not fast enough to be ignored, the forward self-calibration filtering process needs to catch up with the stored time $t_{calib}$ until.
• The above three steps are the implementation process of the fast self-calibration algorithm based on backtracking navigation. At present, most navigation computers are equipped with large-capacity storage elements such as DDR2, SD card, etc., so that the online fast self-calibration algorithm of dual-axis RINS can be realized.

It should be noted that if the filter still does not fully converge after one round of retrospective filtering calibration, the next round of retrospective filtering and calibration can be performed. The reverse filtering process starts from the previous $t_{calib}$, and the subsequent process is the same as before.

5. Experimental Results and Analysis

To verify the effectiveness of the proposed self-calibration method based on the backtracking scheme, we conduct a static test to evaluate the accuracy of the calibration parameters. The dual-axis we utilize in this study is shown in Figure 3. The dual-axis RINS realize the dual-axis turntable with 3D-IMU.

The structure of the dual-axis RINS is shown as follows Figure 4, and the IMU is installed inside the dual-axis turntable to realize the self-contained dual-axis RINS.

The gyroscope we use in the RINS is a laser ring gyro (RLG), the accuracy of the RLG is 0.003 ${}^{\circ }$/h (100 s, $1\sigma$), with a 1 ppm ($1\sigma$) of scale factor repeatability. The quartz accelerometers have an accuracy of 20 $\mathsf{\mu }$g ($1\sigma$). The sampling frequency of the IMU is 1000 Hz. To show the calibration process more clearly, the raw plot of the IMU is shown in Figure 5.

The RINS is fixed in a marble, the algorithm is implemented on the digital signal processer (DSP) chip. We use the method in reference to [10] as a comparison. The self-calibration process lasts 30 mins (the backtracking process also lasts 30 min, with the same data). We use a high-accuracy three-axis turntable to calibrate the IMU parameters as reference [25], this method requires a high-accuracy turntable, the IMU needs to be removed from the dual-axis RINS, as the accuracy of the dual-axis turntable is not high enough (especially horizontal accuracy). The traditional method is described in [15]. The estimation curves of the IMU parameters are shown in Figure 6, Figure 7, Figure 8 and Figure 9, the dotted line part is the reverse filtering process, and the solid line part is the forward filtering process. We can find out all parameters are converged after the backtracking scheme ended.

To verify the estimation accuracy of the parameters, the estimated parameters can be summarized as

Table 4. Estimation results of different methods.

Estimated Parameters	Proposed Method	Traditional Method	Reference Values
${\epsilon }_x$	−0.01345 ${}^{\circ }$/h	−0.03953 ${}^{\circ }$/h	−0.01955 ${}^{\circ }$/h
${\epsilon }_y$	0.012587 ${}^{\circ }$/h	0.04478 ${}^{\circ }$/h	0.01685 ${}^{\circ }$/h
${\epsilon }_z$	0.04521 ${}^{\circ }$/h	0.1023 ${}^{\circ }$/h	0.04002 ${}^{\circ }$/h
${\nabla }_x$	412.23 $\mathsf{\mu }$g	415.22 $\mathsf{\mu }$g	412.75 $\mathsf{\mu }$g
${\nabla }_y$	−812.36 $\mathsf{\mu }$g	−813.56 $\mathsf{\mu }$g	−813.74 $\mathsf{\mu }$g
${\nabla }_z$	694.25 $\mathsf{\mu }$g	691.57 $\mathsf{\mu }$g	695.12 $\mathsf{\mu }$g
$\delta k_{xx}^g$	206,263.25 ${}^{\circ }$/h/pulse	206,269.44 ${}^{\circ }$/h/pulse	206,263.25 ${}^{\circ }$/h/pulse
$\delta k_{yy}^g$	206,269.34 ${}^{\circ }$/h/pulse	206,260.98 ${}^{\circ }$/h/pulse	206,268.87 ${}^{\circ }$/h/pulse
$\delta k_{zz}^g$	206,267.22 ${}^{\circ }$/h/pulse	206,268.74 ${}^{\circ }$/h/pulse	206,267.84 ${}^{\circ }$/h/pulse
$\delta k_{xx}^a$	98,021.66 m/${\mathrm{s}}^2$/pulse	98,020.97 m/${\mathrm{s}}^2$/pulse	98,021.39 m/${\mathrm{s}}^2$/pulse
$\delta k_{yy}^a$	98,015.23 m/${\mathrm{s}}^2$/pulse	98,019.37 m/${\mathrm{s}}^2$/pulse	98,015.54 m/${\mathrm{s}}^2$/pulse
$\delta k_{zz}^a$	98,036.94 m/${\mathrm{s}}^2$/pulse	98,031.52 m/${\mathrm{s}}^2$/pulse	98,036.56 m/${\mathrm{s}}^2$/pulse
$\delta k_{yx}^g$	3.547′	4.125′	3.368′
$\delta k_{zx}^g$	−2.365′	−3.122′	−2.674′
$\delta k_{zy}^g$	11.245′	11.544′	10.941′
$\delta k_{xy}^a$	9.124′	9.426′	9.221′
$\delta k_{xz}^a$	7.586′	7.138′	7.225′
$\delta k_{yx}^a$	1.747′	1.529′	1.596′
$\delta k_{yz}^a$	−5.618′	−5.221′	−5.625′
$\delta k_{zx}^a$	3.027′	3.291′	3.171′
$\delta k_{zy}^a$	6.107′	5.822′	6.128′
$\delta K_{Ax2}$	−5.754 ppm	−7.225 ppm	−5.551 ppm
$\delta K_{Ay2}$	23.485 ppm	29.569 ppm	21.993 ppm
$\delta K_{Az2}$	30.241 ppm	27.226 ppm	30.453 ppm

.

As shown in Table 4, the estimation accuracy of the proposed method is better than the traditional method, especially the calibration parameters of gyros. The errors of gyro biases estimated by the traditional method are 0.02 ${}^{\circ }$/h to 0.03 ${}^{\circ }$/h, using the proposed method, the errors are only within 0.005 ${}^{\circ }$/h. The errors of the gyro scale factors estimated by the traditional method are more than 4 ppm.

Compared with the traditional self-calibration process, in a relatively short time, the method proposed in this paper achieves high estimation accuracy, which can verify the effectiveness of the proposed fast self-calibration method based on the backtracking scheme.

6. Conclusions

In the field of the rotational inertial navigation system, many fields (land and aircraft) have put forward “three-self” performance requirements, namely self-checking, self-alignment, and self-calibration. Although the traditional strapdown INS can realize self-checking and self-alignment, it cannot realize self-calibration. In response to such problems, dual-axis came into research, and the introduction of the dual-axis indexing mechanism fundamentally solved the problem of self-calibration without disassembly of IMU from the dual-axis turntable. Therefore, dual-axis RINS has received extensive attention in the field of land navigation. The accuracy of the calibration parameters can determine the navigation accuracy of RINS. The traditional self-calibration method needs several hours to converge. To shorten the self-calibration time by more than 50 percent, we propose a 39-dimensional online calibration Kalman filtering (KF) model to estimate all calibration parameters. Error relationship between calibration parameter errors and navigation errors are derived, which can be a theoretical guide for the design of calibration rotation path. A backtracking filtering scheme is proposed to shorten the calibration process. Experimental results indicate that the proposed method can shorten the calibration process and improve the calibration accuracy simultaneously compared with the traditional self-calibration method.