Multi-Epoch Kriging-Based 3D Mapping-Aided GNSS and Doppler Measurement Fusion using Factor Graph Optimization

.

several researchers have proposed methods to identify, mitigate, and correct such errors to improve the GNSS performance in urban canyons.
With the digitalization of smart cities, open-source three-dimensional (3D) building models are widely available for use, such as in Hong Kong (Hong Kong Lands Department, n.d.) and New York City (Department of City Planning).Such models can aid in urban positioning by predicting the visibility of satellites and the transmission path of GNSS signals with reflection correction.This type of assisted positioning is known as 3D mapping-aided (3DMA) GNSS positioning (Groves, 2016).3DMA GNSS positioning can be used in a low-cost receiver as a software-based approach to aid in positioning.In a previous study, a 3D building model was utilized to exclude NLOS satellites, and the weighted least-squares (WLS) solution was enhanced (Obst et al., 2012).However, the actual location of the receiver was unknown, and the 3D building model could not implement the actual correction.
To address this problem, a hypothesis position candidate-based approach has been proposed to identify the optimal location with the highest similarity between the predicted and observed measurements.The existing 3DMA GNSS algorithms can be divided into shadow matching (Wang et al., 2013(Wang et al., , 2015) ) and ranging-based 3DMA GNSS, which includes ray-tracing GNSS (Hsu et al., 2016;Miura et al., 2015), likelihood-based ranging 3DMA GNSS (Groves et al., 2020), and skymask 3DMA GNSS (Ng, Zhang, & Hsu, 2020).In addition to studies on using a candidate-based approach for position estimation, a set-based shadow-matching approach based on zonotopes has been proposed to mitigate the multimodal posterior distribution due to symmetries in building geometry (Bhamidipati et al., 2021;Neamati et al., 2022).A study on integrating the shadow-matching concept with direct positioning estimation has been conducted to improve positioning accuracy (Strandjord et al., 2020a(Strandjord et al., , 2020b)).Another study proposed the combination of shadow matching and ray-tracing and the development of a secularity metric to simultaneously predict the likelihood of building reflections and improve positioning (Strandjord, Axelrad, & Mohiuddin, 2020).
The shadow-matching technique uses a 3D building model, satellite ephemeris, and carrier-to-noise ratio (C/N 0 ) values of received measurements to identify the candidate for which the satellite visibility prediction best matches the received satellite data as the receiver location.The C/N 0 value determines the LOS probability for a particular received satellite to mitigate the NLOS reception effect.The observation of any other non-received satellite that appears in the ephemeris data is assumed to correspond to an NLOS satellite.An intelligent classifier (Sun et al., 2019) has been used to identify NLOS-received satellites in shadow matching (Ng, Zhang, & Hsu, 2021).Ranging-based 3DMA GNSS follows a similar principle; however, the matching basis corresponds to the pseudorange.Such approaches model the pseudorange measurements, process the degraded measurement with NLOS correction for each candidate, and compare the values with the observed pseudorange.The candidate with the highest similarity between the modeled and observed pseudoranges is selected as the receiver location.The NLOS correction implemented through the ranging-based 3DMA GNSS may be statistical or geometrical.Likelihood-based ranging 3DMA GNSS uses a statistical approach to remap the NLOS pseudorange error to LOS measurements (Zhong & Groves, 2022).In contrast, ray-tracing GNSS (Hsu et al., 2016;Miura et al., 2015) and skymask 3DMA GNSS (Ng, Zhang, & Hsu, 2020) attempt to identify the reflection point and estimate the reflection delay.The former algorithm validates the reflection over every possible reflector and introduces a high computational load (Hsu et al., 2016;Miura et al., 2015), whereas the latter reduces the computation cost by identifying the reflecting point using an enhanced skymask (Ng, Zhang, & Hsu, 2020), i.e., a skyplot in which building boundaries are combined with building height information.
The two types of 3DMA GNSS positioning algorithms exhibit different positioning performances under different building distributions surrounding the receiver.Because building boundaries are more distinctive in the lateral street direction, shadow matching is excellent for determining the actual side of the street on which the receiver is located.In contrast, ranging-based 3DMA GNSS exhibits high performance in the longitudinal street direction, because the satellite distribution geometry is better along the direction of the street.Considering this complementarity, it is valuable to integrate the 3DMA GNSS techniques to achieve satisfactory performance in all directions (Groves et al., 2020).Another approach to aid in positioning is to use L5 signals in 3DMA GNSS (Ng, Zhang, Luo, & Hsu, 2021).L5-band GNSS signals can mitigate the error caused by signal reflection (Ng, Zhang, Yang et al., 2020).A narrow correlation peak of L5 signals can help eliminate long-delay multipath errors exceeding one chip length.In particular, L5 signals can help mitigate multipath effects and isolate NLOS reception.Therefore, a precise NLOS biased measurement can help enhance the performance of ranging-based 3DMA GNSS positioning (Ng, Zhang, Luo, & Hsu, 2021).
The grid distribution, based on evenly spaced sampling of positioning hypothesis locations, is sometimes flawed, which can lead to two problems: 1) multimodal, local minimum, or drifting issues associated with candidates having a likelihood similar to that of the actual location owing to the high similarity in building geometry, as shown in Figure 1, and 2) a high computational load associated with excessively distributed candidates, especially for candidates far from the actual location.For candidate-based 3DMA GNSS, the required determination time of a single-epoch solution is based on the number of candidates and satellites.For ray-tracing 3DMA GNSS, the number of building surfaces in the model influences the computational load.The combination of these factors corresponds to a massive computation load, which renders the practical and real-time application of 3DMA GNSS challenging.Researchers have attempted to decrease the computational load by using an enhanced skymask: skymask 3DMA GNSS (Ng, Zhang, & Hsu, 2020) determines the reflection point by considering the azimuth angle of the reflecting planes instead of propagating the transmission path over different building surfaces, as in ray-tracing 3DMA GNSS.Moreover, in the Gnome approach, the heavy computational load processes were performed through a cloud-based framework, and the pre-computed information was stored in the device for real-time usage (Liu et al., 2018).Another study accelerated the ray-tracing GNSS method for prac-tical implementation (Ziedan, 2017).Researchers have proposed the improved multipath and NLOS signals identification (iMNSI) algorithm to find an intersect-ing point that is shared between several planes and to discard those satellites whose planes do not intersect at the shared point (Ziedan, 2018).
The objective of the candidate-based 3DMA GNSS method is to identify a position that can minimize the residual between the modeled and observed measurements.The determined position should have a maximum likelihood score.
In other words, one can obtain the 3DMA GNSS solution iteratively instead of using a candidate-based approach, by expressing the likelihood of the particle with a mathematical model.Previously (Ng, Zhang, Wen, & Hsu, 2021), our group attempted to mathematically express the likelihood model of shadow matching based on a set of trigonometric functions and resolved the receiver state by using a nonlinear least-squares (NLS) method, specifically, the Gauss-Newton method.However, the results indicated that this approach could not be applied practically because of its dependence on the initial guess of the receiver location.Nevertheless, the findings highlighted the potential of implementing 3DMA GNSS through an optimization-based approach, which corresponds to a relatively low computational load.
The 3DMA GNSS method can also be enhanced through a multi-epoch approach.Most of the recently proposed 3DMA GNSS algorithms adopt a single-epoch snapshot-based approach to estimate the positioning, and thus, the solutions are not robust.Moreover, most applications require continuous navigation, for instance, when pedestrians hold a smartphone while walking along a street.Temporal links between solution epochs can be considered to enhance 3DMA GNSS.Several research groups have considered the use of filtering techniques to implement 3DMA GNSS with an inter-epoch connection.For example, the Kalman filter (KF) and extended KF (EKF) have been used to recursively update the current state through predictions based on past estimation and the error of current measurements (Groves, 2013b;Zhang et al., 2020).A particle filter has been used to effectively distribute and sample the candidates (Suzuki, 2016;Yozevitch & Moshe, 2015).Moreover, a grid filter has been used to evenly distribute positioning candidates (Zhong & Groves, 2022).These studies have demonstrated that the filtering technique can provide a better and smoother positioning performance.Furthermore, machine learning techniques have been used to estimate the most likely paths on a map, corresponding to optimized position estimation (Ziedan, 2021), where the researchers used machine learning to predict the possibility of a signal status change, corresponding to intelligent signal status estimation.The results indicated that such techniques can help to increase the position estimation accuracy.
In addition to the application of filtering techniques to GNSS, factor graph optimization (FGO) has been performed to optimize states with many constraints (Sünderhauf & Protzel, 2012).Unlike other filtering techniques, FGO represents a batch approach that considers all temporal constraints.Therefore, the optimized results are highly robust.Researchers have integrated GNSS positioning with Doppler measurements to provide optimized states for multi-epoch scenarios, and the FGO codes are available as open source (Wen & Hsu, 2021).Furthermore, FGO has been applied to achieve centimeter-level positioning accuracy through carrier-phase measurements based on GNSS precise point positioning (Watson & Gross, 2018) and real-time kinematics (RTK) (Wen & Hsu, 2021).An optimized GNSS based on FGO can be applied to realize 3DMA-GNSS-based collaborative positioning (Zhang et al., 2021).Overall, the integration of FGO can significantly increase the positioning accuracy.
Considering these aspects, this study focuses on 3DMA GNSS and on determining the receiver location more effectively by iteratively applying an NLS technique.As a key contribution, the iterations are implemented by modeling the cost function (which is a function of the residual between the actual measurement and modeled estimation) based on candidate-based 3DMA GNSS.An ordinary Kriging method (Wackernagel, 2003) is used to model the cost function.The Kriging method corresponds to an optimal or best linear unbiased prediction and is widely used in spatial statistics to interpolate data at unsampled locations.The location-related cost function can be determined to calculate the residual by inputting the location of the sampled positioning candidates and sampled cost to the Kriging system.The Kriging-regressed cost function is a differentiable function that can help clarify the Jacobian matrix for iterative NLS to achieve optimal estimation.In sum, the main benefit of using the Kriging method is the ability to include higher-resolution information on sampled candidates and their costs, thus extending the optimization process.As another benefit of modeling the sampled candidates as a cost function, the spacing between sampling positioning hypothesis locations can be relaxed.This method is termed Kriging-based 3DMA GNSS.The second contribution of this study is the establishment of a multi-epoch approach based on Kriging-based 3DMA GNSS.As the cost function in each epoch is differentiable, batch optimization can be performed by simultaneously optimizing all of the receiver states in the current and historical epochs.In this scenario, FGO can be performed.Thus, we integrate the proposed Kriging-based 3DMA GNSS and Doppler measurements through FGO.The Kriging-modeled cost function can provide higher-resolution information during the integration with Doppler measurements to prevent ambiguous positioning when multiple modes exist.As the measurements are expected to contain many outliers during urban positioning, implementing FGO can increase the overall accuracy and smoothness with an acceptable computational load.Data recorded in Hong Kong and New York City by a low-cost commercial-grade receiver and smartphones are used to evaluate the positioning performance of the proposed framework.The results show that the Kriging-based approach and FGO using Kriging-based 3DMA GNSS can achieve a root mean square (RMS) positioning accuracy of approximately 10 m and 6 m, respectively.
The remainder of this paper is organized as follows.Section 2 provides an overview of the proposed 3DMA GNSS method using FGO.Section 3 describes the modeling of the proposed Kriging-based 3DMA GNSS technique.Section 4 details the Doppler measurement factor considered in FGO.Section 5 describes the integration of the Kriging-based 3DMA GNSS method and Doppler factor using FGO.Section 6 presents the experimental results and analysis.Concluding remarks and the scope for future work are summarized in Section 7.

