RTK positioning vertical protection level optimization method based on improved double-difference subset

Abstract

The invention provides an RTK positioning vertical protection level optimization method based on an improved double-difference subset, and relates to the technical field of satellite navigation positioning. According to the method, firstly, a GNSS inter-station single-difference observation equation and a GNSS station-satellite double-difference observation equation are established, error correction is introduced to reduce errors of a long baseline positioning mode, and to-be-estimated state quantity parameters of a GNSS receiver terminal are updated; then, the GNSS inter-station single-difference observation equation and the GNSS station-satellite double-difference observation equation are processed, and a GNSS station-satellite double-difference observation matrix is obtained; dividing each group of observed quantities in the GNSS station satellite double-difference observation matrix into a subset; and finally, discharging fault data of each subset, and calculating a vertical protection level of a terminal positioning solution of each epoch receiver through an MHSS algorithm. According to the method, ionosphere delay and troposphere delay errors can be reduced, and the precision of observed quantity and the overall reliability of a receiver terminal positioning result are improved.

Description

Technical Field

[0001] The invention relates to the technical field of satellite navigation and positioning, and in particular to an RTK positioning vertical protection level optimization method based on an improved double difference subset. Background Art

[0002] In recent years, autonomous vehicles have shown great potential in various market applications with their high efficiency and low traffic accident rate. At present, the Global Navigation Satellite System (GNSS) is the most important navigation system for ground vehicles, which can provide users with all-day, all-weather and high-precision navigation positioning and timing services. High-precision positioning technologies such as Precise Point Positioning (PPP) technology and Real Time Kinematic (RTK) positioning technology are becoming more and more mature, and RTK positioning technology has received more and more attention in the field of autonomous driving because it can converge quickly under ideal conditions and can achieve dynamic centimeter-level and static millimeter-level positioning accuracy.

[0003] RTK positioning technology uses carrier phase differential to eliminate the error related to the propagation path between the mobile station and the reference station within a certain range (within 10km) to achieve high-precision relative positioning. At present, the research on the integrity monitoring method of carrier phase differential positioning includes the fault detection and identification algorithm of the pseudo-range domain and the carrier phase domain of the GNSS / INS combined positioning, the ambiguity reliability test method based on integrity, the method of monitoring the integrity of the carrier phase using the single-difference carrier phase measurement value, and the RTK integrity assessment method based on the risk probability decomposition of the system end. The above methods have made certain innovative work in the relevant aspects of carrier phase integrity monitoring, but still do not involve the field of subset construction. Therefore, based on the double-difference observation, the present invention proposes a RTK positioning vertical protection level optimization method based on the improved double-difference subset. By constructing a subset based on the double-difference observation equation, extracting the relevant position information in the double-difference design matrix as the input parameter, and bringing it into the Multiple Hypothesis Solution Separation (MHSS) algorithm, multiple fault hypothesis analysis and vertical protection level calculation are performed to provide a stronger positioning accuracy and reliability guarantee for the navigation system of the autonomous vehicle. Summary of the invention

[0004] The technical problem to be solved by the present invention is to provide an RTK positioning vertical protection level optimization method based on an improved double difference subset in view of the shortcomings of the above-mentioned prior art. This method has certain limitations in the application research of current RTK positioning technology and advanced receiver autonomous integrity monitoring technology in theoretical models and technical routes, and there is still a large gap from achieving the goal of integrity monitoring of carrier differential positioning. The method of the present invention studies the carrier differential positioning algorithm and integrity monitoring technology, relies on the global satellite navigation system, combines the Advanced Receiver Autonomous Integrity Monitoring (ARAIM) algorithm of MHSS, monitors the RTK positioning process through RTK positioning vertical protection level calculation, improves positioning accuracy and reliability, and realizes integrity monitoring of RTK positioning.

[0005] In order to solve the above technical problems, the technical solution adopted by the present invention is: an RTK positioning vertical protection level optimization method based on an improved double difference subset, comprising the following steps:

[0006] Step 1: Adjust the existing RTK high-precision positioning model parameters, introduce error correction to reduce the positioning error caused by ionospheric delay and tropospheric delay in the long baseline case, and update the parameters of the state quantity to be estimated of the GNSS receiver terminal;

[0007] Step 1.1: Establish GNSS non-difference observation equations to obtain the functional relationship between the original observations of the GNSS receiver terminal and the unknown parameters in the satellite navigation and positioning system;

[0008] The original observations of the GNSS receiver terminal include pseudorange observations and carrier phase observations;

[0009] At a certain moment, the satellite s observed by the GNSS receiver terminal r is at a frequency f i Pseudorange observations on and carrier phase observation Respectively expressed as:

[0010]

[0011] in, is the geometric distance between satellite s and receiver terminal r, c is the light speed of satellite s, dt r and dt s are the clock errors of the GNSS receiver terminal r and the satellite s, is the ionospheric delay between satellite s and GNSS receiver terminal r, which is proportional to frequency f i related, is the tropospheric delay between the satellite s and the GNSS receiver terminal r, is the wavelength of satellite s, For GNSS receiver terminal r at frequency f i The carrier phase integer ambiguity on , ε and ξ are pseudorange measurement error and carrier phase measurement error respectively;

[0012] Step 1.2: Differencing the GNSS non-differenced observation equation in step 1.1 to obtain the single-differenced observation equation between GNSS stations, thereby eliminating the satellite clock error at the satellite end;

[0013] GNSS receiver terminal r and GNSS receiver terminal b observe the same satellite s at the same time, and then use the observations of satellite s by GNSS receiver terminal r and GNSS receiver terminal b as the single difference between GNSS stations, that is, make a difference to the GNSS non-difference observation equation in step 1.1, and obtain the single difference observation equation between GNSS stations:

[0014]

