A Three-Stage Accelerometer Self-Calibration Technique for Space-Stable Inertial Navigation Systems

Abstract

As a specific force sensor, the tri-axis accelerometer is one of the core instruments in an inertial navigation system (INS). During navigation, its measurement error directly induces constant or alternating navigation errors of the same order of magnitude. Moreover, it also affects the estimation accuracy of gyro drift coefficients during the initial alignment and calibration, which will indirectly result in navigation errors accumulating over time. Calibration can effectively improve measurement accuracy of the accelerometer. Device-level calibration can identify all of the parameters in the error model, and the system-level calibration can accurately estimate part of these parameters. Combining the advantages of both the methods and making full use of the precise angulation of the space-stabilized platform, this paper proposes a three-stage accelerometer self-calibration technique that can be implemented directly in the space-stable INS. The device-level calibration is divided into two steps considering the large amount of parameters. The first step is coarse calibration, which identifies parameters except for the nonlinear terms, and the second step is fine calibration, which not only identifies the nonlinear parameters, but also improves the accuracy of the parameters identified in the first step. The follow-on system-level calibration is carried out on part of the parameters using specific force error and attitude error to further improve the calibration accuracy. Simulation result shows that by using the proposed three-stage calibration technique in the space-stable INS, the estimation accuracy of accelerometer error can reach ${1\times 10}^{-6}\mathrm{g}$ order of magnitude. Experiment results show that after the three-stage calibration, the accuracy of latitude, longitude, and attitude angles has increased by over 45% and the accuracy of velocity has increased by over 22% during navigation.

1. Introduction

Inertial navigation systems (INSs) tend to be the first choice for marine navigation because they are able to output motion parameters completely, autonomously, and continuously [1,2,3]. In order to fulfill the requirement for long-term and high-precision marine navigation, it is essential for INSs to employ high-precision inertial instruments. The gyroscope should be able to work stably for a long time, of which the drift error needs to be small and the constant drift can be precisely compensated. And as for the accelerometer, high stability of the bias and the scale factor is required [4].

All of the inertial instruments must be precisely tested and calibrated before they can work properly to provide accurate measurement in the INS [5,6]. According to the observed quantities, accelerometer calibration can be divided into two types, device-level calibration and system-level calibration. In device-level calibration, the accelerometer is set to measure the gravity of Earth at different positions by using a dividing head, a centrifuge, or a turntable, and error parameters are identified according to the outputs of the accelerometer [7,8,9]. The accelerometer data measured during the multi-position rotation of the turntable are compared with the theoretical value of the local gravity, and then using the least square method, error parameters can be estimated from their difference. The system-level calibration is based on the propagation of error parameters of inertial instruments to the final navigation errors. Error parameters can be precisely identified with a proper rotation scheme that is designed considering the parameter observability [3].

Jiang et al. [9] proposed a device-level calibration method based on the static-base condition, and the precise attitude and attitude rate provided by a tri-axis precision turntable, where the calibration accuracy can reach ${1\times 10}^{-6}\mathrm{g}$ order of magnitude. However, in practical applications, it is difficult to provide an ideal laboratory environment for the implementation of this method. Won et al. [10] and Ye et al. [11] calibrated the accelerometer on the basis that the amplitude of the output vector of the tri-axis accelerometer equals 1 g. Studies [12,13,14,15,16,17,18,19,20,21] presented several system-level calibration methods specific to the microelectromechanical system (MEMS) inertial measuring unit (IMU), but the calibration accuracy needs to be improved. The methods proposed in references [22,23,24,25,26,27,28,29,30] are appropriate for the rotational INSs, and several novel methods are presented for the hybrid INSs in references [31,32]. Gao et al. [22] used position error and velocity error as the observed quantities to estimate the nonlinear error coefficients. They also used velocity error and attitude error to estimate the error coefficients of the gyroscope and the accelerometer [23]. Pan et al. [33] calibrated the nonlinear term in the scale factor based on the change rate of velocity error. Gao et al. [27] and Liu et al. [31] proposed methods to compensate the size effect and the lever arm effect, respectively. As for the turntable rotation strategy, the 6-position [22], 9-position [24,32], 18-position [18,34], and 24-position [35] methods have been proposed. In existing works, calibration methods are rarely designed for the space-stable INSs. In addition, the complementary benefits of device-level and the system-level calibration have not been taken into consideration.

In the space-stable INS, the platform is stabilized in inertial space by the gimbal servo system. With no applied torque, the moment of momentum of the rotor in the free-rotor gyroscope maintains fixed orientation in inertial space, which can ensure the high precision of the system [1]. As there is a stabilized platform in a space-stable INS, self-calibration based on such an own condition becomes the current development trend. This paper attempts to combine the advantages of full-parameter estimation in the device-level calibration and high-precision estimation of partial parameters in the system-level calibration. A three-stage self-calibration technique is proposed for the high-precision accelerometer calibration in practical applications. Firstly, device-level self-calibration is conducted directly in the system. It is divided into two steps considering the large amount of parameters. The first step is coarse calibration which identifies parameters except for the nonlinear terms, and the second step is fine calibration which identifies the nonlinear terms and simultaneously improves the accuracy of parameters identified in the first step. Secondly, the follow-on system-level calibration is carried out on part of the parameters using the errors of specific force and attitude angles to further improve the calibration accuracy. The proposed method is specific to the space-stable INS and has three contributions. First, taking advantage of the high-precision angulation ability of the gimbaled platform, the two-stage device-level calibration can be directly implemented in the system. Second, based on the error propagation of the residual accelerometer model error to navigation errors, the system-level calibration uses specific force error and attitude error to form the observation model. Third, the proposed three-stage calibration method merges the complementary advantages of both the device-level calibration and the system-level calibration to achieve full-parameter calibration of the accelerometer measurement model as well as high-precision estimation of part of these parameters. This can help further improve the performance of the system.

The rest of this paper is arranged as follows: In Section 2, the accelerometer measurement model is introduced. The device-level coarse calibration and device-level fine calibration are presented in Section 3 and Section 4, respectively. The system-level calibration method based on the errors of specific force and attitude angles is proposed in Section 5. Simulation and experiment results are discussed in Section 6. And finally conclusions are drawn in Section 7.

2. Accelerometer Measurement Model

For the convenience of following discussion, several coordinate frames are defined first, and the superscript of a vector will indicate which frame it is in.

$O_n{x}_n{y}_n{z}_n$—the local geographic navigation frame (n-frame), where the $x_n$ axis and the $y_n$ axis are in the local-level plane pointing north and east, respectively, and the $z_n$ axis is angled vertically downward. The gravity vector g has the same orientation with the $z_n$ axis in this coordinate system.

