A New Norm-Observed Calibration Method Based on Improved Differential Evolution Algorithm for SINS

It is vital for a strapdown inertial navigation system (SINS) to be calibrated before normal use. In this paper, a new kind of norm-observed calibration method is proposed. Considering that the norm of the output of accelerometers and gyroscopes can be exactly the norm of local acceleration of gravity and Earth rotation angular velocity, respectively, optimization function about all-parameter calibration and the corresponding 24-position calibration path is established. Diﬀerential evolutionary algorithm (DE) is supposed to be the best option in parameter identiﬁcation due to its strong search and fast convergence abilities. However, the high-dimensional individual vector from calibration error equations restrains the algorithm’s optimum speed and accuracy. To overcome this drawback, improved DE (IDE) optimization is specially designed: First, current “DE/rand/1” and “DE/current-to-best/1” mutation strategies are combined as one with complementary advantages and overall balance during the whole optimization process. Next, with the increase of the evolutionary generation, the mutation factor can adjust itself according to the convergence situation. Multiple identiﬁcation tests prove that our IDE optimization has rapid convergence and high repeatability. Besides, certain motivation of external angular velocity is added to the gyroscope calibration, and a series of dynamic observation paths is formed, further improving the optimization accuracy. The ﬁnal static navigation experiment shows that SINS with calibration parameters solved by IDE has better performance over other identiﬁcation methods, which further explains that our novel method is more accurate and reliable in parameter identiﬁcation.


Introduction
Calibration is a key technology for inertial navigation systems to maintain high accuracy in the long term. Strapdown inertial navigation systems, once mounted on mobile carriers, should avoid repeated disassembly and keep high reliability for a long time in case of emergency usage. erefore, parameters downloaded inside the system need to be precisely identified and remain stable regardless of outside adverse circumstances, which brings a huge challenge for SINS calibration [1][2][3].
Separate calibration is a typical method and has been widely used for decades. Generally, the calibration process contains multiple static positions and continuous rotation tests with the aid of turntables. In each individual test, with the Earth's rotation rate and local gravity acceleration as a reference, a turntable can provide precise excitation of specific forces and angular rates for inertial measurement unit (IMU) [4][5][6]. However, the manufacturing error and installation error do exist in the actual platform, which causes considerable measure errors of specific forces and angular rates and will badly influence the parameter identification accuracy [7]. In [8], mounting errors between the IMU axes and the turntable axes are studied and the compensation is conducted successfully by using a thin-shell algorithm. In [9], a new calibration method based on D-suboptimal design is proposed in order to identify the IMU error model parameters and the turntable errors simultaneously. Although the above methods can restrain navigation errors caused by skew angles, it is rather difficult to eliminate them. Besides external factors that lead to the decrease of calibration accuracy, the effects of parameterfloat, nonlinearity, and coupling are also inevitable during the calibration process [10]. In [11], the accelerometer nonlinear scale factor is added to the error model and finally computed through the estimation of velocity errors. In [12], a general nonlinear model of the IMU output is derived, and online parameter estimation is performed well by means of a transformed unscented Kalman filter.
Later some scholar comes up with the idea of the normobserved calibration method [13]. Under the guidance that the norm of three accelerometer outputs is equal to local gravity and the norm of three gyroscope outputs is equal to the Earth's rotation rate in stationary condition [14], this method has been successfully applied and widely spread on many occasions. In [15], an autocalibration procedure for accelerometers is proposed. It requires no high-precision equipment and allows increasing the accuracy of this kind of sensor. In [16], a modified calibration method that takes advantage of norm observation is utilized to enhance the accuracy and performance of MEMS inertial sensor module. So far, most calibration models contain high-dimensional parameters and are usually computed through nonlinear optimization, which brings a great challenge to algorithm performance, such as convergence, stability, and other characteristics. In [5], the particle swarm optimization (PSO) algorithm is introduced in accelerometer calibration. e experiment shows it can successfully identify the accelerometer nonlinear scale factor and improve the multiposition calibration accuracy. In [17], the proposed PSO algorithm incorporates constriction factor along with time-varying inertial weight for improved results than the traditional PSO scheme. In [18], a novel calibration method by using the optimal artificial fish swarm algorithm is described in detail. Final results show that the system can be recalibrated without disassembly while ensuring accuracy and reliability. In [19], the genetic algorithm is utilized to obtain the inertial platform's optimal rotation trajectory based on the measurement of observability for all parameters.
With the rise of infield inertial sensor calibration, the optimization method needs to shorten the calculation time while ensuring accuracy. Among all the optimization algorithms, a differential evolutionary algorithm is found most suitable for the parameter identification processes, such as calibration of three-axis magnetometers, gyroscopes, and camera calibration. As an optimization heuristic algorithm that searches the global minimum or maximum, it can avoid local solutions of numerical problems. At the same time, it has the characteristics of fast convergence and strong robustness. Most importantly, the DE algorithm is proved not sensitive to initial parameters, which is an important advantage compared with traditional iteration algorithms. In [20], DE algorithm is used for the determination of the scale factor, misalignment, and bias of low-cost gyroscope L3GD20. However, the author in this paper only gave a rough plan and did not verify the final result. In [21], a calibration method using differential evolution is presented to determine the correction parameters of IMU. e final result shows the trend of convergence, and the introduction of root mean square error (RMSE) proves the reliability of DE algorithm identification in this area. But the running speed of the algorithm and the resource usage is not described in this article. With the aid of a nonmagnetic rotation platform, calibration of three-axis magnetometers with DE algorithm is analysed [22]. Compared with the unscented Kalman filter (UKF), recursive least squares (RLS), and genetic algorithm (GA), the DE algorithm has not only the least calibration error but also the best robustness. It is believed that for a three-axis sensor, the number of parameters to be identified is relatively not much. High-end inertial devices such as SINS, Attitude and Heading Reference System (AHRS) usually contain gyroscope, accelerometer, altimeter, magnetometer, and others, which greatly increases the type and number of sensors and complicate the identification process.
For this high-dimensional identification problem, in this paper, a novel norm-observed calibration method based on IDE optimization is proposed. In order to achieve engineering applications, a self-designed calibration scheme is introduced in detail and according to this, real inertial device data is collected and used to verify the performance of our optimization algorithm. After special upgrade and improvement, it is found that IDE greatly enhances the accuracy and repeatability of all-sensor parameter identification compared with the conventional ones. e final static navigation experiment also confirms that point.
is paper is organized as follows: In Section 2, sensor error models and calibration paths are established. Section 3 introduces optimization operation instructions and the algorithm improvement process is explained in detail. Experimental results of the proposed method and the traditional ones are provided in Section 4. Finally, conclusions are drawn in Section 5.

