A Sequential Student’s t-Based Robust Kalman Filter for Multi-GNSS PPP/INS Tightly Coupled Model in the Urban Environment

: The proper stochastic model of a global navigation satellite system (GNSS) makes a signiﬁcant difference on the precise point positioning (PPP)/inertial navigation system (INS) tightly coupled solutions. The empirical Gaussian stochastic model is easily biased by massive gross errors, deteriorating the positioning precisions, especially in the severe GNSS blockage urban environment. In this paper, the distributional characteristics of the gross-error-contaminated observation noise are analyzed by the epoch-differenced (ED) geometry-free (GF) model. The Student’s t distribution is used to express the heavy tails of the gross-error-contaminated observation noise. Consequently, a novel sequential Student’s t-based robust Kalman ﬁlter (SSTRKF) is put forward to adjust the GNSS stochastic model in the urban environment. The proposed method pre-weights all the observations with the a priori residual-derived robust factors. The chi-square test is adopted to distinguish the unreasonable variances. The proposed sequential Student’s t-based Kalman ﬁlter is conducted for the PPP/INS solutions instead of the standard Kalman ﬁlter (KF) when the null hypothesis of the chi-square test is rejected. The results of the road experiments with urban environment simulations reveal that the SSTRKF improves the horizontal and vertical positioning precisions by 57.5% and 62.0% on average compared with the robust Kalman ﬁlter (RKF).


Introduction
Nowadays, GNSS is widely applied for vehicle navigation due to its high accuracy, high availability, and global coverage. Multi-GNSS PPP can achieve a high-accuracy kinematic positioning performance with a sole receiver, attributing to the good satellite geometry. However, high buildings, tunnels, and other typical urban environments can severely impact the observation quality, leading to the biased positioning solutions. The short-term well-performed INS is widely applied to assist environment-dependent PPP in the sophisticated urban environment; furthermore, the PPP/INS tightly coupled system can effectively improve the long-term precise navigation. In the urban environment, the line-of-sight (LOS) multipath effect and the non-line-of-sight (NLOS) reception are two common error sources impacting on the GNSS observation quality. The LOS multipath errors happen when the receiver antenna intercepts both the LOS and multipath satellite signal. This case is mainly attributed to a mixture of reflection and scattering echoes [1]. The LOS multipath errors can reach up to tens of meters in code and a quarter cycle in phase [2,3], which can be modelled by the observation noise variance changes [4]. In the NLOS case, an obstacle completely blocks the LOS satellite signal, but the signal still reaches the antenna mainly through reflection. The NLOS reception is a more serious problem in the urban environment, also leading to significant biases in both code and phase observations [1]. Reference [4] modelled the NLOS reception errors as the observation noise mean value jumps. Hence, due to the impacts of the LOS multipath and NLOS reception

Multi-GNSS PPP/INS Tightly Coupled Model
In the following, the PPP/INS tightly coupled functional and stochastic model are described. We first derive the multi-GNSS UDUC PPP observation model, involving GPS, GLONASS, Galileo, and BDS. Then, the north-oriented INS state model is formulated.