$O_b{x}_b{y}_b{z}_b$—the body-fixed coordinate frame (b-frame). The origin $O_b$ is usually located at the center of the mass of the vehicle. The $x_b$ axis points forward along the longitudinal axis, the $z_b$ axis points to the bottom, and the $y_b$ axis points rightward to form a right-handed coordinate system.

$O_p{x}_p{y}_p{z}_p$—the platform-fixed coordinate system (p-frame). The $z_p$ axis is along the axis of revolution of the platform. And ideally, $x_p$ and $y_p$ are along the revolution axes of the inner gimbal and the middle gimbal, respectively.

$O_a{x}_a{y}_a{z}_a$—the accelerometer coordinate system (a-frame), which is fixed on the tri-axis accelerometer. The $x_a$, $y_a$ and $z_a$ axes are along the three sensitive axes, respectively. This coordinate system is non-orthogonal due to the installation error.

For a static system, the specific force measured in the n-frame should be:

$${\mathit{f}}^n={\left[\begin{array}{ccc}0& 0& -g\end{array}\right]}^{\mathrm{T}},$$

(1)

where g is the scalar gravity. Transforming it to the p-frame yields:

$${\mathit{f}}^p={\mathit{C}}_b^p{\mathit{C}}_n^b{\mathit{f}}^n,$$

(2)

where the direction cosine matrix from the n-frame to the b-frame is:

$${\mathit{C}}_n^b=\left[\begin{array}{ccc}\cos \psi \cos \theta & \sin \psi \cos \theta & -\sin \theta \\ -\sin \psi \cos \varphi +\cos \psi \sin \theta \sin \varphi & \cos \psi \cos \varphi +\sin \psi \sin \theta \sin \varphi & \cos \theta \sin \varphi \\ \sin \psi \sin \varphi +\cos \psi \sin \theta \cos \varphi & -\cos \psi \sin \varphi +\sin \psi \sin \theta \cos \varphi & \cos \theta \cos \varphi \end{array}\right],$$

(3)

and the direction cosine matrix from the b-frame to the p-frame is:

$${\mathit{C}}_b^p=\left[\begin{array}{ccc}-\cos S_t\sin h\cos q-\sin S_t\sin q& \cos S_t\sin h\sin q-\sin S_t\cos q& -\cos S_t\cos h\\ -\sin S_t\sin h\cos q+\cos S_t\sin q& \sin S_t\sin h\sin q+\cos S_t\cos q& -\sin S_t\cos h\\ \cos h\cos q& -\cos h\sin q& -\sin h\end{array}\right].$$

(4)

In Equation (3), $\psi ,\theta ,\mathrm{and}\varphi$ are yaw, pitch, and roll, respectively. In Equation (4), $S_t$ is the rotation angle of the platform’s revolution axis relative to the revolution axis of the inner gimbal, h is the rotation angle of the revolution axis of the middle gimbal relative to that of the outer gimbal, and q is the rotation angle of the revolution axis of the outer gimbal relative to the base. The rotation angle of the revolution axis of the inner gimbal relative to that of the middle gimbal is zero, for the outer gimbal is tracking the rotation angle of the inner gimbal, and for the revolution axes of the platform, the inner gimbal and the middle gimbal are orthogonal.

The geometric relationship between the orthogonal p-frame and the non-orthogonal a-frame is shown in Figure 1. It is assumed that installation error angles are all small quantities. ${\beta }_x,{\beta }_y\mathrm{and}{\beta }_z$ denote the three installation error angles in the x-direction, y-direction, and z-direction, respectively. Each of these three small angles can be decomposed into two rotation angles around the axes of the p-frame, namely ${\beta }_{xy}\mathrm{and}{\beta }_{xz}$, ${\beta }_{yz}\mathrm{and}{\beta }_{yx},$ and ${\beta }_{zx}\mathrm{and}{\beta }_{zy}$.

According to the geometric relationship shown in Figure 1, the direction cosine matrix from the p-frame to the a-frame can be calculated as:

$${\mathit{C}}_p^a=\left[\begin{array}{ccc}1& {\beta }_{xz}& -{\beta }_{xy}\\ -{\beta }_{yz}& 1& {\beta }_{yx}\\ {\beta }_{zy}& -{\beta }_{zx}& 1\end{array}\right].$$

(5)

Thus the specific force measured by the tri-axis accelerometer satisfies an equation, as:

$${\mathit{f}}^a=\mathit{S}{\mathit{F}}^{-1}{\mathit{C}}_p^a{\mathit{C}}_b^p{\mathit{C}}_n^b{\mathit{f}}^n+\mathsf{\delta }\mathit{b}+{\mathit{K}}_{nl}{\left({\mathit{f}}^a-\mathsf{\delta }\mathit{b}\right)}^2,$$

(6)

where SF is the matrix of scale factors, $\mathsf{\delta }\mathit{b}$ is the bias vector, and ${\mathit{K}}_{nl}$ is the matrix of nonlinear coefficients of the scale factors. In this paper, a three-stage calibration technique is proposed, to estimate and compensate all the parameters shown in Equation (6). The first two stages will lay a foundation of full-parameter and full-value estimations of model coefficients. The third stage will then correct part of the coefficients and further improve the calibration accuracy. Although the final calibration accuracy will mainly depend on the novel system-level calibration method, it must work on the basis of the device-level calibration. Therefore, the first two stages of the device-level calibration designed for the space-stable INS will also be introduced to form the integrated three-stage calibration technique.

3. Device-Level Coarse Calibration

For the tri-axis accelerometer, there are totally 18 error coefficients that need to be calibrated: three scale factors, ${SF}_x,{SF}_y\mathrm{and}{SF}_z$; three biases, ${\mathsf{\delta }b}_x,{\mathsf{\delta }b}_y\mathrm{and}{\mathsf{\delta }b}_z$; three nonlinear coefficients of the scale factors, $K_{nlx},K_{nly}\mathrm{and}{K}_{nlz}$; six installation error angles ${\beta }_{xy},{\beta }_{xz}$, ${\beta }_{yz},{\beta }_{yx},{\beta }_{zx}\mathrm{and}{\beta }_{zy}$; and three zero-point offsets of the gimbals’ angular sensors, ${\theta }_{10}$ of the platform axis, ${\theta }_{20}$ of the inner gimbal axis, and ${\theta }_{30}$ of the middle gimbal axis. The platform needs to be rotated to 18 positions.

