IMU Signal Generator Based on Dual Quaternion Interpolation for Integration Simulation.

This paper focuses on the problem of high-update-rate and high accuracy inertial measurement unit signal generation. In order to be in accordance with the vehicle's kinematic and dynamic characteristics as well as the characteristics of pseudorange of post-processed global navigation satellite system and their rate measurements, a novel dual quaternion interpolation and analytic integration algorithm based on actual flight data is proposed. The proposed method can simplify the piecewise analytical expressions of angular rates, angular increments and specific force integral increments. Norm corrections are adopted as constraint conditions to guarantee the accuracy of the signals. Numerical simulations are conducted to validate the method's performance.

1. Introduction

Inertial measurement unit (IMU) signal generators have been used in numerous engineering fields for algorithm testing and verification to reduce experimental costs and shorten developing time under laboratory conditions. It can be used in industrial robots for motion control [1,2,3], path planning of unmanned aerial vehicles (UAV) [4], tightly coupled integrated navigation research [5,6,7] and guidance system simulation [8]. A Global Navigation Satellite System (GNSS) signal simulator is required to generate multiple synchronous signals, one for each satellite in view, and to generate consistent navigation data for the receiver to be able to compute the position. Moreover, a corresponding IMU signal simulator is also necessary to provide gyro and accelerometer measurements for the strapdown inertial navigation system (SINS) and GNSS integration research.

The followings are required for a IMU signal generator to be used in the tightly coupled SINS/GNSS integration simulation:

GNSS’s pseudorange and their rate measurements corresponding to the simulated navigation parameters such as position and velocity as well as satellite ephemeris are required.

The simulated gyro measurements and accelerometer measurements in the trajectory are coincident with the vehicle’s kinematic and dynamic characteristics.

The IMU size effects as well as the level arm effects caused by the displacements of the GNSS receiver antenna from the IMU reference point can be simulated. For this purpose, the analytic angular rate and the angular acceleration data are required.

Various SINS algorithms and update rates can be evaluated under the dynamic environment in trajectories through simulations.

In this paper, an actual-flight-data based IMU signal generation algorithm is presented. Actual GNSS navigation data files (including satellite ephemeris) and observation data files (including GNSS’s pseudorange and their rate measurements) can be obtained through practical flight experiments. The post-processed GNSS navigation solution and a low-cost loosely coupled SINS/GNSS integrated navigation system provide position, velocity and attitude parameters at discrete time intervals (1 Hz). The spline interpolation method is used to generate analytic continuous-time kinematic simulation trajectory. The main contributions of this paper can be summarized as follows: (1) A novel IMU signal generation algorithm is proposed based on dual quaternions and actual flight data. (2) The effectiveness and efficiency of the proposed algorithm are validated by extensive simulation experiments. (3) Compared with state-of-the-art methods, the proposed method generates signals with the highest accuracy in a relatively short computation time.

This paper is organized as follows. Section 2 introduces the related works about the IMU signal generator and dual quaternions. Section 3 provides a dual quaternion representation of angular increments and specific force integral increments for IMU measurements simulation. Section 4 provides two-point piecewise cubic spline interpolation, analytic integration, norm corrections and error analysis of interpolated dual quaternions. Simulation tests for analytic IMU signal generator validation are presented in Section 5. Conclusions are drawn in Section 6.

2. Related Works

For the IMU signal generator, over the decades, various approaches have been developed, which can be classified into three general categories: pure mathematical models [9,10,11,12,13], kinematic or dynamic models [14,15,16,17] and actual flight-data based methods [18,19].

