Self-calibration method for temperature errors in multi-axis rotational inertial navigation system

Jingxuan Ban; Lei Wang; Zengjun Liu; Lixiang Zhang

doi:10.1364/OE.384905

1. Introduction

Compared with strapdown inertial navigation system, rotational inertial navigation system (RINS) is an effective way to realize significant navigation precision improvement using inertial devices of same performance [].

However, inertial sensor’s error still has significant influence on navigation accuracy. Coupled with rotation, navigation errors are magnified by the scale factor error of gyroscope that along with rotary axis and part of installation errors [2]. In order to promote system performance, research in self-calibration should be considered as an important part.

Besides, temperature will affect the performance of inertial sensors. Drifts, biases and scale factor errors are easily affected by working temperature, especially for quartz accelerometer and fiber optic gyroscope (FOG) [3–10]. Even for ring-laser gyroscope (RLG), temperature resulted errors should be considered carefully [11–13]. Moreover, the orientation of FOG input axis will also change over temperature, which leads to the change of installation errors [14,15]. If error parameters were calibrated under constant temperature and remain unchanged, navigation error would increase with working temperature varies. Thus, it is necessary to update error parameters according to real-time temperature. There are three possible solutions to decrease temperature errors: 1) temperature-control method [16]; 2) temperature model identification; 3) temperature-insensitive materials research [17]. Since the last one is device-level solution and the first one increases cost and weight [3], this paper mainly focuses on temperature model identification.

References [3,8,9,11] studied the characteristics of temperature errors directly from inertial sensors’ output with artificial intelligence algorithms, including neural network and support vector machine. But such algorithms are difficult to realize in real-time system, always sacrificed with cost and complexity. In addition, they require large amount of data which will be harmful to time performance.

References [4–6,10] established polynomial model to describe the temperature error for measurements of inertial sensors. References [12,13] employed discrete calibration method to identify temperature model. However, among these references, references [5,6] need two temperature sensors to estimate temperature drift that is not always available in every application. References [4,12,13] ask for temperature stabilization. References [12,13] require keeping temperature stable for more than 3 hours at every constant temperature point.

Compared with traditional methods, this paper proposes a system-level self-calibration method for temperature errors in multi-axis RINS. The method requires one temperature sensor for each inertial device. More important, the method doesn’t need hours to remain temperature stable, and the temperature can change at a constant rate during calibration, which increases production efficiency. Twenty-one error parameters studied in the paper are listed in Table .

Table 1. Error parameters studied in paper.

View Table | View all tables in this article

According to system-level self-calibration theory, a rotation scheme should be designed to excite navigation errors. Multi-position schemes in references [18,19] can be used to estimate all 21 errors within 30 minutes. However, due to the modulation technology, the calibration accuracy of gyroscope drifts can be lower than other error parameters. Thus, rotation scheme can be simplified and self-calibration time can be shortened further.

Self-calibration is a process of optimal estimation. Inertial measurement unit (IMU) error model in the paper is linear and INS error equations are given in differential form. If recursive least squares (RLS) algorithm was adopted, analytical expression of INS error equations should be given [20]. Compared with RLS, Kalman Filter algorithm is briefer. Therefore, the paper applies Kalman Filter algorithm to identify 21 errors of IMU. The system observability analysis is given by piecewise constant system analysis (PWCS) [21,22] and singular value decomposition (SVD) [23,24].

The rest of this paper is organized as follows. Section 2 defines coordinate systems, error parameters and error model of the system. Section 3 introduces Kalman Filter construction, self-calibration rotation scheme and analyzes system observability. Section 4 shows the simulation results, including single self-calibration results at constant temperature and temperature model identification results. Experimental results and analysis are given in Section 5. Finally, the conclusion is given in Section 6.

2. Coordinate system and error model

2.1 Coordinate system definition

Before analyzing installation errors, several coordinate systems must be defined, besides three common coordinate systems called inertial frame (i-frame), navigation frame (n-frame, defined as east-north-up) and body frame (b-frame, defined as right-forward-upward).

The first one is gyroscope measurement frame (mg-frame), which is defined by input axes. Accelerometer measurement frame (ma-frame) is defined similarly. Both ma-frame and mg-frame are nonorthogonal frame. The last one is platform frame (p-frame). Contrary to ma-frame and mg-frame, p-frame is an orthogonal frame. X-axis of p-frame directs to projection of X_mg on O-X_maY_ma plane, Z-axis of p-frame is perpendicular to O-X_maY_ma plane, Y-axis is defined according to right-hand rule. The installation errors of gyroscopes and accelerometers are shown in Fig. . $C_{mg}^p$ and $C_{ma}^p$ describe the rotation matrix from mg-frame and ma-frame to p-frame, respectively.

(1)$$C_{ma}^p = \left[ {\begin{array}{ccc} 1 &{ - \alpha_{ax}^Z} &0\\ {\alpha_{ay}^Z} &1 &0\\ { - \delta_{az}^Y}&{\delta_{az}^X} &1 \end{array}} \right] \qquad \qquad C_{mg}^p = \left[ {\begin{array}{ccc} 1 &0 &{\beta_{gx}^Y}\\ {\alpha_{gy}^Z} &1 &{ - \beta_{gy}^X}\\ { - \delta_{gz}^Y}&{\delta_{gz}^X} &1 \end{array}} \right]$$

Fig. 1. The definition of installation errors. (a) Gyroscope installation errors; (b) Accelerometer installation errors.

Download Full Size | PDF

2.2 IMU error model

IMU error model is described as Eq. (2) and Eq. (3).

(2)$$\delta \omega _{iP}^P = (C_{mg}^P - I + \Delta {K_g})\omega _{im}^m + \varepsilon $$

(3)$$\delta f_{iP}^P = (C_{ma}^P - I + \Delta {K_a})f_{im}^m + \nabla $$

