An Unconventional Multiple Low-Cost IMU and GPS-Integrated Kinematic Positioning and Navigation Method Based on Singer Model

To release the strong dependence of the conventional inertial navigation mechanization on the a priori low-cost inertial measurement unit (IMU) error model, this research applies an unconventional multi-sensor integration strategy to integrate multiple low-cost IMUs and a global positioning system (GPS) for mass-market automotive applications. The unconventional integration strategy utilizes a basic three-dimensional (3D) kinematic trajectory model as the system model to directly estimate navigational parameters, and it allows the measurements from all of the sensors independently participating in measurement updates. However, the less complex kinematic model cannot realize smooth transitions between different motion statuses for the road vehicle with acceleration maneuvers. In this manuscript, we establish a more practical 3D kinematic trajectory model based on a “current” statistical Singer acceleration model to realize smooth transitions for the maneuvering vehicle. In addition, taking advantage of the unconventional strategy, we individually model the systematic errors of each IMU and the measurements of all sensors, in contrast to most existing approaches that adopt the common-mode errors for different sensors of the same design. A real dataset involving a GPS and multiple IMUs is processed to validate the success of the proposed algorithm model under the unconventional integration strategy.


Introduction
With the popularity of intelligent vehicles, positioning and navigation technology is particularly important, because the foremost requirement for all of these safety and assistance systems is an accurate knowledge of the vehicular states (including vehicle position, vehicle velocity, vehicle acceleration, vehicle attitude, etc.) at all times [1]. High-precision sensors can certainly obtain very accurate vehicular states, but their high prices prevent them from being used in general-priced vehicles. Developing affordable sensors is the main research trend for driving applications [2]. On the other hand, the sensor should offer continuous and higher-update-rate observations under all circumstances for driving applications. Thus, considering cost-efficiency and high sampling frequency, a low-cost inertial measurement unit (IMU) is a good option and could autonomously provide navigation information. Furthermore, IMUs mostly contain a gyroscope and accelerometer for all three axes which can offer sufficient navigation information about a sudden acceleration or angular and heading change [3]. However, the accuracy of sensors decreases with the price reduction, and the performance is seriously affected by accumulated bias, drift, and sensor noises [4]. While an automotive-grade gyroscope typically gives drift performance of 1 • /h, a microelectromechanical system (MEMS) gyroscope has a typical performance of 70 • /h [1]. Plainly, low precision is one of the most critical obstacles in the development of a low-cost IMU, which limits its applications such as navigation and guidance [5]. In Figure 1, the time-dependent position vector r  of a point P is uniquely expressed using three unit vectors , , i j k    along three axes , , x y z in frame S .

r xi yj zk
In components, the motion of the moving point at a time instant t with respect to a start time 0 t can be given as shown below.

t z t z t t t z t t t z t t t
where , , , x x x    , , , y y y    , , ,... z z z    are the first, second, and third derivatives of , , x y z with respect to t . It is noteworthy that the time interval In components, the motion of the moving point at a time instant t with respect to a start time t 0 can be given as shown below.
where . x, .. z, ... z , . . . are the first, second, and third derivatives of x, y, z with respect to t. It is noteworthy that the time interval (t − t 0 ) must be short enough in order to require as few terms as where D(t) is the instantaneous rotation matrix from S to S b . According to Equation (5), we can successively derive the velocity vector, acceleration vector, and jerk vector. Moreover, D(t) is coupled with the angular motion. With respect to a rigid body in motion, the basic trajectory parameters are considered in the local navigation frame (S n ), i.e., the position r n nb of the IMU center, the velocity v n nb , the acceleration a n nb , and the jerk j n nb . In essence, v n nb is the derivative of r n nb , and it is also transformed from its opposite number v b nb in the body frame (S b ). Likewise, a n nb and j n nb are the derivatives of v n nb and a n nb , respectively. According to the vector dynamics, the trajectory parameters are directly derived as shown below.
. v n nb = a n nb . a n nb = j n where r n nb , v n nb , a n nb , j n nb , ω n nb are the position, velocity, acceleration, jerk, and angular rate vectors in S n , respectively. v b nb , a b nb , j b nb are the velocity, acceleration, and jerk vectors in S b . C n b is the direction cosine matrix (DCM) from S b to S n . The derivatives of the body velocity components and the body acceleration components are given below.
where j b nbx , j b nby , j b nbz are the body jerk components, ω b nbx , ω b nby , ω b nbz are the body angular rate components, v b nbx , v b nby , v b nbz are the body velocity components, and a b nbx , a b nby , a b nbz are the body acceleration components.