[0015] in, is the single difference operator between GNSS stations, is the satellite s observed by GNSS receiver terminal r and GNSS receiver terminal b at frequency f i The single-difference observations between pseudo-range stations on is the satellite s observed by GNSS receiver terminal r and GNSS receiver terminal b at frequency f i The carrier phase single-difference observation between stations on is the geometric distance difference between GNSS receiver terminal r and GNSS receiver terminal b and satellite s; dt rb is the clock difference between GNSS receiver terminal r and GNSS receiver terminal b; is the ionospheric delay difference between GNSS receiver terminal r and GNSS receiver terminal b and satellite s; is the tropospheric delay difference between GNSS receiver terminal r and GNSS receiver terminal b and satellite s; For GNSS receiver terminal r and GNSS receiver terminal b at frequency f i The inter-station single-difference carrier phase integer ambiguity on ;

[0016] Step 1.3: Based on the single-difference observation equation between GNSS stations, GNSS receiver terminal r and GNSS receiver terminal b simultaneously observe satellite s and satellite k, and establish a long-baseline GNSS station-satellite double-difference observation equation between satellite s and satellite k to eliminate the error of GNSS receiver terminal r, and obtain the state quantity parameters to be estimated of GNSS receiver terminal r in the long-baseline case;

[0017] The long baseline GNSS station double difference observation equation is:

[0018]

[0019] in, is the GNSS station double difference operator, is the satellite s and satellite k observed by receiver terminal r and receiver terminal b at frequency f i The pseudo-range station-satellite double difference observations on; is the satellite s and satellite k observed by receiver terminal r and receiver terminal b at frequency f i The carrier phase station-satellite double difference observations on the is the geometric distance difference between receiver terminal r and receiver terminal b and satellite s and satellite k; is the ionospheric delay difference between receiver terminal r and receiver terminal b and satellite s and satellite k; is the tropospheric delay difference between receiver terminal r and receiver terminal b and satellite s and satellite k; For receiver terminal r and receiver terminal b at frequency f i The carrier phase station-satellite double difference integer ambiguity on the

[0020] When the distance between GNSS receiver terminal r and GNSS receiver terminal b is less than the set threshold, the tropospheric and ionospheric errors are ignored, so the short baseline GNSS station double difference observation equation is rewritten as:

[0021]

[0022] In the short baseline RTK positioning mode, the data of L1 and L2 frequency bands in the global satellite navigation system are received by the GNSS receiver terminal r, and the position parameters of the GNSS receiver terminal r and the integer ambiguity parameters of the single-difference carrier phase between stations are estimated. Then the state quantity x to be estimated of the GNSS receiver terminal r is expressed as:

[0023]

[0024] Among them, r r is the position parameter of the GNSS receiver terminal r, v r is the velocity parameter of the GNSS receiver terminal r, a rb,1 and a rb,2 They are the single-difference carrier phase integer ambiguity parameters between stations in the L1 and L2 frequency bands respectively;

[0025] When the distance between GNSS receiver terminal r and GNSS receiver terminal b is greater than or equal to the set threshold, the tropospheric wet delay and ionospheric delay of the long baseline are taken as unknowns for parameter estimation, and the ionosphere and troposphere state quantities are added to the Kalman filter state quantity, and the time changes of the two are simulated as a random walk model;

[0026] The long baseline GNSS station double difference observation equation is further obtained as:

[0027]

[0028] Then, the estimated state parameter x1 of the GNSS receiver terminal r in the long baseline case is:

[0029]

[0030] Among them, Z T,r and Z T,b are the total delays in the tropospheric direction of GNSS receiver terminal r and GNSS receiver terminal b, respectively, is the inter-station double difference of the ionospheric delay error in the zenith direction on the L1 band, k is the reference satellite, and m is a different satellite system;

[0031] Step 1.4: Use the Kalman filter process to update the estimated state quantity parameters of the GNSS receiver terminal r;

[0032] The covariance matrix corresponding to the state quantities of the ionosphere and troposphere is introduced into the Kalman filter process. In the Kalman filter prediction step, the state estimation value of the GNSS receiver terminal r at the previous moment is used. and its covariance matrix P k-1|k-1 , predict the state estimate of the GNSS receiver terminal r at the current moment and its covariance matrix P k|k-1 :

[0033]

[0034] Among them, Q k is the predicted noise covariance matrix of GNSS receiver terminal r, F k is the state transfer matrix of the GNSS receiver terminal r;

[0035] The predicted noise covariance matrix Q of the GNSS receiver terminal r is k The covariance matrix of the tropospheric error and the ionospheric error is included, and the state transfer matrix F k To spread;

[0036] In the Kalman filter update step, the observation z is used kTo update the parameters of the state quantity to be estimated; the Kalman filter process uses the observation matrix H k and the observation noise covariance matrix R k Update the state of the state to be estimated:

[0037]

[0038] Among them, K k is the Kalman gain, H k T H k The transposed matrix of is the updated state quantity to be estimated;

[0039] The covariance matrix P of the GNSS receiver terminal r after the estimated state quantity is updated k As shown in the following formula:

[0040] P k =(IK k H k ) k|k-1

[0041] Where I is the identity matrix;

[0042] Step 2: Use the GNSS receiver to obtain the GNSS raw observation data of different satellite navigation systems in real time, and screen the GNSS raw observation data, remove abnormal data and conduct quality assessment on the remaining data;

[0043] The different satellite navigation systems include four satellite navigation systems: GPS, BDS, GLONASS and Galileo;

[0044] The GNSS raw observation data includes carrier phase, pseudorange and satellite coordinate information;

[0045] Step 3: Process the GNSS station-to-station single-difference observation equation and the long baseline GNSS station-to-satellite double-difference observation equation given in step 1 to obtain the GNSS station-to-satellite double-difference observation matrix; wherein each set of observations in the GNSS station-to-satellite double-difference observation matrix is ​​a subset, and the GNSS original observation data obtained in step 2 is used for testing and verification;

[0046] Process the single-difference observation equations between GNSS stations and the double-difference observation equations between long-baseline GNSS stations and satellites, and construct the double-difference observation matrix between GNSS stations and satellites;

