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

Authors: Yang Liu; Gongliu Yang; Qingzhong Cai; Lifen Wang
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

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. Differential evolutionary algorithm (DE) is supposed to be the best option in parameter identification 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 identification 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 final static navigation experiment shows that SINS with calibration parameters solved by IDE has better performance over other identification methods, which further explains that our novel method is more accurate and reliable in parameter identification.


---

## Body

## 1. 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. Therefore, 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–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–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 parameter-float, 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 norm-observed 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. The 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. The 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. The final static navigation experiment also confirms that point.This paper is organized as follows: In Section2, 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.


## 2. Error Model of Inertial Sensors and SINS Calibration Path Design

### 2.1. Inertial Sensors Error Model

The 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 axisNg: three-axis gyroscope output in the nonorthogonal frameNa: three-axis accelerometer output in the nonorthogonal frameωb: three-axis angular velocities input in the b framefb: three-axis specific forces input in the b frameKg: diagonal matrix of the gyros’ scale factorsKa: diagonal matrix of the accelerometers’ scale factorsε: three-axis gyroscope constant drift vector∇: three-axis accelerometer constant drift vectorEg: installation error matrix from the b frame to the gyro frameEa: installation error matrix from the b frame to the accelerometer frameThe ideal error model of the gyros can be defined as follows [23]:(1)Ng=KgEgωb+ε.The matrix expansion of (1) is as follows:(2)NgxNgyNgz=Kgx000Kgy000Kgz1EgxyEgxEgyx1EgyzEgzxEgzy1ωxbωybωzb+εxεyεz.Similarly, the ideal error model of the accelerometers can be written as follows:(3)Na=KaEafb+∇.The matrix expansion of (3) is as follows:(4)NaxNayNaz=Kax000Kay000Kaz1EaxyEaxzEayx1EayzEaxxEazy1fxbfybfzb+∇x∇y∇z.


### 2.2. 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),(5)fb2=KaEa−1Na−∇TKaEa−1Na−∇=g2,ωb2=KgEg−1Ng−εTKgEg−1Ng−ε=Ω2,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:(6)f˜b2=K˜aE˜a−1Na−∇˜TK˜aE˜a−1Na−∇˜,ω˜b2=K˜gE˜g−1Ng−ε˜TK˜gE˜g−1Ng−ε˜.Symbols with upper wavy lines represent the actual values. Then the above objective optimization function can be established based on multiple static position observations(7)La=∑m=1Mg2−f˜bm22M−1,Lg=∑m=1MΩ2−ω˜bm22M−1,where M represents the total number of static positions in the calibration path.The objective optimization functionsLa and Lg 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. Thus a 24-position calibration scheme is designed to completely recognize those parameters.The whole calibration path is divided into three phases and every phase contains 8 locating places.Phase 1: assume the SINS’s body frame is theO−XbYbZb right-handed coordinates and axis Xb is fixed towards north (N). In position 1 (POS1), the axis Zb is adjusted towards up (U). Rotating 45° around axis Xb from POS1, we can get POS2. Then, do the same rotation around axis Xb 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(a).Phase 2: AxisYb is fixed towards north. Follow the same rotation operation as Phase 1 around axis Yb 8 times and get 8 positions. Also, keep standing time in each position. The rotation scheme of Phase 2 is shown in Figure 1(b).Phase 3: AxisZb is fixed towards north. Follow the same rotation operation as Phase 1 around axis Zb 8 times and get 8 positions. Also, keep standing time in each position. The rotation scheme of Phase 3 is shown in Figure 1(c).Figure 1
24-position calibration path. (a) Phase 1. (b) Phase 2. (c) Phase 3.
(a)(b)(c)


## 2.1. Inertial Sensors Error Model

