# System-Level Temperature Compensation Method for the RLG-IMU Based on HHO-RVR

Authors: Hao Liang; Yumin Tao; Meijiao Wang; Yu Guo; Xingfa Zhao
Journal: Journal of Sensors
 (2021)

---

## Abstract

The ring laser gyro inertial measurement unit has many systematic error terms and influences each other. These error terms show a complex nonlinear drift that cannot be ignored when the temperature changes, which seriously affects the stability time and output accuracy of the system. In this paper, a system-level temperature modeling and compensation method is proposed based on the relevance vector regression method. First, all temperature-related parameters are modeled; meanwhile, the Harris hawks optimization algorithm is used to optimize each model parameter. Then, the system compensation is modeled to stabilize the system output to the desired temperature. Compared with the least square method, the fitting performance comparison and the system dynamic compensation experiment prove this method’s superiority. The root mean square error, the mean absolute error, theR-squared, and the variance of residual increased by an average of 35.27%, 39.29%, 2.29%, and 30.34%, respectively.


---

## Body

## 1. Introduction

In the inertial sensors, the ring laser gyroscope (RLG) [1–3] has many advantages such as small random drift, wide dynamic range, fast startup speed, and high reliability. It has a high value in the field of inertial navigation. As a representative of strap-down inertial navigation, the ring laser gyroscope inertial measurement unit (RLG-IMU) [4–6], which uses RLG as its core measurement component, is widely used in the field of inertial navigation for its outstanding advantages of small size, high precision, strong environmental adaptability, and low cost. How to improve the performance of the RLG and the RLG-IMU has always been the focus of researchers. After years of research, researchers have done much work on various physical fields that affect the system performance, such as temperature, magnetic field, vibration environment, and how to improve the performance of the RLG-IMU under the influence of various physical fields.Temperature is an important factor affecting the output precision and startup stability time of the RLG-IMU. On the one hand, the operating temperature range of the RLG-IMU is extensive, which requires strong adaptability of the system. On the other hand, many RLG-IMU devices have a temperature control system composed of the heating part, temperature feedback part, and temperature control part. At the same time, the RLG-IMU has many circuit components, which also generate heat. Other components inside the IMU, such as magnetic coils, also generate heat. After the RLG-IMU is started, a large amount of heat will be generated inside the system. The internal temperature changes, and IMU’s output is unstable due to the internal heating part. Such changes in the system’s internal and external temperature affect the RLGs and quartz flexible accelerometers’ performance and the RLG-IMU internal structure’s characteristics. The change of structural characteristics leads to the thermal deformation of inertial devices, which causes a lot of fluctuation and error in output data of the RLG-IMU and limits the startup time and performance of the RLG-IMU. Therefore, it is necessary to study the temperature effect of the IMU system to improve its temperature adaptability and reduce the system stability time.At present, there are two methods to solve the influence of temperature on the IMU system: temperature control and temperature compensation. The method of temperature control makes the system quickly heat up to reach the temperature control point, to make the system stable as soon as possible. After the RLG-IMU is started, the internal temperature control system starts to work. The system temperature generally reaches a relatively stable state after a relatively long time (generally dozens of minutes). At this time, the output of the RLG-IMU is stable, and the output accuracy can meet the requirements. Before that, the output of the IMU was unstable due to the drastic temperature change. However, this method’s problem is that temperature control requires a “long setting time” or called a “long startup time.” The accuracy of the temperature control system will also affect the IMU’s performance, and additional hardware needs to be added, and the system’s power consumption increased. The method of temperature compensation only needs to obtain the error model under the temperature, and the compensation can be started at the startup of the IMU. Theoretically, the system’s output can be stable as soon as it is powered on, but the accuracy of the model establishment will affect the accuracy of the compensation. How to reduce the startup stable time by using the system temperature compensation method is the purpose of this paper. By modeling the system measurement equation’s parameters under different temperatures, the output of the RLG-IMU at different temperatures can be converted to the output at a stable temperature. In this way, the operation of the RLG-IMU will not be affected by the change of temperature so that IMU can meet the requirements without “long setting time.”Several methods have been used for the RLG’s or quartz flexible accelerometer’s zero position temperature compensation, such as stepwise regression [7, 8], artificial neural network [9–12], support vector machine [13–15], and K-mean [16]. The IMU system’s various error terms are very complicated, including the zero position and the influence of the scale factor. The system structure changes will also affect the zero position and scale factor of the sensors. These factors not only are affected by temperature but also have a coupling relationship with each other. Compensating for the RLG or quartz flexible accelerometer zero position alone, without considering the effects of other error terms, cannot cover the system’s error, which will affect the system’s compensation accuracy. Therefore, all drifts must be compensated at the system level. At present, the least square method (LSM) is usually used for system-level compensation in engineering [17–20]. The IMU has more than a dozen parameter items related to temperature. When these parameters change with temperature in a sophisticated nonlinear manner, this traditional modeling method’s fitting accuracy is limited.Machine learning has powerful capabilities in the regression and prediction of complex functions. Based on the Bayesian framework, Tipping proposed the relevance vector machine [21–23]. Compared with the support vector machine, the relevance vector machine uses fewer vectors and has stronger sparsity. Although the training time is longer than that of the support vector machine, the prediction time is much less than that of the support vector machine. The kernel function does not need to meet Mercer’s condition: in a finite input space, the function K is a map. If the kernel matrix is positive semidefinite, then the function K can be a kernel function. In SVM theory, the kernel function must satisfy Mercer’s condition. In relevance vector machine theory, because of the difference between the relevance vector machine and the SVM architecture, the kernel function does not have to satisfy Mercer’s condition, so more kernel functions can be selected. In the case of fewer training data samples, it can ensure excellent generalization ability. The method of using the relevance vector machine for regression is called relevance vector regression (RVR) [24, 25].In order to solve the problem of high-precision compensation of the RLG-IMU system-level temperature, this paper proposes a system-level temperature error model and compensation method for the RLG-IMU based on the RVR. According to the IMU system’s input-output model, all the parameters that affect the output are modeled at the system level and compensated so that the system-level temperature error compensation is more comprehensively achieved. Since the setting of the relevance vector machine kernel function’s width parameter has a severe impact on the regression accuracy, it is necessary to optimize this parameter. The Harris hawks optimization (HHO) [26] is a novel metaheuristic [27, 28] optimization algorithm. Compared with the genetic algorithm [29] and particle swarm optimization algorithm [30], it has fast optimization speed and high precision. In this paper, the HHO algorithm is used to optimize the kernel width parameter in the relevance vector machine to improve the regression accuracy of the model, so the method is called HHO-RVR.For the RLG-IMU system, the influence of temperature on the RLG-IMU system is multifaceted. Although the influence of the temperature on the RLG-IMU is multifaceted, all the effects will be reflected in the measured output pulse at different temperatures. According to the pulse-angular velocity equation and the pulse-apparent acceleration equation, if the objective angular velocity and apparent acceleration are considered real and do not change with temperature, the pulse change caused by temperature is caused by other equation parameters with the change of temperature. This paper is aimed at obtaining the fitting model of the parameters in the equation varying with temperature. When the system works at all temperatures, the system’s output pulse can be converted into a pulse output at a selected temperature. Thus, the output of the system is more stable when the temperature changes. This compensation method does not need to consider every factor affected by temperature change in the system. Therefore, we call this compensation method to be system-level temperature compensation.This paper’s structure is as follows: Section1 introduces the background, current problems, and this paper’s work. Section 2 analyzes the influence of temperature on the IMU system. Section 3 introduces the relevance vector machine regression theory and the HHO optimization process and establishes the system parameters’ temperature and compensation models. Section 4 introduces the experimental methods, experimental results, and analyses to verify the method’s effectiveness and superiority. Section 5 summarizes the whole paper.


## 2. Analysis of Temperature Effects

