Research on Initial Alignment and Self-Calibration of Rotary Strapdown Inertial Navigation Systems

Abstract

The errors of inertial sensors affect the navigation accuracy of the strapdown inertial navigation system (SINS) and are accumulated over time in nature. In order to continuously maintain the high navigation accuracy of vehicles for a long time period, an initial alignment and self-calibration is necessary after the SINS starts. Additionally, the observability analysis is one of the key techniques during the initial alignment and self-calibration process. For marine systems, the observability of inertial sensor errors is extremely low, as their motion states are always slow. Therefore, studying the rotating SINS is urgent. Since traditional analysis methods have their limitations, the global observation analysis method was used in this paper. On the basis of this method, the relationship between the observability and the kinestate of the rotating SINS has been established. After the discussion about the factors that affect the observability in detail, the design principle of the initial alignment and self-calibration rotating scheme, which is appropriate for marine systems, id proposed. With the proposed principle, a novel initial alignment and self-calibration method, named the eight-position rotating scheme, is designed. Simulations and experiments are carried out to verify its performance. The results have shown that compared with other rotating schemes and the static state, the estimated accuracy of the eight-position scheme rotating about axes x and y was the best, and the position error was significantly reduced with this new rotating scheme. The feasibility and effectiveness of the proposed design principle and the rotating scheme were verified.

1. Introduction

In modern marine navigation, the strapdown inertial navigation system (SINS) is widely used, due to its advantages of being more compact and autonomous [1–4]. However, the successive starting error of inertial sensors is one of the most major factors that will affect the system accuracy. Therefore, to maintain high navigation accuracy for a long time period, it is necessary to perform on-line calibration after the SINS starts [5–8]. The core of the SINS on-line calibration is that changing the motion states of the inertial measurement unit (IMU), including the angular motion and the linear motion, can significantly improve the system observability, and then, the optimal estimation of the inertial sensor errors can be obtained [9–11]. For high-speed carriers, such as aircraft, it is obviously easy to change the motion states. However, for sailing ships on the sea surface, during the mooring, the maneuver can be overlooked by approximation; even during sailing, the motor is still unobvious, as the acceleration, de-acceleration and turning of ships are all extremely slow. Therefore, for vessels, it is not easy to change the spatial position of the IMU. In order to deal with the above-mentioned problem, the dual-axial rotating SINS emerged. In the rotating SINS, according to specific rotary schemes, the IMU spatial position can be changed by utilizing a rotating mechanism that has two degrees of freedom, increasing the observability of the inertial sensor errors dramatically [10,12,13].

As we all know, the self-calibration of the rotary SINS can be seen as a system optimal estimation using the Kalman filter, so the observability analysis of the SINS is required before the filter design [14–16]. The traditional observability analysis methods, such as the piece-wise constant system (PWCS) and the singular value decomposition (SVD), are the most commonly used [17–19]. A time-varying system can be approximated by the PWCS method with little loss of the characteristics of the time response. The SVD method is widely used in the SINS, because the filtering process is not needed at all before the observability analysis. However, both of these methods have their limitations. The PWCS method is only able to get the number of the observed system states, while the specific observability degree cannot be determined. For the incompletely observed systems in which the states are coupled, the SVD method cannot explicitly explain the coupling [20–22].

To perfect the above methods' deficiencies, numerous researchers and scholars are dedicating themselves to further researching and improving the methods of observability analysis. One of the most notably modified methods is the global observability analysis method based on the original nonlinear models, proposed in 2011 by Wu et al. [9,23,24]. It should be noted that in most past research, conservative observability concepts, e.g., local observability and linear observability, have extensively been used to locally characterize the estimability. However, the global observability analysis can provide us with comprehensive instructions on whether the estimation is feasible under a given condition or not, as well as how to achieve the state estimates, e.g., by resorting to vehicle maneuvers. Additionally, the global observability analysis method can be used in nonlinear systems. Directly starting from the definition of the observability, this method takes the nonlinear models of the SINS as the basic equations of the observability analysis. After the basic equations are solved, the global observability of the SINS can be analyzed. The global observability analysis method can conduct the observability analysis of the SINS simply, directly and effectively, by avoiding the drawbacks of the traditional methods.

A self-calibration method for dual-axis rotation-modulating INS was introduced in [25], and a calibration rotation sequence was designed. Sun et al. [26] proposed an eight-position calibration method, which can achieve calibration inertial sensor errors when the attitudes were unknown. However, the accuracy of these existing methods was not high, so we proposed a new rotation method for self-calibration and initial alignment.

Based on the global observability analysis, the relationship between resident positions and rotations of the IMU and the relationship between the observability and the kinestate of the rotating SINS have been established in this paper. Accordingly, the design principle of the initial alignment and self-calibration rotating scheme for the SINS was proposed. Besides, by taking the proposed design principle into account, an optimum rotating scheme was designed in this paper. The numerical simulations and experiments showed that the proposed rotary scheme was superior to others, and the position error can be reduced significantly with this new rotary scheme. The rest of the paper is organized as follows. The design principle of the rotating scheme of the SINS is described in Section 2. Section 3 proposes a new rotating scheme based on the design principle. Numerical examples and experiments along with specific analysis are given in Section 4 and Section 5, respectively. Section 6 concludes this manuscript.