[0047] The single-difference observation equation between GNSS stations is rearranged and expressed as: V1 = A1·δx1 using the prior residual matrix V1, the single-difference observation matrix A1 between GNSS stations, and the single-difference state quantity to be estimated between GNSS stations δx1.

[0048]

[0049] Among them, the elements in V1 They are different system frequencies f i The satellite pseudorange and carrier phase residuals are (·) for different satellite systems, G for GPS, C for BDS, The frequencies of different satellite navigation systems are f i The satellite direction cosines, and l, m, n are the components of the direction cosines on the x, y and z axes respectively. The dx, dy, dz contained in the estimated state δx1 are the components of the position correction in the x, y and z axes respectively. c·δt is the influence of the clock error. The influence of clock errors between different satellite systems and the GPS system. The frequency of the jth satellite is f i The fuzziness parameter of

[0050] Then select a reference satellite as the benchmark satellite, subtract the observation value of the benchmark satellite from the observation value of the same frequency of the same satellite navigation system, and obtain the new GNSS station-satellite double difference observation equation: V = A·δx, where V is the GNSS station-satellite double difference residual, A is the GNSS station-satellite double difference observation matrix, δx is the GNSS station-satellite double difference state quantity to be estimated,

[0051]

[0052] Among them, p (·)j is the pseudorange residual of the jth satellite, φ (·)j is the carrier phase residual of the jth satellite, N (·)j is the ambiguity parameter of the jth satellite, l (·)n ,m (·)n ,n (·)n are the direction cosines of satellites in different systems;

[0053] The GNSS station-satellite double-difference observation matrix A contains the geometric relationship information between the GNSS receiver and the satellite. Each set of observations in the GNSS station-satellite double-difference observation matrix is ​​obtained by performing differential operations on the GNSS station-to-station single-difference observations of a reference satellite and an observation satellite. By combining different satellites, multiple different subsets are formed, each of which corresponds to a different satellite configuration.

[0054] Step 4: Exclude the fault data of each subset and calculate the vertical protection level of the receiver terminal positioning solution at each epoch using the MHSS algorithm.

[0055] The beneficial effect of adopting the above technical solution is that: the present invention provides an RTK positioning vertical protection level optimization method based on an improved double-difference subset. Since the station-satellite double-difference observation can eliminate the influence of satellite clock error and receiver clock error in the GNSS system, it can reduce the ionospheric and tropospheric delay errors and improve the accuracy and reliability of the observation. By combining ARAIM with the multi-hypothesis analysis capability of MHSS, it can not only improve the sensitivity of fault detection, but also effectively reduce the influence of a single fault on the positioning results, and enhance the overall reliability of the receiver terminal positioning. BRIEF DESCRIPTION OF THE DRAWINGS

[0056] Figure 1 A flowchart of a method for optimizing the vertical protection level of RTK positioning based on an improved double difference subset provided in an embodiment of the present invention;

[0057] Figure 2 A schematic diagram of double difference observation provided by an embodiment of the present invention;

[0058] Figure 3 A flow chart of calculating the vertical protection level of the positioning solution of a receiver terminal at each epoch based on MHSS provided in an embodiment of the present invention. DETAILED DESCRIPTION

[0059] The specific implementation of the present invention is further described in detail below in conjunction with the accompanying drawings and examples. The following examples are used to illustrate the present invention, but are not intended to limit the scope of the present invention.

[0060] This embodiment provides an RTK positioning vertical protection level optimization method based on an improved double difference subset, such as Figure 1 As shown, the following steps are included:

[0061] Step 1: Adjust the existing RTK high-precision positioning model parameters; introduce error correction to reduce the positioning error caused by ionospheric delay and tropospheric delay factors in the long baseline case, and update the parameters of the state quantity to be estimated of the GNSS receiver terminal, including the following steps:

[0062] This embodiment first analyzes and improves the existing RTK high-precision positioning model, and further improves its accuracy and reliability by finely adjusting the parameters of the existing RTK high-precision positioning model; introduces error correction to reduce the positioning error caused by factors such as ionospheric delay and tropospheric delay in the case of long baseline;

[0063] In order to reduce the positioning error caused by ionospheric delay and tropospheric delay factors in the long baseline case, the ionospheric and tropospheric error correction modules are introduced; these error correction modules are processed by the double difference method, and in the long baseline case, the double difference ionospheric and double difference tropospheric errors are introduced into the error correction module as unknowns, and these errors are estimated using Kalman filtering;

[0064] Step 1.1: Establish GNSS non-difference observation equations to obtain the functional relationship between the original observations of the GNSS receiver terminal and the unknown parameters in the satellite navigation and positioning system;

[0065] GNSS non-differenced observation equations are a set of mathematical equations used to describe the relationship between satellite signals received by satellite signal receivers and factors such as the geometric relationship between satellites and receivers, clock errors, and atmospheric delays. These equations do not rely on differential technology, that is, they do not use the differences in observations between different receivers or different satellites to eliminate errors, but directly use the observations of a single satellite received by a single receiver for calculation.

[0066] GNSS undifferenced observation equations usually include carrier phase observation equations and pseudorange observation equations. Both observation equations take into account factors such as the geometric distance between the satellite and the receiver, clock error, and atmospheric delay (such as ionospheric delay and tropospheric delay).

[0067] The original observations of the GNSS receiver terminal include pseudorange observations and carrier phase observations; the carrier phase observation is the difference between the carrier phase of the satellite signal measured by the receiver and the carrier phase generated by the receiver at the same measurement time. The pseudorange observation is the distance value obtained by multiplying the time required for the satellite signal to propagate from the satellite to the receiver by the speed of light, but due to the existence of various errors, this distance value is not completely equal to the true geometric distance between the satellite and the receiver.

[0068] At a certain moment, the satellite s observed by the GNSS receiver terminal r is at a frequency f i Pseudorange observations on and carrier phase observation Respectively expressed as:

[0069]