As shown in Figure1, a typical RLG-IMU system consists of three orthogonal RLGs and three orthogonal quartz flexible accelerometers mounted on the base. The IMU measures angular velocity in three directions with three RLGs (Gx, Gy, and Gz) and apparent acceleration in three directions with three quartz flexible accelerometers (Ax, Ay, and Az).Figure 1
RLG-IMU diagram.The influence of temperature on IMU mainly includes three aspects:(1)
Impact on the RLGThe laser gyroscope is based on the Sagnac effect principle and uses the optical path difference to measure the rotational angular velocity. The ring laser gyro is essentially an active ring laser. It is a laser source filled with a helium-neon mixture. The effect of temperature on the gyro is comprehensive [31], such as the change of material characteristics of the laser gyroscope, the influence on the length of the resonant cavity and its coplanarity, the change of the gas flow rate, the change of the characteristics of the mirror, the deformation of the capillary, changes of lock-in characteristics, Langmuir flow, and discharge symmetry. These changes are reflected in changes in the scale factor and zero bias [32].
(2)
Impact on the quartz flexible accelerometerWhen the quartz flexible accelerometer is started, the fluctuation of the internal temperature field will cause the differential capacitance sensor to generate a larger current, further increasing the moment coil’s temperature. The increase in temperature affects the magnetic materials and torque coils in the accelerometer, causing changes in the magnetic flux and causing the accelerometer scale factor’s temperature drift. Simultaneously, due to the thermal imbalance in the quartz flexible accelerometer, the arm will produce slight distortion, which will affect the stability of the accelerometer zero bias. Besides, changes in the arm length can affect the scale factor of the accelerometer.(3)
Impact on IMU system structureFor the IMU system, the three RLGs and three quartz flexible accelerometers are fixed on the base and cannot be strictly orthogonal. They have relative installation errors. Due to the generally strong rigidity of the base, the installation errors are relatively stable. However, the geometry of the base is also a small change that will occur with the temperature change, and this change will also be reflected in the parameter drift of the single sensor.From the above analysis, it can be concluded that the influence of temperature on the RLG-IMU system is multifaceted and complex. Temperature affects not only the individual sensor but also the system structure. Structural changes are also reflected in the output of the sensor. It is not straightforward to analyze and compensate for these factors separately. Therefore, this paper uses the system-level compensation method to integrate all the influencing factors to model and compensate directly.


## 3. Relevance Vector Regression and Temperature Compensation Modeling

### 3.1. Function Regression Theory of the Relevance Vector Machine

The relevance vector machine is based on the Bayesian theory. In linear regression problems, the input and output can be described asy=wTx+c, where w is the weight vector, x is the input vector, and c is the offset. When the relationship is nonlinear, it can be expressed as y=wTϕx or yx;w. Given a set of input target data sets: S=xn,tnn=1N,xn∈Rn,tn∈R, N is the number of data samples. Consider that the target set is data samples superimposed with noise, expressed as tn=yxn;w+ϵn, where ϵn is the noise of the Gaussian distribution satisfying the mean of 0 and the variance of σ2, that is, εn~N0,σ2; the conditional probability of the target is [33–35]
(1)ptn∣x=Ntn∣yxn,σ2.tn is independent of each other; then, the likelihood functions of all corresponding data sets are expressed as
(2)pt∣w,σ2=∏i=1NNti∣yxi;w,σ2=2πσ2−N/2exp−12σ2t−Φw2,where t=t1⋯tNT, w=w0⋯wNT, and Φ is the N×N+1 matrix with Φ=ϕx1,ϕx2,⋯,ϕxNT, wherein ϕxn=1,Kxn,x1,Kxn,x2,⋯,Kxn,xNT.Kxi,xj is the kernel function.In order to establish the relationship model between input and output, the weightw in the model is given prior probability.
(3)pw∣α=∏i=0Nwi∣0,αi−1.Hyperparameterαi is used to describe the inverse variance of each wi.The posterior probability of the weight is(4)pw∣t,α,σ2=pt∣w,σ2pw∣αpt∣α,σ2=2π−N+1/2Σ−1/2exp−12w−μTΣ−1w−μ,where the covariance and mean are Σ=σ−2ΦTΦ+A−1, μ=σ−2ΣΦTt, and A=diagα0,α1,⋯,αN.Use the maximum likelihood method to obtain the optimal solutionαMP,σMP2 of pt∣α,σ2.For new inputsx∗, the predicted distribution can be calculated.
(5)pt∗∣t,αMP,σMP2=∫pt∗∣w,σMP2pw∣t,αMP,σMP2dw.Since both of the integrands are Gaussian, this can be calculated.(6)pt∗∣t,αMP,σMP2=Nt∗∣y∗,σ∗2y∗=μTϕx∗,σ∗2=σMP2+ϕx∗TΣϕx∗.In equation (2), the kernel function is used to map the feature vector to the high-dimensional space, which can reduce the calculation complexity. Commonly used are the linear kernel function, polynomial kernel function, sigmoid kernel function, and radial basis kernel function. The linear kernel function is mainly used in the case of linear separability. The application effect of the polynomial kernel function on nonlinear data is not ideal, and there are many parameters to be adjusted, which increases the complexity of the model. The sigmoid kernel function has good performance in dealing with nonlinear data, and two parameters need to be adjusted. The radial basis function is a kind of kernel function with robust localization, which can realize nonlinear mapping. The radial basis function is better than the linear kernel function in nonlinear data processing. The shape of the function is a bell-shaped curve, and only one parameter controls the width of the curve. Compared with the polynomial kernel function and the sigmoid kernel function, the radial basis function has fewer parameters, which greatly reduces the model’s complexity. The radial basis function is the most widely used kernel function. It has good performance in dealing with large and small samples. This paper uses the radial basis kernel function, whose expression is Kxi,xj=exp−xi−xj2/2γ2, where γ denotes the width parameter of the kernel function. Different γ affects the performance of the relevance vector machine, so this parameter must be optimized, and the optimization algorithm based on metaheuristic can significantly improve the optimization efficiency.


### 3.2. Parameter Optimization Based on HHO

HHO is a novel metaheuristic swarm intelligence algorithm. Compared with the genetic algorithm and the particle swarm optimization algorithm, it has better exploration capabilities and convergence speed [26]. When the hawks catch their prey, they adopt different strategies in different stages. The HHO algorithm can find the optimal value by modeling hawks’ behavior in each stage of the predation process.The HHO algorithm includes the exploration phase and the exploitation phase. In the process of predation, with the increase in the moving distance of prey, the physical strength of prey will gradually decrease, which is called escape energy. The algorithm can transfer from the exploration phase to the exploitation phase through the value range of prey’s escape energy. The escape energy is modeled as [26]
(7)E=2E01−tT,where E denotes the escaping energy of the prey, T denotes the maximum number of iterations, t denotes the current number of iterations, and E0 denotes the initial energy. E0 randomly changes inside the interval −1,1 at each iteration.When∣E∣≥1, enter the exploration phase. The hawks inhabit different positions randomly and use two strategies to detect prey. An equal probability q is used to distinguish the two strategies. They perch based on the positions of other hawks and the prey when q<0.5 and perch on random tall trees when q≥0.5. The exploration phase can be expressed as
(8)Xt+1=Xrandt−r1Xrandt−2r2Xt,Xpreyt−Xmt−r3LB+r4UB−LB,q⩾0.5,q<0.5,where Xt+1 denotes the position vector of the next-generation hawks, Xt denotes the position vector of the current-generation hawks, Xpreyt denotes the position of the prey, r1,r2, r3, r4, and q are random values between 0,1 and are updated every generation, LB and UB denote the boundary of the variable, Xrandt denotes a random position selected from the current-generation hawks, and Xmt denotes the average position of the hawks.When∣E∣<1, enter the exploitation phase. At this time, use four siege strategies. r is a random number between 0,1 and indicates the chance for the prey to escape. r<0.5 indicates that the prey can escape, and r≥0.5 indicates that the prey cannot escape successfully [26]. The use of soft and hard besiege strategies depends on the prey’s energy E. When ∣E∣≥0.5, it means that the prey still has high energy. At this time, the Hawks adopt the soft besiege strategy, constantly approaching and chasing the prey to consume the energy of the prey. When ∣E∣<0.5, it indicates that the prey’s energy is weak. At this time, the Hawks can easily catch the prey by using the hard besiege strategy.(1) Soft Besiege. When r≥0.5 and ∣E∣≥0.5, it means that the prey has enough escape energy and jumps randomly. At this time, the hawks encircle it softly to exhaust the prey, expressed as
(9)Xt+1=Xpreyt−Xt−EJXpreyt−Xt,where J is a random number between 0,2, denoting the random jumping energy of prey, which is used to simulate the prey’s activity.(2) Hard Besiege. When r≥0.5 and ∣E∣<0.5, it means that the prey is weak and has low escape energy. At this time, the following equation is used to update the position.
(10)Xt+1=Xpreyt−EXpreyt−Xt.(3) Soft Besiege with Progressive Rapid Dives. When r<0.5 and ∣E∣≥0.5, the prey has enough energy to escape. At this time, equations (11) and (12) are used to update the position, and the Lévy flight is used to simulate the sudden movement of the hawks and the prey. The Lévy flight is expressed as equation (13), which can further enhance the ability of the algorithm to jump out of the local optimum.
(11)Y=Xpreyt−EJXpreyt−Xt,(12)Z=Y+S×LFD,(13)LFx=0.01×u×σv1/β,σ=Γ1+β×sinπβ/2Γ1+β/2×β×2β−1/21/β,where u and v are random values inside 0,1, β is a default constant set to 1.5, and Γz is the Gramma function: Γz=∫0∞tz−1e−tdt.Compare the two update methods, and select the most suitable position update method according to the fitness function value.(14)Xt+1=Y,ifFY<FXt,Z,ifFZ<FXt.(4) Hard Besiege with Progressive Rapid Dives. When r<0.5 and ∣E∣<0.5, it means that the prey does not have enough energy. At this time, equations (15) and (16) are used to update the position and equation (14) is used to select the appropriate update method.
(15)Y=Xpreyt−EJXpreyt−Xmt,(16)Z=Y+S×LFD.The relevance vector machine is used for data fitting and regression. HHO is used to optimize the kernel width parameter of the relevance vector machine.The flow of each HHO-RVR is as follows:Step 1.
Initialize the number of HHO populations and the number of iterations. The hawk’s position in the HHO algorithm represents the kernel width parameter of the relevance vector machine. The position of each hawk in the first-generation group is randomly distributed.Step 2.
Determine the fitness function. The sample data is randomly divided into five parts, 4 of which are used for the relevance vector machine training, and the remaining 1 part is used for verification, and the test accuracy of the relevance vector machine is used as the fitness function of HHO.Step 3.
Calculate the fitness value of each hawk, and use the hawk corresponding to the optimal fitness as the prey.Step 4.
Update the positions of the hawks, according to equations (7)–(16).Step 5.
Perform iterative operations. Repeat steps 3 and 4 until the number of iterations is satisfied and the calculation is stopped. At this time, the prey position is the optimal kernel width parameter.


