An RTK positioning method based on integer constrained Gaussian filtering

Abstract

The invention discloses an RTK positioning method based on integer constrained Gaussian filtering, comprising the following steps: step 1, obtaining mobile station data: obtaining original observation data from a GNSS receiver; step 2, obtaining base station data: obtaining differential observation data of a base station or a virtual CORS station from a communication link; step 3, data preparation: performing single difference and double difference processing on the observation data obtained in the above steps 1 and 2 by calculating the current position and clock difference of a satellite through satellite broadcast ephemeris, so as to obtain a double difference measurement equation after processing; step 4, establishing a system model; step 5, integer constrained progressive Gaussian filtering; step 6, updating a state; the method of the invention proposes an RTK positioning method based on integer constrained progressive Gaussian filtering by constructing a progressive Gaussian filtering framework under integer constraints, combining the advantages of Bayesian reasoning and homology methods, and gradually fusing information by constructing a progressive step length, thereby avoiding the influence of the initial value selection problem on the LAMBDA algorithm, and reducing the influence of a bad floating point solution on ambiguity fixation, thereby further improving the success rate of ambiguity fixation.

Description

Technical Field

[0001] The present invention relates to the technical field of high-precision satellite positioning, and in particular to an RTK positioning method based on integer constrained Gaussian filtering. Background Art

[0002] GNSS-RTK (Global Navigation Satellite System-Real-Time Kinematic) technology has broad application prospects in surveying and mapping, unmanned driving, agriculture and precision engineering due to its high-precision positioning capability. This technology can provide real-time positioning within the centimeter-level accuracy range by using the differential signal between the ground reference station and the mobile terminal. However, although GNSS-RTK technology performs well under ideal conditions, it still faces many challenges when the satellite signal is blocked or interfered, especially in the environment of high-rise buildings in cities, forest-covered areas and areas with dense interference sources, where RTK filtering is prone to divergence. Satellite signal blocking leads to a reduction in the number of effective satellites or an increase in the geometric dilution of precision (GDOP), which in turn affects the positioning accuracy. In addition, the multipath effect of the signal and interference sources (such as radio frequency interference) will lead to an increase in pseudo-range measurement errors. Traditional filtering algorithms are difficult to cope with this uncertainty, causing drastic fluctuations in pseudo-range positioning results. RTK positioning is sensitive to the selection of the initial value of the carrier phase. The above interference will seriously affect the solution accuracy of the filter, and poor floating-point solutions will make it difficult to fix the ambiguity. In addition, the accuracy of GNSS-RTK positioning is greatly affected by the system model, but it is difficult to accurately model the satellite measurement noise in different scenarios, and it is difficult to obtain the statistical characteristics of satellite measurement noise in different scenarios. Currently, there is no technology in GNSS-RTK that can simply and effectively solve the high-precision positioning problem caused by measurement noise uncertainty.

[0003] In order to overcome the problems of low fixation rate and poor robustness of existing GNSS-RTK positioning methods, the present invention proposes an RTK positioning method based on integer constrained Gaussian filtering. By combining Bayesian reasoning with the homotopy method, a progressive Gaussian filtering framework under integer constraints is derived, which effectively improves the fixation rate, reduces the RMS value, and enhances the system robustness to solve the problems existing in the prior art. Summary of the invention

[0004] The present invention solves the problem of low fixation rate and poor robustness of RTK positioning when GNSS signals are interfered or blocked in complex environments. The present invention can improve the fixation rate of ambiguity fixation under the above conditions, thereby improving the accuracy and robustness of RTK positioning.

[0005] The technical solution adopted in the present invention is as follows:

[0006] An RTK positioning method based on integer constrained Gaussian filtering comprises the following steps:

[0007] Step 1, obtain mobile station data: obtain raw observation data from the GNSS receiver, including but not limited to pseudorange, carrier phase, Doppler, signal-to-noise ratio, satellite broadcast ephemeris data, etc.

[0008] Step 2, obtain the base station data: obtain the differential observation data of the base station or virtual CORS station from the communication link, including but not limited to pseudorange, carrier phase, signal-to-noise ratio, reference station location information, etc.

[0009] Step 3, data preparation: Calculate the current position and clock error of the satellite through the satellite broadcast ephemeris, and perform single-difference and double-difference processing on the observation data obtained in the above steps 1 and 2 to obtain the processed double-difference measurement equation.

[0010] Step 4: Establish a system model: Assume that a total of m satellites are observed at time k, and satellite i is used as the reference satellite. The system model is described as follows:

[0011]

