Low-Earth-Orbit Satellites and Robust Theory-Augmented GPS/Inertial-Navigation-System Tight Integration for Vehicle-Borne Positioning

Abstract

Positioning by means of the Global Positioning System (GPS) is a traditional and widely used method. However, its performance is affected by the user environment, such as multi-path effects and poor anti-interference abilities. Therefore, an Inertial Navigation System (INS) has been integrated with GPS to overcome the disadvantages of GPS positioning. INSs do not rely on any external system information and has strong autonomy and independence from the external environment. However, the performance of GPS/INS is visibly degraded in low-observability GPS environments (tall buildings, viaducts, underground tunnels, woods, etc.). Fortunately, with the emergence of Low-Earth-Orbit (LEO) satellites in recent years, the constellation configuration can be extended with the advantages of lower orbits, greater speeds, and richer geometric structures. LEO improves the geometric structure between users and satellites and provides many more observations. Meanwhile, a robust theory approach is applied that can restrain or remove the impact of low-accuracy observations. In this study, we applied LEO data and a robust theory approach to enhance the GPS/INS tight integration. To verify the effectiveness of this method, a set of vehicles and simulated LEO data were analyzed. The results show that robust Kalman filtering (RKF) provides a visible enhancement in the positioning accuracy of GPS/INS integration. This effectively restrains the mutation error and has a smoothing effect on the positioning results. In addition, the addition of LEO data significantly improves the positioning accuracy of a sole GPS and GPS/INS integration. The GPS/LEO/INS integration has the highest positioning accuracy, with Root-Mean-Square Errors (RMSEs) of the north, east, and vertical positions of 2.38 m, 1.94 m, and 2.49 m, respectively, which corresponds to an improvement of 30.21%, 47.43%, and 34.13% compared to sole GPS-based positioning and 8.60%, 17.24%, and 12.14% when compared to the GPS/INS mode. Simultaneously, the simulation results show that LEO and INSs can improve the positioning performance of GPS under GPS-blocked conditions.

1. Introduction

The Global Positioning System (GPS) is a commonly used positioning system that provides users with a high-precision three-dimensional position and velocity [1], owing to its all-weather, real-time, global coverage, and uninterrupted service characteristics. Owing to these advantages, GPS is widely adopted in high-precision measurements, mobile monitoring, vehicle navigation, and search and rescue fields [2,3,4]. In general, GPS positioning technologies can be classified into Single Point Positioning (SPP), Precise Point Positioning (PPP), Real-Time Kinematic (RTK) [5], and PPP-RTK [6,7]. SPP uses pseudo-range observations of GPS and the broadcast ephemeris for positioning and can achieve a few-meter-level positioning accuracy. However, with the user demand for higher positioning accuracy, RTK [8] and PPP [9,10,11] have been proposed. RTK technology eliminates common errors such as receiver and satellite clocks and weakens the orbit error, tropospheric delay, ionospheric delay, and other distance-related errors through the double-difference observation model, to achieve fast centimeter-level positioning [12]. Although it has high accuracy and a short convergence time, the positioning accuracy gradually decreases with an increase in baseline length, which leads to certain limitations in the application range. PPP technology utilizes the precise ephemeris and satellite clock provided by the International GNSS Service (IGS), using only the pseudo-range and carrier phase measurements from a single receiver, to achieve high-precision absolute positioning on a global scale. PPP requires only a single receiver, and its positioning is not limited by the baseline, which makes it more flexible and less costly. However, PPP requires a long initialization time and continuous satellite signals to maintain high accuracy [13].

However, in complex and harsh environments (tall buildings, viaducts, underground tunnels, woods, etc.), the GPS signal, as a form of radio navigation and positioning, is prone to interference and frequent interruptions, resulting in low-quality observation data, insufficiently visible satellites, and even no positioning solutions [14]. For GPS positioning drawbacks in complex and harsh environments, the integration of other technologies with GPS can effectively solve GPS defects.

In light of these drawbacks, GPS is often integrated with the Inertial Navigation System (INS) [15,16,17,18,19,20]. The INS is a completely autonomous navigation system that does not depend on the external environment. It can provide a variety of high-precision navigation parameters independently and autonomously, owing to its autonomy, concealment capability, high sampling rate, and short-term accuracy [21,22]. As a result, INSs can provide navigation information even if GPS is affected. The integration of INSs and GPS [23,24,25] can effectively overcome the shortcomings of GPS and INSs and utilize the advantages of the respective systems to provide users with more reliable and high-precision positioning. However, the errors of INSs accumulate with time, especially for low-cost Micro-Electro-Mechanical Systems (MEMSs) [26]. Several studies have focused on GPS and INS integration. In 2006, Le et al. used 1 h of onboard experimental data to analyze the PPP/INS loose integration and obtained 0.5 m position results in the horizontal direction and 0.6 m in the vertical direction [27]. In 2009, Martell et al. realized PPP/INS tight integration and the results of airborne experiments showed that the positioning accuracy was within 15 cm [28]. In 2017, Zhang et al. proposed a smoothing algorithm for GNSS/SINS tight integration, and experiments showed that the smoothing algorithm improved the average three-dimensional (3D) position Root-Mean-Square Error (RMSE) and the 3D attitude RMSE on average by 65.7% and 70%, respectively [29]. GPS and INS integration has been widely used. However, the integrated system is overly dependent on GPS, which can drastically reduce the accuracy of GPS/INS integration, even though it can provide positioning when the GPS signal is weak.

Due to the strong dependence, LEO data can be introduced to provide many observations. In recent years, with the proposal of a comprehensive PNT (Positioning, Navigation, and Timing), LEO has attracted attention from navigation and communication fields worldwide owing to its advantages, such as rapid speed, strong signal power, and a large number of satellites [30,31,32]. Owing to their unique advantages, LEO satellites can effectively complement and improve GPS services. The tight integration of LEO with GPS and INSs can enhance their accuracy. Ke et al. conducted a preliminary evaluation of LEO-augmented GPS PPP using simulation data. The results show that the convergence time of the GPS/LEO-integrated PPP-AR was reduced by 51.31% and the positioning accuracy was improved by 14.9% compared with that of GPS [33]. Ge et al. simulated LEO satellites with GPS/BDS observations for LEO-augmented GPS/BDS PPP. Their results show that the PPP convergence time can be decreased down to approximately 5 min in most of the global regions. Compared with the 30 s sampling interval, the 1 s sampling interval could further reduce the convergence time [34]. Li et al. analyzed the contribution of LEO to improving the fixed ambiguity of PPP and found that the positioning errors in the east, north, and vertical directions after adding LEO observations improved by 63.2%, 67.2%, and 20.7%, respectively, compared with PPP floating-point solutions [35]. Li et al. found that the greater the number of LEO satellites in LEO-enhanced multi-GNSS PPP, the shorter the convergence time and the higher the positioning accuracy. To achieve a horizontal accuracy of <10 cm, they introduced 60, 96, 192, and 288 LEO satellites, and the convergence time was shortened from 9.6 min to 7.0, 3.2, 2.1, and 1.3 min, respectively [36]. Zhao et al. evaluated the PPP performance of LEO/GNSS integration in harsh environments and showed that the PPP convergence time for LEO/GNSS integration was significantly improved. Compared with a single BDS, the float solution accuracy of the BDS/LEO and BDS/GPS/LEO-integrated PPP improved by 73.77% and 77.43%, respectively [37]. Many scholars have investigated the enhancement in GPS positioning by using simulated LEO data. However, there are few studies on the effect of integrated GPS/INS positioning.