### 3.3. RLG-IMU Temperature Modeling

The physical input of the RLG is the angular velocity, and the measured output is the number of pulses. The computer counts the number of pulses in the sampling time. According to the measurement equation of the system, the angular velocity of the RLG is obtained. The pulse-angular velocity measurement model of the RLGs can be expressed as(17)Ngτ=KgEgωgm+Dg.Ng=NgxNgyNgzT is the vector of the number of pulses output by the RLGs on three orthogonal axes within the sampling time τ. Kg=diagKgx,Kgy,Kgz is the scale factor matrix.Dg=DgxDgyDgzT is the bias matrix. ωgm=ωgxmωgymωgzmT is the angular velocity component vector of the three axes in the IMU measurement coordinate system.
(18)Eg=EgxxEgyxEgzxEgxyEgyyEgzyEgxzEgyzEgzzis the installation error matrix of the RLG relative inertial coordinate system for each axis.In equation (17), the output vector Ng of the RLGs changes with temperature and Eg is very little affected by temperature. Because the influence of temperature on the RLG is multifaceted, the change of the sensor’s output to the measured angular velocity caused by the temperature change is directly reflected in the scale factor matrix Kg and bias matrix Dg, and the structural changes will also be reflected on Kg and Dg, so equation (17) can be expressed as
(19)NgTτ=KgTEgωgm+DgT.The physical input of the quartz flexible accelerometer is the apparent acceleration, and the measured output is the number of pulses. The computer counts the number of pulses in the sampling time. According to the measurement equation of the system, the apparent acceleration of the quartz flexible accelerometer is obtained. The pulse-apparent acceleration measurement model of the quartz flexible accelerometers is expressed as(20)Naτ=KaEaAam+Da.Na=NaxNayNazT is the vector of the number of pulses output by the quartz flexible accelerometers on three orthogonal axes within the sampling time τ. Ka=diagKax,Kay,Kaz is the scale factor matrix of the quartz flexible accelerometers. Since the quartz flexible accelerometer is very sensitive to positive and negative scale factors, it is divided into positive and negative: the positive scale factor matrix is expressed as Kap=diagKapx,Kapy,Kapz and the negative scale factor matrix is expressed as Kan=diagKanx,Kany,Kanz. Da=DaxDayDazT is the bias matrix of the quartz flexible accelerometers. Aam=AaxmAaymAazmT is the apparent acceleration component vector of the three axes in the IMU measurement coordinate system.
(21)Ea=EaxxEayxEazxEaxyEayyEazyEaxzEayzEazzis the installation error matrix of the quartz flexible accelerometer relative inertial coordinate system for each axis. Similarly, equation (20) is expressed as
(22)NaTτ=KaTEaAam+DaT.In equations (19) and (22), the input parameter affected by temperature is called the temperature-related parameter, and there are 15 temperature-related parameters in total.Define a setP=Pi∣i=1⋯15 to denote 15 temperature-related parameters. Make Pi∣i=1⋯6 denote each element in Dg,Kg and Pi∣i=7⋯15 denote each element in Da,Kap,Kan.Since the temperature measurement position of each RLG and quartz flexible accelerometer is different, the temperature measurement curve of each temperature-related parameter will be different. Define a temperature setT=Ti∣i=1⋯15. Ti denotes the temperature change vector of each temperature-related parameter and corresponds one-to-one with Pi; then, the corresponding relationship between set P and temperature can be expressed as PT=PiTi∣i=1⋯15. From this, 15 models can be built as follows:
(23)Pi=Miwi,Ti,i=1,2,⋯,15,where Miwi,Ti is the model corresponding to each temperature-related parameter of the IMU and wi is the optimal width parameter of the kernel function of each model.


### 3.4. System Compensation Model

Ideally, in equations (19) and (22), angular velocity ωgm and apparent acceleration Aam are independent of temperature. First, the output models of angular velocity and apparent acceleration can be obtained according to IMU measurement equations (19) and (22) and the models of temperature-related parameters. The values of the temperature-related parameters at the desired temperature TE are then selected as constants in the angular/apparent acceleration measurement equation to obtain the compensation pulse. The desired temperature TE is a hypothetical stable working state. The purpose is that when the system works at any temperature, the measured value of physical quantity can be converted to the value at the desired temperature so that the system will not be affected by the temperature change. The desired temperature is not unique. It can be set according to the actual needs and requirements. No matter how much the desired temperature is set, the measured value of different temperatures can be converted to the measured value at the desired temperature by compensating for the measurement equation’s parameters to realize the temperature compensation.For the RLGs, according to equation (19), the output model of ωgm is
(24)ωgm=Eg−1NgTgτKg−1Tg−DgTg.Tg corresponding to each RLG is Tgx, Tgy, and Tgz.The compensation equation of RLG output pulse at the desired temperatureTE is
(25)NgCτ=KgTEEgωgm+DgTE.NgC=NgxCNgyCNgzCT is the output pulse of the three RLGs after compensation. KgTE and DgTE are, respectively, the scale factor matrix and bias matrix of the three RLGs at the desired temperature TE.By bringing equation (24) into (25), there is
(26)NgCτ=KgTEEgEg−1NgTgτKg−1Tg−DgTg+DgTE.Corresponding to different RLGs,Tg includes Tgx, Tgy, and Tgz. The data set Tg,Kg can be obtained by temperature experiment. The data set was used to train the relevance vector machine. At the same time, the optimization process introduced in Section 3.2 is used to optimize the kernel width parameter of the relevance vector machine, so as to obtain the optimized KgTg model. Similarly, the model of DgTg can be obtained. Since the temperature models of KgTg and DgTg have been obtained, it is possible to compensate for the output NgTg at all temperatures, thereby compensating for Ng to NgTE at different temperatures.For the quartz flexible accelerometers, according to equation (22), the output model of Aam is
(27)Aam=Ea−1NaTaτKa−1Ta−DaTa.Ta corresponding to each quartz flexible accelerometer is Tax, Tay, and Taz.The compensation equation of quartz flexible accelerometer output pulse at the desired temperatureTE is
(28)NaCτ=KaTEEaAam+DaTE.NaC=NaxCNayCNazCT is the output pulse of the three quartz flexible accelerometers after compensation. KaTE and DaTE are, respectively, the scale factor matrix and bias matrix of the three quartz flexible accelerometers at the desired temperature TE.By bringing equation (27) into (28), there is
(29)NaCτ=KaTEEaEa−1NaTaτKa−1Ta−DaTa+DaTE.Corresponding to different quartz flexible accelerometers,Ta includes Tax, Tay, and Taz. Similarly, the models of KaTa and DaTa can be obtained. Since the temperature models of KaTa and DaTa have been obtained, it is possible to compensate for the output NaTa at all temperatures, thereby compensating for NaTa to NaTE at different temperatures.Equations (26) and (29) are the compensation model of the IMU.