[0012] In the formula, x k =[r k a k ] T is the system state vector at time k; r k is the carrier position; a k is the single difference ambiguity between stations; w k is the process noise; v k To measure noise; is the double difference measurement vector, h(x k ) is a nonlinear measurement function, expressed as:

[0013]

[0014] Step 5: Integer Constrained Progressive Gaussian Filtering

[0015] The asymptotic process of prior distribution-posterior distribution is established, the pseudo time constant τ is introduced, and the integer asymptotic Gaussian filtering framework is given as follows:

[0016]

[0017] The framework includes time update, progressive measurement update and progressive ambiguity update processes.

[0018] Where p(·) represents the probability density function, Z 1:k ={z 1 ,z 2 ,...,z k} is the set of measurements at time k, τ i (i∈[1,N]) is a pseudo-time constant, and τ 0 =0 (indicates that there is no pseudo-time concept here), Where N is the total number of asymptotic steps, Δ i =1 / N is the progressive length of the i-th step;

[0019] Step 5.1, progressive Gaussian filtering

[0020] Based on the Bayesian update process, through the pseudo time constant τ, the following relationship is established:

[0021]

[0022] In the formula,

[0023] Furthermore, the pseudo-time τ i The system position update process is shown as follows:

[0024]

[0025] In the formula, a k represents the fuzziness parameter, p(a k ,τ i |Z 1:k ) represents the normalized distribution formed in its asymptotic process.

[0026] Step 5.2, Asymptotic Gaussian Integer Constraint

[0027] Furthermore, the pseudo-time τ i-1 The asymptotic posterior probability density function at is p(x k ,τ i-1 |Z 1:k ) and obeys Gaussian distribution, then:

[0028]

[0029] In the formula represents the pseudo time τ i-1 The asymptotic posterior estimate of is its corresponding covariance matrix.

[0030] Furthermore, the current asymptotic posterior estimated floating-point ambiguity is rounded to obtain the double-difference asymptotic posterior estimated integer ambiguity. If the asymptotic posterior estimated integer ambiguity satisfies the corresponding conditional test, the integer constraint is successful.

[0031] Step 5.3, integer constrained progressive filtering cutoff condition

[0032] Considering equation (1), the nonlinear measurement z k exist Taylor expansion:

[0033]

[0034] in, represents the estimation error, It means to take the derivative of nonlinear function and ignore the high-order terms:

[0035]

[0036] In the formula, represents the asymptotic posterior estimate residual;

[0037] In step 5.2, if the ambiguity has been successfully fixed, the iteration process can be terminated and step 6 will be executed. However, if the ambiguity has not been successfully fixed, the iteration process will continue until the step length is completed. If the step length has been completed, step 6 will be executed, otherwise jump to step 5.1 to continue the iterative calculation. To prevent the system from overestimating, a cutoff condition is introduced and Among them, ε is a very small constant. If the cutoff condition is met at the same time and Prove that the current solution is the optimal solution under the non-fixed state and execute step 6. Otherwise, jump to step 5.1 to continue iterative calculation.

[0038] Step 6, Update Status

[0039] If the asymptotic Gaussian integer constraint in step 5.2 is successful, the corresponding integer solution can be obtained, otherwise only a floating-point solution can be obtained.

[0040] The beneficial effects of the present invention are as follows:

[0041] The method of the present invention constructs a progressive Gaussian filtering framework under integer constraints, combines the advantages of Bayesian reasoning and homology methods, and proposes an RTK positioning method based on integer constrained progressive Gaussian filtering. By constructing a progressive step length to gradually fuse information, the impact of the initial value selection problem on the LAMBDA algorithm is avoided, while the impact of poor floating-point solutions on ambiguity fixation is reduced, further improving the success rate of ambiguity fixation. The positioning performance can be significantly improved in complex environments such as urban canyons and viaducts. Through the above content, the method of the present invention can improve the proportion of fixed solutions and the positioning accuracy of non-fixed solutions in real-time positioning compared with general methods. BRIEF DESCRIPTION OF THE DRAWINGS

[0042] Figure 1 is the overall system block diagram of the present invention;

[0043] Figure 2It is a flow chart of the RTK positioning method with progressive Gaussian filtering integer constraints of the present invention;

[0044] Figure 3 This is the vehicle test track of Example 1 of the present invention. DETAILED DESCRIPTION

[0045] In order to deepen the understanding of the present invention, the present invention will be further described in detail below in conjunction with examples. The examples are only used to explain the present invention and do not constitute a limitation on the protection scope of the present invention.

