A Parameter Self-Calibration Method for GNSS/INS Deeply Coupled Navigation Systems in Highly Dynamic Environments

The GNSS/INS (Global Navigation Satellite System/Inertial Navigation System) navigation system has been widely discussed in recent years. Because of the unique INS-aided loop structure, the deeply coupled system performs very well in highly dynamic environments. In practice, vehicle maneuvering has a big influence on the performance of IMUs (Inertial Measurement Unit), and determining whether the selected IMUs and receiver parameters satisfy the loop dynamic requirement is still a critical problem for deeply coupled systems. Aiming at this, a new parameter self-calibration method based on the norm principle is proposed which explains the relationship between IMU precision and the velocity error of the system; the method will also provide a detailed solution to calculate the loop steady-state tracking error, so it will eventually make a judgment about the stability of the tracking loop under present system parameter settings. Lastly, a full digital simulation platform is set up, and the results of simulations show good agreement with the proposed method.


Introduction
GNSS/INS deeply coupled technology is an important branch of current research in the field of integrated navigation systems. The theory around the loop assistance mechanism and the navigation performance improvement is quite mature, and the deeply coupled system has been applied in the field already. The deeply coupled system introduces the inertial information into the tracking loop of the receiver, which means that the dynamic stress on the receiver could be highly reduced. Hence, this kind of system performs very well in highly dynamic environments. Around this topic, many research institutions and scholars have been working hard to improve the system performance. In terms of loop theory research, Zhang T. proposed an error analysis method to make the connection between most error sources and the PLL (Phase Lock Loop) tracking loop in 2015 [1]. Liu G. studied the relationship between MEMS (Micro Electro Mechanical Systems) IMU precision and tracking loop stability in 2013, and the bias stability of the IMUs was proposed from the angle of usability [2]. Zeng Q. studied the loop control scheme between loop measurements and INS navigation data in 2016, and two control schemes were derived which were proved to be correct [3]. In 2014, Kirkko-Jaakkola studied the jamming mitigation performance of a deeply coupled GNSS-INS system with a low-cost MEMS IMU and gave actual experimental results instead of simulations [4]. Adeel developed a good structure structure for a tracking loop containing a multicarrier VPLL (Vector PLL) in 2015, and the performance of the proposed architecture was validated by means of tests under different conditions [5].
Moreover, the implementation of actual deeply coupled systems has been under development. In 2014, KIT(Karlsruhe Institute of Technology) set up a highly dynamic and deeply coupled verification prototype; a Simulink software receiver was used to complete the algorithm verification on the Matlab platform for the first time, and a real-time solution for deeply coupled architecture was achieved [6]. In 2015, the Rockwell Collins company launched the NavFire-E Series of highperformance integrated navigation systems, which represent the highest standards of the current practical applications. The highly integrated micro-inertial and anti-jamming receiver technology were used to achieve high output accuracy and high reliability of the system [7].
However, the performance of inertial devices will be affected seriously by highly dynamic environments. The key point is to make sure that the assistance information used to reduce the dynamic stress is reliable. How to select correct inertial devices and receiver parameters to satisfy the loop dynamic requirement is still a problem for the system. However, quite a few articles have studied this: Liu G. discussed the usability of a MEMS IMU in a deeply coupled system, but the gsensitivity error of gyroscope under high dynamic environments was ignored [2]; Kirkko-Jaakkola built an actual experimental platform with a MEMS IMU, but it was not clear how to select the appropriate IMU devices [4]. In this paper, a mathematic model of the deeply coupled system including the g-sensitivity error of gyroscope is introduced. Then, according to the norm principle, the error propagation properties are analyzed by mathematical deduction and this process gives the quantitative velocity error caused by IMU error. In addition, a new method which combines with the tracking error model to determine the loop stability is introduced. It will eventually make a judgment about the stability of the tracking loop under certain parameter settings. Finally, the correctness of the proposed method was verified by several simulations based on a digital simulation platform.