OVERVIEW OF THE PROPOSED SYSTEM
Figure 2 shows an overview of the proposed 3DMA GNSS method using FGO.The system includes a pre-generated skymask database that covers the 3D building models associated with the study area.The 3D building models used in this study were provided by the Hong Kong Lands Department with a level of detail (LoD) of 1 (Hong Kong Lands Department, n.d.) and the Department of Information Technology & Telecommunications' 2014 aerial survey (Department of City Planning) with an LoD of 2 for the testing scenario in Hong Kong and New York City, respectively.3DMA GNSS can be applied only to this coverage area (Ng, Zhang, & Hsu, 2020).After obtaining the pseudorange and Doppler measurements from the receiver, we can estimate the initial position and distribute the candidates.Subsequently, the system can be modeled.All candidates are sampled, and the integrated score for 3DMA GNSS is obtained.Both shadow matching and likelihood-based ranging 3DMA GNSS are applied in the 3DMA GNSS integrated solution.The detailed procedures have been presented elsewhere (Ng, Zhang, Luo, & Hsu, 2021).The determined likelihood scores of candidates can be transformed and normalized to sampled costs, which represent the residuals.Subsequently, we can establish the ordinary Kriging system by inputting the candidates' position and their residuals.The modeled Kriging system and Doppler measurements are input to FGO for batch optimization to establish the FGO using the Kriging-based 3DMA GNSS.For the Kriging-based 3DMA GNSS, only a single-epoch modeled Kriging system must be input for FGO.
The graph problem associated with FGO using the Kriging-based 3DMA GNSS involves two main components: the 3DMA GNSS modeled by a Kriging system and Doppler measurement-related factors.Doppler measurements involve three error factors: Doppler frequency measurements, motion propagation, and motion constraints.These entities are described in Sections 3 and 4, respectively.

