A GNSS multi-system robust adaptive fusion RTK solution method

Abstract

The invention discloses a GNSS multi-system robust adaptive fusion RTK solution method, comprising the following steps: step 1: establishing a conventional RTK double-difference observation equation according to the BDS / GPS / GLONASS / GLILEO multi-system observation data of the base station and the mobile station; step 2: conventional Kalman filter solution; step 3: robust Kalman adaptive filter estimation; step 4: BDS / GPS / GLONASS / GLILEO multi-system fusion Helmert variance component estimation; step 5: BDS / GPS / GLONASS / GLILEO multi-system RTK ambiguity fixation. The invention can effectively reduce the influence of the inaccurate prior random model of GNSS positioning, the abnormal observation and the abnormal dynamic model on the multi-system fusion RTK positioning, improve the positioning accuracy of the GNSS (BDS / GPS / GLONASS / GLILEO) multi-system fusion RTK, and make the RTK positioning result more stable and reliable in complex environments.

Description

Technical Field

[0001] The invention belongs to the field of GNSS satellite navigation and positioning, and in particular relates to a GNSS multi-system robust adaptive fusion RTK solution method. Background Art

[0002] With the continuous development of the global satellite navigation system, the combination of different satellite constellations has gradually increased the number of visible satellites, and the spatial geometric distribution of satellites has also become more perfect. GNSS (BDS / GPS / GLONASS / GLILEO) multi-system fusion has introduced more visible satellites, providing a more reliable guarantee for high-precision positioning. However, due to the differences in the distribution of satellite constellations and satellite orbits of various satellite navigation systems, the navigation positioning performance and observation value accuracy of each system are different. The GNSS multi-system fusion positioning is not accurate and reliable if it only determines the weights of the observation values ​​of each system based on the prior accuracy. At the same time, in cities, canyons and other areas, due to the complex observation environment, satellite signals form multipaths or non-line-of-sight signals due to reflection or scattering of buildings, resulting in an increased probability that the observation values ​​contain gross errors or anomalies, and the carrier motion model under dynamic conditions is difficult to accurately model, making the conventional Kalman filter positioning solution prone to large deviations or even divergence.

[0003] Therefore, for GNSS multi-system fusion positioning, especially in applications in complex environments, how to no longer rely on a priori theories or empirical values ​​to determine the weights of observation values ​​between navigation systems and satellites, and adaptively adjust the weights of observation values ​​between navigation systems and satellites, is of great significance to improving the accuracy and availability of GNSS multi-system fusion RTK positioning. Summary of the invention

[0004] The purpose of the present invention is to provide a GNSS multi-system robust adaptive fusion RTK solution method, which is used to effectively reduce the influence of inaccurate GNSS positioning prior random model, observation anomalies and dynamic model anomalies on multi-system fusion RTK positioning, improve the positioning accuracy of GNSS (BDS / GPS / GLONASS / GLILEO) multi-system fusion RTK, and make the RTK positioning result more stable and reliable in complex environments.

[0005] To achieve the above object, the technical solution of the present invention is:

[0006] A GNSS multi-system robust adaptive fusion RTK solution method, comprising the following steps:

[0007] Step 1: Establish the conventional RTK double-difference observation equation based on the BDS / GPS / GLONASS / GLILEO multi-system observation data of the base station and the rover;

[0008] Step 2: Conventional Kalman filter solution;

[0009] Step 3: Robust Kalman adaptive filter estimation;

[0010] Step 4: Helmert variance component estimation of BDS / GPS / GLONASS / GLILEO multi-system fusion;

[0011] Step 5: BDS / GPS / GLONASS / GLILEO multi-system RTK ambiguity is fixed.

[0012] Further, as a preferred technical solution, the specific process of step 1 is:

[0013] Step 1-1: Establish the basic observation equation based on the undifferenced raw pseudorange and carrier phase observations, as follows:

[0014]

[0015]

[0016] In the formula, the subscripts i and r are the frequency and station of the observed quantity, respectively, and the superscript j is the satellite; are the pseudorange and carrier phase raw observation values, in meters, respectively. is the phase observation value (week); is the geometric distance between the satellite and the station at the time of signal transmission, c is the speed of light in vacuum, dt r and dt j are the receiver clock error and satellite clock error, T j is the tropospheric delay, is the ionospheric delay at frequency i, and are the pseudorange hardware delay bias at the receiver and satellite, and are the phase hardware delay deviations of the receiver and satellite, including the carrier phase hardware delay deviation and the initial phase deviation, λ i and is the carrier phase wavelength and integer ambiguity of frequency i, and are the measurement noise of pseudorange and carrier phase observations, respectively;