## 3.1. Function Regression Theory of the Relevance Vector Machine

The relevance vector machine is based on the Bayesian theory. In linear regression problems, the input and output can be described asy=wTx+c, where w is the weight vector, x is the input vector, and c is the offset. When the relationship is nonlinear, it can be expressed as y=wTϕx or yx;w. Given a set of input target data sets: S=xn,tnn=1N,xn∈Rn,tn∈R, N is the number of data samples. Consider that the target set is data samples superimposed with noise, expressed as tn=yxn;w+ϵn, where ϵn is the noise of the Gaussian distribution satisfying the mean of 0 and the variance of σ2, that is, εn~N0,σ2; the conditional probability of the target is [33–35]
(1)ptn∣x=Ntn∣yxn,σ2.tn is independent of each other; then, the likelihood functions of all corresponding data sets are expressed as
(2)pt∣w,σ2=∏i=1NNti∣yxi;w,σ2=2πσ2−N/2exp−12σ2t−Φw2,where t=t1⋯tNT, w=w0⋯wNT, and Φ is the N×N+1 matrix with Φ=ϕx1,ϕx2,⋯,ϕxNT, wherein ϕxn=1,Kxn,x1,Kxn,x2,⋯,Kxn,xNT.Kxi,xj is the kernel function.In order to establish the relationship model between input and output, the weightw in the model is given prior probability.
(3)pw∣α=∏i=0Nwi∣0,αi−1.Hyperparameterαi is used to describe the inverse variance of each wi.The posterior probability of the weight is(4)pw∣t,α,σ2=pt∣w,σ2pw∣αpt∣α,σ2=2π−N+1/2Σ−1/2exp−12w−μTΣ−1w−μ,where the covariance and mean are Σ=σ−2ΦTΦ+A−1, μ=σ−2ΣΦTt, and A=diagα0,α1,⋯,αN.Use the maximum likelihood method to obtain the optimal solutionαMP,σMP2 of pt∣α,σ2.For new inputsx∗, the predicted distribution can be calculated.
(5)pt∗∣t,αMP,σMP2=∫pt∗∣w,σMP2pw∣t,αMP,σMP2dw.Since both of the integrands are Gaussian, this can be calculated.(6)pt∗∣t,αMP,σMP2=Nt∗∣y∗,σ∗2y∗=μTϕx∗,σ∗2=σMP2+ϕx∗TΣϕx∗.In equation (2), the kernel function is used to map the feature vector to the high-dimensional space, which can reduce the calculation complexity. Commonly used are the linear kernel function, polynomial kernel function, sigmoid kernel function, and radial basis kernel function. The linear kernel function is mainly used in the case of linear separability. The application effect of the polynomial kernel function on nonlinear data is not ideal, and there are many parameters to be adjusted, which increases the complexity of the model. The sigmoid kernel function has good performance in dealing with nonlinear data, and two parameters need to be adjusted. The radial basis function is a kind of kernel function with robust localization, which can realize nonlinear mapping. The radial basis function is better than the linear kernel function in nonlinear data processing. The shape of the function is a bell-shaped curve, and only one parameter controls the width of the curve. Compared with the polynomial kernel function and the sigmoid kernel function, the radial basis function has fewer parameters, which greatly reduces the model’s complexity. The radial basis function is the most widely used kernel function. It has good performance in dealing with large and small samples. This paper uses the radial basis kernel function, whose expression is Kxi,xj=exp−xi−xj2/2γ2, where γ denotes the width parameter of the kernel function. Different γ affects the performance of the relevance vector machine, so this parameter must be optimized, and the optimization algorithm based on metaheuristic can significantly improve the optimization efficiency.


## 3.2. Parameter Optimization Based on HHO

HHO is a novel metaheuristic swarm intelligence algorithm. Compared with the genetic algorithm and the particle swarm optimization algorithm, it has better exploration capabilities and convergence speed [26]. When the hawks catch their prey, they adopt different strategies in different stages. The HHO algorithm can find the optimal value by modeling hawks’ behavior in each stage of the predation process.The HHO algorithm includes the exploration phase and the exploitation phase. In the process of predation, with the increase in the moving distance of prey, the physical strength of prey will gradually decrease, which is called escape energy. The algorithm can transfer from the exploration phase to the exploitation phase through the value range of prey’s escape energy. The escape energy is modeled as [26]
(7)E=2E01−tT,where E denotes the escaping energy of the prey, T denotes the maximum number of iterations, t denotes the current number of iterations, and E0 denotes the initial energy. E0 randomly changes inside the interval −1,1 at each iteration.When∣E∣≥1, enter the exploration phase. The hawks inhabit different positions randomly and use two strategies to detect prey. An equal probability q is used to distinguish the two strategies. They perch based on the positions of other hawks and the prey when q<0.5 and perch on random tall trees when q≥0.5. The exploration phase can be expressed as
(8)Xt+1=Xrandt−r1Xrandt−2r2Xt,Xpreyt−Xmt−r3LB+r4UB−LB,q⩾0.5,q<0.5,where Xt+1 denotes the position vector of the next-generation hawks, Xt denotes the position vector of the current-generation hawks, Xpreyt denotes the position of the prey, r1,r2, r3, r4, and q are random values between 0,1 and are updated every generation, LB and UB denote the boundary of the variable, Xrandt denotes a random position selected from the current-generation hawks, and Xmt denotes the average position of the hawks.When∣E∣<1, enter the exploitation phase. At this time, use four siege strategies. r is a random number between 0,1 and indicates the chance for the prey to escape. r<0.5 indicates that the prey can escape, and r≥0.5 indicates that the prey cannot escape successfully [26]. The use of soft and hard besiege strategies depends on the prey’s energy E. When ∣E∣≥0.5, it means that the prey still has high energy. At this time, the Hawks adopt the soft besiege strategy, constantly approaching and chasing the prey to consume the energy of the prey. When ∣E∣<0.5, it indicates that the prey’s energy is weak. At this time, the Hawks can easily catch the prey by using the hard besiege strategy.(1) Soft Besiege. When r≥0.5 and ∣E∣≥0.5, it means that the prey has enough escape energy and jumps randomly. At this time, the hawks encircle it softly to exhaust the prey, expressed as
(9)Xt+1=Xpreyt−Xt−EJXpreyt−Xt,where J is a random number between 0,2, denoting the random jumping energy of prey, which is used to simulate the prey’s activity.(2) Hard Besiege. When r≥0.5 and ∣E∣<0.5, it means that the prey is weak and has low escape energy. At this time, the following equation is used to update the position.
(10)Xt+1=Xpreyt−EXpreyt−Xt.(3) Soft Besiege with Progressive Rapid Dives. When r<0.5 and ∣E∣≥0.5, the prey has enough energy to escape. At this time, equations (11) and (12) are used to update the position, and the Lévy flight is used to simulate the sudden movement of the hawks and the prey. The Lévy flight is expressed as equation (13), which can further enhance the ability of the algorithm to jump out of the local optimum.
(11)Y=Xpreyt−EJXpreyt−Xt,(12)Z=Y+S×LFD,(13)LFx=0.01×u×σv1/β,σ=Γ1+β×sinπβ/2Γ1+β/2×β×2β−1/21/β,where u and v are random values inside 0,1, β is a default constant set to 1.5, and Γz is the Gramma function: Γz=∫0∞tz−1e−tdt.Compare the two update methods, and select the most suitable position update method according to the fitness function value.(14)Xt+1=Y,ifFY<FXt,Z,ifFZ<FXt.(4) Hard Besiege with Progressive Rapid Dives. When r<0.5 and ∣E∣<0.5, it means that the prey does not have enough energy. At this time, equations (15) and (16) are used to update the position and equation (14) is used to select the appropriate update method.
(15)Y=Xpreyt−EJXpreyt−Xmt,(16)Z=Y+S×LFD.The relevance vector machine is used for data fitting and regression. HHO is used to optimize the kernel width parameter of the relevance vector machine.The flow of each HHO-RVR is as follows:Step 1.
Initialize the number of HHO populations and the number of iterations. The hawk’s position in the HHO algorithm represents the kernel width parameter of the relevance vector machine. The position of each hawk in the first-generation group is randomly distributed.Step 2.
Determine the fitness function. The sample data is randomly divided into five parts, 4 of which are used for the relevance vector machine training, and the remaining 1 part is used for verification, and the test accuracy of the relevance vector machine is used as the fitness function of HHO.Step 3.
Calculate the fitness value of each hawk, and use the hawk corresponding to the optimal fitness as the prey.Step 4.
Update the positions of the hawks, according to equations (7)–(16).Step 5.
Perform iterative operations. Repeat steps 3 and 4 until the number of iterations is satisfied and the calculation is stopped. At this time, the prey position is the optimal kernel width parameter.


