An improved self-calibration method with consideration of inner lever-arm effects for a dual-axis rotational inertial navigation system

As shown in figure 1, the dual-axis RINS mainly consists of the gimbal structure (inner gimbal and outer gimbal) and an IMU. On both sides of the gimbals, a brushless DC motor drives the rotation and a photoelectric encoder is used for collecting the rotation angle. Three fiber optics gyroscopes and three quartz flexible accelerometers are arranged approximately orthogonally on the IMU, together with relevant circuit boards. The device specifications of the dual-axis RINS are listed in table 1.

Zoom In Zoom Out Reset image size

Table 1. Device specifications of the dual-axis RINS.

Characteristics	Parameters
Gyro drift stability	0.02° h−1
Gyro stochastic error	0.001 $^\circ\,\mathrm h^{-\frac{1}{2}}$
Gyro scale factor stability	30 ppm
Accelerometer bias stability	50 ${{\mu }}g$
Accelerometer scale factor stability	50 ppm
Photoelectric encoder accuracy	1''
Sampling frequency	200 Hz

Navigation frame ($n,O-{x_n}{y_n}{z_n}$): the n-frame sets a reference coordinate to operate the navigation algorithm. $O{x_n}$ and $O{y_n}$ are horizontal axes pointing east and north, $O{z_n}$ is a vertical axis pointing skyward.

Body frame ($b,O-{x_b}{y_b}{z_b}$): the b-frame describes the shell of the dual-axis RINS. $O{x_b}$ (the pitch axis) is pointing in the rightward direction, and $\theta$ is the pitch angle. $O{y_b}$ (the roll axis) is a horizontal axis pointing in the forward direction, and $\gamma$ represents the roll angle. $O{z_n}$ (the azimuth axis) constitutes a right-handed orthogonal frame together with $O{x_b}$ and $O{y_b}$, and Ψ is the azimuth angle. $C_n^b$ is a direction cosine matrix describing the relationship between the n-frame and b-frame, denoted by

$$\begin{align} C_n^b &= R(\gamma )R(\theta )R(\psi ) \nonumber\\[4pt] &= \left[\!\!\! \begin{array}{*{20}{c}} {\cos \gamma }&0&{\sin \gamma }\\ 0&1&0\\ { - \sin \gamma }&0&{\cos \gamma } \end{array}\!\!\!\right]\left[\!\!\!\begin{array}{*{20}{c}} 1&0&0\\ 0&{\cos \theta }&{ - \sin \theta }\\ 0&{\sin \theta }&{\cos \theta } \end{array}\!\!\!\right]\left[\!\!\!\begin{array}{*{20}{c}} {\cos \psi }& { - \sin \psi }&0\\ {\sin \psi }& {\cos \psi }&0\\ 0&0&1 \end{array} \!\!\!\right]. \end{align} \tag{ 1 }$$

Platform frame ($p,O-x_py_pz_p$): the p-frame describes the coordinates of the IMU. ${x_g},{y_g},{z_g}$ represent the sensitive axes of the three gyros, and ${x_a},{y_a},{z_a}$ represent the sensitive axes of the three accelerometers. The definition of the orthogonal p-frame can start from the sensitive axis of any inertial component. As shown in figure 2, we first define the sensitive axis of the x-accelerometer as ${x_p}$ of the p-frame, and the plane perpendicular to ${x_p}$ is the $O{y_p}{z_p}$ plane. The orientation of the ${y_p}$ axis is the projection of the y-accelerometer sensitive axis in the $O{y_p}{z_p}$ plane. ${z_p}$ constitutes a right-handed orthogonal frame together with ${x_p}$ and ${y_p}$. By this definition, the orthogonal relationship of the three accelerometers can be described by three installation errors, as shown in figure 2(a). After that, the spatial relationship between each gyro sensitive axis and the corresponding IMU axis needs to be described by two installation errors. Therefore, a total of six installation errors are needed to establish the relationship between the gyro sensitive frame and the p-frame, as shown in figure 2(b).

Zoom In Zoom Out Reset image size

The installation errors are always considered to be small angles. Taking the relationship between ${y_p}$ and ${y_a}$ as an example, the relationship between ${y_p}$ and ${y_a}$ can be described as

$$\begin{align} {y_a} &= {y_p}\cos (\alpha _{ay}^Z) - {x_p}\sin (\alpha _{ay}^Z) \nonumber\\[4pt] &= {y_p}\cos (\alpha _{ay}^Z) - {x_a}\sin (\alpha _{ay}^Z) \nonumber \\[4pt] & \approx {y_p} - {x_a} \cdot \alpha _{ay}^Z \end{align} \tag{ 2 }$$

so

$$\begin{equation}{y_p} = {y_a} + {x_a} \cdot \alpha _{ay}^Z.\end{equation} \tag{ 3 }$$

Similarly, we can find the relationship between $z_{\mathrm p}$ and ${z_a}$, and then the matrix required to transform the accelerometer sensitive axes to the p-frame can be described as

$$\begin{equation}C_a^p{\rm{ = }}\left[ {\begin{array}{*{20}{c}} 1&0&0\\ {\alpha _{ay}^Z} & 1 & 0\\[4pt] { - \delta _{az}^Y} & {\delta _{az}^X} & 1 \end{array}} \right].\end{equation} \tag{ 4 }$$