[0070] in, is the geometric distance between satellite s and receiver terminal r, c is the light speed of satellite s, dt r and dt s are the clock errors of the GNSS receiver terminal r and the satellite s, is the ionospheric delay between satellite s and GNSS receiver terminal r, which is proportional to frequency f i related, is the tropospheric delay between satellite s and GNSS receiver terminal r, λ fi is the wavelength of satellite s, For GNSS receiver terminal r at frequency f i The carrier phase integer ambiguity on , ε and ξ are pseudorange measurement error and carrier phase measurement error respectively;

[0071] Step 1.2: Differencing the GNSS non-differenced observation equation in step 1.1, obtain the single-differenced observation equation between GNSS stations, thereby eliminating the satellite clock error at the satellite s end;

[0072] GNSS receiver terminal r and GNSS receiver terminal b observe the same satellite s at the same time, and then use the observations of satellite s by GNSS receiver terminal r and GNSS receiver terminal b as the single difference between GNSS stations, that is, make a difference to the GNSS non-difference observation equation in step 1.1, and obtain the single difference observation equation between GNSS stations:

[0073]

[0074] in, is the single difference operator between GNSS stations, which indicates the single difference calculation process between GNSS stations. is the satellite s observed by GNSS receiver terminal r and GNSS receiver terminal b at frequency f i The single-difference observations between pseudo-range stations on is the satellite s observed by GNSS receiver terminal r and GNSS receiver terminal b at frequency f i The carrier phase single-difference observation between stations on is the geometric distance difference between GNSS receiver terminal r and GNSS receiver terminal b and satellite s; dt rb is the clock difference between GNSS receiver terminal r and GNSS receiver terminal b; is the ionospheric delay difference between GNSS receiver terminal r and GNSS receiver terminal b and satellite s; is the tropospheric delay difference between GNSS receiver terminal r and GNSS receiver terminal b and satellite s; For GNSS receiver terminal r and GNSS receiver terminal b at frequency f i The inter-station single-difference carrier phase integer ambiguity on ;

[0075] In this embodiment, the pseudo-range observation noise error and carrier phase observation noise error after single difference between GNSS stations are still and It means that the observation noise after single difference between GNSS stations is amplified compared with the observation noise after GNSS non-difference observation. From the above formula, it can be seen that the single difference between GNSS stations can eliminate the error term at the satellite s end.

[0076] Step 1.3: Based on the single-difference observation equation between GNSS stations, GNSS receiver terminal r and GNSS receiver terminal b simultaneously observe satellite s and satellite k, as Figure 2As shown, a long baseline GNSS station-satellite double difference observation equation is established between satellite s and satellite k to eliminate the error of the GNSS receiver terminal r, and the state quantity parameters to be estimated of the GNSS receiver terminal r in the long baseline case are obtained;

[0077] The long baseline GNSS station double difference observation equation is:

[0078]

[0079] in, is the GNSS station-satellite double difference operator, which represents the GNSS station-satellite double difference calculation process. is the satellite s and satellite k observed by receiver terminal r and receiver terminal b at frequency f i The pseudo-range station-satellite double difference observations on; is the satellite s and satellite k observed by receiver terminal r and receiver terminal b at frequency f i The carrier phase station-satellite double difference observations on the is the geometric distance difference between receiver terminal r and receiver terminal b and satellite s and satellite k; is the ionospheric delay difference between receiver terminal r and receiver terminal b and satellite s and satellite k; is the tropospheric delay difference between receiver terminal r and receiver terminal b and satellite s and satellite k; For receiver terminal r and receiver terminal b at frequency f i The carrier phase station-satellite double difference integer ambiguity on the

[0080] From the above formula, we can see that the station-satellite double difference can eliminate the error of the GNSS receiver terminal r;

[0081] When the distance between the GNSS receiver terminal r and the GNSS receiver terminal b is less than the set threshold (20 km), that is, when the vehicle is located in an urban area, the tropospheric and ionospheric errors can also be ignored, so the long baseline GNSS station double difference observation equation is rewritten as:

[0082]

[0083] Since short baseline RTK can eliminate most of the errors between the satellite and the receiver terminal r, in the short baseline RTK positioning mode, the data of the L1 and L2 frequency bands in the global satellite navigation system are received by the GNSS receiver terminal r, and the position parameters of the GNSS receiver terminal r and the integer ambiguity parameters of the single-difference carrier phase between stations are estimated. The state quantity x to be estimated at the GNSS receiver terminal r is expressed as:

[0084]

[0085] Among them, r ris the position parameter of the GNSS receiver terminal r, v r is the velocity parameter of the GNSS receiver terminal r, a rb,1 and a rb,2 They are the single-difference carrier phase integer ambiguity parameters between stations in the L1 and L2 frequency bands respectively;

[0086] The L1 band is one of the most commonly used bands in the global satellite navigation system, mainly used by GPS (Global Positioning System) and GLONASS (Russian Global Navigation Satellite System); the center frequency of the L1 band is 1575.42MHz, with high signal strength and good penetration ability; the navigation signal of the L1 band is suitable for a wide range of applications, including aviation, navigation, vehicle navigation and consumer electronics such as smart phones; the L1 band is characterized by relatively low positioning accuracy and is greatly affected by factors such as multipath effect and atmospheric delay;

[0087] The L2 band is the second core band in the global satellite navigation system, mainly used to improve positioning accuracy and resist interference; the center frequency of the L2 band is 1227.60MHz, which is usually used for military and high-precision applications in GPS and GLONASS systems; the signal of the L2 band has high anti-interference performance and can effectively resist malicious interference and signal shielding attacks; in addition, the L2 band can also achieve higher positioning accuracy by measuring the carrier phase, which is used in scientific research, mapping and measurement.

[0088] When the vehicle is located in a non-urban area, the troposphere and ionosphere of the long baseline (>20km) RTK positioning model are difficult to eliminate. Therefore, the long baseline tropospheric wet delay and ionospheric delay are treated as unknowns for parameter estimation. The ionosphere and troposphere state quantities are added to the Kalman filter state quantity, and their time changes are simulated as a random walk model.