[0017] Step 1-2: According to the RTK double-difference model, the satellite clock error and receiver clock error in step 1-1 can be eliminated, and the influence of satellite orbit error, tropospheric delay error, ionospheric delay error and multipath effect can be weakened, and the carrier phase double-difference observation equation can be established, as follows:

[0018]

[0019] In the formula, superscript j and k represent satellite numbers; subscript r and b represent mobile stations and base stations; is the double difference phase observation value at epoch t; is the distance from satellite s to station r; is the double difference ambiguity; since the reference coordinates are known, the above formula can be transformed into:

[0020]

[0021] Step 1-3: The prior coordinates of the mobile station are estimated using single point positioning, and the approximate estimate X r =(x r ,y r ,z r ) T , let the correction number be δX r =(δx r ,δy r ,δz r ) T , linearizing the final transformed equation in step 1-2 yields the carrier phase double difference error equation:

[0022]

[0023] In the formula,

[0024] Similarly, the pseudorange double difference error equation can be obtained:

[0025]

[0026] in, is the double-difference pseudorange observation value.

[0027] Further, as a preferred technical solution, the specific process of step 2 is:

[0028] Step 2-1: Establish the state equation and observation equation of Kalman filter:

[0029] X k =Φ k,k-1 X k-1 +W k-1

[0030] L k =A k X k +V k

[0031] Where k represents the observation epoch time, X k represents the state vector at time k, Φ k,k-1 represents the transition matrix of the system state from time k-1 to time k, Wk-1 is the dynamic noise vector, L k is the observation vector at time k, A k is the observation equation coefficient matrix, V k is the observation noise vector;

[0032] Step 2-2: The observations are considered to be independent of each other. The dynamic noise and observation noise are zero-mean, uncorrelated white noises, which satisfy the following conditions:

[0033]

[0034] In the formula, Ω k is the variance matrix of process noise, Q k is the variance matrix of the observation noise;

[0035] Step 2-3: Rewrite the state equation and observation equation into error equations respectively:

[0036]

[0037] In the formula, and V k are the predicted values ​​of state parameters and the observed value L k The residual vector of the state parameter is for:

[0038] The covariance matrix of can be obtained:

[0039]

[0040] Construct the objective function according to the least squares criterion:

[0041]

[0042] Where P k and are the observed values ​​L k and state parameter forecast values The weight matrix of

[0043]

[0044]

[0045] Paired Partial derivative and set its partial derivative equal to 0, we have:

[0046]

[0047] Further we can get:

[0048]

[0049] Step 2-4: Substitute the error equation corresponding to the observation equation into formula (1) to obtain:

[0050]

[0051] Furthermore, as a preferred technical solution, the general process of Kalman filter calculation is further derived from the matrix identity transformation and matrix inversion formula:

[0052] Step 2-4-1: Calculate the Kalman filter gain matrix K based on the prediction accuracy of the previous epoch for the current epoch k ,

[0053]

[0054] Step 2-4-2: Update the current epoch observations, i.e. use the weights of the current epoch measurements for weighted correction.

[0055]

[0056]

[0057] Step 2-4-3: Update time, that is, predict the observation value of the next epoch.

[0058] Further, as a preferred technical solution, the specific process of step 3 is:

[0059] Step 3-1: The robust adaptive Kalman filter solution is:

[0060]

[0061] In the formula, For L k The equivalent weight matrix of , which is the adaptive estimate of the observation vector weight matrix, α k is the adaptive factor, 0≤α k ≤1;

[0062] Step 3-2: Determine equivalent weights and adaptive weights.

[0063] Furthermore, as a preferred technical solution, the specific calculation process of the equivalent weight factor is:

[0064] Step 3-2-1: In BDS / GPS / GLONASS / GALILEO fusion positioning, the observation values ​​between systems are independent of each other, and the corresponding equivalent weights are:

[0065]

[0066] Where P i is the original weight matrix, w i is the adaptive weight factor;

[0067] Step 3-2-2: Determine the weight factor using an IGGIII equivalent weight factor or a dual-factor equivalent weight factor function, wherein the function expression based on the IGGIII equivalent weight factor is:

[0068]

[0069] In the formula, is the standardized residual, k1 and k2 are constants, generally k1 = 1.5 ~ 2.0, k2 = 3.0 ~ 8.0;

[0070] The function expression based on the double-factor equivalent weight factor is:

[0071]

[0072]

[0073] In the formula, c is a constant, c = 2.5 ~ 3.0, i and j represent the row and column numbers of the matrix;

[0074] Step 3-2-3: Calculate standardized residuals

[0075]

[0076] Step 3-2-4: For different types of observations, their standardized residuals for:

[0077]

[0078] In the formula, σ MAD =1.483·med|v i -med|v i ||.

[0079] Further, as a preferred technical solution, the specific calculation process of the adaptive weight factor is:

[0080] Adopt two-segment function to calculate adaptive weight factor a k ,Right now:

[0081]

[0082] In the formula, c is a constant, ranging from 1.0 to 2.5.

[0083]

[0084]

[0085] In the formula, is the standardized prediction residual, is the state parameter prediction residual, is the trace of the state parameter prediction residual covariance matrix.

[0086] Further, as a preferred technical solution, the specific process of step 4 is:

[0087] Step 4-1: During the first filtering, the BDS, GPS, GLONASS and GLILEO observations are weighted according to the satellite elevation angle. C :P G :P R :P E ;

[0088] Step 4-2: Perform filtering and solve to obtain V i T P i V i (i=C,G,R,E);

[0089] Step 4-3: Estimation of variance components:

[0090]

[0091] In the formula,

[0092]

[0093] N=A T PA=A C T P C A C +A G T P G A G +A R T P R A R +A E T P E A E =N C +N G +N R +N E

[0094] Among them, n C 、n G 、n R and n Eare the numbers of BDS, GPS, GLONASS and GALILEO observations, respectively;

[0095] Step 4-4: Re-evaluation:

[0096]

[0097] In the formula, c is any constant, optional A value in ;

[0098] Step 4-5: Repeat steps 4-2 to 4-4 until The difference is less than 0.01.

[0099] Furthermore, as a preferred technical solution, when the iteration does not converge or converges distortedly in step 4, the satellites are re-screened or the weights of certain satellites are appropriately reduced before solving the problem.

[0100] Furthermore, as a preferred technical solution, in step 5, the BDS / GPS / GLONASS / GLILEO multi-system RTK ambiguity is fixed by using a least squares ambiguity reduction correlation adjustment method combined with a partial ambiguity fixing strategy.

[0101] Compared with the prior art, the present invention has the following beneficial effects:

[0102] The present invention aims at GNSS multi-system fusion RTK. Firstly, the weight of each satellite observation value is adjusted by adopting the robust Kalman adaptive filtering to suppress the influence of observation anomaly and dynamic model error on RTK valuation parameters. Secondly, the weight of observation values ​​between various satellite navigation systems (BDS / GPS / GLONASS / GLILEO) is adaptively adjusted by using variance component estimation. The weight ratio between observation values ​​of different types or unequal precision is adaptively adjusted through iterative calculation, so that the distribution of weights is more reasonable, and the GNSS multi-system fusion RTK positioning is more accurate and reliable in complex environments. BRIEF DESCRIPTION OF THE DRAWINGS

[0103] Figure 1 It is a schematic diagram of the process structure of the present invention. DETAILED DESCRIPTION

[0104] The following examples are used to illustrate the present invention but are not intended to limit the scope of the present invention.

[0105] Example

[0106] like Figure 1 As shown, a GNSS multi-system robust adaptive fusion RTK solution method described in this embodiment includes the following steps:

[0107] Step 1: Establish the conventional RTK double-difference observation equation based on the BDS / GPS / GLONASS / GLILEO multi-system observation data of the base station and the rover;

[0108] Step 2: Conventional Kalman filter solution;

[0109] Step 3: Robust Kalman adaptive filter estimation;

[0110] Step 4: Helmert variance component estimation of BDS / GPS / GLONASS / GLILEO multi-system fusion;

[0111] Step 5: BDS / GPS / GLONASS / GLILEO multi-system RTK ambiguity is fixed.

[0112] In this embodiment, the specific process of step 1 is:

[0113] Step 1-1: Establish the basic observation equation based on the undifferenced raw pseudorange and carrier phase observations, as follows:

[0114]

[0115]

[0116] In the formula, the subscripts i and r are the frequency and station of the observed quantity, respectively, and the superscript j is the satellite; are the pseudorange and carrier phase raw observation values, in meters, respectively. is the phase observation value (week); is the geometric distance between the satellite and the station at the time of signal transmission, c is the speed of light in vacuum, dt r and dt j are the receiver clock error and satellite clock error, T j is the tropospheric delay, is the ionospheric delay at frequency i, and are the pseudorange hardware delay bias at the receiver and satellite, and are the phase hardware delay deviations of the receiver and satellite, including the carrier phase hardware delay deviation and the initial phase deviation, λ i and is the carrier phase wavelength and integer ambiguity of frequency i, and are the measurement noise of pseudorange and carrier phase observations, respectively;

[0117] Step 1-2: According to the RTK double-difference model, the satellite clock error and receiver clock error in step 1-1 can be eliminated, and the influence of satellite orbit error, tropospheric delay error, ionospheric delay error and multipath effect can be weakened, and the carrier phase double-difference observation equation can be established, as follows:

[0118]

[0119] In the formula, superscript j and k represent satellite numbers; subscript r and b represent mobile stations and base stations; is the double difference phase observation value at epoch t; is the distance from satellite s to station r; is the double difference ambiguity; since the reference coordinates are known, the above formula can be transformed into:

[0120]

[0121] Step 1-3: The prior coordinates of the mobile station are estimated using single point positioning, and the approximate estimate X r =(x r ,y r ,z r ) T , let the correction number be δX r =(δx r ,δy r ,δz r ) T , linearizing the final transformed equation in step 1-2 yields the carrier phase double difference error equation:

[0122]

[0123] In the formula,

[0124] Similarly, the pseudorange double difference error equation can be obtained:

[0125]

[0126] in, is the double-difference pseudorange observation value.

[0127] In this embodiment, the specific process of conventional Kalman filtering solution is:

[0128] Kalman filtering is based on a set of observation sequences L k (k=1,2,...,n) and the system's dynamic model information to solve the state vector estimate. In precision positioning, the state vector consists of the parameters to be estimated. The dynamic model is established based on the correlation between the parameters before and after the epochs to be estimated. The filtered state equation and observation equation are:

[0129] Xk =Φ k,k-1 X k-1 +W k-1

[0130] L k =A k X k +V k

[0131] Where k represents the observation epoch time, X k represents the state vector at time k, Φ k,k-1 represents the transition matrix of the system state from time k-1 to time k, W k-1 is the dynamic noise vector, L k is the observation vector at time k, A k is the observation equation coefficient matrix, V k is the observation noise vector;

[0132] Generally, the observations are considered to be independent of each other, and the dynamic noise and observation noise are zero-mean, uncorrelated white noises, that is, they satisfy the following conditions:

[0133]

[0134] In the formula, Ω k is the variance matrix of process noise, Q k is the variance matrix of the observation noise;

[0135] In order to solve the least squares solution of the state parameters, the Kalman filter is solved according to the free extreme value principle, and the state equation and observation equation are rewritten as error equations respectively:

[0136]

[0137] In the formula, and V k are the predicted values ​​of state parameters and the observed value L k The residual vector of the state parameter is for:

[0138] The covariance matrix of can be obtained:

[0139]

[0140] Construct the objective function according to the least squares criterion:

[0141]

[0142] Where P k and are the observed values ​​L k and state parameter forecast values The weight matrix of

[0143]

[0144]

[0145] Paired Partial derivative and set its partial derivative equal to 0, we have:

[0146]

[0147] Further we can get:

[0148]

[0149] Step 2-4: Substitute the error equation into the above formula to get:

[0150]

[0151] By using the matrix identity transformation and matrix inversion formula, the general process of Kalman filter calculation is further derived as follows:

[0152] Step 2-4-1: Calculate the Kalman filter gain matrix K based on the prediction accuracy of the previous epoch for the current epoch k ,

[0153]

[0154] Step 2-4-2: Update the current epoch observations, i.e. use the weights of the current epoch measurements for weighted correction.

[0155]

[0156]

[0157] Step 2-4-3: Update time, that is, predict the observation value of the next epoch.

[0158] In this embodiment, the specific process of robust Kalman adaptive filtering estimation is as follows:

[0159] In order to control the impact of abnormal dynamic models and observations, the robust estimation principle is adopted. By selecting appropriate equivalent weights to replace the observation noise covariance matrix, the influence of abnormal observation values ​​is reduced or eliminated. The ratio of the covariance matrix of the predicted value and the covariance matrix of the observation noise is adjusted by the adaptive factor to control the impact of abnormal dynamic models on the filter solution. The robust adaptive Kalman filter solution is:

[0160]

[0161] In the formula, For L k The equivalent weight matrix of , which is the adaptive estimate of the observation vector weight matrix, α k is the adaptive factor, 0≤α k ≤1;

[0162] Then determine the equivalent weight factor and the adaptive weight factor. The specific calculation process of the equivalent weight factor is:

[0163] In the BDS / GPS / GLONASS / GALILEO fusion positioning, the observation values ​​between systems are independent of each other, and the corresponding equivalent weights are:

[0164]

[0165] Where P i is the original weight matrix, w i is the adaptive weight factor;

[0166] The weight factor is determined by using an IGGIII equivalent weight factor or a dual-factor equivalent weight factor function, wherein the function expression based on the IGGIII equivalent weight factor is:

[0167]

[0168] In the formula, is the standardized residual, k1 and k2 are constants, generally k1 = 1.5 ~ 2.0, k2 = 3.0 ~ 8.0;

[0169] The function expression based on the double-factor equivalent weight factor is:

[0170]

[0171]

[0172] In the formula, c is a constant, c = 2.5 ~ 3.0, i and j represent the row and column numbers of the matrix;

[0173] Calculate standardized residuals

[0174]

[0175] For GNSS multi-system fusion positioning, due to the differences in the accuracy of observations from different systems, the residual vector of the observations may also have systematic errors, and the variance factor σ v The standard residuals of different types of observations are not consistent, which makes it impossible to accurately distinguish the gross errors from the correct residual information of different types of observations when the adaptive factor is downgraded. for:

[0176]

[0177] In the formula, σ MAD =1.483·med|v i -med|v i ||.

[0178] The specific calculation process of the adaptive weight factor is:

[0179] The dynamic model predicts the new information by using the state estimate of the previous moment and the observation value of the current moment to update the state new information and the weighted gain matrix of the covariance matrix at the current moment. The adaptive factor obtained by the statistics constructed based on the new information residual can better reflect the disturbance of the dynamic system. At the same time, in order to make full use of the dynamic model information, avoid a k When is equal to 0, the adaptive weight factor function adopts a two-segment function, namely:

[0180]

[0181] In the formula, c is a constant, ranging from 1.0 to 2.5.

[0182]

[0183]

[0184] In the formula, is the standardized prediction residual, is the state parameter prediction residual, is the trace of the state parameter prediction residual covariance matrix.

[0185] In this embodiment, the Helmert variance component estimation is to adaptively adjust the weight ratios between observations of different types or unequal precision through iterative calculations, so that the distribution of weights is more reasonable. The specific idea is to first determine the weights of various types of observations based on empirical values, perform adjustment, and estimate the a priori deviations of various types of observations based on the posterior residuals of the observations after adjustment, obtain the weights between the observations, and adjust again until the given test conditions are met. The Helmert variance component estimation is performed based on the established BDS / GPS / GLONASS / GLILEO multi-system double difference error equations and the a priori weight ratios of the observations of each system. The specific calculation process is:

[0186] Step 4-1: During the first filtering, the BDS, GPS, GLONASS and GLILEO observations are weighted according to the satellite elevation angle. C :P G :P R :P E ;

[0187] Step 4-2: Perform filtering and solve to obtain V i T P i V i (i=C,G,R,E);

[0188] Step 4-3: Estimation of variance components:

[0189]

[0190] In the formula,

[0191]

[0192] N=A T PA=A C T P C A C +A G T P G A G +A R T P R A R +A E T P E A E =N C +N G +N R +N E

[0193] Among them, n C 、n G 、n R and n E are the numbers of BDS, GPS, GLONASS and GALILEO observations, respectively;

[0194] Step 4-4: Re-evaluation:

[0195]

[0196] In the formula, c is any constant, optional A value in ;

[0197] Step 4-5: Repeat steps 4-2 to 4-4 until The difference is less than 0.01.