INTEGRATED 3DMA GNSS
According to the literature (Groves et al., 2020), 3DMA GNSS evenly distributes the hypothesis position candidates around the initial position.This study also applies context-based sampling for the candidate distribution (Ng et al., 2022).Subsequently, the simulated measurements are generated to be compared with the actual received measurements for each candidate.Theoretically, the receiver's truth position should obtain the highest similarity, and the candidate distribution should cover the true location to maximize performance.Full implementation of the candidate-based 3DMA GNSS can be found in Ng, Zhang, Luo, & Hsu (2021).
Shadow matching compares the received satellite and visibility predictions for each positioning candidate to locate the position with the highest similarity (Wang et al., 2013(Wang et al., , 2015)).Details regarding the implementation of shadow matching can be found in the literature (Groves et al., 2020).After evaluating the LOS probability of the satellite visibility for candidate j, the shadow-matching likelihood score for the candidate, S j SDM , , can be calculated through the geometrical mean of the visibility consistency of all satellites i I } 1 , P j i : The score, calculated by shadow matching, is integrated with the likelihood score specified by ranging-based 3DMA GNSS.
The primary task for ranging-based 3DMA GNSS is to compare the received pseudorange measurements and modeled pseudorange with NLOS correction.In this study, the likelihood-based ranging 3DMA GNSS is used for NLOS correction.The likelihood score for ranging-based 3DMA GNSS, S j LBR , , is calculated using the NLOS-corrected pseudorange differences, , : where Q is a weighting matrix that contains the uncertainty for each satellite, such that ∑ Q is the summation of the diagonal elements in Q.Details regarding the implementation of likelihood-based ranging 3DMA GNSS can be found in the literature (Groves et al., 2020).
The two scores can be integrated to determine the likelihood score for candidate j: Finally, the location with the maximum likelihood score is the estimated receiver location.The candidate-based 3DMA GNSS can be expressed as an NLS: where x is the state, corresponding to the receiver position in earth-centerearth-fixed (ECEF) coordinates.x * is the optimized state, y is the actual observation provided by the receiver, ŷ is the value estimated by a sophisticated model, and δ ˆ( ) x is the cost function.Thus, the problem can be considered an optimization problem.The optimization can be performed iteratively if the cost function is differentiable.Thus, we propose a method to mathematically describe the cost function based on a Kriging method.

