Optimized Multi-Position Calibration Method with Nonlinear Scale Factor for Inertial Measurement Units

Navigation grade inertial measurement units (IMUs) should be calibrated after Inertial Navigation Systems (INSs) are assembled and be re-calibrated after certain periods of time. The multi-position calibration methods with advantage of not requiring high-precision equipment are widely discussed. However, the existing multi-position calibration methods for IMU are based on the model of linear scale factors. To improve the precision of INS, the nonlinear scale factors should be calibrated accurately. This paper proposes an optimized multi-position calibration method with nonlinear scale factor for IMU, and the optimal calibration motion of IMU has been designed based on the analysis of sensitivity of the cost function to the calibration parameters. Besides, in order to improve the accuracy and robustness of the optimization, an estimation method on initial values is presented to solve the problem of setting initial values for iterative methods. Simulations and experiments show that the proposed method outperforms the calibration method without nonlinear scale factors. The navigation accuracy of INS can be improved by up to 17% in lab conditions and 12% in the moving vehicle experiment, respectively.


Introduction
Inertial Navigation System (INS) is an entirely self-contained system that solves positions of a point by integrating its accelerations [1]. It can provide high-rate attitude, velocity and position information, hence it is widely used as the navigation means of autonomous underwater vehicles, missiles, robots, aircrafts and other autonomous vehicles. Inertial Measurement Unit (IMU), composed of accelerometer and gyroscope, is the essential device of INS that plays a critical role in the precision of INS [2]. IMU calibration is a process of estimating the coefficients that convert the raw outputs of IMU to accelerations and rotation rates, which can be used to reflect the motion status of the vehicle. The applications of IMU calibration can be divided into two types: the factory calibration after the navigation system is assembled and the re-calibration process after certain intervals of time.
Traditionally, the navigation-grade IMU calibration is performed by comparing the IMU outputs with a known reference generated from the high-precision equipment [3]. Apparently, the accuracy of the traditional methods strongly relies on the accuracy of the calibration equipment [4]. Due to the frequent unavailability of costly high-precision calibration equipment, the multi-position calibration method has been widely discussed in recent years [5][6][7][8][9][10][11][12][13][14], which is based on the principle that the norms of outputs of the accelerometer and the gyroscope clusters are equal to the magnitudes of inputs of specific force and rotational velocity respectively. The multi-position calibration method can date back to 1995 when Ferraris F et al conducted their research in which the accelerometers and gyroscopes are calibrated by the local gravity and the geometrical quantities respectively without any rotation table or other velocity standards [5]. The IMU is placed on six static positions to determine the biases and scale factors, but unable to obtain the misalignment errors. Shin and Sheimy extended the methods by estimating the misalignment errors with gravity and earth rotation rate in 2002 [6]. However, the main drawback of this method is that the gyroscope scale factors and misalignment errors cannot be estimated reliably because of the observable problem that the magnitude of the earth rotation rate is very small. The method was modified in 2006 to calibrate the misalignment errors of gyroscopes with a single-axis turntable to provide a strong rotation rate signal by Newton's method [7]. David Jurman et al in 2007 proposed a method which employs constrained Newton optimization procedure for the estimation of 9 parameters including the bias, scale factor errors and misalignment errors [8]. However, this method suffers from the disadvantage of large computation and demand on precise initial conditions. Zhang et al in 2009 transformed the optimal problem to a set of linear equations and proposed a new approach to calibrate the inter-triad misalignment [9]. Although this method do not need any initial values like other iterative method the number of position and rotation clusters for IMU should be more than the number of equations in order to solve the set of equations. Cheuk et al in 2012 proposed a hand motion-based method to calibrate the consumer-grade IMU and utilized the accelerometer and magnetometer as the reference to calibrate the gyroscope in six minutes [10]. Yang et al in 2012 proposed an improved iterative estimation method to derive the scale factors, misalignments, biases and squared coefficients without any orientation information [11]. Cai et al in 2013 presented a calibration method for accelerometer with nonlinear scale factor using the particle swarm optimization to solve the nonlinear equation [12]. Although the results are very promising, the problem of initial values setting still exists. Särkkä et al in 2017 proposed an enhanced multi-position calibration method based on Gauss-Newton method for consumer-grade accelerometers, gyroscopes and magnetometers, and for accelerometers and magnetometers, the direction of reference signals, such as the gravity and the magnetic field of the Earth, are estimated with calibration parameters [13]. Wang et al in 2017 presented a 16-position calibration method for gyroscope's drift coefficients on centrifuge [14].
The calibration process is an optimized estimation for parameters. There are many methods for solving the following unconstrained optimization problem min x∈R n f (x) (1) where f (x) is a smooth cost function. Among them the line search methods [15][16][17] and trust region methods [18][19][20] are the most popular ones. Levenberg-Marquardt (LM) optimization with trust region method makes a good trade-off between the steepest decent method and the Gauss-Newton method which is widely used in nonlinear optimization [21,22]. However, the initial values should be approximately determined before estimation [23]. In a word, although the multi-position calibration method has been widely researched, some problems are still not settled. Firstly, the nonlinear scale factor of gyroscope and accelerometer should be calibrated accurately to improve the precision of INS. IMU has different scale factor for different specific force and angular velocity, due to the change of temperature and magnetic field. Besides, the scale factor often changes as a function of the products of specific force and angular velocity components [12]. The optimization methods used in the literature can be divided into 3 types: (1) Transform the set of nonlinear equations into linear equations; (2) Use the iterative methods and (3) Use the particle swarm optimization. However, the number of position and rotation for IMU should be increased to solve the linear equations in method (1). Rough initial values are needed to obtain global optimal values in method (2). The training process is cumbersome in method (3). Secondly, a simple optimization method without those disadvantages in method (1-3) should be investigated.
The main purpose of this paper is to improve the accuracy of IMU by estimating its nonlinear coefficients. The rest of this paper is organized as follows. The IMU models designed in this paper are Sensors 2019, 19, 3568 3 of 18 discussed in Section 2, including the nonlinear scale factors. The effects of nonlinear scale factors of IMU on the performance of INS are discussed in Section 3. The nonlinear optimization methods based on the optimal motions of IMU and the nonlinear optimization method with initial values estimation are discussed in Section 4. The analysis of simulation and experiment results are presented in Section 5. Finally, the conclusions are concluded in Section 6.

