A Precise Point Positioning Method Based on Factor Graph Optimization

Abstract

The present invention belongs to the technical field of precise point positioning, and particularly relates to a precise point positioning method based on factor graph optimization. First, a PPP observation model is established. Among them, the state variables in the PPP observation model include three-dimensional position, receiver clock error, tropospheric zenith wet delay component, and ambiguity in the carrier phase. Then, the state variables are used as variable nodes of the PPP factor graph optimization model, the GNSS pseudorange observation value and the carrier phase observation value are used as observation constraints to form factor nodes of the PPP factor graph optimization model, and the GNSS pseudorange observation error and the carrier phase observation error are used as edges connecting some variable nodes and the corresponding factor nodes, so as to construct a PPP factor graph optimization model. Finally, optimization calculation is carried out based on the constructed PPP factor graph optimization model to obtain a PPP positioning result with the minimum error. Compared with the traditional Kalman filter, the PPP positioning method based on factor graph optimization can achieve comparable accuracy and eliminate the convergence process in the Kalman filter.

Description

Technical Field

[0001] The present invention belongs to the technical field of precise point positioning, and particularly relates to a precise point positioning method based on factor graph optimization. Background Art

[0002] In the high-precision positioning of the Global Navigation Satellite System (GNSS), Precise Point Positioning (PPP) has attracted extensive attention since its birth. PPP only requires a single GNSS receiver to directly obtain high-precision absolute coordinates globally, and is widely used due to its economic efficiency.

[0003] In 1997, researchers at the Jet Propulsion Laboratory (JPL) of the United States first proposed a single-point positioning method that uses pre-determined high-precision satellite orbits and clock error products, simultaneously performs ionospheric-free combinations on GPS dual-frequency carrier phase and pseudorange observations as combined observations, and applies square root information filtering for parameter estimation. This is the initial precise point positioning model. Experimental results show that it can achieve static positioning accuracy at the centimeter level and dynamic positioning accuracy at the sub-decimeter level. Subsequently, different scholars at home and abroad have successively conducted fruitful research and improvement on precise point positioning technology. Kouba and Héroux further introduced the error correction model in precise point positioning technology and applied the sequential filtering parameter estimation method, with static positioning accuracy reaching the centimeter level. Professors Liu Jingnan and Ye Shirong of Wuhan University conducted in-depth discussions on precise point positioning technology, explained some of the key technologies, and achieved centimeter-level static positioning and decimeter-level dynamic positioning using Kalman filtering; Professor Zhang Xiaohong took the lead in developing the positioning software Trip based on precise point positioning technology in China. The parameter estimation method also uses Kalman filtering, achieving accuracy comparable to that of foreign similar software and being successfully applied in the field of aerial surveying. Currently, there are mainly four methods to accelerate the convergence speed of PPP: fixing non-differenced ambiguities; multi-frequency observations, multi-system observations, and PPP-RTK.

[0004] Summarizing the above research results, it can be seen that all the parameter estimation algorithms in traditional PPP research are Kalman filtering and its derivative forms. Although filtering algorithms have advantages such as high efficiency and low memory occupancy, the filtering-based methods have disadvantages such as being greatly affected by linearization errors, insufficient utilization of observation data, difficulty in supporting plug-and-play, long convergence time, complex processing of delayed observation data, and low positioning accuracy. Summary of the Invention

[0005] The object of the present invention is to provide a precise point positioning method optimized based on a factor graph to solve the problems of long convergence time and low positioning accuracy caused by precise point positioning based on the Kalman filtering algorithm and its derivative forms.

[0006] To solve the above technical problems, the present invention provides a precise point positioning method optimized based on a factor graph, including the following steps:

[0007] 1) Establish a PPP observation model, where the state variables in the PPP observation model include three-dimensional position, receiver clock error, tropospheric zenith wet delay component, and ambiguity in the carrier phase;

[0008] 2) Use the state variables as variable nodes of the PPP factor graph optimization model, use GNSS pseudorange observation values and carrier phase observation values as observation constraints to form factor nodes of the PPP factor graph optimization model, and use GNSS pseudorange observation errors and carrier phase observation errors as edges connecting some variable nodes and corresponding factor nodes, thereby constructing a PPP factor graph optimization model:

[0009] 3) Perform optimization calculations based on the constructed PPP factor graph optimization model to obtain a PPP positioning result with the minimum error.

[0010] The beneficial effects are as follows: The present invention applies factor graph optimization to PPP solution. Specifically, the state variables are used as variable nodes of the PPP factor graph optimization model, GNSS pseudorange observation values and carrier phase observation values are used as observation constraints to form factor nodes of the PPP factor graph optimization model, and GNSS pseudorange observation errors and carrier phase observation errors are used as edges connecting some variable nodes and corresponding factor nodes, thereby constructing a PPP factor graph optimization model, and then performing solution to obtain a better PPP positioning result. Compared with the traditional Kalman filtering, the PPP positioning method based on factor graph optimization can achieve comparable accuracy. Moreover, the batch optimization data processing mode of factor graph optimization can better separate vertical direction errors, improve the strong correlation between vertical direction errors and receiver clock errors and tropospheric zenith wet delays, so that the accuracy in the vertical direction is better in static positioning. In addition, the PPP positioning method based on factor graph optimization can give the optimal estimation of all historical epochs in real time, eliminating the convergence process in Kalman filtering.

[0011] Further, the state variables within the sliding window are expressed as:

[0012] χ = [x0, x1,..., x n , N0, N D ..., N m

[0013] x k ​= [x, y, z, dt r , dT], k ∈ [0, n]

[0014] where χ represents the state variable; N i represents the ambiguity at the i-th frequency, i = 1, 2, …, m, where m represents the number of ambiguities; x, y, z represent the three-dimensional position; dt r , dT represent the clock error and the tropospheric zenith wet delay component respectively; n represents the size of the sliding window.

[0015] Furthermore, the factor node corresponding to the GNSS pseudorange observation value is called the pseudorange factor, and the factor node corresponding to the carrier phase observation value is called the carrier phase factor; the pseudorange factor is a unilateral factor and is only linked to the corresponding factor node x k , and the carrier phase factor is a bilateral factor, with both ends linked to the corresponding factor node x k and the ambiguity N k respectively.

[0016] Furthermore, when performing the optimization calculation in step 3), it is necessary to solve the maximum a posteriori estimation of the state vector according to the observation information and the prior probability, where all the observation values are independent of each other and follow a Gaussian distribution with zero mean.

[0017] Furthermore, when performing the optimization calculation in step 3), it is necessary to consider the inter-epoch constraint factor; if it is static positioning, the change amount between the previous and current epochs is set to 0, and the variance is infinitesimal.

[0018] Furthermore, when performing the optimization calculation in step 3), it is necessary to consider the inter-epoch constraint factor; if it is dynamic positioning, for the receiver clock error, the previous epoch is set as the initial iteration value, and the variance is infinite; for the tropospheric zenith wet delay component, it is set as a random walk model.

[0019] Furthermore, the PPP observation model established in step 1) is:

[0020]

[0021] where P IF represents the dual-frequency ionosphere-free combined pseudorange observation value; represents the geometric distance between the receiver and the satellite; dT represents the tropospheric zenith wet delay component; dt IF represents the receiver clock error; dt IF represents the satellite clock error; Δ represents other error corrections; ε P_IF represents the pseudorange observation noise; Φ IF represents the dual-frequency ionosphere-free combined carrier phase observation value; λ IF represents the wavelength; represents the true floating-point ambiguity estimation; ε Φ_IFIndicates the carrier phase observation noise.

[0022] Furthermore, other error corrections Δ include earth rotation, relativity, antenna phase winding, and antenna phase center deviation. Brief Description of the Drawings

[0023] Figure 1 is a schematic diagram of the PPP factor graph optimization model of the present invention;

[0024] Figure 2 is a three-dimensional error comparison diagram of JFNG station PPP based on Kalman filtering and factor graph optimization of the present invention;

[0025] Figure 3 is a three-dimensional error comparison diagram of PPP for all 10 stations of the present invention;

[0026] Figure 4(a) is a comparison diagram of the carrier phase residual distribution of JFNG station PPP based on Kalman filtering and factor graph optimization of the present invention;

[0027] Figure 4(b) is a comparison diagram of the pseudorange residual distribution of JFNG station PPP based on Kalman filtering and factor graph optimization of the present invention;

