Robust SCKF Filtering Method for MINS/GPS In-Motion Alignment

This paper presents a novel multiple strong tracking adaptive square-root cubature Kalman filter (MSTASCKF) based on the frame of the Sage–Husa filter, employing the multi-fading factor which could automatically adjust the Q value according to the rapidly changing noise in the flight process. This filter can estimate the system noise in real-time during the filtering process and adjust the system noise variance matrix Q so that the filtering accuracy is not significantly reduced with the noise. At the same time, the residual error in the filtering process is used as a measure of the filtering effect, and a multiple fading factor is introduced to adjust the posterior error variance matrix in the filtering process, so that the residual error is always orthogonal and the stability of the filtering is maintained. Finally, a vibration test is designed which simulates the random noise of the short-range guided weapon in flight through the shaking table and adds the noise to the present simulation trajectory for semi-physical simulation. The simulation results show that the proposed filter can significantly reduce the attitude estimation error caused by random vibration.


Introduction
With the development of Micro-Electro-Mechanical System (MEMS) technology, the MEMS inertial measurement unit (MIMU) is gaining attention due to its integrability and miniaturization in several fields such as guided weaponry, unmanned automatic vehicles (UAVs), and robots [1][2][3][4][5][6]. As a dead-reckoning method, the initial alignment accuracy has a vital impact to the final inertial navigation results, while there are only tens of seconds for a tactical weapon. Thus, achieving rapid, accurate, and robust initial alignment is essential to the application of the MIMU. As for the MIMU, it cannot detect the Earth's rotation properly with high gyro noise, so alignment methods for the high-accuracy Inertial Navigation System (INS) are not suitable for it. A widely used method to achieve initial alignment for the MIMU is to use an extra sensor to provide additional information as a reference for initial alignment during the flight process. The in-flight alignment has been investigated as one of the key technologies for guided weapons and UAV applications in MIMUs in the past 10 years.
The numerous methods to solve in-flight alignment for the MIMU could be mainly divided into two categories, namely the optimization-based approach and the Kalman filter approach. The first one usually turns the problem of alignment into a Whaba's problem. To solve the problem, several algorithms have been developed, such as the q-method, Quaternion ESTimator (QUEST), and Recursive QUEST (REQUEST) [7,8]. These methods calculate the optimal quaternion by constructing the optimal observation vectors through different filtering methods [9][10][11][12][13]. The latter method is to obtain the misalignment angle between the navigation frame to the computing frame of the INS by constructing an appropriate model for filtering. Owing to the nonlinearity of the misalignment model and the large noise caused by environment, the nonlinear filtering methods occupy a great proportion of the in-motion alignment methods to achieve stability, accuracy, and robustness.
In order to handle with the strong nonlinear system, several filtering methods based on the Kalman filter have been proposed, such as the extended Kalman filter (EKF) [14,15], unscented Kalman filter (UKF) [16][17][18], Gauss-Hermite quadrature filter (GHQF) [19,20], sparse grid quadrature filter (SGQF) [21][22][23], and cubature Kalman filter (CKF) [24,25]. The EKF approximates the nonlinear model to a linear model through Taylor expansion to deal with nonlinear problems. The advantage of this is that the filter structure is simple and the Kalman filter framework can be directly applied. Fang proposed an innovative adaptive extended Kalman filter method for the aerial alignment of airborne position and azimuth measurement systems [14]. This method improves the real-time performance of the algorithm and reduces the interference from GPS measurement noise. However, EKF can only achieve first-order accuracy and cannot handle strong nonlinear problems. In order to solve this problem, UKF is proposed. UKF assumes that the nonlinear system is a Gaussian probability density function (PDF) distribution and creates a series of sigma points based on the unscented transform rule to approximate the distribution of the system. It can achieve the second-order Taylor expansion polynomial. Saman Mukhtar Siddiqui applied the initialization process of the Central Difference Unscented Kalman filter (CDUKF) to the square-root unscented Kalman filter and verified the advantages of CDUKF in calculation time and accuracy in the on-board experiment to solve the alignment problem of misalignment angles within 30 • [26]. Wang Dingjie proposed an adaptive unscented Kalman filter for a small UAV MEMS navigation system. With an initial error angle of 30 • , the posture is completed in 75 s, which is more accurate than EKF and traditional UKF [27]. However, there is also a problem: For systems with a dimension higher than 3, the weight of the sigma points could be negative, and the filter process will be unstable. As a matter of fact, this problem is caused by the Unscented Transformation (UT) rule. To avoid the negativity of the sigma points, GHQF, SGQF, and CKF have been proposed using different integrated rules rather than the UT rule, which not only solves the mentioned problem but also maintains the accuracy of the second-order Taylor expansion. The GHQF utlizes the Gauss-Hermite quadrature (GHQ) rule to conduct the filter process with the quadrature sampling points that it can be accurate to the arbitrary order Taylor expansion polynomial [23]. However, the drawback of the GHQ rule is that it needs a large number of the sampling points due to the application of the direct tensor product. Thus, the GHQF would also have problems in application in high-dimensional systems. CKF with the cubature rule guarantees the positive weight, and its number of sampling points is 2n. It could also achieve a higher accuracy while applying the high-degree cubature rule with less computational burden [24]. Compared with CKF, the square-root CKF (SCKF) has more advantage not only in convergence speed but also in calculation efficiency [25]. To achieve a better performance of robustness during the in-flight process, we added several methods, such as Sage-Husa noise statistic and fading factors, to the frame of SCKF.
The rest of this paper is organized as follows. An Optimisation-Based initial Alignment (OBA) method based on the Rodrigues parameter for random misalignment angle is introduced in Section 2, where the state model and the measurement equation are used in the simulation. Then, the specific formulation of MSTASCKF is presented in Section 3. In Section 4, a semi-physical simulation is carried out to evaluate the performance of the proposed algorithm under the vibration condition. The alignment results are compared with the conventional SCKF. Finally, the conclusions are drawn in Section 5.

