Novel Process Noise Model for GNSS Kalman Filter Based on Sensitivity Analysis of Covariance with Poor Satellite Geometry

One of the great unsolved GNSS problems is inaccuracy in urban canyons due to Non-Line-Of-Sight (NLOS) signal reception. Owing to several studies about the NLOS signal rejection method, almost all NLOS signals can be excluded from the calculation of the position. However, such precise NLOS rejection would make satellite geometry poor, especially in dense urban environments. This paper points out, through numerical simulations and theoretical analysis, that poor satellite geometry leads to unintentional performance degradation of the Kalman filter with a conventional technique to prevent filter divergence. The conventional technique is to bump up process noise covariance, and causes unnecessary inflation of estimation-error covariance when satellite geometry is poor. We propose a novel choice of process noise covariance based on satellite geometry that can reduce such unnecessary inflation. Numerical and experimental results demonstrate that performance improvement can be achieved by the choice of process noise covariance even for a poor satellite geometry.


Introduction
Global Navigation Satellite System (GNSS) is an absolute positioning sensor widely used in various applications. GNSS can determine the location of a receiver based on trilateration by using pseudo-ranges between GNSS satellites and a receiver. It is wellknown that positioning performance is dramatically degraded because of Non-Line-Of-Sight (NLOS) signal reception in urban environments [1]. NLOS signal reception occurs where tall buildings block Line-Of-Sight (LOS) direct signals from satellites and the signals are received via reflection and/or diffraction. Pseudo-ranges obtained from NLOS signals have tendency to be positive outliers [2], because they are computed from the times of flight of the signals. We, therefore, have to reject NLOS signals from position calculation to improve positioning accuracy in urban environments.
Many researchers have developed NLOS signal rejection methods based on measurement residuals [3][4][5], fisheye-camera images [6,7], 3D building maps with ray-tracing [8,9], and so on. It is possible to reject NLOS signals precisely by applying these methods, and then positioning performance can be improved if a sufficient number of satellites are available. Unfortunately almost all signals in dense urban environments could be NLOS signals. The precise rejection of NLOS signals would greatly reduce the number of satellites available in position calculation, thus cause poor satellite geometry, that is, a biased distribution of visible satellites. Since accuracy of trilateration is highly dependent on the geometric locations of the satellites as seen by a receiver, poor satellite geometry may result in inaccurate positioning. The above mentioned studies have not sufficiently dealt with the problem of poor satellite geometry caused by NLOS rejection. Taking satellite geometry into account in addition to NLOS rejection might yield more accurate positioning.
The estimate of position usually can be obtained by using Kalman filter [10,11] under some assumptions about process model, measurement model, and noise statistics. However, the estimation by Kalman Filter may suffer from its divergence [12,13]. A typical example of filter divergence occurs when the assumptions are inconsistent with actual model. The estimation-error covariance becomes excessively small due to the model errors, the filter gain is therefore small, and eventually the filter begins to ignore new measurements. As a consequence, estimate errors will be biased or diverge. Since these model errors are in general unknown and unbounded, it is difficult to solve the divergence theoretically.
A conventional and simple technique [14] to prevent filter divergence is to bump up process noise covariance to cover the uncertainty of model errors. This technique is easy to implement, that is, we have only to add fictitious process noise covariance to its nominal value. Since this technique is largely heuristic, it is important to discuss its effect not only on divergence but also on filter performance.
In the application of the technique to GNSS positioning [4], bumping up process noise covariance inflates estimation-error covariance in the prediction step. The inflated covariance can be reduced by processing observations right after the inflation, if satellite geometry is good. Consequently, the size of estimation-error covariance is kept to be large enough to avoid the divergence, but not too large. However, if satellite geometry is poor, what would be happened? The inflation cannot be reduced in a certain direction according to the geometry, and the estimation-error covariance would become larger than expected. This unexpectedly large covariance would lead to inaccurate positioning.
A proper choice of the process noise covariance has been of interest to many researchers so far. Some filtering techniques have been developed to have an estimation-error covariance that is insensitive to changes in process noise covariance, through the steadystate analysis of estimation-error covariance [15,16]. Furthermore, also noise-adaptive estimation techniques have been proposed [17][18][19]. However, these techniques cannot handle with preventing the performance degradation caused by poor satellite geometry, because they do not take satellite geometry into account.
There are some approaches that can be used to avoid the performance degradation caused by poor satellite geometry. Constrained Kalman filter techniques can be utilized to mitigate the degradation [20,21]. The techniques automatically detect poor satellite geometry, and estimate only the state variables for which poor satellite geometry would not cause the degradation. However, it is difficult to choose a proper threshold for detecting poor satellite geometry in a heuristic manner. The positioning degradation can also be resolved by integrating other sensors such as IMU (Inertial Measurement Unit) and LiDAR (Light Detection and Ranging) device into GNSS Kalman filter [4]. The observations from the sensors can compensate for the lack of information due to poor satellite geometry. However, adding those sensors is often unacceptable in terms of cost and installation space.
The main contribution of this paper is twofold. First, we clarify the fundamental relationship between process noise and satellite geometry through theoretical analysis. The analysis shows that the inflation of estimation-error covariance by fictitious process noise can be unexpectedly large for poor satellite geometry. Secondly, from the theoretical results, we propose a novel way to choose process noise covariance based on satellite geometry in order to remove the unexpected inflation. Comparing with the approaches mentioned above, our approach is based on a theoretical sensitivity analysis, and does not require other additional sensors.
This paper is constructed as follows. Section 2 illustrates the problem of poor satellite geometry by using motivating examples. We describe theoretical analysis for Kalman filter to show the relationship between process noise covariance and satellite geometry in Section 3. We then introduce the extended Kalman filter for GNSS stationary positioning and derive a novel way to choose process noise covariance based on satellite geometry in Section 4. We present simulation and experimental results in Sections 5 and 6. Section 7 concludes this paper.

