Coupled Integration of CSAC, MIMU, and GNSS for Improved PNT Performance

Positioning, navigation, and timing (PNT) is a strategic key technology widely used in military and civilian applications. Global navigation satellite systems (GNSS) are the most important PNT techniques. However, the vulnerability of GNSS threatens PNT service quality, and integrations with other information are necessary. A chip scale atomic clock (CSAC) provides high-precision frequency and high-accuracy time information in a short time. A micro inertial measurement unit (MIMU) provides a strap-down inertial navigation system (SINS) with rich navigation information, better real-time feed, anti-jamming, and error accumulation. This study explores the coupled integration of CSAC, MIMU, and GNSS to enhance PNT performance. The architecture of coupled integration is designed and degraded when any subsystem fails. A mathematical model for a precise time aiding navigation filter is derived rigorously. The CSAC aids positioning by weighted linear optimization when the visible satellite number is four or larger. By contrast, CSAC converts the GNSS observations to range measurements by “clock coasting” when the visible satellite number is less than four, thereby constraining the error divergence of micro inertial navigation and improving the availability of GNSS signals and the positioning accuracy of the integration. Field vehicle experiments, both in open-sky area and in a harsh environment, show that the integration can improve the positioning probability and accuracy.


Introduction
Positioning, navigation, and timing (PNT) is important for vehicles, aircrafts, robots, and pedestrians, especially for military use. Global navigation satellite systems (GNSS) based on radio signals and inertial navigation systems (INS) based on Newton's law are widely used forms of PNT technology. GNSS can provide automatic, continuous, low-cost, global coverage in all-weather with high-precision positioning and timing. However, the satellites are thousands of kilometers away from the ground and the signal transmission power is limited by the solar power collector. Thus, satellite signals arriving at the ground are extremely weak and susceptible to occlusion and interference, which can cause navigation interruption. The poor dynamic performance of GNSS makes it difficult to provide continuous PNT information to high-speed moving carriers [1][2][3][4]. INS measures acceleration and angular velocity by inertial sensors, and calculates position, velocity, and attitude by an inertial navigation algorithm. INS provides autonomous, high refresh rate, good short-term accuracy, and stable navigation information. INS can work in the air, ground, underground, underwater, and indoors. With the development of micro-electromechanical systems (MEMS) technology, the performance of

Architecture of Coupled Integration
The architecture of CSAC, MIMU, and GNSS coupled integration is described in Figure 1; it includes four blocks, namely, timing, GNSS, INS, and precise time-aiding navigation filters.
The timing block provides timing information. Usually, when the GNSS constellation geometry is good, one pulse per second (1 PPS) signal is synchronized with Coordinated Universal Time (UTC) with an error of only several nanoseconds. Thus, the start time of CSAC is determined by GNSS. One pulse per second signal enters the timing block and uses a feedback loop to restrain error accumulation caused by "clock coasting" of CSAC. This loop also reduces the noise of 1 PPS signal, and it consists of a digital phase detector, regulator, CSAC, and frequency division. If GNSS cannot provide 1 PPS signal, then this loop degrades to clock coasting mode.
In the INS block, carrier motion can be described as r u " xr a, r ωy, which includes three-dimensional acceleration and angular rate. MIMU measures the carrier motion. After sampling and calibration, the measurement can be expressed as followsû k " g 1 pr u k ,ξ k q ξ k "ξḱ`δξ k (1) Initial calibration parameters of MIMU can be determined by laboratory tests or a fast field calibration method. When the navigation filter works, it provides the parameters.
The strap-down inertial navigation algorithm iŝ rḱ " f pr k ,û k q r k " g 2 prḱ , δr k q In the GNSS block, the clock signal is provided by the timing block. The antenna, Radio Frequency (RF) front, down conversion, and analog to Digital Converter (ADC) transfer the satellite RF signals to base band signal. Acquisition of base band signal is aided by the position, velocity, and time information from the navigation filter. ∆ρ and ∆ . ρ are pseudo range increment and pseudo range rate increment, respectively, which are from the navigation filter. The delay lock loop (DLL) used for pseudo random code tracking is aided by ∆ρ and phase lock loop (PLL) used for carrier tracking is aided by ∆ . ρ. These aids can improve the dynamic stress tolerance of GNSS.
Initial calibration parameters of MIMU can be determined by laboratory tests or a fast field calibration method. When the navigation filter works, it provides the parameters.
The strap-down inertial navigation algorithm is In the GNSS block, the clock signal is provided by the timing block. The antenna, Radio Frequency (RF) front, down conversion, and analog to Digital Converter (ADC) transfer the satellite RF signals to base band signal. Acquisition of base band signal is aided by the position, velocity, and time information from the navigation filter.  In the precise time-aiding navigation filter block, precise time, pseudo range, pseudo range rate, and INS output are transmitted to the block. Process and observation models are constructed based on precise time, and the Kalman filter is used to integrate navigation information. The calibrated information is transmitted to sensors and the error parameters are corrected. The PNT information is provided as the integration's output.    In the precise time-aiding navigation filter block, precise time, pseudo range, pseudo range rate, and INS output are transmitted to the block. Process and observation models are constructed based on precise time, and the Kalman filter is used to integrate navigation information. The calibrated information is transmitted to sensors and the error parameters are corrected. The PNT information is provided as the integration's output.
The coupled integration can be degraded because of the failure of subsystems. Table 1 shows degradations of the coupled integration. In the table, " ' " or "ˆ" means that the sensor works or fails, respectively. If all of the subsystems work well, then the integration works in coupled integration mode. If CSAC fails, then the integration works in MIMU/GNSS tightly coupled integration. If MIMU fails, then the integration works in CSAC/GNSS coupled integration mode. If GNSS fails, then the integration works in inertial navigation mode.