"Current" Statistical Singer Acceleration Model
The Singer model assumes that the maneuvering acceleration a(t) obeys a stationary first-order time-correlated process with zero-mean. The correlation function is expressed in the form of exponential decay as where σ 2 a is the variance of the target acceleration, and α is the reciprocal of the maneuver (acceleration) time constant, which is named "maneuvering frequency" and taken generally as an empirical value. Assuming that the probability density function of the acceleration approximately obeys uniform distribution, σ 2 a = a 2 max 3 (1 + 4P max − P 0 ) with the maximum maneuvering acceleration a max and its probability P max , and the non-maneuvering probability P 0 .
Utilizing the correlation function R a (τ), the acceleration a(t) may be expressed in terms of white noise by the Wiener-Kolmogorov whitening procedure [16], as shown below.
. where w a (t) is the white Gaussian noise with the mean of 0 and the variance of 2ασ 2 a . The "current" statistical (CS) model adopts a non-zero mean and modified Rayleigh distribution to characterize the maneuvering acceleration. Specifically, the modified Rayleigh distribution is used to describe the "current" probability density of maneuvering acceleration, and the mean value is the "current" acceleration prediction. The random acceleration still conforms to the first-order time-correlated process. According to the Singer model and the CS model, the body acceleration can be given as where a(t) is the "current" mean of maneuvering acceleration, and it is constant during each sampling period. The differential equation of body acceleration is acquired by consolidating Equation (11), as shown below.
Then, the full-state differential equation in a continuous-time system is as follows:

Attitude Angle Model
Typically, Euler angles (i.e., pitch, roll, and heading) of S b with respect to S n are selected to demonstrate the rotating properties of a 3D object. The differential equations of Euler angles are given in matrix form as follows: wherein P, γ, ψ are the pitch, roll, and heading, respectively, ω b nb is the angular rate in S b , and ω b nbx , ω b nby , ω b nbz are the components of ω b nb in three axes.

Angular Rate Model
For a normally running vehicle with a smooth steering within a small time interval, three components of ω b nb are reasonably treated as independent. In a general way, ω b nbx and ω b nby are both modeled as zero-mean processes by the first-order Markov model, and ω b nbz is modeled as a non-zero-mean random process with random disturbance. More concretely, the zero-mean Singer motion model is adopted to describe the dynamic variation of ω b nbx and ω b nby in the system model, and the modified Singer model is used to express the dynamics of ω b nbz [22,27,28].

The Formulation of the Unconventional Kalman Filter
In a general way, the state and measurement vector ought to be arranged for the system Kalman filter. In this research, the whole value state which describes the kinetic characteristic of the vehicle is one part of the state vector. The systematic errors of each IMU, e.g., the biases, scale factor errors, and so forth, are also modeled individually in the Kalman filter, instead of assuming "the common-mode errors of different sensors of the same design". Therefore, the systematic error is the other part. As for the measurement vector, the raw observables from the GPS and multiple IMUs are considered. Figure 2 illustrates the unconventional integration mechanism in this research. To demonstrate the benefits of the individual modeling for systematic errors and measurements, choosing different low-cost IMUs from different vendors is the best choice. Considering the cost of the experiment and the amount of calculation, it is advisable to integrate three low-cost IMUs and one GPS receiver in the kinematic positioning and navigation system.
In a general way, the state and measurement vector ought to be arranged for the system Kalman filter. In this research, the whole value state which describes the kinetic characteristic of the vehicle is one part of the state vector. The systematic errors of each IMU, e.g., the biases, scale factor errors, and so forth, are also modeled individually in the Kalman filter, instead of assuming "the commonmode errors of different sensors of the same design". Therefore, the systematic error is the other part. As for the measurement vector, the raw observables from the GPS and multiple IMUs are considered. Figure 2 illustrates the unconventional integration mechanism in this research. To demonstrate the benefits of the individual modeling for systematic errors and measurements, choosing different lowcost IMUs from different vendors is the best choice. Considering the cost of the experiment and the amount of calculation, it is advisable to integrate three low-cost IMUs and one GPS receiver in the kinematic positioning and navigation system.