2. Design Principle of the Rotating Scheme

With the design of the initial alignment and self-calibration of the dual-axis rotary SINS, the key is to determine resident positions and to choose rotating axes [9], which will be discussed respectively in this section.

2.1. Determination of Resident Positions

Assume that inertial sensor errors are just constant, then the output of the gyroscopes can be expressed as:

$${\stackrel{\sim }{\mathit{\omega }}}_{ib}^b={\mathit{\omega }}_{ib}^b+\mathit{\varepsilon }$$

(1)

wherein i indicates the inertial coordinate system, while b indicates the body coordinate system; and vectors ${\stackrel{\sim }{\mathit{\omega }}}_{ib}^b$, ${\mathit{\omega }}_{ib}^b$ and ε are the raw output, the theoretical output and the constant drift of the gyroscopes, respectively.

From Equation (1), we can get:

$$\left|{\stackrel{\sim }{\mathit{\omega }}}_{ib}^b-\mathit{\varepsilon }\right|=\left|{\mathit{\omega }}_{ib}^b\right|=\mathbf{\Omega }$$

(2)

Equation (2) can be explained as a spherical surface, whose center is the endpoint of vector ${\stackrel{\sim }{\mathit{\omega }}}_{ib}^b$ with the radius Ω, while ε can be explained as a moving point on this spherical surface. Therefore, the observability of ε is about determining the unique point on this spherical surface.

If there are three spherical surfaces intersected with each other, the number of the intersection points is two. Moreover, these two intersection points are symmetrical about the plane formed by three sphere centers. Therefore, it is not enough if there are only three constraints. If we want to determine the unique point further, another spherical surface whose center is not in the plane should be introduced [13]. Hence, ε can be shown as:

$$\mathit{\varepsilon }={\left(\begin{array}{cc}{\left[\begin{array}{c}2{\left[{\left({\stackrel{\sim }{\mathit{\omega }}}_{ib}^b\right)}_1-{\left({\stackrel{\sim }{\mathit{\omega }}}_{ib}^b\right)}_2\right]}^T\\ 2{\left[{\left({\stackrel{\sim }{\mathit{\omega }}}_{ib}^b\right)}_1-{\left({\stackrel{\sim }{\mathit{\omega }}}_{ib}^b\right)}_3\right]}^T\\ 2{\left[{\left({\stackrel{\sim }{\mathit{\omega }}}_{ib}^b\right)}_1-{\left({\stackrel{\sim }{\mathit{\omega }}}_{ib}^b\right)}_4\right]}^T\end{array}\right]}^T& \left[\begin{array}{c}2{\left[{\left({\stackrel{\sim }{\mathit{\omega }}}_{ib}^b\right)}_1-{\left({\stackrel{\sim }{\mathit{\omega }}}_{ib}^b\right)}_2\right]}^T\\ 2{\left[{\left({\stackrel{\sim }{\mathit{\omega }}}_{ib}^b\right)}_1-{\left({\stackrel{\sim }{\mathit{\omega }}}_{ib}^b\right)}_3\right]}^T\\ 2{\left[{\left({\stackrel{\sim }{\mathit{\omega }}}_{ib}^b\right)}_1-{\left({\stackrel{\sim }{\mathit{\omega }}}_{ib}^b\right)}_4\right]}^T\end{array}\right]\end{array}\right)}^{-1}{\left[\begin{array}{c}2{\left[{\left({\stackrel{\sim }{\mathit{\omega }}}_{ib}^b\right)}_1-{\left({\stackrel{\sim }{\mathit{\omega }}}_{ib}^b\right)}_2\right]}^T\\ 2{\left[{\left({\stackrel{\sim }{\mathit{\omega }}}_{ib}^b\right)}_1-{\left({\stackrel{\sim }{\mathit{\omega }}}_{ib}^b\right)}_3\right]}^T\\ 2{\left[{\left({\stackrel{\sim }{\mathit{\omega }}}_{ib}^b\right)}_1-{\left({\stackrel{\sim }{\mathit{\omega }}}_{ib}^b\right)}_4\right]}^T\end{array}\right]}^T\left[\begin{array}{c}{\left|{\left({\stackrel{\sim }{\mathit{\omega }}}_{ib}^b\right)}_1\right|}^2-{\left|{\left({\stackrel{\sim }{\mathit{\omega }}}_{ib}^b\right)}_2\right|}^2\\ {\left|{\left({\stackrel{\sim }{\mathit{\omega }}}_{ib}^b\right)}_1\right|}^2-{\left|{\left({\stackrel{\sim }{\mathit{\omega }}}_{ib}^b\right)}_3\right|}^2\\ {\left|{\left({\stackrel{\sim }{\mathit{\omega }}}_{ib}^b\right)}_1\right|}^2-{\left|{\left({\stackrel{\sim }{\mathit{\omega }}}_{ib}^b\right)}_4\right|}^2\end{array}\right]$$