[0028] Figure 5 is a three-dimensional error comparison diagram of dynamic PPP of SUTH station based on Kalman filtering and factor graph optimization of the present invention. Detailed Embodiment

[0029] Considering that the factor graph optimization algorithm has the following advantages: 1) It can use historical, current, and future observations within the smoothing window to estimate the system state. Compared with Kalman filtering, which can only use historical and current information, it makes more full use of the observation information; 2) When the degree of equation nonlinearity is high or the initial state is inaccurate, the initial state is continuously modified and iterated during the graph optimization process, and it is less affected by the nonlinearization error; 3) It supports the plug-and-play of sensors. The centralized filtering algorithm requires the time synchronization of observation values. When there is a delay in the observation values, the processing is relatively complex; 4) It can use the sparsity of the matrix for incremental updates to improve the calculation efficiency. Therefore, the present invention applies factor graph optimization to PPP positioning, elaborates on the ambiguity variable factors, as well as the pseudorange and carrier phase observation factors, and verifies the correctness and effectiveness of the new PPP algorithm based on factor graph optimization using the IGS station data distributed globally. The present invention will be described in detail below with reference to the drawings and embodiments.

[0030] Embodiment of the Precise Point Positioning Method Based on Factor Graph Optimization:

[0031] An embodiment of the precise point positioning method based on factor graph optimization of the present invention is as follows:

[0032] Step 1: Establish a PPP observation model. The state variables in the PPP observation model include three-dimensional position, receiver clock bias, tropospheric zenith wet delay component, and ambiguity in the carrier phase. The specific process is as follows:

[0033] In the multi-frequency undifferenced and uncombined model, the carrier phase and pseudorange observation equations are:

[0034]

[0035] In the formula, P i and Φ i (i = 1, 2, 3, …) respectively represent the pseudorange and carrier phase observations of the i-th frequency; ρ represents the geometric distance between the receiver and the satellite; dt r and dt s respectively represent the receiver clock bias and satellite clock bias in meters; dT represents the slant tropospheric delay; I1 represents the ionospheric delay in the line-of-sight direction of the satellite from the B1 station at the first frequency; is defined as the ratio of the first frequency to the i-th frequency, and a i represents the ionospheric delay scaling factor (relative to the first frequency) in the line-of-sight direction of the satellite from the station at the i-th frequency; D i , D i respectively represent the receiver-side hardware delay and satellite-side hardware delay of the pseudorange observation at the i-th frequency; λ i , N i respectively represent the carrier wavelength and ambiguity at the i-th frequency; Pi and ε Φi respectively represent the observation noise of the pseudorange and carrier phase observations at the i-th frequency; Δ represents other error corrections, including earth rotation, relativity, antenna phase winding, antenna phase center deviation, etc.

[0036] When i = 1, 2, multiply the above undifferenced pseudorange and carrier phase observation equations on the left by the transformation matrix to obtain the most commonly used dual-frequency ionosphere-free combined PPP model:

[0037]

[0038] In the formula, P IF represents the dual-frequency ionosphere-free combined pseudorange observation value; Φ IF represents the dual-frequency ionosphere-free combined carrier phase observation value; ε P_IF represents the pseudorange observation noise; ε Φ_IF represents the carrier phase observation noise.

[0039] In precise point positioning (PPP), IGS precise orbits and clock products are generally used to eliminate satellite orbits and clock errors. In the IGS standard model, satellite clock errors include pseudorange hardware delays. Therefore, in the case of using IGS products, D in the above formula IF will be absorbed by dt r , and D IF will be absorbed by dt s to form new satellite clock error parameters and receiver clock error parameters:

[0040]

[0041] At this time, the ambiguity is:

[0042]

[0043] In the formula, λ IF represents the wavelength; represents the true floating-point ambiguity estimate; N IF represents the theoretical value of the ambiguity.

[0044] Substituting the above expressions into Equation (2) gives:

[0045]

[0046] Equation (5) is the classical ionosphere-free combined PPP function model, where the parameters to be estimated include three-dimensional position parameters, receiver clock error, tropospheric zenith wet delay component, and ambiguity parameters in the carrier phase.