where $\delta \omega _{iP}^P = {\left[ {\begin{array}{@{}ccc@{}} {\delta \omega_{iPx}^P}&{\delta \omega_{iPy}^P}&{\delta \omega_{iPz}^P} \end{array}} \right]^T}$ is gyroscope output errors and $\delta f_{iP}^P = {\left[ {\begin{array}{@{}ccc@{}} {\delta f_{iPx}^P}&{\delta f_{iPy}^P}&{\delta f_{iPz}^P} \end{array}} \right]^T}$ is accelerometer output errors. $\omega _{im}^m = {\left[ {\begin{array}{ccc} {\omega_{imx}^m}&{\omega_{imy}^m}&{\omega_{imz}^m} \end{array}} \right]^T}$ and $f_{im}^m = {\left[ {\begin{array}{ccc} {f_{imx}^m}&{f_{imy}^m}&{f_{imz}^m} \end{array}} \right]^T}$ are angular velocity input and acceleration input under measurement frame, respectively. $\varepsilon = {[\begin{array}{ccc} {{\varepsilon _{gx}}}&{{\varepsilon _{gy}}}&{{\varepsilon _{gz}}} \end{array}]^T}$ is gyroscope drift. $\nabla = {[\begin{array}{ccc} {{\nabla _{ax}}}&{{\nabla _{ay}}}&{{\nabla _{az}}} \end{array}]^T}$ is accelerometer bias. $\Delta {K_g} = {\left[ {\begin{array}{ccc} {\Delta {K_{gx}}}&{\Delta {K_{gy}}}&{\Delta {K_{gz}}} \end{array}} \right]^T}$ is gyroscope scale factor error. $\Delta {K_a} = {\left[ {\begin{array}{ccc} {\Delta {K_{ax}}}&{\Delta {K_{ay}}}&{\Delta {K_{az}}} \end{array}} \right]^T}$ is accelerometer scale factor error.

Attitude error equation [25] is obtained by Eq. (4).

(4)$$\dot{\phi } = \phi \times \omega _{in}^n + \delta \omega _{in}^n - C_p^n\delta \omega _{iP}^P$$

where $\phi = {\left[ {\begin{array}{ccc} {{\phi_E}}&{{\phi_N}}&{{\phi_U}} \end{array}} \right]^T}$ denotes attitude error. $C_p^n$ is transformation matrix from p-frame to n-frame. $\omega _{in}^n$ is relative angular velocity between n-frame and i-frame as shown in Eq. (5).

(5)$$\omega _{in}^n = \omega _{ie}^n + \omega _{en}^n = \left[ {\begin{array}{c} { - \frac{{{V_N}}}{{{R_M} + h}}}\\ {{\omega_{ie}}\cos L + \frac{{{V_E}}}{{{R_N} + h}}}\\ {{\omega_{ie}}\sin L + \frac{{{V_E}}}{{{R_N} + h}}\tan L} \end{array}} \right]$$

where ${\omega _{ie}}$ is Earth rotation angular velocity. L denotes latitude. R_M and R_N are radius of curvature in prime vertical and radius of curvature in meridian, respectively.

Velocity error equation [25] is written as Eq. (6).

(6)$$\delta {\dot{V}^n} ={-} {\phi ^n} \times {f^n} - (2\omega _{ie}^n + \omega _{en}^n) \times \delta {V^n} + {V^n} \times (2\delta \omega _{ie}^n + \delta \omega _{en}^n) - C_p^n\delta f_{iP}^P$$

where $\delta {V^n} = {\left[ {\begin{array}{ccc} {\delta {V_E}}&{\delta {V_N}}&{\delta {V_U}} \end{array}} \right]^T}$ denotes velocity error. ${f^n} = {\left[ {\begin{array}{ccc} {{f_E}}&{{f_N}}&{{f_U}} \end{array}} \right]^T}$ is the projection of the body’s specific force in n-frame. ${V^n} = {\left[ {\begin{array}{ccc} {V_E^{}}&{V_N^{}}&{V_U^{}} \end{array}} \right]^T}$ is carrier’s velocity.

$\delta \omega _{ie}^n$ and $\delta \omega _{en}^n$ is given in Eq. (7) and Eq. (8).

(7)$$\delta \omega _{ie}^n = {\left[ {\begin{array}{ccc} 0 &{ - \delta L{\omega_{ie}}\sin L}&{\delta L{\omega_{ie}}\cos L} \end{array}} \right]^T}$$

(8)$$\delta \omega _{en}^n = {\left[ {\begin{array}{ccc} { - \frac{{\delta {V_N}}}{{{R_M} + h}}}&{\frac{{\delta {V_E}}}{{{R_N} + h}}}&{\frac{{\delta {V_E}}}{{{R_N} + h}}\tan L} \end{array}} \right]^T}$$

Position error equation [25] is written as Eq. (9) and Eq. (10).

(9)$$\delta \dot{L} = \frac{{\delta {V_N}}}{{{R_M} + h}}$$

(10)$$\delta \dot{\lambda } = \frac{{\delta {V_E}}}{{{R_N} + h}}\sec L$$

where $\delta L$ and $\delta \lambda$ are latitude error and longitude error, respectively.

3. Self-calibration scheme

3.1 Kalman filter construction

Attitude reference is not always available, which leads to lack of attitude error. Moreover, height calculation of INS is divergent [25]. Therefore, velocity errors, latitude error and longitude error are selected as measurements in this paper.

(11)$$Z = {[{\delta {V_E}\;\delta {V_N}\;\delta {V_U}\;\delta L\;\delta \lambda } ]^T}$$

According to error model in the Section 2, 29 state variables are selected as Eq. (12).

(12)$$\begin{aligned}X &= [{{\phi_E}\;{\phi_N}\;{\phi_U}\;\delta {V_E}} \;\delta {V_N}\;\delta {V_U}\;\delta L\;\delta \lambda \;\;{\varepsilon _{gx}}\;{\varepsilon _{gy}}\;{\varepsilon _{gz}}\;{\nabla _{ax}}\;{\nabla _{ay}}\;{\nabla _{az}}\;{K_{gx}}\;{K_{gy}}\;{K_{gz}}{K_{ax}}\;{K_{ay}}\;{K_{az}}\;\\ &\beta _{gx}^Y\;\alpha _{gy}^Z\;\beta _{gy}^X\;\delta _{gz}^Y\;\delta _{gz}^X\;\alpha _{ax}^Z\;\alpha _{ay}^Z\;\delta _{az}^Y\; {\delta_{az}^X} ]_{29 \times 1}^T \end{aligned}$$

State equations and measurement equations are shown in Eq. (13).

(13)$$\begin{array}{l} \dot{X} = FX + W\\ Z = HX + V \end{array}$$

F is state transition matrix given in Eq. (14).