Table 1. The order of rotation for the 18-position calibration, and the corresponding rotation angles of the gimbals and components of the p-frame gravity vector.

The Order of Rotation	Rotation Angles of Gimbals			Components of the p-Frame Gravity Vector
	${\mathit{\theta }}_1$	${\mathit{\theta }}_2$	${\mathit{\theta }}_3$	${\mathit{g}}_{\mathit{x}}^{\mathit{p}}$	${\mathit{g}}_{\mathit{y}}^{\mathit{p}}$	${\mathit{g}}_{\mathit{z}}^{\mathit{p}}$
1	${\theta }_{10}$	${\theta }_{20}$	${\theta }_{30}$	${-g\theta }_{30}$	${g\theta }_{20}$	$g$
2	$\downarrow$	${\theta }_{20}+\frac{\pi }{4}$	$\downarrow$	$(\frac{\sqrt{2}}{2}{\theta }_{10}{-\theta }_{30})g$	$\frac{\sqrt{2}}{2}g$	$\frac{\sqrt{2}}{2}g$
3	$\downarrow$	${\theta }_{20}+\frac{\pi }{2}$	$\downarrow$	${g(\theta }_{10}{-\theta }_{30})$	$g$	${-g\theta }_{20}$
4	${\theta }_{10}+\frac{\pi }{4}$	$\downarrow$	$\downarrow$	$\frac{\sqrt{2}}{2}g$	$\frac{\sqrt{2}}{2}g$	${-g\theta }_{20}$
5	${\theta }_{10}+\frac{\pi }{2}$	$\downarrow$	$\downarrow$	$g$	${g(\theta }_{30}{-\theta }_{10})$	${-g\theta }_{20}$
6	$\downarrow$	${\theta }_{20}+\frac{\pi }{4}$	$\downarrow$	$\frac{\sqrt{2}}{2}g$	${(\theta }_{30}-\frac{\sqrt{2}}{2}{\theta }_{10})g$	$\frac{\sqrt{2}}{2}g$
7	$\downarrow$	${\theta }_{20}+\pi$	$\downarrow$	${-g\theta }_{20}$	${g\theta }_{30}$	$-g$
8	$\downarrow$	$\downarrow$	${\theta }_{30}-\frac{\pi }{4}$	$-\frac{\sqrt{2}}{2}{g(\theta }_{10}{+\theta }_{20})$	$-\frac{\sqrt{2}}{2}g$	$-\frac{\sqrt{2}}{2}g$
9	$\downarrow$	$\downarrow$	${\theta }_{30}-\frac{\pi }{2}$	${-g\theta }_{10}$	$-g$	${-g\theta }_{30}$
10	${\theta }_{10}+\frac{3\pi }{4}$	$\downarrow$	$\downarrow$	$-\frac{\sqrt{2}}{2}g$	$-\frac{\sqrt{2}}{2}g$	${-g\theta }_{30}$
11	${\theta }_{10}+\pi$	$\downarrow$	$\downarrow$	$-g$	${g\theta }_{10}$	${-g\theta }_{30}$
12	$\downarrow$	$\downarrow$	${\theta }_{30}-\frac{\pi }{4}$	$-\frac{\sqrt{2}}{2}g$	$\frac{\sqrt{2}}{2}{g(\theta }_{10}{+\theta }_{20})$	$-\frac{\sqrt{2}}{2}g$
13	$\downarrow$	${\theta }_{20}+\frac{5\pi }{4}$	${\theta }_{30}$	$(\frac{\sqrt{2}}{2}{\theta }_{10}{+\theta }_{30})g$	$\frac{\sqrt{2}}{2}g$	$-\frac{\sqrt{2}}{2}g$
14	${\theta }_{10}+\frac{3\pi }{4}$	${\theta }_{20}+\frac{3\pi }{2}$	$\downarrow$	$-\frac{\sqrt{2}}{2}g$	$\frac{\sqrt{2}}{2}g$	$\frac{\sqrt{2}}{2}{g\theta }_{20}$
15	${\theta }_{10}+\frac{\pi }{2}$	${\theta }_{20}+\frac{7\pi }{4}$	$\downarrow$	$-\frac{\sqrt{2}}{2}g$	$\frac{\sqrt{2}}{2}{g(\theta }_{30}{+\theta }_{20})$	$\frac{\sqrt{2}}{2}g$
16	$\downarrow$	${\theta }_{20}$	${\theta }_{30}-\frac{\pi }{4}$	$\frac{\sqrt{2}}{2}{g(-\theta }_{10}{+\theta }_{20})$	$-\frac{\sqrt{2}}{2}g$	$\frac{\sqrt{2}}{2}g$
17	${\theta }_{10}+\frac{\pi }{4}$	$\downarrow$	${\theta }_{30}-\frac{\pi }{2}$	$\frac{\sqrt{2}}{2}g$	$-\frac{\sqrt{2}}{2}g$	${g\theta }_{30}$
18	${\theta }_{10}$	$\downarrow$	${\theta }_{30}-\frac{3\pi }{4}$	$\frac{\sqrt{2}}{2}g$	$-\frac{\sqrt{2}}{2}{g(\theta }_{10}{+\theta }_{20})$	$-\frac{\sqrt{2}}{2}g$

lists the order of rotation, and the corresponding rotation angles of the gimbals and components of the p-frame gravity vector [3], where the symbol $\downarrow$ denotes that the angle remains unchanged. Figure 2 shows the directions of the accelerometer’s sensitive axes in detail, where N, W, and U denote the local northern, western and up directions, respectively, and x, y, and z denote the orientations of the tri-axis accelerometer.

During coarse calibration, the nonlinearity of scale factors is ignored. Thus the accelerometer measurement equation shown in Equation (6) can be simplified as:

$${\mathit{f}}^a=\mathit{S}{\mathit{F}}^{-1}{\mathit{C}}_p^a{\mathit{g}}^p+\mathsf{\delta }\mathit{b},$$

(7)

where ${\mathit{C}}_p^a{\mathit{g}}^p$ is the gravity vector sensed by the tri-axis accelerometer.

Each position takes 20 s for the rapid measurement. Based on the outputs of the tri-axis accelerometer at different positions and Table 1, 18 equations can be obtained according to Equation (7) to preliminarily estimate 15 parameters except the nonlinear coefficients. Scale factors and biases can be coarsely estimated by solving the equation set corresponding to the six rotation positions (No. 1, 3, 5, 7, 9, and 11) where each sensitive axis of the tri-axis accelerometer points upward or downward along the $z_n$ axis in turn. The other equations corresponding to the rest of the rotation positions (No. 2, 4, 6, 8, 10, and 12~18), where one of the sensitive axes is parallel to the $x_n$ axis or the $y_n$ axis, and other two sensitive axes have included angles of 45° with the vertical direction, are used to estimate zero-point offsets of the angular sensors and installation error angles.