The 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 axisNg: three-axis gyroscope output in the nonorthogonal frameNa: three-axis accelerometer output in the nonorthogonal frameωb: three-axis angular velocities input in the b framefb: three-axis specific forces input in the b frameKg: diagonal matrix of the gyros’ scale factorsKa: diagonal matrix of the accelerometers’ scale factorsε: three-axis gyroscope constant drift vector∇: three-axis accelerometer constant drift vectorEg: installation error matrix from the b frame to the gyro frameEa: installation error matrix from the b frame to the accelerometer frameThe ideal error model of the gyros can be defined as follows [23]:(1)Ng=KgEgωb+ε.The matrix expansion of (1) is as follows:(2)NgxNgyNgz=Kgx000Kgy000Kgz1EgxyEgxEgyx1EgyzEgzxEgzy1ωxbωybωzb+εxεyεz.Similarly, the ideal error model of the accelerometers can be written as follows:(3)Na=KaEafb+∇.The matrix expansion of (3) is as follows:(4)NaxNayNaz=Kax000Kay000Kaz1EaxyEaxzEayx1EayzEaxxEazy1fxbfybfzb+∇x∇y∇z.


## 2.2. 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),(5)fb2=KaEa−1Na−∇TKaEa−1Na−∇=g2,ωb2=KgEg−1Ng−εTKgEg−1Ng−ε=Ω2,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:(6)f˜b2=K˜aE˜a−1Na−∇˜TK˜aE˜a−1Na−∇˜,ω˜b2=K˜gE˜g−1Ng−ε˜TK˜gE˜g−1Ng−ε˜.Symbols with upper wavy lines represent the actual values. Then the above objective optimization function can be established based on multiple static position observations(7)La=∑m=1Mg2−f˜bm22M−1,Lg=∑m=1MΩ2−ω˜bm22M−1,where M represents the total number of static positions in the calibration path.The objective optimization functionsLa and Lg 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. Thus a 24-position calibration scheme is designed to completely recognize those parameters.The whole calibration path is divided into three phases and every phase contains 8 locating places.Phase 1: assume the SINS’s body frame is theO−XbYbZb right-handed coordinates and axis Xb is fixed towards north (N). In position 1 (POS1), the axis Zb is adjusted towards up (U). Rotating 45° around axis Xb from POS1, we can get POS2. Then, do the same rotation around axis Xb 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(a).Phase 2: AxisYb is fixed towards north. Follow the same rotation operation as Phase 1 around axis Yb 8 times and get 8 positions. Also, keep standing time in each position. The rotation scheme of Phase 2 is shown in Figure 1(b).Phase 3: AxisZb is fixed towards north. Follow the same rotation operation as Phase 1 around axis Zb 8 times and get 8 positions. Also, keep standing time in each position. The rotation scheme of Phase 3 is shown in Figure 1(c).Figure 1
24-position calibration path. (a) Phase 1. (b) Phase 2. (c) Phase 3.
(a)(b)(c)


## 3. Parameter Optimization and Algorithm Improvement Process

### 3.1. Calibration Parameter Identification Process Based on DE Algorithm

The differential evolution algorithm shows great performance in solving optimization problems due to its simple structure, strong searching capability, and robustness [24, 25]. The 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 sizeNP, 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:(8)θi=Kax−iKay−iKaz−iEaxy−iEaxz−iEayx−iEayz−iEazx−iEazy−i∇x−i∇y−i∇z−i,i=1,2,…,NP.For every elementαjj=1,2,…12 in the vector θi, its value should be randomly generated between the lower bound LBjand upper bound UBj.(9)αj=LBj+rand0,1×UBj−LBj,j=1,2,…,12,where rand (0, 1) denotes a random number in the interval [0, 1].According to the results acquired from the traditional separate calibration, we may set boundaries among different elements. See details in Table1.Table 1
Accelerometer calibration parameter boundary settings.
AccelerometersScale factorsConstant biasesInstallation errorsSeparate calibration results13942.35422951−206.592510660.003364200.0019116813345.7248587787.245536070.00271530−0.0020401913062.44709817333.71033323−0.00739421−0.00044851Lower bound13000−400−0.01Upper bound140004000.01The 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:(10)f=1La.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:(11)θMt+1=θxt+Fθyt−θzt,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 ofF 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αji of θi with the corresponding components αjM of θM as follows:(12)αjct+1=αjMt+1,if rand0,1<CR,αjit,otherwise,where αjc 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θit and the trial vector θct+1, the objective function is introduced to identify next-generation members:(13)θit+1=θct+1,iffθct+1>fθit,θit,otherwise.The 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.