(3)

wherein ${\left({\stackrel{\sim }{\mathit{\omega }}}_{ib}^b\right)}_j$ is the output when the gyroscopes are in the j—th (j = 1, 2,…, n) position.

Homoplastically, the zero bias of the accelerometers ∇ can also be determined uniquely. The output vector of the accelerometers is:

$${\stackrel{\sim }{\mathit{f}}}_{ib}^b={\mathit{f}}_{ib}^b+\mathbf{\nabla }$$

(4)

wherein vectors ${\stackrel{\sim }{\mathit{f}}}_{ib}^b$ and ${\mathit{f}}_{ib}^b$ are the raw output and the theoretical output of the accelerometers, respectively.

In a similar way to ε, ∇ can be deduced as:

$$\mathbf{\nabla }={\left(\begin{array}{cc}{\left[\begin{array}{c}2{\left[{\left({\stackrel{\sim }{\mathit{f}}}_{ib}^b\right)}_1-{\left({\stackrel{\sim }{\mathit{f}}}_{ib}^b\right)}_2\right]}^T\\ 2{\left[{\left({\stackrel{\sim }{\mathit{f}}}_{ib}^b\right)}_1-{\left({\stackrel{\sim }{\mathit{f}}}_{ib}^b\right)}_3\right]}^T\\ 2{\left[{\left({\stackrel{\sim }{\mathit{f}}}_{ib}^b\right)}_1-{\left({\stackrel{\sim }{\mathit{f}}}_{ib}^b\right)}_4\right]}^T\end{array}\right]}^T& \left[\begin{array}{c}2{\left[{\left({\stackrel{\sim }{\mathit{f}}}_{ib}^b\right)}_1-{\left({\stackrel{\sim }{\mathit{f}}}_{ib}^b\right)}_2\right]}^T\\ 2{\left[{\left({\stackrel{\sim }{\mathit{f}}}_{ib}^b\right)}_1-{\left({\stackrel{\sim }{\mathit{f}}}_{ib}^b\right)}_3\right]}^T\\ 2{\left[{\left({\stackrel{\sim }{\mathit{f}}}_{ib}^b\right)}_1-{\left({\stackrel{\sim }{\mathit{f}}}_{ib}^b\right)}_4\right]}^T\end{array}\right]\end{array}\right)}^{-1}\begin{array}{cc}{\left[\begin{array}{c}2{\left[{\left({\stackrel{\sim }{\mathit{f}}}_{ib}^b\right)}_1-{\left({\stackrel{\sim }{\mathit{f}}}_{ib}^b\right)}_2\right]}^T\\ 2{\left[{\left({\stackrel{\sim }{\mathit{f}}}_{ib}^b\right)}_1-{\left({\stackrel{\sim }{\mathit{f}}}_{ib}^b\right)}_3\right]}^T\\ 2{\left[{\left({\stackrel{\sim }{\mathit{f}}}_{ib}^b\right)}_1-{\left({\stackrel{\sim }{\mathit{f}}}_{ib}^b\right)}_4\right]}^T\end{array}\right]}^T& \left[\begin{array}{c}{\left|{\left({\stackrel{\sim }{\mathit{f}}}_{ib}^b\right)}_1\right|}^2-{\left|{\left({\stackrel{\sim }{\mathit{f}}}_{ib}^b\right)}_2\right|}^2\\ {\left|{\left({\stackrel{\sim }{\mathit{f}}}_{ib}^b\right)}_1\right|}^2-{\left|{\left({\stackrel{\sim }{\mathit{f}}}_{ib}^b\right)}_3\right|}^2\\ {\left|{\left({\stackrel{\sim }{\mathit{f}}}_{ib}^b\right)}_1\right|}^2-{\left|{\left({\stackrel{\sim }{\mathit{f}}}_{ib}^b\right)}_4\right|}^2\end{array}\right]\end{array}$$

(5)

wherein ${\left({\stackrel{\sim }{\mathit{f}}}_{ib}^b\right)}_j$ is the accelerometer output at the j—th (j = 1, 2,…, n) position.

Based on the above analysis, we can learn how to determine the resident positions in principle. That is, the IMU should stay at four different positions at least. Only if this condition is satisfied can the gyroscope drifts and the accelerometer biases be observed completely.

2.2. The Choices of Rotary Axes

According to the fundamentals of the SINS, the basic equations of the global observability analysis are as follows [24]:

$$\{\begin{array}{c}{\dot{\mathit{C}}}_b^n={\mathit{C}}_b^n\left[\left({\stackrel{\sim }{\mathit{\omega }}}_{ib}^b-\mathit{\varepsilon }\right)\times \right]-\left[\left({\mathit{\omega }}_{ie}^n+{\mathit{\omega }}_{en}^n\right)\times \right]{\mathit{C}}_b^n\\ {\dot{\mathit{\upsilon }}}^n={\mathit{C}}_b^n\left({\stackrel{\sim }{\mathit{f}}}_{ib}^b-\mathbf{\nabla }\right)-\left(2{\mathit{\omega }}_{ie}^n+{\mathit{\omega }}_{en}^n\right)\times {\mathit{\upsilon }}^n+{\mathit{g}}^n\end{array}$$

(6)

wherein ${\mathit{C}}_b^n$ denotes the direction cosine matrix from the b coordinate system to the navigation coordinate system (n); vectors ${\stackrel{\sim }{\mathit{f}}}_{ib}^b$ and ${\stackrel{\sim }{\mathit{\omega }}}_{ib}^b$ are the raw outputs of the accelerometers and the gyroscopes, respectively; ${\mathit{\omega }}_{ie}^n$ indicates the Earth's rotating angular rate vector; ${\mathit{\omega }}_{en}^n$ is the rotating angular rate vector in the n coordinate system relative to the Earth coordinate system (e); gⁿ denotes the gravity vector in the n coordinate system; υⁿ is the transporter's velocity of the n coordinate system and is zero when the transporter is in the stationary base.

Considering that in the stationary base, the IMU is rotating at a uniform velocity around the sensitive axis, Equation (6) can be rewritten as follows:

$$\{\begin{array}{c}{\stackrel{\sim }{\mathit{\omega }}}_{ib}^b={\mathit{\omega }}_{nb}^b+{\mathit{C}}_n^b{\mathit{\omega }}_{ie}^n+\mathit{\varepsilon }\\ {\stackrel{\sim }{\mathit{f}}}_{ib}^b={\mathit{C}}_n^b{\mathit{g}}^n+\mathbf{\nabla }\end{array}$$

(7)

wherein ${\mathit{\omega }}_{nb}^b$ is the rotating angular rate vector in the b coordinate system relative to the n coordinate system.

The differentiation of Equation (7) gives:

$$\{\begin{array}{c}{\dot{\stackrel{\sim }{\mathit{\omega }}}}_{ib}^b={\dot{\mathit{\omega }}}_{nb}^b+{\dot{\mathit{C}}}_n^b{\mathit{\omega }}_{ie}^n={\dot{\mathit{C}}}_n^b{\mathit{\omega }}_{ie}^n=-{\mathit{\omega }}_{nb}^b\times \left({\mathit{C}}_n^b{\mathit{\omega }}_{ie}^n\right)\\ {\dot{\stackrel{\sim }{\mathit{f}}}}_{ib}^b=-{\dot{\mathit{C}}}_n^b{\mathit{g}}^n={\mathit{\omega }}_{nb}^b\left({\mathit{C}}_n^b{\mathit{g}}^n\right)\end{array}$$

(8)

Further, the differentiation of Equation (8) yields:

$$\{\begin{array}{c}{\ddot{\stackrel{\sim }{\mathit{\omega }}}}_{ib}^b=-{\mathit{\omega }}_{nb}^b\times \left(-{\mathit{\omega }}_{nb}^b\times {\mathit{C}}_n^b\right)=-{\mathit{\omega }}_{nb}^b\times {\dot{\stackrel{\sim }{\mathit{\omega }}}}_{ib}^b\\ {\ddot{\stackrel{\sim }{\mathit{f}}}}^b={\mathit{\omega }}_{nb}^b\times \left({\dot{\mathit{C}}}_n^b{\mathit{g}}^n\right)=-{\mathit{\omega }}_{nb}^b\times {\dot{\stackrel{\sim }{\mathit{f}}}}_{ib}^b\end{array}$$

(9)

Since ${\dot{\stackrel{\sim }{\mathit{\omega }}}}_{ib}^b=-{\mathit{\omega }}_{nb}^b\times \left({\mathit{C}}_n^b{\mathit{\omega }}_{ie}^n\right)$ and the property of the vector [9], we know that ${\mathit{C}}_n^b{\mathit{\omega }}_{ie}^n$ is a vector in the plane, which is formed by vectors ${\mathit{\omega }}_{nb}^b$ and ${\mathit{\omega }}_{nb}^b\times {\dot{\stackrel{\sim }{\mathit{\omega }}}}_{ib}^b$.

Therefore, ${\mathit{C}}_n^b{\mathit{\omega }}_{ie}^n$ can be expressed as:

$${\mathit{C}}_n^b{\mathit{\omega }}_{ie}^n=k_1{\mathit{\omega }}_{nb}^b+k_2\left({\mathit{\omega }}_{nb}^b\times {\dot{\stackrel{\sim }{\mathit{\omega }}}}_{ib}^b\right)$$

(10)

wherein both k₁ and k₂ are constant.

Furthermore, due to:

$${\mathit{\omega }}_{nb}^b\perp \left({\mathit{\omega }}_{nb}^b\times {\dot{\stackrel{\sim }{\mathit{\omega }}}}_{ib}^b\right)$$

(11)

we can get:

$${\mathit{\omega }}_{nb}^b\cdot \left({\mathit{\omega }}_{nb}^b\times {\dot{\stackrel{\sim }{\mathit{\omega }}}}_{ib}^b\right)=0$$

(12)

Furthermore, we can obtain:

$$k_2=-\frac{1}{{\left|{\mathit{\omega }}_{nb}^b\right|}^2}$$

(13)

After Equation (10) and:

$${\left|{\mathit{C}}_n^b{\mathit{\omega }}_{ie}^n\right|}^2={\mathbf{\Omega }}^2$$

(14)

We have:

$${\left[k_1{\mathit{\omega }}_{nb}^b+k_2\left({\mathit{\omega }}_{nb}^b\times {\dot{\stackrel{\sim }{\mathit{\omega }}}}_{ib}^b\right)\right]}^2={\mathbf{\Omega }}^2$$

(15)

By expanding the left side of Equation (15), the equation can be expressed as follows:

$${k_1}^2{\left|{\mathit{\omega }}_{nb}^b\right|}^2+2k_1{k}_2{\mathit{\omega }}_{nb}^b\left({\mathit{\omega }}_{nb}^b\times {\dot{\stackrel{\sim }{\mathit{\omega }}}}_{ib}^b\right)+{k_2}^2{\left({\mathit{\omega }}_{nb}^b\times {\dot{\stackrel{\sim }{\mathit{\omega }}}}_{ib}^b\right)}^2={\mathbf{\Omega }}^2$$

(16)

Since:

$${\mathit{\omega }}_{nb}^b\perp \left({\mathit{\omega }}_{nb}^b\times {\dot{\stackrel{\sim }{\mathit{\omega }}}}_{ib}^b\right)$$

(17)

Equation (16) is simplified to:

$${k_1}^2{\left|{\mathit{\omega }}_{nb}^b\right|}^2+{k_2}^2{\left|{\mathit{\omega }}_{nb}^b\right|}^2{\left|{\dot{\stackrel{\sim }{\mathit{\omega }}}}_{ib}^b\right|}^2={\mathbf{\Omega }}^2$$

(18)

From Equations (13) and (18), k₁ is as follows:

$$k_1=\pm \frac{\sqrt{{\left|{\mathit{\omega }}_{nb}^b\right|}^2{\mathbf{\Omega }}^2-{\left|{\dot{\stackrel{\sim }{\mathit{\omega }}}}_{ib}^b\right|}^2}}{{\left|{\mathit{\omega }}_{nb}^b\right|}^2}$$

(19)

According to Equations (10), (13) and (19), the following equation can easily be deduced:

$${\mathit{C}}_n^b{\mathit{\omega }}_{ie}^n=\pm \frac{\sqrt{{\left|{\mathit{\omega }}_{nb}^b\right|}^2{\mathbf{\Omega }}^2-{\left|{\dot{\stackrel{\sim }{\mathit{\omega }}}}_{ib}^b\right|}^2}}{{\left|{\mathit{\omega }}_{nb}^b\right|}^2}{\mathit{\omega }}_{nb}^b-\frac{1}{{\left|{\mathit{\omega }}_{nb}^b\right|}^2}{\mathit{\omega }}_{nb}^b\times {\dot{\stackrel{\sim }{\mathit{\omega }}}}_{ib}^b$$

(20)

Further deduction of Equation (20) gives:

$${\mathit{C}}_n^b{\mathit{\omega }}_{ie}^n={\mathit{C}}_n^b{\mathit{\omega }}_{ie}^n-\frac{{\mathit{\omega }}_{nb}^b\left({\mathit{C}}_n^b{\mathit{\omega }}_{ie}\right)}{{\left|{\mathit{\omega }}_{nb}^b\right|}^2}{\mathit{\omega }}_{nb}^b\pm \frac{{\mathit{\omega }}_{nb}^b\left({\mathit{C}}_n^b{\mathit{\omega }}_{ie}\right)}{{\left|{\mathit{\omega }}_{nb}^b\right|}^2}{\mathit{\omega }}_{nb}^b$$

(21)

Similarly, ${\mathit{C}}_n^b{\mathit{g}}^n$ can also be obtained as follows:

$${\mathit{C}}_n^b{\mathit{g}}^n={\mathit{C}}_n^b{\mathit{g}}^n-\frac{{\mathit{\omega }}_{nb}^b\left({\mathit{C}}_n^b{\mathit{g}}^n\right)}{{\left|{\mathit{\omega }}_{nb}^b\right|}^2}{\mathit{\omega }}_{nb}^b\pm \frac{{\mathit{\omega }}_{nb}^b\left({\mathit{C}}_n^b{\mathit{g}}^n\right)}{{\left|{\mathit{\omega }}_{nb}^b\right|}^2}{\mathit{\omega }}_{nb}^b$$

(22)

From Equations (21) and (22), ${\mathit{C}}_n^b{\mathit{\omega }}_{ie}$ and ${\mathit{C}}_n^b{\mathit{g}}^n$ can obtain two solutions:

$$\{\begin{array}{c}{\left({\mathit{C}}_n^b{\mathit{\omega }}_{ie}^n\right)}_1={\mathit{C}}_n^b{\mathit{\omega }}_{ie}^n\\ {\left({\mathit{C}}_n^b{\mathit{\omega }}_{ie}^n\right)}_2={\mathit{C}}_n^b{\mathit{\omega }}_{ie}^n-2\frac{{\mathit{\omega }}_{nb}^b\left({\mathit{C}}_n^b{\mathit{\omega }}_{ie}^n\right)}{{\left|{\mathit{\omega }}_{nb}^b\right|}^2}{\mathit{\omega }}_{nb}^b\end{array}$$

(23)

$$\{\begin{array}{c}{\left({\mathit{C}}_n^b{\mathit{g}}^n\right)}_1={\mathit{C}}_n^b{\mathit{g}}^n\\ {\left({\mathit{C}}_n^b{\mathit{g}}^n\right)}_2={\mathit{C}}_n^b{\mathit{g}}^n-2\frac{{\mathit{\omega }}_{nb}^b\left({\mathit{C}}_n^b{\mathit{g}}^n\right)}{{\left|{\mathit{\omega }}_{nb}^b\right|}^2}{\mathit{\omega }}_{nb}^b\end{array}$$

(24)

It is easy to know that only when each of ${\mathit{C}}_n^b{\mathit{\omega }}_{ie}^n$ and ${\mathit{C}}_n^b{\mathit{g}}^n$ has one solution uniquely, the ε and ∇ can be absolutely observed. Therefore, the situation that each of ${\mathit{C}}_n^b{\mathit{\omega }}_{ie}^n$ and ${\mathit{C}}_n^b{\mathit{g}}^n$ has two solutions is not acceptable. However, $\frac{{\mathit{\omega }}_{nb}^b\left({\mathit{C}}_n^b{\mathit{\omega }}_{ie}^n\right)}{|{\mathit{\omega }}_{nb}^b{|}^2}{\mathit{\omega }}_{nb}^b$ and $\frac{{\mathit{\omega }}_{nb}^b\left({\mathit{C}}_n^b{\mathit{g}}^n\right)}{|{\mathit{\omega }}_{nb}^b{|}^2}{\mathit{\omega }}_{nb}^b$ can be seen as corresponding errors, while ${\mathit{C}}_n^b{\mathit{\omega }}_{ie}^n$ and ${\mathit{C}}_n^b{\mathit{g}}^n$ as true values. Then, only if the errors are zeros, the solutions are unique and ε and ∇ are absolutely observed.

When ${\mathit{\omega }}_{nb}^b\perp \left({\mathit{C}}_n^b{\mathit{\omega }}_{ie}^n\right)$, that means, when the rotating direction of the IMU is perpendicular to the Earth's rotary axis, the value of $\frac{{\mathit{\omega }}_{nb}^b\left({\mathit{C}}_n^b{\mathit{\omega }}_{ie}^n\right)}{|{\mathit{\omega }}_{nb}^b{|}^2}{\mathit{\omega }}_{nb}^b$ is zero, and in this condition, ε is absolutely observed. When ${\mathit{\omega }}_{nb}^b\perp \left({\mathit{C}}_n^b{\mathit{g}}^n\right)$, that means, when the rotating direction of the IMU is on the horizontal plane, $\frac{{\mathit{\omega }}_{nb}^b\left({\mathit{C}}_n^b{\mathit{g}}^n\right)}{|{\mathit{\omega }}_{nb}^b{|}^2}{\mathit{\omega }}_{nb}^b$ becomes zero, so that ∇ is absolutely observed.

From the above analysis, we can get the choosing principle of the rotary axis in the rotary SINS clearly: when the rotary axes are horizontal axes, the gyroscope drifts and the accelerometer biases are both absolutely observed.

3. New Initial Alignment and Self-Calibration Rotating Scheme

From the analysis in Section 2, the gyroscope drifts, and the accelerometer biases and the misalignment angles are all absolutely observed when the IMU alternately rotates about two horizontal axes and resides at four different positions at least. Therefore, to ensure that all of inertial sensor errors can completely be observed and to minimize the rotation number, a novel eight-position rotating scheme rotating about two horizontal axes was proposed here based on the above rotary principle. In this scheme, the IMU is in rotation and steady alternately. There is a small static state that keeps t_stop after each rotation, wherein t_stop is the residence time at each position. The rotation order in the new scheme is illustrated in Figure 1.

The IMU rotation starts from Position A, and the detailed explanation of this new rotating scheme is as follows:

• A → B: rotate counterclockwise 180° about the x axis from Position A to Position B;

• B → C: rotate counterclockwise 90° about the x axis from Position B to Position C;

• C → D: rotate clockwise 180° about the x axis from Position C to Position D;

• D → A: rotate clockwise 90° about the x axis from Position D to Position A;

• A → E: rotate counterclockwise 90° about the y axis from Position A to Position E;

• E → F: rotate counterclockwise 180° about the x axis from Position E to Position F;

• F → G: rotate counterclockwise 90° about the x axis from Position F to Position G;

• G → H: rotate clockwise 180° about the x axis from Position G to Position H;

• H → E: rotate clockwise 90° about the x axis from Position H to Position E;

• E → A: rotate clockwise 90° about the y axis from Position E to Position A.

In this novel rotation scheme, the IMU rotates about two horizontal axes (x axis and y axis) and resides at eight different positions to satisfy the design principle introduced in the previous section.

4. Simulations and Results

In this section, simulations were introduced to check the feasibility of the proposed rotating scheme. From the discussions in previous sections, the initial alignment and self-calibration of the SINS can be seen as a system optimal estimation based on the Kalman filter, and the block diagram of the Kalman filter is presented as in Figure 2.

Thus, filtering equations of the SINS should first be established. The state equation of the SINS is symbolically as follows:

$$\dot{\mathit{X}}\left(t\right)=\mathit{F}\left(t\right)\mathit{X}\left(t\right)+\mathit{G}\left(t\right)\mathit{W}\left(t\right)$$

(25)

wherein X(t) is the state error vector at t; F(t) is the state-transition matrix; W(t) and G(t) are the noise vector and the associated coefficient matrix, respectively.

The state vector and the noise vector are defined as follows:

$$\mathit{X}\left(t\right)={\left[\mathit{\delta }L\mathit{\delta }\mathit{\lambda }\mathit{\delta }{\mathit{\upsilon }}_E\mathit{\delta }{\mathit{\upsilon }}_N{\mathit{\phi }}_x{\mathit{\phi }}_y{\mathit{\phi }}_z{\mathbf{\nabla }}_x{\mathbf{\nabla }}_y{\mathbf{\nabla }}_z{\mathit{\varepsilon }}_x{\mathit{\varepsilon }}_y{\mathit{\varepsilon }}_z\right]}^T$$

(26)

$$\mathit{W}\left(t\right)={\left[{\mathit{\omega }}_{{\mathbf{\nabla }}_x}{\mathit{\omega }}_{{\mathbf{\nabla }}_y}{\mathit{\omega }}_{{\mathbf{\nabla }}_z}{\mathit{\omega }}_{{\mathit{\varepsilon }}_x}{\mathit{\omega }}_{{\mathit{\varepsilon }}_y}{\mathit{\omega }}_{{\mathit{\varepsilon }}_z}\right]}^T$$

(27)

wherein δL and δλ are the errors of the longitude and the latitude, respectively; δυ_E and δυ_N are the east and north velocity errors; the misalignment angles between the calculated geographical coordinate system and the true geographical coordinate system are denoted by φ_x, φ_y and φ_z; ∇_x, ∇_y and ∇_z are the accelerometer biases in the x, y and z directions, respectively; ω∇_x, ω∇_y and ω∇_z are the corresponding noises of the accelerometers; ε_x, ε_y and ε_z are the gyroscope drifts in the x, y and z directions, respectively; ω_εx, ω_εy and ω_εz are the corresponding noises of the gyroscopes.

Let us take the velocity errors as the measurement vector; the measurement equation is as follows:

$$\mathit{Z}\left(t\right)=\mathit{H}\left(t\right)\mathit{X}\left(t\right)+\mathit{V}\left(t\right)$$

(28)

wherein Z is the measurement vector at t and H and V are the measurement design matrix and the measurement noise, respectively.

The matrix H is given as:

$$\mathit{H}\left(t\right)=\left[{\mathbf{0}}_{2\times 2}{\mathit{I}}_{2\times 2}{\mathbf{0}}_{2\times 9}\right]$$

(29)

Next, simulations were carried out to verify the effectiveness and superiority of the eight-position rotating scheme proposed in this paper. In order to make a better illustration, three groups of simulations and their results were compared below. In one simulation, the SINS was set as stationary. Contrarily, the IMU was set as rotary in the other two simulations. The only difference between these two simulations is in their rotating axes: the first one was rotating about the x axis and y axis proposed in this paper, while the second one was rotating about the y axis and z axis proposed in [26]. The settings of the rotary angles in these two rotary simulations are described in Figures 3 and 4.

The parameter settings of the simulations are shown in

Table 1. Simulation parameters.

Parameters	Parameter Settings
angular rate of rotation: υr	5 (°/s)
angular acceleration of rotation: ar	2.5 (°/s2)
residence time: tstop	500 s
constant drifts of gyroscopes: εx = εy = εz	0.01 (°/h)
zero-biases of accelerometers: ∇x = ∇y = ∇z	9.8 × 10×4 m/s2
noises of gyroscopes: ωεx = ωεy = ωεz	0.001 (°/h)
noises of accelerometers: ω∇x = ω∇y = ω∇z	9.8 × 10−5 m/s2
initial errors of horizontal misalignment angles: φx = φy	0.01°
initial error of azimuth misalignment angle: φz	0.05°
initial longitude: L0	126.67°
initial latitude: λ0	45.78°

. Figures 5, 6, 7, 8, 9 and 10 present the comparisons of the simulation results.

The correlation error curves of the misalignment angles of the three simulations are shown in Figures 5 and 6. When the IMU was set as stationary, the misalignment errors were relatively large. That is, in stationary mode, the misalignment angles cannot be estimated, while those in the rotating simulations can be estimated readily. As can be seen, in rotating simulations, the horizontal angle errors were almost similar, except that the curves had more waves when the IMU rotated about the y axis and z axis. Regarding the azimuth angle error, the simulation with the eight-position scheme rotating about the x axis and y axis can obtain a smaller estimated error.

In Figures 7 and 8, the error curves of the gyroscope drifts were compared. As shown in these figures, in stationary mode, only ε_y can be estimated, but with poor accuracy. The error of ε_y still had about 0.002°/h left. Additionally, ε_z was not restrained in 4 h, while ε_x was not estimated at all. On the other hand, compared with the results of the rotary simulations, the performances with the eight-position scheme rotating around the x axis and y axis were much better. With this scheme, the horizontal gyroscope drifts were convergent and steady within 1 h, while the azimuth gyroscope drift was convergent in 1.5 h. What is more, the convergence rate and the accuracy were both superior to the other rotating simulation.

The estimated results of the accelerometer biases are as in Figures 9 and 10, which clearly show that, when the IMU was in stationary mode, ∇_x, ∇_y and ∇_z cannot be estimated at all. With rotary simulations, the estimated results were approximately the same. Although the bias curves were all convergent within one hour, the superiority still existed in the scheme rotating around the x axis and y axis.

On the basis of the simulation results and analyses, firstly, if the IMU was rotating, a better estimation of the inertial sensor errors can be reached compared with the ones when the IMU was in stationary mode. Secondly, between these two rotating simulations, with the eight-position initial alignment and self-calibration scheme rotating about the x axis and y axis, the estimations of the inertial errors, including misalignment angles, gyroscope drifts and accelerometer biases, were much better. Therefore, the superiority and effectiveness of the rotating scheme were demonstrated evidently. Further, the feasibility of the design principle proposed in this paper was also checked at the same time.

5. Experiments and Results

In this section, experiments were carried on to further verify the superiority of the new eight-position rotary scheme. The SINS made by our lab and the SGT-3 turntable made by AVIC-precision 303 were used, the performance parameters of which are listed in

Table 2. Performance parameters of the turntable.

Parameters	Parameter Settings
load	20∼150 Kg
angular rate of rotation	0.0001∼1000 (°/s)
maximum of rotating angular acceleration	4000 (∼/s2)
positional accuracy	±3″
angular accuracy	±3″

.

The IMU remained static about 0.7 h and then rotated according to the rotary scheme proposed in Section 3. The inertial sensor errors and the misalignment angles were estimated by using the Kalman filter, and the results are shown in Figures 11, 12 and 13.

It can be seen obviously from Figures 11, 12 and 13 that the estimated curves were all stable in 4 h. Thus, in order to verify the estimated results, the estimated results at 4 h were compensated for and the position errors before and after the compensation were compared, as shown in Figure 14.

From Figure 14, we can see that, before compensation, the position error was about 7 nm, while after compensation, the position error was about 3 nm. Therefore, with the new rotary scheme proposed in this paper, the misalignment angles, the gyroscopic drifts and the accelerometer zero-biases can be estimated effectively. Further, the correctness and availability of the new rotary scheme were also verified.

6. Conclusions

In this manuscript, the design principle of a rotating scheme in a dual-axis rotary SINS was proposed in order to improve the slow convergence speed and low estimation accuracy of inertial errors during the initial alignment and self-calibration. On the basis of this principle, a new eight-position rotating scheme was designed. A static simulation and another rotating simulation were used as comparisons, and the only difference between these two rotating simulations was that they have different rotating axes. In these rotary simulations, the IMU periodically rotates about the rotating axes to enhance the observability of the SINS. To further verify and compare, numerical simulations and experiments were carried out. The simulation results showed that better performance can be achieved when the IMU was rotating than when the IMU was static. Although the errors of the inertial sensors and the misalignment angles can effectively be estimated in three rotating simulations, the eight-position rotating scheme rotating about the x axis and y axis had the best performance. Furthermore, the experiment results showed that the position accuracy of the whole SINS can be increased significantly with this new initial alignment and self-calibration rotating scheme. However, the effect of the scale factor error was not taken into account in this study and will further be studied as one of the future tasks.