An Improved Online Self-Calibration Method Utilizing Angular Velocity Observation for Ultra High Accuracy PIGA-Based IMU

In the field of ultra high accuracy inertial measurement unit (IMU), pendulous integrating gyroscopic accelerometer (PIGA) has become a research hot spot due to its high-end performance. However, PIGA is sensitive to angular velocity, and the calibration process of PIGA-based IMU will be very complicated, which makes online self-calibration difficult to implement. To solve the above problems, we proposed an online self-calibration method utilizing angular velocity observation. The main contributions of this study are twofold: (1) An error analysis of PIGA is conducted in this paper, and the error model has also been simplified to suit the self-calibration model. (2) An improved online self-calibration method utilizing angular observation based on a simplified PIGA error model is proposed in this study. Experimental results show that the self-calibration method proposed in this study can improve the PIGA online calibration accuracy effectively (with the accuracy within 0.02 m/s/pulse), which can improve the dynamic accuracy of the PIGA.


Introduction
In the past few decades, pendulous integrating gyroscopic accelerometer (PIGA) has become a research hot spot due to its high-end performance in the field of ultra high accuracy inertial measurement unit (IMU) [1]. The accuracy of gyroscopes has made great progress in the past few decades [2]. However, the accuracy of the accelerometer has a very high correlation between the horizontal attitude and the positioning accuracy of the inertial navigation system (INS), especially in free inertial model [3]. PIGA has many advantages, such as large overload, high accuracy, and high measure range [4]. After a complicated calibration process, PIGA-based IMU can reach high navigation performance [5], which lays the foundation for its utilization in ultra high accuracy INS.
However, PIGA is sensitive to angular velocity, and the calibration process of PIGAbased IMU will be very complicated, which makes online self-calibration difficult to implement. Many researchers have studied the self-calibration method for rotation INS (RINS) and hybrid INS (HINS). In terms of biases estimation of IMU, Fong et al. proposed methods for in-field calibration of IMU without external equipment [6]. In Ref. [7], a gyro-bias calibration method has been proposed for analytic coarse alignment. Han et al. proposed a method for bias calibration of ring laser gyroscopes (RLGs) in single-axis RINS [8]. Furthermore, Li et al. in [9] analyzed the observability of the IMU bias selfcalibration method of single-axis RINS. When the rotational degrees of freedom exceed 1, the calibration of IMU's scale factors can be realized [10]. PIGA contains the angular velocity coupling calibration parameters, with only one rotation axis (only for bias calibration) that cannot meet the self-calibration demand of PIGA.

Kinetics Analysis of PIGA
The internal coordinate system of PIGA is shown as Figure 1: In Figure 1,α denotes the angular velocity of the outer frame assembly about the oy 1 axis,β is the angular velocity of the inner frame assembly about the ox axis relative to the outer frame.φ denotes the angular velocity of the rotor assembly about its axis of rotation relative to the inner frame.
The establishment process of the PIGA internal coordinate system can be expressed as follows: First, it is assumed that the shell coordinate system ox 0 y 0 z 0 of the gyro accelerometer, the outer frame coordinate system ox 1 y 1 z 1 , and the inner frame coordinate system oxyz are completely coincident, and the center of mass of the eccentric mass is located on the z-axis. Secondly, it is assumed that the outer ring component of PIGA is rotated around the axis by an angle, the inner ring component is rotated by α angle around the axis, and the non-orthogonal state of the PIGA internal coordinate system is obtained as shown in the above figure. Under the condition that the rotor is fully dynamically balanced, each axis of the coordinate system ox R y R z R is the main axis of inertia of the rotor, so its inertia product is zero, and the calculation formula of the angular momentum of the rotor is as follows: (1) A = B is the moment of inertia about the two coordinate axes perpendicular to the rotor axis. C is the moment of inertia of the rotor shaft. Since the ox, oy, and oz axes are not the main inertial axes of the inner frame assembly, their inertia products are not equal to zero. The formula for calculating the angular momentum of the inner frame assembly (excluding the rotor) relative to the oxyz axis is: where J xx , J yy , J zz are the moment of inertia of the frame component to the coordinate system oxyz. J xy , J xz and J yz are the inertia product of the inner frame component on the planes to which the x-axis and the y-axis, the x-axis and the z-axis, and the y-axis and the z-axis belong. M x is the total external torque acting on the inner ring shaft, including elastic torque, friction torque, electromagnetic interference torque, etc., which are not related to specific force, and inertial torque, unequal elastic torque, damping torque, etc., which are related to specific force. M x can be described as: where M xr is the sum of the disturbance moments of the inner frame.
Assuming that the coordinates of the center of mass of the PIGA inner ring component are (0, 0, l), the moment generated by the acceleration can be obtained, and M xr can be analyzed and calculated below.
where m is the mass of the inner frame assembly. The expressions of the input-specific force of PIGA in the inner and outer ring coordinate systems can be expressed as: According to the above formula, the inertia moment caused by the input-specific force is:

Simplified PIGA Error Model
After the above analysis of various error sources, to accurately calibrate the nonlinear error term of the PIGA, it is necessary to accurately model the gyroscopes.
It can be seen from the working principle of PIGA that the motion equation of the inner frame shaft of PIGA is deduced according to the Euler equation dH dt + ω × H = M. The precession angle α is much larger than the misalignment angle β, and β is a small angle. Now, let I x = J xx + A, I y = J yy + B, I z = J zz R , H z R = H, and simplify processing, Then cos β = cos 2β = 1, sin β = β, sin 2β = 2β, deriving the output equation of PIGA can obtain:α = ml H a y 1 + ml H β(a x 0 sin α + a z 0 cos α) − ∆φ z H (a x 0 cos α − a z 0 sin α) Since the output equation of PIGA in the ideal state can be expressed asα = ml H a y 1 , that is, in the above formula,α = ml H a y 1 is an effective signal, and the rest can be regarded as error terms, and its physical meaning is as follows: M m H has nothing to do with the input specific force, mainly the error term caused by elastic torque, friction torque and electromagnetic interference torque caused by electromagnetic components. ω y 0 + ω z 1 β is the error term due to the involved motion of the gyro shell. ∆φ z and ∆φ y are the error terms caused by the PIGA outer ring shaft and the inner ring shaft, as well as the out system between the inner ring shaft and the rotor, which are affected by machining and assembly errors.α 2 is the nonlinear error of the gyroscopic accelerometer, which is caused by the centrifugal force of the inner ring assembly as it rotates around the outer ring axis.
Based on Equation (7), the general error model equation of PIGA can be obtained as follows:α = K 0 + K x a x + K y a y + K z a z + K xx a 2 x + K yy a 2 y + K zz a 2 z + K xy a x a y +K xz a x a z + K yz a y a z +∆ xωx +∆ yωy +∆ zωz where K is the static error coefficient, ∆ denotes the dynamic error coefficient, δ is the mixed error coefficient.
The above PIGA error model expression includes the static error model of PIGA, which is the term related to the linear acceleration of PIGA; the dynamic error model, which is the term related to the input angular velocity and angular acceleration; the mixed error model, which is related to the linear acceleration and PIGA.
Since the input shaft and output shaft of PIGA are coincident, the corresponding terms of the cross-axis x-axis and z-axis have little effect on the output. The simplified model of Equation (8) can be written as: It can be seen from the above analysis that K 0 is introduced by unequal elastic torque, friction torque, damping torque, and other torques that have nothing to do with the input specific force; K y is caused by the change of the meter parameters, and the change includes the detection quality and the motor rotor quality. Changes, pendulum length changes, changes in the inertia radius of the motor rotor moment of inertia, and changes in the motor speed, among which the change in mass is small and generally ignored; the quadratic term coefficient K yy is proposed due to the existence ofα 2 , and they are determined by the inertia of the inner frame component. Product, unequal inertia and rotational inertia, etc., are caused by the centrifugal force of the inner frame shaft due to the rotation of the gyro accelerometer around the outer frame shaft under the condition of the verticality error; K ω is caused by the rotation of the output shaft. The error caused by the angular velocity coupling term. Since the coefficient K yy of the gyro accelerometer error model is affected by the interference torque existing on the inner frame axis, when machining the gyro accelerometer, the interference torque existing on the inner frame axis can be utilized to reduce its impact on the scale factors.