Poor Satellite Geometry Caused by NLOS Signal Rejection
This section illustrates the problem caused by poor satellite geometry. Figure 1 is an example of a satellite sky plot superimposed on a fisheye-camera image taken at Kobe city, Japan. As shown in Figure 1, tall buildings surround the location and can block almost all direct signals from satellites. We suppose that satellites overlapped on buildings are NLOS satellites (denoted as red dots), and the others are LOS satellites (denoted as blue dots). If complete NLOS satellite rejection is done, only the five satellites will remain. We can see that remained satellites are located along the northeast direction. If we bump up process noise covariance to avoid the filter divergence, position uncertainty in the direction of southeast will become larger unintentionally. Although it will be explained from a theoretical point of view in Section 3, we will show it through a simulation study in this section.

Motivating Examples for Poor Satellite Geometry
In order to show the problem clearly, we present some simulation results using a simplified model for GNSS, under the assumption that five satellites are stationary and located in poor satellite geometry. A GNSS receiver is also assumed to remain stationary on a horizontal surface. Let the dynamical system and the linear measurement equations be as follows: where x k = (r e , r n , r u ) T is a position vector represented in the ENU (East-North-Up) coordinates, y k is a measurement vector composed of ranges from five satellites to a receiver, H is a 5 × 3 measurement matrix, v k is a measurement noise vector, and w k is a process noise vector. We assume that x 0 ∼ N(0, I), E{w k } = E{v k } = 0, E{w k w T l } = Q(k)δ kl = 0.01Iδ kl (since the position vector x k is stationary, w k and Q(k) should be zero originally. However, such a choice of w k and Q(k) would lead to numerical problems of Kalman filter [11]. To avoid the problems, small w k and Q(k) are usually assumed to exist. In this paper, we suppose that "true process noise" denoted by w k and Q(k) includes this inevitable small noise) and E{a} is the expectation of a random variable a, and δ kl is the Kronecker delta. w k and v k are produced numerically by random number. θ i and ψ i are the elevation angle and the azimuth angle of the ith satellite, respectively. The five satellites are supposed to be located as shown in Figure 2. Their elevation angles are θ = {90, 15, 15, 15, 15} [degree], and their azimuth angles are ψ = {0, 40, 50, 220, 230} [degree]. The satellites are located along the northeast direction as in Figure 1. We can determine whether the satellite geometry (it is well-known that DOP (Dilution of Precision) parameters characterize satellite geometry. In this paper, we attempt to integrate satellite geometry into Kalman filter) is poor or not by evaluating the eigenvalues of H T H; If one of eigenvalues is zero or close to zero, H T H is singular or almost singular, and such a satellite geometry is called poor. The eigenvectors of H T H corresponding to the minimum and maximum eigenvalues coincide with a southeast direction and a northeast direction in this satellite geometry (the locations of the five satellites were chosen so that the directions of the eigenvectors of H T H are intuitively obvious. Although this paper presents the results when the elevation angles θ i for four satellites are 15 degrees (i = 2, . . . 5), similar numerical results are obtained for θ i = 0, 30, or 45 degrees). Simulation results are obtained by Kalman filter that will be described in Section 3, with two parameter settings. For setting (1), the process noise covarianceQ in Kalman filter is chosen as the true one Q. For setting (2),Q is set by adding fictitious noise δqI to Q. Figure 3 shows the estimated positions in a r e -r n plane that were obtained by repeating the filter computation 1000 times with different initial position estimates and different noise realizations. Table 1 shows standard deviations that are computed by projecting the plotted estimate errors to two eigenvectors of H T H, where σ SE and σ NE are standard deviations in the directions of two eigenvectors, southeast direction and northeast direction. The differences between σ SE and σ NE mainly comes from the satellite geometry in both (1) and (2). Since there is no model error in these simulations, the plotted estimation errors spread depending on the satellite geometry and the settings of process and measurement noises, and are consistent with the estimation-error covariance calculated from Kalman filter. However, we should note that the values of σ SE /σ NE in (1) and (2) are quite different-3.48 and 4.22, respectively. This means that the estimation errors in southeast direction are spread more widely than in northeast direction even though the fictitious process noise is added equally in each direction. The fictitious noise to avoid the filter divergence makes the estimation errors worse unintentionally in the southeast direction that corresponds to the eigenvector of the smallest eigenvalue.