Mathematical Model of Precise Time-Aiding Navigation Filter
The precise time-aiding navigation filter is based on the Kalman filter. The integration state model is constructed by inertial sensors and algorithms, and the integration observation model is constructed by GNSS information with precise time aiding. The navigation filter is defined as Λ.

Integration State Model
The state vector X is defined as: where δR " " δϕ δλ δh ı is the error of latitude, longitude, and altitude of MEMS INS, is the error of east, north, and up velocity, Φ " " α β γ ı is the error of pitch, roll, and yaw angles, ∇ " is the three-axis accelerometer bias, " " " ε x ε y ε z ı is the three-axis gyroscope bias, δSF a " " δSF ax δSF ay δSF az ı is the scale factor error of the accelerometers, and δSF g " " δSF gx δSF gy δSF gz ı is the scale factor error of the gyroscopes. Unlike traditional tightly coupled integrated systems, CSAC provides precision time aiding. Clock bias and drift are not regarded as state variables.
Thus, the state equation can be written as: Based on error propagation function of MEMS INS in the navigation coordinate, A can be written in the Appendix.

Integration Observation Model
An observation model is described by the range or velocity error between satellites and the integration device. Suppose that the real position of the device is R 0 " " x 0 y 0 z 0 ı T and the calculated position of the device is R MINS " " x MINS y MINS z MINS ı T , there are N visible satellites named "1, 2¨¨¨N" and the position of the satellite marked "n" is R pnq " " x pnq y pnq z pnq ı T . The position error is: The range between the device and the satellite n is derived based on calculated position.
Meanwhile, the real range is: Taylor expansion of Equation (7) is deformed by neglecting higher order terms.
CSAC clock coasting provides accurate time information in the integration. Thus, the GNSS receiver provides the range between the device and the satellite "n" as follows: The range error of ρ pnq GNSS can be divided into two parts (i.e., ∆ρ pnq c and ε δt u ). ∆ρ pnq c is the error described by factors except the receiver clock noise, which include the satellite clock offset, ionosphere delay, troposphere delay, ephemeris errors, and multipath error. ε δt u is the pseudo range error caused by the receiver clock noise. Hence, If N ě 4, then CSAC aids positioning by weighted linear optimization. The covariance of ∆ " ρ c is: where σ 2 clock is the receiver clock noise variance, and σ 2 URE is the range error covariance of ρ pnq GNSS , which does not include the receiver clock error, I is a NˆN unit matrix, and O " If N ă 4, then CSAC converts the GNSS observations to range measurements by "clock coasting", which can constrain error divergence of micro inertial navigation, as well as improve the GNSS signals availability and positioning accuracy of the integration.
where c is the speed of light, B f 0 is the frequency bias of CSAC, B f 1 is the frequency drift of CSAC, and t´t 0 represents clock coasting time. The range measurement error is: where » -- C LLH ECEF is the transform matrix from earth-centered earth-fixed (ECEF) coordinate system to site-centered coordinate system. Hence, When N satellites are visible, the range measurement error equation is written as: where The derivative of r pnq 0 with time is: .
The derivative of ρ pnq MINS with time is: The derivative of range between device and the satellite "n" from GNSS receiver with time is: where, When the number of visible satellite is N, the time derivative of range error is written as follows: where We then define the observation equation Combining the range error and range rate error to measurement vector, we obtain: The observation matrix is: and the observation noise is:

Experimental Setup
The demo setup of the coupled integration consisted of several parts: SA.45 CSAC, STIM300 MIMU, GPS L1 intermediate frequency signal collector, navigation processor, portable power, and several interfaces. The hardware structure of the demo setup is shown in Figure 2, and an actual photograph is shown in Figure 3. The main parameters of sensors in the demo setup are listed in Table 2.

Experimental Setup
The demo setup of the coupled integration consisted of several parts: SA.45 CSAC, STIM300 MIMU, GPS L1 intermediate frequency signal collector, navigation processor, portable power, and several interfaces. The hardware structure of the demo setup is shown in Figure 2, and an actual photograph is shown in Figure 3. The main parameters of sensors in the demo setup are listed in Table 2.

Experimental Setup
The demo setup of the coupled integration consisted of several parts: SA.45 CSAC, STIM300 MIMU, GPS L1 intermediate frequency signal collector, navigation processor, portable power, and several interfaces. The hardware structure of the demo setup is shown in Figure 2, and an actual photograph is shown in Figure 3. The main parameters of sensors in the demo setup are listed in Table 2.     An optical fiber combined with the inertial navigation system GI7660, which is manufactured by Beijing StarNeto Technology Development Co., Ltd. (Beijing, China), provided the reference trajectory. GI7660 consists of three close-loop optical fiber gyroscopes, three quartz accelerometers, and a mobile mapping grade multi-constellation multi-frequency GNSS receiver. The antenna of the GNSS receiver is HX-CS5601A, which is manufactured by Shenzhen Harxon Antenna Technology Co., Ltd. (Shenzhen, China). The main parameters of GI7660 are shown in Table 3, and an actual photograph is shown in Figure 4. The bias stability of GI7660 gyroscopes is 0.3˝/h, whereas that of STIM300 is 6˝/h with 10 s under average condition. The bias stability of GI7660 accelerometer is 20 µg, whereas that of STIM300 is 70 µg with 10 s under average conditions. Therefore, the inertial performance of GI7660 is much better than that of STIM300. In the tests, GI7660 worked in Real-Time Kinematic (RTK) mode.
The device for comparison consisted of a STIM300 and GPS L1 receiver driven by OCXO. The main parameters of OCXO are shown in Table 4.   The device for comparison consisted of a STIM300 and GPS L1 receiver driven by OCXO. The main parameters of OCXO are shown in Table 4. The demo setup, the reference device, and the compared setup were fixed in a vehicle, as shown in Figure 5, in which the same GNSS antenna is employed (i.e., HX-CS5601A, Shenzhen Harxon The demo setup, the reference device, and the compared setup were fixed in a vehicle, as shown in Figure 5, in which the same GNSS antenna is employed (i.e., HX-CS5601A, Shenzhen Harxon Antenna Technology Co., Ltd., Shenzhen, China).
The navigation result refresh rate of the GI7660 is 10 Hz and was marked as "Ref". The sample rate of the MIMU in the demo setup was 200 Hz. The inertial navigation result refresh rate of the MIMU with the two-sample iteration algorithm was 100 Hz and was marked as "CSAC+MIMU". The positioning result refresh rate of the GPS L1 receiver in the demo setup was 1 Hz and was marked as "CSAC+GNSS". The navigation result refresh rate of the demo with the integration algorithm aided by a precise clock was 100 Hz and was marked as "CSAC+MIMU+GNSS". The navigation result refresh rate of the compared setup was 100 Hz and was marked as "MIMU+GNSS". rate of the MIMU in the demo setup was 200 Hz. The inertial navigation result refresh rate of the MIMU with the two-sample iteration algorithm was 100 Hz and was marked as "CSAC+MIMU". The positioning result refresh rate of the GPS L1 receiver in the demo setup was 1 Hz and was marked as "CSAC+GNSS". The navigation result refresh rate of the demo with the integration algorithm aided by a precise clock was 100 Hz and was marked as "CSAC+MIMU+GNSS". The navigation result refresh rate of the compared setup was 100 Hz and was marked as "MIMU+GNSS".

Open-Sky Route
Field vehicle experiments of open-sky route were conducted in Beiqing Road, Haidian District, Beijing, China. Figure 6 shows a map of the experiment route. The start point is marked by a red

Open-Sky Route
Field vehicle experiments of open-sky route were conducted in Beiqing Road, Haidian District, Beijing, China. Figure 6 shows a map of the experiment route. The start point is marked by a red marker in the map (Figure 6b). The vehicle ran along the direction indicated by the black arrow and back to the start point. The journey was about 5.5 km. Before the experiments, MIMU was calibrated by the fast field calibration method. Then, the setup was fixed in the vehicle and kept stationary for 10 s for initial alignment. The vehicle started and ran along the predetermined route until it reached the terminal. The reference device GI7660 tracked the GPS (L1, L2, and L2C), BD2, GLONASS (L1 and L2), and Galileo (E1) constellation signals. Figure 7 shows the tracked satellite number of GI7660, and the average was 15. Figure 8 shows the reference information of the open-sky route measured by GI7660, including displacements and altitude with time. The reference device GI7660 tracked the GPS (L1, L2, and L2C), BD2, GLONASS (L1 and L2), and Galileo (E1) constellation signals. Figure 7 shows the tracked satellite number of GI7660, and the average was 15. Figure 8 shows the reference information of the open-sky route measured by GI7660, including displacements and altitude with time.    Figure 9 shows the tracked satellite number and space distribution of the GPS L1 receiver in the demo and compared setups. The average number of tracked satellites was less than 7. PRN3, PRN4, PRN16, PRN26, PRN29, and PRN31 were tracked for the entire route.
Sensors 2016, 16,682 13 of 20 Figure 9 shows the tracked satellite number and space distribution of the GPS L1 receiver in the demo and compared setups. The average number of tracked satellites was less than 7. PRN3, PRN4, PRN16, PRN26, PRN29, and PRN31 were tracked for the entire route. The navigation results are compared in Figure 10 and Table 5. The result of the coupled integration navigation of the demo setup "CSAC+MIMU+GNSS" was better than that obtained by any other form. The horizontal positioning accuracy of "CSAC+MIMU+GNSS" was 13% higher than that of "CSAC+GNSS" and 21% higher than that of "MIMU+GNSS". The navigation results are compared in Figure 10 and Table 5. The result of the coupled integration navigation of the demo setup "CSAC+MIMU+GNSS" was better than that obtained by any other form. The horizontal positioning accuracy of "CSAC+MIMU+GNSS" was 13% higher than that of "CSAC+GNSS" and 21% higher than that of "MIMU+GNSS".
Sensors 2016, 16,682 13 of 20 Figure 9 shows the tracked satellite number and space distribution of the GPS L1 receiver in the demo and compared setups. The average number of tracked satellites was less than 7. PRN3, PRN4, PRN16, PRN26, PRN29, and PRN31 were tracked for the entire route. The navigation results are compared in Figure 10 and Table 5. The result of the coupled integration navigation of the demo setup "CSAC+MIMU+GNSS" was better than that obtained by any other form. The horizontal positioning accuracy of "CSAC+MIMU+GNSS" was 13% higher than that of "CSAC+GNSS" and 21% higher than that of "MIMU+GNSS".  A software mask was used to reduce the visible satellite number in the open-sky route. From 100 s to the final, only three satellites could be tracked because of the software mask. The visible satellites were PRN3, PRN16, and PRN29. The comparison of the navigation results by software mask is shown in Figure 11 and Table 6. After application of the software mask, "CSAC+GNSS" could still provide position information, and the positioning errors increased gradually with time. The positioning errors of "MIMU+GNSS" increased quickly with time, reaching >200 m at the end. The positioning errors of "CSAC+MIMU+GNSS" increased slowly and the horizontal error was 9.76 m, which was 14% more accurate than "CSAC+GNSS".  A software mask was used to reduce the visible satellite number in the open-sky route. From 100 s to the final, only three satellites could be tracked because of the software mask. The visible satellites were PRN3, PRN16, and PRN29. The comparison of the navigation results by software mask is shown in Figure 11 and Table 6. After application of the software mask, "CSAC+GNSS" could still provide position information, and the positioning errors increased gradually with time. The positioning errors of "MIMU+GNSS" increased quickly with time, reaching >200 m at the end. The positioning errors of "CSAC+MIMU+GNSS" increased slowly and the horizontal error was 9.76 m, which was 14% more accurate than "CSAC+GNSS".

Harsh Route
Field experiments in harsh environments were conducted, and the route map is shown in Figure 12. The journey was about 11 km and lasted for approximately 1300 s. Figure 13 shows the number of tracked satellites with time for GI7660, the demo, and the compared setups. The average number for GI7660 was 12. The average number of the demo and compared setups was 4.6. At 1030-1113 s, severe occlusion occurred and less than three satellites could be tracked.

Harsh Route
Field experiments in harsh environments were conducted, and the route map is shown in Figure 12. The journey was about 11 km and lasted for approximately 1300 s. Figure 13 shows the number of tracked satellites with time for GI7660, the demo, and the compared setups. The average number for GI7660 was 12. The average number of the demo and compared setups was 4.6. At 1030-1113 s, severe occlusion occurred and less than three satellites could be tracked. Displacements and errors from the harsh route are shown in Figure 14, and navigation results are compared in Table 7. In Table 7 the statistical range in time domain of the "CSAC+GNSS" results was 89% of the whole route and the occlusion area is not included. The positioning probability of CSAC+GNSS was 89%. However, the statistical range in time domain of the "CSAC+MIMU+GNSS" results was 100% of the whole route and the occlusion area is included. The positioning error of the occlusion area expanded. This means that the statistical errors of "CSAC+MIMU+GNSS" are larger than "CSAC+GNSS". The coupled integration positioning probability of the demo setup "CSAC+MIMU+GNSS" was 11% higher than that of "CSAC+GNSS". The horizontal positioning accuracy of "CSAC+MIMU+GNSS" was 24% higher than that of "MIMU+GNSS".  Displacements and errors from the harsh route are shown in Figure 14, and navigation results are compared in Table 7. In Table 7 the statistical range in time domain of the "CSAC+GNSS" results was 89% of the whole route and the occlusion area is not included. The positioning probability of CSAC+GNSS was 89%. However, the statistical range in time domain of the "CSAC+MIMU+GNSS" results was 100% of the whole route and the occlusion area is included. The positioning error of the occlusion area expanded. This means that the statistical errors of "CSAC+MIMU+GNSS" are larger than "CSAC+GNSS". The coupled integration positioning probability of the demo setup "CSAC+MIMU+GNSS" was 11% higher than that of "CSAC+GNSS". The horizontal positioning accuracy of "CSAC+MIMU+GNSS" was 24% higher than that of "MIMU+GNSS".  Displacements and errors from the harsh route are shown in Figure 14, and navigation results are compared in Table 7. In Table 7 the statistical range in time domain of the "CSAC+GNSS" results was 89% of the whole route and the occlusion area is not included. The positioning probability of CSAC+GNSS was 89%. However, the statistical range in time domain of the "CSAC+MIMU+GNSS" results was 100% of the whole route and the occlusion area is included. The positioning error of the occlusion area expanded. This means that the statistical errors of "CSAC+MIMU+GNSS" are larger than "CSAC+GNSS". The coupled integration positioning probability of the demo setup "CSAC+MIMU+GNSS" was 11% higher than that of "CSAC+GNSS". The horizontal positioning accuracy of "CSAC+MIMU+GNSS" was 24% higher than that of "MIMU+GNSS". Displacements and errors from the harsh route are shown in Figure 14, and navigation results are compared in Table 7. In Table 7 the statistical range in time domain of the "CSAC+GNSS" results was 89% of the whole route and the occlusion area is not included. The positioning probability of CSAC+GNSS was 89%. However, the statistical range in time domain of the "CSAC+MIMU+GNSS" results was 100% of the whole route and the occlusion area is included. The positioning error of the occlusion area expanded. This means that the statistical errors of "CSAC+MIMU+GNSS" are larger than "CSAC+GNSS". The coupled integration positioning probability of the demo setup "CSAC+MIMU+GNSS" was 11% higher than that of "CSAC+GNSS". The horizontal positioning accuracy of "CSAC+MIMU+GNSS" was 24% higher than that of "MIMU+GNSS".

Conclusions
This paper discussed the coupled integration of CSAC, MIMU, and GNSS. The architecture of the integration was designed and the mathematical models of precise time aiding navigation filter were derived. If the visible satellite number was four or larger, then CSAC aided GNSS positioning with weighted linear optimization method and integrated with MIMU. If the visible satellite number was less than four, then CSAC converted GNSS observations to ranges by clock coasting and constrained the divergence of the MIMU inertial system. Field vehicle experiments were conducted for both open-sky and harsh routes. Results showed that the coupled integration was more accurate than the traditional techniques, including the tightly coupled GNSS/MIMU. Therefore, the coupled integration of CSAC, MIMU, and GNSS can improve PNT performance.

Acknowledgments:
The authors would like to thank Liu Wei from Beijing Xingyuan Beidou Navigation Technology Co., Ltd. for providing the GNSS test tools. We would also like to thank the reviewers for their beneficial comments and suggestions.
Author Contributions: In this work, the concept was developed by Lin Ma and Zheng You, and the test was developed by Lin Ma and Tianyi Liu. Moreover, Lin Ma and Shuai Shi prepared the final draft and Zheng You guaranteed the critical reading.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
The state matrix of the state equation can be written as:

Conclusions
This paper discussed the coupled integration of CSAC, MIMU, and GNSS. The architecture of the integration was designed and the mathematical models of precise time aiding navigation filter were derived. If the visible satellite number was four or larger, then CSAC aided GNSS positioning with weighted linear optimization method and integrated with MIMU. If the visible satellite number was less than four, then CSAC converted GNSS observations to ranges by clock coasting and constrained the divergence of the MIMU inertial system. Field vehicle experiments were conducted for both open-sky and harsh routes. Results showed that the coupled integration was more accurate than the traditional techniques, including the tightly coupled GNSS/MIMU. Therefore, the coupled integration of CSAC, MIMU, and GNSS can improve PNT performance.
where A ij is a 3ˆ3 matrix, C n b is the transform matrix from carrier coordinate system Ox b y b z b to the navigation coordinate system Ox n y n z n in Figure A1. ϕ, θ, and γ are head, pitch, and roll angles, respectively.  f E , f N , f U are respectively the accelerations of east, north, and up directions in the navigation coordinate, τ δax , τ δay , τ δ az are the negative correlation coefficients of the three-axis accelerometer bias, τ δgx , τ δgy , τ δgz are the negative correlation coefficients of the three-axis gyroscope bias, τ δSFax , τ δSFay , τ δSFaz are the negative correlation coefficients of the three-axis accelerometer scale factor error, τ δSFgx , τ δSFgy , τ δSFgz are the negative correlation coefficients of the three-axis gyroscope scale factor error, and R is the Earth's radius.