Definition of Frames
The frames used in the paper are listed in Table 1. Table 1. The definition of frames.

Symbol Frames
i The orthogonal inertial frame n The orthogonal navigation frame directs east-north-up (ENU) n The computer navigation frame e The earth-fixed frame r The turntable frame a The non-orthogonal frame denoted by accelerometers' sensitivity axes g The non-orthogonal frame denoted by gyroscopes' sensitivity axes p The orthogonal frame defined by a-frame q The orthogonal frame defined by g-frame A six-degree-of-freedom IMU composed of a triple axis accelerometer and a triple axis gyroscope. The accelerometer senses the acceleration along each input axis in a-frame, while gyroscope measures the angular velocity around each input axis in g-frame. Due to the assembly imperfections, each axis of IMU deviates by a small angle from its designed mounting direction. Hence, both a-frame and g-frame are the non-orthogonal frames. The p-frame is defined as follows: axis x p of p-frame coincides with the unit vector x a of a-frame, axis y p is perpendicular to x p in the plane x a y a , while z p , x p and y p together form the right-hand frame. The q-frame is defined as follows: axis x q of q-frame coincides with the unit vector x g of g-frame, axis y q is perpendicular to x q in the plane x g y g , while z q , x q and y q together form the right-hand frame. The definition of p-frame and q-frame are shown in Figure 1, and E ayx , E azx , E azy , E ayx , E azx , E azy are the misalignment errors of IMU. The transformation matrix C q p is defined as the inter-triad misalignment, which represent the transform form p-frame to q-frame. of IMU on the performance of INS are discussed in Section 3. The nonlinear optimization methods based on the optimal motions of IMU and the nonlinear optimization method with initial values estimation are discussed in Section 4. The analysis of simulation and experiment results are presented in Section 5. Finally, the conclusions are concluded in Section 6.

Definition of Frames
The frames used in the paper are listed in Table 1.

Symbol Frames i
The orthogonal inertial frame n The orthogonal navigation frame directs east-north-up (ENU) n′ The computer navigation frame e The earth-fixed frame r The turntable frame a The non-orthogonal frame denoted by accelerometers' sensitivity axes g The non-orthogonal frame denoted by gyroscopes' sensitivity axes p The orthogonal frame defined by a-frame q The orthogonal frame defined by g-frame A six-degree-of-freedom IMU composed of a triple axis accelerometer and a triple axis gyroscope. The accelerometer senses the acceleration along each input axis in a-frame, while gyroscope measures the angular velocity around each input axis in g-frame. Due to the assembly imperfections, each axis of IMU deviates by a small angle from its designed mounting direction. Hence, both a-frame and g-frame are the non-orthogonal frames. The p-frame is defined as follows: axis xp of p-frame coincides with the unit vector xa of a-frame, axis yp is perpendicular to xp in the plane xaya, while zp, xp and yp together form the right-hand frame. The q-frame is defined as follows: axis xq of q-frame coincides with the unit vector xg of g-frame, axis yq is perpendicular to xq in the plane xgyg, while zq, xq and yq together form the right-hand frame. The definition of p-frame and q-frame are shown in Figure 1, and Eayx, Eazx, Eazy, Eayx, Eazx, Eazy are the misalignment errors of IMU. The transformation matrix C q p is defined as the inter-triad misalignment, which represent the transform form p-frame to q-frame.