The State Vector Including the Systematic Errors in Kalman Filter
In this manuscript, the state vector of the constructed multi-sensor integration Kalman filter contains 51 components as follows: where the position vector is expressed by triaxial coordinates in the earth-fixed coordinate system , the accelerometer bias vector

The State Vector Including the Systematic Errors in Kalman Filter
In this manuscript, the state vector of the constructed multi-sensor integration Kalman filter contains 51 components as follows: where the position vector is expressed by triaxial coordinates in the earth-fixed coordinate system , the attitude θ = ( P γ ψ ) T , the body angular rate T , the gyroscope scale factor error s g1 = ( s g1x s g1y s g1z ) T , s g2 = ( s g2x s g2y s g2z ) T , s g3 = ( s g3x s g3y s g3z ) T , and the accelerometer scale factor error s a1 = ( s a1x s a1y s a1z ) T , s a2 = ( s a2x s a2y s a2z ) T , s a3 = ( s a3x s a3y s a3z ) T . Herein, subscripts 1, 2, and 3 are the sequence numbers of the IMUs.

The Discretization of the System Model
On the basis of the differential model in Section 2, the discrete system model for KF can be summarized after omitting the higher-order terms in Taylor series expansion. It is mentionable that the position vector r of a vehicle should be in the earth-fixed coordinate frame (S e ); thus, the position vector r at epoch k + 1 in S e is calculated concretely as follows: • Firstly, the local coordinate increment ∆r n from epoch k to k + 1 in S n is as follows: • Next, ∆r n is transformed from S n to S e as follows: ∆r e = C e n(k) ∆r n . • Finally, we obtain the position vector r at epoch k + 1 as follows: Furthermore, there is a particular discussion for the discretization of the body acceleration model. The discrete state of Equation (13) can be derived as with The discrete acceleration is obtained from Equation (18) as follows: where a is the "current" mean of maneuvering acceleration, which is acquired by calculating the mean value of both sides of Equation (19), as shown below.
where z k is the sequence of all measurements up to the current time, andâ b nb(k) is the acceleration estimate of the previous moment. Equation (20) shows that a(k + 1) is not only related to the current informationâ b nb(k) , but also to the past information a(k). Since the CS model assumes that the maneuvering acceleration at an arbitrary time obeys the modified Rayleigh distribution, the variance of body acceleration can be obtained as follows: As can be seen from Equations (13) and (18), the position vector r and the body velocity vector v b nb are affected by the CS Singer acceleration model. Consequently, combining the 3D kinematic trajectory model in Section 2.1 and the CS Singer acceleration model in Section 2.2, the discrete system equations of kinematic trajectory can be derived as follows: The discrete system equations of the attitude and body angular rate are directly given as follows: The discrete system equations of the systematic errors of three different low-cost IMUs are individually modeled as follows: where i represents the sequence number of the IMU with i = 1, 2, 3, ∆t = t k+1 − t k is the time interval, T x , T y , T z are the time correlation coefficients of the first-order Markov model, w ωx , w ωy , w ωz are the independent white noises for triaxial angular rates, And w bg , w sg , w ba , and w sa are the white-noise vectors for the biases and scale factor errors of gyroscopes and accelerometers, while subscripts 1, 2, and 3 are sequence numbers of three different IMUs. j b nb is treated as process noise for the position vector, velocity vector, and acceleration vector, . ω b nb is the angular acceleration as process noise for the velocity vector, acceleration vector, and attitude vector, and C n b is the DCM. C 3×3 is the coefficient matrix as in Equation (14). C e n , the position cosine matrix, can be obtained through the position information as follows: wherein ϕ and λ are the latitude and longitude, respectively, calculated by the coordinate components of the position vector r in S e .