### 3.2. 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 Figure2.Figure 2
Accelerometer’s convergence process of fitness function (DE/rand/1).It is clear to see that slow convergence in the early stage and poor repeatability among fitness values make the identified calibration parameters unconvincing. Through 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 Table1. All these characteristics require the optimization algorithm must have strong traversal search ability and should not be too time-consuming for practical application.The 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. Therefore, 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:(14)θMt+1=θit+Fθbestt−θit+Fθxt−θyt,where θbestt 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.Figure 3
Accelerometer’s convergence process of fitness function (DE/current-to-best/1).As is shown in Figure3, “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. The integration mutation strategy we designed here is as follows:(15)θMt+1=1−bθxt+bθit+Fθbestt−θit+Fθyt−θzt,where b represents weight distribution factor and it is related to the current evolutional generation:(16)b=2tEG−tEG2.Assuming that the initial value ofb 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. The new results are illustrated in Figure 4.Figure 4
Accelerometer’s Convergence process of fitness function (integration strategy).It can be seen that under the guidance of the new integration strategy, the final fitness value is a lot better than those in Figures2 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. Therefore, 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:(17)Ft=Fmax−Fmax−Fmin11+α1/fmint−1/fmaxt+β1−t/EG,where Fmax,Fmin indicate the upper and lower limits of the mutation factor, fmaxt,fmint 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 self-adjustment. The combined mutation strategy can be further modified as follows:(18)θMt+1=1−bθxt+bθit+Ftθbestt−θit+Ftθyt−θzt.Conduct 3 more optimization process and make a comparison with the results shown in Figure4. The improvement effect after parameter dynamic adaptation is illustrated in Figure 5.Figure 5
Accelerometer’s convergence process of fitness function (integration strategy with parameter dynamic adaptation).The 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. The 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.Table 2
Statistical results of different DE algorithms in accelerometer calibration.
Traditional DEIDE (integration strategy only)IDE (integration strategy and parameter adaptation)Average value ofLa0.0455580.0081710.004854Standard deviation ofLa0.0158416.1852e − 048.5377e − 04


### 3.3. The Application and Correction of IDE Algorithm in Gyroscope Parameter Identification

The 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 Figure1 as the accelerometer. The new individual vector φi can be written as follows:(19)φi=Kgx−iKgy−iKgz−iEgxy−iEgxz−iEgyx−iEgyz−iEgzx−iEgzy−iεx−iεy−iεz−i,i=1,2,…,NP.The new boundary settings are shown in Table3.Table 3
Gyroscope calibration parameters boundary settings.
GyroscopesScale factorsConstant biasesInstallation errorsSeparate calibration results7921.74760617−208.33680363−0.004891320.000960907911.54386049−33.242995330.000581520.000185197924.57007346−54.564757740.004493500.00137873Lower bound7800−500−0.01Upper bound80001000.01The new fitness function is as follows:(20)f=1Lg.Then 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.Figure 6
Gyroscope’s convergence process of the fitness function.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. The key reason for this phenomenon is the ill-fitted data selection. The 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. Then the usage of static observation data causes considerable deviation in scale factors recognition and further makes the installation errors so disordered that some of the values have reached preset thresholds. Therefore, dynamic observation data is needed and a supplementary calibration path is added to the whole program. See the details in Figure7.Figure 7
Gyroscope supplementary calibration process.The main purpose of the extra calibration scheme is to add some certain motivation of external angular velocity to the gyroscopes. As is shown in Figure7, 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:(21)ωx1ωy1ωz1=cosωtsinωt0−sinωtcosωt00010ωiecosϕω+ωiesinϕ,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(22)Ngx1=Kgxωiecosϕsinωt+KgxEgxyωiecosϕcosωt+KgxEgxzω+ωiesinϕ+εx,Ngy1=KgyEgyxωiecosϕsinωt+Kgyωiecosϕcosωt+KgyEgyzω+ωiesinϕ+εy,Ngz1=KgzEgzxωiecosϕsinωt+KgzEgzyωiecosϕcosωt+Kgzω+ωiesinϕ+εz,After the turntable turningn circles, take the sum of both sides in (22) and average the results:(23)Ngx1¯=KgxEgxzω+ωiesinϕ+εx,Ngy1¯=KgyEgyzω+ωiesinϕ+εy,Ngz1¯=Kgzω+ωiesinϕ+εz.Similarly, during negative rotation (ROT2), we can get the following:(24)Ngx2¯=KgxEgxz−ω+ωiesinϕ+εx,Ngy2¯=KgyEgyz−ω+ωiesinϕ+εy,Ngz2¯=Kgz−ω+ωiesinϕ+εz.Adding (23) and (24), then we obtain the following:(25)Ngx1¯−Ngx2¯2ω=KgxEgxz,Ngy1¯−Ngy2¯2ω=KgyEgyz,Ngz1¯−Ngz2¯2ω=Kgz.Take the same steps forX 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:(26)f=1Lgm,Lgm=∑j=1JN¯j+−N¯j−2−K˜gE˜g2ωobj22J−1,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. ωob is the angular velocity of external excitation provided by the turntable. In our calibration scheme, the norm of ωob is set as 10°/s, 20°/s, and 30°/s, respectively. Use the improved fitness function to rerun the parameter optimization process. The convergence process is shown in Figure 8 and the final statistical results are listed in Table 4.Figure 8
Gyroscope’s convergence process of the new fitness function.Table 4
Statistical results of different DE algorithms in gyroscope calibration.
IDE with static calibration pathIDE with dynamic calibration pathAverage value ofLg/Lgm0.49350.003075Standard deviation ofLg/Lgm4.0192e − 067.8156e − 04It can be seen from Table4 that the high repeatability of the objective function occurs when the static calibration path is used. This 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. The 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.