Nonlinear Model of Accelerometers
In this paper, the nonlinear model of accelerometers is established inspired by reference [12], and the accelerometer model is expressed as

Nonlinear Model of Accelerometers
In this paper, the nonlinear model of accelerometers is established inspired by reference [12], and the accelerometer model is expressed as where f p is the true value of the specific force in p-frame,f p andf a are the raw output of accelerometers in p-frame and a-frame, respectively; K a1 = diag(K ax1 , K ay1 , K az1 ), K a2 = diag(K ax2 , K ay2 , K az2 ) and K a3 = diag(K ax3 , K ay3 , K az3 ) are the linear, the second-order and the third-order scale factor of accelerometers, respectively; a and v p are the bias and the noise of accelerometers, respectively; N a is the output of accelerometers in pulses per second, N a(2) denotes the vector with the square of each element in N a and N a(3) denotes the vector with the cube of each element in N a : C p a is the transformation matrix from a-frame to p-frame, and can be written as where E ayx , E azx , E azy are the misalignment errors of accelerometers.

Nonlinear Model of Gyroscopes
Similarly, the gyroscope model can be expressed as where ω q iq is the true value of gyroscope output in q-frame, K g1 = diag(K gx1 , K gy1 , K gz1 ) and K g2 = diag(K gx2 , K gy2 , K gz2 ) are the linear and the second-order scale factor of gyroscopes, respectively; ε g and v g are the bias and the noise of gyroscopes, respectively; N g is the raw output of gyroscopes in pulses per second, N g(2) denotes the vector with the square of each element in N g : C q g is the transformation matrix from g-frame to q-frame, and can be written as where E gyx , E gzx , E gzy are the misalignment errors of gyroscopes.
The task of IMU calibration is to estimate the linear scale factors, nonlinear scale factors, biases and the misalignments of IMU. According to (2) and (6), the raw output of IMU can be calibrated in the orthogonal p-frame and q-frame from non-orthogonal a-frame and g-frame, respectively. In order to ensure the accuracy of alignment and navigation, it is necessary to calibrate the accelerometer and gyroscope in the same frame. Hence, the calibration method of the inter-triad misalignment in Reference [9] is utilized in this paper.

The Effects of Nonlinear Scale Factors
For high-precision navigation grade IMU with the precision of 0.01 mg and 0.01 • /h, the nonlinear scale factors should not be ignored. This section discuses about the effects of Nonlinear Scale Factors of IMU on the precision of INS. The well-known differential equations of navigation error for INS are described as [24]: . φ n = φ n × ω n ie + ω n en + δω n ie + δω n en − δω n ib (10) where V n = [V E V N V U ] T denotes the velocity of vehicle in n-frame, f n = [f E f N f U ] T denotes the accelerometers output in n-frame and φ = [φ E φ N φ U ] T denotes the misalignment angle between n-frame and n -frame. δ[·] denotes the error of the vector [·]. The rotation rate of the earth in n-frame is where L denotes the latitude of the vehicle, and the rotation rate of n-frame related to e-frame is where R M and R N denote the meridian radius of the earth. In order to analysis the effects of nonlinear scale factors of IMU only, substitute Equations (2), (6) into Equations (9), (10) and we have, . φ n = φ n × ω n ie + ω n en + δω n ie + δω n en − C n b K g2 N g(2) Hence, if the nonlinear scale factors are not calibrated, the velocity and attitude errors of INS are increased as shown in Equations (13) and (14). There exists a constant error in the differential equations caused by the nonlinear scale factors of IMU. Besides, when the vehicle is in the quasi-static state, the navigation error equations can be simplified as It can be seen that as time grows, the error caused by the nonlinear scale factors will be integrated into the velocity and attitude errors of INS. In summary, the nonlinear scale factors of IMU determines the growth speed of navigation errors. Hence, in order to improve the precision of INS, the nonlinear scale factors should be calibrated accurately.