Multi-GNSS PPP/INS Tightly Coupled Observation Model
After applying the known corrections of the precise satellite products, i.e., orbits, clock offsets, and satellite differential code biases, the linearized UDUC PPP observation equations of code P k,r and phase L k,r for m available satellites read [20,21] P k,r = (e 2 ⊗ A k,r )p e k,r + e 2 ⊗ g k,r T s k,r + e 2m δt G k,r,IF + (e 2 ⊗ Λ m )Σ k + U ι k + ε P k,r (1) L k,r = (e 2 ⊗ A k,r )p e k,r + e 2 ⊗ g k,r T s k,r + e 2m δt G k,r,IF + (e 2 ⊗ Λ m )Σ k − U ι k − a r + ε L k,r (2) where the subscripts k and r and the superscript s refer to epoch, receiver, and satellite, respectively. In the following, the superscripts G, R, E, and C represent GPS, GLONASS, Galileo, and BDS, respectively. P k,r = (P 1 ) T (P 2 ) T T is the vector of the dual-frequency multi-GNSS code observations with P j = P G k,r,j T P R k,r,j T P E k,r,j T P C k,r,j T T (j = 1, 2), where the subscript j denotes the frequency. P sys k,r,j = P 1 k,r,j · · · P m sys k,r,j T is the jth frequency code observation vector of the m sys satellites for system sys. Note that the dual-frequency multi-GNSS UDUC phase observation vector L k,r can be organized in a same structure as P k,r . A k,r is the design matrix to the earth-centered earthfixed (ECEF) coordinate parameters p e k,r with e representing the ECEF frame (e-frame).
is the mapping function vector of the zenith tropospheric delay (ZTD) T s k,r where g sys k,r = g 1 k,r · · · g m sys k,r T . To eliminate the rank deficiency between the clock offset errors and the code hardware delays for four systems, two processing steps are introduced in new reparameterization: (i) The receiver clock error is first parameterized as δt G k,r,IF , absorbing the ionospheric-free (IF) combination of the GPS code hardware delay. (ii) The inter-system code biases (ISCB) vector Σ k = η T k,GR η k,GE η k,GC T is defined referring to the δt G k,r,IF , where η k,GR = δt R k,r,IF − δt G k,r,IF e m R , η k,GE = δt E k,r,IF − δt G k,r,IF , and η k,GC = δt C k,r,IF − δt G k,r,IF . Note that GLONASS is a frequency division multiple access (FDMA) system, leading to the satellite-specific receiver code hardware delays absorbed into the clock errors δt R k,r,IF . e m is an m-column vector with all elements being 1, while I m is an m-dimension identity matrix.
is the parameterized ionospheric delay vector, where ι sys k = ι 1 k,r · · · ι m sys k,r T absorbs the code hardware delay residuals.
is the float ambiguity vector, where a sys k = a 1 k,r · · · a m sys k,r T absorbs the code and phase hardware delays. The code and phase noise vectors are ε P k,r and ε L k,r . Moreover, ⊗ denotes the Kronecker product. The PPP/INS tightly coupled state parameters is organized as where δφ k , δv n k , δp n k , ∇ k , and ε k are the attitude error vector, velocity and position error vectors under the navigation frame (n-frame) and the accelerometer and gyroscope bias vectors, respectively. The observation functional model can be constructed as follows where H y,k is the design matrix to x k . Here, ρ INS,k is the INS-derived geometry distance vector. G k = A k C e n denotes the INS position error design matrix in the navigation frame (n-frame), where C e n is the rotation matrix from the n-frame to the e-frame. The empirical multi-GNSS Gaussian stochastic model is constructed as Q k = Q m ⊗ Q 0 with Q m = blkdiag(Q P , Q L ) and Q 0 = blkdiag(Q G , Q R , Q E , Q C ) defining the elevationdependent dispersion for the four systems [22]. Moreover, Q P = diag σ 2 P 1 , σ 2 P 2 and Q L = diag σ 2 L 1 , σ 2 L 2 describe the code and phase precisions.

North-Oriented INS State Model
The PPP/INS tightly coupled state functional model is based on the north-oriented INS error function as follows [23] where the points upon δφ, δv n k , and δp n k denote the derivative; the superscript n represents the n-frame; and δ denotes the error. ω n in , ω n ie , and ω n ib are the n-frame, e-frame, and body frame (b-frame) angular rate relative to the inertial frame (i-frame), projected to the n-frame. ω n en can be expressed in the same way above. f n is the specific force in the n-frame. Moreover, ∇ k , and ε k are contained in δf n and δω n ib , respectively. M pv and M pp are the coefficient matrices of δv n and δp in the position error function, which read where ϕ and h are the latitude and altitude, respectively. v E and v N are the east and north velocity. R M and R N represent the principal radii of curvature along the meridional and prime-vertical normal sections, respectively.

PPP/INS Observation Noise Distributional Characteristic Analysis
The PPP/INS observation noise distributional characteristic is analyzed as follows. First, the raw GNSS observation noise is extracted and analyzed by the ED GF model. Then, we analyze the distributional characteristics of the gross-error-contaminated observation noise.