Similarly, the matrix required to transform the gyro sensitive axes to the p-frame can be described as

$$\begin{equation}C_g^p = \left[ {\begin{array}{*{20}{c}} 1 & { - \alpha _{gx}^Z} & {\beta _{gx}^Y} \\[4pt] {\alpha _{gy}^Z} & 1 & { - \beta _{gy}^X} \\[4pt] { - \delta _{gz}^Y} & {\delta _{gz}^X} & 1 \end{array}} \right]\end{equation} \tag{ 5 }$$

where, $\alpha _{ay}^Z,\delta _{az}^X,\delta _{az}^Y,\alpha _{gx}^Z,\alpha _{gy}^Z,\beta _{gx}^Y,\beta _{gy}^X,\delta _{gz}^X,\delta _{gz}^Y$ represent the installation errors in the IMU. The naming principle is as follows: taking $\beta _{gx}^Y$ as an example, the subscript indicates it is a misalignment error of the x-gyro, and the Y superscript means the angle is generated along the ${y_p}$ axis.

The error modeling for the dual-axis RINS can be mainly divided into two categories: one is the traditional IMU sensor errors and the other is the inner lever-arm parameters. The traditional IMU sensor errors include gyro drifts, accelerometer biases, gyro and accelerometer scale factor errors and installation errors. The relationship between the angular rates of the gyro sensitive axes and the p-frame can be described as

