A Self-Calibration Method for Dual-axis Rotational Inertial Navigation Systems

2 MATHEMATICAL MODEL OF ROTATIONAL INS

2.3 Mathematical Model of Rotational INS

For this model of rotational INS, let dθ be the non-orthogonal between the rotation axes. As shown in Figure 1, the vector of the outer rotation axis in its initial position can be expressed in the s-frame as $x_r^s=\left[\begin{array}{ccc}\cos \left(d\theta \right)& 0& -\sin \left(d\theta \right)\end{array}\right]$, and the rotation matrix $R_x^s\left({\theta }_x\right)$ around this non-orthogonal outer axis can be represented as:

$R_x^s\left({\theta }_x\right)=\left[\begin{array}{ccc}1& d\theta \sin \left({\theta }_x\right)& d\theta \left(\cos ({\theta }_x\right)-1)\\ -d\theta \sin \left({\theta }_x\right)& \cos \left({\theta }_x\right)& -\sin \left({\theta }_x\right)\\ d\theta \left(\cos ({\theta }_x\right)-1)& \sin \left({\theta }_x\right)& \cos \left({\theta }_x\right)\end{array}\right]$5

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 1

Definition of the non-orthogonal angle

The corresponding rotation matrix around the inner axis is given by:

$R_z^s\left({\theta }_z\right)=\left[\begin{array}{ccc}\cos \left({\theta }_z\right)& -\sin \left({\theta }_z\right)& 0\\ \sin \left({\theta }_z\right)& \cos \left({\theta }_z\right)& 0\\ 0& 0& 1\end{array}\right]$6

where θx and θz in the equations above denote the angles of the outer and inner rotation axes, respectively, relative to their initial angles.

The transformation matrix from the s-frame to the body frame can be expressed as:

$C_s^b=R_x^s\left({\theta }_x\right)R_z^s\left({\theta }_z\right)$7

The vehicle’s attitude angles α0, β0, and γ0 represent the rotation of the vehicle frame relative to the navigation frame. To align the outer gimbal axis with the vehicle’s x-axis, the Euler angle sequence is defined as z-y-x. The corresponding transformation from the b-frame to the n-frame is illustrated in Figure 2, where the colors of the rotation angles represent the colors of the respective rotation axes. Mathematically, this transformation is defined as:

$C_b^n=R_z^b\left(\gamma 0\right)R_y^b\left(\beta 0\right)R_x^b\left(\alpha 0\right)$8

where $R_z^b\left(\gamma 0\right)$, $R_y^b\left(\beta 0\right)$, and $R_x^b\left(\alpha 0\right)$ are the rotation matrices around the x-, y-, and z-axes of the body frame, respectively.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 2

Definition of the attitude angles between the navigation frame and body frame

When the inner and outer frames of the gimbals rotate at angular rates ω_x and ω_z, respectively, the s-frame angular velocity and specific force can be expressed as:

${\omega }_{is}^s=C_b^s\left(C_n^b\left({\omega }_{ie}^n+{\omega }_{en}^n\right)+{\omega }_x^b\right)+{\omega }_z^s$9

$f^s=C_b^s{C}_n^b(a^n+(2{\omega }_{ie}^n+{\omega }_{en}^n)\times v^n-g^n)$10

where ${\omega }_{ie}^n$ and ${\omega }_{en}^n$ are rotation rates of the Earth and the vehicle, respectively, projected onto the n-frame; ${\omega }_x^b={\left[\begin{array}{ccc}{\omega }_x& 0& 0\end{array}\right]}^T$ and ${\omega }_z^s={\left[\begin{array}{ccc}0& 0& {\omega }_z\end{array}\right]}^T$; gⁿ is the acceleration due to gravity, and vⁿ is the velocity of the vehicle. When the vehicle is stationary, ${\omega }_{en}^n=0$ and vⁿ = 0. By substituting Equations (9) and (10) into Equations (1) and (2), the outputs of the inertial sensors in the stationary state can be expressed as:

$N_g=(I+K_g)(I+E_g)(C_b^s(C_n^b{\omega }_{ie}^n+{\omega }_x^b)+{\omega }_z^s)+\epsilon$11

$N_a=-\left(I+K_a\right)\left(I+E_a\right)C_b^s{C}_n^b{g}^n+\nabla$12

3 EFFECTS OF NON-ORTHOGONAL ANGLES ON SYSTEM ERRORS

3.2 Effects of Non-orthogonality on Navigation

The effect of non-orthogonal rotation errors on navigation performance depends on the coordinate system used for computation and navigation. When navigation calculations are performed in the body frame, the IMU outputs must first be transformed into the body coordinate system. The transformation equations for angular velocity are:

${\omega }_{ib}^b=C_s^b\left({\omega }_{is}^s-{\omega }_z^s\right)-{\omega }_x^s$30

${\tilde{\omega }}_{ib}^b={\tilde{C}}_s^b({\omega }_{is}^s-{\omega }_z^s)-{\omega }_x^s$31

The resulting error introduced by the gimbal’s non-orthogonality into the body frame’s angular velocity output is:

$\delta {\omega }_{ib}^b={\tilde{\omega }}_{ib}^b-{\omega }_{ib}^b=\delta C_s^b\left({\omega }_{is}^s-{\omega }_z^s\right)=\delta C_s^b{\omega }_{ib}^s$32

where $\delta C_s^b={\tilde{C}}_s^b-C_s^b$ represents the rotation error matrix caused by non-orthogonality of the gimbal axes. The corresponding non-orthogonal error in the accelerometer outputs (represented in the body frame) is:

$f^b={\tilde{f}}^b-f^b=\delta C_s^b{f}^s$33

These expressions indicate that IMU output errors in the body frame are functions of both the gimbal-induced transformation errors and body-frame inertial inputs. In particular, higher rotation rates amplify the effect of non-orthogonal errors. Note that these derivations assume that the IMU and gimbal coordinate systems are perfectly aligned; any misalignment would cause the non-orthogonal errors to introduce more complex and pronounced effects (Deng et al., 2017).

To mitigate gimbal-induced navigation errors, modern navigation frameworks compute state estimates directly in the IMU coordinate system and then transform the estimated attitude into the body frame. In this approach, gimbal non-orthogonality only affects attitude estimates. The transformation from the IMU frame to the navigation frame is given by:

${\tilde{C}}_b^n=C_s^n{\tilde{C}}_b^s=C_s^n{R}_z\left({\theta }_z\right)R_x\left({\theta }_x\right)$34

The corresponding attitude error matrix becomes:

$\delta C_b^n=C_n^b{\tilde{C}}_b^n=C_s^b{\tilde{C}}_b^s=\left[\begin{array}{ccc}1& -d\theta \sin \left({\theta }_x\right)& 2d\theta {\sin }^2\left(\frac{{\theta }_x}{2}\right)\cos \left({\theta }_z\right)\\ d\theta \sin \left({\theta }_x\right)& 1& 0\\ 2d\theta {\sin }^2\left(\frac{{\theta }_x}{2}\right)\cos \left({\theta }_z\right)& 0& 1\end{array}\right]$35

The estimated attitude error vector can then be derived from this transformation matrix:

$att={\left[\begin{array}{ccc}\arctan \left(\frac{C_n^b\left(3,2\right)}{C_n^b\left(3,3\right)}\right)& \mathrm{asin}\left(-C_n^b\left(3,1\right)\right)& \arctan \left(C_n^b\left(2,1\right)/C_n^b\left(1,1\right)\right)\end{array}\right]}^T$ 36

The attitude error introduced by the non-orthogonal angle is therefore:

$\delta att={[0\hspace{0.17em}-2d\theta {\sin }^2(\frac{{\theta }_x}{2}\left)\cos \right({\theta }_z)\hspace{0.17em}d\theta \sin ({\theta }_x\left)\right]}^T$37

To evaluate the practical effect of gimbal non-orthogonality, we rotated a gimbal mechanism following a standard 16-sequence rotation scheme (Yuan et al., 2012). For this analysis, we neglected the IMU’s inherent systematic errors to isolate the effects of the gimbal. For a non-orthogonal angle error of 20 arcseconds, the induced attitude errors are illustrated in Figure 3. The results demonstrate that gimbal-induced attitude errors depend solely on angular displacement, with the maximum error magnitude approximately equal to the non-orthogonal angle. Furthermore, the carrier’s attitude errors are periodic, corresponding to the gimbal’s periodic rotation. Performing navigation directly in the IMU coordinate system thus limits the propagation of non-orthogonal errors into velocity and acceleration estimates, effectively confining their influence to attitude determination.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 3

Effect of gimbal non-orthogonality on carrier attitude errors

4 MULTI-POSITION CALIBRATION METHOD

Our proposed calibration procedure to account for non-orthogonality consists of three key steps. First, the outputs of each axis (in multiple positions) are used to calibrate the accelerometer error parameters, gyroscope biases, and the non-orthogonal angles between rotation axes. Second, the left gyroscope error parameters are calibrated using turntable rotations. Third, the measured parameters are used to adjust the IMU outputs via reverse error correction.

In the case of an unknown initial attitude, initial calibration is performed with 12 different sensor orientations, shown in Figure 4. Each of the three IMU axes is placed in two symmetric positions, one pointing upward and one pointing downward relative to the inner gimbal axes. All positions follow the inner frame symmetry.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 4

Multi-position calibration scheme

Because the body frame’s attitude cannot be determined, it cannot be assumed to be negligible. However, α0 is the initial angle between the outer gimbal axis and the vehicle body’s x-axis, and the outer axis angle can be manually adjusted to keep α0 within a small range. A rough estimate of α0 can therefore be obtained by rotating the outer axis by 90° and using:

$\tan \left(\alpha 0\right)\approx \frac{N_z\left(\theta x={90}^{\circ },\theta z=0\right)}{N_z(\theta x=0,\theta z=0)}=\frac{g\cos \left(\beta 0\right)\sin \left(\alpha 0\right)+{\nabla }_z}{g\cos \left(\beta 0\right)\cos \left(\alpha 0\right)+{\nabla }_z}$38

The accuracy of this estimate depends on the accelerometer bias and is lowest when α0 → 0. Assuming α0 = 0 and an accelerometer bias on the order of 1 mg, the estimate for α0 from Equation (38) is 0.057°. In contrast, β0 represents the angle between the outer gimbal axis and the horizontal plane. Because the outer gimbal is fixed to the carrier, β0 cannot be manually adjusted. In this setup, α0, dθ, and E_aij can therefore be considered small quantities, while β0 cannot.

The initial solution for the systematic error obtained through the least squares estimation method follows Equations (24)–(29) in Section 3.1. The compensated output of the accelerometers using the rough measured error parameters is as follows:

${\left((I+{\tilde{K}}_a)(I+{\tilde{E}}_a)\right)}^{-1}{N}_a=\frac{1+d\theta \hspace{0.17em}\tan \left(\beta 0\right)}{\cos \left(\beta 0\right)}{f}^s+{\left((I+{\tilde{K}}_a)(I+{\tilde{E}}_a)\right)}^{-1}\nabla$39

After this rough calibration, the accelerometer outputs are proportional to the true s-frame inputs, except for bias errors. However, due to inner gimbal misalignment, the x- and y-axis outputs vary as the inner gimbal rotates. As shown in Equation (20), this effect on f^s is maximized in symmetric positions.

Let L denote the difference between roughly calibrated accelerometer outputs under symmetric inner frame rotations. The expression for L is as follows:

$L=\left[\begin{array}{c}{L}_x\\ L_y\end{array}\right]=\left[\begin{array}{cc}{I}_{2\times 2}& 0\end{array}\right]{\left((I+{\tilde{K}}_a)(I+{\tilde{E}}_a)\right)}^{-1}(N_a({\theta }_x,{\theta }_z)-N_a({\theta }_x,{\theta }_z+\pi ))$40

With the attitudes used for the calibration above, we obtain:

$L^T=XA$41

In Equation (41), A is a function of the gimbal angles θx and θz:

$A={\left[\begin{array}{ccc}-2\hspace{0.17em}\cos \left({\theta }_z\right)& 2\hspace{0.17em}\sin \left({\theta }_z\right)\cos \left({\theta }_x\right)& \left(1-\cos \left({\theta }_x\right)\right)\cos \left({\theta }_z\right)\\ 2\hspace{0.17em}\sin \left({\theta }_z\right)& 2\hspace{0.17em}\cos \left({\theta }_z\right)\cos \left({\theta }_x\right)& \left(1-\cos \left({\theta }_x\right)\right)\sin \left({\theta }_z\right)\end{array}\right]}^T$42

X represents the parameters to be measured:

$X=\left[\begin{array}{ccc}tan\left(\beta 0\right)& \alpha 0& d\theta \end{array}\right]$43

The least-squares solution for the components of matrix X is therefore:

$X=L^T{A}^T{\left(A{A}^T\right)}^{-1}$44

After α0, β0, and dθ have been estimated, they can be substituted into Equations (21)–(26) to obtain refined estimates of the accelerometer scale factors, installation errors, and bias errors, thus completing the calibration.

5 SIMULATION RESULTS AND ANALYSIS

We conducted a series of simulations to validate the effectiveness of the proposed calibration scheme. The attitude of the carrier was set at [0.1° 0.1° 20°], and the initial location was set at 120.664265°E, 30.686651°N. The true values of all error parameters are shown in Table 2. The random walk of accelerometer velocity was set at 20ug/√Hz, and the angular random walk of the gyroscope was 0.02°/√h.

The simulation framework comprises the following steps:

1. Initialization: The initial attitude of the carrier is used to derive the static angular velocity and acceleration components in the turntable coordinate system.

2. Rotational sequence design: A rotation sequence is then designed. The turntable angles and rotational velocities are substituted into Equations (11) and (12) to compute the angular velocity and acceleration of the inner frame.

3. Error propagation: The IMU error parameters are then used to calculate the actual IMU outputs.

To assess the proposed scheme’s ability to estimate non-orthogonal angle errors, the value of the non-orthogonal angle was gradually increased from 1” to 40”. For each value, 200 Monte Carlo simulation experiments were conducted, and calibration accuracy was assessed using the relative standard deviation (RSD). Simulation results are shown in Figure 5. Notably, due to the limitations of device precision, the RSD falls below 10% once the non-orthogonal angle exceeds 7”.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 5

Relative standard deviation error of the non-orthogonal angle estimated via simulation

To further validate the proposed calibration method, we simulated and compared three different calibration schemes:

Scheme 1: The proposed method, which calibrates IMU errors and nonorthogonality errors of the rotary table.

Scheme 2: Traditional six-position calibration, which calibrates IMU errors but does not account for rotary table non-orthogonality.

Scheme 3: System-level calibration for IMU errors, which involves a discrete calibration approach for the rotation matrices between (i) the IMU and the inner axes of the rotary table and (ii) the inner and outer axes of the rotary table.