(14)$$F = {\left[ {\begin{array}{ccccccc} {F{1_{3 \times 8}}}&{C_p^n}&{{0_{3 \times 3}}}&{F{2_{3 \times 3}}}&{{0_{3 \times 3}}}&{F{3_{3 \times 5}}}&{{0_{3 \times 4}}}\\ {F{4_{3 \times 8}}}&{{0_{3 \times 3}}}&{C_p^n}&{{0_{3 \times 3}}}&{F{5_{3 \times 3}}}&{{0_{3 \times 5}}}&{F{6_{3 \times 4}}}\\ {F{7_{2 \times 8}}}&{{0_{2 \times 3}}}&{{0_{2 \times 3}}}&{{0_{2 \times 3}}}&{{0_{2 \times 3}}}&{{0_{2 \times 5}}}&{{0_{2 \times 4}}}\\ {{0_{21 \times 8}}}&{{0_{21 \times 3}}}&{{0_{21 \times 3}}}&{{0_{21 \times 3}}}&{{0_{21 \times 3}}}&{{0_{21 \times 5}}}&{{0_{21 \times 4}}} \end{array}} \right]_{29 \times 29}}$$

F1, F4 and F7 are derived from attitude error equation, velocity error equation and position error equation.

F2 indicates the relationship between gyroscope scale factor errors and attitude errors. F5 expresses velocity errors lead by accelerometer scale factor errors. They are written in Eq. (15).

(15)$$F2 = \left[ {\begin{array}{@{}ccc@{}} {C_p^n({1,1} )\omega_{imx}^m}&{C_p^n({1,2} )\omega_{imy}^m}&{C_p^n({1,3} )\omega_{imz}^m}\\ {C_p^n({2,1} )\omega_{imx}^m}&{C_p^n({2,2} )\omega_{imy}^m}&{C_p^n({2,3} )\omega_{imz}^m}\\ {C_p^n({3,1} )\omega_{imx}^m}&{C_p^n({3,2} )\omega_{imy}^m}&{C_p^n({3,3} )\omega_{imz}^m} \end{array}} \right]\,F5 = \left[ {\begin{array}{@{}ccc@{}} {C_p^n({1,1} )f_{imx}^m}&{C_p^n({1,2} )f_{imy}^m}&{C_p^n({1,3} )f_{imz}^m}\\ {C_p^n({2,1} )f_{imx}^m}&{C_p^n({2,2} )f_{imy}^m}&{C_p^n({2,3} )f_{imz}^m}\\ {C_p^n({3,1} )f_{imx}^m}&{C_p^n({3,2} )f_{imy}^m}&{C_p^n({3,3} )f_{imz}^m} \end{array}} \right]$$

F3 and F6 derived from attitude errors and velocity errors caused by installation errors. They are given in Eq. (16).

(16)$$F3 = {\left[ {\begin{array}{@{}ccc@{}} {C_p^n({1,1} )\omega_{imz}^m} &{C_p^n({2,1} )\omega_{imz}^m}&{C_p^n({3,1} )\omega_{imz}^m}\\ {C_p^n({1,2} )\omega_{imx}^m} &{C_p^n({2,2} )\omega_{imx}^m} &{C_p^n({3,2} )\omega_{imx}^m}\\ { - C_p^n({1,2} )\omega_{imz}^m} &{ - C_p^n({2,2} )\omega_{imz}^m} &{ - C_p^n({3,2} )\omega_{imz}^m}\\ { - C_p^n({1,3} )\omega_{imx}^m} &{ - C_p^n({2,3} )\omega_{imx}^m} &{ - C_p^n({3,3} )\omega_{imx}^m}\\ {C_p^n({1,3} )\omega_{imy}^m} &{C_p^n({2,3} )\omega_{imy}^m} &{C_p^n({3,3} )\omega_{imy}^m} \end{array}} \right]^T} \,F6 = {\left[ {\begin{array}{@{}ccc@{}} { - C_p^n({1,1} )f_{imy}^m} &{ - C_p^n({2,1} )f_{imy}^m} &{ - C_p^n({3,1} )f_{imy}^m}\\ {C_p^n({1,2} )f_{imx}^m} &{C_p^n({2,2} )f_{imx}^m} &{C_p^n({3,2} )f_{imx}^m}\\ { - C_p^n({1,3} )\omega_{imx}^m} &{ - C_p^n({2,3} )\omega_{imx}^m} &{ - C_p^n({3,3} )\omega_{imx}^m}\\ {C_p^n({1,3} )\omega_{imy}^m} &{C_p^n({2,3} )\omega_{imy}^m} &{C_p^n({3,3} )\omega_{imy}^m} \end{array}} \right]^T}$$

V and W are measurement noise and state noise, respectively, which are regarded as white noises. H is measurement matrix.

(17)$$H = {\left[ {\begin{array}{ccc} {{0_{\textrm{5} \times \textrm{3}}}}&{{I_{5 \times 5}}}&{{0_{5 \times \textrm{21}}}} \end{array}} \right]_{5 \times 29}}$$

where ${I_{5 \times 5}}$ denotes the identity matrix.

3.2 Design of rotation scheme

An eight-step rotation scheme was designed to excite all 21 errors, with rotation angular velocity is 4°/s. The calibration costs 585 seconds. The rotation scheme is shown in Fig. 2.

Fig. 2. Rotation Scheme.

Download Full Size | PDF

Taking step 1 and step 2 as an example to illustrate how the rotation scheme works. Suppose the initial attitude of IMU is zero, and rotation angular velocity is ${\omega _r}$. Then, ideal input of gyroscope is given in Eq. (18).

(18)$$\omega _{ip}^p = C_b^p\left[ {\begin{array}{c} 0\\ {{\omega_{ie}}\cos L + {\omega_r}}\\ {{\omega_{ie}}\sin L} \end{array}} \right] = C_b^p\left[ {\begin{array}{c} 0\\ {{\omega_N} + {\omega_r}}\\ {{\omega_U}} \end{array}} \right]$$

(19)$$C_b^p = \left[ {\begin{array}{ccc} {\cos {\omega_r}t} &0 &{ - \sin {\omega_r}t}\\ 0 &1 &0\\ {\sin {\omega_r}t} &0 &{\cos {\omega_r}t} \end{array}} \right]$$

Equivalent drifts in navigation frame are written in Eq. (20).