## 3.3. RLG-IMU Temperature Modeling

The physical input of the RLG is the angular velocity, and the measured output is the number of pulses. The computer counts the number of pulses in the sampling time. According to the measurement equation of the system, the angular velocity of the RLG is obtained. The pulse-angular velocity measurement model of the RLGs can be expressed as(17)Ngτ=KgEgωgm+Dg.Ng=NgxNgyNgzT is the vector of the number of pulses output by the RLGs on three orthogonal axes within the sampling time τ. Kg=diagKgx,Kgy,Kgz is the scale factor matrix.Dg=DgxDgyDgzT is the bias matrix. ωgm=ωgxmωgymωgzmT is the angular velocity component vector of the three axes in the IMU measurement coordinate system.
(18)Eg=EgxxEgyxEgzxEgxyEgyyEgzyEgxzEgyzEgzzis the installation error matrix of the RLG relative inertial coordinate system for each axis.In equation (17), the output vector Ng of the RLGs changes with temperature and Eg is very little affected by temperature. Because the influence of temperature on the RLG is multifaceted, the change of the sensor’s output to the measured angular velocity caused by the temperature change is directly reflected in the scale factor matrix Kg and bias matrix Dg, and the structural changes will also be reflected on Kg and Dg, so equation (17) can be expressed as
(19)NgTτ=KgTEgωgm+DgT.The physical input of the quartz flexible accelerometer is the apparent acceleration, and the measured output is the number of pulses. The computer counts the number of pulses in the sampling time. According to the measurement equation of the system, the apparent acceleration of the quartz flexible accelerometer is obtained. The pulse-apparent acceleration measurement model of the quartz flexible accelerometers is expressed as(20)Naτ=KaEaAam+Da.Na=NaxNayNazT is the vector of the number of pulses output by the quartz flexible accelerometers on three orthogonal axes within the sampling time τ. Ka=diagKax,Kay,Kaz is the scale factor matrix of the quartz flexible accelerometers. Since the quartz flexible accelerometer is very sensitive to positive and negative scale factors, it is divided into positive and negative: the positive scale factor matrix is expressed as Kap=diagKapx,Kapy,Kapz and the negative scale factor matrix is expressed as Kan=diagKanx,Kany,Kanz. Da=DaxDayDazT is the bias matrix of the quartz flexible accelerometers. Aam=AaxmAaymAazmT is the apparent acceleration component vector of the three axes in the IMU measurement coordinate system.
(21)Ea=EaxxEayxEazxEaxyEayyEazyEaxzEayzEazzis the installation error matrix of the quartz flexible accelerometer relative inertial coordinate system for each axis. Similarly, equation (20) is expressed as
(22)NaTτ=KaTEaAam+DaT.In equations (19) and (22), the input parameter affected by temperature is called the temperature-related parameter, and there are 15 temperature-related parameters in total.Define a setP=Pi∣i=1⋯15 to denote 15 temperature-related parameters. Make Pi∣i=1⋯6 denote each element in Dg,Kg and Pi∣i=7⋯15 denote each element in Da,Kap,Kan.Since the temperature measurement position of each RLG and quartz flexible accelerometer is different, the temperature measurement curve of each temperature-related parameter will be different. Define a temperature setT=Ti∣i=1⋯15. Ti denotes the temperature change vector of each temperature-related parameter and corresponds one-to-one with Pi; then, the corresponding relationship between set P and temperature can be expressed as PT=PiTi∣i=1⋯15. From this, 15 models can be built as follows:
(23)Pi=Miwi,Ti,i=1,2,⋯,15,where Miwi,Ti is the model corresponding to each temperature-related parameter of the IMU and wi is the optimal width parameter of the kernel function of each model.


## 3.4. System Compensation Model

Ideally, in equations (19) and (22), angular velocity ωgm and apparent acceleration Aam are independent of temperature. First, the output models of angular velocity and apparent acceleration can be obtained according to IMU measurement equations (19) and (22) and the models of temperature-related parameters. The values of the temperature-related parameters at the desired temperature TE are then selected as constants in the angular/apparent acceleration measurement equation to obtain the compensation pulse. The desired temperature TE is a hypothetical stable working state. The purpose is that when the system works at any temperature, the measured value of physical quantity can be converted to the value at the desired temperature so that the system will not be affected by the temperature change. The desired temperature is not unique. It can be set according to the actual needs and requirements. No matter how much the desired temperature is set, the measured value of different temperatures can be converted to the measured value at the desired temperature by compensating for the measurement equation’s parameters to realize the temperature compensation.For the RLGs, according to equation (19), the output model of ωgm is
(24)ωgm=Eg−1NgTgτKg−1Tg−DgTg.Tg corresponding to each RLG is Tgx, Tgy, and Tgz.The compensation equation of RLG output pulse at the desired temperatureTE is
(25)NgCτ=KgTEEgωgm+DgTE.NgC=NgxCNgyCNgzCT is the output pulse of the three RLGs after compensation. KgTE and DgTE are, respectively, the scale factor matrix and bias matrix of the three RLGs at the desired temperature TE.By bringing equation (24) into (25), there is
(26)NgCτ=KgTEEgEg−1NgTgτKg−1Tg−DgTg+DgTE.Corresponding to different RLGs,Tg includes Tgx, Tgy, and Tgz. The data set Tg,Kg can be obtained by temperature experiment. The data set was used to train the relevance vector machine. At the same time, the optimization process introduced in Section 3.2 is used to optimize the kernel width parameter of the relevance vector machine, so as to obtain the optimized KgTg model. Similarly, the model of DgTg can be obtained. Since the temperature models of KgTg and DgTg have been obtained, it is possible to compensate for the output NgTg at all temperatures, thereby compensating for Ng to NgTE at different temperatures.For the quartz flexible accelerometers, according to equation (22), the output model of Aam is
(27)Aam=Ea−1NaTaτKa−1Ta−DaTa.Ta corresponding to each quartz flexible accelerometer is Tax, Tay, and Taz.The compensation equation of quartz flexible accelerometer output pulse at the desired temperatureTE is
(28)NaCτ=KaTEEaAam+DaTE.NaC=NaxCNayCNazCT is the output pulse of the three quartz flexible accelerometers after compensation. KaTE and DaTE are, respectively, the scale factor matrix and bias matrix of the three quartz flexible accelerometers at the desired temperature TE.By bringing equation (27) into (28), there is
(29)NaCτ=KaTEEaEa−1NaTaτKa−1Ta−DaTa+DaTE.Corresponding to different quartz flexible accelerometers,Ta includes Tax, Tay, and Taz. Similarly, the models of KaTa and DaTa can be obtained. Since the temperature models of KaTa and DaTa have been obtained, it is possible to compensate for the output NaTa at all temperatures, thereby compensating for NaTa to NaTE at different temperatures.Equations (26) and (29) are the compensation model of the IMU.


## 4. Experiment and Analysis

### 4.1. Data Acquisition and Modeling