## 3.1. Calibration Parameter Identification Process Based on DE Algorithm

The differential evolution algorithm shows great performance in solving optimization problems due to its simple structure, strong searching capability, and robustness [24, 25]. The 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 sizeNP, 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:(8)θi=Kax−iKay−iKaz−iEaxy−iEaxz−iEayx−iEayz−iEazx−iEazy−i∇x−i∇y−i∇z−i,i=1,2,…,NP.For every elementαjj=1,2,…12 in the vector θi, its value should be randomly generated between the lower bound LBjand upper bound UBj.(9)αj=LBj+rand0,1×UBj−LBj,j=1,2,…,12,where rand (0, 1) denotes a random number in the interval [0, 1].According to the results acquired from the traditional separate calibration, we may set boundaries among different elements. See details in Table1.Table 1
Accelerometer calibration parameter boundary settings.
AccelerometersScale factorsConstant biasesInstallation errorsSeparate calibration results13942.35422951−206.592510660.003364200.0019116813345.7248587787.245536070.00271530−0.0020401913062.44709817333.71033323−0.00739421−0.00044851Lower bound13000−400−0.01Upper bound140004000.01The 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:(10)f=1La.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:(11)θMt+1=θxt+Fθyt−θzt,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 ofF 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αji of θi with the corresponding components αjM of θM as follows:(12)αjct+1=αjMt+1,if rand0,1<CR,αjit,otherwise,where αjc 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θit and the trial vector θct+1, the objective function is introduced to identify next-generation members:(13)θit+1=θct+1,iffθct+1>fθit,θit,otherwise.The 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.


## 3.2. 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 Figure2.Figure 2
Accelerometer’s convergence process of fitness function (DE/rand/1).It is clear to see that slow convergence in the early stage and poor repeatability among fitness values make the identified calibration parameters unconvincing. Through 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 Table1. All these characteristics require the optimization algorithm must have strong traversal search ability and should not be too time-consuming for practical application.The 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. Therefore, 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:(14)θMt+1=θit+Fθbestt−θit+Fθxt−θyt,where θbestt 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.Figure 3
Accelerometer’s convergence process of fitness function (DE/current-to-best/1).As is shown in Figure3, “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. The integration mutation strategy we designed here is as follows:(15)θMt+1=1−bθxt+bθit+Fθbestt−θit+Fθyt−θzt,where b represents weight distribution factor and it is related to the current evolutional generation:(16)b=2tEG−tEG2.Assuming that the initial value ofb 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. The new results are illustrated in Figure 4.Figure 4
Accelerometer’s Convergence process of fitness function (integration strategy).It can be seen that under the guidance of the new integration strategy, the final fitness value is a lot better than those in Figures2 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. Therefore, 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:(17)Ft=Fmax−Fmax−Fmin11+α1/fmint−1/fmaxt+β1−t/EG,where Fmax,Fmin indicate the upper and lower limits of the mutation factor, fmaxt,fmint 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 self-adjustment. The combined mutation strategy can be further modified as follows:(18)θMt+1=1−bθxt+bθit+Ftθbestt−θit+Ftθyt−θzt.Conduct 3 more optimization process and make a comparison with the results shown in Figure4. The improvement effect after parameter dynamic adaptation is illustrated in Figure 5.Figure 5
Accelerometer’s convergence process of fitness function (integration strategy with parameter dynamic adaptation).The 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. The 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.Table 2
Statistical results of different DE algorithms in accelerometer calibration.
Traditional DEIDE (integration strategy only)IDE (integration strategy and parameter adaptation)Average value ofLa0.0455580.0081710.004854Standard deviation ofLa0.0158416.1852e − 048.5377e − 04