In this study, an LEO-enhanced GPS/INS tight integration algorithm was modeled. Based on the previous GPS/INS tight integration model, LEO data were added to realize a tight integration of GPS/LEO/INS. The remainder of this paper is organized as follows. The mathematical model of the tight integration algorithm is presented in Section 2, where the observation equations, state equations, robust Kalman filtering model, and algorithmic framework are given in detail. The experimental scheme and analysis are presented in Section 3, where the tight integration algorithm is experimentally validated using on-board data and simulated LEO data. At the same time, the applicability of the tight integration algorithm in a weak observation environment is simulated. A discussion of the experimental validation results is presented in Section 4. Finally, the conclusions are presented in Section 5.

2. GPS/LEO/INS Mathematical Models

GPS/LEO/INS tight integration uses a dual-frequency pseudo-range and Doppler to form an innovation vector in the Kalman filter. The carrier’s position, velocity, attitude, bias of gyroscopes and accelerometers, scale factor of gyroscopes and accelerometers, and receiver clock offset and clock drift of GPS and LEO were selected as the estimated parameters. In the following equations, the GPS, INS, and LEO observations are denoted as G, I, and L, respectively.

The Kalman filter state equation for the tight integration of GPS/LEO/INS is expressed as [38,39]

$$\dot{\mathit{X}}\left(\mathit{t}\right)=\mathit{F}\left(\mathit{t}\right)\mathit{X}\left(\mathit{t}\right)+\mathit{G}\left(\mathit{t}\right)\mathit{w}\left(\mathit{t}\right)$$

(1)

where $\mathit{X}\left(\mathit{t}\right)$ represents the system state vector, $\mathit{F}\left(\mathit{t}\right)$ represents the system state transition matrix, $\mathit{G}\left(\mathit{t}\right)$ represents the system noise driving matrix, and $\mathit{w}\left(\mathit{t}\right)$ represents the system noise vector that follows a Gaussian normal distribution with a mean value of zero.

The observation innovation equation for the tight integration of GPS/LEO/INS is expressed as

$${\mathit{Z}}_{\mathit{k}}={\mathit{H}}_{\mathit{k}}{\mathit{X}}_{\mathit{k}}+{\mathit{v}}_{\mathit{k}}$$

(2)

where ${\mathit{H}}_{\mathit{k}}$ represents the design coefficient matrix of the k epoch, and ${\mathit{v}}_{\mathit{k}}$ represents the observed noise, which follows a Gaussian normal distribution with a mean value of zero.

2.1. Algorithm Structure

The GPS/LEO/INS tight integration algorithm is mainly composed of a GPS/INS model, LEO/INS model, and GPS/LEO/INS model, which is illustrated in Figure 1. The raw observations (pseudo-range observations or Doppler observations) from GPS and LEO were used to integrate with the pseudo-range or pseudo-range rate predicted by the INS.

2.2. Observation Function of GPS/LEO/INS Tight Integration

The observation innovation vector of GPS/LEO/INS tight integration was obtained by subtracting the dual-frequency pseudo-range and Doppler observations of GPS and LEO ($\mathit{O}$) from the pseudo-range and Doppler values predicted by the INS (${\tilde{\mathit{O}}}_{\mathit{I}\mathit{N}\mathit{S}}$):

$${\mathit{Z}}_{\mathit{k}}=\mathit{O}-{\tilde{\mathit{O}}}_{\mathit{I}\mathit{N}\mathit{S}}=\left[\begin{array}{c}{\mathit{Z}}_{{\mathit{P}}_{\mathit{P}\mathit{C}}}\\ {\mathit{Z}}_{{\dot{\mathit{P}}}_{\mathit{D}\mathit{C}}}\end{array}\right]$$

(3)

To obtain INS-predicted values, INS mechanization was required to update the receiver position and velocity. However, owing to the different reference centers of the INS and GPS/INS antenna, a lever-arm correction was required. Accordingly, the linearized function ${\mathit{Z}}_{\mathit{k}}$ after considering the lever-arm correction can be expressed as [22,40]

$$\mathit{\delta }{\mathit{Z}}_{{\mathit{P}}_{\mathit{PC}}^{\mathit{G}}}={\mathit{C}}_1\left(\mathit{\delta }{\mathit{p}}_{\mathit{INS}}^{\mathit{n}}+\left({\mathit{C}}_{\mathit{b}}^{\mathit{n}}{\mathit{\iota }}_{\mathit{G}}^{\mathit{b}}\times \right)\mathit{\delta }\mathbf{\psi }\right)+\mathit{c}\mathit{\delta }{\mathit{t}}_{\mathit{G}}$$

(4)

$$\mathit{\delta }{\mathit{Z}}_{{\dot{\mathit{P}}}_{\mathit{DC}}^{\mathit{G}}}={\mathit{C}}_2{\mathit{D}}^{-1}\mathit{\delta }{\mathit{p}}_{\mathit{INS}}^{\mathit{n}}+{\mathit{C}}_{\mathit{n}}^{\mathit{e}}{\mathit{H}}_{\mathit{\Psi }}\mathit{\delta }\mathit{\Psi }+{\mathit{C}}_{\mathit{n}}^{\mathit{e}}\mathit{\delta }{\mathit{v}}_{\mathit{INS}}^{\mathit{n}}+{\mathit{C}}_{\mathit{n}}^{\mathit{e}}{\mathit{C}}_{\mathit{b}}^{\mathit{n}}\left({\mathit{\iota }}_{\mathit{G}}^{\mathit{b}}\right)\mathit{\delta }{\mathit{w}}_{\mathit{ib}}^{\mathit{b}}+\mathit{c}\mathit{\delta }{\dot{\mathit{t}}}_{\mathit{G}}$$

(5)

$$\mathit{\delta }{\mathit{Z}}_{{\mathit{P}}_{\mathit{PC}}^{\mathit{L}}}={\mathit{C}}_1\left(\mathit{\delta }{\mathit{p}}_{\mathit{INS}}^{\mathit{n}}+\left({\mathit{C}}_{\mathit{b}}^{\mathit{n}}{\mathit{\iota }}_{\mathit{L}}^{\mathit{b}}\times \right)\mathit{\delta }\mathbf{\psi }\right)+\mathit{c}\mathit{\delta }{\mathit{t}}_{\mathit{L}}$$

(6)