(20)$$\left\{ {\begin{array}{{l}} \begin{array}{l} \varepsilon _x^n = \frac{{\left( {{\Delta }{K_{gz}} - {\Delta }{K_{gx}}} \right)}}{2}{\omega _U}\sin 2{\omega _r}t + \frac{{\left( {\beta _{gx}^Y + \delta _{gz}^Y} \right)}}{2}{\omega _U} + \frac{{\left( {\beta _{gx}^Y - \delta _{gz}^Y} \right)}}{2}{\omega _U}\cos 2{\omega _r}t\\ \qquad+ \delta _{gz}^X\left( {{\omega _U} + {\omega _r}} \right)\sin {\omega _r}t + {\varepsilon _{gx}}\cos {\omega _r}t + {\varepsilon _{gz}}\sin {\omega _r}t \end{array}\\ {\varepsilon _y^n = - \alpha _{gy}^Z{\omega _U}\sin {\omega _r}t + {\Delta }{K_{gy}}\left( {{\omega _N} + {\omega _r}} \right) + \beta _{gy}^X{\omega _U}\cos {\omega _r}t + {\varepsilon _{gy}}}\\ \begin{array}{l} \varepsilon _z^n = \frac{{\left( {{\Delta }{K_{gx}} + {\Delta }{K_{gz}}} \right)}}{2}{\omega _U} + \frac{{\left( {{\Delta }{K_{gz}} - {\Delta }{K_{gx}}} \right)}}{2}{\omega _U}\cos 2{\omega _r}t + \frac{{\left( {\delta _{gz}^Y - \beta _{gx}^Y} \right)}}{2}{\omega _U}\sin 2{\omega _r}t\\ \qquad+ \delta _{gz}^X\left( {{\omega _U} + {\omega _r}} \right)\cos {\omega _r}t - {\varepsilon _{gx}}\sin {\omega _r}t + {\varepsilon _{gz}}\cos {\omega _r}t \end{array} \end{array}} \right.$$

Since ${\omega _U} < < {\omega _r}$ and ${\omega _N} < < {\omega _r}$, integration of such terms, including ${\omega _U}\sin 2{\omega _r}t$, ${\omega _U}\cos 2{\omega _r}t$, ${\varepsilon _{gx}}\cos {\omega _r}t$, ${\varepsilon _{gz}}\sin {\omega _r}t$, ${\varepsilon _{gx}}\sin {\omega _r}t$ and ${\varepsilon _{gz}}\cos {\omega _r}t$, are small enough to be ignored. Thus, Eq. (20) can be reduced to Eq. (21).

(21)$${\varepsilon ^n} = \left[ {\begin{array}{c} {\varepsilon_x^n}\\ {\varepsilon_y^n}\\ {\varepsilon_z^n} \end{array}} \right]\textrm{ = }\left[ {\begin{array}{c} {\frac{{({\beta_{gx}^Y + \delta_{gz}^Y} )}}{2}{\omega_U} + \delta_{gz}^X{\omega_r}\sin {\omega_r}t}\\ {\Delta {K_{gy}}{\omega_r} + {\varepsilon_{gy}}}\\ {\frac{{({\Delta {K_{gx}} + \Delta {K_{gz}}} )}}{2}{\omega_U} + \delta_{gz}^X{\omega_r}\cos {\omega_r}t} \end{array}} \right]$$

From Eq. (21), $\delta _{gz}^X$, $\Delta {K_{gy}}$ and ${\varepsilon _{gy}}$ will be identified at the end of step 2. The analysis for accelerometer related parameters is similarly to the analysis for error parameters of gyroscope.

3.3 Observability analysis

During rotation, the state transition matrix is changing and the system is a time-varying system. But in small time interval, state transition matrix can be regarded as constant and such system is a PWCS. Therefore, observability analysis of time-varying system can be replaced by the analysis of PWCS [22–24].

The observability is obtained by calculating the rank of stripped observability matrix (SOM), expressed as Eq. (22).

(22)$${R_k} = {\left[ {\begin{array}{cccc} {{r_1}}&{{r_2}} & \cdots &{{r_k}} \end{array}} \right]^T},\,{r_j} = {\left[ {\begin{array}{cccc} H &{H{F_j}} & \cdots &{HF_j^{29 - 1}} \end{array}} \right]^T}$$

where j = 1, 2, …, k. H is shown in Eq. (17) and remains constant during calibration. F is shown in Eq. (14) and F_j is state transition matrix in jth interval. Only if the rank of R_k is 29, the system is completely observable [22]. The simulation frequency is 200 Hz. The rank at the end of each step is listed in Table 2. The rank increases to 29 after step 4.

Furthermore, SVD method is used to analysis observability [23,24]. A larger singular value means corresponding state variable has higher observability [24]. The singular value of SOM is defined as Eq. (23).

(23)$${R_k} = U\sum {V^T}\,\left( {\sum \textrm{ = }diag({{S_{r \times 1}},{0_{({n - r} )\times 1}}} )\;,S = diag({{\sigma_1},{\sigma_2}, \cdots {\sigma_r}} )} \right)$$

where R_k is a $m \times 29({m = 5 \times j} )$ matrix, ${\sigma _i}({i = 1,2, \cdots r} )$ is singular value, r is the rank of R_k. $U = [{{u_1},{u_2}, \cdots ,{u_m}} ]$ is a $m \times m$ matrix. $V = [{{v_1},{v_2}, \cdots ,{v_{29}}} ]$ is a $29 \times 29$ matrix.

Table 2. Rank of R_k.

View Table | View all tables in this article

For every singular value ${\sigma _i}$, if the kth value of ${{({u_i^T\tilde{Z}{v_i}} )} \mathord{\left/ {\vphantom {{({u_i^T\tilde{Z}{v_i}} )} {{\sigma_i}}}} \right.} {{\sigma _i}}}$ is largest, ${\sigma _i}$ is the singular value that corresponds to the kth state variable [24]. $\tilde{Z}$ is a $m \times 1$ matrix consists of measurements in Eq. (11). $\tilde{Z}$ is given in Eq. (24).

(24)$$\tilde{Z} = [{{z_1};{z_2}; \cdots ;{z_j}} ]\,\,{z_j} = [{\underbrace{{{Z_j};{Z_j}; \cdots ;{Z_j}}}_{{29 \; elements}}} ]$$

Table 3 lists singular value at the end of step 2, step 5 and step 8. Different steps increase observability for different error parameters, and all error parameters have significantly singular value in the end.

Table 3. Singular value for error parameters.

View Table | View all tables in this article

For example, during step 1 and step 2, the singular value of $\Delta {K_{gy}}$, $\Delta {K_{ax}}$ and $\Delta {K_{az}}$ are apparently larger than that of other scale factor errors. Then, IMU rotates along axis X_p, the singular value of $\Delta {K_{gx}}$, $\Delta {K_{ay}}$ and $\Delta {K_{az}}$ are increased. Finally, as the IMU rotates along axis Z_p, the observability of $\Delta {K_{gz}}$ is increased, almost equal to $\Delta {K_{gx}}$ and $\Delta {K_{gy}}$.

4. Simulation results and analysis

4.1 Self-calibration results

Simulation results are shown to verify the rotation scheme. The initial velocity and attitude are zero. The initial location is 116°E, 39°N. The simulation frequency is 200 Hz. IMU errors and alignment errors are listed in Table 4. The simulation is conducted fifty times. The standard deviation (STD) of fifty simulations are given in Table 5.

Table 4. Self-calibration simulation conditions.

View Table | View all tables in this article

Table 5. STD of fifty simulations.

View Table | View all tables in this article

The STD of installation errors and scale factor errors are less than 1 arc sec and 2 ppm. The STD of accelerometer biases is inside of 1.5ug. Thus, these errors are identified with high repeatability. Compared with their ideal value, the STD of gyroscope drifts are higher. But in RINS, because of modulation, estimation accuracy for gyroscope is a relatively low requirement.

According to the self-calibration results and observability analysis, the designed rotation scheme and the Kalman filter construction is verified. Therefore, with Kalman filter algorithm, the rotation scheme can be repeated to calibrate the temperature model.

4.2 Temperature self-calibration results

To obtain temperature model, rotation scheme mentioned in Section 3 is repeated ten times. For every error parameter, there will be ten different values at different temperature in one simulation. The initial velocity and attitude are zero. The initial location is 116°E, 39°N. The simulation frequency is 200Hz. The initial attitude errors and IMU random errors are listed in Table 4.

Temperature range is −20${\circ}{C}$ to +60${\circ}{C}$ with increasing rate is 12${\circ}{C}$/h. Six inertial sensors’ temperature make no difference. The temperature is described as Eq. (25).

(25) $$T ={-} 20.0 + 0.0034t$$

where t represents for time (sec) and T represents for temperature (${\circ}{C}$).

Temperature model is described by a second-order polynomial equation shown in Eq. (26).

(26)$$x = {a_0} + {a_1}T + {a_2}{T^2}$$

where a₀, a₁ and a₂ are coefficients, x is parameter listed in Table . The units of coefficients and the value of coefficients are shown in Table 6 to Table 9.

Table 6. Units of coefficients.

View Table | View all tables in this article

Table 7. Coefficients of drifts and biases

View Table | View all tables in this article

Table 8. Coefficients of scale factor errors.

View Table | View all tables in this article

Table 9. Coefficients of installation errors.

View Table | View all tables in this article

Figure 3 shows the residual error of ${\varepsilon _{gx}}$. Residual error still varies linearly with temperature, from −0.01deg/h to 0.01deg/h. In other words, the temperature model cannot reflect parameter’s real change under different temperature. The possible reason is estimated value isn’t equivalent to real value of mean temperature during one calibration procedure.

Fig. 3. Simulation results of x-gyro drift. Including the comparison between ideal value and estimated value, and the residual error of ${\varepsilon _{gx}}$.

Download Full Size | PDF

Iterative calibration method is used to correct temperature model. The iterative calculation does not need to conduct experiments repeatedly. The method aims at decreasing the residual errors for the same set of data. Thus, multiple iterations will not affect time performance seriously. After compensated with estimated parameters, residual error is estimated and old temperature model is modified. Repeating the process until residual error meets the demand. Figure 4 shows the change of estimated x-gyroscope drift and its residual error after iterative calibration.

Fig. 4. The change of estimated x-gyroscope drift and its residual error after iterative calibration.

Download Full Size | PDF

According to Fig. 4, iterative calibration decreases residual error effectively. After second iterative calibration, residual error of x-gyroscope drift declines to maximum of 1.2 × 10⁻³ deg/h. Moreover, computational costs are measured by MATLAB, as listed in Table 10.

Table 10. Computational costs of iterative calibration.

View Table | View all tables in this article

Maximum residual errors after four times iterative calibration are given in Table 11 and Table 12. Compared with Table 5, the residual errors after compensation are equivalent to the standard deviations of simulation results under constant temperature. Therefore, the temperature model is accuracy and repeatability of the model is acceptable.

Table 11. Maximum residual error of drifts, biases and scale factor errors.

View Table | View all tables in this article

Table 12. Maximum residual error of installation errors.

View Table | View all tables in this article

5. Experimental results and analysis

Taking advantages of temperature elevation since system power on, experiments are conducted to verify the method. The experiments are based on tri-axis RINS with FOG. IMU data was stored at 200Hz and processed offline. The location is 116.3468°E, 39.9788°N.

The RINS consists of FOG based IMU and three gimbals, including inner gimbal, middle gimbal and outer gimbal, as shown in Fig. 5. Each gimbal is driven by a torque motor and an angular encoder. The characteristics of RINS are given in Table 13.

Fig. 5. Experimental equipment.

Download Full Size | PDF

Table 13. Characteristics of RINS.

View Table | View all tables in this article

The experiments are conducted as follows: 1) system power on; 2) three gimbals back to initial position; 3) initial alignment; 4) gimbals rotate according to designed rotation scheme; 5) repeat 4) at least ten times; 6) system power off.