[0089] In this embodiment, in the long baseline RTK positioning mode, the state changes of the ionosphere and troposphere are simulated as random walk processes, and their covariance matrices are also introduced into the Kalman filter prediction; in the prediction and update steps of the Kalman filter, the covariance matrix of the error is propagated through the state transfer matrix, and the influence of the ionosphere and troposphere error covariance is fully considered to achieve more accurate error estimation and elimination;

[0090] The long baseline GNSS station double difference observation equation is further obtained as:

[0091]

[0092] Then, the estimated state parameter x1 of the GNSS receiver terminal r in the long baseline case is:

[0093]

[0094] Among them, Z T,r and Z T,b are the total delays in the tropospheric direction of GNSS receiver terminal r and GNSS receiver terminal b, respectively, is the inter-station double difference of the ionospheric delay error in the zenith direction on the L1 band, k is the reference satellite, and m is a different satellite system;

[0095] Step 1.4: Use the Kalman filter process to update the estimated state quantity parameters of the GNSS receiver terminal r;

[0096] The magnitude of the tropospheric error is related to the baseline length and the elevation difference between the GNSS receiver terminal r and the GNSS receiver terminal b. The constraint and random walk noise of the troposphere in the zenith direction in the L1 band are defined as:

[0097] P trop =(log(1+D·5·10 -4 )·0.05+H·5·10 -5 )

[0098]

[0099] Among them, P trop is the covariance matrix corresponding to the tropospheric state quantity, Q trop is the process noise matrix corresponding to the tropospheric state quantity, D is the baseline length, H is the elevation difference between receiver terminal r and receiver terminal b, and △t is the epoch interval;

[0100] The magnitude of the ionospheric error is not only related to the baseline length and the elevation difference between the GNSS receiver terminal r and the GNSS receiver terminal b, but is also closely related to the latitude position. Therefore, the ionospheric constraint and random walk noise are defined as:

[0101] P iono =(D·5·10 -6 exp((90-lat) / 50-1)) / sin(El)

[0102]

[0103] Among them, P iono is the covariance matrix corresponding to the ionospheric state quantity, Q iono is the process noise matrix corresponding to the ionospheric state quantity, lat is the average latitude of the receiver terminal b and the GNSS receiver terminal r, and El is the altitude angle of the observed satellite;

[0104] The correlation covariance matrix of the ionosphere and troposphere is introduced into the Kalman filter process. In the Kalman filter prediction step, the state estimation value of the GNSS receiver terminal r at the previous moment is used. and its covariance matrix P k-1|k-1 , predict the state estimate of the GNSS receiver terminal r at the current moment and its covariance matrix P k|k-1 :

[0105]

[0106] Among them, Q k is the predicted noise covariance matrix of GNSS receiver terminal r, F k is the state transfer matrix of the GNSS receiver terminal r;

[0107] The predicted noise covariance matrix Q of the GNSS receiver terminal r is k The covariance matrix Q including tropospheric error and ionospheric error trop and Q iono , and through the state transfer matrix F k To spread;

[0108] In the Kalman filter update step, the observation z is used k To update the parameters of the state quantity to be estimated; the Kalman filter process uses the observation matrix H k and the observation noise covariance matrix R k To make a status update:

[0109]

[0110] Among them, K k is the Kalman gain, H k T H k The transposed matrix of is the updated state quantity to be estimated;

[0111] Kalman gain K k The calculation process of fully considers the influence of the ionosphere and troposphere error covariance matrix;

[0112] The covariance matrix P of the GNSS receiver terminal r after the estimated state quantity is updated k As shown in the following formula:

[0113] P k =(IK k H k ) k|k-1

[0114] Where I is the identity matrix;

[0115] In this update equation, the covariance matrix P k As each observation update gradually decreases, it shows that the influence of ionospheric and tropospheric errors is gradually eliminated; in long baseline RTK positioning, the use of this RTK high-precision positioning model can effectively reduce the positioning error caused by ionospheric and tropospheric errors, making the final positioning solution more accurate and reliable;

[0116] This optimized positioning function model can improve the accuracy and stability under different baseline conditions. First, the double-difference observation equation performs differential operations between the receiver and the satellite, effectively eliminating most of the system errors on the satellite and receiver sides, such as satellite clock errors and receiver clock errors; at the same time, through the modeling and estimation of tropospheric and ionospheric errors, under long baseline conditions, even if these errors cannot be completely eliminated, they can be corrected in real time by adding state parameters and using Kalman filtering. This optimization model adaptively adjusts the random walk model of ionospheric and tropospheric delay errors, taking into account key factors such as baseline length, elevation difference and latitude position, making the model more robust and adaptable. In long-baseline RTK positioning, the use of this model can effectively reduce the positioning error caused by ionospheric and tropospheric errors, making the final positioning solution more accurate and reliable;

[0117] Step 2: Use GNSS receivers to obtain GNSS raw observation data of GPS, BDS, GLONASS and Galileo satellite systems in real time;

[0118] The GNSS raw observation data includes carrier phase, pseudorange and satellite coordinate information, which is the basis for subsequent processing and integrity monitoring. To ensure the integrity and reliability of the data, the GNSS raw observation data is screened, abnormal data is removed and the remaining data is evaluated for quality.

[0119] Step 3: Process the GNSS station-to-station single-difference observation equation and the long baseline GNSS station-to-satellite double-difference observation equation given in step 1 to obtain the GNSS station-to-satellite double-difference observation matrix; wherein each set of observations in the GNSS station-to-satellite double-difference observation matrix is ​​a subset, and the GNSS original observation data obtained in step 2 is used for testing and verification;

[0120] Process the single-difference observation equations between GNSS stations and the double-difference observation equations between long-baseline GNSS stations and satellites, and construct the double-difference observation matrix between GNSS stations and satellites;

[0121] The single-difference observation equation between GNSS stations is rearranged and expressed as: V1 = A1·δx1 using the prior residual matrix V1, the single-difference observation matrix A1 between GNSS stations, and the single-difference state quantity to be estimated between GNSS stations δx1.

