{'fraction_non_alphanumeric': 0.030517305545217716, 'fraction_numerical': 0.007294380349832527, 'mean_word_length': 4.987299465240642, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Evaluation of Kinematic Precise Point Positioning Convergence with an Incremental Graph Optimizer

Estimation techniques to precisely localize a kinematic platform with GNSS observables can be broadly partitioned into two categories: differential, or undifferenced. The differential techniques (e.g., real-time kinematic (RTK)) have several attractive properties, such as correlated error mitigation and fast convergence; however, to support a differential processing scheme, an infrastructure of reference stations within a proximity of the platform must be in place to construct observation corrections. This infrastructure requirement makes differential processing techniques infeasible in many locations. To mitigate the need for additional receivers within proximity of the platform, the precise point positioning (PPP) method utilizes accurate orbit and clock models to localize the platform. The autonomy of PPP from local reference stations make it an attractive processing scheme for several applications; however, a current disadvantage of PPP is the slow positioning convergence when compared to differential techniques. In this paper, we evaluate the convergence properties of PPP with an incremental graph optimization scheme (Incremental Smoothing and Mapping (iSAM2)), which allows for real-time filtering and smoothing. The characterization is first conducted through a Monte Carlo analysis within a simulation environment, which allows for the variations of parameters, such as atmospheric conditions, satellite geometry, and intensity of multipath. Then, an example collected data set is utilized to validate the trends presented in the simulation study.


I. INTRODUCTION
The ability to precisely localize a platform is of paramount importance to a myriad of fields (e.g., augmented reality [1], autonomous navigation [2], and natural hazard monitoring [3]).
To facilitate the precise localization of the platform, several navigation aids can be utilized (e.g., vision [4], lidar [5], inertial [6]). One navigation aid that is commonly utilized for terrestrial applications is a global navigation satellite system (GNSS) receiver. The signals propagated by a GNSS satellite provides the algorithm with information that allows for accurate, global localization of the platform.
One commonly used methodology for processing GNSS signals is the Precise Point Positioning (PPP) approach [7]. The PPP algorithm utilizes the dual-frequency undifferenced GNSS observables, which allows the technique to operate without the need of external reference stations. The undifferenced observations are use along with precise GNSS orbit and clock bias products to mitigate the errors removed through observation differencing [8]. The orbit and clock products that enable the PPP method to achieve decimeter level positioning can be broadcast to an end-user in real-time (e.g., L-band [9] and Iridium modem link [10]).
Real-time kinematic-PPP (kPPP) provides similar positioning performance when compared to traditional differential GPS (DGPS) (i.e., Real Time Kinematic) for dynamic platforms [11]. The comparable positioning accuracy without the need for a nearby static GPS reference station makes it an attractive processing formulation. However, it has been noted in several studies that the PPP formulation has a longer convergence period than comparable differential techniques [12], [13].
In an attempt to decrease the initial convergence period, there has been a plethora of research into augmenting the PPP approach with additional information sources. One of the most commonly utilized augmentation sources for traditional single constellation PPP is additional GNSS observables. One example of this type of augmentation is the incorporation of multiple constellation observations [14], [15]. Another example of this type of augmentation is the PPP-RTK formulation [16] which provides faster convergence by enabling integer ambiguity resolution [17]. Another well studied form of PPP augmentation is the tightly coupled PPP inertial navigation (INS) formulation [18], which has also been shown to decrease the initial convergence period of PPP [19]. However, all PPP augmented methods require additional infrastructure (e.g., a network of reference stations, or additional sensors on-board), which can be prohibitive for many applications.
Another method to decrease the convergence period of PPP is to utilize a novel optimization framework. This has the potential to provide a benefit over the previously discussed PPP methods because all the previously provided methods utilizing the same underlying optimization framework (i.e., a variant of the Kalman filter [20]- [22]). Where this framework estimates the desired states by marginalizing all prior information and propagating with dynamic models to the next time step. For PPP, where a subset of the desired states are not observable over a single epoch (i.e., the carrier-phase ambiguity states), this may not be the best framework.
In this paper we evaluate the convergence properties of PPP utilizing an incremental graph optimization framework that allows for real-time smoothing. This work relies upon advances made within the robotics community on efficient, real-time smoothing. Where research into smoothing, within the robotics community, has been dominated by graph based methodologies since the seminal paper on the subject was published in 1997 [23]. When [23] was published, graphbased smoothing was not widely utilized due to computation complexity of solving the initial formulation. However, quickly thereafter, methods were proposed to greatly reduce complexity through the utilization of factor graphs [24]. The √ SAM formulation as presented in [25] was particularly influential as it provided connections between the factor graph formulation and sparse linear algebra. The idea of batch factor graph optimization was later extended to an incremental inference framework in [26], [27]. The work presented in [27] provides a frame-work to conduct real-time [28], non-linear graph based filter and smoothing.
The rest of this paper is organized in the following manner. First, the technical approach will be discussed. The technical discussion will provide an overview of factor graph optimization for GNSS optimization, and the ability to incrementally updating using the Bayes tree data structure. Next the discussion will shift to the evaluation of the algorithm with both simulated and collected datasets. Finally, some concluding remarks and future work will be discussed.

A. Factor Graphs
The ability to conduct accurate and efficient inference is at the center of all navigation algorithms. One way to represent the inference problem is with a probabilistic graphical model [29], which can take several forms. One convenient graphical model for conducting state estimation is the factor graph [24]. At a fundamental level, the factor graph provides a convenient framework for factorizing a function of several variables into smaller subsets. More explicitly, this model provides a useful framework for factorizing the posterior distribution, which allows for efficient calculation of the state vector that maximizes the a posteriori distribution. The factorization is represents as a bipartite graph, G = (F , X , E), where there are two types of vertices: the states to be estimated, X , and the probabilistic constraints applied to those states, F . An edge E only exists between a state vertex and a factor vertex if that factor is a constraint applied on that time-step. An example factor graph is depicted in Fig.1, where X n represents the states to be estimate at time-step n, ψ p,n−1 represents prior information about the estimated states at time-step n − 1, ψ b,n represents the motion model of the platform from time-step n− 1 to n, and ψ l represents the constraint applied to the state by a measurement (e.g., a GNSS pseudorange observable). As previously mentioned, the factor graph provides a factorization over the posterior distribution, p(X|Z). Thus, we can easily calculate the state vector that maximizes the posterior (MAP) by finding the state vector that maximizes the product of factors, as depicted in Eq. 1.X = argmax For a through discussion on factor graph based state estimation the reader is refered to [30]. The optimization problem presented in Eq. 1 can be reduced to non-linear least-squares formulation if Gaussian noise is assumed, as provided in Eq. 2.
Now that a general discussion of the factor graph framework has been provided, we can proceed by constructing GNSS specific factors. For this work, we will detail the construction of two factors: the GNSS observation factor, and the carrierphase bias factor.

B. Constructing the GNSS Observation Factor
To allow autonomy of the PPP approach from local reference stations, the undifferenced dual-frequency GNSS observables are utilized. Due to the undifferenced nature of the observations, the PPP processing technique must incorporate GNSS error mitigation models -these models provide corrects for the corrupting sources that would be mitigated through observation differencing -to provide an accurate positioning solution. The sources that corrupt a GNSS observation can be segregated into three partitions: the error contributed by the propagation medium, the error contributed by the control segment, and the error contributed by the user.
To begin constructing our measurement model, the method implemented to mitigate the propagation medium errors are discussed. The error attributed to the propagation medium is composed of delay due to the ionosphere and the delay due to the troposphere. To mitigate the ionospheric delay, we leverage the dispersive nature of the medium, and a linear combination of the GPS L 1 and L 2 frequencies (1575.42 MHz and 1227.60 MHz, respectively) is formed to produce ionospheric-free (IF) pseudorange and carrier phase measurements [8]. The IF combination of an observable, O j , can be seen in Eq. 3.
To mitigate the error attributed to the troposphere, both the wet and the dry component of the troposphere must be modeled, as shown in Eq. 4. For this study, the Hopfield model [31] is used to model the dry component of the troposphere. To compensate for the wet delay -the wet component only accounts for approximately 10% of the total troposphere error -and the residual error of the dry delay model, a stochastic random variable is added to the state vector.
To mitigate the error attributed to the control segment, the PPP approach utilizes orbit and clock corrections. These global corrections are generated through a network of reference stations.
Finally, a discussion on the user error segment is provided. The user error segment is composed of two sources: multipath error, and receiver thermal noise error. For this study, no methods were implemented to explicitly model the user segment error; however, as noted in [32], the magnitude or the user error is proportional to the elevation angle between the platform and the satellite so, within this evaluation, the uncertainty in the observation is scaled by the elevation angle. It should be noted that PPP observational models for moving platforms typically include corrections for relativistic effects (i.e. from the GPS broadcast correction), receiver and satellite antenna phase center variation, and carrier-phase wind-up; however, these effects were neglected within this simulation study. Additionally, dynamic platform generally couple inertial information with the GNSS observables to mitigate uncertainty in the platforms dynamic model [19].
Utilizing the provided error mitigation techniques, the PPP observation model can be constructed. The pseudorange and carrier-phase measurements are modeled as shown in Eq. 5 and Eq. 6, respectively: where, R j = ||x s − x u || is the geometric range between the platform and the j th satellite, δt u is the receiver's clock bias, δt s is the satellite's clock bias, T z is the tropospheric delay in the zenith direction, M d (el j ) is a user to satellite elevation angle dependent mapping function, δ Rel. is the correction attributed to relativistic effect [33], δ P.C. is the correction attributed to the offset between the satellite's center of mass and the phase center of the antenna [34], δ D.C.B is the differential code bias correction [31], δ W.U. is the correction attributed to the windup effect on the phase observables [35], λ IF is the wavelength corresponding to the IF combination, and N IF is phase ambiguity. In Eqs. 5 and 6 the remaining unmodelled error sources are indicated with ǫ. To implement the provided observation model in software, the open-source library GPSTk [36] is utilized.
Using the PPP observation model, we can construct a GNSS constraint for the factor graph [37]. To begin, we note that the GNSS observations are providing a set of likelihood constraint, L(O|X), on the optimization process. If the assumption is made that the state and measurement noise models are Gaussian, then this constraint can be incorporated into factor graph through the mahalanobis distance, as provided in Eq. 7, where z is the observed measurement,ẑ is the estimated measurement -calculated using Eq's 5 and 6 -and Σ is the uncertainty in the observation.
C. Incorporating the Carrier-Phase Ambiguity States There are several way in which the carrier-phase ambiguity states can be incorporated into the factor graph. One such way is to incorporate a new carrier-phase ambiguity state for each epoch. Consecutive carrier-phase ambiguity states can be constrained by a process noise update. The measurement Jacobian associated with this graph construction is represented in Fig. 2.A. From Fig 2.A, we can see a measurement Jacobian that is more densely populated than desired.
To construct less densly populated measurement Jacobian (i.e., a more efficient optimization scheme), we can leverage the knowledge that the true carrier-phase ambiguity value for a given satellite within a continually tracked phase-arc is a constant value. Due to this property of the true ambiguity value, the carrier-phase ambiguity factor can be represented as a random constant variable. Where, initially, a single factor is added for each satellite, and a new factor is added only if the there is a cycle-slip or if a new satellite is tracked. By treating the carrier-phase bias factor in this manner (i.e., like a "landmark" variable in traditional pose-graph SLAM [38]), and utilizing the Bayes tree based optimizer, an efficient real-time smoothing formulation for GNSS signal processing is presented. The measurement Jacobian associated with this graph construction is represented in Fig. 2.B. From Fig. 2.B, we see a less densely populated measurement Jacobian, as desired.

D. Incremental Factor Graph Inference
The formulation presented in the previous sections provides an efficient estimation of X MAP when all of the information is provided a priori. However, it is generally the case that information is arriving sequentially, and it is desired to incrementally provide state estimates. The ability to provide an incremental estimator lies in the capability of the optimizer to reuse prior computations. A well studied technique of computation reuse, for state estimation, is to employ QR-factorization to update the previous matrix factorization [22], [26]; however, this technique only works for linearized systems.
To overcome this limitation, the Incremental Smoothing and Mapping (iSAM2) formulation was developed [27]. The iSAM2 formulation allows for incremental inference over linear or non-linear objective functions through the utilization of a novel graphical model, the Bayes tree [39]. To provide insight into this formulation, specifically for GNSS applications, a simple GNSS example will be presented. Where it will be shown how to convert the GNSS factor graph into a Bayes tree. Additionally, a discussion will be provided on how the Bayes tree graphical models allows for efficient inference.
To begin our discussion, a factor graph that represents the GNSS inference problem is presented in Fig. 3. With this factor graph, it is desired to estimate the states {X, B}. In this formulation, X represents the position, troposphere, and receiver clock bias states, as provided in Eq. 8. Additionally, the vertices B represents the carrier-phase bias states.

Prior Factors Between Factors Pseudorange Factors
Carrier Phase Factors Utilizing the factor graph presented in Fig. 3, we can begin the process of converting a factor graph into a Bayes tree. To do this, we must take an intermediate step and construct a Bayes net. The Bayes net can be constructed from the factor graph using a variable elimination game [40]. For our specific example, the Bayes net is provided in Fig. 4, if the bias states are eliminated first then the positing states (i.e., the variable elimination is B 1 , B 2 , X 1 , X 2 , X 3 ). It should be noted that if the elimination ordering is varied, the resultant Bayes net will change and this can have a substantial impact on the run-time of the optimizer [41].
Utilizing the previously constructed Bayes net, we can now construct the Bayes tree. The Bayes tree is constructed to take advantage of the clique structure within the Bayes net. That is, by re-writing the Bayes net we are left with a directed tree structure, where the vertices in the graph represent cliques in the original Bayes net. For our GNSS example, the constructed Bayes tree is provided in Fig. 5.
The tree structure present in the Bayes tree plays a pivotal role in the ability of the data structure to provide an efficient incremental inference engine. The primary advantage of the tree structure is in the idea that only local sections (i.e., a branch in the Bayes tree structure) of the data structure needs to be re-linearized when new constraints are added to the graph. Where re-linearization can be conducted by converting a subset of the Bayes tree back into a factor graph adding the new constraint [30].

III. EXPERIMENTAL SETUP
To conduct an analysis of PPP convergence, a simulation environment was constructed. For synthetic observations generation, the SatNav-3.04 Toolbox [42] is utilized, which provides a Matlab environment for generating dual-frequency pseudorange and carrier-phase observations for a specified trajectory. For this evaluation, four trajectories of varying dynamic were created -an example flight trajectory is provided in Fig. 6. To tailor the toolbox for an evaluation of kinematic PPP airborne positioning, several minor modifications were made, as discussed in [43]. For example, toolbox was modified to include attitude dependent satellite masking and carrierphase breaks. That is, when a satellite is obscured or nearly obscured due to a change in platform attitude, it is masked from view and the the potential of a carrier-phase breaks is increased. Additionally, for a PPP analysis, a methodology of constructing an orbit and clock model is required. This error model is constructed by differing JPL's International GNSS Service (IGS) submission with European Center for Orbit Determination (CODE) submission. For a more detailed discussion on the simulation environment, the reader is directed to [43]. To evaluate the positioning performance of the PPP incremental graph optimizer, a Monte Carlo style experiment was implemented. Specifically, one hundred datasets were generated where several parameters, which are known to adversely affect GNSS positioning performance, were randomly initialized for each flight -see Table I for additional information on varied parameters. To provide a reference positioning solution, a traditional extended Kalman filter (EKF) was utilized, where the specific EKF formulation details are provided in [43]. To provide a fair comparison, the same stochastic models were implemented for both estimators. The specific stochastic models utilized for the comparison are provided in Table II.

A. Example Simulated Flight Evaluation
To begin an evaluation of the PPP incremental graph optimizer, a single data set, which is representative of all datasets simulated for this study, will be analyzed. As a starting point, the residual sum of squares (RSOS) positioning error is utilized to evaluated the performance of both estimators, as shown in Fig. 7. From Fig. 7, it can be readily seen that the incremental graph optimizer more quickly converges -when compared to the EKF -to a steady-state value. Additionally, it should be noted that the both estimators converge to approximately the same value. The RSOS positioning error statistics for both estimator are provided in Table III. From Table III it should be noted that the incremental graph optimizer outperforms the EKF with respect to all metrics provided (e.g., the incremental graph optimizer provides a 25 cm error reduction with respect to the mean RSOS positioning error). To continue an analysis of this example data set, it can be seen in Fig. 7 that the most substantial positioning error reduction attributed to the incremental graph optimizer occurs during the first several minutes of the flight (i.e., during the PPP convergence period). To provide insight into the accelerated convergence rate of the incremental optimizer, next, an evaluation of both estimators ability to correctly resolve the phase bias states is provided in Fig. 8. Where it can be seen that the incremental graph optimizer provides a substantial carrier-phase bias estimation error when compared to the EKF.

B. Positioning Performance Over All Simulated Flights
Now, the evaluation shifts from a single flight to the performance of both estimators over all simulated data sets. As with the previous evaluation of a single flight, we will utilize the RSOS positioning error as the metric of comparison.
To begin our evaluation, the cumulative distribution function (CDF) of the RSOS positioning error for both estimators is evaluated, as provided in Fig. 9. From Fig. 9 it can be noted that there is a considerable shift to the left for the CDF of the RSOS positioning error of the incremental graph optimizer for  To confirm that RSOS positioning error seen in Fig. 9 for the incremental graph optimizer is occurring during the initial convergence period, a CDF of the RSOS positioning error for both optimizers during the first 15 minutes of each data set is provided in Fig. 10. As indicated by the right shift in Fig. 10 of the EKF RSOS positioning error line, the incremental graph optimizer provides a more accurate positioning solution during the initial convergence period. The specific RSOS positioning error reduction during the initial convergence period can be seen in Table V.

C. Example Evaluation With Collected Data
Finally, to verify the positioning performance benefits noted in the simulation study, a similar analysis is conducted on  an example collected dataset. The dataset to be evaluated was collected on-board a small, fixed-wing Unmanned Aerial Vehicle (UAV). This UAV testbed (Phastball) -as depicted in Fig. 11 -was developed at West Virginia University as a research platform [44]. The Phastball is equipped with a NovAtel OEM-615 dualfrequency GNSS receiver, which provides 10 Hz GNSS observables over the duration of the flight. The flight profile is depicted in Fig. 12. A second OEM-615 NovAtel GNSS receiver was placed near the runway to allow for a postprocessed RTK solution, where RTKLIB [45] was utilized to generated the reference solution. Fig. 11: Phastball research platform [44] in flight over the West Virginia University Jackson's Mill airfield.
Utilizing this dataset, the PPP incremental graph optimizer is evaluated against a Kalman filter based PPP approach, where both estimators are given the same inital conditions (i.e., both estimators are provided the same initial error covariance, the same measurement noise model, and the same process noise model). In Fig. 13 the 3D RSOS positioning error for both estimators is provided. The result presented in Fig. 13 follow the trend provided by the simulation study (i.e., the incremental graph optimizer provides faster positioning error convergence than that provided by the Kalman filter). This is further validated by looking at Table VI, where a substantial median 3D RSOS positioning error reduction is granted by the PPP incremental graph optimizer when compared to the Kalman filter.  The desire to precisely localize a platform is of paramount importance to a myriad of fields. This desire has lead to a plethora of research into precise GNSS localization due to its ability to provide a precise and globally consistent solution. One of the most commonly utilized GNSS formulations is the precise point positioning (PPP) technique due to its autonomy from local reference stations. However, it has been noted in several studies that PPP has a relatively long initial convergence period when compared to differential techniques.
To reduce the convergence time of PPP, this paper proposes the use of recent advances in real-time smoothing made within the robotics community. Specifically, this paper makes connections between GNSS localization and incremental posegraph optimization. The connection between the two fields lies in the ability to treat phase bias states as "landmark" nodes in the graph. By treating the phase bias state in this manner, and utilizing a Bayes tree based optimizer, efficient smoothing of the position states can be conducted in real-time.
To the quantify the benefit of this formulation, a Monte-Carlo style experiment was conducted within a simulation environment. Utilizing the simulated data, the incremental graph optimization was evaluated along with a traditional EKF-PPP formulation. Through this evaluation, it was found that the incremental graph optimization technique provided a substantial RSOS positioning error reduction during the initial PPP convergence period when compared to a traditional EKF formulation. This finding was also validated with an evaluation of a short duration dataset collected with a fixed-wing UAV.
Finally, to allow for external validation and collaboration, all software developed for this evaluation has been released publicly at github.com/wvu-navLab. Included with the source code are several example datasets.

VI. FUTURE WORK
In this evaluation, the only comparison solution was generated by a traditional EKF-PPP formulation where phase biases are estimated as floating parameters. However, there are several additional formulation that are known to provide faster convergence rates (e.g., integer ambiguity enabled PPP). With that in mind, there is a need to evaluate the incremental graph optimizer against other state of the art formulations, and to leverage these techniques within the graph.