Settings
(1) We devote Section 3 to reveal the theoretical relationship among the inflation of estimationerror covariance, process noise covariance, and satellite geometry. Then, in Section 4, we will describe how to reduce the unnecessary inflation of estimation-error covariance.

Sensitivity Analysis of Estimation-Error Covariance
We consider a linear system and Kalman filter for it: where we use the same notations as in Equations (1) and (2), and variables in the filter are listed in Table 2. We suppose thatR(k) andQ(k) are different from their true ones R(k) and Q(k).
First of all, we derive the difference equation between P(k − 1|k − 1) and P(k|k). Since is a diagonal matrix whose diagonal elements are the eigenvalues of H T (k)R −1 (k)H(k). Now we make the following assumptions: Variable Definition Posterior estimate P(k|k − 1) Prior estimation-error covaraiance P(k|k) Posterior estimation-error covaraiance Filter gain matrix k Time step We drop the notation (k) for simplicity except for P(k − 1|k − 1) and P(k|k). Using Assumption 1 and substituting H TR−1 H = GΛG T into Equation (8), we can rewrite it as Multiplying G T and G to both sides of Equation (10), we obtain From Equation (11) and Assumption 2, G T P(k|k)G is also diagonal. Then the iith element of G T P(k|k)G denoted asP i,k can be written as whereλ i is the iith element ofΛ, andq i is the iith element of G TQ G. Equation (12) is the difference equation of estimation-error covariance from k − 1 to k, and shows thatP i,k depends onq i andλ i , that is, process noise covariance and satellite geometry. By using the true process noise covariance Q and a fictitious process noise δq, we representQ asQ = Q + δq. It is assumed that each component inQ is diagonalized by G, that is, G T QG and G T δqG are diagonal matrices. Then, Equation (12) can be rewritten as: where q i and δq i are the iith element of G T QG and G T δqG, respectively. Now we define the variation ofP i,k due to δq i as ∆P i . From Equation (13), ∆P i is approximated as the following equation.
From Equation (14), we examine the sensitivity of ∆P i with respect toλ i . Since H TR−1 H is a positive semidefinite matrix,λ i ≥ 0 and That is, ∆P i (λ i ) decreases monotonically asλ i goes from zero to infinity. At the limits For good satellite geometry, there is no eigenvalueλ i close to zero. The inflations of P i,k that are caused by the fictitious noise δq i between time steps k − 1 and k are suppressed for all i by the observation at time step k, according to Equation (14). By choosing an appropriate value of δq i through trial and error, we can keep the covariance P(k|k) reasonably large to avoid the filter divergence even for the system with model errors. On the other hand, at least one of eigenvalues is zero or close to zero for poor satellite geometry. If we choose the size of δq i for good satellite geometry, it may be too large for poor satellite geometry and cause unintentionally large inflation ofP i,k . Too large covarianceP i,k leads to widely distributed estimation errors along the corresponding direction and degrades the filter performance.
In the next section, we will propose a new way to choose the fictitious process noise δq i based on satellite geometry H T (k)H(k) to reduce the inflation. To derive the choice, we make an additional assumption: Assumption 3.R(k) = r k I, where r k is a positive scalar.