Modeling 3DMA GNSS using an Ordinary Kriging Method
Kriging methods are commonly used in spatial statistics when only a limited number of sample points are considered (Matheron, 1963).We use a Kriging method to mathematically model the cost function of the optimization problem described in Equation (4), as illustrated in Figure 3.The 3DMA GNSS method calculates the candidate likelihood score and transforms it into a sampled cost.After the samples are obtained, the candidate locations are used to interpolate and model the cost function δ ˆ( ) x based on an ordinary Kriging method (Wackernagel, 2003).Because the Kriging method regresses the parameters of an assumed function based on the sampled cost data, the function is differentiable.In other words, the Jacobian matrix of the cost function can be determined to perform iterative NLS optimization.
Although several types of Kriging methods are available, the ordinary Kriging method is used because it does not require the mean and covariance function of the sampled data to form the prediction function (Webster & Oliver, 2007).Thus, the formation of our cost function can be simplified.Furthermore, the candidate-based 3DMA GNSS with the Kriging method is modeled to calculate the Jacobian matrix for the NLS and FGO and to predict the unknown values at an unsampled location from the observed data.In this manner, we can relax the separation of the candidates during sampling.
At the position of interest, the ordinary Kriging method can be expressed as follows: where δ ˆ( ) x is the predicted cost at a position state x to be estimated, ω ω is the ordinary Kriging weighting vector with the Lagrange multiplier, and T 1 2 , , , , 1 is the normalized sampled cost at candidate j.In the proposed Kriging-based 3DMA GNSS, we aim to minimize the cost function.Therefore, we transform the likelihood score to a cost (which is also known as a residual between the measurement and modeled estimation) for the candidate as follows: Subsequently, the costs  δ δ are normalized and used as the sample for Kriging interpolation.
The ordinary Kriging weighting vector ω ω is calculated by considering the distance between the estimated position and all sampled candidates, which represents the spatial covariance of the sampled candidates.The calculation of ω ω can be expressed in matrix form as follows: where λ is the Lagrange multiplier and ⋅ is the Euclidean norm distance between two positions.γ ( ) h is a semivariogram model, which is set as an exponential model in this study: Here, h is the Euclidean norm distance between two positions, and a and b are constants determined by a semivariogram analysis, representing a range parameter and a sill value, respectively (Wackernagel, 2003;Webster & Oliver, 2007).
The semivariogram analysis estimates the average dissimilarity γ * * ( ) h of all sampled candidates: where denotes the set of all pairs of candidates within the set of lag distance h * .For example, h * * { , , , } } 5 10 h k (in meters), and the union x in the sampled candidates.The selected lag distance depends on the separations of candidates and is related to the final number of candidates to support the Kriging interpolation.For a larger lag distance, fewer candidates support the interpolation.In this case, the The range parameter a is the distance h at which the semivariogram first exceeds the sill value.The pseudocode for determining the variogram constant is presented as Algorithm 1. for each element in sampled candidates, as After the constants a and b have been determined for the model, the ordinary Kriging weighting ω ω can be derived: Therefore, the Kriging method used to estimate the cost (residual) at the position of interest can be expressed as in Equation ( 5).The cost function of 3DMA GNSS based on the Kriging method can be expressed as follows: The block matrix Γ and sampled residual  δ δ are obtained by the sampled candidates' positions associated with the candidate-based 3DMA GNSS, which are known terms in this function.
Thus, we can express the error factor of the proposed Kriging-based 3DMA GNSS at epoch t as follows: where and y is a single column array.σ σ t,3 2 DMA is a diagonal variance matrix of the Kriging-based 3DMA GNSS at the x-, y-, and z-axes.The variance at each axis is taken by the weighted distance variation between the Kriging-based 3DMA GNSS solution and all sampled candidates.The weighting is given by the renormalized sampled cost of the corresponding candidate, such that: where c 1 is an empirically determined tuning factor for the 3DMA GNSS error factor.
The modeled error function can be solved as the Kriging-based 3DMA GNSS by applying the Levenberg-Marquardt (LM) algorithm using the Google Ceres Solver (Agarwal & Mierle, 2012).Note that the Kriging-based 3DMA GNSS can only achieve a result identical to that of the candidate-based method, as it only models the sampled candidates mathematically.In other words, if multimodal behavior or drifting occurs, such as in the scenarios shown in Figure 1, a similar degradation should be reflected in the estimated position.Moreover, the error factor can be used in batch optimization, as discussed in the case of FGO using Kriging-based 3DMA GNSS.

DOPPLER MEASUREMENTS
Doppler frequencies and velocity connect two consecutive GNSS position solution epochs.This study integrates received Doppler frequencies as tightly coupled with the Kriging-based 3DMA GNSS.FGO optimizes the error factor formed by each available Doppler measurement separately.The estimated velocity also introduces two error factors connecting two epochs in the FGO process: motion propagation and motion constraint.