Multi-GNSS Observation Noise Extraction
Two types of GF models are formulated to extract the phase and code observation noise [24] GF s 1,k,r = L s k,r,1 − L s k,r,2 = a s r,2 − a s r,1 + (µ 2 − 1)ι s k,r + δB k,r,12 + δb s k,12 + ε GF s 1,k,r GF s 2,k,r,j = P s k,r,j − L s k,r, where δB k,r,12 = B k,r,1 − B k,r,2 and δb k,r,12 = b k,r,1 − b k,r,2 denote the inter-frequency receiver and satellite phase hardware delay differences, respectively. D k,r,j is the frequency-specific receiver code hardware delay. ε GF s 1,k,r and ε GF s 2,k,r,j are the unmodelled error residuals of GF s 1,k,r and GF s 2,k,r,j , involving the LOS multipath errors, the NLOS reception errors, and the GF observation noise. Moreover, the phase hardware delays and unmodelled error residuals are evidently smaller than those of code [25]. Thus, Equation (8) can be simplified as All cycle slips are detected and repaired by applying the TurboEdit method in the post processing [26,27]. Hence, the hardware delays and the phase ambiguities can be eliminated in the epoch difference GF combinations ∆GF s 1,k,r and ∆GF s 2,k,r,j , which read where only the ionospheric delay residuals and the unmodelled error residuals are contained in ∆GF s 1,k,r and ∆GF s 2,k,r,j . Moreover, satellites with large a priori residuals (an empirical threshold 30 m) are rejected to guarantee the observations to be gross-error-free. Afterwards, 1Hz sampling static multi-GNSS data from the IGS MGEX GAMG station (7:00:00-7:59:59, 21 February 2022) is applied to analyze the observation noise. To eliminate the influence of the multipath errors as much as possible, the cut-off elevation is set to 35 • . Moreover, the ionospheric delay residuals in ∆GF s 1,k,r and ∆GF s 2,k,r,j are ignored in this case, as the stability of ED ionospheric delays within high-rate (1Hz) data [28,29]. Eventually, the phase and code observation noise are extracted by ∆GF s 1,k,r and ∆GF s 2,k,r,j . Considering that the observation precisions are satellite-type-specific, the distributional characteristics of GPS, GLONASS, Galileo, BDS geostationary earth orbit (GEO), BDS inclined geosynchronous orbit (IGSO), and BDS medium earth orbit (MEO) observation noise are analyzed afterwards in Figures 1 and 2. In the following, N(mean, cov) denotes the Gaussian probability density functions (PDF) with the mean vector mean and the covariance matrix cov. Moreover, var(· · · ) is the variance derivation operation. Figure 1 indicates that the phase observation noise can be well fitted by the zero-mean-Gaussian distribution in absence of the gross errors. Moreover, the same conclusion is obtained for the code observation noise by analyzing Figure 2. Note that the Gaussian characteristic can also be extracted from ∆GF s 2,k,r,2 , which is not further elaborated here. In this case, the Gaussian-based KF can properly handle the PPP/INS solutions. The distributional characteristics of the gross-error-contaminated observations are discussed in the next section.
Eventually, the phase and code observation noise are extracted by ∆ , , and ∆ , , , . Considering that the observation precisions are satellite-type-specific, the distributional characteristics of GPS, GLONASS, Galileo, BDS geostationary earth orbit (GEO), BDS inclined geosynchronous orbit (IGSO), and BDS medium earth orbit (MEO) observation noise are analyzed afterwards in Figures 1 and 2. In the following, N , denotes the Gaussian probability density functions (PDF) with the mean vector and the covariance matrix . Moreover, var ⋯ is the variance derivation operation.  Figure 1 indicates that the phase observation noise can be well fitted by the zeromean-Gaussian distribution in absence of the gross errors. Moreover, the same conclusion is obtained for the code observation noise by analyzing Figure 2. Note that the Gaussian characteristic can also be extracted from ∆ , , , , which is not further elaborated here. In this case, the Gaussian-based KF can properly handle the PPP/INS solutions. The distributional characteristics of the gross-error-contaminated observations are discussed in the next section.