Mathematical Model of a Deeply Coupled Navigation System
The core concept of the deeply coupled system aided by an IMU is to extract the doppler frequency and code phase deviation based on the inertial information and satellite ephemeris, and this can help to control the code/carrier NCO (Numerically Controlled Oscillator) in the tracking loop. This structure can ensure the system maintains high performance of the tracking loop in narrow bandwidth conditions. The overall architecture of the system is shown in Figure 1.  The performance of the inertial devices is obviously affected by highly dynamic circumstances, so the IMUs' error model should be established first. Due to the smaller and lighter tendency of the deeply coupled system, this article mainly discusses the MEMS device.
The error model of gyroscope and accelerometer are given as follows: The performance of the inertial devices is obviously affected by highly dynamic circumstances, so the IMUs' error model should be established first. Due to the smaller and lighter tendency of the deeply coupled system, this article mainly discusses the MEMS device.
The error model of gyroscope and accelerometer are given as follows: where ω b ib is the angular velocity measurement of the gyroscope, i denotes the inertial coordinate system, b denotes the body coordinate system, b g is the gyroscope bias, g s is the matrix of g-sensitivity coefficients, f b ib is the specific force, η g is the measurement error of gyroscope, f b ib is the specific force measurement of the accelerometer, and b a is the accelerometer bias [8,9].
In order to analyze the relationship between the IMU errors and the INS navigation errors in the deeply coupled system, the state variables related to the inertial system in the integrated navigation filter could be expressed as follows: where δR e are the error vectors of the 3D position, δV e are the error vectors of 3D velocity, and δψ e are the error vectors of 3D attitude; e denotes the ECEF coordinate system; δω b are the error vectors of the gyroscopes; and δA b are the error vectors of the accelerometers in the body coordinate system. δT are the clock error vectors including clock bias error δT bias and clock drift error δT dri f t [10,11]. The equation of state is given as follows: where G denotes the parameter of gravity, µ = 3.986005 × 10 14 (m 3 /s 2 ), which represents the product of the gravitation constant and the mass of earth; x y z denote the position vectors in ECEF coordinates, r = x 2 + y 2 + z 2 , which is the distance to the center of the Earth; Ω e is the matrix of the Earth's rotational angular velocity, ω e = 7.292115 × 10 −5 (rad/s); F e is the matrix of specific force in the ECEF coordinate system; C e b is the direction cosine matrix between the body coordinate system and the ECEF coordinate system; Λ is the clock error coefficient matrix; ω T = 0 ω f T is the clock error vectors; and ω f is the clock drift error vectors.
To make sure that the system stability could be analyzed, Wagner J. F. studied the observability of the integrated system [12]. According to his conclusion, when the system consists of three antennas which form a triangle of sufficient area, observability would be guaranteed under all circumstances. Thus, this paper continues on the premise that the observability of the system is guaranteed by the system hardware.

Parameter Self-Calibration Method for Loops in Deeply Coupled Systems Based on Norm Analysis
According to the analysis in the above section, the position and velocity of the vehicle calculated by INS are directly involved in the feedback to the receiver tracking loop, and the carrier frequency error is mainly related to the velocity error, while the code phase error is only related to the position error. When the closed-loop Kalman filter is run to the steady state, the estimation error is convergent, and the accuracy of the filtering result is mainly determined by the covariance matrix of the state noise. Therefore, what needs to be considered is the INS error in one calculation cycle. Because the code phase error caused by the position error is much smaller than the carrier frequency error and leads to a quite small influence on the system, this article focuses on the discussion of the INS velocity error analysis.

Norm Analysis to the IMU Error Propagation Properties
Before the norm analysis referring to the INS error, the definition of the two kinds of norms should be made clear [13]. If A is a vector, the Euclid norm, 2 ) 1/2 . Both of these two norms represent the meaning of spatial distance and satisfy the following axioms: (1) The constant coefficients theorem is described as follows: Let A be an matrix of m × m scalars, X be an m × n matrix of unknowns, and B be an m × n matrix of given functions. Then the differential equation with initial value X(t 0 ) = X 0 has the solution where t 0 denotes the initial moment. According to the theorem above, Equation (4) could be expressed briefly as the following differential equation: and the approximate analytical solution of Equation (7) is where ∆T denotes the time interval between two consecutive estimations. From Equation (4), the velocity and position differential equation could be expressed as follows: From Equation (11), the position error δR e in one time interval should be δV e ∆T. Thus, Equation (12) could be expressed as follows: Combining with Equation (10), the velocity error caused by the accelerometer error is Thus, from Equation (14), the velocity error caused by the accelerometer error is and it should be noted that the small quantities of high order (∆T ≥4 ) in this step have been ignored. Both sides of Equation (16) could be expressed in the norm form as follows: According to the norm property, Equation (17) could be transformed into the following expression: where C e b is the unitary matrix, so C e b F = tr(C b e C e b )) = √ 3. Equation (18) could be transformed into the following expression: because Ω e F = √ 2ω e and G F = √ 6µ r 2 , according to the definition of F-norm. Also according to Equation (4), the velocity error caused by the gyro error could be concluded as follows with the same analysis idea above: where both sides of Equation (20) could be expressed in the norm form as follows: From Equations (19) and (21), the relationship between the IMU errors and the related velocity error is clear. It could be concluded that the INS velocity error mainly depends on four factors: the gyro error δω b , the accelerometer error δA b , the time interval ∆T, and the vehicle maneuvering which can be expressed as f 2 x + f 2 y + f 2 z . Through the mathematical deduction above, the quantitative velocity error can be calculated directly when the IMUs' precision and the dynamic requirement for the system are determined. This will do great favors for the stability judgement of the tracking loop in the next step.