Estimated values of zero-point offsets, ${\hat{\theta }}_{20}$ and ${\hat{\theta }}_{30}$, are used to coarsely correct the orientation of the p-frame, and to make the directions of the middle gimbal axis and the inner gimbal axis approach the directions of axes of the p-frame. Concretely, the inner gimbal axis and the middle gimbal axis are inversely rotated ${\hat{\theta }}_{20}$ and ${\hat{\theta }}_{30}$, respectively, and then coarse calibration is re-executed. This procedure will repeat until the re-estimated ${\hat{\theta }}_{20}$ and ${\hat{\theta }}_{30}$ decrease to negligibly small quantities.

4. Device-Level Fine Calibration

Fine calibration will be executed on the basis of coarse calibration. The orientation of the platform has been corrected and new data will be measured in more accurate rotation positions. In order to improve the data collection accuracy, the measurement time of each position is extended to 10 min. Moreover, to improve the data processing accuracy, the estimation error of coarse calibration is estimated instead of the full values of error parameters.

According to Equation (6), the ideal specific force in the p-frame satisfies:

$${\mathit{f}}^p={\mathit{C}}_a^p\mathit{S}F\left({\mathit{f}}^a-\mathsf{\delta }\mathit{b}-{\mathit{K}}_{nl}{\left({\mathit{f}}^a-\mathsf{\delta }\mathit{b}\right)}^2\right).$$

(8)

If the result of coarse calibration and outputs of the tri-axis accelerometer during fine calibration ($f^a$) have been obtained, the coarse estimation of $f^p$ can be calculated as:

$${\widehat{\mathit{f}}}^p={\widehat{\mathit{C}}}_a^p\widehat{\mathit{S}\mathit{F}}\left({\mathit{f}}^a-\widehat{\mathsf{\delta }\mathit{b}}\right),$$

(9)

where “hatted” symbols denote the estimated values.

In order to reduce the requirement for position control accuracy of the platform during fine calibration, the observed variable is not the direct difference between the coarsely estimated value and the ideal value of the specific force at one certain position, but the difference between the squared amplitudes of them. Subtracting the squared amplitude of Equation (9) from that of Equation (8), and then taking half of the result yields the observed variable, as:

$$\begin{array}{cc}\hfill \frac{1}{2}\left({\left({\mathit{f}}^p\right)}^{\mathrm{T}}{\mathit{f}}^p-{\left({\widehat{\mathit{f}}}^p\right)}^{\mathrm{T}}{\widehat{\mathit{f}}}^p\right)& ={\left(\widehat{\mathit{S}\mathit{F}}\left({\mathit{f}}^a-\mathsf{\delta }\mathit{b}\right)\right)}^2{\widehat{\mathit{S}\mathit{F}}}^{-1}\Delta \mathit{S}\mathit{F}\hfill \\ & +\frac{1}{2}{\left(\widehat{\mathit{S}\mathit{F}}\left({\mathit{f}}^a-\mathsf{\delta }\mathit{b}\right)\right)}^{\mathrm{T}}\left({\left(\Delta {\mathit{C}}_p^a\right)}^{\mathrm{T}}+\Delta {\mathit{C}}_p^a\right)\left(\widehat{\mathit{S}\mathit{F}}\left({\mathit{f}}^a-\mathsf{\delta }\mathit{b}\right)\right)\hfill \\ & -{\left(\widehat{\mathit{S}\mathit{F}}\left({\mathit{f}}^a-\mathsf{\delta }\mathit{b}\right)\right)}^{\mathrm{T}}\widehat{\mathit{S}\mathit{F}}\widetilde{\mathsf{\delta }\mathit{b}}-{\left(\widehat{\mathit{S}\mathit{F}}\left({\mathit{f}}^a-\mathsf{\delta }\mathit{b}\right)\right)}^{\mathrm{T}}\widehat{\mathit{S}\mathit{F}}{\mathit{K}}_{\mathrm{nl}}\left({\mathit{f}}^a-\mathsf{\delta }\mathit{b}\right)\left({\mathit{f}}^a-\mathsf{\delta }\mathit{b}\right)\hfill \end{array},$$

(10)

where:

$$\frac{1}{2}\left({\left({\mathit{f}}^p\right)}^{\mathrm{T}}{\mathit{f}}^p-{\left({\widehat{\mathit{f}}}^p\right)}^{\mathrm{T}}{\widehat{\mathit{f}}}^p\right)=\frac{1}{2}\left({\left(f^p\right)}^2-{\left({\widehat{f}}^p\right)}^2\right)\approx f^p\left(f^p-{\widehat{f}}^p\right),$$

(11)

$$\Delta \mathit{S}\mathit{F}=\mathit{S}\mathit{F}-\widehat{\mathit{S}\mathit{F}}=\left[\begin{array}{ccc}\Delta S{F}_x& 0& 0\\ 0& \Delta S{F}_y& 0\\ 0& 0& \Delta S{F}_z\end{array}\right],$$

(12)

$$\widetilde{\delta \mathit{b}}=\delta \mathit{b}-\widehat{\delta \mathit{b}}=\left[\begin{array}{ccc}\widetilde{\delta b_x}& 0& 0\\ 0& \widetilde{\delta b_y}& 0\\ 0& 0& \widetilde{\delta b_z}\end{array}\right],$$

(13)

$$\Delta {\mathit{C}}_p^a=\left[\begin{array}{ccc}0& -\Delta {\beta }_{xz}& \Delta {\beta }_{xy}\\ \Delta {\beta }_{yz}& 0& -\Delta {\beta }_{yx}\\ -\Delta {\beta }_{zy}& \Delta {\beta }_{zx}& 0\end{array}\right],\mathrm{and}$$

(14)

$${\widehat{\mathit{f}}}^a=\widehat{\mathit{S}\mathit{F}}\left({\mathit{f}}^a-\delta \mathit{b}\right)=\left(\begin{array}{c}{\widehat{f}}_x^a\\ {\widehat{f}}_y^a\\ {\widehat{f}}_z^a\end{array}\right).$$

(15)

According to Equations (10)–(15), the observation equation at position j can be written as:

$$\begin{array}{c}[\begin{array}{ccccccccc}{\left({\widehat{f}}_{xj}^a\right)}^2& {\left({\widehat{f}}_{yj}^a\right)}^2& {\left({\widehat{f}}_{zj}^a\right)}^2& -{\widehat{f}}_{xj}^a& -{\widehat{f}}_{yj}^a& -{\widehat{f}}_{zj}^a& -{\left({\widehat{f}}_{xj}^a\right)}^3& -{\left({\widehat{f}}_{yj}^a\right)}^3& -{\left({\widehat{f}}_{zj}^a\right)}^3\end{array}\hfill \\ \begin{array}{ccc}{\widehat{f}}_{xj}^a{\widehat{f}}_{yj}^a& {\widehat{f}}_{xj}^a{\widehat{f}}_{zj}^a& {\widehat{f}}_{zj}^a{\widehat{f}}_{yj}^a\end{array}]\mathrm{d}\mathit{x}=f^p\left(f^p-{\widehat{f}}_j^p\right)\end{array},$$

(16)

where the unknown error vector is:

$$\begin{array}{cc}\hfill d\mathit{x}=& [\begin{array}{cccccc}\frac{\widehat{\Delta S{F}_x}}{\widehat{S{F}_x}}& \frac{\widehat{\Delta S{F}_y}}{\widehat{S{F}_y}}& \frac{\widehat{\Delta S{F}_z}}{\widehat{S{F}_z}}& \widehat{S{F}_x}{\widetilde{\Delta b}}_x& \widehat{S{F}_y}{\widetilde{\Delta b}}_y& \widehat{S{F}_z}{\widetilde{\Delta b}}_z\end{array}\hfill \\ & \begin{array}{cccccc}\frac{K_{\mathrm{nlx}}}{\widehat{S{F}_x}}& \frac{K_{\mathrm{nly}}}{\widehat{S{F}_y}}& \frac{K_{\mathrm{nlz}}}{\widehat{S{F}_z}}& \Delta {\beta }_{yz}-\Delta {\beta }_{xz}& \Delta {\beta }_{xy}-\Delta {\beta }_{zy}& \Delta {\beta }_{zx}-\Delta {\beta }_{yx}\end{array}{]}^{\mathrm{T}}\hfill \end{array}.$$

(17)

Based on the data collected at the 18 positions, an $18\times 12$-dimension observation equation can be formed based on Equation (16), and the unknown variables shown in Equation (17) can be solved by the least squares method. Finally, using scale factors to convert the solved $\mathrm{d}\mathit{x}$ and adding the result to the coarse estimation will yield a full-value fine estimation of 12 error coefficients. It should be noted that ${\beta }_{xz}\mathrm{and}{\beta }_{yz}$, ${\beta }_{zy}\mathrm{and}{\beta }_{xy},\mathrm{and}{\beta }_{yx}\mathrm{and}{\beta }_{zx}$ are linearly correlated in pairs and cannot be decomposed from each other. Therefore the fine calibration can only improve the accuracy of ${\mathit{f}}^p$, but it cannot calibrate the rotation characteristics of the a-frame.

5. System-Level Calibration Based on Specific Force Error and Attitude Error

Errors of inertial instruments will induce the navigational errors of INSs, including the positional error, the velocity error and the attitude error. Inversely, the independently observed navigation errors such as the specific force error, the attitude error and the position error can be used to calibrate the main error parameters in the system. As the position error is affected by both the gyroscope error and accelerometer error, it is generally difficult to directly extract the accelerometer-induced term from the position error. On the other hand, the specific force error and the horizontal attitude error are mainly caused by the accelerometer error. Therefore this paper employs the specific force and horizontal attitude angles to correct part of the accelerometer’s error coefficients. The prerequisite for accurate observation in such a system-level calibration is that the system should be installed on a static and horizontal base.

5.1. Calibration Using Specific Force Error

Under the ideal condition, that there is no installation error between the base and the stabilized platform, and that the measurement of gimbals’ rotation angle is accurate, there is an approximate relationship:

$${\mathit{C}}_b^p{\mathit{C}}_n^b\approx \left[\begin{array}{ccc}-\cos S_t\sin L& -\sin S_t& -\cos S_t\cos L\\ -\sin S_t\sin L& \cos S_t& -\sin S_t\cos L\\ \cos L& 0& -\sin L\end{array}\right],$$

(18)

where L denotes the latitude. On the static and horizontal base, the local scalar gravity g can be accurately calculated. Then, according to Equations (1) and (18), the ideal specific force vector in the p-frame can be written as:

$$\frac{{\mathit{f}}^p}{g}=\left(\begin{array}{c}{f}_x^P/g\\ f_y^P/g\\ f_z^P/g\end{array}\right)={\mathit{C}}_b^p{\mathit{C}}_n^b\left(\begin{array}{c}0\\ 0\\ -1\end{array}\right)=\left(\begin{array}{c}\cos S_t\cos L\\ \sin S_t\cos L\\ \sin L\end{array}\right).$$

(19)

Similar to the device-level calibration method presented in the previous section, the difference of the squared amplitudes of the measured and the ideal specific force vectors is chosen as the observed variable, as:

$$z=\frac{1}{2}\left[{\left|\frac{{\widehat{\mathit{f}}}^p}{g}\right|}^2-{\left|\frac{{\mathit{f}}^p}{g}\right|}^2\right]\approx {\left(\frac{{\mathit{f}}^p}{g}\right)}^{\mathrm{T}}\frac{\Delta {\mathit{f}}^p}{g}.$$

(20)

A 5-dimensional state vector is constructed as $\mathit{x}={\left[\begin{array}{ccccc}{x}_1& x_2& x_3& x_4& x_5\end{array}\right]}^{\mathrm{T}}$, where:

$$x_1=\left(\frac{\Delta S{F}_x}{2S{F}_x}+\frac{\Delta S{F}_y}{2S{F}_y}\right){\cos }^{2}L-\left(\frac{{\widetilde{\delta b}}_z}{g}-\frac{\Delta S{F}_z}{S{F}_z}\sin L\right)\sin L,$$

(21)

$$x_2=\frac{{\widetilde{\delta b}}_x}{g}-\Delta {\beta }_{xy}\sin L+\Delta {\beta }_{zy}\sin L,$$

(22)

$$x_3=\frac{{\widetilde{\delta b}}_y}{g}+\Delta {\beta }_{yx}\sin L-\Delta {\beta }_{zx}\sin L,$$

(23)

$$x_4=\frac{\Delta S{F}_x}{2S{F}_x}-\frac{\Delta S{F}_y}{2S{F}_y},\mathrm{and}$$

(24)

$$x_5=\frac{\Delta {\beta }_{yz}-\Delta {\beta }_{xz}}{2}.$$