[0122]

[0123] Among them, the elements in V1 They are different system frequencies f i The satellite pseudorange and carrier phase residuals are (·) for different satellite systems, G for GPS, C for BDS, They are different system frequencies f i The satellite direction cosines, and l, m, n are the components of the direction cosines on the x, y and z axes respectively. The dx, dy, dz contained in the estimated state δx1 are the components of the position correction in the x, y and z axes respectively. c·δt is the influence of the clock error. The influence of clock errors between different satellite systems and the GPS system. The frequency of the jth satellite is f i The fuzziness parameter of

[0124] Then select a reference satellite as the benchmark satellite, subtract the observation value of the benchmark satellite from the observation value of the same system and the same frequency, and get the new GNSS station-satellite double difference observation equation: V = A·δx, where V is the GNSS station-satellite double difference residual, A is the GNSS station-satellite double difference observation matrix, δx is the GNSS station-satellite double difference state quantity to be estimated,

[0125]

[0126] Among them, p (·)j is the pseudorange residual of the jth satellite, φ (·)j is the carrier phase residual of the jth satellite, N (·)j is the ambiguity parameter of the jth satellite, l (·)n ,m (·)n ,n (·)n are the direction cosines of satellites in different systems;

[0127] The GNSS station-satellite double-difference observation matrix A contains the geometric relationship information between the GNSS receiver and the satellite, preparing for subsequent integrity monitoring;

[0128] Each set of observations in the GNSS station-satellite double-difference observation matrix is ​​obtained by performing a differential operation between the reference satellite and the GNSS station-to-station single-difference observations of a certain observation satellite. By combining different satellites, multiple different subsets are formed, each of which corresponds to a different satellite configuration. By analyzing these subsets, the reliability and robustness of the positioning results can be effectively evaluated. When an abnormality occurs in a subset, the problem satellite can be identified and isolated by elimination or weighted adjustment to ensure the overall integrity of the positioning system. This method not only improves the system's anti-interference ability, but also enhances the accuracy and stability of positioning.

[0129] Step 4: Exclude the fault data of each subset and calculate the vertical protection level of the terminal positioning solution of the receiver at each epoch by using the MHSS algorithm, such as Figure 3 As shown;

[0130] For each subset, the MHSS algorithm considers different satellite failure assumptions and calculates the positioning solution and corresponding protection level (PL) of the receiver terminal r at each moment. The protection level is a key indicator used to quantify the maximum error range of the position solution under a specific confidence level, and to calculate and analyze the user's vertical protection level (VPL). The MHSS method uses solution separation technology to calculate the solution separation values ​​of each subset under different failures. These values ​​are used to evaluate the error performance of the GNSS system under various possible failure scenarios. The MHSS algorithm identifies potential abnormal satellites or faulty subsets based on the solution separation values. If the protection level of a subset exceeds the safety threshold, it means that the subset may be affected by errors or failures. The algorithm will eliminate or reduce the weight of the subset to ensure that the final positioning result meets the integrity requirements.

[0131] The linear expression of the station-satellite double difference observation equation is as follows:

[0132] L=AX+△

[0133] Where L is the difference between the observed value and the linearized predicted value, X is the m×1 dimensional state vector to be estimated, which consists of positioning parameters and receiver clock error parameters, A is the n×m dimensional design matrix, which is also the observation matrix in the MHSS algorithm, and △ is the n×1 dimensional observation error vector; when the observed value is normal, △ obeys the distribution E n (0, ∆ω), in is the variance of each observation, i.e. the measurement precision of different observations;

[0134] Based on the linearized observation equation, the following analysis is performed:

[0135] Where n is the total number of observations;

[0136] Assume H 1,i If it holds, then the i-th observation value contains a fault, that is, after removing the i-th observation value from the formula, the observation equation is obtained:

[0137] L(i)=A(i)X+△(i)

[0138] Among them, i is the serial number assigned when reading the observation value;

[0139] The satellite numbers of non-reference satellites are stored in correspondence with the corresponding observation numbers for subsequent elimination. Using the least squares method, the estimated value of X can be obtained by substituting it into the observation equation of the uneliminated observations. (0)

[0140] x (0) =S0L

[0141] Where S0=(A T W URA A) -1 A T W URA , W URA is a diagonal matrix with the jth diagonal element being is the variance of multipath and receiver terminal errors, is the variance of the tropospheric residual;

[0142] Correspondingly, the estimated value of X x can be obtained by removing the new observation equation after the i-th observation value is removed. (i)

[0143] x (i) =S (i) L

[0144] Among them, S (i) =(A T M i W URA A) -1 A T M i W URA , M i is the identity matrix whose (i,i)th element is 0;

[0145] Combine x (0) and x (i) , we get the test statistic d to determine whether the i-th double difference observation value is faulty. (i) , the specific formula is:

[0146] d (i) =|x z(i) -x z(0) |

[0147] Among them, x z(i) is the parameter to be estimated in the vertical direction after removing the i-th double difference observation, x z(0) is the parameter to be estimated in the vertical direction estimated using all observations;

[0148] And find the detection threshold D (i) for:

[0149]

[0150] Among them, K ffd,j The false alarm rate P FA The upper 1-P value of the standard normal distribution obtained for the significance level FA / 2 Quantile value; z is the vertical direction, and △S (i) =S (i) -S0,b norm is the continuity deviation value;

[0151] By comparing the test statistic d (i) With the detection threshold D (i) , to determine whether the i-th observation value has a fault; when the detection statistic d (i) Exceeding the detection threshold D (i) When , it is considered that the observed value has a fault; otherwise, it is considered that there is no fault; for the observed value with a fault, the observed value is eliminated and re-tested until all the observed values ​​meet the fault detection conditions, thus completing the troubleshooting of the observed data;

[0152] The calculation of the user vertical protection level is divided into two cases: no fault and fault. When there is no fault, the user protection level is:

[0153]

[0154] Among them, VPL0 is the user protection level when there is no fault, K md,0 The error rate P md,0 The upper 1-P value of the standard normal distribution obtained for the significance level md,0 / 2 Quantile value, P0=(A T W URA A) -1 , b max is the integrity deviation value;

[0155] If the fault cannot be eliminated and the i-th satellite still contains potential faults, the user protection level is:

[0156]

[0157] Among them, VPL (i) is the user protection level when the i-th satellite still contains potential failures, K md,i The error rate P md,i The upper 1-P value of the standard normal distribution obtained for the significance level md,i / 2 Quantile value; P (i) =(A T M i W URA A) -1 ;

[0158] The final user protection level VPL depends on the worst case under all fault assumptions:

[0159] VPL=max{VPL0,max VPL (i)}

[0160] Through this process, the multi-hypothesis analysis capability of MHSS can be utilized to not only improve the sensitivity of fault detection, but also effectively reduce the impact of a single fault on the positioning results, thereby enhancing the overall reliability of the system.

[0161] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, rather than to limit it. Although the present invention has been described in detail with reference to the aforementioned embodiments, those skilled in the art should understand that they can still modify the technical solutions described in the aforementioned embodiments, or make equivalent replacements for some or all of the technical features therein. However, these modifications or replacements do not cause the essence of the corresponding technical solutions to deviate from the scope defined by the claims of the present invention.

Claims

1. A method for optimizing the vertical protection level of RTK positioning based on an improved double difference subset, characterized in that: The following steps are involved: Step 1: Adjust the existing RTK high-precision positioning model parameters, introduce error correction to reduce the positioning error caused by ionospheric delay and tropospheric delay in the long baseline case, and update the parameters of the state quantity to be estimated of the GNSS receiver terminal; Step 2: Use the GNSS receiver to obtain the GNSS raw observation data of different satellite navigation systems in real time, and screen the GNSS raw observation data, remove abnormal data and conduct quality assessment on the remaining data; Step 3: Process the GNSS station-to-station single-difference observation equation and the long baseline GNSS station-to-satellite double-difference observation equation given in step 1 to obtain the GNSS station-to-satellite double-difference observation matrix; wherein each set of observations in the GNSS station-to-satellite double-difference observation matrix is ​​a subset, and the GNSS original observation data obtained in step 2 is used for testing and verification; Step 4: Exclude the fault data of each subset and calculate the vertical protection level of the receiver terminal positioning solution at each epoch using the MHSS algorithm.

2. The method for optimizing the vertical protection level of RTK positioning based on an improved double difference subset according to claim 1, characterized in that: The step 1 comprises: Step 1.1: Establish GNSS non-difference observation equations to obtain the functional relationship between the original observations of the GNSS receiver terminal and the unknown parameters in the satellite navigation and positioning system; Step 1.2: Differencing the GNSS non-differenced observation equation in step 1.1 to obtain the single-differenced observation equation between GNSS stations, thereby eliminating the satellite clock error at the satellite end; Step 1.3: Based on the single-difference observation equation between GNSS stations, GNSS receiver terminal r and GNSS receiver terminal b simultaneously observe satellite s and satellite k, and establish a long-baseline GNSS station-satellite double-difference observation equation between satellite s and satellite k to eliminate the error of GNSS receiver terminal r, and obtain the state quantity parameters to be estimated of GNSS receiver terminal r in the long-baseline case; Step 1.4: Use the Kalman filter process to update the estimated state parameters of the GNSS receiver terminal r.

3. The method for optimizing the vertical protection level of RTK positioning based on the improved double difference subset according to claim 2, characterized in that: The original observations of the GNSS receiver terminal include pseudorange observations and carrier phase observations; At a certain moment, the satellite s observed by the GNSS receiver terminal r is at a frequency f i Pseudorange observations on and carrier phase observation Respectively expressed as: in, is the geometric distance between satellite s and receiver terminal r, c is the light speed of satellite s, dt r and dt s are the clock errors of the GNSS receiver terminal r and the satellite s, is the ionospheric delay between satellite s and GNSS receiver terminal r, which is proportional to frequency f i related, is the tropospheric delay between the satellite s and the GNSS receiver terminal r, is the wavelength of satellite s, For GNSS receiver terminal r at frequency f i The carrier phase integer ambiguity on , ε and ξ are the pseudorange measurement error and carrier phase measurement error respectively.

4. The method for optimizing the vertical protection level of RTK positioning based on the improved double difference subset according to claim 3 is characterized in that: The step 1.2 comprises: GNSS receiver terminal r and GNSS receiver terminal b observe the same satellite s at the same time, and then use the observations of satellite s by GNSS receiver terminal r and GNSS receiver terminal b as the single difference between GNSS stations, that is, make a difference to the GNSS non-difference observation equation in step 1.1, and obtain the single difference observation equation between GNSS stations: in, is the single difference operator between GNSS stations, is the satellite s observed by GNSS receiver terminal r and GNSS receiver terminal b at frequency f i The single-difference observations between pseudo-range stations on is the satellite s observed by GNSS receiver terminal r and GNSS receiver terminal b at frequency f i The carrier phase single-difference observation between stations on is the geometric distance difference between GNSS receiver terminal r and GNSS receiver terminal b and satellite s; dt rb is the clock difference between GNSS receiver terminal r and GNSS receiver terminal b; is the ionospheric delay difference between GNSS receiver terminal r and GNSS receiver terminal b and satellite s; is the tropospheric delay difference between GNSS receiver terminal r and GNSS receiver terminal b and satellite s; For GNSS receiver terminal r and GNSS receiver terminal b at frequency f i The inter-station single-difference carrier phase integer ambiguity on .