Parameter Self-Calibration Method for the Tracking Loop in a Deeply Coupled System
A method based on the norm property was proposed to analyze the relationship between the IMU errors and the system velocity error above; the next step is to make the judgement on the stability of the tracking loop in deeply coupled systems in highly dynamic environments.
For the deeply coupled system, the total tracking error of carrier loop could be given as follows: θ e = 360 a n n R (n) where σ tPLL is the standard deviation of phase error caused by the thermal noise, σ v is the standard deviation of phase error caused by the jitter of the reference oscillator frequency, σ A is the standard deviation of phase error caused by the frequency drift of the Allan crystal oscillator, and θ e is the dynamic stress error caused by the relative movement. Further, B L is the equivalent bandwidth of the loop, C/N 0 is the carrier-to-noise ratio, T coh is the coherent integration time, c is the speed of light, λ is the carrier wavelength, σ A (τ) is the Allan standard deviation of the frequency, and a n is filter parameters of the loop; n is the loop order, a 1 = 0.25, a 2 = 0.53, a 3 = 0.7845 [8].
It needs to be emphasized that R (n) = d n R dt n , where R (n) is the nth-order derivative of the distance between the vehicle and the satellite.
Because the second-order locked loop is sensitive to the accelerated movement and the third-order loop is sensitive to the jerk movement, the velocity error should be transformed to the equivalent relative movement. Assuming that the INS-aided velocity error in one loop-aided cycle is δv e , the (n − 1)th-order derivative is as follows: where R (n) represents the nth-order derivative of the position error vectors. In order to keep the tracking loop stable, 3 times the variance of the phase measurement error should be less than 45 • [14]. According to Equations (22), (25) and (26), the following inequation could be inferred: Furthermore, according to Equations (19) and (21), the (n − 1)th-order derivative of the INS-aided velocity error is where σ v is about 2 • , n is the loop order, a is the sum acceleration of the vehicle, and ∆T is the time interval between two INS-aided moments. So far, all the key parameters of the receiver are involved in Equations (27) and (28), and the whole deeply coupled system can be configured to satisfy the engineering requirement in consideration of the vehicle maneuvering. If the IMU precision related to δA b /δω b and the dynamic requirement of the deeply coupled system related to a are defined, it can be directly determined whether the system can guarantee the stability of the tracking loop or not. The entire parameter self-calibration method is proposed above; the correctness and advantages will be proved next.

Highly Dynamic Simulations and Results
In order to prove the correctness and advantages of the proposed method, a full digital simulation platform based on Matlab was constructed. The orbital distribution and operation cycle of GPS (Global Position System) is generated according to actual GPS parameters; the IMUs' information is selected from the data book of Analog Devices Inc. In this part, the norm-based analysis method is verified first, and after combining with three different kinds of MEMS gyros' parameters, the actual result of the parameter self-calibration method is put forward.