$$\begin{align} \omega_{ip}^{p}&= (C_{g}^{p}+K_{g})\omega^{g}+\varepsilon^{p}\nonumber\\[4pt] &= \left[ {\begin{array}{*{20}{c}} {1 + {K_{gx}}}&{ - \alpha _{gx}^Z}&{\beta _{gx}^Y}\\[4pt] {\alpha _{gy}^Z}&{1 + {K_{gy}}}&{ - \beta _{gy}^X}\\[4pt] { - \delta _{gz}^Y}&{\delta _{gz}^X}&{1 + {K_{gz}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {\omega _x^g}\\[4pt] {\omega _y^g}\\[4pt] {\omega _z^g} \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} {\varepsilon _x^p}\\[4pt] {\varepsilon _y^p}\\[4pt] {\varepsilon _z^p} \end{array}} \right]\end{align} \tag{ 6 }$$

where $\omega _{ip}^p$ is the angular rate vector in the p-frame, ${{K_g} = diag\left[{\begin{array}{@{}lll@{}} {{K_{gx}}} & {{K_{gy}}} & {{K_{gz}}} \end{array}}\right]}$ is the gyro scale factor error matrix and ${\omega ^g} = {\left[\!\!\!\!\!{\begin{array}{*{20}{c}} {\omega _x^g} & {\omega _y^g} & {\omega _z^g} \end{array}}\!\!\!\!\!\right]^T}$ is the angular rate vector in the gyro sensitive axes. ${\varepsilon ^p} = {\left[\!\!\!\!\! {\begin{array}{*{20}{c}} {\varepsilon _x^p} & {\varepsilon _y^p} & {\varepsilon _z^p} \end{array}} \!\!\!\!\!\right]^T}$ is the gyro constant drift vector in the p-frame. Similarly, the relationship between the accelerations of the accelerometer sensitive axes and the p-frame could be described as

$$\begin{align} f_{ip}^{\,\,\,p} &= (C_a^p + {K_a}){f^{\,\,g}} + \Delta _{}^p \nonumber \\[4pt] &= \left[ {\begin{array}{*{20}{c}} {1 + {K_{ax}}} & 0 & 0 \\ [2pt] {\alpha _{ay}^Z} & {1 + {K_{ay}}} & 0 \\[4pt] { - \delta _{az}^Y} & {\delta _{az}^X} & {1 + {K_{az}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {f_x^{\,\,\,g}} \\[4pt] {f_y^{\,\,\,g}} \\[4pt] {f_z^{\,\,\,g}} \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} {\Delta _x^p} \\[4pt] {\Delta _y^p} \\[4pt] {\Delta _z^p} \end{array}} \right] \end{align} \tag{ 7 }$$

where $f_{ip}^{\,\,p}$ is the acceleration vector in the p-frame, ${{K_a} = diag\left[\!\!\!\!\!{\begin{array}{*{20}{c}} {{K_{ax}}} & {{K_{ay}}} & {{K_{az}}}\end{array}}\!\!\!\!\!\right]}$ is the accelerometer scale factor error matrix and ${f^{\,\,g}} = {\left[\!\! {\begin{array}{*{20}{c}} {f_x^{\,\,g}} & {f_y^{\,\,g}} & {f_z^{\,\,g}} \end{array}} \!\!\!\!\!\right]^T}$ is the acceleration vector in the accelerometer sensitive axes. ${{\Delta ^p} = {\left[{\begin{array}{@{}lll@{}} {\Delta _x^p} & {\Delta _y^p} & {\Delta _z^p} \end{array}}\right]^T}}$ is the accelerometer constant bias vector in the p-frame. From equations (6) and (7), the error model of traditional IMU sensor errors could be extended as follows:

$$\begin{equation}\Delta \omega _{ip}^{\,p} = \left[ {\begin{array}{*{20}{c}} {\Delta \omega _{ipx}^p} \\[4pt] {\Delta \omega _{ipy}^p} \\[4pt] {\Delta \omega _{ipz}^p} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{K_{gx}}\omega _x^g - \alpha _{gx}^Z\omega _y^g + \beta _{gx}^Y\omega _z^g + \varepsilon _x^p} \\[4pt] {{K_{gy}}\omega _y^g + \alpha _{gy}^Z\omega _x^g - \beta _{gy}^X\omega _z^g + \varepsilon _y^p} \\[4pt] {{K_{gz}}\omega _z^g - \delta _{gz}^Y\omega _x^g + \delta _{gz}^X\omega _y^g + \varepsilon _z^p} \end{array}} \right]\end{equation} \tag{ 8 }$$

$$\begin{equation}\Delta f_{ip}^{\,\,p}=\left[\begin{array}{c}\Delta f_{ipx}^{\,\,p}\\[4pt] \Delta f_{ipy}^{\,\,p}\\[4pt] \Delta f_{ipz}^{\,\,p}\end{array}\right]=\left[\begin{array}{c}K_{ax\,\,}f_x^{\,\,g}+\Delta_x^p\\[4pt] K_{ay\,\,}f_y^{\,\,g}+\alpha_{ay\,\,}^Zf_x^{\,\,g}+\Delta_y^p\\[4pt] K_{az\,\,}f_z^{\,\,g}-\delta_{az\,\,}^Yf_x^g+\delta_{az\,\,}^Xf_y^g+\Delta_z^p\end{array}\right].\end{equation} \tag{ 9 }$$

In actual assembly, the accelerometer sensitive points cannot completely coincide with the IMU center. The spatial relationship between the three accelerometers and the IMU center is shown in figure 3.

Zoom In Zoom Out Reset image size

In order to describe the spatial relationship between the accelerometer measurement center and the IMU center, three-dimensional coordinates for each accelerometer are required. For example, $(r_{xx}^{\,p},r_{xy}^{\,p},r_{xz}^{\,p})$ shows the coordinates of the x-accelerometer in the p-frame, the first subscript x representing the 'x-accelerometer' and the second subscript representing the direction of the inner lever-arm parameters. Similarly, $(r_{yx}^{\,p},r_{yy}^{\,p},r_{yz}^{\,p})$ and $(r_{zx}^{\,p},r_{zy}^{\,p},r_{zz}^{\,p})$ represent the coordinates of the inner lever-arm parameters for the y- and z-accelerometers.

According to the rigid motion theorem [12], when the IMU rotates, harmful accelerations caused by the inner lever-arm effect could be extended as

$$\begin{align} \Delta f_{ip\_r}^{\,p} &= \dot \omega _{ip}^p \times {r^{\,p}} + \omega _{ip}^p \times (\omega _{ip}^p \times {r^{\,p}}) \nonumber \\[2pt] &= \left[ {\begin{array}{*{20}{c}} {( - \omega {{_{ipy}^p}^2} - \omega {{_{ipz}^p}^2})r_{xx}^{\,p} + ( - \dot \omega _{ipz}^p + \omega _{ipx}^p\omega _{ipy}^p)r_{xy}^{\,p} + (\dot \omega _{ipy}^p + \omega _{ipx}^p\omega _{ipz}^p)r_{xz}^p} \\[4pt] {(\dot \omega _{ipz}^p + \omega _{ipx}^p\omega _{ipy}^p)r_{yx}^{\,p} + ( - \omega {{_{ipz}^p}^2} - \omega {{_{ipx}^p}^2})r_{yy}^{\,p} + ( - \dot \omega _{ipx}^p + \omega _{ipy}^p\omega _{ipz}^p)r_{yz}^{\,p}} \\[4pt] {( - \dot \omega _{ipy}^p + \omega _{ipx}^p\omega _{ipz}^p)r_{zx}^{\,p} + (\dot \omega _{ipx}^p + \omega _{ipy}^p\omega _{ipz}^p)r_{zy}^{\,p} + ( - \omega {{_{ipx}^p}^2} - \omega {{_{ipy}^p}^2})r_{zz}^{\,p}} \end{array}} \right] \end{align} \tag{ 10 }$$

where $\Delta f_{ip\_r}^{\,\,p}$ is the harmful acceleration vector and $\dot \omega _{ip}^p$ represents the angular acceleration vector. In this equation, $\Delta f_{ip\_r}^{\,\,p}$ mainly includes two components: $\dot \omega _{ip}^p \times {r^p}$ is the tangential acceleration caused by $r_{xy}^{\,p},r_{xz}^p,r_{yx}^{\,p},r_{yz}^{\,p},r_{zx}^{\,p},r_{zy}^{\,p}$, and $\omega _{ip}^p \times (\omega _{ip}^p \times {r^p})$ is the centripetal acceleration caused by $r_{xx}^{\,p},r_{yy}^{\,p},r_{zz}^{\,p}$. It can be seen that $r_{xx}^{\,p},r_{yy}^{\,p},r_{zz}^{\,p}$ are relative to the angular rate, which can be activaed by any kind of rotation, while $r_{xy}^{\,p},r_{xz}^{\,p},r_{yx}^{\,p},r_{yz}^{\,p},r_{zx}^{\,p},r_{zy}^{\,p}$ are relative to angular accelerations and can only be activacted by variable angular rate rotation.

Based on the published research [15], the principles of rotation scheme design in self-calibration can be summarized under three aspects. Firstly, all parameters should be fully activated by adopting the rotation scheme, and the observable degrees of state variables should be as high as possible. Secondly, the rotation scheme should have simplicity for motor control. Thirdly, the calibration time should be as short as possible, because rapid self-calibration could improve the carrier maneuverability in military applications.

In the traditional self-calibration method, gyro scale factor errors and installation errors could be separated from gyro drifts when the rotation scheme contains changeable angular rate input, such as continuous bi-directional rotation. Similarly, accelerometer scale factor errors and installation errors could be separated from accelerometer biases when changeable acceleration input exists. The gravitational acceleration is an important incentive for static base calibration, so upward and downward movement is always involved in rotation schemes. In summary, bi-directional rotation with a uniform angular rate along the IMU with three directions pointing to the horizontal axis sequentially is considered to be the most effective rotation scheme [15]. Based on these principles, figure 4 provides a typical rotation scheme (Scheme A) for the dual-axis RINS, and the details are described as follows.

Zoom In Zoom Out Reset image size

Step 1: Inner gimbal locks at zero, while the outer gimbal rotates bi-directionally with a uniform angular rate of 6° s⁻¹ and repeats this process for five cycles. In this step, the IMU rotates along the ${y_p}$ axis in the horizontal plane.

Step 2: Inner gimbal rotates to 90° and then locks. After that, the outer gimbal rotates bi-directionally with a uniform angular rate of 6° s⁻¹ and repeats this process for five cycles. In this step, the IMU rotates along the ${x_p}$ axis in the horizontal plane.

Step 3: Inner gimbal rotates back to zero and then the outer gimbal rotates to 90°. After that, the inner gimbal rotates bi-directionally with a uniform angular rate of 6° s⁻¹ and repeats this process for five cycles. In this step, the IMU rotates along the ${z_p}$ axis in the horizontal plane.

By adopting Scheme A, all IMU traditional sensor errors could be activated in sequence. However, refering to equation (8), uniform rotational movement may lead to the inner lever-arm errors caused by $r_{xx}^{\,p},r_{yy}^{\,p},r_{zz}^{\,p}$, which will bring non-ignorable errors to velocity observations. If inner lever-arm parameters are not established in state variables of the Kalman filter, the observation errors will be incorrectly assigned to the other state variables during the estimating process, thereby damaging the calibration accuracy of traditional IMU sensor errors.

If we extend inner lever-arm parameters as the state variables, the calibration accuracy of IMU sensor errors will be improved unquestionably. However, $r_{xy}^{\,p},r_{xz}^p,r_{yx}^{\,p},r_{yz}^{\,p},r_{zx}^{\,p},r_{zy}^{\,p}$ cannot be fully excited by this rotation scheme. In order to obtain enough angular acceleration incentives, we improve the triangle wave into a sine wave (shown in figure 5, and referred to as Scheme B), which can provide variable angular rate to each axis of the IMU in sequence. By adopting Scheme B, the estimation of traditional IMU sensor errors will not be affected, and the observable degrees of $r_{xy}^{\,p},r_{xz}^{\,p},r_{yx}^{\,p},r_{yz}^{\,p},r_{zx}^{\,p},r_{zy}^{\,p}$ could be greatly increased as well.

Zoom In Zoom Out Reset image size

In addition, in [1] it is declared that inner lever-arm parameters can be calibrated when two gimbals simultaneously rotate with uniform angular rate for a tri-axis RINS [1]. The calibration results achieve a level of satisfaction, but the rotation scheme seems to lack sufficient explanation. In this paper, we try to analyze the error excitation mechanism of this rotation strategy. Actually, when two gimbals rotate with constant speed simultaneously, two gyros perpendicular to the inside rotation axis will respectively be sensitive to sine and cosine angular rate, and another gyro coincident with the inside rotation axis will be sensitive to uniform angular rate. Figure 6 provides the angular rates measured by the gyros when the inner and outer gimbals rotate simultaneously with uniform speed.

Zoom In Zoom Out Reset image size

Table 2 shows the angular rate types measured by the three gyros when two of the gimbals rotate simultaneously in a tri-axis RINS.

It can be seen that the IMU's three gyros can measure at least one variable angular rate input when the whole rotation scheme is finished, and thereby all inner lever-arm parameters can be effectively activated. However, due to the dismantling of the middle gimbal, a dual-axis RINS can only achieve simultaneous rotation of the inner and outer gimbals. Therefore, this rotation strategy is only suitable for the tri-axis RINS and it cannot be generalized to the dual-axis RINS, while the proposed Scheme B can be operated on both the dual-axis RINS and the tri-axis RINS, extending the universal applicability of the proposed self-calibration method.

When the IMU rotates along a fixed horizontal axis, most of the installation errors can be activated without any coupling errors, but there still remain exceptions. For instance, when the IMU rotates along the ${y_p}$ axis, $\alpha _{ay}^Z$ could be activated in the northern velocity error as

$$\begin{equation}\Delta {\dot V_N} = \Delta {\dot V_E}\sin (\alpha _{ay}^Z) \approx \alpha _{ay}^Zg\sin ({\omega _s}t)\end{equation} \tag{ 11 }$$

where ω_s is the angular rate, and the installation error like $\alpha _{gx}^Z$ is considered to be a small angle. $\Delta {\dot V_E}$ and $\Delta {\dot V_N}$ are the acceleration errors in the eastern and northern directions, respectively. At the same time, $\alpha _{gx}^Z$ could be excited in the northern velocity error as

$$\begin{equation}\left\{ \begin{array}{*{20}{l}} {\Delta {\omega _E} = \alpha _{gx}^Z{\omega _s}\cos \left( {{\omega _s}t} \right)}\\ {\Delta {{\dot V}_N} = g\Delta {\phi _E} = g\mathop \smallint \nolimits \Delta {\omega _E}dt = \alpha _{gx}^Zg\sin \left( {{\omega _s}t} \right).} \end{array}\right.\end{equation} \tag{ 12 }$$

From equations (11) and (12), it is obvious that $\alpha _{gx}^Z$ and $\alpha _{ay}^Z$ could be excited in the northern velocity with the same form. That is to say, if we take the northern velocity error as the observation with which to calibrate $\alpha _{gx}^Z$ and $\alpha _{ay}^Z$, these two parameters are coupled under this rotation. Nevertheless, when the IMU rotates along the ${z_p}$ axis in step 3, $\alpha _{ay}^Z$ could be activated independently, and these two variables could then be separated. Similarly, the groups of $\beta _{gx}^Y,\delta _{az}^X$ and $\beta _{gy}^X,\delta _{az}^Y$ are both coupled during the rotation along the ${z_p}$ axis, while $\delta _{az}^X$ and $\delta _{az}^Y$ could be singly activated when the IMU rotates along the ${x_p}$ and ${y_p}$ axes. In section 4.1, the simulation curves intuitively illustrate the decoupling of the relevant variables. As a result, all installation errors can be decoupled when the rotation scheme finishes.

Self-calibration is carried out through establishing the relationship between navigation errors and state variables. By observing navigation errors, the state variables can realize the optimal estimation under the Kalman filter algorithm. In the traditional self-calibration method, the velocity and attitude error equations of the RINS are used to construct the state equations of the Kalman filter:

$$\begin{align} \delta {{\dot V}^n} &= {f^{\,\,n}} \times {\phi ^n} - (2\omega _{ie}^n + \omega _{en}^n) \times \delta {V^n} + C_p^n\Delta f_{ip}^{\,\,p} \nonumber \\[2pt] {{\dot \phi }^n} &= - \omega _{ie}^n \times {\phi ^n} + (\delta \omega _{ie}^n + \delta \omega _{en}^n) + C_p^n\Delta \omega _{ip}^p \end{align} \tag{ 13 }$$

where $\delta {V^n} = {\left[ {\begin{array}{@{}lll@{}} {\delta {V_E}} & {\delta {V_N}} & {\delta {V_U}} \end{array}} \right]^T}$ represents the vector of velocity errors, ${f^{\,\,n}} = {\left[ {\begin{array}{@{}lll@{}} {{f_E}} & {{f_N}} & {{f_U}} \end{array}} \right]^T}$ denotes the vector of accelerations in the n-frame and ${\phi ^n} = {[\begin{array}{@{}lll@{}} {{\phi _E}} & {{\phi _N}} & {{\phi _U}} \end{array}]^T}$ is the vector of misalignment angles between the ideal n-frame and the computational navigation frame. $\omega _{ie}^{}$ is the angular rate vector of Earth. $\Delta f_{ip}^{\,\,p},\Delta \omega _{ip}^p$ denote the acceleration errors and angular rate errors caused by traditional IMU sensor errors respectively, which can be extended as equations (8) and (9). In previous research [8, 9, 15], the filter state vector is usually established by extending all the sensor errors as a large augmented vector, as follows,

$$\begin{equation}\begin{array}{l}X = \left[\delta V_E\,\,\delta V_N\,\,\delta V_U\,\,\phi_E\,\,\phi_N\,\,\phi_U\,\,\varepsilon_x\,\,\varepsilon_y\,\,\varepsilon_z\,\,\Delta_x\,\,\Delta_y\,\,\Delta_z\,\,\right.\\[4pt] \qquad K_{gx}\,\,K_{gy}\,\,K_{gz}\,\,K_{ax}\,\,K_{ay}\,\,K_{az}\\[4pt] \qquad \left.\alpha_{gx}^Z\,\,\alpha_{gy}^Z\,\,\beta_{gx}^Y\,\,\beta_{gy}^X\,\,\delta_{gz}^X\,\,\delta_{gz}^Y\,\,\alpha_{ay}^Z\,\,\delta_{az}^X\,\,\delta_{az}^Y\right]^T.\end{array}\end{equation} \tag{ 14 }$$

There is no doubt that this establishment can achieve the calibration of traditional IMU sensor errors, but it is flawed. On the one hand, if only one Kalman filter is established, the state vector dimension may burden the computer with heavy matrix calculations, which is not conducive to the real-time performance of the algorithm. On the other hand, the rotation strategy in this paper contains the transition processes of switching the IMU rotary axis by a 90° rotation. Due to unavoidable error approximation of the Kalman filter, the transition processes may cause instability of the state's convergencein, resulting calibration errors. In this paper, we design three separate Kalman filters for each rotation step; the filter state variables in each step are selected according to the law of error incentives, which can not only reduce the burden of calculation but can also ensure the stability of the convergence process.

The traditional self-calibration method without consideration of the inner lever-arm effect is named Method A, while Method B denotes the improved self-calibration method with consideration of the inner lever-arm effect. The differences between Method A and Method B are listed in table 3, including state vectors in three rotation steps, state equations and rotation scheme.

Table 3. The establishment of Method A and Method B.

	Method A	Method B
State vector in step 1	$\begin{gathered} {X_1} = [\delta {V_E}\delta {V_N}\delta {V_U}{\phi _E}{\phi _N}{\phi _U} \hfill\\ \qquad {\varepsilon _x}{\varepsilon _y}{\varepsilon _z}{\Delta _x}{\Delta _z}{K_{gy}}{K_{ax}}{K_{az}} \hfill\\ \delta _{gz}^X\,\, (\alpha _{gx}^Z - \alpha _{ay}^Z)\,\, \delta _{az}^Y{]^T} \end{gathered}$	$\begin{gathered} {X_1} = [\delta {V_E}\delta {V_N}\delta {V_U}{\phi _E}{\phi _N}{\phi _U} \hfill \\ \qquad {\varepsilon _x}{\varepsilon _y}{\varepsilon _z}{\Delta _x}{\Delta _z}{K_{gy}}{K_{ax}}{K_{az}} \hfill \\ \qquad \delta _{gz}^X{\text{ (}}\alpha _{gx}^Z-\alpha _{ay}^Z{\text{) }}\delta _{az}^Yr_{xx}^{\,p}r_{xz}^{\,p}r_{zx}^{\,p}r_{zz}^{\,p}{]^T} \hfill \\ \end{gathered}$
State vector in step 2	$\begin{gathered} {X_2} = [\delta {V_E}\delta {V_N}\delta {V_U}{\phi _E}{\phi _N}{\phi _U} \hfill \\ \qquad {\varepsilon _x}{\varepsilon _y}{\varepsilon _z}{\Delta _y}{\Delta _z}{K_{gx}}{K_{ay}}{K_{az}} \hfill \\ \qquad \alpha _{gy}^Z\delta _{gz}^Y\delta _{az}^X{]^T} \hfill \\ \end{gathered}$	$\begin{gathered} {X_2} = [\delta {V_E}\delta {V_N}\delta {V_U}{\phi _E}{\phi _N}{\phi _U} \hfill \\ \qquad {\varepsilon _x}{\varepsilon _y}{\varepsilon _z}{\Delta _y}{\Delta _z}{K_{gx}}{K_{ay}}{K_{az}} \hfill \\ \qquad \alpha _{gy}^Z\delta _{gz}^Y\delta _{az}^Xr_{yy}^{\,p}r_{yz}^{\,p}r_{zy}^{\,p}r_{zz}^{\,p}{]^T} \hfill \\ \end{gathered}$
State vector in step 3	$\begin{gathered} {X_{\text{3}}} = [\delta {V_E}\delta {V_N}\delta {V_U}{\phi _E}{\phi _N}{\phi _U} \hfill \\ \qquad {\varepsilon _x}{\varepsilon _y}{\varepsilon _z}{\Delta _x}\;{\Delta _y}{K_{gz}}{K_{ax}}{K_{ay}} \hfill \\ \qquad {\text{ (}}\beta _{gx}^Y-\delta _{az}^X{\text{) (}}\beta _{gy}^X-\delta _{az}^Y)\alpha _{ay}^Z{]^T} \hfill \\ \end{gathered}$	$\begin{gathered} {X_{\text{3}}} = [\delta {V_E}\delta {V_N}\delta {V_U}{\phi _E}{\phi _N}{\phi _U} \hfill \\ \qquad {\varepsilon _x}{\varepsilon _y}{\varepsilon _z}{\Delta _x}\;{\Delta _y}{K_{gz}}{K_{ax}}{K_{ay}} \hfill \\ \qquad {\text{ (}}\beta _{gx}^Y-\delta _{az}^X{\text{) (}}\beta _{gy}^X-\delta _{az}^Y)\alpha _{ay}^Zr_{xx}^{\,p}r_{xy}^{\,p}r_{yx}^{\,p}r_{yy}^{\,p}{]^T} \hfill \\ \end{gathered}$
State equations	$\begin{gathered} \delta {{\dot V}^{\,\,n}} = {f^{\,\,n}} \times {\phi ^n} - (2\omega _{ie}^n + \omega _{en}^n) \times \delta {V^{\,\,n}} + C_p^{\,\,n}\Delta f_{ip}^{\,\,p} \hfill \\ {{\dot \phi }^n} = - \omega _{ie}^n \times {\phi ^n} + (\delta \omega _{ie}^n + \delta \omega _{en}^n) + C_p^n\Delta \omega _{ip}^p \hfill \\ \end{gathered}$	$\begin{gathered} \delta {{\dot V}^{\,\,n}} = {f^{\,\,n}} \times {\phi ^n} - (2\omega _{ie}^n + \omega _{en}^n) \times \delta {V^{\,\,n}} + C_p^n\Delta f_{ip}^{\,\,p} + C_p^{\,\,n}\Delta f_{ip\_r}^{\,\,p} \hfill \\ {{\dot \phi }^n} = - \omega _{ie}^n \times {\phi ^n} + (\delta \omega _{ie}^n + \delta \omega _{en}^n) + C_p^n\Delta \omega _{ip}^p \hfill \\ \end{gathered}$
Rotation scheme	Bi-directional rotation scheme with	Bi-directional rotation scheme with
	uniform angular rate (figure 4)	variable angular rate (figure 5)

As shown in table 3, inner lever-arm errors are not considered in state vectors and state equations in Method A, and the bi-directional rotation scheme with uniform angular rate (figure 4) is adopted. In the improved Method B, the corresponding inner lever-arm errors are included in the state vectors according to the excitation sequence. Simultaneously, the inner lever-arm errors are added to the velocity error model of the state equations. In order to fully stimulate the inner lever-arm errors, the bi-directional rotation scheme with variable angular rate (figure 5) is adopted. According to the analysis in section 3.2, when the IMU rotates along the ${y_p}$ axis in step 1, $\alpha _{gx}^Z$ and $\alpha _{ay}^Z$ need to be considered as a united variable. When the IMU rotates along the ${x_p}$ axis in step 2, all filter states can be separated during the estimation process. In the final step, step 3, $\beta _{gx}^Y,\delta _{az}^X$ and $\beta _{gy}^X,\delta _{az}^Y$ are both coupled, and need to be considered as united variables.

In each step of the calibration process, the velocity errors in three directions are chosen as observations, denoted by

$$\begin{equation}Z={\lbrack\delta V_E\delta V_N\delta V_U\rbrack}^T.\end{equation} \tag{ 15 }$$

The filter model including state equations and observation equations can be described as

$$\begin{equation}\begin{gathered} \dot X = FX + W \hfill \\ Z = HX + V \hfill \\ \end{gathered} \end{equation} \tag{ 16 }$$

where $F$ denotes the state transition matrix, $H$ stands for the observation matrix and $G$ denotes the driving noise matrix. $W$ and $V$ represent system noise and measurement noise respectively, which are both zero-mean Gaussian white-noise processes. The elements in each matrix are provided in the appendix at the end of this article. From equations (13) and (15), the state equations and measurement equations can be approximated as linear, so the linear Kalman filter algorithm could be adopted to estimate the state variables [9, 10]. The proposed Kalman filter is suitable for the three steps of the self-calibration method. In each step, the difference is that the initial values of the covariance matrix should be set according to the error parameters of current excitations (refer to table 3), while the corresponding initial values in the state covariance matrix are set based on the average level of each state variable under the statistical characteristics.

Simulations are conducted to compare the performance of Method A with that of Method B. The simulation parameters are determined according to the practical characteristics in table 1. The frequency of navigation calculation is set as 200 Hz, and the frequency of the Kalman filter is 1 Hz. The given values of the relative sensor errors are listed in table 4. The convergence curves from adopting Method A and Method B are provided in figures 7(a)–(c), and the estimation curves of the inner lever-arm parameters with Method B are shown in figure 7(d). Simulations are performed 20 times for both methods, and the estimation errors of relevant parameters are listed in table 4.

Table 4. Estimation errors of the self-calibration with Method A and Method B.

		Estimation errors
Error parameters	Given values	Method A	Method B
${\varepsilon _x}$ (° h−1)	0.02	0.0182	0.0039
${\varepsilon _y}$ (° h−1)	0.02	0.0071	−0.0016
${\varepsilon _z}$ (° h−1)	0.02	0.0086	0.0065
${\Delta_x}$ (${\mu }g$)	50	42.38	1.61
${\Delta _y}$ ($\mu g$)	50	51.56	2.47
${\Delta _z}$ (${\mu }g$)	50	38.59	−2.26
${K_{gx}}$ (ppm)	50	−17.54	−0.40
${K_{gy}}$ (ppm)	50	−15.68	1.67
${K_{gz}}$ (ppm)	50	16.71	−0.44
${K_{ax}}$ (ppm)	50	−2.64	2.40
${K_{ay}}$ (ppm)	50	−2.88	−1.59
${K_{az}}$(ppm)	50	5.21	−1.69
$\alpha _{gx}^Z$ ('')	50	6.82	0.04
$\alpha _{gy}^Z$ ('')	50	−0.43	−0.07
$\beta _{gx}^Y$ ('')	50	2.18	−0.47
$\beta _{gy}^X$ ('')	50	−5.01	−1.45
$\delta _{gz}^X$ ('')	50	−1.46	−0.07
$\delta _{gz}^Y$ ('')	50	−1.41	−0.01
$\alpha _{ay}^Z$('')	50	−10.32	−1.04
$\delta _{az}^X$ ('')	50	15.29	−0.52
$\delta _{az}^Y$ ('')	50	−13.69	−1.50
$r_{xx}^{\,p}$ (mm)	50	–	2.54
$r_{xy}^{\,p}$ (mm)	50	–	0.82
$r_{xz}^{\,p}$ (mm)	50	–	−3.22
$r_{yx}^{\,p}$ (mm)	50	–	−1.78
$r_{yy}^{\,p}$ (mm)	50	–	−1.19
$r_{yz}^{\,p}$ (mm)	50	–	2.35
$r_{zx}^{\,p}$ (mm)	50	–	0.80
$r_{zy}^{\,p}$ (mm)	50	–	−0.56
$r_{zz}^{\,p}$ (mm)	50	–	−2.74

Zoom In Zoom Out Reset image size

We can learn from figures 7(a)–(c) that the estimation curves of traditional IMU sensor errors using Method B can approximately reach the given value, while the estimation results of certain state variables using Method A has constant errors, demonstrating that the observation errors caused by the inner lever-arm effect will be wrongly assigned to other IMU sensor errors if the inner lever-arm errors are not considered in the model. Table 4 indicates that the maximum estimation error of IMU sensor errors is less than 5% of the given values. Furthermore, figure 7(d) indicates that the estimation of inner lever-arm parameters could approximately reach the given values using Method B as well, and the maximum estimation error is less than 4 mm. In summary, the proposed Method B can not only improve the estimation accuracy of IMU sensor errors, but can also realize a satisfactory calibration of inner lever-arm parameters.

In military applications, the dual-axis RINS can operate the self-calibration process before navigation without disassembling from the vehicle, which can greatly improve the carrier maneuverability. In this section, a long-term vehicle navigation experiment is designed according to actual situations, and the vehicle experimental environment is shown in figure 8. The dual-axis RINS is fixed on the mounting bracket with the system's body frame approximately coincident with the vehicle's right-forward-upward coordinate. A GPS antenna is installed on top of the vehicle to provide velocity and position references, with a velocity accuracy of 0.1 m s⁻¹ and a position accuracy of 10 m in single-point location mode. The data output frequency of the GPS receiver is 10 Hz. The RINS is powered by a 27 V DC lead battery, and relevant data are collected by a laptop inside the vehicle.

Zoom In Zoom Out Reset image size

The trajectory of the vehicle experiment is shown in figure 9. In order to activate more errors, mount climbing and descending, S-curves, acceleration and deceleration are all involved. During navigation, the rotation modulation strategy mainly contains bi-directional rotation along the IMU vertical axis, periodically changing direction along the IMU's horizontal axis. Before the vehicle is driven, self-calibration Methods A and B are executed once in sequence, and the calibration results of the error parameters are stored. After loading the calibration parameters obtained by Method B into the dual-axis RINS, the vehicle experiment is started. During driving, the real-time navigation results are recorded, and the raw data, including the output of gyroscopes, accelerometers and encoders and GPS data, are collected simultaneously. After the vehicle experiment is completed, the calibration results obtained using Method A are loaded for offline navigation processing.

Zoom In Zoom Out Reset image size

The navigation errors in figure 10 are obtained through subtraction between the dual-axis RINS and the GPS, separately compensating for the self-calibration results with Methods A and B, followed by the statistical results of maximum navigation errors and circular eprobable (CEP) (table 5). In order to reduce the spikes introduced by the GPS signal, the navigation results provided in figure 10 are averaged every 10 s.

Zoom In Zoom Out Reset image size

In table 5, CEP is an important index to characterize the positioning performance of an INS. Equations (17) and (18) provide the CEP calculation method at a 50% circle probability radius:

$$\begin{equation}\Delta P({T_i}) = \frac{1}{{{T_i}}}\sqrt {{{(\Delta {P_e}({T_i}))}^2} + {{(\Delta {P_n}({T_i}))}^2}} \end{equation} \tag{ 17 }$$

$$\begin{equation}CEP=0.83\sqrt{\frac1n\sum_{i=1}^n{(\Delta P{(T_i)}^2)}}\end{equation} \tag{ 18 }$$

where ${T_i}$ is a navigation timestamp sequence, such as 10 s, 20 s, 30 s, etc and $\Delta {P_e}({T_i})$ and $\Delta {P_n}({T_i})$ represent the eastern and northern positioning error sequence at ${T_i}$. $\Delta P({T_i})$ represents the radial error rate at ${T_i}$. $n$ represents the number of sampling points throughout the vehicle experiment. From figure 10, we can see that the navigation performance has been greatly enhanced by adopting the improved self-calibration method. The maximum velocity and positioning errors in table 5 have been reduced by approximately half as well. The CEP of the vehicle navigation has reduced from 0.33 n mile h⁻¹ to 0.25 n mile h⁻¹.