The methods based on pure mathematical models mainly depend on the vehicle’s straightforward manoeuvre, such as climbing up, straight flight, turning flight and diving. Moreover, the integration of these manoeuvres can be simulated. Two types of highly straightforward models are commonly used. The first one is the linear acceleration that has phase step with 50 g, wherein the acceleration is a constant. The other one is the circular motion or turn that has radial acceleration with 50 g; furthermore, all the velocities and derivations follow the sinusoidal waveform. Although the form is adopted, all of them are highly straightforward. A trajectory generated by these methods can neither reflect the complex kinematics and dynamics of a vehicle nor be conducive to combining with actual flight experiments to produce consistent simulation trajectory using satellite signals and multi-sensor measurement data such as those from image sensors. In kinematic and dynamic model based methods, classical kinematic or dynamic equations are used to generate inertial sensor measurements under flight simulation conditions, and the impact of the pneumatic environment of a vehicle can also be considered. The drawback is that it is generally challenging to predetermine the position and attitude of the vehicle precisely in the simulated trajectory. In addition, physical and mathematical models are often excessively complex to be described, and the degree of approximation between the generated trajectory and the vehicle’s actual motion depends on the accuracy of the kinematic models. Based on actual flight data, the incremental outputs of the inertial sensors are inversely deduced from the attitude, velocity and positioning information using a numerical strapdown inertial navigation system (SINS) attitude and velocity updating algorithm, which can be referred to as an inverse SINS algorithm. Generally, this algorithm requires high update rate of data of navigational parameters (100 Hz–1000 Hz). Because the inverse SINS algorithm is used, the accuracy influences of the strapdown inertial navigation algorithm and update rates are not considered. It is not feasible to simulate error factors such as size effects and lever arm effects because of the lack of the angular rate and the angular acceleration data in numerical simulations using the inversed SINS algorithm.

Dual quaternions (DQ) can provide a uniform representation of the rotational and translational motion of a rigid body in space. It was proposed by Clifford in 1873 and based on the dual quaternion algebra to describe the spiral variation. In 1992, dual quaternion was first used in the inertial navigation field by Branets and Shmyglevsky [20]. They analysed its feasibility and drew the outline of various coordinates’ relationships in SINS. In 2004, Y.X Wu elaborated upon strapdown inertial navigation system algorithms based on the dual quaternion and demonstrated the consistency of the differential equations between dual quaternion algorithms and traditional algorithms [21]. Compared with the general methods, the superior accuracy performance of the dual quaternion has been proved [21,22]. In [21], the author analysed the algorithm errors and proved that the algorithm based on DQ was superior to conventional ones in terms of the accuracy of the translational components (velocity vectors and position vectors) from both the analytic forms and the numerical simulation perspectives. Moreover, it was indicated that the new algorithm was a better choice than all conventional algorithms for applications in high-accuracy navigation systems and high-manoeuvre applications. In [22], the accuracy advantage of the algorithm based on DQ was analysed, and the analytic comparisons indicated that the thrust velocity solution as well as the attitude and gravitational velocity solution were superior to the traditional navigation algorithms in terms of accuracy.

In this study, a novel dual quaternion interpolation and analytic integration algorithm is used to generate high update rate and high accuracy IMU signals for actual-flight-data based integrated navigation simulation, in which the piecewise analytical expressions of angular increments and specific force integral increments in the body frame are simplified. A method to solve the norm constraint problems in quaternion and dual quaternion interpolation is presented.

3. Dual Quaternion Representation of Angular Increments and Specific Force Integral Increments

In this section, the relationship between dual quaternion and angular rates in the body frame and that between dual quaternion and specific force in the body frame are presented for ease of reference and for providing the foundation for the dual quaternion interpolation.

3.1. Rotation and Translation Expressed as a Dual Quaternion

Rotation and translation between the inertial frame (i frame) and the body frame (b frame) can be expressed as a dual quaternion in the following manner [21] where, is the dual quaternion, and ε2 = 0. The real part q is the general quaternion to facilitate the representation of rotation, p is the dual part of , and represents the quaternion or dual quaternion multiplication. The translation component t can be defined as (see Appendix A for details) where, is the velocity vector of frame b relative to frame i expressed in frame c, is the gravitational acceleration vector expressed in frame c.

The rotational and translational kinematics can be expressed as a differential equation of the dual quaternion as in which where, is a dual vector and generally referred to as twist, × represents the vector product, and represents the angular rate vector of frame c2 relative to frame c1.

Using and the fact that , where * represents the conjugate operator, the following is obtained: where is the specific force vector expressed in frame c, is the acceleration vector of frame b relative to frame i expressed in frame c.

Substituting Equation (5) in Equation (4), we obtain where the twist is a combination of the angular velocity and the specific force in frame b.

Substituting Equation (6) in Equation (3), the twist is obtained as follows:

3.2. Integral Increments of Angular Rates and Specific Force in Body Frame

Integrating Equation (7) over the sampling time interval , the angular rate increments and the specific force integral increments at time can be expressed as where considering and are quaternions with zero scalar parts.