[0198] When the iteration does not converge or converges distortedly in step 4, the reason for analysis is that the number of observations of each type is unbalanced. If the number of observations of one type is too small, the weight ratio determined by the Helmert variance component estimation is easily distorted; second, in the process of observation data preprocessing, if the residuals of various observations differ too much, the iteration process may not converge or the number of iterations is too many, and the S matrix will be a singular matrix and unable to locate the solution. Therefore, during the solution process, it is necessary to pay attention to the number of satellites of each type in the epoch to decide whether to use the Helmert variance component estimation; when the iteration process diverges, it is necessary to analyze the residuals of each satellite, re-screen the satellites or appropriately reduce the weights of some satellites, and then solve.

[0199] In the present embodiment, BDS / GPS / GLONASS / GLILEO multi-system RTK ambiguity is fixed and preferably adopts LAMBDA (Least-square Ambiguity Decorrelation Adjustment, least squares ambiguity descending correlation adjustment method) method to fix in conjunction with partial ambiguity fixing strategy. LAMBDA (Least-square Ambiguity Decorrelation Adjustment, least squares ambiguity descending correlation adjustment method) method belongs to prior art, and it is easy to realize for those skilled in the art that it is fixed in conjunction with partial ambiguity fixing strategy, and the present embodiment is only to provide such a preferred mode.

[0200] Although the present invention has been described in detail above by general description and specific embodiments, it is obvious to those skilled in the art that some modifications or improvements can be made to the present invention. Therefore, these modifications or improvements made without departing from the spirit of the present invention all belong to the scope of protection claimed by the present invention.

Claims

1. A GNSS multi-system robust adaptive fusion RTK solution method, characterized in that: The following steps are involved: Step 1: Establish the conventional RTK double-difference observation equation based on the BDS / GPS / GLONASS / GLILEO multi-system observation data of the base station and the rover; Step 2: Conventional Kalman filter solution; Step 3: Robust Kalman adaptive filter estimation; Step 4: Helmert variance component estimation of BDS / GPS / GLONASS / GLILEO multi-system fusion; Step 5: BDS / GPS / GLONASS / GLILEO multi-system RTK ambiguity fixation; The specific process of step 1 is as follows: Step 1-1: Establish the basic observation equation based on the undifferenced raw pseudorange and carrier phase observations, as follows: In the formula, the subscripts i and r are the frequency and station of the observed quantity, respectively, and the superscript j is the satellite; are the pseudorange and carrier phase raw observation values, in meters, respectively. is the phase observation value, in weeks; is the geometric distance between the satellite and the station at the time of signal transmission, c is the speed of light in vacuum, dt r and dt j are the receiver clock error and satellite clock error respectively, T j is the tropospheric delay, is the ionospheric delay at frequency i, b Pi,r and are the pseudorange hardware delay bias at the receiver and satellite, respectively, and b φi,r and are the phase hardware delay deviations of the receiver and satellite, including the carrier phase hardware delay deviation and the initial phase deviation, λ i and is the carrier phase wavelength and integer ambiguity of frequency i, and are the measurement noise of pseudorange and carrier phase observations, respectively; Step 1-2: According to the RTK double-difference model, the satellite clock error and receiver clock error in step 1-1 can be eliminated, and the influence of satellite orbit error, tropospheric delay error, ionospheric delay error and multipath effect can be weakened, and the carrier phase double-difference observation equation can be established, as follows: In the formula, superscript j and k represent satellite numbers; subscript r and b represent mobile stations and base stations; is the double difference phase observation value at epoch t; is the distance from satellite s to station r; is the double difference ambiguity; since the reference coordinates are known, the above formula can be transformed into: Step 1-3: The prior coordinates of the mobile station are estimated using single point positioning, and the approximate estimate X r =(x r ,y r ,z r ) T , let the correction number be δX r =(δx r ,δy r ,δz r ) T , linearizing the final transformed equation in step 1-2 yields the carrier phase double difference error equation: In the formula, Similarly, the pseudorange double difference error equation can be obtained: in, is the double-difference pseudorange observation value; The specific process of step 2 is: Step 2-1: Establish the state equation and observation equation of Kalman filter: X k =Φ k,k-1 X k-1 +W k-1 L k =A k X k +V k Where k represents the observation epoch time, X k represents the state vector at time k, Φ k,k-1 represents the transition matrix of the system state from time k-1 to time k, W k-1 is the dynamic noise vector, L k is the observation vector at time k, A k is the observation equation coefficient matrix, V k is the observation noise vector; Step 2-2: The observations are considered to be independent of each other. The dynamic noise and observation noise are zero-mean, uncorrelated white noises, which satisfy the following conditions: In the formula, Ω k is the variance matrix of process noise, Q k is the variance matrix of the observation noise; Step 2-3: Rewrite the state equation and observation equation into error equations respectively: In the formula, and V k are the predicted values ​​of state parameters and the observed value L k The residual vector of the state parameter is for: The covariance matrix of is: Construct the objective function according to the least squares criterion: Where P k and are the observed values ​​L k and state parameter forecast values The weight matrix of Paired Partial derivative and set its partial derivative equal to 0, we have: Further we can get: Step 2-4: Substitute the error equation corresponding to the observation equation into formula (1) to obtain: By using the matrix identity transformation and matrix inversion formula, the Kalman filter calculation process is further derived as follows: Step 2-4-1: Calculate the Kalman filter gain matrix K based on the prediction accuracy of the previous epoch for the current epoch k , Step 2-4-2: Update the current epoch observations, i.e. use the weights of the current epoch measurements for weighted correction. Step 2-4-3: Update time, i.e. predict the observation value of the next epoch; The specific process of step 3 is as follows: Step 3-1: The robust adaptive Kalman filter solution is: In the formula, For L k The equivalent weight matrix of , which is the adaptive estimate of the observation vector weight matrix, α k is the adaptive factor, 0≤α k ≤1; Step 3-2: Determine equivalent weight factors and adaptive weight factors; The specific process of step 4 is as follows: Step 4-1: During the first filtering, the BDS, GPS, GLONASS and GLILEO observations are weighted according to the satellite elevation angle. C :P G :P R :P E ; Step 4-2: Perform filtering and solve to obtain V i T P i V i , i=C,G,R,E; Step 4-3: Estimation of variance components: In the formula, N=A T PA=A C T P C A C +A G T P G A G +A R T P R A R +A E T P E A E =N C +N G +N R +N E Among them, n C 、n G 、n R and n E are the numbers of BDS, GPS, GLONASS and GALILEO observations, respectively; Step 4-4: Re-evaluation: In the formula, c is any constant; Step 4-5: Repeat steps 4-2 to 4-4 until The difference is less than 0.