From Assumption 3, H T (k)R −1 (k)H(k) can be rewritten as
where Λ(k) is diagonal. From Equations (13) and (14), we obtaiñ where λ i is the iith element of Λ(k). It should be noted that λ i is an eigenvalue of H T (k)H(k) and, unlikeλ i , represents the satellite geometry independently ofR(k). Moreover, we can compute a steady-state solution of Equation (12) (or Equation (16)), under the assumptions thatR(k), H(k) andQ(k) are time-invariant, and G T P(0|0)G is diagonal. Denoting the steady values ofP i,k andP i,k−1 asP i and substituting them into Equation (12), we can obtain a steady-state solution: The solution in Equation (18) is equivalent to the solution of an algebraic Riccati equation of estimation-error covariance [22]. It should be noted that the solution in Equation (18) is consistent with the results in Section 2.2. The values ofP i that are calculated withλ i andq i for the motivating example coincide with the standard deviations in Table 1. Even though a steady-state solution is well-known, the analysis in this section revealed how the covarianceP i,k is inflated at each time step depending on the fictitious noise δq i and the satellite geometry λ i (orλ i ). By using the analytical result in Equation (17), we will propose a novel process noise model in the next section to avoid an unintentionally large inflation at each time step.

Extended Kalman Filter
In this section, we describe the extended Kalman filter for GNSS stationary positioning. Since the pseudo-range and Doppler frequency measurement equations are non-linear, the extended Kalman filter (referred as EKF) is generally used. We define the state vector as is an inter-system bias (ISB) vector between GPS time and other GNSS's time. We suppose to use a multi-GNSS single frequency receiver, for example, FURUNO GN-8720, and t s consists of three elements that are GPS-QZSS, GPS-GLONASS, and GPS-Galileo ISBs.
The dynamical system and the measurement equations are given as follows: where h (j) (x k ) and h (j+M) (x k ) are the jth and the j + Mth elements of h(x k ), r (j) (k) andṙ (j) (k) are the position and velocity vectors of the jth satellite, and t m is the ISB of the GNSS system to which the jth satellite belongs. M is the number of satellites. The EKF for this system consists of the computations outlined from Equation (5) to Equation (9), but with Equation (7) replaced bŷ and H(k) computed as The EKF may often produce a biased estimate and diverge due to unmodeled measurement errors, dynamics errors, and non-linear effects of Equation (21) [23]. To cover the uncertainties by a fictitious process noise, we usually choose the process noise covariance in the EKF asQ For the EKF in this paper, we set Q to 0.01I to avoid the numerical problems as in Section 2.2. ThisQ makes gain matrix K(k) large from Equations (6) and (9), thus the EKF updates the estimatex by overweighting current observations. It should be noted that a largerR cannot be effective in preventing filter divergence, since a largerR makes K(k) smaller. We chooseR so that it represents actual measurement noise level. In this paper,R is determined based on the Signal-Noise-Ratios (SNRs) of the receiver FURUNO GN-8720 as follows [2]:R = diag(σ 2 1 (s 1 ), . . . , σ 2 M (s M ), µ 2 1 (s 1 ), . . . , µ 2 M (s M )), where s k , σ k and µ k are the SNR, the standard deviations of pseudo-range noise and Doppler frequency noise for the kth satellite signal.