The experiments were conducted four times. Figure 6 shows part of calibration results. The variation with temperature for x-accelerometer bias, x-gyroscope scale factor and x-accelerometer is severer than installation errors and gyroscope drifts. Gyroscope drifts do not obviously present the trend changing with temperature. Since the calibration process doesn’t last for long time, the calibration accuracy for gyroscope drifts is lower than other error parameters, that is a possible explanation to gyroscope drifts’ calibration results. Moreover, the small range of temperature leads to the minor change for installation errors.

Fig. 6. Self-calibration results.

Download Full Size | PDF

Three temperature models are identified from the last three experiments. The average coefficients of three models are listed in Table 14 to Table 16. Since the error parameters are different at every time after system power on, a₀ is set as zero.

Table 14. Coefficients of drifts and biases.

View Table | View all tables in this article

Table 15. Coefficients of scale factor errors.

View Table | View all tables in this article

Table 16. Coefficients of installation errors.

View Table | View all tables in this article

After compensated with the average temperature model, the residual errors are shown in Fig. 7. Compared with Fig. 6, variation with temperature in Fig. 7 has been significantly reduced. STD of residual errors are listed in Table 17 and Table 18.

Fig. 7. Residual errors of temperature model.

Download Full Size | PDF

Table 17. STD of residual errors (drifts, biases and scale factor errors).

View Table | View all tables in this article

Table 18. STD of residual errors (installation errors).

View Table | View all tables in this article