[0046] It should be noted that, in this article, relational terms such as first and second, etc. are only used to distinguish one entity or operation from another entity or operation, and do not necessarily require or imply any such actual relationship or order between these entities or operations. Moreover, the terms "include", "comprise" or any other variants thereof are intended to cover non-exclusive inclusion, so that a process, method, article or device including a series of elements includes not only those elements, but also other elements not explicitly listed, or also includes elements inherent to such process, method, article or device.

[0047] Embodiment 1

[0048] The present invention proposes an RTK positioning method based on integer constrained Gaussian filtering, and the specific steps of the method are as follows:

[0049] Step 1: Obtain raw observation data from the GNSS receiver, including but not limited to pseudorange, carrier phase, Doppler, and signal-to-noise ratio.

[0050]

[0051] Where P and Φ represent pseudorange observation and carrier phase observation respectively, and their subscripts represent satellite number i and receiver number r respectively; represents the Euclidean distance between the satellite and the ground, where r i and r r They represent the three-dimensional coordinate vectors of the satellite and the receiver respectively. The satellite position is calculated using the ephemeris; c represents the speed of light; dt r and dt i Represent the clock differences of the receiver and the satellite respectively; It represents the tropospheric delay error between satellite and ground; represents the ionospheric delay error between the satellite and the ground; λ and a represent the wavelength and ambiguity of the carrier phase observation respectively; ε P and ε Φ represents the unmodeled error term, including multipath effects and receiver measurement noise.

[0052] Step 2: Obtain differential observation data of the reference station or virtual CORS station from the communication link, including but not limited to pseudorange, carrier phase, Doppler, and signal-to-noise ratio.

[0053] Step 3, performing single-difference and double-difference processing on the observation data obtained in the above steps 1 and 2 to obtain the processed double-difference measurement equation.

[0054]

[0055] In the formula, and They represent the double difference measurement of pseudorange and carrier phase, and represents the single-difference ambiguity, and represents the unmodeled error after double difference processing.

[0056] Step 4: Establish a system model. Assume that m satellites are observed at time k, and satellite i is used as the reference satellite. The system model is described as follows:

[0057]

[0058] In the formula, x k =[r k a k ] T is the system state vector at time k; r k is the carrier position; a k is the single difference ambiguity between stations; w k is the process noise; v k To measure noise; is the double difference measurement vector, h(x k ) is a nonlinear measurement function, expressed as:

[0059]

[0060] Step 5: Integer Constrained Progressive Gaussian Filtering

[0061] The asymptotic process of prior distribution-posterior distribution is established, the pseudo-time constant τ is introduced by using Bayes' theorem, and the integer asymptotic Gaussian filtering framework is given as follows:

[0062]

[0063] This framework includes time update, gradual measurement update and gradual ambiguity update processes. Where p(·) represents the probability density function, Z 1:k ={z 1 ,z 2 ,...,z k} is the set of measurements at time k, τi (i∈[1,N]) is a pseudo-time constant, and τ 0 =0, Where N is the total number of asymptotic steps, Δ i =1 / N is the progressive length of the i-th step; p(z k ,Δ i |x k ) is the asymptotic likelihood function, which is defined as:

[0064]

[0065] Step 5.1, progressive Gaussian filtering

[0066] Based on the Bayesian update process, through the pseudo time constant τ, the recursive relationship can be established as follows:

[0067]

[0068] In the formula, C k (τ i-1 )=1 / ∫p(z k ,Δ i |x k )p(x k ,τ i-1 |Z 1:k-1 )dx k Then the pseudo time τ i The asymptotic posterior probability density function at is:

[0069]

[0070] In the formula,

[0071] Correspondingly, the pseudo-time τ i The system position update process is shown as follows:

[0072]

[0073] In the formula, a k represents the fuzziness parameter, p(a k ,τ i |Z 1:k ) represents the normalized distribution formed in its asymptotic process.

[0074] Step 5.2, Asymptotic Gaussian Integer Constraint

[0075] Assume the pseudo time τ i-1 The asymptotic posterior probability density function at is p(x k ,τ i-1 |Z 1:k ) and obeys Gaussian distribution, then:

[0076]

[0077] In the formula represents the pseudo time τ i-1 The asymptotic posterior estimate of As its corresponding covariance matrix, the single difference ambiguity parameters in the state vector are converted into double difference ambiguity parameters using the conversion matrix G, and we can get:

[0078]

[0079] In the formula, represents the state vector after conversion, represents the position vector, represents the double difference ambiguity vector, and the matrix G is the conversion matrix from single difference to double difference, which is expressed as follows:

[0080]

[0081] Status The corresponding covariance matrix is ​​expressed as:

[0082]

[0083] The corresponding state floating point solution has been obtained And the corresponding variance covariance matrix Now we constrain it to integers. The joint probability density function of position and ambiguity is expressed as follows