$$\mathit{\delta }{\mathit{Z}}_{{\dot{\mathit{P}}}_{\mathit{DC}}^{\mathit{L}}}={\mathit{C}}_2{\mathit{D}}^{-1}\mathit{\delta }{\mathit{p}}_{\mathit{INS}}^{\mathit{n}}+{\mathit{C}}_{\mathit{n}}^{\mathit{e}}{\mathit{H}}_{\mathit{\Psi }}\mathit{\delta }\mathit{\Psi }+{\mathit{C}}_{\mathit{n}}^{\mathit{e}}\mathit{\delta }{\mathit{v}}_{\mathit{INS}}^{\mathit{n}}+{\mathit{C}}_{\mathit{n}}^{\mathit{e}}{\mathit{C}}_{\mathit{b}}^{\mathit{n}}\left({\mathit{\iota }}_{\mathit{L}}^{\mathit{b}}\right)\mathit{\delta }{\mathit{w}}_{\mathit{ib}}^{\mathit{b}}+\mathit{c}\mathit{\delta }{\dot{\mathit{t}}}_{\mathit{L}}$$

(7)

where

$${\mathit{D}}^{-1}=diag\left(\frac{1}{R_M+h},\frac{1}{R_N+h},-1\right)$$

(8)

$${\mathit{C}}_1=\left[\begin{array}{ccc}\frac{-\left(R_N+h\right)cosLsinB}{R_M+h}& -sinL& -cosBcosL\\ \frac{-\left(R_N+h\right)sinBsinL}{R_M+h}& cosL& -cosBcosL\\ \frac{\left[R_N\left(1-{\lambda }^2\right)\right]cosBcosL}{R_M+h}& 0& -sinB\end{array}\right]$$

(9)

$${\mathit{C}}_2=\left({\mathit{C}}_{\mathit{b}}^{\mathit{n}}{\mathit{\iota }}_{\mathit{G}}^{\mathit{b}}\right)\times$$

(10)

where n, b, i, e, and p represent the navigation frame, body frame, inertial frame, Earth-centered Earth-fixed coordinate system, and platform frame, respectively; $\mathit{\delta }{\mathit{p}}_{\mathit{I}\mathit{N}\mathit{S}}^{\mathit{n}}$, ${\mathit{v}}_{\mathit{INS}}^{\mathit{n}}$, and $\mathit{\delta }\mathit{\Psi }$ represent the position, velocity, and attitude correction in the n-frame, respectively; ${\mathit{\iota }}_{\mathit{G}}^{\mathit{b}}$ represents the value of the lever-arm in the b-frame for GPS; ${\mathit{C}}_{\mathit{k}}^{\mathit{j}}$ (k = n and b, j = e and n) represents the rotation matrix from the j-system to the k-system; $\mathit{\delta }{\mathit{t}}_{\mathit{G}}$ and $\mathit{\delta }{\dot{\mathit{t}}}_{\mathit{G}}$ represent the receiver clock offset and clock drift of GPS, respectively; ${\mathit{w}}_{\mathit{ib}}^{\mathit{b}}$ represents the gyroscope error; ${\mathit{C}}_1$ represents the transformation matrix from e-frame coordinates to the n-frame; ${\mathit{C}}_2$ represents the transformation matrix related to the compensation of the lever-arm error; $R_M$ and $R_N$ represent the radius of curvature of the meridian circle and the radius of curvature in the prime vertical where the carrier is located, respectively; $\mathit{\lambda }$ represents the eccentricity of the meridian ellipse; and $B$, $L$, and $h$ are geodetic latitude, longitude, and height, respectively.

By combining Equations (4)–(7), the observation innovation vector for the GPS/LEO/INS tight integration can be expressed as

$${\mathit{Z}}_{\mathit{k}}=\left[\begin{array}{c}{\mathit{Z}}_{{\mathit{P}}_{\mathit{PC}}^{\mathit{G}}}\\ {\mathit{Z}}_{{\dot{\mathit{P}}}_{\mathit{DC}}^{\mathit{G}}}\\ {\mathit{Z}}_{{\mathit{P}}_{\mathit{PC}}^{\mathit{L}}}\\ {\mathit{Z}}_{{\dot{\mathit{P}}}_{\mathit{DC}}^{\mathit{L}}}\end{array}\right]=\left[\begin{array}{c}\mathit{\delta }{\mathit{Z}}_{{\mathit{P}}_{\mathit{PC}}^{\mathit{G}}}\\ \mathit{\delta }{\mathit{Z}}_{{\dot{\mathit{P}}}_{\mathit{DC}}^{\mathit{G}}}\\ \mathit{\delta }{\mathit{Z}}_{{\mathit{P}}_{\mathit{PC}}^{\mathit{L}}}\\ \mathit{\delta }{\mathit{Z}}_{{\dot{\mathit{P}}}_{\mathit{DC}}^{\mathit{L}}}\end{array}\right]$$

(11)

The design coefficient matrix of the GPS/LEO/INS tight integration can be expressed as follows:

$${\mathit{H}}_{\mathit{k}}=\left[\begin{array}{cc}{\mathit{H}}_{\mathit{k}}^{\mathit{G}}& {\mathit{H}}_{\mathit{k}}^{\mathit{L}}\end{array}\right]$$

(12)

$${\mathit{H}}_{\mathit{k}}^{\mathit{G}}=\left[\begin{array}{c}{\mathit{H}}_{{\mathit{P}}_{\mathit{PC}}}^{\mathit{G}}\\ {\mathit{H}}_{{\dot{\mathit{P}}}_{\mathit{DC}}}^{\mathit{G}}\end{array}\right]=\left[\begin{array}{ccccccccc}{\mathit{H}}_1^{\mathit{G}}& 0& {\mathit{H}}_2^{\mathit{G}}& 0& 0& 0& 0& {\mathit{H}}_3^{\mathit{G}}& 0\\ {\mathit{H}}_4^{\mathit{G}}& {\mathit{H}}_5^{\mathit{G}}& {\mathit{H}}_6^{\mathit{G}}& 0& {\mathit{H}}_7^{\mathit{G}}& 0& {\mathit{H}}_8^{\mathit{G}}& 0& {\mathit{H}}_3^{\mathit{G}}\end{array}\right]$$

(13)

$${\mathit{H}}_{\mathit{k}}^{\mathit{L}}=\left[\begin{array}{c}{\mathit{H}}_{{\mathit{P}}_{\mathit{PC}}}^{\mathit{L}}\\ {\mathit{H}}_{{\dot{\mathit{P}}}_{\mathit{DC}}}^{\mathit{L}}\end{array}\right]=\left[\begin{array}{ccccccccc}{\mathit{H}}_1^{\mathit{L}}& 0& {\mathit{H}}_2^{\mathit{L}}& 0& 0& 0& 0& {\mathit{H}}_3^{\mathit{L}}& 0\\ {\mathit{H}}_4^{\mathit{L}}& {\mathit{H}}_5^{\mathit{L}}& {\mathit{H}}_6^{\mathit{L}}& 0& {\mathit{H}}_7^{\mathit{L}}& 0& {\mathit{H}}_8^{\mathit{L}}& 0& {\mathit{H}}_3^{\mathit{L}}\end{array}\right]$$

(14)

By combining Equations (13) and (14), the final design coefficient matrix of the GPS/LEO/INS tight integration can be expressed as