4. Piecewise Constraint Interpolation and Analytic Integration Algorithm

In this section, a novel IMU signal generation algorithm is developed based on the dual quaternion interpolation, which is suitable for implementation in navigation fields. The two-point piecewise cubic spline interpolation and the analytic integration are used for the dual quaternion, and angular increments and specific force increments can be obtained. Additionally, the error analysis and methods for accuracy enhancement are also presented.

An algorithm flow chart is illustrated in Figure 1 that provides the detailed steps to streamline the algorithm proposed in this paper. The proposed method provides a necessary basis for the offline test of the performance of SINS and its fusion with other sensors’ signals to verify the performance of the integrated navigation algorithm. By using this algorithm, we can get IMU signal data with any frequency. If the sensor (IMU) outputs and all the flight data can be collected, the proposed algorithm is not necessary.

Figure 1

Algorithm flow chart.

The algorithm can be divided into the following three parts:

Attitude Part

Convert Euler angles (roll, pitch and yaw) from the navigation frame n to quaternion in frame i.

Set the interpolation interval . The cubic spline interpolation is adopted for the quaternion. Combined with the constraint conditions, the discrete quaternion and its derivative are obtained.

Using as the new interpolation interval, the polynomial of are obtained.

Calculate the angular rate and its integral increments.

Use average approximation to modify the values of angular increments.

Position Part

Convert the position information (latitude, longitude and altitude) from the navigation frame n to the inertial frame i.

Set the interpolation interval . The cubic spline interpolation is adopted for the position and velocity. Then, the discrete velocity values can be obtained.

Gravity values in the Earth-centered inertial (ECI) coordinates are calculated according to the gravity formula.

Set the interpolation interval . The cubic spline interpolation is adopted for the gravity values. Then, discrete gravity integral values can also be obtained.

Using the discrete quaternion and its derivative, discrete velocity values and its derivative, discrete gravity and its integral values, discrete and its derivative can be obtained.

Using as the new interpolation interval, the polynomials of are obtained.

Combined with the polynomials of , the specific force and its integral increments are obtained.

Use average approximation to modify the values of specific force integral increments.

Solution Part

Based on the two sub-sample error compensation algorithms, the inertial sensor simulated samples and are used to update pure SINS. The attitude and the position results are compared with the actual flight data input.

In this algorithm, two sample time intervals are involved. The first one is the interpolation time interval, which is the interval between the interpolation points selected from the original input data. The second one is the integration time interval, which is the interval for the integration of angular rate and specific force.

4.1. Dual Quaternion Calculation at Discrete Time Based on Actual Flight Data

The actual flight data including navigational parameters such as position , velocity and quaternion at discrete time can be obtained through a SINS/GNSS loosely integrated navigation system in a UAV.

The dual quaternion and its derivative at discrete time can be calculated based on the following equations: where, is the function of . The gravitation integration is calculated as follows: can be obtained by the piecewise cubic spline interpolation in the time period , .

Differentiate Equation (14) twice, and let , and . in Equation (12) can be calculated as in Equation (12) can be obtained by the interpolation.

The piecewise cubic spline interpolation for the quaternion is expressed by [23] where are the coefficients of the polynomial, and they vary across the interpolation intervals. The subscript m (or in the following content) represents the time point.

represents the piecewise cubic spline quaternion in the time interval

Considering the constraint condition , the quaternion interpolation and its derivatives can be rectified as follows: where, Let , considering where are the four elements of the quaternion, and are the derivatives of the four elements of quaternion.

4.2. Two-Point Piecewise Cubic Spline Interpolation for Dual Quaternions

In order to obtain the analytic integrations of Equations (8) and (9), the cubic spline interpolation for the dual quaternion is expressed by where, , subject to the following constraints of two-point boundary conditions is the dual quaternion interpolation in the time interval τ ∈ [0, T] without normalization, is the quaternion interpolation in the time interval τ ∈ [0, T] without normalization, and is the dual part interpolation of in the time interval τ ∈ [0, T].

4.3. Analytic Integration of Dual Quaternions

The twist, which is defined in the dual quaternion interpolations, is expressed by

By substituting Equations (24) and (25) into Equation (32), integrating Equation (32) over the sampling time interval and supposing , the following equations can be determined. represents the integration frequency.

4.4. Norm Corrections of Dual Quaternion