To comprehensively evaluate each scheme’s calibration performance, 200 Monte Carlo simulations were conducted for Schemes 1 and 2. Calibration accuracy was measured as the bias between the ideal value and the mean calibration result across all simulations, while the calibration uncertainty was quantified as the corresponding the standard deviation. For Scheme 3, the bias between the ideal value and a single two-hour calibration result was used as the performance metric. This was done because system-level calibration requires longer computational time, but the higher noise resistance results in high repeatability. To ensure errors were comparable among schemes, the installation errors obtained from Scheme 3 (which were initially referenced to the IMU coordinate system) were transformed into the rotary table coordinate system using known rotation matrices between the IMU and the inner/outer axes of the rotary table.

Figure 6 shows the convergence of the estimated gyroscope error parameters over the Kalman filtering process. The zero bias and scale factor converge rapidly, while the installation misalignment converges more slowly and with less stability. The true and calibrated errors for zero bias and scale factor are summarized in Table 2, whereas those for installation misalignment and non-orthogonality angle are provided in Table 3.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 6

Estimation of gyroscope errors around (a) zero bias, (b) scale factor, and (c) installation errors calculated via Kalman filtering

As shown in Tables 2 and 3, the proposed scheme (Scheme 1) achieves higher estimation accuracy for the z-axis accelerometer zero bias than the other two schemes. The traditional discrete calibration scheme (Scheme 2) exhibits lower calibration precision for the accelerometer scale factor, primarily due to its failure to compensate for rotary table non-leveling errors. While these non-leveling and non-orthogonality errors also affect the precision of installation misalignment errors, the proposed scheme still outperforms the system-level calibration approach (Scheme 3). The lower performance of Scheme 3 is due to the complex nature of error propagation inherent in system calibration, which results in lower observability of installation errors.

In summary, the traditional discrete calibration scheme (Scheme 2) suffers significant accuracy loss even when the turntable’s leveling error and axis non-orthogonality are relatively small. Leveling errors can be suppressed by introducing additional inner-loop symmetric positions, but residual errors require further compensation. System-level calibration (Scheme 3) performs well for accelerometer scale factors but not for installation errors, largely because this approach relies on velocity and position errors as observables. Although system-level calibration is also more tolerant to measurement noise, its effectiveness depends on correct noise parameters, and it does not yield a substantial improvement in scale factor accuracy. In contrast, the discrete calibration scheme proposed here (Scheme 1) improves overall calibration precision and reduces measurement uncertainty by explicitly calibrating non-orthogonal angle errors and compensating for other error sources. The proposed approach also avoids the need to empirically estimate device noise parameters, thus enhancing its practical robustness and ease of implementation.

6 EXPERIMENTAL RESULTS AND ANALYSIS

Zheng et al. (2023) conducted a detailed analysis of how sensor errors with varying precision affect navigation accuracy. They identified gyroscope zero bias and installation as the primary contributors to heading angle errors, while accelerometer zero bias and installation errors were the main sources of velocity errors. The effects of installation and scale factor errors were state-dependent, with systematic errors scaling proportionally with the intensity of the motion.

To empirically demonstrate the effectiveness of the proposed calibration scheme, different calibration approaches were experimentally tested by fixing a fiber-optic INS to a small-scale dual-axis reorientation mechanism. The system structure is shown in Figure 7. Because the proposed self-calibration method does not require the carrier’s attitude information, the turntable was not leveled or aligned with the north direction.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 7

experimental setup with a static base

Prior to the experiment, the performance of the IMU and the angular stability of the rotary table were assessed using Allan variance analysis. The results of this analysis are shown in Figure 8 for the gyroscope, accelerometer, and the inner/ outer axes of the rotary table. The 10-second average standard deviation of the gyroscope and accelerometer are approximately 0.08°/h and 4 μg, respectively. For the rotary table, the standard deviation declined from 0.0028 arcseconds over 10 seconds to 0.0004 arcseconds over 100 seconds, demonstrating excellent long-term angular position stability. Overall, these results confirm that the vibrations of the rotary table have negligible effects on the 100-second averages of the gyroscope and accelerometer outputs.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 8

Allan variance analysis about the (a) gyroscope, (b) accelerometer, and (c) rotary table