From Table 17 and Table 18, the temperature models are more suitable for scale factor errors, installation errors and accelerometer biases. The STD of them are less than 4ppm, 1arc sec and 2ug, which means these models have higher repeatability and accuracy. However, the STD of gyroscope is inside of 0.005deg/h. A comparison of Fig. 6 and Fig. 7 indicates that the identified temperature model still works with gyroscope drifts, while lower accuracy estimation for gyroscope drifts is still a disadvantage to the method.

6. Conclusion

A self-calibration method for temperature errors in RINS is proposed in this paper. Simulations under variable temperature are performed, with the designed rotation scheme is repeated ten times. After temperature errors are reduced by the proposed method, the maximum residual errors are equivalent to the calibration errors under constant temperature. Therefore, the method is effective in decreasing temperature resulted errors for inertial sensors. Real implementations with FOG based RINS are carried out to demonstrate the method. Experimental results indicate that the identified temperature model works with all error parameters. Compared with traditional methods, the proposed method does not require temperature stabilization and calibration time is shortened.

However, the proposed method also has disadvantages. Since the rotation scheme does not last long, the method has a low accuracy for gyroscope drifts estimation. But for RINS, gyroscope drifts are modulated. Therefore, calibration precision requirement of gyroscope drifts is lower than installation errors and scale factor errors.

Funding

Aeronautical Science Foundation of China (NO. 20175851030).

Disclosures

The authors declare no conflicts of interest.

References

1. G. Chang, J. Xu, A. Li, and K. Cheng, “Error Analysis and Simulation of the Dual-Axis Rotation-Dwell Autocompensating Strapdown Inertial Navigation System,” in 2010 International Conference on Measuring Technology and Mechatronics Automation, (2010), pp. 124–127.

2. K. Li, Y. Chen, and L. Wang, “Online self-calibration research of single-axis rotational inertial navigation system,” Measurement 129, 633–641 (2018). [CrossRef]

3. X. Chen and C. Shen, “Study on temperature error processing technique for fiber optic gyroscope,” Optik 124(9), 784–792 (2013). [CrossRef]

4. C. Shen and X. Chen, “Analysis and modeling for fiber-optic gyroscope scale factor based on environment temperature,” Appl. Opt. 51(14), 2541–2547 (2012). [CrossRef]

5. X. Li, C. Zhang, Z. He, and Z. Zhong, “Temperature errors of IFOG and its compensation in engineering application,” in 2009 9th International Conference on Electronic Measurement & Instruments, (IEEE, 2009), 2-230-232-234.

6. C.-L. Zhou, Q. Zhang, S.-H. Yan, L. Gao, and G.-C. Wang, “Modeling and compensation for temperature errors of interferometric fiber optic gyroscope,” in 2009 International Conference on Information and Automation, (IEEE, 2009), 1443–1446.

7. A. M. Kurbatov and R. A. Kurbatov, “Temperature characteristics of fiber-optic gyroscope sensing coils,” J. Commun. Technol. Electron. 58(7), 745–752 (2013). [CrossRef]

8. Y. Zhang, Y. Guo, C. Li, Y. Wang, and Z. Wang, “A new open-loop fiber optic gyro error compensation method based on angular velocity error modeling,” Sensors 15(3), 4899–4912 (2015). [CrossRef]

9. Y. Pan, L. Li, C. Ren, and H. Luo, “Study on the compensation for a quartz accelerometer based on a wavelet neural network,” Meas. Sci. Technol. 21(10), 105202 (2010). [CrossRef]

10. J.-M. Gao, K.-B. Zhang, F.-B. Chen, and H.-B. Yang, “Temperature characteristics and error compensation for quartz flexible accelerometer,” Int. J. Autom. Comput. 12(5), 540–550 (2015). [CrossRef]

11. G. Wei, G. Li, Y. Wu, and X. Long, “Application of Least Squares-Support Vector Machine in system-level temperature compensation of ring laser gyroscope,” Measurement 44(10), 1898–1903 (2011). [CrossRef]

12. Q.-L. Xu, B. Han, and Z.-L. Deng, “Calibration Method for Laser Strapdown Inertial Measurement Unit In All Environment Temperature [J],” Journal of Chinese Inertial Technology6 (2004).

13. F. Jiao and Z. Chu, “An improved six-position hybrid calibration for RLG POS in full temperature,” in 2012 8th IEEE International Symposium on Instrumentation and Control Technology (ISICT) Proceedings, (IEEE, 2012), 246–250.

14. Y. Zhi-huai, Z. Wen-long, Z. Xiao-ya, and Z. Shu-ying, “Modeling and compensation method for input axis misalignment angle of FOG with continuous rotation,” Journal of Chinese Inertial Technology, 686–689 (2017).

15. Y. Mingye, S. Zhangqi, Z. Xueliang, and C. Yuzhong, “Temperature characteristics of input axis misalignment angle of fiber optic gyro,” J. of Nation. Univ. of Def. Techno. 36(3), 46–50 (2014). [CrossRef]

16. Q. Yang, R. Zhang, and H. Li, “Economical High(-)Low Temperature and Heading Rotation Test Method for the Evaluation and Optimization of the Temperature Control System for High-Precision Platform Inertial Navigation Systems,” Sensors 18(11), 3967 (2018). [CrossRef]

17. H. Yang, W. Huang, S. Jiao, X. Sun, W. Hong, and L. Qiao, “Temperature independent polarization-maintaining photonic crystal fiber with regular pentagon air hole distribution,” Optik 185, 390–396 (2019). [CrossRef]

18. Q. Ren, B. Wang, Z. Deng, and M. Fu, “A multi-position self-calibration method for dual-axis rotational inertial navigation system,” Sens. Actuators, A 219, 24–31 (2014). [CrossRef]

19. B. Liu, S. Wei, G. Su, J. Wang, and J. Lu, “An Improved Fast Self-Calibration Method for Hybrid Inertial Navigation System under Stationary Condition,” Sensors 18(5), 1303 (2018). [CrossRef]

20. P. Gao, K. Li, L. Wang, and Z. Liu, “A self-calibration method for tri-axis rotational inertial navigation system,” Meas. Sci. Technol.27(2016).

21. X. Wang and G. Shen, “A fast and accurate initial alignment method for strapdown inertial navigation system on stationary base,” J. Control Theory Appl. 3(2), 145–149 (2005). [CrossRef]

22. Q. Chong, C. Jiabin, H. Yongqiang, S. Chunlei, and W. Mingjie, “Online calibration method based on dual-axis rotation-modulating laser gyro SINS,” in 2016 Chinese Control and Decision Conference (CCDC), (IEEE, 2016), 3311–3315.