Error Model of PIGAs in Filtering
As the core component of the strategic weapon IMU system, the gyro accelerometer's calibration level directly affects and determines the actual combat effectiveness of the weapon [18]. In the actual system use, when the rocket engine is turned off during the transition from power flight to ballistic flight, the gyro accelerometer inertial group is required to control the cut-off speed and initial pitch, and heading of the missile. Setup, simulation analysis, etc., consist of only two PIGAs and a quartz accelerometer.
The IMU accelerometer inertial group in this section consists of two gyro accelerometers and a quartz accelerometer (hereinafter referred to as the accelerometer component). The calibration parameters include the accelerometer bias, scale factor, installation error angle, quadratic term coefficient, and angle Rate coupling term coefficients.
Considering that the accelerometers in the IMU are not strictly orthogonally installed, in the inertial navigation solution, the output of each inertial device must be projection is in the same orthogonal coordinate system, which is named the IMU coordinate system (m-frame) in this paper.
In this work, to convert the measurement information of the accelerometer from the oblique coordinate system to the IMU coordinate system, an orthogonal coordinate system (o-frame) is defined by the sensitive axis of the accelerometer. The x o axis is consistent with the sensitive axis x a of the accelerometer, and the y o axis is in the plane formed by the sensitive axes x a and y a of the accelerometer and differs from the y a axis by a small angle β yz , the z o axis is the accelerometer sensitive axis z a rotated by a small angle β zx around the x a axis, and then rotated by a small angle β zy around the y a axis, as shown in Figure 2. The installation error matrix of the accelerometer can be expressed as: where β ij represents the deflection angle of the i-axis of the accelerometer sensitive axis around the j-axis of the o-frame.
Since the installation error between the gyro component and the accelerometer component is a small angle, the coordinate transformation matrix from the system to the IMU coordinate system can be written as: where the installation error angle η i between the gyro component and the accelerometer component is the Euler angle concerning the i-axis.
The installation error matrix of the accelerometer can be defined by 6 small angles, as shown in Figure 3, and in the IMU coordinate frame: The unit vectors measured by the accelerometer components in the IMU are ox a , oy a , and oz a , respectively, and the order of magnitude of the cubic term of the instrument is very small, so the actual accelerometer inertial group input and output model is simplified as: Neglecting the accelerometer noise term, we can obtain: where f m is the specific force measured by the accelerometer.

43-Dimensional Kalman Filtering Model for Self-Calibration
In the process of self-calibration, there is a lever arm at the observation point of the velocity position and the rotation center of the IMU, which is called the outer lever arm in this paper. If the rotation center of the turntable does not coincide with the IMU, the observation of the velocity and position will have errors along with the rotation of the indexing mechanism. Therefore, it is necessary to analyze the effect of the outer lever arm.
Assuming that the outer lever arms between the IMU and the turntable are δl b x , δl b y and δl b z , the velocity and position observations of the IMU can be written as: where R M and R N are the earth radius parameter, h is the height. It should be noted that in Equation (15), The non-orthogonal angles of the gyro and accelerometer are usually only a few tens of arc minutes, so the formula can be written as: In this paper, a 43-dimensional Kalman filter is designed to estimate the error and calibration parameters of the IMU. The state quantities in the error equation include the scale error, zero bias, inner lever arm, outer lever arm, and gyro acceleration of the gyroscope and accelerometer. The state quantity of the filtering method proposed in this paper can be written as: where δp = δL δλ δh T , X g is the gyroscope calibration error parameters. X a represents the accelerometer calibration error parameters. δl bT is the outer lever arm vector, and δr bT is the inner lever arm vector. The gyro error vector and accelerometer error vector are shown in Equation (19).
The state function of KF can be described as: According to the previous error analysis of the IMU, the state transition matrix F can be obtained as: The elements of the F are shown as follows: Here, N g i is the output of the three gyroscopes, and N a i denotes the output of the three PIGAs. According to the analysis of the inner lever arm and the analysis of the time asynchronous error of the gyro accelerometer, the element matrices F 27 and F 26 in the state transition matrix can be obtained: The outer lever arm effect is mainly reflected in the observation equation in the Kalman filter. During the rotation of the indexing mechanism, the attitude error cannot be obtained in real time. However, during the rotation of the indexing mechanism, after compensating the outer lever arm, the observed speed and position are both 0.
Therefore, the observation transition matrix H can be written as: The elements of the H are shown as follows: Based on the previous analysis, the rotation path of the self-calibration utilizing angular velocity observation can be designed in Table 1: Table 1. Rotation path of self-calibration process.  The entire rotation path includes 18 rotation stages with a duration of 0.5 h; in the rotation stages 1-18, each rotation stage has a duration of 90 s, including the rotation movement and parking with an angular velocity of 20, where the purpose of rotation stages 1-10 is to excite and decouple all system-level self-calibration parameters; rotation stages 11-19 are mainly designed for rotation path design principle, whose purpose is to make all self-calibration parameters are fully estimated, especially the gyro bias error term. In addition, the initial alignment of the static base is performed for 120 s of coarse alignment and 180 s of fine alignment before self-calibration.