The IMU used in the experiment uses three RLGs and three quartz flexible accelerometers as measurement sensors. The scale factor of the selected RLG on the IMU system is approximately 2.14 pulses/arcsec. The scale factor of the selected accelerometer on the IMU system is approximately2.4e+4pulses/s.g0. The RLG-IMU is fixed on a two-axis rotating platform in the temperature chamber. Set the temperature chamber to -5°C, and keep it for 4 hours to ensure the temperature inside the IMU and the chamber temperature are consistent. The temperature chamber is set to 7000 minutes from -5°C to 65°C. In the process of temperature rise in the temperature chamber, the IMU continuously carries out rotation calibration work in the temperature chamber. The IMU uses the systematic calibration method, and the rotation speed is 20 deg/s. Various errors will be excited by the continuous rotation of multiple positions. Using the Kalman filter, various error parameters of the RLGs and accelerometers can be calculated, that is, Kg, Dg, Eg, Ka, Da, and Ea in Equations (17) and (20). Each calibration takes approximately 40 minutes. A total of 175 calibrations have been completed within the set temperature range and time. The IMU temperature change is caused not only by the temperature change in the temperature chamber but also by the heating part inside the IMU. The temperature sensor measurement shows that the temperature change between every two calibrations is about 0.4°C. The temperature change is tiny. Therefore, the temperature is considered to be stable in one calibration process. Through the experimental process described above, 175 sets of calibration data of the IMU can be obtained over the entire experimental temperature range. After offline calibration data processing, 175 sets of temperature-related parameters at 175 temperature points can be obtained.Through experiments in the temperature chamber, the data setPi of each temperature-related parameter and the corresponding temperature set Ti are obtained. There are 15 such sets in total. It can be expressed as S=Sii=115,Si=Tij,Pijj=1175. Each Si is used as the relevance vector machine training data. Use the HHO-RVR process in Section 3.2 to train each relevance vector machine. The kernel width parameter of the relevance vector machine is used as the input variable of the HHO. The total population size is set to 30, the max number of iteration is set to 500, and the value range is set to 0.01,200. The kernel width parameter with the highest prediction accuracy in the range of values is found using the 5-fold cross-validation method. Each temperature-related parameter is modeled according to this process. The set of models M=Mii=115 in equation (23) can be obtained.For comparison, the least square method is used to model the 15 temperature-related parameters. The model equation is expressed as(30)P=AT,where P=P1tP2t⋯P15tT is the output vector of the temperature-related parameters,
(31)A=a1,0a1,1a1,2a1,3a2,0a2,1a2,2a2,3⋯⋯⋯⋯a15,0a15,1a15,2a15,3is the coefficient matrix of the equation, and Τ=1tt2t3T is the temperature input vector.Fifteen temperature-related parameters were modeled by the RVR method and the least square method, respectively. In the RVR model, because the radial basis function is used as the kernel function, there is only one kernel width parameter. According to equation (30), the least square method of a cubic fitting polynomial is used. Each fitting polynomial needs to determine four coefficients: the row vector of matrix A in equation (30). The kernel width parameters, the number of relevance vectors, the proportion of relevance vectors of the 15 RVR models, and the coefficients of 15 least square fitting polynomials are listed in Table 1.Table 1
Model parameters.
The temperature-related parameterParameters of the RVR modelsCoefficients of the least square methodKernel function width parameter (γ)Number of relevance vectorsThe proportion of relevance vectors (%)a0a1a2a3P13231.71−2.11E−08−2.17E−063.69E−042.15E−02P24.4174−1.10E−075.82E−06−7.12E−055.32E−02P310.6952.86−1.06E−078.23E−06−8.18E−05−5.50E−02P46.6752.863.64E−11−2.88E−099.01E−082.14E+00P511.2352.862.21E−10−2.26E−088.34E−072.14E+00P625.6521.14−1.48E−112.02E−094.57E−082.14E+00P73252.864.68E−10−6.57E−071.67E−04−1.13E−02P87.7184.573.72E−09−2.46E−074.04E−066.17E−03P911.0895.143.10E−09−6.43E−076.02E−057.01E−03P1028.9252.86−2.50E−051.11E−02−1.16E+002.40E+04P1119.4274−2.51E−051.24E−02−4.87E−012.42E+04P1227.7452.86−1.32E−059.90E−03−9.13E−012.38E+04P1331.0963.43−2.55E−051.12E−02−1.16E+002.40E+04P1419.474−2.55E−051.25E−02−4.93E−012.42E+04P1515.4284.57−1.59E−051.02E−02−9.19E−012.38E+04


### 4.2. Analysis of the Regression Performance of the Models

To compare the fitting and regression performance of HHO-RVR and LSM, HHO-RVR and LSM were used to model the temperature value of 15 temperature-related parameters, with a total of 30 models. Three indicators are used to measure the fitting and regression performance of the two modeling methods:(1) Root Mean Square Error [36]. It is used to measure the error and dispersion between the regression value and the real value. The smaller the value, the smaller the estimation error and the smaller the error dispersion. The equation is as follows:
(32)fX,Φ=1m∑i=1mΦxi−yi2.(2) Mean Absolute Error [37]. It is used to measure the average absolute errors. Reflect the actual situation of the predicted value error. The smaller the value, the smaller the estimated error. The equation is as follows:
(33)fX,Φ=1m∑i=1mΦxi−yi.(3) R-Squared. It is used to measure the fitting accuracy of the model; the closer to 1, the higher the fitting accuracy. The equation is as follows:
(34)fX,Φ=1−∑i=1mΦxi−yi2∑i=1mymean−yi2.In equations (32), (33), and (34), X=xi∣i=1⋯m is the sample input, Φ is the regression model, Φxi is the regression output, yi is the sample output value, ymean is the mean of the sample output, and m is the number of samples.The results of the three indicators are shown in Table2. Compared with the least square method, all the regression test performances of the 15 models were improved by using the HHO-RVR method. The root mean square error value increased by a maximum of 80.12% and an average of 35.27%. The mean absolute error value increased by a maximum of 83.51% and an average of 39.29%. The R-squared value increased by a maximum of 9.15% and an average of 2.29%. The results show that the RVR has better performance than the least square method in terms of error, error dispersion, and fitting accuracy.Table 2
Performance comparison.
The temperature-related parameterRoot mean square errorMean absolute errorR-squaredThe least square methodHHO-RVRImprovement (%)The least square methodHHO-RVRImprovement (%)The least square methodHHO-RVRImprovement (%)P16.1037E−045.7205E−046.284.7865E−044.5685E−044.569.6287E−019.6739E−010.47P28.5874E−046.7484E−0421.426.8025E−045.4583E−0419.769.0215E−019.3957E−014.15P36.6615E−045.4146E−0418.725.3673E−044.1850E−0422.037.8803E−018.5995E−019.13P49.9417E−079.7066E−072.377.9998E−077.6307E−074.614.3971E−014.6591E−015.96P51.4264E−061.0149E−0628.851.0347E−067.8771E−0723.879.0700E−019.5292E−015.06P68.3064E−078.2362E−070.846.4871E−076.4206E−071.038.9077E−018.9261E−010.21P71.4118E−058.0060E−0643.291.0219E−056.3037E−0638.319.9997E−019.9999E−010.002P81.9559E−054.3197E−0677.911.5416E−053.3696E−0678.149.1229E−019.9572E−019.15P92.3627E−055.1660E−0678.141.9734E−054.0025E−0679.729.9843E−019.9992E−010.15P108.5537E−027.4064E−0213.417.0884E−025.5104E−0222.269.9995E−019.9996E−010.001P118.6871E−021.7273E−0280.127.5096E−021.2384E−0283.519.9975E−019.9999E−010.02P126.4535E−025.9523E−027.775.2699E−023.5745E−0232.179.9991E−019.9992E−010.001P136.7920E−024.1029E−0239.595.5552E−022.6900E−0251.589.9997E−019.9999E−010.002P147.5613E−021.9322E−0274.456.5216E−021.3913E−0278.679.9980E−019.9999E−010.02P155.8027E−023.7193E−0235.904.8269E−022.4574E−0249.099.9993E−019.9997E−010.004As shown in Figure2, the original data is plotted against the data calculated by regression using HHO-RVR and LSM models, respectively. Observe the goodness of fit between the regression curve using both methods and the original data curve. The horizontal axis is temperature, and the vertical axis is the value of each temperature-related parameter. Since some of the curves do not show any significant differences from the graph, nine representative curves are selected for comparison. According to the order of the figures, they are Dax (the bias of the accelerometer x), Daz (the bias of the accelerometer z), Kany (the negative scale factor of the accelerometer y), Kgx (the scale factor of the gyro x), Dgz (the bias of the gyro z), Dgx (the bias of the gyro x), Dgy (the bias of the gyro y), Kgy (the scale factor of the gyro y), and Day (the bias of the accelerometer y). Figures Dax, Daz and Kany show that when the original data trend is relatively simple and the fluctuation and dispersion are small, the performance difference between the fitting accuracy of the two methods is not very obvious. Figures Kgx, Dgz, Dgx, and Dgy show that the data are characterized by large dispersion, obvious fluctuations, and complex nonlinearity. In this case, the curve fitted by the RVR can accurately fit the complex fluctuations of the original data. However, the least square method can only fit the trend of the whole data and cannot accurately represent the complex fluctuation of the original data, thus affecting the performance of regression and prediction. Figures Kgy and Day show that when the original data dispersion is small, but the data change trend is obvious, the RVR fits the data more accurately than the least square method.Figure 2
Comparison curve of two regression models.
(a)
Dax(b)
Daz(c)
Kany(d)
Kgx(e)
Dgz(f)
Dgx(g)
Dgy(h)
Kgy(i)
Day


### 4.3. System Dynamic Compensation