$${\mathit{H}}_{\mathit{k}}=\left[\begin{array}{c}{\mathit{H}}_{{\mathit{P}}_{\mathit{PC}}}^{\mathit{G}}\\ {\mathit{H}}_{{\dot{\mathit{P}}}_{\mathit{DC}}}^{\mathit{G}}\\ {\mathit{H}}_{{\mathit{P}}_{\mathit{PC}}}^{\mathit{L}}\\ {\mathit{H}}_{{\dot{\mathit{P}}}_{\mathit{DC}}}^{\mathit{L}}\end{array}\right]=\left[\begin{array}{ccccccccccc}{\mathit{H}}_1^{\mathit{G}}& 0& {\mathit{H}}_2^{\mathit{G}}& 0& 0& 0& 0& {\mathit{H}}_3^{\mathit{G}}& 0& 0& 0\\ {\mathit{H}}_4^{\mathit{G}}& {\mathit{H}}_5^{\mathit{G}}& {\mathit{H}}_6^{\mathit{G}}& 0& {\mathit{H}}_7^{\mathit{G}}& 0& {\mathit{H}}_8^{\mathit{G}}& 0& {\mathit{H}}_3^{\mathit{G}}& 0& 0\\ {\mathit{H}}_1^{\mathit{L}}& 0& {\mathit{H}}_2^{\mathit{L}}& 0& 0& 0& 0& 0& 0& {\mathit{H}}_3^{\mathit{L}}& 0\\ {\mathit{H}}_4^{\mathit{L}}& {\mathit{H}}_5^{\mathit{L}}& {\mathit{H}}_6^{\mathit{L}}& 0& {\mathit{H}}_7^{\mathit{L}}& 0& {\mathit{H}}_8^{\mathit{L}}& 0& 0& 0& {\mathit{H}}_3^{\mathit{L}}\end{array}\right]$$

(15)

with