5. The method for optimizing the vertical protection level of RTK positioning based on improved double difference subset according to claim 4 is characterized in that: The long baseline GNSS station double difference observation equation in step 1.3 is: in, is the GNSS station double difference operator, is the satellite s and satellite k observed by receiver terminal r and receiver terminal b at frequency f i The pseudo-range station-satellite double difference observations on; is the satellite s and satellite k observed by receiver terminal r and receiver terminal b at frequency f i The carrier phase station-satellite double difference observations on the is the geometric distance difference between receiver terminal r and receiver terminal b and satellite s and satellite k; is the ionospheric delay difference between receiver terminal r and receiver terminal b and satellite s and satellite k; is the tropospheric delay difference between receiver terminal r and receiver terminal b and satellite s and satellite k; For receiver terminal r and receiver terminal b at frequency f i The carrier phase station-satellite double difference integer ambiguity on the When the distance between GNSS receiver terminal r and GNSS receiver terminal b is less than the set threshold, the tropospheric and ionospheric errors are ignored, so the short baseline GNSS station double difference observation equation is rewritten as: In the short baseline RTK positioning mode, the data of L1 and L2 frequency bands in the global satellite navigation system are received by the GNSS receiver terminal r, and the position parameters of the GNSS receiver terminal r and the integer ambiguity parameters of the single-difference carrier phase between stations are estimated. Then the state quantity x to be estimated of the GNSS receiver terminal r is expressed as: Among them, r r is the position parameter of the GNSS receiver terminal r, v r is the velocity parameter of the GNSS receiver terminal r, a rb,1 and a rb,2 They are the single-difference carrier phase integer ambiguity parameters between stations in the L1 and L2 frequency bands respectively; When the distance between GNSS receiver terminal r and GNSS receiver terminal b is greater than or equal to the set threshold, the tropospheric wet delay and ionospheric delay of the long baseline are taken as unknowns for parameter estimation, and the ionosphere and troposphere state quantities are added to the Kalman filter state quantity, and the time changes of the two are simulated as a random walk model; The long baseline GNSS station double difference observation equation is further obtained as: Then, the estimated state parameter x1 of the GNSS receiver terminal r in the long baseline case is: Among them, Z T,r and Z T,b are the total delays in the tropospheric direction of GNSS receiver terminal r and GNSS receiver terminal b, respectively, is the inter-station double difference of the ionospheric delay error in the zenith direction on the L1 band, k is the reference satellite, and m is a different satellite system.

6. The method for optimizing the vertical protection level of RTK positioning based on improved double difference subsets according to claim 5, characterized in that: The step 1.4 comprises: The covariance matrix corresponding to the state quantities of the ionosphere and troposphere is introduced into the Kalman filter process. In the Kalman filter prediction step, the state estimation value of the GNSS receiver terminal r at the previous moment is used. and its covariance matrix P k-1|k-1 , predict the state estimate of the GNSS receiver terminal r at the current moment and its covariance matrix P k|k-1 : Among them, Q k is the predicted noise covariance matrix of GNSS receiver terminal r, F k is the state transfer matrix of the GNSS receiver terminal r; The predicted noise covariance matrix Q of the GNSS receiver terminal r is k The covariance matrix of the tropospheric error and the ionospheric error is included, and the state transfer matrix F k To spread; In the Kalman filter update step, the observation z is used k To update the parameters of the state quantity to be estimated; the Kalman filter process uses the observation matrix H k and the observation noise covariance matrix R k Update the state of the state to be estimated: Among them, K k is the Kalman gain, H k T H k The transposed matrix of is the updated state quantity to be estimated; The covariance matrix P of the GNSS receiver terminal r after the estimated state quantity is updated k As shown in the following formula: P k =(I-K k H k )P k|k-1 Where I is the identity matrix.

7. The method for optimizing the vertical protection level of RTK positioning based on improved double difference subset according to claim 1, characterized in that: The different satellite navigation systems described in step 2 include four satellite navigation systems: GPS, BDS, GLONASS and Galileo; the GNSS raw observation data includes carrier phase, pseudorange and satellite coordinate information.

8. The method for optimizing the vertical protection level of RTK positioning based on improved double difference subsets according to claim 6, characterized in that: The step 3 comprises: Process the single-difference observation equations between GNSS stations and the double-difference observation equations between long-baseline GNSS stations and satellites, and construct the double-difference observation matrix between GNSS stations and satellites; The single-difference observation equation between GNSS stations is rearranged and expressed as: V1 = A1·δx1 using the prior residual matrix V1, the single-difference observation matrix A1 between GNSS stations, and the single-difference state quantity to be estimated between GNSS stations δx1. Among them, the elements in V1 They are different system frequencies f i The satellite pseudorange and carrier phase residuals are (·) for different satellite systems, G for GPS, C for BDS, The frequencies of different satellite navigation systems are f i The satellite direction cosines, and l, m, n are the components of the direction cosines on the x, y and z axes respectively. The dx, dy, dz contained in the estimated state δx1 are the components of the position correction in the x, y and z axes respectively. c·δt is the influence of the clock error. The influence of clock errors between different satellite systems and the GPS system. The frequency of the jth satellite is f i The fuzziness parameter of Then select a reference satellite as the benchmark satellite, subtract the observation value of the benchmark satellite from the observation value of the same frequency of the same satellite navigation system, and obtain the new GNSS station-satellite double difference observation equation: V = A·δx, where V is the GNSS station-satellite double difference residual, A is the GNSS station-satellite double difference observation matrix, δx is the GNSS station-satellite double difference state quantity to be estimated, Among them, p (·)j is the pseudorange residual of the jth satellite, φ (·)j is the carrier phase residual of the jth satellite, N (·)j is the ambiguity parameter of the jth satellite, l (·)n ,m (·)n ,n (·)n are the direction cosines of satellites in different systems; The GNSS station-satellite double-difference observation matrix A contains the geometric relationship information between the GNSS receiver and the satellite. Each set of observations in the GNSS station-satellite double-difference observation matrix is ​​obtained by performing a differential operation between the GNSS station-to-station single-difference observations of a reference satellite and an observation satellite. By combining different satellites, multiple different subsets are formed, and each subset corresponds to a different satellite configuration.

Images