The Measurement Model of IMU
Generally, the IMU raw outputs include specific forces from three orthogonal accelerometers and angular rates from three orthogonal gyroscopes. We should derive three groups of measurement equations because there are three IMUs in the integration system. Figure 3 shows the structure diagram of IMUs, where Figure 3a is the general view of the vehicle, and Figure 3b represents the partial enlarged details. C is the DCM. 3 3 C × is the coefficient matrix as in Equation (14).
wherein ϕ and λ are the latitude and longitude, respectively, calculated by the coordinate components of the position vector r in e S .

The Measurement Model of IMU
Generally, the IMU raw outputs include specific forces from three orthogonal accelerometers and angular rates from three orthogonal gyroscopes. We should derive three groups of measurement equations because there are three IMUs in the integration system. Figure 3 shows the structure diagram of IMUs, where Figure 3a is the general view of the vehicle, and Figure 3b represents the partial enlarged details.  The measurement equations for low-cost IMUs could be simplified depending on specific needs [15,19]. Considering that IMUs cannot be located at the same point on the body, the measurements from different IMUs must be transformed to the same reference frame so as to perform the fusion algorithm [29]. Here, the central IMU is selected as the reference. Three groups of measurement equations for angular rate and specific force are respectively derived as shown below based on the particular structure. The measurement equations for low-cost IMUs could be simplified depending on specific needs [15,19]. Considering that IMUs cannot be located at the same point on the body, the measurements from different IMUs must be transformed to the same reference frame so as to perform the fusion algorithm [29]. Here, the central IMU is selected as the reference. Three groups of measurement equations for angular rate and specific force are respectively derived as shown below based on the particular structure.
where g n is the local gravity vector in S n , ω b nb , a b nb are the rotation rate and acceleration vectors of S b with respect to S n , S g , S a are 3 × 3 scale factor error matrices for gyros and accelerometers, b g , b a have the same meanings as mentioned in Section 3.1, and ∆ g, ∆ a are Gaussian white noises for the angular rate vector and specific force vector. The lever arm parameters of the remaining IMUs relative to the central one are r 2 = [−0.5, 0, 0] T and r 3 = [0.5, 0, 0] T .

The Measurement Model of the GPS
In this study, not only do the raw observables of IMUs participate in measurement updates, but so do those of the GPS. The GPS, as another sensor distinct from the IMU but with equal status in the system, offers two types of observables (carrier phase and pseudo-range). However, only the pseudo-range is adopted to complete the specific navigation task. The pseudo-range observation equation is generally given as where j = 1, 2, · · · , n denotes the j − th satellite, and ρ

The Implementation of the Kalman Filter
On the basis of the system model proposed above, the KF is structured straightforwardly. The state Equations (22)-(32) and measurement Equations (33)-(39) are separately generalized by the discrete nonlinear system and measurement models as shown below.
where X k is the state vector as determined in Section 3.1, Z k is the measurement vector as described in Sections 3.3 and 3.4, f() and h() are nonlinear mathematical functions, and they are constructed from Equations (22)-(32) and Equations (33)-(39), respectively, B k and Γ k are coefficient matrices, u k is the system input, u k = a(k), ∆ k is the measurement noise vector, and W k is the process noise vector including the jerk vector, the derivative of the angular rate vector, and so forth. Specifically, B k , Γ k , and W k are given as As is well known, the Kalman filter performs the prediction of the state vector through two information update processes: time update and measurement update. Specifically, the one-step prediction and the variance propagation of the state vector proceed from epoch k to k + 1 during the time update.X where P k is the mean square error matrix of the state vector, and Q k is the a priori variance-covariance matrix of the process noise vector. Φ k+1,k is the Jacobian matrix of the nonlinear system model in Equation (40) (i, j represent the serial numbers; f [i] represents the i-th system equation; X [j] represents the j-th component of state vector X), as follows: Subsequently, the measurement update is performed as follows when the measurements at epoch k + 1 are available, thus accomplishing the state estimation: where K k is the Kalman filter gain matrix, And R k is the a priori variance-covariance matrix of the measurement noise vector. H k is the Jacobian matrix of the nonlinear measurement model in Equation (41) .