[0047] Step 2: Introduce the factor graph optimization method widely used in SLAM to construct the precise point positioning model, construct the PPP factor graph optimization model, combine the optimization problem with graph theory. The variable nodes in the PPP factor graph optimization model are the state variables composed of the parameters to be estimated in Step 1, the factor nodes are GNSS pseudorange observations and carrier phase observations, and the edges connecting some variable nodes and the corresponding factor nodes are GNSS pseudorange observation errors and carrier phase observation errors. By adjusting and optimizing the state variables to best satisfy the constraints of the edges, the error is minimized, thus obtaining the PPP positioning result with the minimum error (the corresponding error function is essentially Equation (5)). The specific process is as follows:

[0048] Among them, define the state variables of precise point positioning within the sliding window as:

[0049]

[0050] In the formula, n represents the size of the sliding window; m represents the number of ambiguities; x, y, z, dt r , dT represent position, receiver clock error, and tropospheric zenith wet delay component respectively.

[0051] When obtaining GNSS pseudorange and carrier phase observations, its estimation model can be modeled as a factor graph optimization problem. Essentially, it is to solve the maximum a posteriori estimation of the state vector according to the observation information and prior probability, combining Bayesian methods. Assume that all observations are independent and follow a Gaussian distribution with zero mean:

[0052]

[0053] In the formula, {r p , H p} represents the prior information of the state variable, z k represents the k-th carrier phase or pseudorange observation, r represents the residual function constructed according to the carrier phase or pseudorange observation model, and |…| P represents the Mahalanobis distance. It can be seen from the above formula that adding carrier phase or pseudorange observations, that is, increasing the factor nodes of the corresponding observation constraints, the constructed non-linear least squares problem can be iteratively solved by the Levenberg-Marquardt algorithm. The model of PPP positioning based on factor graph optimization is as Figure 1 shown.

[0054] First, the prior information factor in the figure can be obtained by standard single-point positioning solution and constitutes the prior constraint factor. Second, the pseudorange observations and carrier phase observations respectively constitute the pseudorange factor and the carrier phase factor. According to formula (5), the pseudorange factor is a unilateral factor, which is only linked to the position, clock error, and tropospheric variable x k at time k, while the carrier phase factor is a bilateral factor, with both ends linked to x k and the ambiguity variable N m . The core of PPP positioning based on factor graph optimization is the processing of carrier phase ambiguity. When the carrier phase is continuously observed and no cycle slips occur, the ambiguity theoretically remains unchanged. Therefore, multiple carrier phase observation factors can be linked to the same ambiguity variable. When cycle slips occur or new satellites appear, a new ambiguity variable needs to be added. Finally, the construction of the inter-epoch constraint factor is considered for the position, receiver clock error, and tropospheric wet delay respectively. For static positioning, the position remains unchanged between epochs. Therefore, the change amount can be set to zero and the variance can be set to infinitely small. For dynamic positioning, the positions between epochs are different, and the previous epoch is only used as the initial iteration value, with the variance set to infinitely large. For the receiver clock error, its change is difficult to model and is generally set as white noise. Therefore, the previous epoch is only used as the initial iteration value, with the variance set to infinitely large. For the tropospheric wet delay, it is generally set as a random walk model, so the variance can be set to an appropriate noise. After completing the construction of the PPP factor graph, write the corresponding solution program based on the current open-source computing library GTSAM and perform incremental update calculations to obtain the PPP result.

[0055] The method of the present invention will be applied to specific examples below to illustrate the effectiveness of the method of the present invention.

[0056] PPP solution experiment analysis was carried out based on the data of the global IGS observation network. 10 stations were randomly selected, with 24-hour observation data on August 4, 2021 (day of year 216), a sampling rate of 30 s, and each station had 2,880 epochs. In the experiment process, multi-system precise products of GFZ were used, including precise orbits, precise clock offsets, and Earth rotation parameter products. At the same time, error corrections such as antenna phase center, tropospheric dry delay, phase winding, and ocean tide and solid tide corrections were made according to the standard models recommended by IGS. In the PPP solution process, the traditional Kalman filter and factor graph optimization algorithm were used respectively, denoted by Kalman and Factor graph optimization respectively.