[0084]

[0085] Therefore, the integer ambiguity is obtained by solving the following integer least squares problem:

[0086]

[0087] In the formula, represents the double difference integer ambiguity vector to be solved. This optimization problem can be solved by the LAMBDA algorithm. First, the double difference ambiguity covariance matrix Perform descending correlation processing and obtain the optimal integer solution vector through Ratio test. It can be simplified as follows:

[0088]

[0089] If the corresponding Ratio test value is met, the LAMBDA algorithm solves successfully, otherwise the LAMBDA algorithm fails to solve and the integer constraint fails.

[0090] Step 5.3, integer constrained progressive filtering cutoff condition

[0091] Considering equation (1), the nonlinear measurement z k exist Taylor expansion:

[0092]

[0093] in represents the estimation error, Derivative of nonlinear function, ignoring high-order terms, we can get

[0094]

[0095] In the formula, represents the asymptotic posterior estimate residual;

[0096] In step 5.2, when the ambiguity has been successfully fixed, the iteration process can be terminated, that is, Execute step 6. [·] indicates the rounding process. However, if the ambiguity cannot be fixed successfully, the iterative process will continue until the step length is completed. If the step length has been completed, execute step 6, otherwise jump to step 5.1 to continue iterative calculation. To prevent the system from overestimating, the cutoff condition is introduced and Among them, ε is a very small constant. If the cutoff condition is met at the same time and Prove that the current solution is the optimal solution under the non-fixed state and execute step 6. Otherwise, jump to step 5.1 to continue iterative calculation.

[0097] Step 6, Update Status

[0098] If the LAMBDA algorithm solves successfully in step 5.2, the corresponding double difference integer solution can be obtained The corresponding single-difference integer ambiguity vector It can be recovered by solving the double difference integer ambiguity vector. And it can be obtained by the double difference integer solution Modified integer estimate and the corresponding covariance matrix P k|k,AR Otherwise only floating point blur is obtained. and the corresponding floating point solution and the corresponding covariance matrix P k|k .

[0099] The test was conducted on June 28, 2023 in Pudong New Area, Shanghai. It used the original FSS-NAV619 navigation equipment with a sampling frequency set to 1Hz. The RTK module was a UM982 board that supports BDS and GPS full-band, GLONASSL1 / L2 band, GalileoE1 / E5b band, and can record RTK raw observation data. The sampling rate of the GNSS receiver was 1Hz, and a total of 1351s of measured data were collected. The satellite cutoff elevation angle was set to 15 degrees. The progressive step length (total number of progressive steps) was set to N=10. The reference value uses the dual-frequency multi-GNSS-RTK / INS data processing results of the Wuhan Maipu space-time integrated navigation system data processing software, with horizontal and vertical accuracies of 0.02m and 0.03m respectively. The experimental trajectory is as follows Figure 3 As shown in the figure, the experimental environment includes open satellite observation environment, as well as complex observation environments such as viaducts and urban roads. In the first half of the experiment, a small number of satellites were severely blocked, and in the second half, signal blockages were frequent. The number of visible satellites varied relatively widely, with the minimum number of visible satellites being only 3.

[0100] In order to verify the improvement of the positioning system performance by the proposed method, the proposed method is compared with the extended Kalman filter (EKF), iterative EKF (IEKF), robust EKF (REKF) and iterative REKF (IREKF), and the root mean square error (RMSE) and ambiguity fixation rate in the horizontal and vertical coordinate directions are calculated respectively, and recorded in Table 1.

[0101] Table 1 Example - RMSE and fixation rate

[0102] method EKF IREKF IEKF REKF Proposed method promote RMSE-level (m) 0.8718 1.0223 0.9600 0.9075 0.6696 23.19％ RMSE-vertical (m) 0.2959 0.6174 0.5150 0.2871 0.2062 28.18％ Fixing rate (%) 50.33 61.58 58.48 71.13 90.38 19.25％

[0103] Obviously, the method proposed in the present invention maximizes the use of the effective information carried by the data and has significant performance improvements compared to other methods.

[0104] Finally, it should be noted that the above is only a preferred embodiment of the present invention and is not intended to limit the present invention. Although the present invention has been described in detail with reference to the aforementioned embodiments, those skilled in the art can still modify the technical solutions described in the aforementioned embodiments or replace some of the technical features therein by equivalents. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present invention shall be included in the protection scope of the present invention.

Claims