Nonlinear Optimization of Accelerometer
Theoretically, the norm of accelerometer output in p-frame is equal to the specific force in n-frame, that is: where G is the gravity vector. According to (2) and (6), the nonlinear cost function of accelerometer can be expressed as: where N a i is the output of accelerometer in the ith position. In order to find the optimized position cluster and identify all parameters of accelerometer, the sensitivity of measurement to the parameters should be maximized. Define Take the partial derivatives of (19) with respect to the linear scale factor Ka 1 x, and ignore the little quantity items with misalignment errors: Similarly, take the partial derivatives of (19) with respect to other parameters: Based on Equations (20)- (25), the sensitivity of measurement to the parameters can be analyzed. Firstly, the sensitivity to linear scale factor, nonlinear scale factors and bias of the i-axis (i = x, y, z) accelerometer are maximized, when the sensitive direction of the i-axis (i = x, y, z)) accelerometer is parallel to the gravity vector. Secondly, the sensitivity to misalignment errors E aij (i = y, j = x or i = z, j = x, y) is maximized, when the gravity vector is in the plane formed by the i-axis and the j-axis accelerometer with an angle of 45 • or 135 • between the axis and the gravity vector. Hence, considering the sensitivity of measurement to parameters, the optimal position clusters to estimate the 15 parameters in (2) are shown in Figure 2.

Nonlinear Optimization of Gyroscope
The output of gyroscope contains the angular rate of SINS and the earth's rotation rate, that is: where ω q iq is the true value of gyroscope output, ω q ie is the earth's rotation rate, ω q r is the true value of the rotation rate that SINS is sensitive to. Different from the optimization of accelerometer, the earth's rotation rate is too small for gyroscope calibration of accurate scale factors and misalignment [6]. Hence, clockwise and counter-clockwise rotation of SINS is utilized to estimate the scale factors and misalignment of gyroscope firstly.

Nonlinear Optimization of Gyroscope
The output of gyroscope contains the angular rate of SINS and the earth's rotation rate, that is: Where q iq ω is the true value of gyroscope output, q ie ω is the earth's rotation rate, q r ω is the true value of the rotation rate that SINS is sensitive to. Different from the optimization of accelerometer, the earth's rotation rate is too small for gyroscope calibration of accurate scale factors and misalignment [6]. Hence, clockwise and counter-clockwise rotation of SINS is utilized to estimate the scale factors and misalignment of gyroscope firstly. Take (6) in (26), and integrate the angular rate over the clockwise rotation period t 1 , we have: Where N(τ) g+ is the output of gyroscope during clockwise rotation. Similarly, when SINS rotates in counter-clockwise during the period of t 2 : Where N(τ) gis the output of gyroscope during counter-clockwise rotation. Subtract (28) from (27), and make t= t 1 = t 2 : is the rotation angle of SINS in a period of t, hence the nonlinear cost function of gyroscope can be expressed as: Take (6) in (26), and integrate the angular rate over the clockwise rotation period t 1 , we have: where N(τ) g+ is the output of gyroscope during clockwise rotation. Similarly, when SINS rotates in counter-clockwise during the period of t 2 : where N(τ) g− is the output of gyroscope during counter-clockwise rotation. Subtract (28) from (27), and make t = t 1 = t 2 : where θ = t 0 ω q r (τ) dτ is the rotation angle of SINS in a period of t, hence the nonlinear cost function of gyroscope can be expressed as: Similar to accelerometer, in order to find the optimized rotation cluster and then identify the scale factors and the misalignment errors of gyroscope, the sensitivity of measurement to the parameters is maximized. Define Take partial derivatives of (31) with respect to the scale factors and misalignment errors: According to equations (32)-(34), the sensitivity to scale factors of the i-axis (i = x, y, z) gyroscope are maximized, when the sensitive direction of the i-axis (i = x, y, z) gyroscope is parallel to the rotation vector. Besides, the sensitivity to misalignment errors E gij (i = y, j = x or i = z, j = x, y) is maximized, when the rotation vector is in the plane formed by the i-axis and the j-axis gyroscope with an angle of 45 • or 135 • between the axis and the rotation vector. Hence, the optimization of rotation clusters is to estimate the 9 parameters: the scale factors and misalignment errors in (6) are shown in Figure 3. Similar to accelerometer, in order to find the optimized rotation cluster and then identify the scale factors and the misalignment errors of gyroscope, the sensitivity of measurement to the parameters is maximized. Define Where Nsum=[Nsumx, Nsumy, Nsumz] T , and define Take partial derivatives of (31) with respect to the scale factors and misalignment errors: According to equations (32)-(34), the sensitivity to scale factors of the i-axis (i=x, y, z) gyroscope are maximized, when the sensitive direction of the i-axis (i=x, y, z) gyroscope is parallel to the rotation vector. Besides, the sensitivity to misalignment errors Egij (i=y, j=x or i=z, j=x, y) is maximized, when the rotation vector is in the plane formed by the i-axis and the j-axis gyroscope with an angle of 45° or 135° between the axis and the rotation vector. Hence, the optimization of rotation clusters is to estimate the 9 parameters: the scale factors and misalignment errors in (6) are shown in Figure 3. The bias of gyroscope can be estimated using the position clusters redundantly shown in Figure 2, with the cost function:  The bias of gyroscope can be estimated using the position clusters redundantly shown in Figure 2, with the cost function: where C q g , K g1 and K g2 are estimated in (30), N g i is the output of gyroscope in the ith position.