The norm of the dual quaternion representing rotation and translation between the inertial frame and the body frame is defined as in which, , and is the square of the dual quaternion norm. Generally, the cubic spline interpolation and do not satisfy these constraints, which results in undesirable fluctuation errors in and .

To solve these problems, norm corrections of the dual quaternion interpolation is brought into the analytic integrations as follows:

Substituting Equation (36) into Equation (7), the corrected twist is derived as where

Substituting Equation (38) into Equation (37), the angular rates and specific force can be expressed as is scalar and can be neglected in the vector computation for and .

Considering the fact that sampling time is generally marginal (5 ms or less) over the sampling time interval , the following approximations can be adopted:

Based on the above equations, the integral angular rate increments in Equation (8) and specific force integral increments in Equation (9) can be rewritten as

By this procedure, simulated IMU signals including gyro measurements and accelerometer measurements can be calculated by the dual quaternion interpolation and the analytic integration.

As and can be obtained analytically, the error sources such as the size effects and the lever arm effects can also be simulated [24]. For example, when , and (the displacements of the x-, y- and z- axis accelerometers, respectively, from the IMU reference point expressed in frame b ) are specified, the size effect caused by the displacements of the accelerometers from the IMU reference point can be simulated. Similarly, when (the displacement of the GNSS receiver antenna from the IMU reference point expressed in frame b ) is specified, the lever arm effect caused by the displacement of the GNSS receiver antenna from the IMU reference point can also be simulated.

4.5. Error Analysis and Methods for Accuracy Enhancement

The error caused by approximations in Equations (40) and (41) can be expressed as: where is defined in Equation (41), is the error vector of the angular rate, and is the error vector of the specific force.

For a low dynamic vehicle trajectory, amplitudes of and are generally negligible; therefore, errors in Equations (44) and (45) can be neglected. For a high dynamic vehicle trajectory with large angular rates and specific force, a smaller sampling integration time interval is to be selected to reduce the errors.

In terms of the computational complexity, this algorithm only needs a series of sequential calculations including interpolation, differentiation, multiplication and integration for each point. When there are n points, the number of calculation steps is proportional to n. Therefore, the computational complexity is .

5. Simulation Tests for Analytic IMU Signal Generator Performance Validation

This section describes the numerical simulations to validate the proposed algorithm. We choose the quaternion error and the position error as the test metrics to illustrate the performance of the algorithm. The quaternion error means the error between the quaternion of the original data we input and the quaternion of the results we obtained. The position error means the error between the position in ECI coordinates of the original data we input and the position in ECI coordinates of the results we obtained. The results of quaternion and the position in ECI coordinates are obtained by using the strapdown inertial navigation solution with the angular increment and the specific force increment obtained by the proposed algorithm.

The actual flight data including position, velocity, attitude and time as well as GNSS’s pseudoranges and their rate measurements have been procured from post-processed GNSS navigation solutions and a low-cost loosely coupled SINS/GNSS integrated navigation system in a UAV. The information of the attitude and the position are used for the algorithm proposed in this paper to generate the IMU signal. The update rate of flight data is 1 Hz and the total duration is 1200 s. The flight profile includes climbing, turning, cruise and a few flying manoeuvres, which are illustrated in Figure 2.

Figure 2

Actual flight trajectory.

We choose 1 s as the analytic interpolation interval, i.e., 1, which is consistent with the sampling frequency of the satellite signals. The simulated angular rate integral increments and the specific force integral increments are generated based on the dual quaternion interpolations at a frequency of 400 Hz. That is to say, 0.0025 s. It is a common sampling frequency of SINS in engineering. The strapdown inertial navigation simulation is performed to verify the accuracy of the simulated gyro and accelerometer measurements. The coning and sculling compensation algorithms are used for the attitude and velocity update, while the trapezoidal integration method is used for the position update [24]. The strapdown inertial navigation simulation results are compared with the flight trajectory at time without considering the initial alignment errors and instrument errors such as gyro drifts, accelerometer biases and scale factor errors. The attitude and position error comparisons of the simulation results are illustrated in Figure 3 and Figure 4, respectively.

Figure 3

Curve of the quaternion error (based on DQ).

Figure 4

Curve of the position error (based on DQ).