01.

2. A GNSS multi-system robust adaptive fusion RTK solution method according to claim 1, characterized in that: The specific calculation process of the equivalent weight factor is: Step 3-2-1: In BDS / GPS / GLONASS / GALILEO fusion positioning, the observation values ​​between systems are independent of each other, and the corresponding equivalent weights are: Where P i is the original weight matrix, w i is the adaptive weight factor; Step 3-2-2: Determine the weight factor using an IGGIII equivalent weight factor or a dual-factor equivalent weight factor function, wherein the function expression based on the IGGIII equivalent weight factor is: In the formula, is the standardized residual, k1 and k2 are constants, k1 = 1.5 ~ 2.0, k2 = 3.0 ~ 8.0; The function expression based on the double-factor equivalent weight factor is: In the formula, c is a constant, c = 2.5 ~ 3.0, i and j represent the row and column numbers of the matrix; Step 3-2-3: Calculate standardized residuals Step 3-2-4: For different types of observations, their standardized residuals for: In the formula, 3. A GNSS multi-system robust adaptive fusion RTK solution method according to claim 1, characterized in that: The specific calculation process of the adaptive weight factor is: Adopt two-segment function to calculate adaptive weight factor a k ,Right now: In the formula, c is a constant, ranging from 1.0 to 2.

5. In the formula, is the standardized prediction residual, is the state parameter prediction residual, is the trace of the state parameter prediction residual covariance matrix.

4. The GNSS multi-system robust adaptive fusion RTK solution method according to claim 1, characterized in that: When the iteration does not converge or converges distortedly in step 4, re-screen the satellites or appropriately reduce the weights of some satellites before solving the problem again.

5. The GNSS multi-system robust adaptive fusion RTK solution method according to claim 1, characterized in that: In step 5, the BDS / GPS / GLONASS / GLILEO multi-system RTK ambiguity is fixed by using the least squares ambiguity reduction correlation adjustment method combined with a partial ambiguity fixing strategy.

Images