[0057] Figure 2 The PPP three-dimensional error comparison of the JFNG station based on Kalman filter and factor graph optimization is given. It can be seen from the figure that there is an obvious accuracy convergence process in the PPP results based on the Kalman filter. This is because in the initial stage, the estimation of the carrier phase ambiguity depends on the accuracy of the pseudorange observations and the estimation is not accurate. As the number of observations accumulates, the ambiguity remains unchanged when no cycle slips occur, and the estimation accuracy gradually improves until it converges to the best accuracy. However, there is no accuracy convergence process in the PPP results of the factor graph optimization. This is because the factor graph optimization algorithm does not only estimate the current epoch state when processing the data of the current epoch, but optimizes all historical and current state variables. When the accuracy of the ambiguity parameter estimation becomes higher and higher with the accumulation of the number of observations, substituting it into the historical observation equation can improve the accuracy of the historical state variables. In this sense, the factor graph optimization algorithm is equivalent to the least squares batch processing algorithm. It should be noted that the factor graph algorithm can fuse all existing historical and current observations to give the optimal estimates of the current and historical epochs.

[0058] To further compare the accuracies of the two algorithms, Figure 3 the result comparison of all stations is given, Figure 3In the figure, two columns are grouped together. Among them, the left column in a group uses Kalman filtering, and the right column uses factor graph optimization. It can be seen from the figure that the mean absolute deviations of PPP in the northeast celestial direction based on Kalman filtering are (0.004m, 0.001m, 0.034m), while the corresponding deviations of factor graph optimization are (0.002m, 0.001m, 0.020m). It can be seen that both the factor graph optimization and Kalman filtering methods can obtain PPP results with millimeter level in the horizontal direction and centimeter level in the vertical direction. However, the factor graph optimization method is significantly more accurate than the Kalman filtering method in the elevation direction. The possible reason is that the vertical component of the position error in PPP positioning has a strong correlation with the receiver clock error and the wet delay in the tropospheric zenith direction. The factor graph optimization method can decouple them using all the existing observation data, while Kalman filtering can only estimate them in a sequential processing manner. Relatively speaking, the factor graph optimization method can better separate the vertical direction error.

[0059] Figures 4(a) and 4(b) show the comparison of PPP residuals of the JFNG station based on Kalman filtering and factor graph optimization. Among them, the darker black color represents factor graph optimization, and the lighter gray color represents Kalman filtering. Theoretically speaking, the smaller the residual and the more concentrated it is near zero and closer to the normal distribution, the more accurate the adjustment system is. It is not difficult to see from the figure that the PPP carrier phase residuals based on factor graph optimization are significantly smaller than the PPP residuals based on Kalman filtering and are also closer to zero. This also proves from another angle that PPP based on factor graph optimization can obtain higher-precision results. It should be noted that no similar phenomenon is found for the pseudorange residuals. This is because, for PPP, the weight of the carrier phase observation value is much larger than that of the pseudorange observation value (at least 100 times), and its contribution to the final PPP accuracy is also greater. The noise of the pseudorange observation value is larger, reaching the order of several meters. The improvement of PPP based on factor graph optimization for accuracy is only at the centimeter level. Therefore, the change of the pseudorange residuals is not obvious.

[0060] In order to further evaluate the performance of PPP based on Kalman filtering and factor graph optimization in a dynamic environment, static IGS station data was selected and solved using the dynamic mode. Figure 5The dynamic PPP three-dimensional errors resolved by two methods at the SUTH station are given respectively. It can be seen from the figure that there is an obvious accuracy convergence process in the PPP results based on Kalman filtering, while there is no accuracy convergence process in the factor graph optimized PPP results, which is consistent with the phenomenon of static PPP. In terms of accuracy, when the convergence is completed (with the criterion that the three-dimensional errors of 10 consecutive epochs are less than 10 cm), the PPP based on factor graph optimization is roughly equivalent to that based on Kalman filtering. The above results once again prove the effectiveness of the PPP method based on factor graph optimization.