Distributional Characteristic Analysis on Gross-Error-Contaminated Observation Noise
The gross error simulations were introduced into the raw observations. According to the research by reference [1], in a typical city urban environment (Shanghai-Lujiazui area), 95% of the LOS multipath error lifetime is shorter than 13.6 s when the rover velocity is 0~3 km/h. Hence, the 10 s LOS multipath errors were introduced into the G08 raw  [4] as Moreover, reference [1] reveals that most LOS multipath errors are less than 350m, while the observations with much too large multipath errors are always eliminated by the a priori residual test (e.g., an empirical threshold 30 m) in PPP data pre-processing. Hence, a typical multipath standard deviation 15 m is selected to simulate the LOS multipath errors, which can pass the priori residual test. The phase tracking errors can be up to a quarter cycle (maximum 0.064 m for GPS L5, Galileo E5a, and BDS B2a and minimum 0.047 m for the last satellite in GLONASS L1 among all the systems). Consequently, introducing the 0.01 m phase variance jumps into the phase observations is a reasonable choice.
As for the NLOS reception errors, 95% of their lifetime is shorter than 68.2 s with 0~3 km/h rover velocity [1]. Hence, the NLOS reception errors can be modelled as the 60 s constant mean value jumps in the code and phase. However, ∆GF s 1,k,r and ∆GF s 2,k,r,j cannot extract the NLOS reception errors modelled as constant mean value jumps which are eliminated by the epoch-differenced strategy. Hence, only the LOS multipath errors are simulated to accomplish the analysis. The ∆GF s 1,k,r PDF of G08, Student's t distribution PDF (DOF = 3), and N 0, var ∆GF s 1,k,r are shown in Figure 3.
ens. 2022, 14, x FOR PEER REVIEW As for the NLOS reception errors, 95% of their lifetime is shorter than km/h rover velocity [1]. Hence, the NLOS reception errors can be mode constant mean value jumps in the code and phase. However, ∆ , , and not extract the NLOS reception errors modelled as constant mean value ju eliminated by the epoch-differenced strategy. Hence, only the LOS multi simulated to accomplish the analysis. The ∆ , , PDF of G08, Student' PDF (DOF 3), and N 0, var ∆ , , are shown in Figure 3. , and Student's t PDF (DOF 3) Figure 3 demonstrates that the distributional characteristics of the g taminated observation noise can be divided into the Gaussian part and When fitting the heavy tail, the Student's t distribution with a proper DO preferred to the Gaussian distribution. However, the Gaussian part is o fitted by the Gaussian distribution than the Student's t distribution. Henc and Student's t combined processing strategy will outperform to any sing based strategy. Eventually, if the NLOS reception errors were introduced bution would obtain a heavier tail due to the mean value jump model.  Figure 3 demonstrates that the distributional characteristics of the gross-error-contaminated observation noise can be divided into the Gaussian part and the heavy tail. When fitting the heavy tail, the Student's t distribution with a proper DOF parameter is preferred to the Gaussian distribution. However, the Gaussian part is obviously better fitted by the Gaussian distribution than the Student's t distribution. Hence, the Gaussian and Student's t combined processing strategy will outperform to any single distribution-based strategy. Eventually, if the NLOS reception errors were introduced, the true distribution would obtain a heavier tail due to the mean value jump model.

Sequential Student's t-Based Robust Kalman Filter
First, the robust Kalman filter and the Student's t based Kalman filter (STKF) are introduced briefly. Then, the SSTRKF is derived by incorporating the RKF and the sequential Student's t Kalman filter. Eventually, we provide the SSTRKF processing strategy.

Robust Kalman Filter Based on IGG-III Function
The RKF tunes the gross-error-contaminated observation noise in the standard Kalman filter by setting the robust factors as follows where Q k is the new stochastic model with Q k,i = r k,i Q k,i . Q k,i represents the ith element in the main diagonal of the raw observation stochastic model Q k . r k,i is the ith observation robust factor, derived by the IGG-III function [30] where k 0 and k 1 are presetting constants, chosen from 1.0 to 2.5 and from 3.5 to 8.0, respec- where v k,i and Q v k,i are the ith observation a priori residual and the corresponding variance, respectively. Note that Q v k,i is the ith element in the main diagonal of the a priori residual covariance matrix Q v k = H y,k Q k H T y,k + Q k . Here, Q k is the covariance matrix of the one-step prediction x k . After redetermining the stochastic model of the observation system, the modified Kalman gain matrix can be written as The RKF state estimators and the corresponding covariance matrices will be obtained after updating the Kalman gain matrix. As shown in Equation (14), the robust factors vary with v k,i which is determined by Q v k,i in Equation (15). Meanwhile Q v k,i is significantly susceptible to the worse PPP/INS positioning solutions with poor satellite geometry.

Student's t-Based Kalman Filter
The Student's t distribution is applied to fit the gross-error-contaminated observation noise. The STKF can be described as follows [18,19] where St y k ; H y,k x k , Q k , ϑ k denotes the Student's t PDF with mean vector H y,k x k , scale matrix Q k , and DOF ϑ k . ξ k is an auxiliary random variable obeying Gamma distribution G ξ k ; ϑ k 2 , ϑ k 2 . The STKF adjusts the observation stochastic model mainly by ξ k . Moreover, ϑ k is also modelled as Gamma distribution with shape parameters p k and rate parameter q k . The estimation problem can be solved by the VB method as follows [19] {q(x k ), q(ξ k ), q(ϑ k )} = arg min KLD(q(x k ), q(ξ k ), q(ϑ k )||p(x k , ξ k , ϑ k |y 1:k )) (19)  where q(· · · ) and p(· · · ) denote the approximate posterior and posterior PDF, respectively. KLD(· · · ) is the Kullback-Leibler divergence. Equations (17) and (18) shows that the STKF simultaneously adjusts the stochastic model of every observation with solely one parameter ξ k , which can mis-weight the gross-error-free observations. Moreover, the gross-errorcontaminated observation noise has diverse stochastic characteristics, requiring separately processing. It is noticed that the computation burden of the STKF is relatively heavy as Equation (19) is solved iteratively. Although the Student's t distribution degenerates to be Gaussian as the DOF tends to be infinite, conducting STKF with high-dimension matrices in every epoch can lead to disastrous computation burdens. Hence, the SSTRKF is proposed in the next section.