Alignment Model for Random Misalignment Angle
In this section, a dynamic alignment model in the form of a random misalignment angle [28] is introduced.
The coordinate systems are defined as follows: Body frame (b-frame): The origin of the coordinate system is at the center of mass of the carrier, the Y axis is the longitudinal axis of the carrier, the X axis points to the right of the carrier and is perpendicular to the Y axis, and the Z axis and X and Y axes form a right-handed coordinate system. Navigation frame (n-frame): The origin of the coordinate system is at the center of mass of the carrier, the X axis points to the geographic north direction, the Y axis points to the geographic east direction, and the Z axis points to the up direction.
Inertial navigation frame (i n -frame): An inertial coordinate system, which coincides with the initial navigation frame.
Inertial body frame (i b -frame): The inertial coordinate system, which coincides with the body frame at the initial time.
The attitude matrix can be decomposed into three parts based on the chain rule as follows: where C n in (t) C ib b (t) denotes the time-varying transformation matrix from the i n -frame to the n-frame and the transformation matrix from the b-frame to the i b -frame, respectively. C in ib is a constant transformation matrix from the i b -frame to the i n -frame. From Equation (1), we can see that the attitude matrix C n b (t) can be divided into three parts: C n in (t) represents the time-varying attitude due to the rotation of the n-frame caused by the earth rotation rate and the vehicle's transport rate related to earth, which can be calculated with the change of latitude and the longitude.
C ib b (t) describes the time-varying attitude due to the rotation of the vehicle from the b-frame to the i b -frame, which can be computed with the angular rate measured by the gyro. It is easy to derive the recursive form: represents the rotation of the b-frame from time t k−1 to time t k . It is appropriate to assume the rotation is tiny while the sampling rate of INS is high.
As for the last part, C in ib can be derived by solving the attitude determination problem using the relationship between vectors: where the V in and V ib are the integration of the specific force in the i n -frame and the i b -frame, respectively. We assume the specific force detected by accelerometer at time t k in the i b -frame is given by following expression: where the δV ib (t k ) is the integration error ofV ib (t k ). Equation (4) could be rewritten as: The attitude matrix C in ib can be written in the form of the Rodrigues vector [29] as: where the l× is the skew-symmetric matrix of the Rodrigues vector. Substituting Equation (7) to Equation (6) leads to: where the w t k is the inertial sensor error.
Defining the sum and difference as follows: We could obtain: Equation (10) is the measurement equation of our filtering model. Once the optimal estimation of the Rodrigues vector l is obtained, we can derive C in ib from Equation (7). Considering the bias and the random noise of the gyro and the accelerometer, the output of gyro and accelerometer can be derived as: Referring to the error model of INS under n-frame, the error equation of the attitude and velocity can be obtained: The gradient of the Rodrigues vector is zero, since the attitude matrix C in ib is a constant matrix.
The gyroscope drifts ε b and the accelerometer bias ∇ b are irrelevant with time, so we can derive: The state vector was selected as follows:

Mstasckf Method
In this section, we demonstrate the algorithm design ideas and the overall process of the algorithm. The algorithm was based on the frame of square-root Kalman filter and the Sage-Husa noise estimator. The aim to our work was to add a criterion to the Sage-Husa noise estimator to ensure the calculated system error variance remains positive, which could make the filter stable. Besides that, a multiple fading factor was implemented to make the residual orthogonal, which could accelerate system convergence and improve the accuracy.