In order to observe the system’s dynamic performance after temperature compensation, the IMU dynamic data is used for verification, which can better characterize the compensation performance when the IMU is in motion.Take 50°C as the desired compensation temperature. First, the temperature compensation model of the temperature-related parameters is obtained in Section4.1, and the gyro pulse data at the corresponding temperature are brought into equations (26) and (29) to get the compensated pulse data at all temperature points. Secondly, the compensated data are used for calibration operations to obtain the temperature-related parameters of compensated dynamic data. Finally, each temperature-related parameter calculated after compensation is subtracted from the temperature-related parameter at the desired temperature. Two methods are used to model and regress the temperature-related parameters, and the temperature-related parameters are compensated to the desired temperature. The compensated value is subtracted from the measured value to obtain the compensation residual. The smaller the fluctuation of the residual curve, the smaller the residual variance is and the higher the compensation accuracy will be.Figure3 shows the residuals of 15 temperature-related parameters compensated by RVR and the least square method. Because the magnitudes of Kgx, Kgy, Kgz, Dgx, Dgy, and Dgz are small and the trend is not apparent, the two curves overlap more. From the curve of 9 parameters of Kapx, Kapy, Kapz, Kanx, Kany, Kanz, Dax, Day, and Daz, it can be seen that the temperature-related parameters compensated by the RVR method have smaller fluctuations and are more stable than those compensated by the least square method. Figure Day is a typical curve. It can be seen that the residual curve after the LSM compensation has a large fluctuation. The residual curve of the HHO-RVR compensation is very stable. It shows that the original data has a strong nonlinearity. The fitting performance of the LSM is not as good as that of the HHO-RVR, which leads to a large fluctuation of residual after compensation. Other similar figures also illustrate this conclusion.Figure 3
Comparison of residuals after the temperature-related parameter compensation.
(a)
Kgx(b)
Kgy(c)
Kgz(d)
Dgx(e)
Dgy(f)
Dgz(g)
Kapx(h)
Kapy(i)
Kapz(j)
Kanx(k)
Kany(l)
Kanz(m)
Dax(n)
Day(o)
DazAs shown in Table3, the variances of the residuals are compared in a table. It can be seen that the residual variance after compensation using the HHO-RVR method is lower than that using the least square method, indicating that the performance of using RVR to compensate for the temperature-related parameters is more superior.Table 3
Variance comparison data for residuals.
The temperature-related parameterOriginal dataHHO-RVRThe least square methodImprovement (%)P13.16773E−036.42336E−046.94344E−047.49P22.74528E−038.32613E−049.54550E−0412.77P31.44687E−035.85158E−047.18217E−0418.53P41.32818E−061.72603E−061.73323E−060.42P54.67740E−061.91785E−062.11570E−069.35P62.51331E−062.00860E−062.03553E−061.32P72.56482E−039.96659E−061.71954E−0542.04P86.60438E−055.72778E−062.01836E−0571.62P95.95731E−046.61036E−062.44531E−0572.97P101.19157E+017.68137E−028.15592E−025.82P115.44946E+002.00643E−028.83737E−0277.30P126.83441E+006.02220E−026.64657E−029.39P131.19019E+014.94948E−026.28477E−0221.25P145.39584E+002.22470E−027.57132E−0270.62P156.72286E+003.89916E−025.99048E−0234.91


## 4.1. Data Acquisition and Modeling

The IMU used in the experiment uses three RLGs and three quartz flexible accelerometers as measurement sensors. The scale factor of the selected RLG on the IMU system is approximately 2.14 pulses/arcsec. The scale factor of the selected accelerometer on the IMU system is approximately2.4e+4pulses/s.g0. The RLG-IMU is fixed on a two-axis rotating platform in the temperature chamber. Set the temperature chamber to -5°C, and keep it for 4 hours to ensure the temperature inside the IMU and the chamber temperature are consistent. The temperature chamber is set to 7000 minutes from -5°C to 65°C. In the process of temperature rise in the temperature chamber, the IMU continuously carries out rotation calibration work in the temperature chamber. The IMU uses the systematic calibration method, and the rotation speed is 20 deg/s. Various errors will be excited by the continuous rotation of multiple positions. Using the Kalman filter, various error parameters of the RLGs and accelerometers can be calculated, that is, Kg, Dg, Eg, Ka, Da, and Ea in Equations (17) and (20). Each calibration takes approximately 40 minutes. A total of 175 calibrations have been completed within the set temperature range and time. The IMU temperature change is caused not only by the temperature change in the temperature chamber but also by the heating part inside the IMU. The temperature sensor measurement shows that the temperature change between every two calibrations is about 0.4°C. The temperature change is tiny. Therefore, the temperature is considered to be stable in one calibration process. Through the experimental process described above, 175 sets of calibration data of the IMU can be obtained over the entire experimental temperature range. After offline calibration data processing, 175 sets of temperature-related parameters at 175 temperature points can be obtained.Through experiments in the temperature chamber, the data setPi of each temperature-related parameter and the corresponding temperature set Ti are obtained. There are 15 such sets in total. It can be expressed as S=Sii=115,Si=Tij,Pijj=1175. Each Si is used as the relevance vector machine training data. Use the HHO-RVR process in Section 3.2 to train each relevance vector machine. The kernel width parameter of the relevance vector machine is used as the input variable of the HHO. The total population size is set to 30, the max number of iteration is set to 500, and the value range is set to 0.01,200. The kernel width parameter with the highest prediction accuracy in the range of values is found using the 5-fold cross-validation method. Each temperature-related parameter is modeled according to this process. The set of models M=Mii=115 in equation (23) can be obtained.For comparison, the least square method is used to model the 15 temperature-related parameters. The model equation is expressed as(30)P=AT,where P=P1tP2t⋯P15tT is the output vector of the temperature-related parameters,
(31)A=a1,0a1,1a1,2a1,3a2,0a2,1a2,2a2,3⋯⋯⋯⋯a15,0a15,1a15,2a15,3is the coefficient matrix of the equation, and Τ=1tt2t3T is the temperature input vector.Fifteen temperature-related parameters were modeled by the RVR method and the least square method, respectively. In the RVR model, because the radial basis function is used as the kernel function, there is only one kernel width parameter. According to equation (30), the least square method of a cubic fitting polynomial is used. Each fitting polynomial needs to determine four coefficients: the row vector of matrix A in equation (30). The kernel width parameters, the number of relevance vectors, the proportion of relevance vectors of the 15 RVR models, and the coefficients of 15 least square fitting polynomials are listed in Table 1.Table 1
Model parameters.
The temperature-related parameterParameters of the RVR modelsCoefficients of the least square methodKernel function width parameter (γ)Number of relevance vectorsThe proportion of relevance vectors (%)a0a1a2a3P13231.71−2.11E−08−2.17E−063.69E−042.15E−02P24.4174−1.10E−075.82E−06−7.12E−055.32E−02P310.6952.86−1.06E−078.23E−06−8.18E−05−5.50E−02P46.6752.863.64E−11−2.88E−099.01E−082.14E+00P511.2352.862.21E−10−2.26E−088.34E−072.14E+00P625.6521.14−1.48E−112.02E−094.57E−082.14E+00P73252.864.68E−10−6.57E−071.67E−04−1.13E−02P87.7184.573.72E−09−2.46E−074.04E−066.17E−03P911.0895.143.10E−09−6.43E−076.02E−057.01E−03P1028.9252.86−2.50E−051.11E−02−1.16E+002.40E+04P1119.4274−2.51E−051.24E−02−4.87E−012.42E+04P1227.7452.86−1.32E−059.90E−03−9.13E−012.38E+04P1331.0963.43−2.55E−051.12E−02−1.16E+002.40E+04P1419.474−2.55E−051.25E−02−4.93E−012.42E+04P1515.4284.57−1.59E−051.02E−02−9.19E−012.38E+04


## 4.2. Analysis of the Regression Performance of the Models