The Nonlinear Optimization with Initial Values Estimation
The nonlinear unconstrained optimization problems described in Equations (18) f i (x) 2 (36) The Levenberg-Marquardt (LM) algorithm, one of the most popular algorithms for the solution of nonlinear least squares problems, is used in this paper. To implement LM algorithm, Ceres Solver, an open source C++ library to handle large complex optimization problems, is used. For cost function f (x) that is strongly convex and twice differentiable, the iterative sequence using LM algorithm will be where x (i) (i = k, k+1) is the parameters vector at the ith iteration, H is the Hessian matrix of f (x), J is the Jacobian matrix of f (x) and β is the blending factor that determines the mix between steepest descent and Newton-Raphson [25].
where K 0 a1 and K 0 g1 are the rough linear scale factor matrix of accelerometer and gyroscope respectively. Hence, the initial values of linear scale factors can be approximately determined.

Analysis of Simulation Results
Monte Carlo simulations are conducted to verify the proposed method. Using MATLAB, the IMU outputs with random errors are generated in a uniform distribution as the true values. The calibration is conducted using the proposed method, and the calibration errors can be calculated.
The simulation conditions are set as: the turntable angle errors are 3", 3 and 3 • (1σ), respectively to verify the relationship between the accuracy of proposed method and the accuracy of turntable. The random bias of accelerometer and gyroscope are 100 µg (1σ) and 0.  Table 2. It should be pointed out that the distribution of all parameters shown in Table 2 are set based on the IMU parameters in real experiments.  Table 2 shows that: firstly, the mean values and the root mean squares of calibration errors in Monte Calo experiments are far less than the calibration parameters which means that all calibration parameters of accelerometers can be calibrated accurately. The maximum calibration errors of bias, linear scale factor, second-order scale factor, third-order scale factor and misalignment error are 0.045 µg, −7.49 × 10 −7 µg/pulse, 2.55 × 10 −3 f g/pulse 2 , 1.96 × 10 −2 zg/pulse 3 and −6.00 × 10 −4 , respectively.
Secondly, whatever the turntable errors are, the calibration results of the proposed method are almost the same, and the calibration errors of parameters are of the same order of magnitude. Hence, the turntable angle errors have no effects on the calibration errors of proposed method, which means that the proposed method can calibrate all parameters of accelerometers, including the nonlinear scale factors, without relying on the error of the turntable.
The statistical results of 500 Monte Carlo simulations of gyroscopes are shown in Table 3. Similar to accelerometers, all calibration parameters of gyroscopes can be calibrated accurately without relying on the angle errors of the turntable. The maximum calibration errors of bias, linear scale factor, second-order scale factor and misalignment error are 9.88 × 10 −10 • /s, −3.38 × 10 −16 • /s/pulse, −6.21 × 10 −20 • /s/pulse 2 and 7.95 × 10 −10 , respectively. The above simulation results verified the correctness of establishing calibration model of IMU with nonlinear scale factor. The calibration results shown in Tables 2 and 3   The calibration results shown in Table 2 and Table 3 are with the initial values estimation proposed in 4.3. To verify that the proposed method can help the LM algorithm converge to the global optimal values, 500 Monte Carlo simulations are carried out, and the calibration errors without and with the initial values estimation of gyroscopes are shown in Figure 4 and Figure 5 respectively.   The calibration results shown in Table 2 and Table 3 are with the initial values estimation proposed in 4.3. To verify that the proposed method can help the LM algorithm converge to the global optimal values, 500 Monte Carlo simulations are carried out, and the calibration errors without and with the initial values estimation of gyroscopes are shown in Figure 4 and Figure 5 respectively.   As shown in Figure 4, the 34th and the 469th Monte Carlo simulation converge to the local optimal values leading to the unacceptable calibration errors because of the incorrect initial values of LM algorithm. Otherwise, the calibration errors with the initial values estimation proposed in 4.3 are acceptable as shown in Figure 5. Figures 4 and 5 shows that estimating the linear scale factors of IMU firstly and make them as the initial values of the optimization can improve the accuracy and the robustness of calibration. The above simulation results showed that the proposed calibration method not only efficiently identified the nonlinear scale factors of IMU without relying on the accuracy of the turntable, but also improved the accuracy and robustness of the calibration with the initial values estimation.