Rotation Angle along I/O Axis
The self-calibration process is designed as follows: As shown in Figure 4, the designed IMU consists of fiber optic gyroscopes (FOGs) and PIGAs, and the angular velocity observation is based on the FOG's angular velocity output. Utilizing the error models of FOG and PIGA, the state equation and measurement equation of KF can be derived, substitute the two equations into KF's time update and measurement update, the calibration results of the PIGA-based IMU can be obtained.

Experimental Results and Analysis
For the purpose of verifying the feasibility and effectiveness of the proposed online self-calibration method utilizing angular velocity observation for ultra high accuracy PIGAbased IMU. A self-calibration test is conducted to evaluate the accuracy of the calibration parameters. The circuit design of the embedded calculation and collect module is shown in Figure 5.
The flow chart of the inertial navigation signal is shown in Figure 5. The FOG transmits information such as uncalibrated angular increment, mechanical frequency jitter, and amplitude jitter to the Field Programmable Gate Array (FPGA) through the 3.3 V TTL level, and the current signal of the accelerometer passes through the I/F. After the module is converted into a frequency signal, it is transmitted to the FPGA. It should be noted that the FOG has a delay of 4 ms due to the low-pass filtering process. Therefore, the clock phase of the FPGA sampling signal to the DSP motion control module needs to be shifted to the left by 4 ms phase relative to the sampling signal to the FOG. In this way, the encoder angle information transmitted by the DSP motion control module to the FPGA through the serial port is synchronized with the output information of the IMU in time, and no related errors will be caused during the attitude demodulation process. After processing the relevant information, the FPGA transmits it to the DSP navigation module through EMIF, and also gives a 200 Hz square wave signal to the GPIO port of the DSP as the solution cycle (timed interrupt). Under the condition of large airborne dynamics, the operation rate of 200 Hz cannot meet the accuracy requirements. Therefore, in the process of the 4 k sampling of the IMU, the signal will not be accumulated, but will be latched and sent to the FPGA through EMIF. During the process of calibration, the output of PIGAs and FOGs are shown in Figure 6. It can be seen from Figure 6 that the PIGA's output is related to the rotation process, which verifies the error model derived in Section II. The accuracy of FOG we utilize in this designed IMU is 0.002 • /h (10 s, 1σ), with a 10 ppm (1σ) of scale factor repeatability. The designed RINS is fixed in the marble, and the algorithm is implemented on a digital signal processor (DSP) chip. We use the method in ref. [16] as a comparison, and the self-calibration process lasts 30 min. In addition, we use a high-precision threeaxis turntable to calibrate the IMU parameters as a reference. This method requires a high-precision turntable, and the IMU needs to be removed from the dual-axis RINS. For example, the accuracy of the dual-axis turntable is not high (especially horizontal accuracy). Traditional methods are described in [13]. The estimated curves of the IMU parameters are shown in Figures 7-9:   Particularly, we draw the angular velocity sensitivity curve of PIGA separately in Figure 10 to perform a separate analysis of its convergence. It can be seen from Figure 7 to Figure 9 that the errors added to the model do not affect the convergence of the IMU bias, scale factors, and installation angles. In Figure 10, The PIGA's angular velocity coupling factors start to converge when the turntable rotates. To better discuss the experimental effect of this method, the estimated parameters are summarized in Table 2: As shown in Table 2, the estimation accuracy of the proposed method is better than the traditional method, especially the calibration parameters of gyros. In addition, we find that the estimation results of PIGA's angular velocity coupling factors are very close to the results utilizing high accuracy offline calibration method. The errors of gyro biases estimated by the traditional method are 0.012 • /h to 0.023 • /h, using the proposed method, the errors are only within 0.003 • /h. The errors of the FOGs' scale factors estimated by the traditional method are more than 15 ppm.
The self-calibration experiment results show that the propsoed method can not only estimate the PIGA's angular velocity coupling factors, but also improve the gyroscope calibration parameters when utilzing the PIGA-based IMU.

Conclusions
Here, we propose an online self-calibration method utilizing angular velocity observation. Experimental results indicate that the proposed method accurately estimates the PIGA's angular velocity coupling factors and improves the calibration accuracy by up to 0.02 m/s/pulse simultaneously, compared with the traditional self-calibration method (with an accuracy of 0.2 m/s/pulse) for PIGA-based IMU. After compensating for the PIGA's angular velocity coupling factors, the navigation dynamic accuracy can be greatly improved. The self-calibration method also simplifies the calibration process and calibration implementation conditions, which make it possible to perform online-calibration without disassembling it and returning it to the factory for calibration.
There are still some error mechanisms that are understudied. In the future, research on decoupling the calibration of angular velocity and acceleration coupling coefficients will be carried out to improve the accuracy of PIGA continuously.