It can be observed in Figure 3 and Figure 4 that the attitude quaternion error and position error accumulate over time, which is consistent with the inertial navigation’s error characteristics (the attitude error and the position error accumulate with time). In Figure 4, the fluctuation is mainly owing to the use of the average approximation and the computing error, and it can be neglected. The quaternion error is smaller than , and the position error is smaller than 0.021 m. The simulated angular rate increments and specific force integral increments are highly satisfactory and show high accuracy because the errors in Figure 3 and Figure 4 are substantially smaller than those caused by misalignments or bias and scale factor errors of gyros and accelerometers.

In order to further illustrate the algorithm’s performance, the inverse SINS algorithm in [18] and the trajectory generation algorithm based on quaternion [19] are used with the same actual flight data. For the inverse SINS algorithm, the interpolation interval is also set as 1 s and the integrated frequency is 400 Hz. The attitude and position error comparisons of the simulation results are illustrated in Figure 5 and Figure 6, respectively.

Figure 5

Curve of the quaternion error (based on the inverse SINS).

Figure 6

Curve of the position error (based on the inverse SINS).

According to the trajectory generation algorithm based on the quaternion [19], the following can be obtained: where is the function of quaternion . Using the first-order Taylor expansion for , Equations (46) and (47) can be rewritten as

With the same actual flight data, the error result can be obtained by MATLAB simulation. The attitude and the position errors are illustrated in Figure 7 and Figure 8, respectively.

Figure 7

Curve of the quaternion error (based on the quaternion).

Figure 8

Curve of position error (based on quaternion).

In the strapdown inertial navigation system, the attitude accuracy is highly critical for the whole system’s calculation. Through the simulation, it can be observed that the dual quaternion interpolation algorithm developed in this study shows the highest attitude accuracy, with an error of . The error in the quaternion interpolation algorithm is lower than that of the inverse SINS algorithm in [18]. Then, the position errors of the three algorithms are compared. From Figure 4, Figure 6 and Figure 8, it can be observed that the dual quaternion interpolation algorithm developed in this study shows the highest accuracy also, with a position norm error of m. The quaternion interpolation algorithm is the second, with a position norm error of h m. The performance of the inverse SINS algorithm is the worst, with a position norm error of m. Furthermore, high order Taylor expansion approximation can be used in Equations (46) and (47) to enhance attitude and position accuracy and the algorithm’s performance. However, when comparing Equation (42) with Equation (48) and Equation (43) with Equation (49), it is observed that Equations (48) and (49) show higher order than Equations (42) and (43), although only the first-order Taylor expansion approximation is used. An increase in the order of the equations increases the time cost.

The computation time of the three algorithms is shown in Table 1. From the data in the table we can see that the computation time of the three algorithms is similar for the dataset containing 1200 data points. The algorithm based on the inverse SINS takes the most time and its accuracy is the worst among the three algorithms. Although the proposed algorithm takes slightly more time than the algorithm based on quaternion, the attitude accuracy obtained by the former one is better than that of the latter one. The longer time is acceptable.

Table 1

The computation time of the three algorithms.

Algorithm Proposed in This Paper	Algorithm Based on Inverse SINS	Algorithm Based on Quaternion
784.74 s	813.14 s	750.92 s

Therefore, it can also be concluded that the analytic IMU signal generation algorithm, which is based on the actual flight data of a UAV and dual quaternion interpolations, show high accuracy. Moreover, it is capable of satisfying the requirements of dynamic simulations of SINS/GNSS integrated navigation research studies, in which simulated IMU signals are to be consistent with GNSS’s pseudoranges and their rate measurements in an actual flight test.

6. Conclusions

This paper proposed a new analytical IMU signal generation algorithm based on the dual quaternion interpolation. Compared with other algorithms, the new algorithm provides a concise analytical form for angular rate increments and specific force integral increments, and has higher accuracy than other state-of-the-art algorithms. The simulation results indicate that the quaternion error is smaller than , and the position error is smaller than 0.021 m. The accuracy of the proposed algorithm is capable of satisfying the requirements of dynamic simulations of the integrated navigation.

The applications can be extended to other SINS/multi-sensor integrated navigation or rigid body motion control simulations. When the position, velocity and attitude data at discrete time epoches are given, a continuous smooth six-degrees-of-freedom (6DOF) kinematic or dynamic trajectory can always be generated through the dual quaternion interpolation.