(25)

Thus the observation equation is:

$$z=\mathit{H}x+v.$$

(26)

Considering Equation (19), the observation matrix is:

$$\mathit{H}=\left[\begin{array}{ccccc}\mathit{1}& -f_x^p/g& -f_y^p/g& \left({\left(f_x^p\right)}^2-{\left(f_y^p\right)}^2\right)/g^2& 2f_x^p{f}_y^p/g^2\end{array}\right].$$

(27)

The five components of $\mathit{x}$ are actually the mean value, the sine and the cosine fundamental waves, and the sine and the cosine second harmonics of the Earth’s rotation rate in the specific force error. Assuming that the observation noise is a zero-mean Gaussian white noise, $\mathit{x}$ can be estimated by a Kalman filter.

5.2. Calibration Using Attitude Error

On the static and horizontal base, the attitude error of the INS with both vertical damping and horizontal velocity damping mainly consists of a constant value and a 24 h long period component, which is induced by the gyroscope error, the accelerometer error, and the rotation angle errors of platform’s gimbals. Using the space-stable mechanization, the azimuth error will accumulate over time due to the gyroscope error, while the horizontal attitude errors are mainly caused by the accelerometer error and have no such accumulation. Therefore, the horizontal attitude errors are chosen as the observed quantities for accelerometer calibration.

Neglecting the short-period misalignment angle of gyro case rotation, the attitude error on the static and horizontal base has an expression as:

$$\left(\begin{array}{c}\mathsf{\delta }\varphi \\ \mathsf{\delta }\theta \\ \mathsf{\delta }\psi \end{array}\right)=\left(\begin{array}{c}\frac{\Delta f_y^e}{g}\\ \frac{\Delta f_x^e\sin L-\Delta f_z^e\cos L}{g}\\ -\frac{\Delta S_1}{\cos L}-\frac{\Delta f_y^e}{g}\tan L\end{array}\right)-\left(\begin{array}{c}\mathsf{\delta }{S}_t\cos L+\mathsf{\delta }{\kappa }_x^b\cos \psi -\mathsf{\delta }{\kappa }_y^b\sin \psi \\ -\mathsf{\delta }h+\mathsf{\delta }{\kappa }_x^b\sin \psi +\mathsf{\delta }{\kappa }_y^b\cos \psi \\ \mathsf{\delta }q-\mathsf{\delta }{S}_t\sin L+\mathsf{\delta }{\kappa }_z^b\end{array}\right),$$

(28)

where ${\mathsf{\delta }\kappa }_x^b$, ${\mathsf{\delta }\kappa }_y^b$, and ${\mathsf{\delta }\kappa }_z^b$ are the transverse tilt error, the longitudinal tilt error, and the azimuth error of the base, respectively, and the superscript e denotes the earth-fixed frame (e-frame).

Substituting the e-frame specific force error into Equation (28) yields the expressions of pitch error and roll error. Both of them consist of the mean value, the fundamental wave, and the second harmonic of the Earth’s rotation rate. The five Fourier coefficients of the pitch error are:

$$\overline{\mathsf{\delta }\theta }=-\left(-\mathsf{\delta }h+\mathsf{\delta }{\kappa }_x^b\sin \psi +\mathsf{\delta }{\kappa }_y^b\cos \psi \right)+\left[\frac{{\nabla }_z}{g}-\left(\frac{\Delta S{F}_z}{S{F}_z}-\frac{\Delta S{F}_x}{2S{F}_x}-\frac{\Delta S{F}_y}{2S{F}_y}\right)\sin L\right]\cos L,$$

(29)

$$\widetilde{\mathsf{\delta }{\theta }_{1c}}=-\left(\frac{{\nabla }_x}{g}\sin L-\Delta {\beta }_{xy}{\sin }^{2}L-\Delta {\beta }_{zy}{\cos }^{2}L\right),$$

(30)

$$\widetilde{\mathsf{\delta }{\theta }_{1s}}=-\left(\frac{{\nabla }_y}{g}\sin L+\Delta {\beta }_{yx}{\sin }^{2}L+\Delta {\beta }_{zx}{\cos }^{2}L\right),$$

(31)

$$\widetilde{\mathsf{\delta }{\theta }_{2c}}=\left(\frac{\Delta S{F}_x}{2S{F}_x}-\frac{\Delta S{F}_y}{2S{F}_y}\right)\sin L\cos L,\mathrm{and}$$

(32)

$$\widetilde{\mathsf{\delta }{\theta }_{2s}}=\left(\frac{\Delta {\beta }_{yz}}{2}-\frac{\Delta {\beta }_{xz}}{2}\right)\sin L\cos L.$$

(33)

The other five Fourier coefficients of the roll error are:

$$\overline{\mathsf{\delta }\varphi }=-\left(\mathsf{\delta }{S}_t\cos L+\mathsf{\delta }{\kappa }_x^b\cos \psi -\mathsf{\delta }{\kappa }_y^b\sin \psi \right)+\left(\frac{\Delta {\beta }_{yz}}{2}+\frac{\Delta {\beta }_{xz}}{2}\right)\cos L,$$

(34)

$$\widetilde{\mathsf{\delta }{\varphi }_{1s}}=\frac{{\nabla }_x}{g}-\Delta {\beta }_{xy}\sin L,$$

(35)

$$\widetilde{\mathsf{\delta }{\varphi }_{1c}}=-\frac{{\nabla }_y}{g}-\Delta {\beta }_{yx}\sin L,$$

(36)

$$\widetilde{\mathsf{\delta }{\varphi }_{2s}}=-\left(\frac{\Delta S{F}_x}{2S{F}_x}-\frac{\Delta S{F}_y}{2S{F}_y}\right)\cos L,\mathrm{and}$$

(37)

$$\widetilde{\mathsf{\delta }{\varphi }_{2c}}=\left(\frac{\Delta {\beta }_{yz}}{2}-\frac{\Delta {\beta }_{xz}}{2}\right)\cos L.$$

(38)

The subscripts indicate which component each expression is (1 for the fundamental waves, 2 for the second harmonics, s for the sine functions, and c for the cosine functions). After the error curves of pitch and roll are obtained, these coefficients can be extracted by Fourier expansion.

5.3. Parameter Separation Algorithm for Combined Error Coefficients of Accerelometer

In Section 5.1 and Section 5.2, 15 combined coefficients (Equations (21)–(25), (29)–(33), and (34)–(38)) are estimated by static Kalman filter and Fourier expansion. Twelve of them are selected and divided into three groups to extract the residual error parameters of the tri-axis accelerometer. Since only the fundamental waves and second harmonics of pitch and roll errors are concerned, the requirement for the horizontality of the base need not be very strict.