$$\left\{\begin{array}{c}{\mathit{H}}_1=\mathit{A}{\mathit{C}}_1\\ {\mathit{H}}_2={\mathit{H}}_1\left({\mathit{C}}_{\mathit{b}}^{\mathit{n}}{\mathit{\iota }}^{\mathit{b}}\times \right)\\ {\mathit{H}}_3={\left[\begin{array}{ccc}1& \cdots & 1\end{array}\right]}^{\mathit{T}}\\ {\mathit{H}}_4=\mathit{A}{\mathit{D}}^{-1}{\mathit{C}}_2\\ {\mathit{H}}_5=\mathit{A}{\mathit{C}}_{\mathit{n}}^{\mathit{e}}\\ {\mathit{H}}_6={\mathit{H}}_4\left({\mathit{C}}_{\mathit{b}}^{\mathit{n}}{\mathit{\iota }}^{\mathit{b}}\times \right)-{\mathit{H}}_5\left[\left({\mathit{w}}_{\mathit{en}}^{\mathit{n}}\times +{\mathit{w}}_{\mathit{ie}}^{\mathit{n}}\times \right){\mathit{C}}_{\mathit{b}}^{\mathit{n}}\left({\mathit{\iota }}^{\mathit{b}}\times \right)+{\mathit{C}}_{\mathit{b}}^{\mathit{n}}\left({\mathit{\iota }}^{\mathit{b}}\times {\mathit{w}}_{\mathit{ib}}^{\mathit{b}}\right)\times \right]\\ {\mathit{H}}_7=-{\mathit{H}}_5{\mathit{C}}_{\mathit{b}}^{\mathit{n}}\left({\mathit{\iota }}^{\mathit{b}}\times \right)\\ {\mathit{H}}_8={\mathit{H}}_7\mathit{diag}\left({\mathit{w}}_{\mathit{ib}}^{\mathit{b}}\right)\end{array}\right.$$

(16)

where $\mathit{A}$ covers the direction cosine of the satellite–receiver vector; ${\mathit{C}}_{\mathit{b}}^{\mathit{n}}$ and ${\mathit{C}}_{\mathit{n}}^{\mathit{e}}$ represent the conversion matrix from b-frame to n-frame and the conversion matrix from n-frame to e-frame, respectively; ${\mathit{w}}_{\mathit{ib}}^{\mathit{b}}$ represents the angular rate of the gyroscope output of the b-frame to the i-frame, sensed in the b-frame; ${\mathit{w}}_{\mathit{en}}^{\mathit{n}}$ and ${\mathit{w}}_{\mathit{ie}}^{\mathit{n}}$ represent the angular rate of the n-frame to the e-frame and the e-frame to i-frame, sensed in the n-frame, respectively.

2.3. State Function

The state vector of INS is expressed as

$${\mathit{X}}_1\left(\mathit{t}\right)=\left[\mathit{\delta }{\mathit{p}}_{\mathit{I}}\mathit{\delta }{\mathit{v}}_{\mathit{I}}\mathit{\delta }\mathit{\Psi }\mathit{\delta }{\mathit{B}}_{\mathit{a}}\mathit{\delta }{\mathit{B}}_{\mathit{g}}\mathit{\delta }{\mathit{S}}_{\mathit{a}}\mathit{\delta }{\mathit{S}}_{\mathit{g}}\right]$$

(17)

where $\mathit{\delta }{\mathit{p}}_{\mathit{I}}$ represents position correction; $\mathit{\delta }{\mathit{v}}_{\mathit{I}}$ represents velocity correction; $\mathit{\delta }\mathit{\Psi }$ represents attitude correction; $\mathit{\delta }{\mathit{B}}_{\mathit{a}}$ and $\mathit{\delta }{\mathit{B}}_{\mathit{g}}$ represent the bias corrections for the accelerometer and gyroscope, respectively; and $\mathit{\delta }{\mathit{S}}_{\mathit{a}}$ and $\mathit{\delta }{\mathit{S}}_{\mathit{g}}$ represent the scale factor corrections for the accelerometer and gyroscope, respectively.

The state transfer matrix of the INS is expressed as

$${\mathit{F}}_{\mathit{I}}\left(\mathit{t}\right)=\left[\begin{array}{ccccccc}{\mathit{F}}_1& {\mathit{F}}_2& 0& 0& 0& 0& 0\\ {\mathit{F}}_3& {\mathit{F}}_4& {\mathit{F}}_5& {\mathit{F}}_6& 0& {\mathit{F}}_7& 0\\ 0& 0& {\mathit{F}}_8& 0& {\mathit{F}}_9& 0& {\mathit{F}}_{10}\\ 0& 0& 0& {\mathit{F}}_{11}& 0& 0& 0\\ 0& 0& 0& 0& {\mathit{F}}_{12}& 0& 0\\ 0& 0& 0& 0& 0& {\mathit{F}}_{13}& 0\\ 0& 0& 0& 0& 0& 0& {\mathit{F}}_{14}\end{array}\right]$$

(18)

with

$$\left\{\begin{array}{c}\begin{array}{c}{\mathit{F}}_1={\left[\mathit{I}-\left({\mathit{w}}_{\mathit{en}}^{\mathit{n}}\times \right)\mathbf{\Delta }\mathit{t}\right]}_{3\times 3}\\ {\mathit{F}}_2={\left[\mathit{I}\mathbf{\Delta }\mathit{t}\right]}_{3\times 3}\\ {\mathit{F}}_3=\mathit{diag}{\left(\frac{-g\Delta t}{R_M+h},\frac{-g\Delta t}{R_N+h},\frac{-g\Delta t}{2\left(\sqrt{R_N{R}_M}+h\right)}\right)}_{3\times 3}\\ {\mathit{F}}_4={\left[\mathit{I}-\left(\left({\mathit{w}}_{\mathit{ie}}^{\mathit{n}}+{\mathit{w}}_{\mathit{in}}^{\mathit{n}}\right)\mathbf{\Delta }\mathit{t}\right)\times \right]}_{3\times 3}\\ {\mathit{F}}_5={\left[\left({\mathit{f}}^{\mathit{n}}\times \right)\mathbf{\Delta }\mathit{t}\right]}_{3\times 3}\\ {\mathit{F}}_6={\left[{\mathit{C}}_{\mathit{b}}^{\mathit{p}}\mathbf{\Delta }\mathit{t}\right]}_{3\times 3}\\ {\mathit{F}}_7={\left[{\mathit{C}}_{\mathit{b}}^{\mathit{p}}{\mathit{f}}^{\mathit{b}}\mathbf{\Delta }\mathit{t}\right]}_{3\times 3}\end{array}\\ \begin{array}{l}\hspace{1em}{\mathit{F}}_8={\left[\mathit{I}-\left({\mathit{w}}_{\mathit{in}}^{\mathit{n}}\times \right)\mathbf{\Delta }\mathit{t}\right]}_{3\times 3}\\ \hspace{1em}{\mathit{F}}_9={\left[-{\mathit{C}}_{\mathit{b}}^{\mathit{n}}\mathbf{\Delta }\mathit{t}\right]}_{3\times 3}\\ \hspace{1em}{\mathit{F}}_{10}={\left[{\mathit{C}}_{\mathit{b}}^{\mathit{p}}{\mathit{w}}_{\mathit{ib}}^{\mathit{b}}\mathbf{\Delta }\mathit{t}\right]}_{3\times 3}\\ \hspace{1em}{\mathit{F}}_{11}=\mathit{diag}{\left(e^{-\Delta t/T_{B_a}},e^{-\Delta t/T_{B_a}},e^{-\Delta t/T_{B_a}}\right)}_{3\times 3}\\ \hspace{1em}{\mathit{F}}_{12}=\mathit{diag}{\left(e^{-\Delta t/T_{S_a}},e^{-\Delta t/T_{S_a}},e^{-\Delta t/T_{S_a}}\right)}_{3\times 3}\\ \hspace{1em}{\mathit{F}}_{13}=\mathit{diag}{\left(e^{-\Delta t/T_{B_g}},e^{-\Delta t/T_{B_g}},e^{-\Delta t/T_{B_g}}\right)}_{3\times 3}\\ \hspace{1em}{\mathit{F}}_{14}=\mathit{diag}{\left(e^{-\Delta t/T_{S_g}},e^{-\Delta t/T_{S_g}},e^{-\Delta t/T_{S_g}}\right)}_{3\times 3}\end{array}\end{array}\right.$$

(19)

The GPS receiver clock usually introduces a clock error which is translated into a ranging error. The state vector of the GPS receiver clock is expressed as

$${\mathit{X}}_G\left(\mathit{t}\right)={\left[\delta t_G\hspace{1em}\delta {\dot{t}}_G\right]}^T$$

(20)

where $\delta t_G$ and $\delta {\dot{t}}_G$ represent the GPS receiver clock offset and clock drift, respectively.

Similarly, the state transfer matrix of the GPS receiver clock is expressed as

$${\mathit{F}}_G\left(\mathit{t}\right)=\left[\begin{array}{cc}1& \Delta t\\ 0& 1\end{array}\right]$$

(21)

Similarly, the state vector of the LEO receiver clock is expressed as

$${\mathit{X}}_L\left(\mathit{t}\right)={\left[\mathit{\delta }{t}_L\hspace{1em}\mathit{\delta }{\dot{t}}_L\right]}^T$$

(22)

where $\mathit{\delta }{t}_L$ and $\mathit{\delta }{\dot{t}}_L$ represent the LEO receiver clock offset and clock drift, respectively.

The state transfer matrix of the LEO receiver clock is expressed as

$${\mathit{F}}_{\mathit{L}}\left(\mathit{t}\right)=\left[\begin{array}{cc}1& \Delta t\\ 0& 1\end{array}\right]$$

(23)

Therefore, by combining Equations (17), (20), and (22), we obtain a state function for tight integration of the GPS/LEO/INS:

$$\left[\begin{array}{c}{\dot{\mathit{X}}}_{\mathit{I}}(\mathit{t})\\ {\dot{\mathit{X}}}_{\mathit{G}}(\mathit{t})\\ {\dot{\mathit{X}}}_{\mathit{L}}(\mathit{t})\end{array}\right]=\left[\begin{array}{ccc}{\mathit{F}}_{\mathit{I}}(\mathit{t})& 0& 0\\ 0& {\mathit{F}}_{\mathit{G}}(\mathit{t})& 0\\ 0& 0& {\mathit{F}}_{\mathit{L}}(\mathit{t})\end{array}\right]\left[\begin{array}{c}{\mathit{X}}_{\mathit{I}}(\mathit{t})\\ {\mathit{X}}_{\mathit{G}}(\mathit{t})\\ {\mathit{X}}_{\mathit{L}}(\mathit{t})\end{array}\right]+\left[\begin{array}{ccc}{\mathit{G}}_{\mathit{I}}(\mathit{t})& 0& 0\\ 0& {\mathit{G}}_{\mathit{G}}(\mathit{t})& 0\\ 0& 0& {\mathit{G}}_{\mathit{L}}(\mathit{t})\end{array}\right]\left[\begin{array}{c}{\mathit{w}}_{\mathit{I}}(\mathit{t})\\ {\mathit{w}}_{\mathit{G}}(\mathit{t})\\ {\mathit{w}}_{\mathit{L}}(\mathit{t})\end{array}\right]$$

(24)

with

$$\mathit{X}\left(\mathit{t}\right)={\left[{\mathit{X}}_{\mathit{I}}(\mathit{t}){\mathit{X}}_{\mathit{G}}(\mathit{t}){\mathit{X}}_{\mathit{L}}(\mathit{t})\right]}^{\mathit{T}}=\left[\mathit{\delta }{\mathit{p}}_{\mathit{I}}\mathit{\delta }{\mathit{v}}_{\mathit{I}}\mathit{\delta }\mathit{\Psi }\mathit{\delta }{\mathit{B}}_{\mathit{a}}\mathit{\delta }{\mathit{B}}_{\mathit{g}}\mathit{\delta }{\mathit{S}}_{\mathit{a}}\mathit{\delta }{\mathit{S}}_{\mathit{g}}\mathit{\delta }{t}_G\mathit{\delta }{\dot{t}}_G\mathit{\delta }{t}_L\mathit{\delta }{\dot{t}}_L\right]$$

(25)

In general, the states of the clock offset and clock drift of GPS and LEO can be modeled as random walk process by means of the dynamic models:

$$\mathit{\delta }{t}_k=\mathit{\delta }{t}_{k-1}+\mathit{\delta }{\dot{t}}_{k-1}\Delta t+v_{k-1}$$

(26)

$$\mathit{\delta }{\dot{t}}_k=\mathit{\delta }{\dot{t}}_{k-1}+{\dot{v}}_{k-1}$$

(27)

where $\mathit{\delta }{t}_k$ and $\mathit{\delta }{\dot{\mathit{t}}}_k$ represent the receiver clock offset and clock drift, respectively; $v_{k-1}$ and ${\dot{v}}_{k-1}$ represent the driving white noise of the receiver clock offset and clock drift, respectively; and $\Delta t$ represents the interval between adjacent epochs. According to the law of variance–covariance propagation, the prior variance can be expressed as

$${\mathit{\sigma }}_{v_{k-1}}^{\mathbf{2}}=\mathbf{2}c{h}_{\mathbf{0}}\mathbf{\Delta }t$$

(28)

$${\mathit{\sigma }}_{{\dot{v}}_{k-\mathbf{1}}}^{\mathbf{2}}=\mathbf{8}c{h}_{\mathbf{2}}{\pi }^2\mathbf{\Delta }t$$

(29)

where c represents the speed of light in vacuum and the values of $h_0$ and $h_2$ are related to the model selected by the receiver crystal oscillator.

In general, accelerometer and gyroscope zero-bias errors and scale factor errors are modeled as first-order Gauss–Markov processes, expressed as

$$S_k=e^{-\Delta t/T}{S}_{k-1}+v_{S_{k-1}},v_{S_{k-1}}\sim \left(0,2{\sigma }^{2\Delta t}/T\right)$$

(30)

$$B_k=e^{-\Delta t/T}{B}_{k-1}+v_{B_{k-1}},v_{B_{k-1}}\sim \left(0,2{\sigma }^{2\Delta t}/T\right)$$

(31)

where $\Delta t$ and $T$ represent the INS interval and correlation time of bias, respectively; $v_{S_{k-1}}$ and $v_{B_{k-1}}$ represent the driving white noise of scale factor and bias, respectively. The prior variance was determined by the hardware performance of the INS.

2.4. Robust Kalman Filter

In this study, to limit the effect of low-quality observations and anomalous observations on parameter estimation, ICG-III [41] was used:

$${\gamma }_{ii}=\left\{\begin{array}{c}1\\ \frac{k_0}{\left|v_i/\mathit{\delta }\right|}{\left(\frac{k_1-\left|v_i/\mathit{\delta }\right|}{k_1-k_0}\right)}^2,\\ 0\end{array}\right.\begin{array}{c}\left|v_i/\mathit{\delta }\right|\le k_0\\ k_0\text{<}\left|v_i/\mathit{\delta }\right|<k_1\\ \left|v_i/\mathit{\delta }\right|\ge k_0\end{array}$$

(32)

where $v_i$ represents the observed residuals; $\delta$ represents the variance of the observed residuals; ${\gamma }_{ii}$ represents the resistance factors calculated by equivalent weight functions; $k_0$ and $k_1$ are constants, for which $k_0$ is usually 1.0–1.5 and $k_1$ is usually 2.5–8.0.

Based on the factor of ${\gamma }_{ii}$, ${\mathit{R}}_{\mathit{k}}$ can be expressed as

$${\mathit{R}}_{\mathit{k}}=\left[\begin{array}{ccc}{\mathit{\gamma }}_{11}{\mathit{\sigma }}_1^2& \cdots & {\mathit{\gamma }}_{1m}{\mathit{\sigma }}_{1m}\\ ⋮& \ddots & ⋮\\ {\mathit{\gamma }}_{m1}{\mathit{\sigma }}_{m1}& \cdots & {\mathit{\gamma }}_{\mathit{mm}}{\mathit{\sigma }}_m^2\end{array}\right]$$

(33)

$${\gamma }_{\mathit{ij}}=\sqrt{{\gamma }_{\mathit{ii}}{\gamma }_{\mathit{jj}}}$$

(34)

where ${\gamma }_{\mathit{ij}}$ is the covariance expansion factor.

Such resistance factors are applied in the Kalman filter update phase by

$${\mathit{K}}_{\mathit{k}}={\mathit{P}}_{\mathit{k},\mathit{k}-1}{\mathit{H}}_{\mathit{k}}^{\mathit{T}}\left({\mathit{H}}_{\mathit{k}}{\mathit{P}}_{\mathit{k},\mathit{k}-1}{\mathit{H}}_{\mathit{k}}^{\mathit{T}}+{\mathit{R}}_{\mathit{k}}\right)$$

(35)

$${\mathit{X}}_{\mathit{k}}={\mathit{X}}_{\mathit{k},\mathit{k}-1}+{\mathit{K}}_{\mathit{k}}\left({\mathit{Z}}_{\mathit{k}}-{\mathit{H}}_{\mathit{k}}{\mathit{X}}_{\mathit{k},\mathit{k}-1}\right)$$

(36)

$${\mathit{P}}_{\mathit{k}}=\left(\mathit{I}-{\mathit{K}}_{\mathit{k}}{\mathit{H}}_{\mathit{k}}\right){\mathit{P}}_{\mathit{k},\mathit{k}-1}{\left(\mathit{I}-{\mathit{K}}_{\mathit{k}}{\mathit{H}}_{\mathit{k}}\right)}^{\mathit{T}}+{\mathit{K}}_{\mathit{k}}{R}_{\mathit{k}}{\mathit{K}}_{\mathit{k}}^{\mathit{T}}$$

(37)

$${\mathit{P}}_{\mathit{k},\mathit{k}-1}={\mathit{\Phi }}_{\mathit{k},\mathit{k}-1}{\mathit{X}}_{\mathit{k}-1}{\mathit{\Phi }}_{\mathit{k},\mathit{k}-1}^{\mathit{T}}+{\mathit{Q}}_{\mathit{k}-1}$$

(38)

where ${\mathit{X}}_{\mathit{k}}$ represents the estimated state vector at $t_k$, ${\mathit{X}}_{\mathit{k}-1}$ represents the estimate at $t_{k-1}$, ${\mathit{P}}_{\mathit{k},\mathit{k}-1}$ represents the predicted variance–covariance matrix, and ${\mathit{P}}_{\mathit{k}}$ represents the estimated variance–covariance matrix.

3. Experimental Tests and Evaluations

To present the performance of our method, a set of experiment tests and two simulation tests are arranged and analyzed.

3.1. Data Collection

A field test, based on a vehicle being outfitted with a GPS receiver and a tactical-grade IMU (POS320), was arranged around complex urban environments (tall buildings, large trees, etc.) in Beijing on 17 December 2020, with the corresponding trajectory in Figure 2. The details of the IMU are shown in

Table 1. Technical parameters of the POS320.

IMU	Sampling Rate	Bias		Random Walk
	Hz	Gyro. °/h	Acc. mGal	Angular °/s/√ℎ	Velocity m/s/√ℎ
POS320	200	0.5	25	0.05	0.1

. Since there are no available LEO navigation satellites operated in orbit, we simulated two commonly used low-orbit constellations, namely the inclined orbit constellation and the polar orbit constellation, with a total of 160 satellites providing global coverage. The orbital altitude of the LEO satellites is 970 km. There are 70 satellites in six polar orbits with an orbital inclination of 90 degrees and 90 satellites in 10 inclined orbits with an orbital inclination of 60 degrees. The simulated LEO observations mainly include the satellite clock error, tidal error, Earth rotation, relativistic effect, tropospheric delay, phase winding, and other simulated errors. The sampling rates of the GPS, LEO, and IMU tests were 1 Hz, 1 Hz, and 200 Hz, respectively. The GPS RTK/INS tight integration algorithm provided by the NovAtel Inertial Explorer (IE) software serves to generate the reference values of positions.

3.2. Data Quality Analysis

Based on those measured GPS observations, INS data, and simulated LEO data, six modes, namely GPS, LEO, GPS/LEO integration, GPS/INS tight integration, LEO/INS tight integration, and GPS/LEO/INS tight integration, were analyzed. GPS and LEO modes indicate single-system positioning, and GPS/LEO, GPS/INS, LEO/INS, and GPS/LEO/INS modes indicate integrated positioning between GPS, LEO, or INSs. Since the data quality is one of the most important and direct factors for positioning accuracy, the data quality in terms of available GPS satellite number, LEO satellite number, the Position Dilution of Precision (PDOP), multipath noise, and signal-to-noise ratio (SNR) [42,43] are analyzed in this section.

The PDOP is the spatial geometric intensity factor of satellite distribution. Generally, when the satellite distribution is better, the PDOP value is smaller, and a value generally less than 3 is the ideal state. The number of available GPS and LEO satellites is shown in Figure 3, and the corresponding PDOP is shown in Figure 4. As is shown, the availability of GPS satellites during this test is low, which leads to some epochs with fewer than four satellites and even no available satellites. Relatively, more than five LEO satellites are available during such a period. According to the statistics, the numbers of available satellites on average for GPS, LEO, and GPS/LEO are 7.0, 9.5, and 16.5, respectively, and the corresponding PDOP is 3.4, 2.0, and 1.1. It is clearly visible that both the number of available satellites has increased and the geometry structure has improved. After applying LEO satellites, the smaller the corresponding PDOP, the better the satellite distribution, and theoretically the higher the positioning accuracy.

In addition, the multipath error and SNR are direct indicators for the quality assessment of GPS observations. The multipath error is related to the satellite positions, satellite altitude angle, and user environments. Those mentioned issues could lead to pseudo-range multipath errors with a maximum code width of 0.5, which has a significant impact on positioning accuracy and is difficult to eliminate [43]. As shown in Figure 5 and Figure 6, the pseudo-range multipath on GPS L1-frequency and L2-frequency observations for those observed satellites are within ±3.8 m, with a corresponding RMSE of 0.396 m and 0.489 m. Analyzing the SNR, it can be seen that L2 has a lower SNR than L1, and the lower SNR makes L2 more susceptible to multipath errors.

In view of the above analysis of PDOP, multipath, and SNR, more observations and better spatial geometry can be provided in the integrated mode, which can improve its positioning performance.

3.3. Enhancements in LEO and Robust Theory on GPS/INS Tight Integration

The RMSE was used to evaluate the deviation between the estimated and reference values and the Cumulative Distribution Function (CDF) was used to estimate the probability distribution of the positioning error. Therefore, RMSE and CDF can be used to evaluate positioning errors and their distribution characteristics effectively, respectively:

$$\mathit{RMSE}=\sqrt{\frac{1}{N}{\displaystyle \sum _{i=1}^N{\left(X_i-\widehat{X}\right)}^2}}$$

(39)

$$\mathit{CDF}\left(Y\right)=P\left(Y\le y\right)$$

(40)

where $X_i$ represents the estimated value; $\widehat{X}$ represents the reference value; $N$ and $P$ represent the number of estimates and the probability, respectively; and $Y$ and $y$ represent the random variable and the independent variable, respectively.

3.3.1. Enhancements in Robust Theory on GPS/INS Tight Integration

The n-frame is the basic coordinate system used in INS algorithms. The uniform transformation of the positioning results to the navigation coordinate system can maintain consistency. So, the position differences between the reference values and the solutions calculated by the models of GPS/INS tight integration (GPS/INS), robust GPS/INS tight integration (GPS/INS-RKF), GPS/LEO/INS tight integration (GPS/LEO/INS), and robust GPS/LEO/INS tight integration (GPS/LEO/INS-RKF) were first projected in the n-frame and then analyzed. According to the results and statistics presented in Figure 7 and Figure 8, the performance of GPS/INS and GPS/LEO/INS could be improved after applying the robust algorithm, especially for the GPS/INS tight integration mode. The statistical RMSE values showed that the position results of the GPS/INS-RKF in the north, east, and vertical directions were improved to 2.61 m, 2.34 m, and 2.83 m, respectively, with corresponding improvement percentages of 44.23%, 64.81%, and 72.36% compared to those of GPS/INS. The position results of the GPS/LEO/INS-RKF mode in the three directions were 2.38 m, 1.94 m, and 2.49 m, respectively, with corresponding improvements of 8.11%, 1.52%, and 2.73%. Meanwhile, it can be observed that the GPS/INS mode delivers the worst positioning results compared to other modes, due to the significant enhancement provided by the robust theory and LEO observations. For the GPS/LEO/INS mode, the introduction of LEO can improve its positioning results to a certain extent, so the robust theory has no significant enhancement effect on it than the GPS/INS mode. Since the impact of the robust theory is significant, the following assessments are all based on the robust theory-augmented models.

3.3.2. Enhancements in LEO on GPS/INS Tight Integration

In this section, four modes, namely robust GPS (GPS), robust GPS/LEO integration (GPS/LEO), robust GPS/INS tight integration (GPS/INS), and robust GPS/LEO/INS tight integration (GPS/LEO/INS), were analyzed. According to the results in Figure 9, the GPS-only solutions were significantly affected by the available GPS satellites. After obtaining enhancements from LEO satellites and INSs, the positioning results were significantly improved. By comparing the solutions of GPS/LEO and GPS/LEO/INS with those without LEO augmentation, the average enhancements provided by LEO were about 11.81%, 16.47%, and 10.57% in the north, east, and vertical directions, respectively. Based on the CDF of the position offsets for these four modes in Figure 10, the percentages for those epochs with position offsets within 5 m in the north direction were 92.6%, 98.5%, 93.5%, and 94.7%, respectively, for the GPS, GPS/LEO, GPS/INS, and GPS/LEO/INS modes. Such percentages in the east direction were 93.8%, 94.2%, 95.8%, and 97.4%, and they were 90.3%, 94.4%, 90.5%, and 93.0% in the vertical direction. Given the above CDF trends, the integrated GPS/LEO, GPS/INS and GPS/LEO/INS modes offer better stability and higher accuracy than the single GPS mode.

In general, LEO satellites improve the geometry between the user and the satellite (GPS and LEO satellites) and provide more observational data. Accordingly, the stability, continuity, and accuracy of GPS and GPS/INS were significantly improved after adding LEO enhancements. With the completion of the LEO constellations in the future, the positioning accuracy of the current GPS/INS integrated navigation can be further improved, which is of great significance for integrated navigation models both for vehicle and airborne applications.

3.4. Performance under GPS with Low-Observability Conditions

The results above have proven that both LEO satellites and robust theory are effective in upgrading the performance of GPS and GPS/INS tight integration. In this section, we further evaluate the performance of the presented method under the low observability of GPS satellites, i.e., in complex and harsh urban environments.

In harsh environments, GPS signals may be partially or completely lost, which leads to poor positioning performance. To address this, we simulated the partial GPS loss condition based on the raw GPS observations, and these simulated data were processed by using GPS/LEO and GPS/LEO/INS modes, aiming to analyze the enhancement effects of LEO and INSs under weak GPS conditions. Figure 11 shows the results for the conditions with 1, 2, 3, and 4 available GPS satellites integrated with LEO satellites, which contributed significantly to the GPS positioning accuracy under weak GPS conditions. Along with the decreasing number of available GPS satellites, the fluctuations in the positioning offsets in the north, east, and vertical directions were kept small and tended to stabilize in general, which is due to the augmentation of LEO satellites. According to the RMSE for the position in the GPS/LEO mode using 1, 2, 3, and 4 available GPS satellites in Figure 12, the RMSEs of GPS/LEO using four available GPS satellites were 3.90 m, 3.66 m, and 4.10 m in the north, east, and vertical directions, respectively, and such values when using one available GPS satellite were 4.30 m, 5.54 m, and 4.96 m in the three directions. Significantly, the impact of the low observability of GPS satellites is constrained visibly after introducing the LEO satellites, even when only one GPS satellite is available.

Figure 13 shows the results for 1, 2, 3, and 4 available GPS satellites under GPS/LEO/INS tight integration. Under weak GPS conditions, LEO and INSs contributed significantly to the GPS positioning results with smooth fluctuations during the whole period. Similar to the conclusions in the above test, there is no visible position accuracy change while the number of available GPS satellites decreases from 4 to 1. According to the RMSE in Figure 14, the positioning RMSE was 3.64 m, 2.16 m, and 2.79 m in the north, east, and vertical directions when using one available GPS satellite, and the RMSEs of these using four GPS satellites were 2.38 m, 1.94 m, and 2.49 m, respectively. For a pseudo-range-based vehicle positioning method, such a positioning accuracy loss is acceptable. The results in this section show that LEO satellites are helpful in upgrading the pseudo-range-based positioning.

3.5. Performance under LEO with Low-Observability Conditions

Similar to GPS, LEO signals may also be partially or completely lost in complex and harsh environments. To present the performance of our method under such conditions, we simulated partial LEO outages. Figure 15 shows the results for the conditions with 1, 2, 3, and 4 available LEO satellites to integrate with GPS satellites, which contributed significantly to the LEO positioning accuracy under weak LEO conditions. Along with the decreasing number of available LEO satellites, the fluctuations in the positioning offsets in the north, east, and vertical directions were kept small and tended to stabilize in general, which is due to the augmentation of GPS satellites. According to the RMSE for the position in the GPS/LEO mode using 1, 2, 3, and 4 available LEO satellites in Figure 16, the RMSEs of GPS/LEO using 1 available LEO satellite were 3.99 m, 4.12 m, and 3.94 m in the north, east, and vertical directions, respectively, and such values when using 4 available LEO satellites were 2.4 m, 3.41 m, and 3.50 m in the north, east, and vertical directions, respectively. Significantly, with the increase in available LEO satellites, the performance of GPS satellite positioning has improved to a certain extent.

Figure 17 shows the results for 1, 2, 3, and 4 available LEO satellites under GPS/LEO/INS tight integration. Under weak LEO conditions, GPS and INSs contributed significantly to the LEO positioning results with smooth fluctuations during the whole period. Similar to the conclusions in the above test, there was no visible position accuracy change while the number of available LEO satellites decreased from 4 to 1. According to the RMSE in Figure 18, the positioning RMSE was 2.635 m, 2.463 m, and 2.739 m in the north, east, and vertical directions when using 1 available LEO satellite, respectively, while the RMSEs of these when using 4 available LEO satellites were 2.519 m, 2.05 m, and 2.540 m, respectively. So the results in this section show that enhancing the integrated navigation approach with LEO observations has explicitly improved vehicle navigation in harsh environments and allows for more stable and extensive applications.

4. Discussion

The above field tests and simulation tests have proven the contributions of our method in upgrading the positioning accuracy for vehicles under harsh environments. As a summary, we list the RMSE values of position offsets for the GPS, GPS/LEO, GPS/INS, GPS/INS-RKF, GPS/LEO/INS, and GPS/LEO/INS-RKF modes in Figure 8 and Figure 19. It can be found that positioning accuracy can be improved by using the robust theory and LEO satellites. In general, the position improvement percentages provided by the robust algorithm are 26.17%, 33.17%, and 37.55%, and such percentages provided by LEO are 11.81%, 16.47%, and 10.57% in the north, east, and vertical directions, respectively. In addition, the INS brings 20.55%, 37.10%, and 26.37% improvements, respectively. The reasons are as follows: (1) The robust theory in Section 2.4 has the capability to restrain or remove the impact of low-accuracy observations and to smooth the position time series. (2) The LEO satellites have lower orbits, faster speeds, and stronger signals. They have two benefits: (a) they improve the geometry structure between users and satellites (GPS and LEO), and (b) they provide many more observations that can be used in parameter estimation. (3) Due to the introduction of LEO, the enhancement in the GPS/LEO/INS mode was smaller than that in the GPS/LEO mode, but the enhancement in the results in the GPS/LEO/INS mode was very large. The results in this mode from the introduction of an INS are due to the enhancements from the INS [27,28,29].

5. Conclusions

The meter-level positioning method is widely needed for vehicle-related applications at present. Currently, the GPS pseudo-range-based SPP is widely used. However, such SPP has visible drawbacks, especially for the continuity and reliability in complex and harsh environments. To solve two such problems, this paper presented a robust theory and LEO-enhanced GPS/INS tight integration method. To validate such a method, a set of vehicle data and simulated LEO data were used. The corresponding results illustrated the following: (1) The robust theory can control low-quality observations and abnormal observations, with position improvements on average about 32.29%. (2) The LEO satellites can improve the geometry structure and provide many more observations. The position improvements on average are about 12.95%. (3) The GPS/LEO/INS tight integration has the highest positioning accuracy, with RMSEs of 2.38 m, 1.94 m, and 2.49 m in the north, east, and vertical directions, respectively, with corresponding improvement percentages of 30.21%, 47.43%, and 34.13% compared to those of GPS, and such percentages are 8.60%, 17.24%, and 12.14% compared to those of GPS/INS, respectively. Based on our works in this paper, we will focus on precise positioning by using the GNSS, LEO, and INS measurements in the future. At the same time, it is necessary to focus on the differences between LEO and GNSS systems, such as LEO’s data quality control, cycle slip detection, and error correction.