Analysis of Experiment Results
In order to verify the proposed calibration method in practice, the calibration experiments, the la navigation experiments and the moving vehicle experiments are carried out based on 1-axis Rotational Inertial Navigation System (RINS). The RINS consists of three fiber optic gyroscopes, three quartz accelerometers with the precision of 0.01 • /h and 0.01 mg respectively, a rotating mechanism with the rotation axis pointing to vertical direction and a core control processor based on Digital Signal Processor (DSP).

Calibration Experiments results
The calibration experiments are conducted using three methods: (1) The Traditional calibration method based on the high-precision 3-axis turntable shown in Figure 6 with about 5"angle errors, whose accuracy relies on the accuracy of the turntable.
(2) The Multi-position calibration method in Reference [9] based on the low-precision 2-axis turntable shown in Figure 7 with about 5 angle errors and the rotating mechanism of RINS.
(3) The proposed method in this paper based on the turntable in method 2 and the rotating mechanism of RINS.
It should be pointed out that the rotating mechanism of RINS is used to help the RINS complete the rotation shown in Figure 3 on the 2-axis turntable.
acceptable as shown in Figure 5. Figure 4 and Figure 5 shows that estimating the linear scale factors of IMU firstly and make them as the initial values of the optimization can improve the accuracy and the robustness of calibration. The above simulation results showed that the proposed calibration method not only efficiently identified the nonlinear scale factors of IMU without relying on the accuracy of the turntable, but also improved the accuracy and robustness of the calibration with the initial values estimation.

Analysis of Experiment Results
In order to verify the proposed calibration method in practice, the calibration experiments, the la navigation experiments and the moving vehicle experiments are carried out based on 1-axis Rotational Inertial Navigation System (RINS). The RINS consists of three fiber optic gyroscopes, three quartz accelerometers with the precision of 0.01 °/h and 0.01 mg respectively, a rotating mechanism with the rotation axis pointing to vertical direction and a core control processor based on Digital Signal Processor (DSP).

Calibration Experiments results
The calibration experiments are conducted using three methods: (1) The Traditional calibration method based on the high-precision 3-axis turntable shown in Figure 6 with about 5″angle errors, whose accuracy relies on the accuracy of the turntable.
(2) The Multi-position calibration method in Reference [9] based on the low-precision 2-axis turntable shown in Figure 7 with about 5′angle errors and the rotating mechanism of RINS.
(3) The proposed method in this paper based on the turntable in method 2 and the rotating mechanism of RINS.
It should be pointed out that the rotating mechanism of RINS is used to help the RINS complete the rotation shown in Figure 3 on the 2-axis turntable.  Method 1 calibrates the IMU in the frame of turntable r-frame, while method 2 and 3 firstly calibrates the accelerometers and gyroscopes in p-frame and q-frame respectively, and then calibrates the inter-triad misalignment, making the accelerometers and gyroscopes calibrated to the same orthogonal frame q-frame [9]. The inter-triad misalignment C q p , the transformation matrix C r a and r can be expressed by Euler angles α1, β1, γ1, α2, β2, γ2, α3, β3 and γ3, respectively. Method 1 calibrates the IMU in the frame of turntable r-frame, while method 2 and 3 firstly calibrates the accelerometers and gyroscopes in p-frame and q-frame respectively, and then calibrates the inter-triad misalignment, making the accelerometers and gyroscopes calibrated to the same orthogonal frame q-frame [9]. The inter-triad misalignment C q p , the transformation matrix C r a and C r g can be expressed by Euler angles α 1 , β 1 , γ 1 , α 2 , β 2 , γ 2 , α 3 , β 3 and γ 3 , respectively.
The calibration results of method 1, 2 and 3 for RINS are shown in Table 4. It is obvious that method 3 proposed in the paper can calibrate the nonlinear scale factors of IMU compared with Method 2. Method 1 calibrates the IMU by the comparison of the IMU outputs with the turntable outputs, hence the results of calibration rely on the precision of the turntable. However, the multi-position calibration methods calibrate the IMU in a reference frame instead of the turntable frame. Therefore, the precision of the turntable has no effects on the calibration results.