## 3.3. The Application and Correction of IDE Algorithm in Gyroscope Parameter Identification

The 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 Figure1 as the accelerometer. The new individual vector φi can be written as follows:(19)φi=Kgx−iKgy−iKgz−iEgxy−iEgxz−iEgyx−iEgyz−iEgzx−iEgzy−iεx−iεy−iεz−i,i=1,2,…,NP.The new boundary settings are shown in Table3.Table 3
Gyroscope calibration parameters boundary settings.
GyroscopesScale factorsConstant biasesInstallation errorsSeparate calibration results7921.74760617−208.33680363−0.004891320.000960907911.54386049−33.242995330.000581520.000185197924.57007346−54.564757740.004493500.00137873Lower bound7800−500−0.01Upper bound80001000.01The new fitness function is as follows:(20)f=1Lg.Then 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.Figure 6
Gyroscope’s convergence process of the fitness function.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. The key reason for this phenomenon is the ill-fitted data selection. The 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. Then the usage of static observation data causes considerable deviation in scale factors recognition and further makes the installation errors so disordered that some of the values have reached preset thresholds. Therefore, dynamic observation data is needed and a supplementary calibration path is added to the whole program. See the details in Figure7.Figure 7
Gyroscope supplementary calibration process.The main purpose of the extra calibration scheme is to add some certain motivation of external angular velocity to the gyroscopes. As is shown in Figure7, 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:(21)ωx1ωy1ωz1=cosωtsinωt0−sinωtcosωt00010ωiecosϕω+ωiesinϕ,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(22)Ngx1=Kgxωiecosϕsinωt+KgxEgxyωiecosϕcosωt+KgxEgxzω+ωiesinϕ+εx,Ngy1=KgyEgyxωiecosϕsinωt+Kgyωiecosϕcosωt+KgyEgyzω+ωiesinϕ+εy,Ngz1=KgzEgzxωiecosϕsinωt+KgzEgzyωiecosϕcosωt+Kgzω+ωiesinϕ+εz,After the turntable turningn circles, take the sum of both sides in (22) and average the results:(23)Ngx1¯=KgxEgxzω+ωiesinϕ+εx,Ngy1¯=KgyEgyzω+ωiesinϕ+εy,Ngz1¯=Kgzω+ωiesinϕ+εz.Similarly, during negative rotation (ROT2), we can get the following:(24)Ngx2¯=KgxEgxz−ω+ωiesinϕ+εx,Ngy2¯=KgyEgyz−ω+ωiesinϕ+εy,Ngz2¯=Kgz−ω+ωiesinϕ+εz.Adding (23) and (24), then we obtain the following:(25)Ngx1¯−Ngx2¯2ω=KgxEgxz,Ngy1¯−Ngy2¯2ω=KgyEgyz,Ngz1¯−Ngz2¯2ω=Kgz.Take the same steps forX 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:(26)f=1Lgm,Lgm=∑j=1JN¯j+−N¯j−2−K˜gE˜g2ωobj22J−1,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. ωob is the angular velocity of external excitation provided by the turntable. In our calibration scheme, the norm of ωob is set as 10°/s, 20°/s, and 30°/s, respectively. Use the improved fitness function to rerun the parameter optimization process. The convergence process is shown in Figure 8 and the final statistical results are listed in Table 4.Figure 8
Gyroscope’s convergence process of the new fitness function.Table 4
Statistical results of different DE algorithms in gyroscope calibration.
IDE with static calibration pathIDE with dynamic calibration pathAverage value ofLg/Lgm0.49350.003075Standard deviation ofLg/Lgm4.0192e − 067.8156e − 04It can be seen from Table4 that the high repeatability of the objective function occurs when the static calibration path is used. This 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. The 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.