Each calibration attitude was maintained for 100 seconds, with a gimbal rotation speed of 4°/s. As with the simulation results, the calibration parameters exhibited only minor variation, so only representative parameters are listed in Table 4.

After calibration, we assessed the practical effectiveness of our proposed scheme by conducting a six-hour navigation test under rotation modulation conditions. Of the potential rotation schemes, we chose a 16-sequence rotation that is widely used in practical applications (Ji et al., 2013; Yuan et al., 2012; Zhou et al., 2008). The INS was then subjected to periodic rotations following the procedure described by Yuan et al. (2012). Compared to compact modulated INS systems, our experimental setup exhibited significantly larger internal lever arm and time synchronization errors. These errors were therefore calibrated and compensated prior to the navigation test. The results from the six-hour navigation test are shown in Figures 9–11. As shown in the figures, the proposed scheme (Scheme 1) demonstrates higher navigation accuracy than the traditional approach (Scheme 2).

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 9

Attitude errors under rotation modulation

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 10

Velocity errors under rotation modulation

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 11

Position errors under rotation modulation

The navigation results presented in Figures 9–11 further demonstrate that the proposed calibration scheme significantly enhances navigation accuracy under dynamic conditions like rotation modulation. Consistent with the simulation results, the attitude errors exhibit periodic oscillations corresponding to the rotation cycle. These oscillations arise because the alternative calibration schemes (Schemes 2 and 3) do not compensate for non-orthogonality errors, highlighting the importance of correcting for these errors.

Titterton and Weston (2004) evaluated how equivalent bias errors from gyroscopes and accelerometers during static-base inertial navigation propagate in a manner that induces these periodic oscillations. These oscillations include the

Schuler period $T_s=\frac{2\pi }{{\omega }_s}=84.4\hspace{0.17em}\min ,$ the Earth’s rotation period $T_e=\frac{2\pi }{{\omega }_{ie}}=24h,$, and the Foucault period $T_f=\frac{2\pi }{{\omega }_{ie}\sin L}$. The Schuler oscillations are caused by attitude misalignment between the true and computed geographic frames. The two horizontal components of the attitude error, combined with gravity, form a second-order negative feedback loop. This feedback loop produces stable oscillations that are immune to dynamic accelerations due to Schuler tuning. The Earth rotation-induced tuning. The Earth rotation-induced oscillations result from coupling between tri-axial attitude errors and latitude estimation errors, which introduce components of Earth’s rotation into the inertial solution. Specifically, errors in the north- and zenith-pointing gyroscopes induce oscillations in the heading angle errors, constant biases in the east velocity, and latitude errors. In contrast, errors in the east-pointing gyroscope introduce a constant heading bias, while gyroscope bias in the north direction leads to oscillations in the north velocity errors. The Foucault oscillations stem from incomplete compensation of parasitic accelerations.

The amplitude of the oscillatory navigation errors thus serves as a critical indicator of calibration quality, and Table 5 accordingly shows these amplitudes for the three tested calibration schemes. As shown in the table, the calibration scheme proposed here (Scheme 1) offers the most stable attitude estimation. Over the six-hour test, the roll, pitch, and yaw errors under Scheme 1 remain low and increase only gradually with time, indicating that these scheme effectively compensates for gyroscope bias. The velocity and position errors correspondingly exhibit minimal periodic oscillation. Most notably, the north position error remains below 1 km throughout the experiment, which reflects not only adequate compensation for accelerometer bias but also effective suppression of the Schuler resonance amplification.

In contrast, Schemes 2 and 3 exhibit more pronounced oscillations during the six-hour navigation test. Over time, both velocity and position errors begin to show increasingly periodic fluctuations with growing amplitudes. Of the two schemes, Scheme 3 exhibits smaller velocity oscillations than Scheme 2, suggesting that it more effectively compensates for the equivalent accelerometer bias. However, both schemes fall short of the performance achieved by Scheme 1.