[0061] In summary, starting from the GNSS raw observation equation, the present invention constructs a PPP factor graph optimization model according to Bayes' theorem, details the ambiguity variable factor and the pseudo-range and carrier phase observation factors therein, and verifies the correctness and effectiveness of the new PPP algorithm based on factor graph optimization based on the data of IGS stations distributed globally. Compared with the traditional Kalman filtering, the PPP algorithm based on factor graph optimization can achieve comparable accuracy, and has better accuracy in the vertical direction for static positioning. The possible reason is that the data processing mode of batch optimization by factor graph optimization can better separate the vertical direction errors and improve the strong correlation between the vertical direction errors, receiver clock errors and tropospheric zenith wet delays. Another significant feature of the PPP algorithm based on factor graph optimization is that it can give the optimal estimation of all historical epochs in real time, eliminating the convergence process in Kalman filtering. In traditional positioning and attitude systems, such as the current commercial software Inertial Explorer, repeated forward and backward smoothing is required to obtain continuous high-precision pose estimation. The PPP algorithm based on factor graph optimization does not require repeated smoothing and directly gives the optimal estimation of the position, and is expected to replace the data processing method in traditional positioning and attitude, and be widely used in fields such as mobile mapping, photogrammetry, and unmanned driving.

Claims

1. A precise point positioning method based on factor graph optimization, characterized in that It includes the following steps: 1) Establish an ionosphere-free combined PPP observation model, where the state variables in the PPP observation model include three-dimensional position, receiver clock error, tropospheric zenith wet delay component, and ambiguity in the carrier phase; 2) Use the state variables as the variable nodes of the PPP factor graph optimization model, use GNSS pseudorange observation values and carrier phase observation values as the factor nodes of the PPP factor graph optimization model, and use GNSS pseudorange observation errors and carrier phase observation errors as the edges connecting some variable nodes and the corresponding factor nodes, so as to construct a PPP factor graph optimization model: 3) Perform optimization calculations based on the constructed PPP factor graph optimization model to obtain the PPP positioning result with the minimum error.

2. The precise point positioning method optimized based on factor graph according to claim 1, characterized in that, The state variables within the sliding window are expressed as: χ = [x0, x1, …, x n , N0, N1, …, N m ​ x k = [x, y, z, dt r , dT], k ∈ [0, n] where χ represents the state variable; N i represents the ambiguity at the i-th frequency, i = 1, 2, …, m, where m represents the number of ambiguities; x, y, z represent the three-dimensional position; dt r , dT represent the clock error and the tropospheric zenith wet delay component, respectively; n represents the size of the sliding window.

3. The precise point positioning method optimized based on factor graph according to claim 2, characterized in that The factor node corresponding to the GNSS pseudorange observation value is called the pseudorange factor, and the factor node corresponding to the carrier phase observation value is called the carrier phase factor; the pseudorange factor is a unilateral factor and is only linked to the corresponding factor node x k The carrier phase factor is a bilateral factor, and both ends are respectively linked to the corresponding factor node x k and the ambiguity N k linked.

4. The precise point positioning method optimized based on factor graph according to claim 1, wherein, When performing the optimization calculation in step 3), it is necessary to solve the maximum a posteriori estimation of the state vector by combining the Bayesian method according to the observation information and prior probability, where all the observation values are independent of each other and all follow a Gaussian distribution with zero mean.

5. The precise point positioning method optimized based on factor graph according to claim 1, characterized in that When performing the optimization calculation in step 3), it is necessary to consider the epoch-to-epoch constraint factor; if it is static positioning, set the change amount between the previous and current epochs to 0 and the variance to infinitesimal.

6. The precise point positioning method optimized based on factor graph according to claim 1, characterized in that When performing the optimization calculation in step 3), it is necessary to consider the epoch-to-epoch constraint factor; if it is dynamic positioning, for the receiver clock error, set the previous epoch as the initial iteration value and the variance to infinite; for the tropospheric zenith wet delay component, set it as a random walk model.

7. The precise point positioning method optimized based on factor graph according to claim 1, wherein The ionosphere-free combined PPP observation model established in step 1) is: In the formula, P IF represents the dual-frequency ionosphere-free combined pseudorange observation value; ρ represents the geometric distance between the receiver and the satellite; dT represents the tropospheric zenith wet delay component; dt IF represents the receiver clock error; dt IF represents the satellite clock error; Δ represents other error corrections; ε P_IF represents the pseudorange observation noise; Φ IF represents the dual-frequency ionosphere-free combined carrier phase observation value; λ IF represents the wavelength; represents the true floating ambiguity estimate; ε Φ_IF represents the carrier phase observation noise.

8. The precise point positioning method optimized based on factor graph according to claim 7, characterized in that Other error corrections Δ include earth rotation, relativity, antenna phase winding, and antenna phase center deviation.

Images