Sage-Husa Noise Estimator
The Sage-Husa filtering method was originally proposed by P.A Sage and G.W Husa in 1969 [30]. The core idea of the algorithm is to construct a real-time filtered noise statistic based on maximum a posteriori information (MAP) to count the system noise and measurement noise in the filtering process. For a nonlinear system as shown in (18): where the w k−1 and v k are the system noise and measurement noise, respectively. f (*) and h(*) denote the nonlinear state model and measurement model, respectively.
Definite the system's posterior error variance matrix at the initial moment as: where the x 0 is the mean value of the initial state vector x 0 . At time t k , assuming the system noise q k , the measurement noise r k , the system noise covariance matrix Q k , and the measurement noise covariance matrix R k are unknown. The state quantity set and quantity measurement set from the initial time to k time are, respectively, It can be seen from the Bayes formula that the probability of estimating the state vector at time t k , the system noise, the measurement noise, and its noise variance matrix from the measurement sequence from the initial time to time t k can be expressed as: where p(Z k ) has nothing to do with noise estimation. Then, (20) can be reduced to the problem of solving the maximum value of probability density as follows: Considering that Cov(q k , r k ) = 0 and each measurement can be regarded as independent.
When the noise changes with time, the filter estimates and measurements at the more recent time can better reflect the current noise situation. Therefore, it is necessary to adjust the weights of the filter estimates and measurements at different times in the noise statistics process to obtain more accurate noise statistics.
Define the weighting coefficient set 1, 2, . . . , k, b as a constant value between 0 to 1. Hence, we can derive a discrete noise estimator for time-varying noise as:q where Φ klk−1 is the one-step prediction matrix from time k − 1 to time k and H k is the measurement model of the system. v k = z k − H kXk|k−1 − r k is the residual of the measurement updating process.

SCKF with Sage-Husa Noise Estimator
There is a problem when applying the Sage-Husa noise estimator to the SCKF frame. Under the framework of square-root volume Kalman filtering, it is necessary to use the Cholesky decomposition form of the system noise variance matrix and the measurement noise variance matrix to perform filtering iterations during the filtering process [24], which requires us to always obtain positive definite in the noise estimator. Otherwise, the filter will be unstable.
Define the second part of (29) and (28) as: When the initial system noise covariance is positive, the sign ofQ k is only related to the sign of J q , and the same situation occurs forR k . To ensure the stability of the filter, we used the equation as follows to computeQ k andR k in the condition J q < 0 or J r < 0.
Unfortunately, when the covariance matrix is computed with (32) or (33), the result is no longer unbiased. Thus, to make up to the accuracy loss caused by the bias, a multiple fading factor was introduced.

Multiple Fading Factor Calculation
The multiple fading factors method is based on the orthogonality principle of the residual. The computation process of multiple fading factors under SCKF could be summarized as follows, and readers can refer to the methodology described by the authors of [31] for more analysis and additional derivations.
Definite the multiple fading factor sequence as: where n is the dimension of the state vector.

MSTASCKF Update
Step 1: Initialization where S k−1 is the lower triangle matrix of the result of the Cholesky decomposition.
Step 2: Prediction In the prediction step, MSTASCKF computes the mean x k|k−1 and S (l) k|k−1 , the squareroot form of the associated covariance without the correction of the multiple fading factors, using the cubature rule, in which each cubature point has the same weight.
Obtain the cubature points: Propagate the cubature points through the state model: Estimate the predicted state: In the proposed algorithm, the weights of the cubature points are equal, which is Square-root factor of the predicted error covariance: where S Q k−1 is the square-root of the system noise Q k−1 at time k − 1 and qr() is the QR decomposition function to obtain the square-root factor of the matrix.
Step 3: Calculating the multiple fading factor Obtain the cubature points: Propagate the cubature points through measurement model: Obtain the predicted measurement: Obtain the multiple fading factor with (35)-(38) and adjust S k|k−1 with multiple fading factors: The l in the upper right corner of the symbol indicates that this quantity has not been corrected by multiple fading factors.
Step 4: Measurement updating Create cubature points: Propagate the cubature points through measurement model: Obtain the Kalman gain: Update the state estimation: Obtain the square-root of noise covariance estimation: Step 5: Estimation of the noise covariance Obtain J q and J r .
Compute the system noise covariance and the measurement noise covariance based on the sign of J q and J r using (28), (29) or (32), and (33).

Semi-Physical Simulation Results and Discussion
In order to assess the performance of the MSTASCKF during the in-flight alignment for GPS/MIMU under vibration, a semi-physical simulation is described in this section. A trajectory maintaining the IMU data and GPS data of typical short-range guided weapon was created, and an experiment to simulate the flight vibration was carried out to obtain simulated vibration noise data. Then, the vibration noise was added to the IMU data to simulate the vibration in-flight. The performance of MSTASCKF was compared with the traditional SCKF, which validated the efficiency of the two strategies added in the MSTASCKF.