## 4. Experimental Results and Algorithm Performance Comparison

### 4.1. Experimental Condition Setting

The system’s configuration is shown in Figure9(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. The 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.Figure 9
Configuration of the test system. (a) SINS. (b) Turntable. (c) Power supply and acquisition equipment.
(a)(b)(c)


### 4.2. 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. The 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 Subsections2.2 and 3.3. Finally, statistical results of parameter identification are shown in Table 5 and navigation results are illustrated in Figure 10.Table 5
Statistical results of different methods in SINS calibration.
Comparison in accelerometer parameter identificationGAPSOIDESuccess rate (%)73.3386.6796.67Average value ofLa0.15130.033870.004854Standard deviation ofLa0.10960.047838.5377e − 04Comparison in gyroscope parameter identificationSuccess rate (%)6083.33100Average value ofLgm0.083600.050490.003075Standard deviation ofLgm0.094560.038067.8156e-04PS: GA represents genetic algorithm and PSO represents particle swarm optimization.Figure 10
Comparison of navigation results. (a) East velocity errors. (b) North velocity errors. (c) Latitude errors. (d) Longitude errors.
(a)(b)(c)(d)It is clear that the IDE method has the best performance in Table5. 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. Therefore, parameters gained from IDE optimization calibration are proved to be accurate and show better navigation performance.


### 4.3. 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 isM, and the population size is generally 3–8 times the dimension, taking 5 M. The 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:(27)ONMO5NM2ONM.The overall time complexity of the traditional algorithm is as follows:(28)ONM+5NM2+NM=ONM2.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:(29)ONM2O5NM2ONM.And time complexity of the new algorithm is as follows:(30)ONM2+5NM2+NM=ONM2.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. Therefore, the new method adds a certain storage unit.It can be seen that the proposed IDE algorithm represents the compromise between accuracy and efficiency. The 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.


## 4.1. Experimental Condition Setting

The system’s configuration is shown in Figure9(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. The 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.Figure 9
Configuration of the test system. (a) SINS. (b) Turntable. (c) Power supply and acquisition equipment.
(a)(b)(c)


## 4.2. 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. The 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 Subsections2.2 and 3.3. Finally, statistical results of parameter identification are shown in Table 5 and navigation results are illustrated in Figure 10.Table 5
Statistical results of different methods in SINS calibration.
Comparison in accelerometer parameter identificationGAPSOIDESuccess rate (%)73.3386.6796.67Average value ofLa0.15130.033870.004854Standard deviation ofLa0.10960.047838.5377e − 04Comparison in gyroscope parameter identificationSuccess rate (%)6083.33100Average value ofLgm0.083600.050490.003075Standard deviation ofLgm0.094560.038067.8156e-04PS: GA represents genetic algorithm and PSO represents particle swarm optimization.Figure 10
Comparison of navigation results. (a) East velocity errors. (b) North velocity errors. (c) Latitude errors. (d) Longitude errors.
(a)(b)(c)(d)It is clear that the IDE method has the best performance in Table5. 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. Therefore, parameters gained from IDE optimization calibration are proved to be accurate and show better navigation performance.


## 4.3. 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 isM, and the population size is generally 3–8 times the dimension, taking 5 M. The 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:(27)ONMO5NM2ONM.The overall time complexity of the traditional algorithm is as follows:(28)ONM+5NM2+NM=ONM2.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:(29)ONM2O5NM2ONM.And time complexity of the new algorithm is as follows:(30)ONM2+5NM2+NM=ONM2.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. Therefore, the new method adds a certain storage unit.It can be seen that the proposed IDE algorithm represents the compromise between accuracy and efficiency. The 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.


## 5. 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. The 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.



---

*Source: 6684119-2021-01-29.xml*