Doppler Measurement Modeling
The received Doppler frequency measurement,  D t i , is compared with the modeled value, ˆ, i t D repeated for all received satellites, and forms the error factor for FGO to estimate the receiver's velocity at time t, v t ecef , .The modeled Doppler frequency for the i-th satellite can be expressed as follows: where λ i is the wavelength of the measurement, v t i is the velocity of the satellite provided by the ephemeris, a i is the LOS unit vector from the satellite to the receiver, c is the speed of light, ⋅ δt t i is the satellite's clock drift provided by the ephemeris, and ⋅ δ t t u , is the clock drift of the receiver to the corresponding constellation.
Therefore, the error factor of the Doppler frequency for the i-th satellite at epoch t can be formed: where V W is the variance of the Doppler frequency measurement of the i-th satellite at time t.τ t i ,1 is the uncertainty of the L1-band measurement, as in Equation ( 7) in Ng, Zhang, Luo, & Hsu (2021).
Doppler measurements can also provide information on the clock drift of the receiver.Clock drift should remain stable over a time period.Therefore, the stable clock drift factor can minimize the clock drift of the available satellite constellations, ⋅ δ δt t u , , between two consecutive epochs, which can be expressed as follows:

Motion Propagation
Motion propagation provides a connection in the temporal dimension based on the receiver's estimated positions and velocity.Here, the estimated positions of two epochs, x t and x t +1 , are used, and the error based on the estimated velocity, v t , is calculated.The error factor can be obtained as follows: where ∆t is the time difference between epoch t and t + 1 and is a diagonal covariance matrix associated with the velocity v t ecef , at the x-, y-, and z-axes.c 4 is an empirically determined tuning factor for the motion propagation factor.

Constant-Velocity Motion Constraint
This motion constraint assumes no variation in velocity between two epochs (Li et al., 2018).In other words, the users' motions assume a constant velocity with negligible acceleration.As a result, the optimized trajectory can become smoother.The motion constraint factor minimizes the error between the position change between two epochs and the averaged estimated velocity, which can be expressed as follows: is the average diagonal covariance matrix at time t and t + 1. c 5 is an empirically determined tuning factor for the motion constraint factor.

FGO USING KRIGING-BASED 3DMA GNSS
FGO is performed to integrate the error factors associated with the proposed Kriging-based 3DMA GNSS and Doppler measurements.FGO is used instead of other filtering techniques, such as KF or EKF methods, because of the robust estimation yielded by FGO (Wen et al., 2021), as FGO implements multiple iterations through the estimation process.Furthermore, FGO can accurately reflect the time correlation between measurements and states when measurements do not follow the Gaussian noise assumption, based on a batch of historical data.Therefore, this study uses FGO to optimize the state set as a batch optimization framework, denoted as tightly coupled FGO using Kriging-based 3DMA GNSS (FGO-Krg-Doppler), in the evaluation of this study.A corresponding graph structure is shown in Figure 4.
The state set of the receiver, χ χ, for FGO is expressed as follows: , is a vector on the clock drift of the available satellite constellations.The initial value for the position is provided by the Kriging-based 3DMA GNSS.Meanwhile, the velocities and clock drift are initialized via a least-squares approach using Doppler measurements (Hofmann-Wellenhof et al., 2007).
Therefore, the cost function for the position estimation of the Kriging-based 3DMA GNSS using FGO is formulated as follows: where χ χ * denotes the optimal estimation of the state set.The empirically determined tuning factors in the uncertainty values for each factor, in Equations ( 13) and ( 15) -( 18), are set as c ,

Experimental Setup
Three designed experiments were conducted and are presented in this section: one pedestrian case in New York City, a vehicular case in Hong Kong, and a pedestrian case in Hong Kong.Details of the experiments are presented in Table 1, and trajectories are shown in Figure 5.The average skymask elevation angle, SM, is calculated by averaging the average elevation angle of the skymask at ground truth over the whole experimental period as follows: where SM t az is the skymask elevation angle at the corresponding azimuth angle az of the ground truth at epoch t and T is the number of epochs in the experiment.Here, a larger value of SM represents a more complex environment for the experiment.
The open-source code RTKLIB (Takasu, 2009) is used for GNSS data processing, and the Google Ceres Solver (Agarwal & Mierle, 2012) is used for the NLS and FGO processes.The optimization frame corresponds to the ECEF coordinate system.The ground truth for vehicular cases is given by NovAtel synchronous position, attitude and navigation (SPAN)-CPT, which provides centimeter-level accuracy (Hsu et al., 2021).A reference trajectory for pedestrian walking scenarios is obtained via interpolation by postprocessing.The pedestrian subjects who collected the data were equipped with a smartphone to record the device's fused location output during the experiment.Smartphone-outputted locations are used to calculate the longitudinal speed and project to the vector between the starting and ending locations as reference locations.