Sequential Student's t-Based Robust Kalman Filter
The proposed SSTRKF first pre-weights the observations by Equation (14), which can give the proper observation noise covariance matrix with ideal satellite geometry and gross-error-free observation noise. Hence, after being pre-weighted by the IGG-III function, we conduct the chi-square test satellite-by-satellite as follows where v s k and Q s v,k = H s y,k Q k H s y,k T + Q s k are the a priori residual vector and the corresponding IGG-III-modified covariance matrix of the satellite s, respectively. th is the presetting chi-square threshold. Furthermore, H 0 and H 1 represent the null and alternative hypothesis, respectively. If the H 0 case is accepted, the IGG-III-derived variances of the satellite s are reasonable. In contrast, the observations are sorted into the unreasonable variance group when accepting H 1 . The reasonable variance group obeys the Gaussian assumption, whereas the unreasonable variance group is modelled by the Student's t distribution. Hence, we first process the reasonable group by Equation (16).
The unreasonable group is sequentially processed due to the diverse stochastic characteristics of the gross-error-contaminated code and phase observations. Hence, the estimation problem of the unreasonable group can be presented as p x k,i y 1:k−1 , y k,N , y k,ST,1: where y k,ST,1:i = y k,ST,1 · · · y k,ST,i is the observation vector with the first i observations in the unreasonable group. y k,N denotes the reasonable variance group observation vector. Note that the subscript i represents the ith observation in the unreasonable group. Q k,ST,i is the ith unreasonable observation variance, determined by Equation (13). Then, the estimation problem can be solved using the VB method as Equation (19).  k,i are the tth iteration estimator and the corresponding covariance matrix, respectively. dim is the observation dimension. tr[· · · ] is the trace of a matrix. Furthermore, ξ k can be calculated through the expectation of the Gamma distribution as is exploited to be the posterior PDF of the DOF ϑ k,i , resulting in where ψ(· · · ) denotes the digamma function. The expectation E (t+1) [ϑ k,i ] is obtained as Then, the posterior PDF of x k,i is modelled by the Gaussian PDF N x k,i ;x . Thus, we obtain the measurement update process, given as More derivation details refer to references [18] and [19]. To intuitively demonstrate the process of SSTRKF, the flow chart is shown in Figure 4. Moreover, we elaborate the processing steps in the kinematic applications: (i) After the time update, the IGG-III pre-weighted observation stochastic model is calculated by Equations (12)- (14). (ii) The chi-square test is applied to sort the available satellites into the reasonable and unreasonable group. (iii) The measurement update with the Kalman gain obtained by Equation (16) is conducted for the reasonable group. (iv) The estimator and the corresponding covariance matrix of the reasonable group are assigned to x k,1 and Q k,1 , respectively. (v) p