Inertial Sensors Error
Model. e symbols involved in this paper are defined as follows (the bold symbols in the formula are all represented as vectors): b: orthogonal body frame with right-front-up axis N g : three-axis gyroscope output in the nonorthogonal frame N a : three-axis accelerometer output in the nonorthogonal frame ω b : three-axis angular velocities input in the b frame f b : three-axis specific forces input in the b frame K g : diagonal matrix of the gyros' scale factors K a : diagonal matrix of the accelerometers' scale factors ε: three-axis gyroscope constant drift vector ∇: three-axis accelerometer constant drift vector E g : installation error matrix from the b frame to the gyro frame E a : installation error matrix from the b frame to the accelerometer frame e ideal error model of the gyros can be defined as follows [23]: (1) e matrix expansion of (1) is as follows: Similarly, the ideal error model of the accelerometers can be written as follows: e matrix expansion of (3) is as follows:

Principle of Norm-Observed Calibration and Its Path
Design. It is already known that the norm of SINS input angular velocities and specific forces are fixed values in static condition. According to (1) and (3), where g represents local gravitational acceleration and Ω is the Earth's rotation angular velocity. Practically, the calibration process contains measurement noises and identification errors. So, according to the actual situation, we have the following equations: Symbols with upper wavy lines represent the actual values. en the above objective optimization function can be established based on multiple static position observations where M represents the total number of static positions in the calibration path. e objective optimization functions L a and L g takes the form similar to the standard deviation and the optimal solution must be reached at the moment when the objective function is minimized.
In order to ensure that the parameters are fully stimulated, the optimal arrangement of static positions is also necessary. According to (2) and (4), each of them has 12 accelerometer feature parameters that need to be calibrated, that is, 3 scale factors, 6 installation errors, and 3 constant biases. us a 24-position calibration scheme is designed to completely recognize those parameters. e whole calibration path is divided into three phases and every phase contains 8 locating places.
right-handed coordinates and axis X b is fixed towards north (N). In position 1 (POS1), the axis Z b is adjusted towards up (U). Rotating 45°around axis X b from POS1, we can get POS2. en, do the same rotation around axis X b until reaching POS8. In each position, the turntable stops around 1 minute for SINS to collect enough data. All 8 locating places are illustrated in Figure 1