Impact of Increasing the Candidate Sampling Separation
One of the limitations of the candidate-based 3DMA GNSS is the unnecessary computational load for candidate sampling.Increasing the separation between FIGURE 4 Structure of the factor graph of the proposed framework candidates while maintaining the sampling radius is one approach for improving the effectiveness of the 3DMA GNSS.The Kriging-based approach expresses discrete candidates as a continuous function.With this approach, the cost value can be interpolated at an unsampled location even if the gap becomes larger if and only if the separation increment does not destroy the gradient within the sampling space.This section uses the vehicular case in Hong Kong (Experiment 2) to evaluate the relationship between candidate sampling separation, positioning accuracy, and computational load.The vehicular case in Hong Kong covers various levels of urbanization complexity, which is better for describing general cases.Figure 6 shows a plot of candidate separations versus the additional positioning error introduced by the Kriging-based 3DMA GNSS as well as the average computational time.Note that the Kriging-based 3DMA GNSS is assumed to achieve a performance identical to that of the candidate-based 3DMA GNSS.Therefore, this evaluation uses the positioning solution of the candidate-based 3DMA GNSS as the initialization and an evaluation baseline with a distribution radius of 40 m and a separation of 2 m.
The results show that the separation distance is proportional to the positioning accuracy of the Kriging-based 3DMA GNSS, as it becomes ambiguous for the continuous function to model and interpolate the overall gradient of the sampled space.Meanwhile, the computational time is inversely proportional to the separation because of the reduction in sampled candidates when the gap increases for the same distribution radius.It can be observed that for the same sampled candidate set (2-m separation), the Kriging-based approach introduces an additional positioning error of approximately 0.1 m compared with the candidate-based approach.One of  the reasons for this result is the solution obtained when a multimodal case occurs.The candidate-based approach attempts to find a weighted average position, and the solution tends to fall in the area between multiple local optima.In contrast, the Kriging-based approach resolves a solution located at either one of the local minimum areas.For the case in which the optimal area is located far from the actual locations, the positioning error of the Kriging-based approach is higher than that of the candidate-based approach with weighted averaging.
As the separation increases, the positioning error continues increasing, whereas the computation time is reduced.When the separation is increased to 4 m, the Kriging-based 3DMA GNSS introduced an additional positioning error of 0.2 m in total, while the computation time for one epoch decreased by nearly one half.Saturation occurs at approximately 8-to 10-m separations, where an additional positioning error of approximately 0.5 m and a duration of 0.4 s are required to calculate the solution.Regarding the balance between positioning accuracy and computational efficiency, a separation of 4 m should be optimal as the setup for the positioning performance evaluation.In the following section, the sampling setup of a 40-m distribution radius with a 4-m separation is used for the Kriging-based 3DMA GNSS.

Evaluation of the Positioning Performance
This section provides postprocessed results and evaluations on the performance of the proposed algorithm based on the setting determined in the previous section.The evaluation is aimed at comparing the proposed algorithms and the following conventional solutions:

Device solution (National Marine Electronics Association [NMEA]):
Position solution outputted by the receiver using the algorithm developed by the manufacturer, for example, the NMEA solution.2. Weighted least-squares (WLS): Conventional WLS method (Takasu, 2009), which uses the pseudorange to estimate the receiver location.The statistics of the postprocessed results for each algorithm are summarized in Table 2.As observed for single-epoch positioning, 3DMA GNSS can provide better positioning than the conventional WLS.Meanwhile, Kriging-based 3DMA GNSS is assumed to achieve a performance identical to that of the candidate-based 3DMA GNSS in Section 3.1, which is confirmed by our experimental results.Positioning accuracy and robustness are improved after the multi-epoch positions are connected.The FGO results outperform EKF results in all experiments.For the integration using FGO, improvement can be found after the 3DMA GNSS error factor is replaced with the Kriging-modeled function compared with the case in which the position is loosely coupled.As the Doppler frequency measurements are further integrated, the position accuracy is improved compared with the results obtained using the Doppler velocity.
The first designed experiment is a pedestrian case in a deep urban area in New York City, for which Figure 7 shows the positioning results and statistics.Both Cand-3DMA and Krg-3DMA achieve an average positioning accuracy of within 20 m.After the Doppler measurements have been integrated, at least a 3-m improvement can be found.The integration of Kriging-based 3DMA GNSS and Doppler measurements can suppress the average positioning error to 10 m.
The candidate-based and Kriging-based 3DMA GNSS approaches should provide a similar positioning performance, whereas the latter approach only expresses the sampled candidates with a mathematical model and solves the model with the LM method.The overall positioning result shows that the assumption is correct in this scenario, and Kriging-based 3DMA GNSS introduces a positioning RMSE that is 0.2 m greater overall than that of the candidate-based 3DMA GNSS.However, two typical limitations of candidate-based 3DMA GNSS, i.e., multiple modes and drifting, cannot be solved by the Kriging-based 3DMA GNSS alone.Two snapshot epochs are selected for further discussion, as shown in Figure 8. Kriging-based 3DMA GNSS shows an error identical to that of the candidate-based 3DMA GNSS.
After the information from Doppler measurements has been integrated, the multimodal issue, shown in Figure 8(a), can be solved by some algorithms.Large position errors can still be found in EKF-cand-vel and FGO-cand-vel because of the domination of the error from the single-epoch solution.In contrast, improvement can be found in FGO-Krg-vel and FGO-Krg-Doppler because of the higher-resolution information from 3DMA GNSS and Doppler measurements that are input to the optimization.This result shows that the Krg-3DMA GNSS improves the integrated solution.In comparison, Cand-3DMA determines the solution by averaging candidates and their scores, and Krg-3DMA follows an optimization approach to converge to a minimal-cost location.Krg-3DMA is more sensitive to multimodal issues, with a significant bias.Therefore, when the velocity constraint is added, this issue can be significantly improved to reveal the benefits of the Krg-3DMA method in determining the correct minimal-cost area.However, integrating the temporal domain information cannot solve the limitation of the drifting case in Figure 8(b) because the sampled candidates consistently provide wrong information.Therefore, these candidates with a small sampled cost are concentrated and cause a small estimated variance based on Equation ( 13), which cannot bound the actual positioning error.The wrong information deteriorates the optimization, resulting in the incorrect result that the 3DMA GNSS is accurate and causing deviation in the determined solution.Interestingly, although FGO-Krg-vel and FGO-Krg-Doppler are both affected by drifting at the epoch chosen here, these methods prove to be beneficial by providing higher-resolution information to the optimization.The drifting is obviously slower than that for FGO-cand-vel, especially in the earlier period.
Another important indicator that influences the optimization performance is the estimated variance of 3DMA GNSS.The estimated variance should bound the actual positioning error to provide optimal performance.However, especially for the drifting case, the estimated variance is small whereas the positioning error is large, as the candidates with low cost are concentrated and consistently provide wrong information.Figure 9 shows the actual positioning error and estimated standard deviation (STD) of single-epoch 3DMA GNSS and the FGO-integrated algorithms.
In some of the periods, the actual positioning error is large while the estimated variance is small; such cases indicate drifting, as shown in the two yellow frames in Figure 9.The optimization process of FGO-Krg-vel wrongly trusts the 3DMA GNSS and causes the result to be driven by the error of 3DMA GNSS as well.In contrast, tightly integrating Doppler measurements can help in some cases.FGO-Krg-Doppler provides stronger constraints on the temporal domain by directly using the Doppler measurements, and a better positioning performance can be achieved, as shown in the two yellow framed periods in Figure 9.
However, tight coupling of the Doppler measurements fails to constrain the position in some cases.For example, in the period in the purple frame in Figure 9, the estimated STD is similar for the Kriging-based 3DMA GNSS and FGO-integrated solution.Two FGO results successfully suppress the positioning error but show similar positioning error trends.This result demonstrates the importance of variance estimation for integration and optimization, especially for the identification of drifting cases.Further work is needed to detect outliers in the 3DMA GNSS.
The second experiment is a vehicular case with various urban environments in Hong Kong, for which the positioning results and statistics are shown in Figure 10.This experiment begins in a relatively open-sky environment, continues into a deep urban area, and then returns to the relatively open-sky region as a closed loop.The positioning error is reduced by nearly one half with the aid of the 3D building model, from an average of 11.75 m for WLS to approximately 6.5 m.For the temporally connected algorithm, both the positioning accuracy and robustness are improved by at least 1 m.For the loosely coupled approaches, FGO outperforms EKF by approximately 1 m in terms of accuracy and deviation.A slight positioning improvement of 0.6 m can be found after tightly integrating the Kriging-based 3DMA GNSS and Doppler measurements.However, the receiver output solution still outperforms the proposed algorithm by approximately 1 m on average in this case.
Two cases are selected for a detailed comparison between each 3DMA GNSS method, i.e., the EKF and FGO performance, for a stationary case in which the vehicle is waiting for a traffic light (blue framed period in Figure 10) and a moving case in which the vehicle makes a right turn (green framed period in Figure 10), as shown in Figures 11 and 12, respectively.In the first case study, in which the vehicle is waiting for a traffic light, the vehicle remained stationary while waiting for the traffic light during epochs 11 to 64.Thus, any changes in position or velocity during these instances are caused by unmodeled noise.Between the epochs of approximately 20 s to 40 s, EKF-cand-vel, FGO-cand-vel, and FGO-Krg-vel are affected by the position error of the single-epoch 3DMA GNSS.For FGO-Krg-Doppler, the Doppler measurements are tightly coupled in optimization; thus, this method can perform better in de-weighing wrong or outlier measurements to provide stronger constraints in the zero-velocity case.
The second case is a moving case in which the vehicle is turning right.EKF introduced a large error caused by drifting of the 3DMA GNSS result.Among the three FGO approaches, the FGO-cand-vel provided an accuracy similar to that of FGO-Krg-Doppler before the turn is made but showed the worst performance after epoch 20.This result shows that only loosely coupling the position and velocity cannot provide a constraint that is sufficiently strong to achieve satisfactory accuracy.After the position error factor is replaced with a Kriging-modeled factor, the position error can be suppressed as the vehicle drives straight after making the turn in epoch 25; however, a relatively large positioning error can be found when the vehicle moves out of the turn between epoch 20 and 25.Integrating the Kriging-based 3DMA GNSS and Doppler measurements achieves a smaller overall positioning error, and the whole trajectory shows a substantial amount of overlapping with the reference trajectory.Some relatively large errors are introduced during the turn (epoch 18) and after the turn (epoch 28), which primarily contribute to the longitudinal side of the driving direction.The final designed experiment is a pedestrian case in an urban environment in Hong Kong, as shown in Figure 13.The walking trajectory is a straight line that runs parallel to a building.The average positioning error of the single-epoch 3DMA GNSS is approximately 8.5 m.After the Doppler measurements are integrated, the positioning error is reduced to 3.4 m on average.
The positioning statistics indicate that the two FGO approaches using the Kriging-based 3DMA GNSS (FGO-Krg-vel and FGO-Krg-Doppler) can provide better positioning accuracies, especially between epochs 30 and 70 and after epoch 120.However, there was a large positioning error around epoch 110 for FGO-Krg-vel.This poor performance is dominated by drifting of the lowest sampled cost area at the longitudinal side of the direction of motion.