Experiments and Discussions
The road experiments were conducted with SPAN-CPT IMU and a Trimble receiver in March 2019, Xuzhou, China. The dataset is available in the open-access software Ginav [31]. The experiment duration is 1721 s with 100 Hz IMU data and 1 Hz GNSS data. Note that we applied the RTK/INS tightly coupled solutions as the reference to evaluate the navigation performance. The precise clock and orbit products were obtained from the MGEX dataprocessing center (http://www.igs.org/mgex/products, accessed on 31 October 2021), and the satellite and receiver antenna parameters were gained from igs14.atx. In real-time navigation, the real-time clock and orbit products can be used to support PPP/INS. Figure 5 is the 2D experiment trajectory, containing two laps on a same road. We select three segments of the road to conduct the six severe urban environment simulations, presented in red in Figure 5. To simulate the blockage urban environment as severe as possible, the cut-off elevation was set to 90 • and 60 • for the satellites with azimuth 285 •~3 60 • and 0 •~2 85 • , respectively. This environment can be always found in the urban environment, such as the narrow streets surrounded by high buildings. Moreover, the LOS multipath and the NLOS reception errors were introduced into an available satellite and the blocked satellite with the highest elevation, respectively. According to the analysis in Section 3.2, the lifetime and simulation model of the LOS multipath and the NLOS reception errors are summarized in Table 1. It is noticed that the elevation of the satellite contaminated by the NLOS reception errors is higher than 50 • , which tends to lead the NLOS reception errors lower than 50 m [1]. Therefore, the value of the NLOS reception errors was set to 15 m with the same consideration as the LOS multipath errors in Section 3.2. Figure 5 is the 2D experiment trajectory, containing two laps on a same road. We select three segments of the road to conduct the six severe urban environment simulations, presented in red in Figure 5. To simulate the blockage urban environment as severe as possible, the cut-off elevation was set to 90° and 60° for the satellites with azimuth 285°~360° and 0°~285°, respectively. This environment can be always found in the urban environment, such as the narrow streets surrounded by high buildings. Moreover, the LOS multipath and the NLOS reception errors were introduced into an available satellite and the blocked satellite with the highest elevation, respectively. According to the analysis in Section 3.2, the lifetime and simulation model of the LOS multipath and the NLOS reception errors are summarized in Table 1. It is noticed that the elevation of the satellite contaminated by the NLOS reception errors is higher than 50°, which tends to lead the NLOS reception errors lower than 50 m [1]. Therefore, the value of the NLOS reception errors was set to 15 m with the same consideration as the LOS multipath errors in Section 3.2.    Model P s k,r,j = P s k,r,j + 15m L s k,r,j = L s k,r,j + 0.01m The available satellite (AS) number of the multi-GNSS (GPS/GLONASS/Galileo/BDS), GPS (G), GLONASS (R), Galileo (E), BDS (C) is shown in Figure 6. During the experiment, the average AS number of multi-GNSS, G, R, E, and C are 18.9, 5.8, 2.5, 3.0, and 7.5, respectively. The first urban environment simulation was from 80 s to 125 s with 6, 1, 0, 1, and 4 average AS number for multi-GNSS, G, R, E, and C. During the 320 s~350 s, we conducted the second urban environment simulation, where the AS number drops to 4.3, 1, 0, 0, and 3.3 for multi-GNSS, G, R, E, and C on average, respectively. As for the third urban environment simulation during 540 s~560 s, the multi-GNSS, G, R, E, and C AS number maintain 4, 1, 0, 0, and 3, respectively. The fourth, fifth and sixth urban environment simulations were from 945 s to 1000 s, from 1227 s to 1260 s, and from 1445 s to 1465 s, respectively, where 4, 1, 0, 0, and 3 satellites for multi-GNSS, G, R, E, and C are available on average. Figure 6 indicates that the AS number reduces significantly during the urban environment simulations. Moreover, the shadow areas in Figures 6-9 represent the urban environment simulations.
To reveal the impacts of the LOS and NLOS multipath errors on the positioning solutions in the urban blockage environment, we first conducted three schemes with the KF, summarized in Table 2. Figure 7 depicts the positioning errors of the three schemes, where the Scheme 1 achieves the best precisions among the three schemes. When suffering the blockage environment, the positioning errors (blue line in Figure 7) increases due to the barren observation conditions. Moreover, the positioning errors diverge dramatically when the LOS multipath and NLOS reception errors are introduced into the observations, such as the red line in Figure 7. Therefore, the gross errors induced by the LOS multipath effect and NLOS reception are rather serious problems for the PPP/INS tightly coupled system, especially in the urban severe blockage environment. the urban environment simulations. Moreover, the shadow areas in Figures 6-9 represent the urban environment simulations.
To reveal the impacts of the LOS and NLOS multipath errors on the positioning solutions in the urban blockage environment, we first conducted three schemes with the KF, summarized in Table 2. Figure 7 depicts the positioning errors of the three schemes, where the Scheme 1 achieves the best precisions among the three schemes. When suffering the blockage environment, the positioning errors (blue line in Figure 7) increases due to the barren observation conditions. Moreover, the positioning errors diverge dramatically when the LOS multipath and NLOS reception errors are introduced into the observations, such as the red line in Figure 7. Therefore, the gross errors induced by the LOS multipath effect and NLOS reception are rather serious problems for the PPP/INS tightly coupled system, especially in the urban severe blockage environment.  Table 2. Three simulation schemes to analyze the impacts of the urban environment.

The LOS Multipath and NLOS Reception Errors
Blockage Environment Yes Yes The chi-square test threshold was set to 13.277 with the significance level and the DOF being 1% and 4, respectively. Figure 8 gives the chi-square statistics during the road experiment. The chi-square statistics average 0.23 when neither the LOS multipath errors nor the NLOS reception errors are introduced. Moreover, only 15 in total 31,360 chi-square statistics exceed the chi-square threshold in this case, indicating that the IGG-III scheme can properly adjust the observation noise stochastic model in this case. On the contrary, the average chi-square statistics of the satellites impacted by the LOS multipath and NLOS reception errors are 41.6, 68.1, 89.3, 74.0, 52.8, and 72.0 for the six urban environment simulations, respectively. The chi-square statistics of the gross-error-contaminated satellites evidently exceed the chi-square threshold during the six urban environment simulations. Consequently, the IGG-III-based RKF fails to give proper variances for the gross-errorcontaminated observations. Moreover, the SSTRKF avoids an unnecessary computation burden from processing every observation with the sequential Student's t-based Kalman filter by dividing the observations into the reasonable and unreasonable variance group. The chi-square test threshold was set to 13.277 with the significance level and the DOF being 1% and 4, respectively. Figure 8 gives the chi-square statistics during the road experiment. The chi-square statistics average 0.23 when neither the LOS multipath errors nor the NLOS reception errors are introduced. Moreover, only 15 in total 31,360 chi-square statistics exceed the chi-square threshold in this case, indicating that the IGG-III scheme can properly adjust the observation noise stochastic model in this case. On the contrary, the average chi-square statistics of the satellites impacted by the LOS multipath and NLOS reception errors are 41.6, 68.1, 89.3, 74.0, 52.8, and 72.0 for the six urban environment simulations, respectively. The chi-square statistics of the gross-error-contaminated satellites evidently exceed the chi-square threshold during the six urban environment simulations. Consequently, the IGG-III-based RKF fails to give proper variances for the gross-errorcontaminated observations. Moreover, the SSTRKF avoids an unnecessary computation burden from processing every observation with the sequential Student's t-based Kalman filter by dividing the observations into the reasonable and unreasonable variance group. In the proposed SSTRKF, the sequential Student's t-based Kalman filter is activated to process the observations that cannot be appropriately weighted by the IGG-III model. According to Section 4.3, , and , were initialized to 0.3 and 1 before processing every observation, respectively. The maximum iteration number of the SSTRKF was set to 10 [17][18][19]. To explore the superiority of the proposed SSTRKF, the KF, the RKF, the SSTRKF, and the STKF were applied to the multi-GNSS PPP/INS tightly coupled system in the Scheme 3 setting. precision improvement are achieved by conducting the SSRTKF. Moreover, the SSTRKF improves the positioning solutions slightly in the simulation 1. As shown in Figure 6, the observation condition of simulation 1 is better than the other five simulations, contributing to well-independent positioning. Hence, all three schemes can properly identify and handle the gross-error-contaminated satellites. According to the analysis above, the proposed SSTRKF evidently outperforms to the KF and RKF during the urban environment simulations, especially in the severe blockage environment. Figure 9. Positioning errors of the three filters. Figure 9. Positioning errors of the three filters. Table 2. Three simulation schemes to analyze the impacts of the urban environment.

The LOS Multipath and NLOS Reception Errors
Blockage Environment The chi-square test threshold was set to 13.277 with the significance level and the DOF being 1% and 4, respectively. Figure 8 gives the chi-square statistics during the road experiment. The chi-square statistics average 0.23 when neither the LOS multipath errors nor the NLOS reception errors are introduced. Moreover, only 15 in total 31,360 chi-square statistics exceed the chi-square threshold in this case, indicating that the IGG-III scheme can properly adjust the observation noise stochastic model in this case. On the contrary, the average chi-square statistics of the satellites impacted by the LOS multipath and NLOS reception errors are 41.6, 68.1, 89.3, 74.0, 52.8, and 72.0 for the six urban environment simulations, respectively. The chi-square statistics of the gross-error-contaminated satellites evidently exceed the chi-square threshold during the six urban environment simulations. Consequently, the IGG-III-based RKF fails to give proper variances for the gross-errorcontaminated observations. Moreover, the SSTRKF avoids an unnecessary computation burden from processing every observation with the sequential Student's t-based Kalman filter by dividing the observations into the reasonable and unreasonable variance group.
In the proposed SSTRKF, the sequential Student's t-based Kalman filter is activated to process the observations that cannot be appropriately weighted by the IGG-III model.
According to Section 4.3, p (0) k,i and q (0) k,i were initialized to 0.3 and 1 before processing every observation, respectively. The maximum iteration number of the SSTRKF was set to 10 [17][18][19]. To explore the superiority of the proposed SSTRKF, the KF, the RKF, the SSTRKF, and the STKF were applied to the multi-GNSS PPP/INS tightly coupled system in the Scheme 3 setting. Figure 9 shows the KF, the RKF, the STKF, and the SSTRKF positioning errors during the road experiments. The STKF simultaneously tunes the observation noise covariance matrix with solely one factor, leading to the severely biased PPP stochastic model. Hence, the STKF errors diverge rapidly. It is evident that the positioning precisions are impacted by the urban environment simulation. The root mean square (RMS) positioning errors during the six urban environment simulations are collected in Figure 10. Compared with the KF, the proposed SSTRKF improves the horizontal and vertical precisions by 65.7% and 72.1% on average, respectively. Moreover, the SSTRKF achieves 57.5% and 62.0% average precision improvement relative to the RKF in the horizontal and vertical components, respectively. During simulation 2, although the RKF ameliorates the positioning precisions compared with the KF, the RKF positioning errors still diverge to 1.78m in maximum. As for the SSTRKF, the positioning precisions are improved by 29.4% and 32.6% relative to the RKF in the horizontal and vertical components during simulation 2. Moreover, the RKF impairs the positioning precisions compared with the KF during simulation 3, indicating that the RKF fails to derive proper observation variances with the IGG-III function. In contrast, during simulation 3, the SSTRKF suppresses the error divergence during simulation 3, reducing the horizontal and vertical positioning errors by 84.9% and 71.0% compared with the RKF and by 83.9% and 68.6% compared with the KF, respectively. In simulation 4, compared with RKF and KF, the horizontal and vertical PPP/INS precisions are improved by 80.3% and 95.3% and by 87.5% and 96.3% after applying the SSTRKF. As for simulation 5, the horizontal and vertical positioning errors of RKF and KF are reduced by 72.6% and 81.8% and by 54.2% and 76.7% with the SSTRKF. During simulation 6, the SSTRKF improves by 78.5% and 61.5% the horizontal and vertical positioning precisions compared with the RKF. As for KF, 89.3% and 76.7% horizontal and vertical precision improvement are achieved by conducting the SSRTKF. Moreover, the SSTRKF improves the positioning solutions slightly in the simulation 1. As shown in Figure 6, the observation condition of simulation 1 is better than the other five simulations, contributing to well-independent positioning. Hence, all three schemes can properly identify and handle the gross-error-contaminated satellites. According to the analysis above, the proposed SSTRKF evidently outperforms to the KF and RKF during the urban environment simulations, especially in the severe blockage environment.

Conclusions
The standard Kalman filter is optimal for PPP/INS tightly coupled only when the observation noise obeys the Gaussian assumption. However, the LOS multipath and NLOS reception gross errors in the urban environment lead to heavy tails in the observation noise. In the urban severe blockage environment, the heavy tails cannot be properly fitted by the IGG-III function. Hence, in this paper, we apply the Student's t distribution to fit the heavy tails of the gross-error-contaminated observation noise. Moreover, the SSTRKF is proposed by sequentially conducting the RKF and the sequential Student's tbased Kalman filter. The research findings and conclusions are summarized as follows.
1. The GNSS phase and code observation noise obey the Gaussian assumption in the absence of the LOS multipath and NLOS reception errors. Moreover, the Student's t distribution can fit the heavy tails of the gross-error-contaminated observation noise. 2. The proposed SSTRKF can adjust the IGG-III function-derived unreasonable variances through the chi-square test and the sequential Student's t-based Kalman filter, respectively. 3. The numerical comparisons have validated our proposed SSTRKF for the gross-error-contaminated observations. Compared with the RKF, the proposed SSTRKF improves the horizontal and vertical positioning precisions by 57.5% and 62.0% on average during the urban environment simulations. Consequently, the proposed SSTRKF is superior to the KF and the RKF in the urban environment.
Author Contributions: Figure 10. RMS positioning errors of the three schemes during the six simulations.

Conclusions
The standard Kalman filter is optimal for PPP/INS tightly coupled only when the observation noise obeys the Gaussian assumption. However, the LOS multipath and NLOS reception gross errors in the urban environment lead to heavy tails in the observation noise. In the urban severe blockage environment, the heavy tails cannot be properly fitted by the IGG-III function. Hence, in this paper, we apply the Student's t distribution to fit the heavy tails of the gross-error-contaminated observation noise. Moreover, the SSTRKF is proposed by sequentially conducting the RKF and the sequential Student's t-based Kalman filter. The research findings and conclusions are summarized as follows.

1.
The GNSS phase and code observation noise obey the Gaussian assumption in the absence of the LOS multipath and NLOS reception errors. Moreover, the Student's t distribution can fit the heavy tails of the gross-error-contaminated observation noise.

2.
The proposed SSTRKF can adjust the IGG-III function-derived unreasonable variances through the chi-square test and the sequential Student's t-based Kalman filter, respectively. 3.
The numerical comparisons have validated our proposed SSTRKF for the gross-errorcontaminated observations. Compared with the RKF, the proposed SSTRKF improves the horizontal and vertical positioning precisions by 57.5% and 62.0% on average during the urban environment simulations. Consequently, the proposed SSTRKF is superior to the KF and the RKF in the urban environment.