Calibration Parameter Identification Process Based on DE
Algorithm. e differential evolution algorithm shows great performance in solving optimization problems due to its simple structure, strong searching capability, and robustness [24,25]. e operation of this method normally contains 4 main stages, that is, initialization, mutation, crossover, and selection [26,27]. Counting on the algorithm's global search ability and rapid convergence characteristic, we decide to put DE into accurate identification of calibration parameters.
First, determine initial population size NP, evolutional generation EG, and objective function f to be maximized. Take the accelerometer calibration as an example, and there are 12 feature parameters that need to be calibrated, so the individual vector θ i is created:

Mathematical Problems in Engineering
For every element α j (j � 1, 2, . . . 12) in the vector θ i , its value should be randomly generated between the lower bound LB j and upper bound UB j .
According to the results acquired from the traditional separate calibration, we may set boundaries among different elements. See details in Table 1.
e population size is generally 3-8 times the dimension, so we set NP � 60. Other settings include EG � 1200 and objective function is set as follows: Secondly, in the mutation operation, "DE/rand/1" mutation strategy is initially selected [28]. And mutant vector at the next generation is generated by three different existing individuals: where θ x , θ y , θ z are distinct vectors uniformly chosen from the set θ 1 , θ 2 , . . . , θ NP , θ M is the mutant vector, t and t + 1 represent certain generation and its next, F is the mutation factor.
In this optimization process, the initial setting of F is set as 1.2, which normally ranges on the interval of (0, 2].
Next, in order to increase the diversity of the group, crossover operation is conducted by replacing some components α i j of θ i with the corresponding components α M j of θ M as follows: where α c j stands for elements in the vector θ C generated after the crossover operation, CR denotes crossover probability and is set as 0.6.
At last, to select the better one from the parent vector θ i (t) and the trial vector θ c (t + 1), the objective function is introduced to identify next-generation members: e latest generation members will be treated as parent vector and the mutation, crossover, and selection operation will start a new cycle until it reaches maximum generation or falls within preset convergence interval. In this way, an optimal solution will be obtained.

Improvement of Conventional DE Algorithm and Optimization Strategy of Key Parameters.
Follow the instructions in the previous subsection and perform 5 times optimization processes, and the convergence trend of the fitness function is shown in Figure 2.
It is clear to see that slow convergence in the early stage and poor repeatability among fitness values make the identified calibration parameters unconvincing. rough comprehensive analysis, we know that the whole calibration process needs to identify a total of 12 parameters, which means a high-dimensional optimization problem and will absolutely complicate the solving process. Besides, there are huge differences in the magnitude of parameters in Table 1. All these characteristics require the optimization algorithm must have strong traversal search ability and should not be too time-consuming for practical application. e chosen "DE/rand/1" strategy has very good global search capabilities [28], but local search is insufficient, and population diversity will become worse in the later stages of evolution. erefore, the current results are of poor repeatability and stability and it will take a lot more time to get a global optimum solution. "DE/current-to-best/1" is another frequently used mutation strategy and its basic configuration can be written as follows: where θ best (t) is the best vector at the current generation t.
Replace (11) with (14) and restart the optimization process, then we can get new convergence progress in Figure 3.
As is shown in Figure 3, "DE/current-to-best/1" usually has a faster convergence rate and performs well when solving unimodal problems. However, when solving multimodal problems, it is more likely to fall into a local optimum, leading to premature convergence. And the fitness function value obtained by this method is much smaller than the one solved by "DE/rand/1".
Since the above two schemes have their own advantages and disadvantages, they can be used in combination to achieve a nearly perfect equilibrium strategy. e integration mutation strategy we designed here is as follows: where b represents weight distribution factor and it is related to the current evolutional generation: Assuming that the initial value of b is 0, at this moment, the integration mutation strategy can be simplified as "DE/rand/1" in (11). With the increase of evolutionary generation, the value of b is also increasing. When b is equal to 1, the evolution strategy degenerates to "DE/current-to-best/1" in (14). In this way, the improved combination strategy can ensure the global ergodic feature in the early stage of evolution, and in the later stage of evolution, it can improve the convergence speed and the local search performance.
Put (15) into the mutation operation and restart the optimization process.
e new results are illustrated in Figure 4.
It can be seen that under the guidance of the new integration strategy, the final fitness value is a lot better than those in Figures 2 and 3. Besides, in the later stages of evolution, all of the results have converged to a stable state. We can transform the fitness function values into the objective optimization function in (7) and get the value range [0.007656, 0.008999], which shows that the final results have acceptable repeatability.
Back to (15), the mutation factor F remains a fixed value throughout the evolution process. Actually, the selection of this parameter has a great influence on the final evolutionary effect. It is generally believed that a smaller scaling factor F value is likely to cause the algorithm to converge locally, while a larger F value will reduce the convergence speed and even make the algorithm difficult to converge. erefore, F should have the ability to dynamically adjust itself during the entire evolution process: In the early stage of differential evolution algorithm, in order to avoid local optimization, the mutation factor should be larger. In the later stage, in order to retain elite individuals and improve search efficiency, the mutation factor should be smaller. Based on the above ideas, the adaptive design of the mutation factor is as follows:  where F max , F min indicate the upper and lower limits of the mutation factor, f max (t), f min (t) represent maximum and minimum values of fitness function at generation t, α, β are different weight factors. After many trials, finally, we make α � 0.1, β � 0.15. Now that mutation factor has the ability of selfadjustment. e combined mutation strategy can be further modified as follows:

Mathematical Problems in Engineering
Conduct 3 more optimization process and make a comparison with the results shown in Figure 4. e improvement effect after parameter dynamic adaptation is illustrated in Figure 5. e 3 new tests (test01, test02, and test03) are in use of the mutation strategy in (18). As is shown in Figure 5, the introduction of parameter dynamic adaptation further improves the identification accuracy. e corresponding function values in (7) of the new tests are in the range of [0.004162, 0.005808], which is smaller than those of the old tests. Final statistical results are listed in Table 2.

e Application and Correction of IDE Algorithm in Gyroscope Parameter Identification.
e improved DE algorithm has been successful in the identification of accelerometer parameters. Now IDE is applied in the parameter identification process of the gyroscope. What should be noted is that gyro data is collected from the same calibration path in Figure 1 as the accelerometer. e new individual vector φ i can be written as follows: e new boundary settings are shown in Table 3. e new fitness function is as follows: en following the same instructions as the accelerometer's identification and using the latest mutation strategy in (18), the gyroscope parameter identification results are illustrated in Figure 6.
With the same mutation strategy applied, the gyroscope's final fitness values are obviously smaller than the accelerometer's ones, which means the poor recognition process and substandard accuracy. e key reason for this phenomenon is the ill-fitted data selection. e calibration scheme is a kind of stationary method. During the 24-position static observations, the turntable stays still and the only input for gyros is the Earth's rotation angular velocity. It cannot create sufficient external excitation to be sensed. en the usage of static observation data causes      Mathematical Problems in Engineering e main purpose of the extra calibration scheme is to add some certain motivation of external angular velocity to the gyroscopes. As is shown in Figure 7, each axis is conducted with continuous positive and negative rotation. Take Z axis rotations (ROT1 and ROT2 in Figure 7) for example, the input of three-axis gyros ω x1 , ω y1 , ω z1 in positive rotation (ROT1) can be expressed as follows: where ω ie is the Earth's rotation angular velocity, ω is the turntable's rotation angular velocity, ϕ stands for local latitude, and t is the turntable rotation time. Substituting (21) into (2), we can get N gx1 � K gx ω ie cos ϕ sin ωt + K gx E gxy ω ie cos ϕ cos ωt + K gx E gxz ω + ω ie sin ϕ + ε x , N gy1 � K gy E gyx ω ie cos ϕ sin ωt + K gy ω ie cos ϕ cos ωt + K gy E gyz ω + ω ie sin ϕ + ε y , After the turntable turning n circles, take the sum of both sides in (22) and average the results: Similarly, during negative rotation (ROT2), we can get the following: Figure 7: Gyroscope supplementary calibration process.

Mathematical Problems in Engineering
Adding (23) and (24), then we obtain the following: Take the same steps for X and Y axis; then scale factors and installation errors of three-axis gyros can all be received. Due to the introduction of additional paths, the fitness function is no longer applicable and according to the derived solution in (25), the modified fitness function can be written as follows: where J represents the total number of symmetrical rotations, N + j , N − j stand for average outputs during positive and negative rotation along the same axis. ω b o is the angular velocity of external excitation provided by the turntable. In our calibration scheme, the norm of ω b o is set as 10°/s, 20°/s, and 30°/s, respectively. Use the improved fitness function to rerun the parameter optimization process. e convergence process is shown in Figure 8 and the final statistical results are listed in Table 4.
It can be seen from Table 4 that the high repeatability of the objective function occurs when the static calibration path is used. is is the result of falling into a local extreme. After adjusting the calibration path, the final convergence of the gyro parameter optimization is maintained at the same level as the accelerometer, which proves that the path selection and algorithm performance are equally important for accurate calibration.
What needs to be added is there are still constant biases waiting to be recognized. e already obtained scale factors and installation errors are now input parameters and substitute them into the static objective function in (7). At this time, the remaining optimization process only needs to identify the constant drifts of the three-axis gyroscope and is easy to complete.