Road Test and Results
The proposed algorithm model under the unconventional integration strategy was adopted to process the navigation data from multiple low-cost IMUs and a GPS-integrated system on a land vehicle. Several road tests were performed by our ground-based vehicle navigation system with a Harxon Mini Survey Antenna GPS500 and three IMUs (FOS/05M, ADIS16405BMLZ, and Crossbow IMU440CA), as shown in Figure 4. Table 1 shows the performance parameters of the three IMUs. The experimental results are given in this section from one of our road tests. The used dataset was collected in Harbin, China, from which an 8-min data fragment of the data source was selected for demonstration purposes. Figure 5a reveals the environment where the measurements were made. As can be seen from Figure 5a, the chosen test environment was an urban highway, with the Songhua River and buildings nearby. Figure 5b,c depict the trajectory and velocity profile of the test vehicle.

Further Insights into the Unconventional Integration Mechanism
The distinctions between the conventional integration mechanism (left side) and the unconventional counterpart (right side) are shown in Figure 6: (1) the embedment of the prediction of the kinematic states, (2) the introduction of IMU measurements, and (3) the navigation parameters directly used in the state vector.
The acceleration vector and angular rate vector are only regarded as inputs of the mechanization in the conventional error state-based Kalman filter, while the realistic measurements of the IMU directly participate in the measurement update process in the unconventional Kalman filter. In principle, the Kalman filter is equivalent to a sequential least square with a time-variant state vector and the process noises [30]. That is to say, the system model is In order to obtain the ground truth as the reference, we equipped a high-grade Fiber-Optic Gyroscope (FOG) inertial navigation system (INS) (gyroscope: constant drift less than 0.01 • /h, random noise less than 0.001 • /h; accelerometer: constant bias less than 100 µg, random noise less than 10 µg) developed by our research group on the test vehicle, as shown in Figure 4a. Hence, the benchmarks for attitude, velocity, and position could be provided by fusing the measurements of the high-grade INS and the GPS.