23. X.-X. Yang and Y.-M. Yin, “Discussions on observability and its applications in SINS,” Journal of Chinese Inertial Technology 20, 406–409 (2012).

24. J. Wu, G. Sun, R. Zhang, and Y. Zhang, “Analysis on observability of INS transfer alignment based on SVD method,” Journal of Chinese Inertial Technology 13, 26–30 (2005).

25. Q. Yongyuan, Inertial navigation, 1st ed. (Science Press, Beijing, China, 2006), pp. 249–258.

Error Parameters	Symbol	Error Parameters	Symbol
Gyroscope drifts	$[\begin{array}{ccc} ε_{g x} & ε_{g y} & ε_{g z} \end{array}]$	Accelerometer biases	$[\begin{array}{ccc} \nabla_{a x} & \nabla_{a y} & \nabla_{a z} \end{array}]$
Gyroscope scale factor errors	$[\begin{array}{ccc} Δ k_{g x} & Δ k_{g y} & Δ k_{g z} \end{array}]$	Accelerometer scale factor errors	$[\begin{array}{ccc} Δ k_{a x} & Δ k_{a y} & Δ k_{a z} \end{array}]$
Gyroscope installation errors	$[\begin{array}{ccccc} β_{g x}^{Y} & α_{g y}^{Z} & β_{g y}^{X} & δ_{g z}^{Y} & δ_{g z}^{X} \end{array}]$	Accelerometer installation errors	$[\begin{array}{cccc} α_{a x}^{Z} & α_{a y}^{Z} & δ_{a z}^{Y} & δ_{a z}^{X} \end{array}]$

Error parameters	Singular value
Error parameters	At the end of step 2	At the end of step 5	At the end of step 8
$[\begin{array}{ccc} ε_{g x} & ε_{g y} & ε_{g z} \end{array}]$	$[\begin{array}{ccc} 34 & 48 & 34 \end{array}]$	$[\begin{array}{ccc} 61 & 61 & 48 \end{array}]$	$[\begin{array}{ccc} 72 & 71 & 69 \end{array}]$
$[\begin{array}{ccc} \nabla_{a x} & \nabla_{a y} & \nabla_{a z} \end{array}]$	$[\begin{array}{ccc} 5 & 0 & 5 \end{array}]$	$[\begin{array}{ccc} 5 & 7 & 7 \end{array}]$	$[\begin{array}{ccc} 7 & 9 & 8 \end{array}]$
$[\begin{array}{ccc} Δ k_{g x} & Δ k_{g y} & Δ k_{g z} \end{array}]$	$[\begin{array}{ccc} 0 & 3 & 0 \end{array}]$	$[\begin{array}{ccc} 3 & 3 & 0 \end{array}]$	$[\begin{array}{ccc} 4 & 3 & 3 \end{array}]$
$[\begin{array}{ccc} Δ k_{a x} & Δ k_{a y} & Δ k_{a z} \end{array}]$	$[\begin{array}{ccc} 34 & 0 & 34 \end{array}]$	$[\begin{array}{ccc} 34 & 34 & 51 \end{array}]$	$[\begin{array}{ccc} 48 & 50 & 52 \end{array}]$
$[\begin{array}{ccccc} β_{g x}^{Y} & α_{g y}^{Z} & β_{g y}^{X} & δ_{g z}^{Y} & δ_{g z}^{X} \end{array}]$	$[\begin{array}{ccccc} 0 & 0 & 0 & 0 & 2 \end{array}]$	$[\begin{array}{ccccc} 1 & 2 & 1 & 2 & 2 \end{array}]$	$[\begin{array}{ccccc} 3 & 2 & 3 & 3 & 2 \end{array}]$
$[\begin{array}{cccc} α_{a x}^{Z} & α_{a y}^{Z} & δ_{a z}^{Y} & δ_{a z}^{X} \end{array}]$	$[\begin{array}{cccc} 0 & 34 & 21 & 0 \end{array}]$	$[\begin{array}{cccc} 34 & 34 & 23 & 34 \end{array}]$	$[\begin{array}{cccc} 49 & 43 & 43 & 49 \end{array}]$

Error parameters	Symbol	Ideal values	Unit
Gyroscope drifts	$[\begin{array}{ccc} ε_{g x} & ε_{g y} & ε_{g z} \end{array}]$	$[\begin{array}{ccc} 0.0 4 & 0.0 4 & 0.0 6 \end{array}]$	deg/h
Accelerometer biases	$[\begin{array}{ccc} \nabla_{a x} & \nabla_{a y} & \nabla_{a z} \end{array}]$	$[\begin{array}{ccc} 400 & 880 & 1760 \end{array}]$	ug
Gyroscope scale factor errors	$[\begin{array}{ccc} Δ k_{g x} & Δ k_{g y} & Δ k_{g z} \end{array}]$	$[\begin{array}{ccc} 90 & 45 & 100 \end{array}]$	ppm
Accelerometer scale factor errors	$[\begin{array}{ccc} Δ k_{a x} & Δ k_{a y} & Δ k_{a z} \end{array}]$	$[\begin{array}{ccc} 1500 & 1600 & 1700 \end{array}]$	ppm
Gyroscope installation errors	$[\begin{array}{ccccc} β_{g x}^{Y} & α_{g y}^{Z} & β_{g y}^{X} & δ_{g z}^{Y} & δ_{g z}^{X} \end{array}]$	$[\begin{array}{ccccc} 40 & 60 & 430 & 120 & 460 \end{array}]$	arc sec
Accelerometer installation errors	$[\begin{array}{cccc} α_{a x}^{Z} & α_{a y}^{Z} & δ_{a z}^{Y} & δ_{a z}^{X} \end{array}]$	$[\begin{array}{cccc} 100 & 480 & 60 & 250 \end{array}]$	arc sec
Gyroscope random drifts	$[\begin{array}{ccc} ε_{r x} & ε_{r y} & ε_{r z} \end{array}]$	$[\begin{array}{ccc} 0.0 1 & 0.01 & 0.0 1 \end{array}]$	deg/h
Accelerometer random biases	$[\begin{array}{ccc} \nabla_{r x} & \nabla_{r y} & \nabla_{r z} \end{array}]$	$[\begin{array}{ccc} 50 & 50 & 50 \end{array}]$	ug
Alignment errors	$[\begin{array}{ccc} Δ θ & Δ γ & Δ ψ \end{array}]$	$[\begin{array}{ccc} 10 & - 10 & 120 \end{array}]$	arc sec