Simulation Conditions and Track Settings
According to the data book of ADI Co., Ltd. product, three typical gyroscopes are selected; they are the high-precision gyroscope ADIS16490, the middle-precision gyroscope ADIS16448, and the low-precision gyroscope ADIS16300 (Analog Devices, Norwood, MA. USA). The main error parameters of the gyroscopes are shown in Table 1, and the bias errors of accelerometers are set as 16 mg/20 mg/60 mg. The parameters of the GPS system are as follows: the carrier-to-noise ratio (CNR) is 50 dB-Hz, the bandwidth of the carrier tracking loop is 20 Hz, the coherent integration time is 1 ms, the Allan standard deviation of the frequency is 10 −6 Hz, and the carrier wavelength of GPS L1 is 0.19 m. For the deeply coupled system, the calculation cycle of INS and the cycle of INS-aided loop control are both 10 ms, and the filter cycle time is 0.1 s. A typical dynamic trajectory of an X-43 hypersonic aircraft is set as the reference track from running take-off until level flight [15]. The initial position is in Nanjing where the coordinate is (18.80 • , 32.03 • , 10.00 m), the initial velocity is 0 m/s, and the initial heading angle is 90 • . The 3D trajectory is as shown in Figure 2, and each stage of the dynamic trajectory is as shown in Table 2.

Simulation of Inertial Calculation
In this part, the correctness of the proposed norm analysis method will be proved from the aspect of inertial calculation. The calculated velocity error by norm analysis is compared with the Monte Carlo simulation result of the inertial calculation in one cycle.

Norm Analysis Simulation
In order to prove the correctness of the norm analysis method above, the first 10 ms of the acceleration of the X-43 are picked up; these 10 ms give the greatest dynamic movement in the whole trajectory. The initial velocity of the vehicle is set as (2146.4 m/s, 0, 0). Ten Monte Carlo simulations were conducted and the root-mean-square errors (RMSEs) of three velocity errors were compared with the calculated velocity error from Equation (28). A histogram of the RMSE is shown in Figure 3, and the detailed statistical analysis of Figure 3 is shown in Table 3.

Simulation of Inertial Calculation
In this part, the correctness of the proposed norm analysis method will be proved from the aspect of inertial calculation. The calculated velocity error by norm analysis is compared with the Monte Carlo simulation result of the inertial calculation in one cycle.

Norm Analysis Simulation
In order to prove the correctness of the norm analysis method above, the first 10 ms of the acceleration of the X-43 are picked up; these 10 ms give the greatest dynamic movement in the whole trajectory. The initial velocity of the vehicle is set as (2146.4 m/s, 0, 0). Ten Monte Carlo simulations were conducted and the root-mean-square errors (RMSEs) of three velocity errors were compared with the calculated velocity error from Equation (28). A histogram of the RMSE is shown in Figure 3, and the detailed statistical analysis of Figure 3 is shown in Table 3. Sensors 2018, 18, x 9 of 14  Table 3. RMSEs of the velocity errors in one cycle.  Figure 3 and Table 3, the velocity error calculated by Equation (28) matches well with the RMSE of the velocity errors from the Monte Carlo simulation. Hence, the proposed norm analysis method is correct and the calculated velocity error can be trusted to work in the loop stability analysis. Further, the calculation times of the norm analysis method and the Monte Carlo method were compared by Matlab platform as performed on a PC equipped with an I7-7700 CPU; the result shows that the norm analysis method only takes 0.01 s to get the velocity error while the Monte Carlo method needs 0.42 s. So, from this part, we can see that the proposed norm analysis method is more concise, efficient, quantitative, and precise.

Stability Simulation under Highly Dynamic Circumstances
With the help of the simulation above, the velocity error caused by the inertial calculation has been quantized, and the stability of the tracking loop is related to Equation (27) directly. So, in this part, the correctness of Equation (27) will be proved.
The parameters of both the INS and GNSS systems are set as the same as before; it should be mentioned that the Kalman Filter (KF) parameters should be configured according to the system and measurement parameters, and this article will not discuss the setting of KF parameters. The structure of the carrier loop is determined as the third-order phase-locked loop aided by the second-order frequency-locked loop, so n in Equation (27) Table 4.  Table 3. RMSEs of the velocity errors in one cycle. From Figure 3 and Table 3, the velocity error calculated by Equation (28) matches well with the RMSE of the velocity errors from the Monte Carlo simulation. Hence, the proposed norm analysis method is correct and the calculated velocity error can be trusted to work in the loop stability analysis. Further, the calculation times of the norm analysis method and the Monte Carlo method were compared by Matlab platform as performed on a PC equipped with an I7-7700 CPU; the result shows that the norm analysis method only takes 0.01 s to get the velocity error while the Monte Carlo method needs 0.42 s. So, from this part, we can see that the proposed norm analysis method is more concise, efficient, quantitative, and precise.