Further Insights into the Unconventional Integration Mechanism
The distinctions between the conventional integration mechanism (left side) and the unconventional counterpart (right side) are shown in Figure 6: (1) the embedment of the prediction of the kinematic states, (2) the introduction of IMU measurements, and (3) the navigation parameters directly used in the state vector.
The acceleration vector and angular rate vector are only regarded as inputs of the mechanization in the conventional error state-based Kalman filter, while the realistic measurements of the IMU directly participate in the measurement update process in the unconventional Kalman filter. In principle, the Kalman filter is equivalent to a sequential least square with a time-variant state vector and the process noises [30]. That is to say, the system model is composed of a group of virtual measurements for the state vector and reflects the connections between the state vectors from epoch to epoch. These extra virtual measurements mean that the measurement redundancies in the unconventional Kalman filter for the body acceleration vector a b nb and the body angular rate vector ω b nb are evidently better than those in the conventional Kalman filter. Therefore, the accuracy of a b nb and ω b nb in the unconventional KF for multiple low-cost IMUs and a GPS-integrated system will be undoubtedly improved because they are not only predicted but also measured. The prediction for a b nb and ω b nb in the system model can be used as a rigorous reference to check on the performance of the IMU without increasing the complexity of the filtering structure. Furthermore, the system equations of the unconventional KF can act as the dynamic constraints for the navigation parameters, for example, assuming that v b nbz = 0 and/or v b nbx = 0. Generally speaking, the most significant feature of the unconventional integration mechanism lies in the improvement of the overall measurement redundancy of the system through the system model based on the 3D kinematic trajectory model. Generally speaking, the most significant feature of the unconventional integration mechanism lies in the improvement of the overall measurement redundancy of the system through the system model based on the 3D kinematic trajectory model. As discussed above, it is expected that the accuracy of the navigation solutions shall be improved by using the rigorous trajectory model as the system model in the unconventional KF. As the components of the state vector, the acceleration and angular rate vectors in the body frame ( b S ) will also profit from the unconventional KF time updates for the same reason. For better visual effects, Figures 7 and 8 show the comparisons between the raw IMU outputs and filtered signals.
As can be seen from Figure 7, the overall variation of the IMU raw angular rates is slightly larger than that of the filtered angular rates. However, it is apparent from Figure 8 that the overall fluctuation of the filtered accelerations is significantly lower than that of the IMU raw specific forces. The improvement is due to the introduction of the novel system model in the unconventional Kalman filter. It is worth mentioning that the raw output of the accelerometer in the vertical direction includes acceleration due to gravity, while the corresponding filtered acceleration does not. That is why there is a large difference (about 10 2 m/s ) between the raw and filtered acceleration values along the vertical axis in Figure 8; here, we only compare the fluctuating ranges of the two curves. As discussed above, it is expected that the accuracy of the navigation solutions shall be improved by using the rigorous trajectory model as the system model in the unconventional KF. As the components of the state vector, the acceleration and angular rate vectors in the body frame (S b ) will also profit from the unconventional KF time updates for the same reason. For better visual effects, Figures 7 and 8 show the comparisons between the raw IMU outputs and filtered signals.
As can be seen from Figure 7, the overall variation of the IMU raw angular rates is slightly larger than that of the filtered angular rates. However, it is apparent from Figure 8 that the overall fluctuation of the filtered accelerations is significantly lower than that of the IMU raw specific forces. The improvement is due to the introduction of the novel system model in the unconventional Kalman filter. It is worth mentioning that the raw output of the accelerometer in the vertical direction includes acceleration due to gravity, while the corresponding filtered acceleration does not. That is why there is a large difference (about 10 m/s 2 ) between the raw and filtered acceleration values along the vertical axis in Figure 8; here, we only compare the fluctuating ranges of the two curves.
The comparison results from Figures 7 and 8 confirm the fact that the influences of the IMU measurement noises on the final navigation solutions are effectively mitigated due to the participation of the IMU outputs in the KF measurement updates.               Figure 9 shows the solution accuracy comparisons for 3D position using the basic 3D kinematic trajectory model or the practical model based on the "current" statistical Singer acceleration model (CSSAM). The overall 3D position accuracy using the proposed algorithm model is under 5 m , which is a great improvement compared with the accuracy (10 m ) using the basic one. Figure 10 shows that the estimation errors for the velocity state vector in the three axial directions of the body frame using the proposed algorithm model are within 0.05 0. 15 − ∼ / m s , 0.2 ± / m s , and 0.2 ± / m s , respectively. The velocity error in the direction of travel is obviously reduced compared with that using the basic model. Figure 11 shows that the estimation errors for acceleration state vector in the three axial directions of the body frame are within 0.5 ± 2 / m s , 1 ± 2 / m s , and 2 ± 2 / m s , respectively. As can be seen, the accuracy of the acceleration estimation is apparently improved by using the proposed algorithm model. As shown in Figure 12, the accuracies for attitude (pitch, roll,   Figure 9 shows the solution accuracy comparisons for 3D position using the basic 3D kinematic trajectory model or the practical model based on the "current" statistical Singer acceleration model (CSSAM). The overall 3D position accuracy using the proposed algorithm model is under 5 m , which is a great improvement compared with the accuracy (10 m ) using the basic one. Figure 10 shows that the estimation errors for the velocity state vector in the three axial directions of the body frame using the proposed algorithm model are within 0.05 0.15 − ∼ / m s , 0.2 ± / m s , and 0.2 ± / m s , respectively. The velocity error in the direction of travel is obviously reduced compared with that using the basic model. Figure 11 shows that the estimation errors for acceleration state vector in the three axial directions of the body frame are within 0.5 ± 2 / m s , 1 ± 2 / m s , and 2 ± 2 / m s , respectively. As can be seen, the accuracy of the acceleration estimation is apparently improved by using the proposed algorithm model. As shown in Figure 12, the accuracies for attitude (pitch, roll,  Figure 9 shows the solution accuracy comparisons for 3D position using the basic 3D kinematic trajectory model or the practical model based on the "current" statistical Singer acceleration model (CSSAM). The overall 3D position accuracy using the proposed algorithm model is under 5 m, which is a great improvement compared with the accuracy (10 m) using the basic one. Figure 10 shows that the estimation errors for the velocity state vector in the three axial directions of the body frame using the proposed algorithm model are within −0.05 ∼ 0.15 m/s, ±0.2 m/s, and ±0.2 m/s, respectively. The velocity error in the direction of travel is obviously reduced compared with that using the basic model. Figure 11 shows that the estimation errors for acceleration state vector in the three axial directions of the body frame are within ±0.5 m/s 2 , ±1 m/s 2 , and ±2 m/s 2 , respectively. As can be seen, the accuracy of the acceleration estimation is apparently improved by using the proposed algorithm model. As shown in Figure 12, the accuracies for attitude (pitch, roll, and heading) are within ±2 • , −0.2 ∼ 0.6 • , and −1 ∼ 3 • , respectively. It is to be observed that the pitch angle has an apparent fluctuation from 3-4.5 min. The trajectory curve in Figure 5b illustrates that the vehicle experienced ups and downs in the test duration, which could affect the estimation of attitude angles (especially for pitch angle). It also indicates that the attitude estimation model still needs improvement so that it can be applied to more complex motion forms.

Verification of the Proposed Algorithm Model under the Unconventional Integration Strategy
Apparently, the navigation parameters during the 8-min vehicle motion were estimated within acceptable ranges by adopting the proposed algorithm model under the unconventional integration strategy. Furthermore, the reason why the 3D position error using the basic trajectory model did not diverge is that the output of the central accelerometer with real-time calibration was utilized to regulate the acceleration estimation error when significant acceleration maneuvers occurred.

IMU Systematic Error Estimation
Since the systematic errors of these IMU arrays, i.e., biases and scale factor errors of gyroscopes and accelerometers, are individually modeled in Kalman filtering, the systematic error estimations of each IMU can be obtained separately, as illustrated in Figures 13-16. As an example, this manuscript only shows the systematic errors of the central IMU due to space limitations. The common existing approach under the conventional inertial navigation mechanization adopts a set of common shared errors for all IMUs and, sometimes, these error parameters are from initial calibration results or technical specifications [7,28]. However, in fact, the a priori error model defined for a static low-cost MEMS inertial sensor needs to be checked and compensated for in the dynamic working environment as the vibration on a low-cost IMU might cause significant changes in its scale factors and noise level compared to those in the static case [31]. and heading) are within 2 ±°, 0.2 0.6°−  , and 1 3 −° , respectively. It is to be observed that the pitch angle has an apparent fluctuation from 3-4.5 min. The trajectory curve in Figure 5b illustrates that the vehicle experienced ups and downs in the test duration, which could affect the estimation of attitude angles (especially for pitch angle). It also indicates that the attitude estimation model still needs improvement so that it can be applied to more complex motion forms.
Apparently, the navigation parameters during the 8-min vehicle motion were estimated within acceptable ranges by adopting the proposed algorithm model under the unconventional integration strategy. Furthermore, the reason why the 3D position error using the basic trajectory model did not diverge is that the output of the central accelerometer with real-time calibration was utilized to regulate the acceleration estimation error when significant acceleration maneuvers occurred.

IMU Systematic Error Estimation
Since the systematic errors of these IMU arrays, i.e., biases and scale factor errors of gyroscopes and accelerometers, are individually modeled in Kalman filtering, the systematic error estimations of each IMU can be obtained separately, as illustrated in Figures 13-16. As an example, this manuscript only shows the systematic errors of the central IMU due to space limitations. The common existing approach under the conventional inertial navigation mechanization adopts a set of common shared errors for all IMUs and, sometimes, these error parameters are from initial calibration results or technical specifications [7,28]. However, in fact, the a priori error model defined for a static low-cost MEMS inertial sensor needs to be checked and compensated for in the dynamic working environment as the vibration on a low-cost IMU might cause significant changes in its scale factors and noise level compared to those in the static case [31].       As can be seen from Figures 13-16, in contrast to the a priori constant systematic error with the common existing approach, the systematic error estimations with the individual modeling method vary with time, which conforms better with the real situation. In other words, the individual modeling method under the unconventional integration strategy can adjust the systematic error estimations of each IMU according to the real-time measurements from each IMU. This technique can firstly verify if those IMUs really share the same systematic errors quantitatively, even if they are the common errors physically; then, they can be modeled either separately or combined, laying the groundwork for future research such as auto calibration and fault detection. While applying the individual modeling method for the IMU array, in spite of its characteristic properties, one must confront the following problem: high computation load caused by high-rate measurement updates in the Kalman filter, e.g., with IMUs whose measurement rate may be at 100 Hz. Although the modern computation capability is considerably improved, one cannot stand such a high measurement update rate, especially one so unreasonable, which renders it useless in real time. How to appropriately reduce the high-rate KF measurement updates without compromising the valuable information embedded in the high-rate measurements will become a topic for further study.

Conclusions
This research fused the information from one GPS receiver and three low-cost IMUs by applying an unconventional multi-sensor integration strategy. The enhanced and improved parts involved establishing a more practical 3D kinematic trajectory model based on the "current" statistical Singer acceleration model as the core of the system model, and individually modeling the measurements and systematic errors of these IMU arrays in Kalman filtering. The processing results of the experimental data demonstrated the success of the proposed algorithm model under the unconventional integration strategy with satisfactory solution performance and reliability. Future work will involve developing a more precise fusion algorithm by using the carrier phase information from the GPS.  As can be seen from Figures 13-16, in contrast to the a priori constant systematic error with the common existing approach, the systematic error estimations with the individual modeling method vary with time, which conforms better with the real situation. In other words, the individual modeling method under the unconventional integration strategy can adjust the systematic error estimations of each IMU according to the real-time measurements from each IMU. This technique can firstly verify if those IMUs really share the same systematic errors quantitatively, even if they are the common errors physically; then, they can be modeled either separately or combined, laying the groundwork for future research such as auto calibration and fault detection. While applying the individual modeling method for the IMU array, in spite of its characteristic properties, one must confront the following problem: high computation load caused by high-rate measurement updates in the Kalman filter, e.g., with IMUs whose measurement rate may be at 100 Hz. Although the modern computation capability is considerably improved, one cannot stand such a high measurement update rate, especially one so unreasonable, which renders it useless in real time. How to appropriately reduce the high-rate KF measurement updates without compromising the valuable information embedded in the high-rate measurements will become a topic for further study.

Conclusions
This research fused the information from one GPS receiver and three low-cost IMUs by applying an unconventional multi-sensor integration strategy. The enhanced and improved parts involved establishing a more practical 3D kinematic trajectory model based on the "current" statistical Singer acceleration model as the core of the system model, and individually modeling the measurements and systematic errors of these IMU arrays in Kalman filtering. The processing results of the experimental data demonstrated the success of the proposed algorithm model under the unconventional integration strategy with satisfactory solution performance and reliability. Future work will involve developing a more precise fusion algorithm by using the carrier phase information from the GPS.