Error Parameters	STD	Unit
$[\begin{array}{ccc} ε_{g x} & ε_{g y} & ε_{g z} \end{array}]$	$[\begin{array}{ccc} 0.0 021 & 0.0 014 & 0.0 015 \end{array}]$	deg/h
$[\begin{array}{ccc} \nabla_{a x} & \nabla_{a y} & \nabla_{a z} \end{array}]$	$[\begin{array}{ccc} 0 .34 & 1 .04 & 0 .92 \end{array}]$	ug
$[\begin{array}{ccc} Δ k_{g x} & Δ k_{g y} & Δ k_{g z} \end{array}]$	$[\begin{array}{ccc} 0 .34 & 0 .30 & 0 .64 \end{array}]$	ppm
$[\begin{array}{ccc} Δ k_{a x} & Δ k_{a y} & Δ k_{a z} \end{array}]$	$[\begin{array}{ccc} 1 .77 & 0 .61 & 0 .81 \end{array}]$	ppm
$[\begin{array}{ccccc} β_{g x}^{Y} & α_{g y}^{Z} & β_{g y}^{X} & δ_{g z}^{Y} & δ_{g z}^{X} \end{array}]$	$[\begin{array}{ccccc} 0 .18 & 0 .35 & 0 .59 & 0 .55 & 0 .25 \end{array}]$	arc sec
$[\begin{array}{cccc} α_{a x}^{Z} & α_{a y}^{Z} & δ_{a z}^{Y} & δ_{a z}^{X} \end{array}]$	$[\begin{array}{cccc} 0 .44 & 0 .64 & 0 .43 & 0 .60 \end{array}]$	arc sec

Coefficients	Gyroscope drifts	Accelerometer biases	Scale factor errors	Installation errors
$a_{0}$	deg/h	ug	ppm	arc sec
$a_{1}$	(deg/h)/ $^{\circ} C$	ug / $^{\circ} C$	ppm / $^{\circ} C$	(arc sec)/ $^{\circ} C$
$a_{2}$	(deg/h)/( $^{\circ} C$ )²	ug /( $^{\circ} C$ )²	ppm/( $^{\circ} C$ )²	(arc sec)/( $^{\circ} C$ )²

Self-calibration method for temperature errors in multi-axis rotational inertial navigation system

Abstract

1. Introduction

2. Coordinate system and error model

2.1 Coordinate system definition

2.2 IMU error model

3. Self-calibration scheme

3.1 Kalman filter construction

3.2 Design of rotation scheme

3.3 Observability analysis

4. Simulation results and analysis

4.1 Self-calibration results

4.2 Temperature self-calibration results

5. Experimental results and analysis

6. Conclusion

Funding

Disclosures

References

Cited By

Figures (7)

Tables (18)

Equations (26)

Optics Express

No.	Error parameters
No.	$[\begin{array}{ccccc} β_{g x}^{Y} & α_{g y}^{Z} & β_{g y}^{X} & δ_{g z}^{Y} & δ_{g z}^{X} \end{array}]$	$[\begin{array}{cccc} α_{a x}^{Z} & α_{a y}^{Z} & δ_{a z}^{Y} & δ_{a z}^{X} \end{array}]$
1	$[\begin{array}{ccccc} 0.31 & 0.24 & 0.37 & 0.013 & 0.080 \end{array}]$	$[\begin{array}{cccc} 0.34 & 0.86 & 0.51 & 0.55 \end{array}]$
2	$[\begin{array}{ccccc} 0.060 & 0.25 & 0.41 & 0.064 & 0.15 \end{array}]$	$[\begin{array}{cccc} 0.31 & 1.1 & 0.30 & 0.40 \end{array}]$
3	$[\begin{array}{ccccc} 0.21 & 0.17 & 0.25 & 0.044 & 0.24 \end{array}]$	$[\begin{array}{cccc} 0.24 & 1.0 & 0.42 & 0.48 \end{array}]$
4	$[\begin{array}{ccccc} 0.22 & 0.23 & 0.30 & 0.078 & 0.21 \end{array}]$	$[\begin{array}{cccc} 0.35 & 0.78 & 0.22 & 0.46 \end{array}]$
5	$[\begin{array}{ccccc} 0.13 & 0.26 & 0.29 & 0.19 & 0.25 \end{array}]$	$[\begin{array}{cccc} 0.48 & 1.0 & 0.11 & 0.49 \end{array}]$

Parameters		Value
Gyroscope	drift stability	≤0.05 deg/h
Gyroscope	scale factor repeatability	≤20 ppm
Accelerometer	bias repeatability	≤60 ug
Accelerometer	scale factor repeatability	≤100 ppm
Encoder	resolution	2.47 arc sec
Encoder	precision	≤36 arc sec

No.	Error parameters
No.	$[\begin{array}{ccccc} β_{g x}^{Y} & α_{g y}^{Z} & β_{g y}^{X} & δ_{g z}^{Y} & δ_{g z}^{X} \end{array}]$	$[\begin{array}{cccc} α_{a x}^{Z} & α_{a y}^{Z} & δ_{a z}^{Y} & δ_{a z}^{X} \end{array}]$
1	$[\begin{array}{ccccc} 0.48 & 0.81 & 0.65 & 0.14 & 0.34 \end{array}]$	$[\begin{array}{cccc} 0.79 & 0.26 & 0.38 & 0.51 \end{array}]$
2	$[\begin{array}{ccccc} 0.45 & 0.74 & 0.88 & 0.082 & 0.40 \end{array}]$	$[\begin{array}{cccc} 0.70 & 0.18 & 0.31 & 0.76 \end{array}]$
3	$[\begin{array}{ccccc} 0.39 & 0.67 & 0.54 & 0.14 & 0.16 \end{array}]$	$[\begin{array}{cccc} 0.53 & 0.13 & 0.35 & 0.39 \end{array}]$
4	$[\begin{array}{ccccc} 0.35 & 0.59 & 0.29 & 0.24 & 0.21 \end{array}]$	$[\begin{array}{cccc} 0.60 & 0.14 & 0.39 & 0.34 \end{array}]$

Step number	1	2	3	4	5	6	7	8
Rank(R_k)	24	26	26	29	29	29	29	29

Number of iterations	Computational cost (unit: sec)
1	604.48
2	1266.39
3	2123.80
4	3285.10

Step number	1	2	3	4	5	6	7	8
Rank(R_k)	24	26	26	29	29	29	29	29

Number of iterations	Computational cost (unit: sec)
1	604.48
2	1266.39
3	2123.80
4	3285.10