Stability Simulation under Highly Dynamic Circumstances
With the help of the simulation above, the velocity error caused by the inertial calculation has been quantized, and the stability of the tracking loop is related to Equation (27) directly. So, in this part, the correctness of Equation (27) will be proved.
The parameters of both the INS and GNSS systems are set as the same as before; it should be mentioned that the Kalman Filter (KF) parameters should be configured according to the system and measurement parameters, and this article will not discuss the setting of KF parameters. The structure of the carrier loop is determined as the third-order phase-locked loop aided by the second-order frequency-locked loop, so n in Equation (27) Figure 4 shows that the high-precision IMU helps maintain tracking through almost all the process. Only at the times of 265 s and 270 s does there exist a phase detecting error saltation, because the acceleration changed from 19 m/s/s to 300 m/s/s suddenly and the tracking loop cannot adjust to the great jerk instantly. During the 5 s acceleration of X-43, the high-precision IMU can help maintain tracking, and the middle-precision IMU performs worse than the high-precision IMU, but the lowprecision IMU cannot guarantee the stability of the tracking loop even when the vehicle is under a low dynamic environment. Further, according to Table 4, the calculated phase errors show that both high-precision and middle-precision IMUs could satisfy the stability condition of Equation (27) but the low-precision IMU could not; this result agrees well with the phenomenon observed in Figure 4.    Figure 4 shows that the high-precision IMU helps maintain tracking through almost all the process. Only at the times of 265 s and 270 s does there exist a phase detecting error saltation, because the acceleration changed from 19 m/s/s to 300 m/s/s suddenly and the tracking loop cannot adjust to the great jerk instantly. During the 5 s acceleration of X-43, the high-precision IMU can help maintain tracking, and the middle-precision IMU performs worse than the high-precision IMU, but the low-precision IMU cannot guarantee the stability of the tracking loop even when the vehicle is under a low dynamic environment. Further, according to Table 4, the calculated phase errors show that both high-precision and middle-precision IMUs could satisfy the stability condition of Equation (27) but the low-precision IMU could not; this result agrees well with the phenomenon observed in Figure 4.
The navigation results aided by the three kinds of IMUs are further compared in Figures 5 and 6. The detailed statistical analysis of Figures 5 and 6 is shown in Table 5.  Table 5. RMSEs of the deeply coupled system errors aided by the three kinds of IMUs.

RMSE of System Errors High-Precision Mid-Precision Low-Precision
East-position (m)   Table 5. RMSEs of the deeply coupled system errors aided by the three kinds of IMUs.

RMSE of System Errors High-Precision Mid-Precision Low-Precision
East-position (m)   Table 5 show that with the help of the high-precision IMU, the deeply coupled navigation system performs very well in highly dynamic environments. The middleprecision IMU can guarantee the navigation precision of the system but performs worse than the high-precision IMU. However, the low-precision IMU causes the system to lose positioning ability throughout the whole process.
In all, the simulation results above have proved that Equation (27) could make a correct judgment about the stability of the tracking loop in consideration of the IMUs' precision, the system parameters, and the vehicle maneuvering.

Conclusions
A new mathematical model of a deeply coupled system including the g-sensitivity error of the gyroscope was proposed. According to the system model, the error propagation properties were analyzed by mathematical deduction from the perspective of the norm principle. Then, an inequality considering most of the deeply coupled navigation system parameters was proposed to judge the  Table 5 show that with the help of the high-precision IMU, the deeply coupled navigation system performs very well in highly dynamic environments. The middle-precision IMU can guarantee the navigation precision of the system but performs worse than the high-precision IMU. However, the low-precision IMU causes the system to lose positioning ability throughout the whole process.
In all, the simulation results above have proved that Equation (27) could make a correct judgment about the stability of the tracking loop in consideration of the IMUs' precision, the system parameters, and the vehicle maneuvering.

Conclusions
A new mathematical model of a deeply coupled system including the g-sensitivity error of the gyroscope was proposed. According to the system model, the error propagation properties were analyzed by mathematical deduction from the perspective of the norm principle. Then, an inequality considering most of the deeply coupled navigation system parameters was proposed to judge the stability of the tracking loop. All the above theories were verified by a full digital simulation platform, and the simulation results proved the correctness of the relevant conclusions.
The proposed method can provide basic guidance for the design of deeply coupled systems. The stability judging process is reliable and flexible, and so can provide new ideas for the engineering realization of deeply coupled navigation systems. Also, the actual effect of the theory in practical systems will be checked in the near future.