The first group consists of three equations, which are the coefficients of the sine and cosine second harmonics of Earth’s rotation rate about $\frac{{\Delta SF}_x}{{2SF}_x}-\frac{{\Delta SF}_y}{{2SF}_y}$. They are the linearly correlated Equations (24), (32), and (37). Thus $\frac{{\Delta SF}_x}{{2SF}_x}-\frac{{\Delta SF}_y}{{2SF}_y}$ can be solved using their equally weighted average, as:

$$\widehat{\frac{\Delta S{F}_x}{2S{F}_x}}-\widehat{\frac{\Delta S{F}_y}{2S{F}_y}}=\frac{1}{3}\left(x_4-\frac{{\varphi }_{\mathrm{N}}{}_{2\mathrm{s}}}{\cos L}+\frac{{\varphi }_{\mathrm{E}}{}_{2\mathrm{c}}}{\cos L\sin L}\right).$$

(39)

Similarly, the second group consists of Equations (25), (33), and (38), and $\frac{{\Delta \beta }_{yz}{-\Delta \beta }_{xz}}{2}$ can be solved as:

$$\frac{\widehat{\Delta {\beta }_{yz}}-\widehat{\Delta {\beta }_{xz}}}{2}=\frac{1}{3}\left(x_5+\frac{{\varphi }_{\mathrm{N}}{}_{2\mathrm{c}}}{\cos L}+\frac{{\varphi }_{\mathrm{E}}{}_{2\mathrm{s}}}{\cos L\sin L}\right).$$

(40)

Another six parameters ($\frac{{𝛻}_x}{g}$, $\frac{{𝛻}_y}{g}$, ${\Delta \beta }_{yx}$, ${\Delta \beta }_{zx}$, ${\Delta \beta }_{xy}$, ${\Delta \beta }_{zy}$) can be solved using the least square method based on the remaining six combined coefficients, according to:

$$\left[\begin{array}{cccccc}1& 0& 0& 0& -\sin L& 0\\ 0& -1& -\sin L& 0& 0& 0\\ -\sin L& 0& 0& 0& {\sin }^{2}L& {\cos }^{2}L\\ 0& -\sin L& -{\sin }^{2}L& -{\cos }^{2}L& 0& 0\\ 1& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0\end{array}\right]\left[\begin{array}{c}{𝛻}_x/g\\ {𝛻}_y/g\\ \Delta {\beta }_{yx}\\ \Delta {\beta }_{zx}\\ \Delta {\beta }_{xy}\\ \Delta {\beta }_{zy}\end{array}\right]=\left[\begin{array}{c}\widetilde{\mathsf{\delta }{\varphi }_{1s}}\\ \widetilde{\mathsf{\delta }{\varphi }_{1\mathrm{c}}}\\ \widetilde{\mathsf{\delta }{\theta }_{1\mathrm{c}}}\\ \widetilde{\mathsf{\delta }{\theta }_{1\mathrm{s}}}\\ \widehat{x_2}\\ \widehat{x_3}\end{array}\right].$$

(41)

The six parameters in Equation (41) can be directly compensated by substituting them into the accelerometer’s error model. Also, $\frac{{\Delta \beta }_{yz}{-\Delta \beta }_{xz}}{2}$ and $\frac{{\Delta SF}_x}{{2SF}_x}-\frac{{\Delta SF}_y}{{2SF}_y}$ will be indirectly compensated by deducting the error terms about them from the n-frame specific force expression.

6. Simulation and Experiment

6.1. Simulation Result

As the final accuracy of the three-stage calibration mainly depends on the performance of the last stage, the simulation of the system-level calibration is conducted first. In the simulation, the location of the static base is set at 40° N. The non-orthogonal installation error of the tri-axis accelerometer is quite a stable mechanical error, and after the 18-position two-stage device-level calibration, the residuals are generally less than 10″. Thus, the residuals of accelerometer’s error parameters are assigned as

Table 2. Residuals of scale factor errors and biases.

$\mathbf{\Delta }{\mathit{SF}}_{\mathit{x}}/{\mathit{SF}}_{\mathit{x}}$	$\mathbf{\Delta }{\mathit{SF}}_{\mathit{y}}/{\mathit{SF}}_{\mathit{y}}$	$\mathbf{\Delta }{\mathit{SF}}_{\mathit{z}}/{\mathit{SF}}_{\mathit{z}}$	${\mathbf{\nabla }}_{\mathit{x}}$	${\mathbf{\nabla }}_{\mathit{y}}$	${\mathbf{\nabla }}_{\mathit{z}}$
${10}^{-4}$	${2\times 10}^{-4}$	${-10}^{-4}$	${-10}^{-3}{\mathrm{m}/\mathrm{s}}^2$	${5\times 10}^{-4}{\mathrm{m}/\mathrm{s}}^2$	${10}^{-3}{\mathrm{m}/\mathrm{s}}^2$

and

Table 3. Residuals of installation error angles.

${\mathbf{\Delta }\mathit{\beta }}_{\mathit{xz}}$	${\mathbf{\Delta }\mathit{\beta }}_{\mathit{yz}}$	${\mathbf{\Delta }\mathit{\beta }}_{\mathit{zy}}$	${\mathbf{\Delta }\mathit{\beta }}_{\mathit{xy}}$	${\mathbf{\Delta }\mathit{\beta }}_{\mathit{yx}}$	${\mathbf{\Delta }\mathit{\beta }}_{\mathit{zx}}$
5″	10″	−5″	−10″	5″	−5″

. In addition, a zero-mean random noise with an amplitude of ${5\times 10}^{-3}{\mathrm{m}/\mathrm{s}}^2$ is added to each sensitive axis of the tri-axis accelerometer.

Truth values of the system-level combined coefficients are calculated according to the error parameters shown in Table 2 and Table 3, and are listed in

Table 4. Simulation result of system-level calibration of the combined accelerometer error coefficients (unit: rad).