Process Noise Covariance Based on Satellite Geometry
In practical applications of EKF, it is not easy to find an appropriate choice of the fictitious noise δq in Equation (25). It is often chosen and fixed through a trial and error process, so that a required performance of EKF is obtained. If the satellite geometry is always good and all the eigenvalues λ i are equally large, such a choice of δq could work well. In this paper, we call it the conventional choice, and consider only δq in the following form: where δq is a positive scalar. However, if the satellite geometry is poor, the conventional choice of δq would cause an unintentional inflation of P as described in Sections 2 and 3.
In this section, we propose a novel way to choose the fictitious noise δq based on satellite geometry by assuming that Equation (16) approximately holds in our EKF. The basic idea is to choose δq so that ∆P i in Equation (17) is less than or equal to a specified value c i at each time step k in order to suppress the unintentional inflation. The fictitious noise is chosen as follows.
where the matrix Q is diagonal, and Q i is the ith diagonal element of Q, and δq is a positive scalar for the conventional choice. From Equations (17) and (28), if δq i ≤ δq, ∆P i satisfies ∆P i ≈ c i even for a small λ i . If δq i > δq, ∆P i < c i holds approximately even for δq i = δq. Consequently, we can expect that the unintentional inflation of P is suppressed, because ∆P i is limited to the maximum value of c i , regardless of satellite geometry λ i . For the proposed choice of δq, we can keep an appropriate inflation of P by the parameters c i . We will determine c i through simulation study in Section 5. Moreover, the steady-state solutionP i in Equation (18) does not hold for δq in Equation (28), becauseq i in Equation (12) includesP i,k−1 . By substituting Equation (28) into Equation (16), we can obtain the following equation: By solving the above equation numerically, we can obtain an approximate steady-state solutionP i .
Strictly speaking, Equation (16) does not hold for the EKF in this section, because the three assumptions in Section 3 may be unsatisfied, that is to say, (1) Φ(k) = I, (2) G T (k)P(k − 1|k − 1)G(k) is not diagonal, and (3)R = r k I. For Assumption 3,we choosẽ R(k) as diagonal and give its diagonal elements based on SNRs in this paper. Since SNRs of LOS signals usually have similar values, we can expect thatR(k) ≈ r k I and Assumption 3 is approximately satisfied. For Assumptions 1 and 2, we cannot justify them theoretically. Although Assumption 1 is necessary to show that both P(k|k) and P(k − 1|k − 1) are diagonalized by G(k), they are almost diagonalized for the experimental results shown in Section 6. Therefore, Assumption 2 is approximately satisfied in the experimental results, and then we do not need Assumption 1 to show it. However, for Φ = I, Equation (11) itself does not hold. In this paper, we suppose that the effects of Φ with F in Equation (23) are small, because only one nondiagonal element is nonzero. We also note that, in Assumption 2,Q can be chosen with δq in Equation (28) such that G T (k)QG(k) is diagonal. As a result, we believe that Equation (16) approximately holds in our EKF.
We finally summarize the proposed noise model. Employing the fictitious noise δq in Equation (28),Q varies depending on satellite geometry, and ∆P i can be insensitive to satellite geometry. We, therefore, call the model satellite-geometry adaptive model. Although a noise-adaptive model introduced in Section 1 is well known,Q is estimated based on recent observations, and an accurate estimation would be difficult for poor satellite geometry. To our best knowledge, our model is the first process noise model that changes with H T H, i.e., satellite geometry.