Navigation Experiments in lab results
Install the RINS on the 3-axis turntable shown in Figure 6, and perform the self-alignment and navigation experiments in 2 states: (1) The quasi-static state that keeps the turntable in a fixed angle; (2) The swing state that enables the turntable swing along 3 axes in the condition shown in Table 5, where Heading, Pitch, Roll denotes the rotation axis of the turntable.

Heading Pitch Roll
Swing frequency (Hz) 2 2 8 Swing amplitude ( • ) 3 3 5 The position errors in state (1) and (2) are shown in Figures 8 and 9 respectively. Method 1, 2 and 3 are the calibration methods described in Section 5.2.1. It is obvious that Method 1 and Method 2 have the similar accuracy on the position errors, while Method 3 proposed in this paper leads to higher precision of navigation for RINS.

Navigation Experiments in lab results
Install the RINS on the 3-axis turntable shown in Figure 6, and perform the self-alignment and navigation experiments in 2 states: (1) The quasi-static state that keeps the turntable in a fixed angle; (2) The swing state that enables the turntable swing along 3 axes in the condition shown in Table 5, where Heading, Pitch, Roll denotes the rotation axis of the turntable. The position errors in state (1) and (2) are shown in Figure 8 and Figure 9 respectively. Method 1, 2 and 3 are the calibration methods described in 4.2.1. It is obvious that Method 1 and Method 2 have the similar accuracy on the position errors, while Method 3 proposed in this paper leads to higher precision of navigation for RINS.   Install the RINS on the 3-axis turntable shown in Figure 6, and perform the self-alignment and navigation experiments in 2 states: (1) The quasi-static state that keeps the turntable in a fixed angle; (2) The swing state that enables the turntable swing along 3 axes in the condition shown in Table 5, where Heading, Pitch, Roll denotes the rotation axis of the turntable. The position errors in state (1) and (2) are shown in Figure 8 and Figure 9 respectively. Method 1, 2 and 3 are the calibration methods described in 4.2.1. It is obvious that Method 1 and Method 2 have the similar accuracy on the position errors, while Method 3 proposed in this paper leads to higher precision of navigation for RINS.    (16). Besides, it is concluded that compared with the navigation accuracy under quasi-static conditions, the navigation accuracy under the dynamic conditions can be more improved by the proposed method. It is because that the gyroscope only senses the earth rotation rate in quasi-static state, while it also senses the rotation rate shown in Table 5 in swing state. Hence, the non-linearity of gyroscope's scale factor is more significant in swing state than in quasi-static stare, and it contributes more to the position error in swing state. the dynamic conditions can be more improved by the proposed method. It is because that the gyroscope only senses the earth rotation rate in quasi-static state, while it also senses the rotation rate shown in Table 5 in swing state. Hence, the non-linearity of gyroscope's scale factor is more significant in swing state than in quasi-static stare, and it contributes more to the position error in swing state.

Moving Vehicle Experiment Results
To verify the calibration method based on the navigation errors of RINS, the moving vehicle experiments are carried out 3 times using the same calibration results to compare the results of three different calibration methods. The navigation experiment vehicle is shown in Figure 11, which is a human operated automobile equipped with a GPS receiver, 1-axis RINS and a computer for data visualization. The 1-axis RINS is installed inside the vehicle as shown in Figure 12. The precision of GPS is 1m, as the reference for navigation. The trajectory of the vehicle is shown in Figure 13, and the route includes movements of turning, uphill, downhill, acceleration and deceleration within 3.2 hours.

Moving Vehicle Experiment Results
To verify the calibration method based on the navigation errors of RINS, the moving vehicle experiments are carried out 3 times using the same calibration results to compare the results of three different calibration methods. The navigation experiment vehicle is shown in Figure 11, which is a human operated automobile equipped with a GPS receiver, 1-axis RINS and a computer for data visualization. The 1-axis RINS is installed inside the vehicle as shown in Figure 12. The precision of GPS is 1m, as the reference for navigation. The trajectory of the vehicle is shown in Figure 13, and the route includes movements of turning, uphill, downhill, acceleration and deceleration within 3.2 hours. the dynamic conditions can be more improved by the proposed method. It is because that the gyroscope only senses the earth rotation rate in quasi-static state, while it also senses the rotation rate shown in Table 5 in swing state. Hence, the non-linearity of gyroscope's scale factor is more significant in swing state than in quasi-static stare, and it contributes more to the position error in swing state.