CONCLUSIONS AND FUTURE WORK
This study makes two main contributions.First, we modeled the candidate-based 3DMA GNSS via an ordinary Kriging method, allowing the receiver location to be iteratively estimated using the NLS approach.The Kriging system was established, and the LM method was used to determine the receiver location; this method is denoted the Kriging-based 3DMA GNSS method.Second, we integrated the Kriging-based 3DMA GNSS and Doppler measurements through FGO.In this manner, multi-epoch states could be optimized through FGO to obtain a more robust position solution.According to the postprocessing results of actual data recorded in Hong Kong and New York City, the performance of the Kriging-based 3DMA GNSS is comparable to that of candidate-based 3DMA GNSS in most cases.Moreover, the FGO method using Kriging-based 3DMA GNSS can optimize the states for multiple epochs, and a smoother trajectory can be obtained.The proposed FGO using Kriging-based 3DMA GNSS was compared with FGO and EKF using candidate-based 3DMA GNSS with loose coupling.The tightly coupled FGO using Kriging-based 3DMA GNSS achieves superior positioning results with better positioning accuracy and robustness, giving an average positioning RMSE of approximately 6 m over all three designed experiments.A case study showed that the proposed method can suppress the multimodal issue of 3DMA GNSS to improve positioning performance.Moreover, the proposed method uses fewer position candidates and thus corresponds to a decreased computational load, and the optimal setup of the separation is set as 4 m.To improve the candidate sampling efficiency, the uncertainty of the previous epoch's solution and velocity can be considered to determine the sampling area in the future.
However, the results show that drifting of the candidates is still a limitation of 3DMA GNSS, where the adaptive variance estimation in Equation ( 13) cannot detect the occurrence of drifting.Identification of these outlier epochs remains to be explored.For the multimodal limitation in the position domain, inputting the samples as multiple hypotheses into FGO may better represent the occurrence of a Gaussian mixture.Additionally, the proposed Kriging-based 3DMA GNSS still relies on the fundamental basis of the candidate-based 3DMA GNSS; thus, future work will focus on studying the sampling area to achieve optimal performance for the Kriging method.Moreover, mathematically expressing the 3DMA GNSS algorithm or even 3D building models as a differentiable function input to FGO can maximize the contribution of 3DMA GNSS.
To enhance the positioning performance in urban canyons, it is necessary to introduce carrier-phase measurements in the algorithm.The introduction of carrier-phase measurements can aid in positioning, for instance, by integrating time-difference carrier-phase (TDCP) measurements with RTK GNSS.TDCP assumes that no cycle slip occurs between consecutive epochs, and thus, no ambiguity resolution is required.In other words, when TDCP measurements are used, no reference station is required for positioning, and position estimation is not relative to a reference station.We previously demonstrated that the use of a 3D building model can enhance the GNSS RTK in urban canyons (Ng & Hsu, 2021).In a future study, the benefits of introducing carrier-phase measurements in the proposed Kriging-based 3DMA GNSS using FGO can be evaluated.

a c k n o w l e d g e m e n t s
This research was supported by the University Grants Committee of Hong Kong under the Research Impact Fund for project R5009-21: "Reliable Multiagent Collaborative Global Navigation Satellite System Positioning for Intelligent Transportation Systems." The authors thank John-Ross Rizzo, MD, and his team from the Grossman School of Medicine, New York University for conducting experiments in New York City and for sharing the raw GNSS data for this work.

FIGURE 1
FIGURE 1 Example of typical limitations of candidate-based 3DMA GNSS for (a) multiple modes/local minima and (b) drifting

FIGURE 2
FIGURE 2 System architecture of the proposed multi-epoch Kriging-based 3DMA GNSS using FGO

FIGURE 3
FIGURE 3 Process flow for calculating the cost function in the proposed Kriging-based 3DMA GNSS (a) Likelihood score of the candidate-based 3DMA GNSS; (b) sampled cost of the candidates calculated from the likelihood score (Equation (6)); and (c) interpolated cost function based on the ordinary Kriging method.Eq.: equation.

FIGURE 6
FIGURE 6 Average additional positioning error introduced by Kriging-based 3DMA GNSS and average computation time for one epoch solution Note that the comparison baseline of the positioning error is the comparison between the candidate-based 3DMA GNSS and the proposed Kriging-based 3DMA GNSS.

FIGURE 9
FIGURE 9 Positioning statistic of Experiment 1 on positioning error and the STD estimated from Equation (13)

TABLE 1
Environment of the Designed Experiments Note: Exp: experiment; GPS: global positioning system; GLONASS: globalnaya navigazionnaya sputnikovaya sistema or global navigation satellite system

TABLE 2
Positioning Result Statistics