Vibration Experiment
Guided weapons generally experience storage, transportation, and flight environments after their production cycle. During flight, the high-speed turbulent flow field formed around the projectile due to high-speed flight induces complex vibration, overload, high temperatures, and other environments [32]. The vibration environment has a greater impact on the structure of the projectile and may even cause deformation of the projectile or shedding of solder joints on the circuit board. The vibration of guided weapons during flight can be divided into two categories: 0-2000 HZ, low-frequency noise excited by vibration mainly transmitted by mechanical structures, and 10-10,000 HZ noise excited by frequency range. Among them, vibration excitation noise has an obvious effect on the mass-spring structure, and excitation noise mainly has an obvious effect on the plate and shell structure. Therefore, the main focus of this paper was the impact characteristics of the frequency range 0-2000 HZ by vibration excitation noise on the output of the inertial navigation system.
In order to obtain the output data of MIMU under the vibration environment, the following random vibration test was designed. The power spectrum density and frequency range of the vibration are as shown in the Figure 1. The experimental platform is shown in Figure 2, which is composed of a MIMU (MTI-1 from the XSENS company), a vibration table, and the tooling for the vibration table. The specification of the MIMU is listed in Table 1.   Since the MTI-1 MIMU cannot be sensitive to the Earth's rotation angular rate, its gyro output under static conditions can be regarded as zero bias and random noise. Vibration noise can be obtained by subtracting the average value of the data at rest when the vibration table is turned on. Thus, the process of the experiment was arranged as in Figure 3 to obtain both the static output and the vibration output of the three axes. The processed vibration noise data in the experiment are shown in Figure 4.

Simulation Condition
Generally, the flight time of short-range guided weapons is within 60 s, and the flight distance is within 20 km. According to the trajectory characteristics of short-range guided weapons, the trajectory attitude and speed used in the semi-physical simulation are shown in the Table 2. The attitude and the velocity of the trajectory are shown in Figure 5. The trajectory was set for 30 s, and the initial attitude angle was set as ψ = 60 • , θ = 45 • , γ = 0 • , representing the yaw, pitch, and roll, respectively. The yaw change was positive from north to east and negative from west to north, the pitch change was positive from the horizontal section of the projectile upward and negative from the downward direction, the hour hand was positive and the counterclockwise was negative, and the projectile speed was the projectile speed at the end of the current time period. In the simulation, the three-axis gyro data were all added with a 100 • /h bias. The speed error of the GPS was 0.1 m/s. Thus, the initial state covariance was set as follows.

Simulation Results under Vibration Noise
The in-flight alignment progresses were calculated using MSTASCKF and traditional SCKF, respectively. The attitude error is shown in Figure 6. The specific value of the estimated attitude error and its RMS of the proposed method and SCKF is shown in Table 3.  We can see that the proposed method had a better RMS performance than the SCKF. The RMS values of the three attitudes of MSTASCKF and SCKF at the end of the alignment process were 1.13 • 0.14 • 5.35 • and 2.5 • 0.99 • 6.12 • , respectively. The RMS of the yaw and roll was much larger than the pitch. This was caused by the estimation scheme, since the MSTASCKF and SCKF had the same trend. The first estimation of roll and pitch had such a large deviation that the subsequent RMS value was affected.
In fact, the main metrics of alignment performance was the final accuracy at the end of the alignment process rather than the accuracy during the process. The attitude error is shown in the Table 4 when the alignment time was limited to 10 s and 20 s. These data could also show that the MSTASCKF had an advantage in the convergence speed and the estimation stability.

Discussion
In this section, a semi-physical simulation for in-flight alignment under vibration noise was carried out to evaluate the performance of the proposed algorithm. Compared with conventional SCKF, it is apparent that the proposed method had a better performance not only in the convergence speed in the beginning of the process but also in the accuracy and steady state. For the proposed filter, the results are as follows.
We proposed a robust SCKF algorithm, which could significantly increase the in-flight alignment accuracy under the vibration noise.
Compared with conventional SCKF filter, the estimation performance was better. The attitude error was smaller and the RMS was also smaller. The convergence time of the proposed method was no longer than 10 s. Thus, the proposed method can meet the rapidity and accuracy requirements of in-flight alignment.

Conclusions
The rapid and accurate in-flight alignment has a vital impact on the performance of the guided weapon applied MIMU. To improve the in-flight alignment performance of MIMU/GPS system under flight vibration, this paper proposed a MSTASCKF filter based on the Sage-Husa noise estimator and multiple fading factors. The filter solved the instability problem of SCKF caused by the Sage-Husa estimator and improved the accuracy. The main work of this paper was as follows.
An experimental vibration environmental test was conducted to maintain the simulated vibration noise data during the flight and a semi-physical simulation was carried out to evaluate the performance of the proposed algorithm. The filtering method had an advantage in the convergence speed, accuracy, and stability compared to SCKF.