Experimental Results and Algorithm
Performance Comparison

Experimental Condition Setting.
e system's configuration is shown in Figure 9(a). It should be noted that all the calibration paths of the above optimization tests are generated by the turntable in Figure 9(b) according to manual control instructions. e original calibration data of the accelerometers and gyroscopes are collected by the data acquisition device in Figure 9(c). Now in order to test navigation performance under different calibration parameters, a static navigation experiment is conducted under the same experimental environment and equipment.

Comparison of Experiment
Results. Since completely accurate calibration parameters are not available in advance, the accuracy of the identified parameters can be decided by the degree of divergence of the static navigation error. In this static navigation experiment, first preheat SINS long enough to eliminate changes of the ambient temperature on the effects of navigation accuracy. Next, download local latitude, longitude, and initial attitudes into the systems, respectively, and make preparations for navigation solution. During the formal navigation phase, raw data from inertial components are collected and several navigation solutions are operated inside with parameters gained from different identification methods.
To ensure the fairness of the comparative test, first, in the process of parameter optimization, the fitness functions used are all the same. e number of repeated optimization and iteration remains the same. Next, in each navigation solution with acquired calibration parameters, raw data are all measured and collected under the same calibration path we proposed in Subsections 2.2 and 3.3. Finally, statistical results of parameter identification are shown in Table 5 and navigation results are illustrated in Figure 10.
It is clear that the IDE method has the best performance in Table 5. Experiment results in Figure 10 show that the system calibrated by our proposed method IDE has the smallest velocity error divergence during the static navigation process: less than 2 m/s velocity error in the north and 10 m/s in the east. erefore, parameters gained from IDE optimization calibration are proved to be accurate and show better navigation performance.

Cost Effectiveness and Limitations of IDE Algorithm.
Now we discuss the time and space complexity of the algorithm as follows: Suppose the number of iterations is N, the dimension of the solution is M, and the population size is generally 3-8 times the dimension, taking 5 M. e standard differential evolution algorithm includes three steps of mutation, crossover, and selection. According to the number of times the basic sentence is executed, the time complexity of each step can be determined as follows: e overall time complexity of the traditional algorithm is as follows: As for our IDE method, due to the addition of the population optimal value calculation in the mutation process, the time complexity of each step of the algorithm becomes as follows: And time complexity of the new algorithm is as follows: In a word, the improved DE algorithm and the standard DE algorithm have the same time complexity, and there is only a constant difference in running time, which shows that the new method can be used efficiently.
In terms of space complexity, the storage space required during the execution of the algorithm includes three parts: the space occupied by the algorithm program, the storage space occupied by the input initial data, and the extra space required during the execution of the algorithm. Compared with the standard DE algorithm, the improved DE algorithm adds the determination of the best individual in group history and the adjustment process of the variation factor. erefore, the new method adds a certain storage unit.
It can be seen that the proposed IDE algorithm represents the compromise between accuracy and efficiency. e improved method sacrifices a certain convergence speed but in exchange for better optimization accuracy. However, compared with system-level calibration where the Kalman filter is used to estimate parameters, our algorithm has a slight advantage in time consumption. For further application in the field environment and fast calibration occasions, there are still improvements and modifications that need to be studied.  PS: GA represents genetic algorithm and PSO represents particle swarm optimization.

Conclusions
To overcome the impacts of turntable errors for traditional separate calibration, in this paper, a novel norm-observed calibration based on an improved differential evolution algorithm is introduced. First error models of inertial sensors are established and according to this, fitness functions for the identification of accelerometer and gyroscope parameters are built, respectively. Secondly, the 24-position calibration path is designed for traditional DE optimization. After the improvement of the integration mutation strategy and parameter dynamic adaptation, our IDE proves to be more efficient and accurate in finding a global optimum for the calibration problem. Aiming at the degeneration of the gyroscope parameter recognition, the analysis concludes that small motivation of external angular velocity occurs during static observation. So in the optimization process of gyroscope calibration, raw data is changed into a dynamic one and accordingly, its fitness function is modified. e final static navigation experiment shows that recalibrated parameters by IDE indeed improve the navigation performance. In a word, the norm-observed calibration method of SINS proposed in this paper is an effective way in parameter identification, which has great advancement in the SINS calibration field and much potential in future practical use, such as the temperature drift modeling and compensation, the alignment of different coordinate systems during camera shooting, and extraction and decomposition of various harmonic components in gravity measurement.

Data Availability
e data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this paper.