Simulation Study
We present here some simulation results of the proposed process noise model to demonstrate its effectiveness and to determine an appropriate value of c i . The model is applied to the Kalman filter for the simplified system in Section 2.2, where all the assumptions in Section 3 are satisfied and Equation (16) holds.
We consider the following three settings of the process noiseQ. Setting (1): No fictitious process noise:Q = Q (30) Setting (2): Conventional choice of fictitious process noise: Setting (3): Proposed choice of fictitious process noise: The first and second settings are the same as the settings in Figure 3a,b, where δq is chosen as δq = 1 in Equation (31). In this paper, we assume that, for Equation (32), all the c i are the same, that is, c i = c for ∀i. Numerical simulations for setting (3) were performed with different values of c: (3-a) with c = 1.0 to (3-f) with c = 0.01 as shown in Table 3. The standard deviations σ SE , σ NE , and their ratios for all the settings are summarized in Table 3, where the results for settings (1) and (2) in Table 1 are shown again. The Root Mean Squared Error (RMSE) d is also computed for the three-dimensional position to show the positioning accuracy in the table.
For determining an appropriate value of c, we should pay attention to two points: avoiding the filter divergence and suppressing the unintentional inflation of P along the eigenvectors corresponding to small eigenvalues λ i . For the first point, σ NE , that is the standard deviation along the direction with a large eigenvalue, is increased to 0.80 in setting (2), comparing with setting (1). A similar magnitude of σ NE would be necessary even in setting (3) for the divergence avoidance. For the second point, we check the ratio of σ SE to σ NE . In setting (1), the ratio comes only from the satellite geometry, because there is no model error or fictitious noise. Therefore, we suppose that the ratio in setting (1) is desirable even in setting (3) for suppressing the inflation. From these two points, we can choose setting (3-c) c = 0.36 for the proposed noise covariance model in this paper. Moreover, it should be noted that the values of standard deviations in (3-a) to (3-f) are consistent with the steady-state solutions computed from Equation (29). Table 3. Standard deviations, their ratios, and RMSEs with different process noise settings. The effectiveness of the proposed noise model is verified by comparing the simulation results in settings (2) and (3-c). The estimated positions in a r e -r n plane for the two settings are plotted in Figure 4a,b, where the filter computation was repeated 1000 times for each setting as in Section 2.2 and Figure 4a is the same as Figure 3b. Ellipsoids in these figures are the 1σ contours of probability density function (Gaussian) that are calculated from Equation (18) for Figure 4a and from Equation (29) for Figure 4b. The distribution of horizontal position errors obtained by numerical simulations coincides with the ellipsoid obtained theoretically for each setting. We can see that the estimation errors in setting (3-c) are reduced in the direction of southeast, that is, the direction of the minimum eigenvalue of H T H. In this numerical results, the proposed noise model suppresses the unintentional inflation of σ SE by about 20%. It should be noted that these results show that both precision and accuracy of positioning in the horizontal plane are improved in setting (3-c), because the mean value of estimation errors is almost zero. We can also check the positioning accuracy that includes the vertical errors, by comparing the RMSEs for settings (2) and (3-c) in Table 3.

Experimental Results of Stationary Positioning
This section presents experimental results to show the effectiveness of the proposed process noise model by applying it to actual data obtained in an urban canyon. The data was collected at the time and location shown in Table 4, and Figure 5 is a fisheye-camera image with a satellite sky plot at the beginning of data collection. We used a multi-GNSS single frequency receiver, FURUNO GN-8720, that can receive the signals of GPS, QZSS L1C/A, GLONASS L1OF, and Galileo E1, and acquired about ten minutes the data that includes pseudo-ranges, Doppler frequencies, navigation messages, and so on. The receiver was fixed at the location that corresponds to the center of Figure 5 during the experiment, and the pseudo-ranges and Doppler frequencies were used as the measurements for EKF. In the figure, blue dots are LOS satellites, and red dots are NLOS satellites, as in Section 2.1. By detecting NLOS satellites from the fisheye-camera images, the signals from them were excluded from the measurements. The LOS satellites are quite close to each other, and the satellite geometry can be considered as poor, because the minimum eigenvalue of the LOS satellite geometry is always about 0.01 in the data (although multipath signals may occur via diffraction/reflection even for LOS satellites, they would be largely attenuated compared to LOS direct signals. In this paper, we do not consider the effects of multipath signals by assuming that they are sufficiently small). Based on the numerical results shown in Section 5, we examine positioning performance of the EKF in Section 4.1 with two parameter settings for process noise covariance: Setting (2) with δq = 1 and Setting (3-c) with c = 0.36. The EKFs with the conventional and the proposed noise settings were implemented in a laptop computer as post-processing programs, by customizing the EKF program implemented in FURUNO GN-8720. The EKFs compute the estimated position of the receiver by using only the measurements from LOS satellites.
In order to obtain an ensemble average of estimation errors, we make N measurement data sets whose length is l, from the data acquired for about ten minutes, by shifting the start point by one step. The EKF calculation is performed for each of the N data sets, Y 1 = {y 1 , y 2 , · · · , y l }, Y 2 = {y 2 , y 3 , · · · , y l+1 }, · · · , Y N = {y N , y N+1 , · · · , y l+N−1 }, with the initial conditions shown in Table 4. The length l of each data set is chosen such that the changes inx and P become sufficiently small in the EKF as the time step k approaches l. We focus on the estimation error in position at the final step,r(l|l) − r, where the true position r is known as in Table 4 , and denote it asr m for mth data set.