Combined Coefficients	Reference Values	Estimated Values	Errors
$x_1$	$1.1228\times {10}^{-4}$	$1.1256\times {10}^{-4}$	$2.8\times {10}^{-7}$
$x_2$	$-8.645\times {10}^{-5}$	$-8.612\times {10}^{-5}$	$3.3\times {10}^{-7}$
$x_3$	$8.217\times {10}^{-5}$	$8.206\times {10}^{-5}$	$-1.1\times {10}^{-7}$
$x_4$	$-5.000\times {10}^{-5}$	$-4.957\times {10}^{-5}$	$4.3\times {10}^{-7}$
$x_5$	$1.212\times {10}^{-5}$	$1.185\times {10}^{-5}$	$-2.7\times {10}^{-7}$
$\mathsf{\delta }{\varphi }_{1s}$	$-1.3319\times {10}^{-4}$	$-1.3343\times {10}^{-4}$	$-2.4\times {10}^{-7}$
$\mathsf{\delta }{\varphi }_{1c}$	$-3.542\times {10}^{-5}$	$-3.575\times {10}^{-5}$	$-3.3\times {10}^{-7}$
$\mathsf{\delta }{\varphi }_{2s}$	$3.830\times {10}^{-5}$	$3.848\times {10}^{-5}$	$1.8\times {10}^{-7}$
$\mathsf{\delta }{\varphi }_{2c}$	$-9.28\times {10}^{-6}$	$-9.47\times {10}^{-6}$	$-1.9\times {10}^{-7}$
${\mathsf{\delta }\theta }_{1s}$	$-5.521\times {10}^{-5}$	$-5.536\times {10}^{-5}$	$-1.5\times {10}^{-7}$
${\mathsf{\delta }\theta }_{1c}$	$1.3629\times {10}^{-4}$	$1.3695\times {10}^{-4}$	$6.6\times {10}^{-7}$
${\mathsf{\delta }\theta }_{2s}$	$-5.97\times {10}^{-6}$	$-5.76\times {10}^{-6}$	$2.1\times {10}^{-7}$
${\mathsf{\delta }\theta }_{2c}$	$-2.462\times {10}^{-5}$	$-2.484\times {10}^{-5}$	$-2.2\times {10}^{-7}$

as reference values. The combined coefficients are estimated by system-level calibration based on specific force error and attitude error, and are also listed in Table 4. Comparison shows that the estimation accuracy of each combined coefficient is better than ${10}^{-6}\mathrm{rad}$, which is equivalent to an accelerometer calibration accuracy that is better than ${10}^{-6}$ g. It indicates that the method can be used for the accelerometer calibration of high-precision INSs.

6.2. Experiment Results

In order to validate the proposed three-stage accelerometer calibration method, a static-base navigation test has been conducted using a high-precision space-stable INS. The INS is in the navigation grade class and has been introduced in reference [36]. Two dual-axis free-rotor gyros and a tri-axis quartz accelerometer are mounted on a gimbaled platform. According to the outputs of the free-rotor gyros, the platform is servo-controlled to be stabilized in the inertial space. The gyros have an ultra-low drift of less than 0.001°/h. The measurement error of the accelerometer is of a ${1\times 10}^{-5}$ g order of magnitude, and its stability is better than ${1\times 10}^{-5}$ g per year. The static measurement accuracy of gimbals’ rotation angles can reach the level of 1″. The position error of this navigation system can reach the order of magnitude of 1 n mile per day. The data used in the experiment are collected when the system is installed on a static and horizontal base.

Navigation errors before calibration, after the original two-stage 18-position coarse and fine self-calibration, and after the entire three-stage calibration are compared in Figure 3, Figure 4 and Figure 5 and

Table 5. Maximums of navigation errors.

Navigation Errors	Before Calibration	Two-Stage Calibration	Improvement after Two-Stage Calibration	Three-Stage Calibration	Improvement after Three-Stage Calibration
Longitude	0.84	0.58	30.95%	0.33	60.71%
Latitude	0.86	0.60	30.23%	0.26	69.77%
Northern Velocity	0.85	0.70	17.65%	0.56	34.12%
Eastern Velocity	0.87	0.78	10.34%	0.67	22.99%
Roll	0.89	0.50	43.82%	0.38	57.30%
Pitch	0.95	0.58	38.95%	0.43	54.74%
Yaw	0.80	0.54	32.50%	0.44	45.00%

. Since the system employs vertical damping, height error and vertical velocity error can be neglected and are not presented. The output position of the system is compared with the reference position of the base to obtain the position error. The output velocity and horizontal attitude angles are actually the velocity error and horizontal attitude error, for the base is static and horizontal. All experiment results here have been normalized.

The system has employed gyro case rotation techniques to decrease the gyro drift. Although its influence on the platform’s motion has been deducted during navigation calculating, there are still 8-minute residuals in the velocity and attitude error curves. As the cycle of these residuals are much shorter than the test time, they are severely squeezed in the error curves shown in Figure 4 and Figure 5.

Accelerometer parameters will change after the high-precision INS works for a period of time, which will consequently lead to the increase of navigation errors. Therefore, on-site and online calibration is necessary to improve the navigation precision. Latitude error, longitude error, and attitude errors are reduced by over 30% after the two-stage calibration, and by over 45% after the three-stage calibration. Since the velocity error mainly consists of Schuler oscillation, and the error term is related to the gyro case rotation residual, the improvement of velocity accuracy is not conspicuous. The velocity error has decreased by over 10% after the two-stage calibration, and by over 22% after the three-stage calibration.

7. Conclusions

A three-stage tri-axis accelerometer self-calibration technique is proposed with a thorough theoretical and mathematical analysis. Compared with existing calibration methods, the proposed three-stage calibration method has the following features and advantages: (1) Taking full advantage of high-precision angulation of the gimbaled platform in the space-stable INS, the two-stage device-level calibration is implemented directly on the platform of the space-stable INS. Thus the entire three-stage self-calibration can be implemented in the INS after the tri-axis accelerometer has been mounted on the platform. Besides, during the device-level calibration, zero-point offsets of the gimbals’ angular sensors are simultaneously estimated and used to correct the platform’s orientation. This will help improve the system’s navigation performance. (2) Different from the existing calibration methods which usually estimate error parameters of both accelerometer and gyroscope simultaneously, the proposed system-level calibration is designed based on error propagation to choose observed quantities only resulting from accelerometer measurement error. This can get rid of the influence of gyro drifts. (3) Combining device-level calibration and system-level calibration, the proposed method can achieve both full-parameter estimation of the accelerometer measurement model and high-precision correction for part of the parameters, while existing calibration methods normally concentrate on one of the two aspects. Both simulation and experiment results validate the proposed method. The simulation result shows that using the proposed three-stage calibration technique, the estimation accuracy of accelerometer error is better than ${1\times 10}^{-6}$ g. Experiment results show that after implementing the three-stage calibration, during navigation, the accuracy of latitude, longitude, and attitude angles has increased by over 45%, and the accuracy of velocity has increased by over 22%.