To compare the fitting and regression performance of HHO-RVR and LSM, HHO-RVR and LSM were used to model the temperature value of 15 temperature-related parameters, with a total of 30 models. Three indicators are used to measure the fitting and regression performance of the two modeling methods:(1) Root Mean Square Error [36]. It is used to measure the error and dispersion between the regression value and the real value. The smaller the value, the smaller the estimation error and the smaller the error dispersion. The equation is as follows:
(32)fX,Φ=1m∑i=1mΦxi−yi2.(2) Mean Absolute Error [37]. It is used to measure the average absolute errors. Reflect the actual situation of the predicted value error. The smaller the value, the smaller the estimated error. The equation is as follows:
(33)fX,Φ=1m∑i=1mΦxi−yi.(3) R-Squared. It is used to measure the fitting accuracy of the model; the closer to 1, the higher the fitting accuracy. The equation is as follows:
(34)fX,Φ=1−∑i=1mΦxi−yi2∑i=1mymean−yi2.In equations (32), (33), and (34), X=xi∣i=1⋯m is the sample input, Φ is the regression model, Φxi is the regression output, yi is the sample output value, ymean is the mean of the sample output, and m is the number of samples.The results of the three indicators are shown in Table2. Compared with the least square method, all the regression test performances of the 15 models were improved by using the HHO-RVR method. The root mean square error value increased by a maximum of 80.12% and an average of 35.27%. The mean absolute error value increased by a maximum of 83.51% and an average of 39.29%. The R-squared value increased by a maximum of 9.15% and an average of 2.29%. The results show that the RVR has better performance than the least square method in terms of error, error dispersion, and fitting accuracy.Table 2
Performance comparison.
The temperature-related parameterRoot mean square errorMean absolute errorR-squaredThe least square methodHHO-RVRImprovement (%)The least square methodHHO-RVRImprovement (%)The least square methodHHO-RVRImprovement (%)P16.1037E−045.7205E−046.284.7865E−044.5685E−044.569.6287E−019.6739E−010.47P28.5874E−046.7484E−0421.426.8025E−045.4583E−0419.769.0215E−019.3957E−014.15P36.6615E−045.4146E−0418.725.3673E−044.1850E−0422.037.8803E−018.5995E−019.13P49.9417E−079.7066E−072.377.9998E−077.6307E−074.614.3971E−014.6591E−015.96P51.4264E−061.0149E−0628.851.0347E−067.8771E−0723.879.0700E−019.5292E−015.06P68.3064E−078.2362E−070.846.4871E−076.4206E−071.038.9077E−018.9261E−010.21P71.4118E−058.0060E−0643.291.0219E−056.3037E−0638.319.9997E−019.9999E−010.002P81.9559E−054.3197E−0677.911.5416E−053.3696E−0678.149.1229E−019.9572E−019.15P92.3627E−055.1660E−0678.141.9734E−054.0025E−0679.729.9843E−019.9992E−010.15P108.5537E−027.4064E−0213.417.0884E−025.5104E−0222.269.9995E−019.9996E−010.001P118.6871E−021.7273E−0280.127.5096E−021.2384E−0283.519.9975E−019.9999E−010.02P126.4535E−025.9523E−027.775.2699E−023.5745E−0232.179.9991E−019.9992E−010.001P136.7920E−024.1029E−0239.595.5552E−022.6900E−0251.589.9997E−019.9999E−010.002P147.5613E−021.9322E−0274.456.5216E−021.3913E−0278.679.9980E−019.9999E−010.02P155.8027E−023.7193E−0235.904.8269E−022.4574E−0249.099.9993E−019.9997E−010.004As shown in Figure2, the original data is plotted against the data calculated by regression using HHO-RVR and LSM models, respectively. Observe the goodness of fit between the regression curve using both methods and the original data curve. The horizontal axis is temperature, and the vertical axis is the value of each temperature-related parameter. Since some of the curves do not show any significant differences from the graph, nine representative curves are selected for comparison. According to the order of the figures, they are Dax (the bias of the accelerometer x), Daz (the bias of the accelerometer z), Kany (the negative scale factor of the accelerometer y), Kgx (the scale factor of the gyro x), Dgz (the bias of the gyro z), Dgx (the bias of the gyro x), Dgy (the bias of the gyro y), Kgy (the scale factor of the gyro y), and Day (the bias of the accelerometer y). Figures Dax, Daz and Kany show that when the original data trend is relatively simple and the fluctuation and dispersion are small, the performance difference between the fitting accuracy of the two methods is not very obvious. Figures Kgx, Dgz, Dgx, and Dgy show that the data are characterized by large dispersion, obvious fluctuations, and complex nonlinearity. In this case, the curve fitted by the RVR can accurately fit the complex fluctuations of the original data. However, the least square method can only fit the trend of the whole data and cannot accurately represent the complex fluctuation of the original data, thus affecting the performance of regression and prediction. Figures Kgy and Day show that when the original data dispersion is small, but the data change trend is obvious, the RVR fits the data more accurately than the least square method.Figure 2
Comparison curve of two regression models.
(a)
Dax(b)
Daz(c)
Kany(d)
Kgx(e)
Dgz(f)
Dgx(g)
Dgy(h)
Kgy(i)
Day


## 4.3. System Dynamic Compensation

In order to observe the system’s dynamic performance after temperature compensation, the IMU dynamic data is used for verification, which can better characterize the compensation performance when the IMU is in motion.Take 50°C as the desired compensation temperature. First, the temperature compensation model of the temperature-related parameters is obtained in Section4.1, and the gyro pulse data at the corresponding temperature are brought into equations (26) and (29) to get the compensated pulse data at all temperature points. Secondly, the compensated data are used for calibration operations to obtain the temperature-related parameters of compensated dynamic data. Finally, each temperature-related parameter calculated after compensation is subtracted from the temperature-related parameter at the desired temperature. Two methods are used to model and regress the temperature-related parameters, and the temperature-related parameters are compensated to the desired temperature. The compensated value is subtracted from the measured value to obtain the compensation residual. The smaller the fluctuation of the residual curve, the smaller the residual variance is and the higher the compensation accuracy will be.Figure3 shows the residuals of 15 temperature-related parameters compensated by RVR and the least square method. Because the magnitudes of Kgx, Kgy, Kgz, Dgx, Dgy, and Dgz are small and the trend is not apparent, the two curves overlap more. From the curve of 9 parameters of Kapx, Kapy, Kapz, Kanx, Kany, Kanz, Dax, Day, and Daz, it can be seen that the temperature-related parameters compensated by the RVR method have smaller fluctuations and are more stable than those compensated by the least square method. Figure Day is a typical curve. It can be seen that the residual curve after the LSM compensation has a large fluctuation. The residual curve of the HHO-RVR compensation is very stable. It shows that the original data has a strong nonlinearity. The fitting performance of the LSM is not as good as that of the HHO-RVR, which leads to a large fluctuation of residual after compensation. Other similar figures also illustrate this conclusion.Figure 3
Comparison of residuals after the temperature-related parameter compensation.
(a)
Kgx(b)
Kgy(c)
Kgz(d)
Dgx(e)
Dgy(f)
Dgz(g)
Kapx(h)
Kapy(i)
Kapz(j)
Kanx(k)
Kany(l)
Kanz(m)
Dax(n)
Day(o)
DazAs shown in Table3, the variances of the residuals are compared in a table. It can be seen that the residual variance after compensation using the HHO-RVR method is lower than that using the least square method, indicating that the performance of using RVR to compensate for the temperature-related parameters is more superior.Table 3
Variance comparison data for residuals.
The temperature-related parameterOriginal dataHHO-RVRThe least square methodImprovement (%)P13.16773E−036.42336E−046.94344E−047.49P22.74528E−038.32613E−049.54550E−0412.77P31.44687E−035.85158E−047.18217E−0418.53P41.32818E−061.72603E−061.73323E−060.42P54.67740E−061.91785E−062.11570E−069.35P62.51331E−062.00860E−062.03553E−061.32P72.56482E−039.96659E−061.71954E−0542.04P86.60438E−055.72778E−062.01836E−0571.62P95.95731E−046.61036E−062.44531E−0572.97P101.19157E+017.68137E−028.15592E−025.82P115.44946E+002.00643E−028.83737E−0277.30P126.83441E+006.02220E−026.64657E−029.39P131.19019E+014.94948E−026.28477E−0221.25P145.39584E+002.22470E−027.57132E−0270.62P156.72286E+003.89916E−025.99048E−0234.91


## 5. Conclusion

In order to compensate for the influence of temperature on the RLG-IMU system and reduce the system startup stability time, an RLG-IMU system-level temperature compensation method is proposed. The RVR was used to model the 15 temperature-related parameters of the RLG-IMU, and the HHO algorithm was used to optimize the width parameters of the RVR kernel function. A system compensation model was established. The test results show that the fitting accuracy and compensation accuracy are more excellent than those of the least square method. In terms of fitting and regression performance, the root mean square error value increased by an average of 35.27%, the mean absolute error value increased by an average of 39.29%, and theR-squared value increased by an average of 2.29%. In terms of system compensation performance, the variance of the residual after compensation increased by 30.34% on average. The system-level compensation method can effectively reduce the startup stability time of the RLG-IMU, which has high application value and practical engineering value.



---

*Source: 6613574-2021-02-13.xml*