Configurations and Parameters Remarks
The minimum and maximum standard deviations σ min and σ max of the estimation errorsr m are obtained by calculating the minimum and maximum eigenvalues of Σ. Table 5 summarizes the values of σ min , σ max , σ max /σ min and d for the two settings. The following two points can be seen in Table 5: (a) σ max and σ max /σ min in setting (3-c) are less than the ones in setting (2), while the values of σ min in both settings are almost the same, and (b) RMSE d in setting (3-c) is also less than the one in setting (2). The first point (a) would indicate that the proposed noise model can suppress the unintentional inflation of estimation-error covariance even for actual data in an urban canyon. From the second point (b), the positioning accuracy with the proposed noise model is also improved compared to the conventional noise model. The estimation errorsr m for the two settings are shown in Figure 6a,b, where the errorsr m are projected to the plane spanned by the unit eigenvectors e min and e max that correspond to σ min and σ max , respectively. The ellipsoid in each figure is the 1σ contour of probability density function that is supposed to be Gaussian with the standard deviations σ min and σ max in Table 5.  Although the results in this section demonstrate that the suppression of covariance inflation is achieved by the proposed noise model, it should be noted that the eigenvector corresponding to σ max is almost along the altitude direction. The standard deviations σ min and σ max in Table 5 are also quite different from the ones calculated from Equations (18) and (29). These differences from the simulation results in Section 5 may be due to the following two reasons. First, the satellite geometry in Figure 5 is largely different from the one for the simulation results, because there is no LOS satellite with a low elevation angle in Figure 5. The uncertainty in the altitude direction would be large due to the clock bias t b for such a satellite geometry (if the elevation angles of all the satellites are close to 90 degrees, the errors in altitude estimation and clock bias estimation would be indistinguishable. Through numerical simulations based on the simplified simulation model introduced in Section 2.2 for the same satellite geometry as in Figure 5, we can see that the uncertainty in the altitude direction with clock bias estimation is much larger than the one without it). Secondly, Assumption 1 is not satisfied in the EKF in Section 4.1. Equations (18) and (29) do not hold without the assumption. The non-diagonal component of F in Equation (23) would be non-negligible, especially because the position estimation in the altitude direction is highly dependent on the clock bias t b .
The above inconsistency between the simulation results and the experimental results could be resolved by extending the proposed process noise model to the systems without the assumptions in Section 3. Our preliminary results indicate that theoretical extension of the sensitivity analysis is possible, and it will allow us to develop a process noise model in a more consistent manner for the systems without the assumptions. Further improvement of the estimation performance may also be possible by changing the choice of c or δq itself. The extension and improvement will be presented in our next paper.

Conclusions
This paper pointed out that a fictitious process noise to prevent the filter divergence can degrade the filter performance for a poor satellite geometry due to an unintentional inflation of estimation-error covariance. A sensitivity analysis of estimation-error covariance by fictitious process noise and satellite geometry was performed under some assumptions, and we proposed a novel model of the fictitious process noise based on satellite geometry in order to suppress the unintentional inflation of estimation-error covariance. The effectiveness of the proposed noise model was shown via simulation and experimental results. Although the sensitivity analysis and the proposed model derived from it are based on the assumptions, the approach in this paper can be extended to the systems without the assumptions. The extension and further improvement of the proposed approach will be presented in our future work.

Data Availability Statement:
The data presented in this study are partially available on request from the corresponding author.

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

Abbreviations
The following abbreviations are used in this manuscript: GNSS Global Navigation Satellite System NLOS Non-Line-Of-Sight LOS Line-Of-Sight EKF Extended Kalman Filter SNR Signal to Noise Ratio RMSE Root Mean Squared Error