Moving Vehicle Experiment Results
To verify the calibration method based on the navigation errors of RINS, the moving vehicle experiments are carried out 3 times using the same calibration results to compare the results of three different calibration methods. The navigation experiment vehicle is shown in Figure 11, which is a human operated automobile equipped with a GPS receiver, 1-axis RINS and a computer for data visualization. The 1-axis RINS is installed inside the vehicle as shown in Figure 12. The precision of GPS is 1m, as the reference for navigation. The trajectory of the vehicle is shown in Figure 13, and the route includes movements of turning, uphill, downhill, acceleration and deceleration within 3.2 hours.  The maximum position errors of the three field experiments are shown in Table 6. Because the results of the three experiments are similar, we take the first experiment result for example for detailed analysis in this paper. The trajectories of the GPS output and the RINS navigation result using the parameters of methods 1-3 are shown in Figure 13. The position errors of three RINS navigation results are compared with the GPS output in Figure 14. The field experiment results show The maximum position errors of the three field experiments are shown in Table 6. Because the results of the three experiments are similar, we take the first experiment result for example for detailed analysis in this paper. The trajectories of the GPS output and the RINS navigation result using the parameters of methods 1-3 are shown in Figure 13. The position errors of three RINS navigation results are compared with the GPS output in Figure 14. The field experiment results show that the precision of method 1 is approximately equal to that of method 2. As the nonlinear scale factors can be accurately calibrated, the navigation results using the parameters of method 3 are better than both method 1 and method 2. As shown in Table 6, the position error of method 3 in moving vehicle navigation experiments can be decreased by 12.19% maximally compared to that of method 1. Besides, the results of three field experiments show that the maximum position error can be reduced by an average of 11% with the calibration and compensation of nonlinear scale factor of IMU. Different from the lab experiments, the accelerometer senses the acceleration of the vehicle in addition to the gravity acceleration in the field experiments. Therefore, the non-linearity of accelerometer's scale factor is more significant in field experiments than that in lab experiments. Hence, similar to the results of navigation experiments in lab conditions, it can be concluded that the position precision is also improved in field condition using the proposed method, which can estimate the IMU nonlinear scale factors accurately without high-precision turntable.   2 The decreased percentage of maximum position error between method 3 and method 1 Figure 13. Trajectories of the navigation results.

Conclusions
This paper presents a study on the optimized calibration method with nonlinear scale factor for IMU. The effects of nonlinear scale factors of IMU are analyzed, and it proved that the nonlinear scale factors should not be ignored in order to improve the accuracy of navigation for high-precision INS. A nonlinear optimization model of IMU is established, and the optimized calibration motion of IMU is designed based on the analysis of sensitivity of the cost function to the calibration parameters. To solve the nonlinear optimization problems and obtain the global optimal values, LM algorithm of Ceres Solver is used, and in addition, the model for estimating the initial values of nonlinear optimization is established to improve the accuracy and robustness of the optimization. Finally, simulations and experiments are conducted to test the performance of the proposed method. The results of navigation experiments based on the traditional calibration method, the multi-position

Conclusions
This paper presents a study on the optimized calibration method with nonlinear scale factor for IMU. The effects of nonlinear scale factors of IMU are analyzed, and it proved that the nonlinear scale factors should not be ignored in order to improve the accuracy of navigation for high-precision INS. A nonlinear optimization model of IMU is established, and the optimized calibration motion of IMU is designed based on the analysis of sensitivity of the cost function to the calibration parameters. To solve the nonlinear optimization problems and obtain the global optimal values, LM algorithm of Ceres Solver is used, and in addition, the model for estimating the initial values of nonlinear optimization is established to improve the accuracy and robustness of the optimization. Finally, simulations and experiments are conducted to test the performance of the proposed method. The results of navigation experiments based on the traditional calibration method, the multi-position calibration method without the nonlinear scale factors and the proposed calibration method are compared. It is shown that in the calibration of nonlinear scale factor for IMU without high-precision turntable, the position precisions can be improved by up to 17% in the lab conditions and 12% in the moving vehicle experiment respectively. It is concluded that the traditional calibration method and the multi-position calibration method without the nonlinear scale factors have the similar accuracy, while the proposed method with the nonlinear scale factors leads to higher precision of navigation for INS.
Author Contributions: X.C. and Z.W. proposed the initial idea and conceived the experiments. Z.W. performed the experiments and wrote the paper. X.C. and J.F. reviewed and edited the manuscript. All authors read and approved this manuscript.
Funding: This research was funded by the National Natural Science Foundation of China (No. 61773116).

Conflicts of Interest:
The authors declare no conflict of interest.