1. An RTK positioning method based on integer constrained Gaussian filtering, characterized in that: The following steps are involved: Step 1, obtain mobile station data: obtain raw observation data from the GNSS receiver, including pseudorange, carrier phase, Doppler, signal-to-noise ratio, and satellite broadcast ephemeris data; Step 2, obtain the base station data: obtain the differential observation data of the base station or virtual CORS station from the communication link, including pseudorange, carrier phase, signal-to-noise ratio, and reference station location information; Step 3, data preparation: calculate the satellite's current position, clock error, clock drift, and ionosphere estimation value through satellite broadcast ephemeris, and correct the error; perform single-difference and double-difference processing on the observation data obtained in the above steps 1 and 2 to obtain the processed double-difference measurement equation; Step 4: Establish a system model: Assume that a total of m satellites are observed at time k, and satellite i is used as the reference satellite. The system model is described as follows: In the formula, x k =[r k a k ] T is the system state vector at time k; r k is the carrier position; a k is the single difference ambiguity between stations; F k|k-1 is the state transfer matrix; w k is the process noise; v k To measure noise; is the double difference measurement vector; is the corresponding double difference pseudorange; is the corresponding double difference carrier; h(x k ) is a nonlinear measurement function, expressed as: Where λ is the carrier wavelength of the corresponding frequency; is the single-difference ambiguity between reference satellite stations; is the single difference ambiguity between satellite stations m; Step 5, integer constrained progressive Gaussian filtering: The asymptotic process of prior distribution-posterior distribution is established, the pseudo time constant τ is introduced, and the integer asymptotic Gaussian filtering framework is given as follows: In the formula, p(·) represents the probability density function, Z 1:k ={z1,z2,...,z k } is the set of measurements at time k, τ i (i∈[1,N]) is a pseudo time constant, and τ0=0, Where N is the total number of asymptotic steps, Δ i =1 / N is the progressive length of the i-th step; Step 6, update the status.

2. The RTK positioning method based on integer constrained Gaussian filtering according to claim 1, characterized in that: The specific process of step 5 is as follows: Step 5.1, progressive Gaussian filtering Based on the Bayesian update process, through the pseudo time constant τ, the following relationship is established: In the formula, Pseudo-time τ i The location update process is shown as follows: In the formula, p(a k ,τ i |Z 1:k ) represents the normalized distribution formed in its asymptotic process; Step 5.2, Asymptotic Gaussian Integer Constraint Pseudo-time τ i-1 The asymptotic posterior probability density function at is p(x k ,τ i-1 |Z 1:k ) and obeys Gaussian distribution, then: In the formula represents the pseudo time τ i-1 The asymptotic posterior estimate of is its corresponding covariance matrix; Indicates that it satisfies the normal distribution; The single-difference ambiguity parameters in the current asymptotic posterior estimation are converted into double-difference ambiguity parameters using the conversion matrix, and the double-difference ambiguity is obtained. After the calculation is completed, the double difference ambiguity The double difference asymptotic posterior estimated integer ambiguity is obtained by rounding; if the asymptotic posterior estimated integer ambiguity satisfies the corresponding condition test, the integer constraint is successful; Step 5.3, integer constrained progressive filtering cutoff condition Combining equation (1), the nonlinear measurement z k exist Taylor expansion: in represents the estimation error, Represents the derivative of a nonlinear function; Ignoring higher-order terms, we get: In the formula, represents the asymptotic posterior estimate residual; When the ambiguity has been successfully fixed in step 5.2, the iteration process can be terminated and step 6 can be executed; If the ambiguity cannot be fixed successfully, the iterative process will continue until the step length is completed. If the step length has been completed, execute step 6, otherwise jump to step 5.1 to continue iterative calculation; To prevent overestimation, a cutoff condition is introduced and Among them, ε is a custom threshold; if the cutoff conditions are met at the same time and To prove that the solution is the optimal solution under the non-fixed state, execute step 6, otherwise jump to step 5.1 to continue iterative calculation.

3. The RTK positioning method based on integer constrained Gaussian filtering according to claim 1, characterized in that: The step 6 is specifically as follows: if the asymptotic Gaussian integer constraint in step 5.2 is successful, the corresponding integer solution can be obtained, otherwise only a floating point solution can be obtained.

4. The RTK positioning method based on integer constrained Gaussian filtering according to claim 1, characterized in that: The errors in step 3 include the earth's rotation error and the relativistic effect.

5. The RTK positioning method based on integer constrained Gaussian filtering according to claim 1, characterized in that: In step 4, the satellite with the highest elevation angle is used as the reference satellite, and the other satellites are subtracted from the reference satellite once to obtain an inter-satellite single difference. The observation data of the same satellite between the reference station and the mobile station are subtracted to obtain an inter-station single difference. The inter-satellite single difference is subtracted from the inter-station single difference to obtain a double difference observation value.

Images