High-precision precise point positioning method based on BDS-3

Abstract

The invention provides a high-precision precise point positioning method based on BDS-3, and relates to the technical field of wireless communication. The method comprises the following steps: firstly, introducing a traditional ionosphere-free combination model to eliminate the influence of ionosphere delay on positioning precision, and eliminating or weakening a satellite orbit error and a clock error error by using satellite orbit and clock error data provided by an IGS; secondly, constructing a troposphere correction model of adaptive weighted wet delay for eliminating troposphere delay errors; thirdly, establishing a self-adaptive random model, and reasonably weighting pseudo-range and carrier phase observed quantity; and finally, estimating the parameters to be estimated through extended Kalman filtering so as to obtain a positioning calculation result. The method has obvious advantages in the aspect of improving the positioning precision, can meet the requirement of high-precision positioning, can realize centimeter-level positioning precision no matter in a dynamic mode or a static mode, and can even reach millimeter-level positioning precision in the static mode in some observation stations and directions.

Description

Technical Field

[0001] The present invention relates to the field of wireless communication technologies, and particularly to a high-precision precise point positioning method based on BDS-3. Background Art

[0002] As one of the global satellite navigation systems, the BeiDou Navigation Satellite System (BDS) is widely used in multiple fields such as transportation, surveying and mapping, agriculture, and disaster monitoring due to its high-precision positioning ability. The Precise Point Positioning (PPP) technology is one of the core positioning methods in the BeiDou system, which has advantages such as not relying on a reference station, reducing labor intensity and production costs. Through a single receiver, the PPP technology can achieve high-precision static and dynamic positioning globally, solving the problem that traditional relative positioning methods cannot be applied in the case of long distances or lack of reference stations. However, since the PPP technology does not rely on a reference station for error correction, its positioning accuracy is affected by various error sources, especially the interference of factors such as satellite clock errors, orbit errors, and atmospheric delays. In order to improve the positioning accuracy, an effective error correction model must be established to reduce the influence of these errors on the positioning result. In addition, the PPP positioning solution depends on two types of observables, pseudorange and carrier phase, and their accuracy and stability are directly affected by the weighting of the observables. Reasonably allocating the weights of the observables and effectively eliminating various errors will directly improve its positioning accuracy and reliability. Therefore, in-depth research on the PPP technology in the BeiDou system and optimizing the error correction model have important theoretical significance and practical value for improving the positioning accuracy of the BeiDou system, expanding its application scope, and enhancing its global service capabilities. Summary of the Invention

[0003] The technical problem to be solved by the present invention is to provide a high-precision precise point positioning method based on BDS-3 in view of the above-mentioned prior art deficiencies.

[0004] To solve the above technical problem, the technical solution adopted by the present invention is as follows:

[0005] A high-precision precise point positioning method based on BDS-3. First, introduce the traditional ionosphere-free combination model to eliminate the influence of ionospheric delay on positioning accuracy, and use the satellite orbit and clock data provided by IGS to eliminate or weaken satellite orbit errors and clock errors. Secondly, construct an adaptive weighted wet delay tropospheric correction model to eliminate tropospheric delay errors. Then, establish an adaptive stochastic model to reasonably weight the pseudorange and carrier phase observables. Finally, estimate the parameters to be estimated through the extended Kalman filter to obtain the positioning solution result.

[0006] Furthermore, in the ionosphere-free combined model, by constructing the phase combination and pseudorange combination of the ionosphere-free delay, the observations of two frequencies are weighted and combined to cancel the interference of the ionosphere delay. The observations of the ionosphere-free model phase and pseudorange are respectively:

[0007]

[0008]

[0009] Wherein, and are respectively the pseudorange and carrier phase observations of the ionosphere-free combination; the superscript s represents the satellite, the subscript r represents the station position, and the subscript IF represents the ionosphere-free combination; is the geometric distance between the satellite position and the station position; c is the speed of light in vacuum; dt r is the receiver clock error; dt s is the satellite clock error; is the tropospheric delay of the signal propagation path; is the noise corresponding to the ionosphere-free delay pseudorange combination of station r and satellite s; is the noise corresponding to the ionosphere-free delay phase combination of station r and satellite s; is the phase ambiguity.

[0010] Furthermore, the tropospheric correction model of the adaptive weighted wet delay is specifically as follows:

[0011] The tropospheric delay on the signal propagation path is expressed as:

[0012]

[0013] Wherein, Z H and Z W are respectively the tropospheric dry delay and the tropospheric wet delay; and are respectively the projection functions corresponding to the dry delay component and the wet delay component, as follows:

[0014]

[0015] Wherein, h is the vertical distance of the station; is the elevation angle in the direction of the line connecting the station and the satellite; a ht , b ht , c ht are the constants determined by the NMF model; a d , b d , c d are obtained by interpolating the NMF dry component coefficients; a W , bW , c W Obtained by interpolating the NMF wet component coefficient; introducing the northward horizontal gradient G of the wet delay N and the eastward horizontal gradient G E On this basis, further dynamically adjust the weight of the wet delay gradient, and the projection function of the wet delay Is further corrected to:

[0016]

[0017] Among them, Is the azimuth angle in the direction of the line connecting the station and the satellite; W h (RH) is the humidity weight function, and its expression is:

[0018]

[0019] Among them, α and β are adjustment parameters used to flexibly control the change range of the humidity weight; RH represents humidity;

[0020] In the humidity weight function W h (RH), when the humidity is high, that is, RH>80%, increase the weight of the wet delay gradient, and the adjustment factor α is set to 0.05 - 0.15: If the humidity is high but not extremely wet, that is, the humidity is around 80%, the adjustment factor is set to 0.05; if the humidity is higher than 80%, the influence of the wet delay on signal propagation is greater, and stronger correction is needed to compensate for the delay change caused by humidity, and the higher adjustment factor α = 0.05 + 0.1×[(RH - 80) / 20] is used;

[0021] When the humidity is low, that is, RH<30%, reduce the weight of the wet delay gradient, and the adjustment factor β is set to 0.05 - 0.1; when the humidity is very low, that is, below 20%, the influence of the wet delay on signal propagation is extremely small, and a smaller adjustment factor is selected, β is set to 0.05; if the humidity is on the low side, that is, between 20% - 30%, although the influence of the wet delay is small, it still exists. At this time, appropriately increase the adjustment factor, β = 0.05 + 0.05×[(30 - RH) / 30];

[0022] Under medium humidity conditions, that is, 30≤RH≤80%, keep the weight as the reference value 1.

[0023] After synthesizing the above components, the tropospheric delay on the final signal propagation path Is calculated by the following formula:

[0024]

[0025] Among them, Z T Represents the tropospheric delay in the zenith direction;

[0026] Obtained by the projection function model NMF The dry component delay Z in the zenith direction H Calculated by the Saastamoinen model according to the standard meteorological parameters; the tropospheric delay Z in the zenith direction T And the northward horizontal gradient G of the wet delay N And the eastward horizontal gradient G E Are regarded as unknown parameters for estimation;

[0027] After processing the tropospheric delay error, linearize equations (1) and (2) to obtain the simplified observation equation, and the specific expression is:

[0028]

[0029] Where, L Φ And L P Are the difference vectors between the phase and pseudorange observation values and the calculated values; A is the design matrix; X is the parameter to be estimated, including the station position r, the distance error cdt corresponding to the receiver clock error r The tropospheric delay Z in the zenith direction T And the horizontal gradient G N G E ; Y is the phase ambiguity; e is the identity matrix; v Φ And v P Are the residuals of the phase and pseudorange observations;

[0030] After solving the observation equation (9), the parameters to be estimated for joint positioning using pseudorange and phase observations are obtained as follows:

[0031]

[0032] Where, And Respectively represent the estimated values of X and Y; Is the corresponding weight matrix of the pseudorange and phase observation values; w Φ And w P Are determined according to the PPP stochastic model.

[0033] Furthermore, for the establishment of the adaptive stochastic model, first evaluate the observation accuracy of GEO, IGSO, and MEO satellites through the residual values, and then combine the satellite elevation angle model to construct the adaptive stochastic model; the specific method is as follows:

[0034] In the ionosphere-free combination model, subtract the ionosphere-free phase and pseudorange combinations to obtain:

[0035]

[0036] Combined with the standard deviations of the phase and pseudorange observations, calculate the delay standard deviation of the ionosphere-free phase and pseudorange combination difference to quantify the magnitude and distribution of the residual errors in the difference; given the phase standard deviation and the pseudorange standard deviation in the case, through the error propagation law, synthesize the errors of the phase and pseudorange, and thus calculate the root mean square error corresponding to their difference as:

[0037]

[0038] where, is the root mean square error of the ionosphere-free phase and pseudorange combination difference corresponding to station r and satellite s;

[0039] In the case where no cycle slips occur in u observation epochs, the ionosphere-free combination can effectively reduce the influence of the ionospheric error; at this time, the difference between the phase and pseudorange observations, that is, the residual, is caused by the system noise and other error sources, and the residual related to this difference is described by the following expression:

[0040]

[0041] where, is the residual corresponding to the ionosphere-free delay phase and pseudorange combination difference;

[0042] Quantify the overall accuracy of the satellite observation data by calculating the root mean square error value; for the residual related to the ionosphere-free combination difference, the root mean square error value is used to evaluate the error magnitude of the phase and pseudorange difference; for the ionosphere-free combination, the root mean square error of the residual is calculated by the following expression:

[0043]

[0044] Calculate the average root mean square errors of GEO, IGSO, and MEO satellites respectively, and measure the accuracy of each type of satellite data based on these errors, as shown in the following formula:

[0045]

[0046] where, t represents the type of Beidou satellite, that is, GEO, IGSO, or MEO; σ t ′ is the average root mean square error value of one of the three types of Beidou satellites; k represents the number of observations of the t-th type of satellite;

[0047] According to Equation (15), if the root mean square error value of MEO is the largest, that is, the observation accuracy of MEO satellites is the lowest, the variance of the stochastic model is expressed as:

[0048]

[0049] Among them, σ is solved according to the elevation angle, as shown in the following formula:

[0050] σ 2 = a 2 + b 2 / (sinE) 2 (17)

[0051] Among them, both a and b are set to 0.003m; E represents the elevation angle;

[0052] According to Equation (15), if the observation accuracy of the IGSO satellite is the lowest, then in the above formula (16), when t = IGSO, when t = GEO, MEO,

[0053] According to Equation (15), if the observation accuracy of the GEO satellite is the lowest, then in the above formula (16), when t = GEO, when t = IGSO, MEO,

[0054] By analyzing the observation value accuracies of different types of satellites GEO, IGSO, and MEO, different weights are set for each type of satellite to optimize the stochastic model in the positioning solution process; the adaptive stochastic model combined with the observation value variance factor and variance is:

[0055]

[0056] Among them, represents the stochastic model, and its subscripts Φ and P are the phase observation value and the pseudorange observation value respectively; the variance factor R r is the weight ratio of the phase observation value to the pseudorange observation value, R r = σ b / σ c σ b σ c are the root mean square error values of the carrier phase observation value residual v b and the pseudorange observation value residual v c respectively;

[0057] Thus, the corresponding weight matrices of the pseudorange and phase observation values in Equation (10) are:

[0058]

[0059] Furthermore, the specific method for estimating the parameters to be estimated by extended Kalman filtering is:

[0060] After obtaining the weight matrix of the phase pseudorange observation value through formula (19), in the PPP mode, the extended Kalman filter is introduced to estimate the parameters to be estimated in formula (10). It can be seen from formula (10) that the vector to be estimated is:

[0061]

[0062] The observation data required for the extended Kalman filter estimation are: the phase observation Φ IF and the pseudorange observation P IF , and these data can provide the relationship between the system state and the observation, so as to determine the expression of the observation vector as:

[0063]

[0064] In the process of extended Kalman filter estimation, by updating the state vector and covariance matrix to be estimated, the state estimation is continuously corrected and improved according to the current observation information at each epoch; at epoch t k the state vector to be estimated and the corresponding covariance matrix P k can be estimated using the observation vector y k at the same moment. The model is:

[0065]

[0066] In the formula, P k (+) is the covariance matrix updated by the observation data at epoch t k ; P k (-) is the covariance matrix before the update of the observation data at epoch t k ; K k is the Kalman gain matrix at epoch t k , which controls the balance between prediction and measurement; h() is the observation matrix, representing the relationship between the state variable and the observation value; H() is the partial derivative matrix; R k is the covariance matrix of the measurement error at epoch t k ; and are the state vectors before and after the update at epoch t k respectively;

[0067] H(x) is the partial derivative matrix of the observation model with respect to the state vector, describing the influence of the change of the state vector on the observation value and reflecting the linearized relationship between the state and the observation; the form of H(x) is:

[0068]

[0069] Among them, each term represents the partial derivative of the observed value with respect to the state variable, specifically as follows:

[0070]

[0071] In the formula, is the unit vector in the direction from the station to the satellite, s = 1, 2, … m; the superscript T represents the transpose;

[0072]

[0073] Among them, is the projection function of the tropospheric wet delay in the zenith direction.

[0074] Furthermore, the extended Kalman filtering method updates the state vector and covariance matrix at each stage, so as to more accurately estimate the current state of the system; specifically, the formulas for updating the state vector and the associated covariance matrix are described as:

[0075]

[0076] Among them, is the system noise transfer matrix from the epoch t k to the epoch t k+1 ; is the system noise covariance matrix from the epoch t k to the epoch t k+1 .

[0077] Furthermore, in the positioning of the high-precision precise point positioning method based on BDS-3, in the static mode, the station position is assumed to be a fixed value for estimation, and the receiver clock error is regarded as white noise and its mean value is taken as 0 to reduce the influence of the clock error on the positioning accuracy; while in the dynamic mode, the station position is estimated through a random walk process, and the position estimate is updated in real time to meet the positioning requirements in the dynamic environment.

[0078] The beneficial effects of adopting the above technical solutions are as follows: The high-precision precise point positioning method based on BDS-3 provided by the present invention has obvious advantages in improving the positioning accuracy, can meet the requirements of high-precision positioning, and can achieve centimeter-level positioning accuracy both in the dynamic mode and the static mode. And in some stations and directions, the positioning accuracy in the static mode can even reach the millimeter level. BRIEF DESCRIPTION OF THE DRAWINGS

[0079] Figure 1 is the schematic diagram of the high-precision precise point positioning method based on BDS-3 provided by the embodiment of the present invention;

[0080] Figure 2Statistical chart of positioning errors in the E direction of three precise point positioning methods at URUM, SGOC, CAS1, and WUH2 stations in the static mode provided by the embodiments of the present invention;

[0081] Figure 3 Statistical chart of positioning errors in the N direction of three precise point positioning methods at URUM, SGOC, CAS1, and WUH2 stations in the static mode provided by the embodiments of the present invention;

[0082] Figure 4 Statistical chart of positioning errors in the U direction of three precise point positioning methods at URUM, SGOC, CAS1, and WUH2 stations in the static mode provided by the embodiments of the present invention;

[0083] Figure 5 Error comparison chart of three dynamic precise point positioning methods at the URUM station provided by the embodiments of the present invention;

[0084] Figure 6 Error comparison chart of three dynamic precise point positioning methods at the SGOC station provided by the embodiments of the present invention;

[0085] Figure 7 Error comparison chart of three dynamic precise point positioning methods at the WUH2 station provided by the embodiments of the present invention;

[0086] Figure 8 Error comparison chart of three dynamic precise point positioning methods at the CAS1 station provided by the embodiments of the present invention. Detailed implementation manners

[0087] The following combines the accompanying drawings and embodiments to further describe in detail the specific implementation manners of the present invention. The following embodiments are used to illustrate the present invention but are not used to limit the scope of the present invention.

[0088] This embodiment provides a high-precision precise point positioning method based on BDS-3. In this method, aiming at the problems of uneven distribution of tropospheric wet delay and humidity differences under different conditions, an adaptive weighted wet delay tropospheric correction model is constructed. Based on the ZTD parameter estimation model, this model introduces the northward and eastward horizontal gradients of the wet delay, corrects the projection function of the wet delay, and constructs a humidity weight function to further dynamically adjust the projection function of the wet delay. In addition, in this method, considering that the traditional stochastic model fails to fully reflect the differences in the observation accuracies of different types of BDS-3 satellites, an adaptive stochastic model is established based on the accuracy evaluation results of BDS-3 observation data to further improve the positioning accuracy of the high-precision precise point positioning method based on BDS-3. As Figure 1 shown, the method of this embodiment is described as follows.

[0089] Step 1: The high-precision precise point positioning method based on BDS-3 introduces the traditional ionosphere-free combination model, which can effectively eliminate the influence of ionospheric delay on positioning accuracy during signal propagation. In this model, by constructing the ionosphere-free delay phase combination and pseudorange combination, the observations of two frequencies are weighted and combined to offset the interference of ionospheric delay. The observations of the ionosphere-free model phase and pseudorange are as follows:

[0090]

[0091]

[0092] In the formula, and are the pseudorange and carrier phase observations of the ionosphere-free combination respectively; the superscript s represents the satellite, the subscript r represents the station coordinates, and the subscript IF represents the ionosphere-free combination; is the geometric distance between the satellite position and the station position; c is the speed of light in vacuum; dt r is the receiver clock error; dt s is the satellite clock error; is the tropospheric delay of the signal propagation path; is the noise corresponding to the ionosphere-free delay pseudorange combination of station r and satellite s; is the noise corresponding to the ionosphere-free delay phase combination of station r and satellite s; is the phase ambiguity.

[0093] It can be seen from Equation (1) and Equation (2) that although the traditional ionosphere-free model eliminates the influence of the ionosphere, there are still satellite clock error dt s , receiver clock error dt r , observation noise and as well as the influence of tropospheric delay error and so on. Therefore, in the data processing of precise point positioning, the satellite orbit and clock data provided by IGS are adopted to eliminate or weaken the satellite orbit error and clock error. Thus, eliminating the influence of tropospheric delay error becomes a key issue.

[0094] Step 2: Considering the non-uniform distribution of the tropospheric wet delay in the tropospheric delay error and the influence brought by the humidity difference under different climate conditions, based on the ZTD parameter estimation model, an adaptive weighted wet delay tropospheric correction model is constructed.

[0095] The tropospheric delay on the signal propagation path is expressed as:

[0096]

[0097] In the formula, Z H and Z W are the tropospheric dry delay and the tropospheric wet delay respectively; and are the projection functions corresponding to the dry and wet delay components, as shown in the following two formulas:

[0098]

[0099]

[0100] where h is the vertical distance of the station; is the elevation angle in the direction of the line connecting the station and the satellite; a ht , b ht , c ht are constants determined by the NMF model; a d , b d , c d are obtained by interpolating the NMF dry component coefficients; a W , b W , c W are obtained by interpolating the NMF wet component coefficients.

[0101] In order to correct the uneven distribution of the tropospheric wet delay in different directions and its tropospheric delay distribution under different climatic conditions, on the basis of introducing the northward horizontal gradient G N and the eastward horizontal gradient G E of the wet delay, the weight of the wet delay gradient is further adjusted dynamically. For this, the projection function of the wet delay is further corrected to:

[0102]

[0103] where is the azimuth angle in the direction of the line connecting the station and the satellite; W h (RH) is the humidity weight function, and its expression is:

[0104]

[0105] In the formula, α and β are adjustment parameters used to flexibly control the variation range of the humidity weight. When the humidity RH is high (RH > 80%), the weight of the wet delay gradient is increased, and the adjustment factor α is set to 0.05 - 0.15; if the humidity is high but not extremely wet (humidity around 80%), the adjustment factor α is set to 0.05. If the humidity is higher than 80% and the influence of the wet delay on signal propagation is large, stronger correction is needed to compensate for the delay change caused by humidity, and a higher adjustment factor is used, α = 0.05 + 0.1×[(RH - 80) / 20]. When the humidity RH is low (RH < 30%), the weight of the wet delay gradient is decreased, and the adjustment factor β is set to 0.05 - 0.1; when the humidity is very low (such as below 20%), the influence of the wet delay on signal propagation is extremely small, and excessive correction may introduce unnecessary errors, so a smaller adjustment factor is selected, β is set to 0.05; if the humidity is on the low side (between 20% - 30%), although the influence of the wet delay is small, it still exists. At this time, the adjustment factor can be appropriately increased, β = 0.05 + 0.05×[(30 - RH) / 30], to better correct the error caused by humidity. Under medium humidity conditions (30 ≤ RH ≤ 80%), the weight is maintained at the reference value of 1.

[0106] After synthesizing the above components, the tropospheric delay on the final signal propagation path can be calculated by the following formula:

[0107]

[0108] where Z T represents the tropospheric delay in the zenith direction.

[0109] During the process of high-precision precise point positioning based on BDS-3, obtained by the projection function model NMF, the dry component delay Z in the zenith direction H is calculated by the Saastamoinen model according to standard meteorological parameters. The tropospheric delay Z in the zenith direction T and the northward horizontal gradient G of the wet delay N and the eastward horizontal gradient G E are regarded as unknown parameters for estimation in order to more accurately correct the influence of the tropospheric delay on the positioning result.

[0110] After processing the tropospheric delay error, linearization is performed on equations (1) and (2) to simplify the process of solving using the observation equations. Through linearization, the simplified observation equations are obtained, and the specific expressions are:

[0111]

[0112] where L Φ and L P are the difference vectors between the phase and pseudorange observations and the calculated values; A is the design matrix; X is the parameter to be estimated, including the station position r, the distance error cdt corresponding to the receiver clock error r , the tropospheric delay Z in the zenith direction T , and the horizontal gradient G N , G E ; Y is the phase ambiguity; e is the identity matrix; v Φ and v P are the residuals of the phase and pseudorange observations.

[0113] After solving the observation equation (9), the parameters to be estimated for joint positioning using pseudorange and phase observations can be obtained:

[0114]

[0115] where and represent the estimated values of X and Y respectively; is the corresponding weight matrix of the pseudorange and phase observations, and w Φ and w P are determined according to the PPP stochastic model.

[0116] Step 3: As can be seen from Equation (10), the key to solving the parameter to be estimated lies in accurately calculating the corresponding weight w of the pseudorange observation P and the corresponding weight w of the phase observation Φ . Therefore, it is crucial to establish a reasonable stochastic model. Considering the limitation that the traditional stochastic model fails to fully reflect the differences in the observation accuracies of different satellites, based on the accuracy evaluation results of BDS-3 observation data, an adaptive stochastic model is proposed for weighting the pseudorange and carrier phase observations. This method first evaluates the observation accuracies of GEO, IGSO, and MEO satellites through the residual values, and then constructs an adaptive stochastic model in combination with the satellite elevation angle model.

[0117] In the ionosphere-free combination model, by subtracting the phase and pseudorange combinations, the influence of the ionospheric delay can be effectively eliminated, thereby obtaining a relatively pure expression of the observation value difference, which helps to more accurately evaluate various observation errors. The subtraction of the ionosphere-free phase and pseudorange combinations can obtain:

[0118]

[0119] Furthermore, by combining the standard deviations of the phase and pseudorange observations, the delay standard deviation of the ionosphere-free phase and pseudorange combination difference can be calculated to quantify the magnitude and distribution of the residual errors in the difference. Given the phase standard deviation and the pseudorange standard deviation In this case, through the error propagation law, the errors of the phase and the pseudorange are synthesized, and the root mean square error corresponding to their difference can be calculated as follows:

[0120]

[0121] In the case where no cycle slips occur in u observation epochs, the ionosphere-free combination can effectively reduce the influence of the ionospheric error. At this time, the difference (i.e., the residual) between the phase and the pseudorange observations is usually caused by system noise and other error sources. The residual related to this difference can be described by the following expression:

[0122]

[0123] where is the residual corresponding to the difference between the ionosphere-free delay phase and the pseudorange combination.

[0124] The overall accuracy of satellite observation data is quantified by calculating the root mean square error value. For the residual related to the ionosphere-free combination difference, the root mean square error value can be used to evaluate the error magnitude of the phase and pseudorange difference. For the ionosphere-free combination, the root mean square error of the residual can be calculated by the following expression:

[0125]

[0126] Next, the average root mean square errors of the GEO, IGSO, and MEO satellites are calculated respectively, and the accuracy of each type of satellite data is measured based on these errors:

[0127]

[0128] In the formula, t represents the type of Beidou satellite, i.e., GEO, IGSO, or MEO; σ t ′ is the average root mean square error value of one of the three types of Beidou satellites; k represents the number of observations of the t-th type of satellite.

[0129] According to Equation (15), if the root mean square error value of MEO is the largest, that is, the observation accuracy of MEO satellites is the lowest, the variance of the stochastic model can be expressed as:

[0130]

[0131] where σ is solved according to the elevation angle, as shown in the following formula:

[0132] σ 2 = a 2 + b 2 / (sinE) 2 (17)

[0133] Among them, both a and b are set to 0.003 m; E represents the elevation angle.

[0134] According to Equation (15), if the observation accuracy of the IGSO satellite is the lowest, then in the above formula (16), when t = IGSO, when t = GEO, MEO,

[0135] According to Equation (15), if the observation accuracy of the GEO satellite is the lowest, then in the above formula (16), when t = GEO, when t = IGSO, MEO,

[0136] By analyzing the observation value accuracies of different types of satellites GEO, IGSO, and MEO, different weights are set for each type of satellite to optimize the stochastic model in the positioning solution process. The adaptive stochastic model combined with the observation value variance factor and variance is:

[0137]

[0138] Among them, represents the stochastic model, and its subscripts Φ and P are the phase observation value and the pseudorange observation value respectively; the variance factor R r is the weight ratio of the phase observation value to the pseudorange observation value, R r = σ b / σ c σ b σ c are the root mean square error values of the carrier phase observation value residual v b and the pseudorange observation value residual v c respectively.

[0139] Thus, the corresponding weight matrices of the pseudorange and phase observation values in Equation (10) are:

[0140]

[0141] Step 4: After obtaining the weight matrix of the phase pseudorange observation value through Equation (19), in the PPP mode, the extended Kalman filter (EKF) is introduced to estimate the parameters to be estimated in Equation (10). It can be seen from Equation (10) that the vector to be estimated is:

[0142]

[0143] The observation data required for EKF filtering estimation are: the phase observation quantity Φ IF and the pseudorange observation quantity P IF, these data can provide the relationship between the system state and the observed quantities, thereby determining the expression of the observation vector as:

[0144]

[0145] In the EKF filtering estimation process, by updating the state vector to be estimated and the covariance matrix, the state estimation can be continuously corrected and improved according to the current observation information at each epoch. At time t k The state vector to be estimated at epoch and the corresponding covariance matrix P k can be estimated using the observation vector y k at the same moment. The model is:

[0146]

[0147] where P k (+) is the covariance matrix updated at time t k for the observed data; P k (-) is the covariance matrix before updating the observed data at time t k ; K k is the Kalman gain matrix at time t k , which controls the balance between prediction and measurement; h(·) is the observation matrix, representing the relationship between the state variable and the observed value; H(·) is the partial derivative matrix; R k is the covariance matrix of the measurement error at time t k ; and are the state vector before updating and the state vector after updating at time t k respectively.

[0148] H(x) is the partial derivative of the observation model with respect to the state vector. It describes the influence of the change in the state vector on the observed value and reflects the linearized relationship between the state and the observation. Specifically, the form of H(x) is:

[0149]

[0150] In the partial derivative matrix H(x), each term represents the partial derivative of the observed value with respect to the state variable, specifically as follows:

[0151]

[0152] where is the unit vector in the direction of the satellite from the measuring station.

[0153]

[0154]

[0155] In the formula, is the projection function of the tropospheric wet delay in the zenith direction, and the meanings of other symbols are the same as those in formula (6).

[0156] The calculation and determination of each term in H(x) depend on the parameters in the state vector x to be estimated. For example, M T the elements in the matrix (the projection function of the tropospheric wet delay in the zenith direction) is related to parameters such as the tropospheric delay in the zenith direction and the tropospheric horizontal gradient. Because the calculation of the tropospheric delay and the determination of the projection function involve these tropospheric-related parameters, and these parameters are included in the state vector x to be estimated.

[0157] The EKF filtering method updates the state vector and covariance matrix at each stage, so as to more accurately estimate the current state of the system. Specifically, the formulas for updating the state vector and the associated covariance matrix are described as:

[0158]

[0159] In the formula, is the system noise transfer matrix from the epoch time t k to the time t k+1 ; is the system noise covariance matrix from the time t k to the time t k+1 .

[0160] In the high-precision precise point positioning based on BDS-3, in the static mode, the station position is usually assumed to be a fixed value for estimation, and the receiver clock offset is regarded as white noise and its mean value is taken as 0 to reduce the influence of the clock offset on the positioning accuracy. In the dynamic mode, the station position is estimated through a random walk process, and the position estimate is updated in real time to meet the positioning requirements in the dynamic environment.

[0161] To verify the effectiveness and superiority of the high-precision precise point positioning method based on BDS-3 in this embodiment, a series of positioning simulation experiments are carried out, aiming to evaluate the ability of this technology to accurately calculate the station position coordinates. Multiple stations in the MGEX network are selected for the experiment, the observation time is 24 hours on December 1, 2023, the epoch interval is 30 seconds, and a total of 2880 epochs of data are included. To achieve high-precision PPP positioning, this experiment uses the 30-second precise ephemeris and satellite clock offset products provided by a university IGS data center, and uses the BDS-3 system for positioning solution. The specific configuration is as follows: the satellite elevation angle is set to 15°, the BDS-3B1 / B3 ionospheric elimination combined observation value is selected, and the receiver sampling interval is set to 30 seconds.

[0162] To evaluate the positioning performance of the high-precision precise point positioning method based on BDS-3, it is compared and analyzed with two currently widely used positioning methods: one is the precise point positioning method that determines weights by combining the tropospheric delay model based on the parameter estimation method with the experience ratio (Weighting according to experience ratio, WH-ER), which is called the precise point positioning WH-ER method; the other is the precise point positioning method that determines weights according to the residual ratio of the tropospheric delay combined with phase and pseudo-range observations based on the parameter estimation method (The weight is determined according to the residual ratio of phase and pseudo-range observations, WH-RPO), which is called the precise point positioning WH-RPO method. Simulation experiments in dynamic and static modes are carried out respectively for the high-precision precise point positioning method based on BDS-3, the precise point positioning WH-ER method, and the precise point positioning WH-RPO method. By comparing the positioning accuracies of different methods, the positioning performances of each method are evaluated.

[0163] Next, the positioning accuracy analysis of the high-precision positioning method based on BDS-3 in the static mode is carried out.

[0164] Figure 2 、 Figure 3 and Figure 4 show the statistical charts of the positioning errors of the three precise point positioning methods at the URUM, SGOC, CAS1, and WUH2 stations in the E, N, and U directions in the static mode. It can be intuitively seen from Figure 2 that the positioning errors of the high-precision precise point positioning method based on BDS-3 in the E, N, and U directions are all at the centimeter level, meeting the standards of positioning requirements and verifying the effectiveness of this method. At the same time, compared with the precise point positioning WH-ER method and the WH-RPO method, this method shows obvious advantages in the positioning errors in each direction. Although there is a slight increase in errors at some stations and in specific directions (such as the E direction at the URUM station), exceeding the precise point positioning WH-RPO method, generally speaking, the positioning errors of the high-precision precise point positioning method based on BDS-3 are at the lowest level in all directions. Therefore, the precise point positioning method based on BDS-3 has obvious advantages in positioning accuracy compared with traditional methods, can provide higher-precision positioning results, and meet the requirements of high-precision positioning.

[0165] Tables 1 and 2 present the improvements in the E, N, and U directions and the three-dimensional average positioning accuracy of the high-precision precise point positioning method based on BDS-3 compared to the WH-ER method and the WH-RPO method of precise point positioning in the static mode.

[0166] Table 1 Improvements in positioning accuracy of the high-precision precise point positioning method compared to the WH-ER method of precise point positioning

[0167]

[0168] As can be seen from Table 1, the high-precision precise point positioning method based on BDS-3 shows improvements in the positioning accuracy in the E, N, and U directions at all stations compared to the WH-ER method of precise point positioning. Especially at the SGOC station, the improvements in the positioning accuracy in the E, N, and U directions are all relatively significant, with an average improvement of 3.40 cm, while the average positioning improvements at other stations are 3.17 cm (URUM), 2.46 cm (CAS1), and 2.82 cm (WUH2) respectively. These results further verify the advantages of the high-precision precise point positioning method in improving positioning accuracy compared to the WH-ER method. Table 2 Average improvements in positioning accuracy of the high-precision precise point positioning method compared to the WH-RPO method of precise point positioning

[0169]

[0170]

[0171] As can be seen from Table 2, the high-precision precise point positioning method based on BDS-3 shows varying degrees of improvements in the positioning accuracy at each station compared to the WH-RPO method of precise point positioning. Although there is a slight decrease in accuracy in the E direction (-0.17 cm), the accuracy in the N and U directions is significantly improved, and the overall average positioning accuracy is improved by 2.01 cm. Other stations, such as SGOC, CAS1, and WUH2, also show varying degrees of improvements, among which the SGOC station has the largest average improvement, reaching 2.38 cm. These results show that the high-precision precise point positioning method of this embodiment has relatively consistent accuracy advantages compared to the WH-RPO method.

[0172] Next, the positioning accuracy analysis of the high-precision positioning method based on BDS-3 in the dynamic mode is carried out.

[0173] Figure 5 、 Figure 6 、 Figure 7 and Figure 8It shows the error comparison of three dynamic precise point positioning methods at URUM, SGOC, WUH2 and CAS1 stations. It can be clearly seen from the figure that the error fluctuation of the positioning curve of the high-precision precise point positioning method based on BDS-3 is significantly smaller than that of the other two methods. Especially at the initial epoch of URUM and SGOC stations, the initial epoch of WUH2 station, and the epoch around 1600, the positioning error fluctuation of the high-precision precise point positioning method based on BDS-3 is small, showing a more stable characteristic. In addition, at the CAS1 station, the positioning error fluctuation range of the high-precision precise point positioning method based on BDS-3 is significantly lower than that of the precise point positioning WH-ER method and the precise point positioning WH-RPO method, further verifying the superiority of this method in dynamic positioning.

[0174] Table 3 shows the positioning accuracy of the three dynamic precise point positioning methods in different directions. Through statistical analysis, the following conclusions can be drawn. At the URUM station, the three-dimensional average positioning errors of the dynamic high-precision precise point positioning method are 3.70 cm and 2.42 cm lower than those of the dynamic precise point positioning WH-RPO method and the WH-ER method respectively, which fully shows that the dynamic high-precision method has significant advantages in positioning accuracy. For the SGOC station, the improvement amplitude of the dynamic high-precision precise point positioning is more significant, which are 3.68 cm and 2.31 cm respectively, further verifying the superiority of this method in improving positioning accuracy. At the same time, at the WUH2 station, the dynamic high-precision precise point positioning method is 4.99 cm and 1.58 cm higher than the precise point positioning WH-RPO method and the WH-ER method respectively, and this result also shows the excellent performance of this method in improving positioning accuracy. Based on the above results, it can be concluded that the dynamic high-precision precise point positioning method has higher positioning accuracy than the precise point positioning WH-RPO method and the WH-ER method, showing its obvious advantages in accuracy improvement.

[0175] Table 3 Positioning accuracy of dynamic precise point positioning

[0176]

[0177]

[0178] In summary, whether in the dynamic mode or the static mode, the high-precision precise point positioning method based on BDS-3 can achieve centimeter-level positioning accuracy, and in some stations and directions, the positioning accuracy in the static mode can even reach the millimeter level. This result fully verifies the effectiveness of the method. In addition, compared with the precise point positioning WH-RPO method and WH-ER method, the high-precision precise point positioning method shows better positioning accuracy in all experiments, further demonstrating its advantages over traditional positioning methods. Therefore, based on the experimental results, it can be concluded that the high-precision precise point positioning method based on BDS-3 has obvious advantages in improving positioning accuracy and can meet the requirements of high-precision positioning.

[0179] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and are not intended to limit them; although the present invention has been described in detail with reference to the foregoing embodiments, those of ordinary skill in the art should understand that they can still modify the technical solutions described in the foregoing embodiments, or perform equivalent replacements for some or all of the technical features; and 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 high-precision precise single-point positioning method based on BDS-3, characterized by: Firstly, the traditional ionospheric-free combination model is introduced to eliminate the influence of ionospheric delay on positioning accuracy. The satellite orbit and clock error data provided by IGS are used to eliminate or weaken the satellite orbit error and clock error. Secondly, an adaptive weighted wet delay tropospheric correction model is constructed to eliminate the tropospheric delay error. Then, an adaptive random model is established to reasonably weight the pseudorange and carrier phase observations. Finally, the parameters to be estimated are estimated through the extended Kalman filter to obtain the positioning solution result.

2. According to claim 1, a high-precision precise single-point positioning method based on BDS-3 is characterized in that: In the ionosphere-free combination model, the observation values ​​of the two frequencies are weightedly combined by constructing a phase combination and a pseudorange combination without ionosphere delay, thereby offsetting the interference of ionosphere delay; the observation values ​​of the phase and pseudorange of the ionosphere-free model are respectively: in, and are the pseudorange and carrier phase observation values ​​of the ionosphere-free combination, respectively; the superscript s represents the satellite, the subscript r represents the station location, and the subscript IF represents the ionosphere-free combination; is the geometric distance between the satellite position and the station position; c is the speed of light in vacuum; dt r is the receiver clock error; dt s is the satellite clock error; is the tropospheric delay of the signal propagation path; is the corresponding noise of the pseudorange combination without ionospheric delay for station r and satellite s; is the corresponding noise of the phase combination without ionospheric delay at station r and satellite s; is the phase ambiguity.

3. The high-precision precise single-point positioning method based on BDS-3 according to claim 2 is characterized in that: The tropospheric correction model of the adaptive weighted wet delay is as follows: Tropospheric delay along the signal propagation path It is expressed as: Among them, Z H and Z W They are the tropospheric dry delay and the tropospheric wet delay respectively; and They are the projection functions corresponding to the dry delay component and the wet delay component, as shown in the following two equations: Where h is the vertical distance of the measuring station; is the altitude angle of the line connecting the station and the satellite; a ht , b ht 、c ht A constant determined by the NMF model; d , b d 、c d Obtained by interpolation of NMF dry component coefficients; a W , b W 、c W Obtained by interpolation of NMF wet component coefficients; introducing the north-direction horizontal gradient G of wet delay N and the eastward horizontal gradient G E Based on this, the weight of the wet delay gradient is further adjusted dynamically, and the projection function of the wet delay Further corrected to: in, is the azimuth of the line connecting the station and the satellite; W h (RH) is the humidity weight function, and its expression is: Among them, α and β are adjustment parameters used to flexibly control the change range of humidity weight; RH represents humidity; The humidity weight function W h In (RH), when the humidity is high, that is, RH>80%, the weight of the wet delay gradient is increased, and the adjustment factor α is set to 0.05-0.15: If the humidity is high but not extremely humid, that is, the humidity is around 80%, the adjustment factor is set to 0.05; if the humidity is higher than 80%, the wet delay has a greater impact on signal propagation, and a stronger correction is required to compensate for the delay change caused by humidity, and a higher adjustment factor α=0.05+0.1×[(RH-80) / 20] is used; When the humidity is low, that is, RH<30%, the weight of the wet delay gradient is reduced, and the adjustment factor β is set to 0.05-0.1; when the humidity is very low, that is, below 20%, the effect of wet delay on signal propagation is minimal, and a smaller adjustment factor β is selected, and β is set to 0.05; if the humidity is low, that is, between 20% and 30%, the effect of wet delay is small but still exists, and the adjustment factor is appropriately increased at this time, β=0.05+0.05×[(30-RH) / 30]; Under moderate humidity conditions, i.e. 30 ≤ RH ≤ 80%, keep the weight at the base value of 1; After combining the above components, the final tropospheric delay on the signal propagation path is Calculated by the following formula: Among them, Z T represents the tropospheric delay in the zenith direction; Obtained by the projection function model NMF, The dry component delay Z in the zenith direction H Calculated by the Saastamoinen model based on standard meteorological parameters; zenith tropospheric delay Z T and the northerly horizontal gradient of the wet delay G N and the eastward horizontal gradient G E are considered as unknown parameters for estimation; After processing the tropospheric delay error, equations (1) and (2) are linearized to obtain the simplified observation equations. The specific expression is: Among them, L Φ and L P is the difference vector between the phase and pseudorange observations and the calculated values; A is the design matrix; X is the parameter to be estimated, including the station position r, the distance error cdt corresponding to the receiver clock error r , zenith tropospheric delay Z T and the horizontal gradient G N , G E ; Y is the phase ambiguity; e is the unit matrix; v Φ and v P is the residual of the phase and pseudorange observations; After solving the observation equation (9), the estimated parameters for joint positioning using pseudorange and phase observations are obtained as follows: x=(B T WB) -1 B T WL (10) in, and Represent the estimated values ​​of X and Y respectively; is the corresponding weight matrix of pseudorange and phase observations; w Φ and w P Determined based on the PPP stochastic model.

4. The high-precision precise single-point positioning method based on BDS-3 according to claim 3 is characterized in that: The adaptive random model is established by first evaluating the observation accuracy of GEO, IGSO and MEO satellites through residual values, and then combining the altitude angle model of the satellite to construct an adaptive random model; the specific method is as follows: In the ionosphere-free combination model, the ionosphere-free phase and pseudorange combination are subtracted to obtain: Combining the standard deviations of phase and pseudorange observations, the delay standard deviation of the ionospheric-free phase and pseudorange combination difference is calculated to quantify the size and distribution of the residual error in the difference. and pseudorange standard deviation In the case of , the error propagation law is used to synthesize the phase and pseudorange errors, so that the root mean square error corresponding to their difference is calculated as: in, is the root mean square error of the difference between the ionosphere-free phase and pseudorange combination of station r and satellite s; When no cycle slip occurs in u observation epochs, the ionospheric-free combination can effectively reduce the influence of ionospheric errors. At this time, the difference between the phase and pseudorange observations, i.e., the residual, is caused by system noise and other error sources. The residual related to the difference is described by the following expression: in, is the residual corresponding to the difference of the phase and pseudorange combination without ionospheric delay; The overall accuracy of satellite observation data is quantified by calculating the root mean square error value; for the residuals related to the ionosphere-free combination difference, the root mean square error value is used to evaluate the error size of the phase and pseudorange difference; for the ionosphere-free combination, the root mean square error of the residual is calculated using the following expression: The average root mean square error of satellite GEO, IGSO and MEO is calculated respectively, and the accuracy of each type of satellite data is measured based on these errors, as shown in the following formula: Where t represents the type of BeiDou satellite, i.e. GEO, IGSO or MEO; σ t ′ is the average root mean square error value of one of the three types of Beidou satellites; k represents the number of observations of the tth type of satellite; According to formula (15), if the root mean square error of MEO is the largest, that is, the observation accuracy of MEO satellite is the lowest, the variance of the random model is expressed as: in, σ is solved according to the altitude angle as follows: s 2 =a 2 +b 2 / (sinE) 2 (17) Where a and b are both set to 0.003m; E represents the altitude angle; According to formula (15), if the observation accuracy of the IGSO satellite is the lowest, then in the above formula (16), when t = IGSO, When t=GEO, MEO, According to formula (15), if the observation accuracy of GEO satellite is the lowest, then in the above formula (16), when t = GEO, When t=IGSO, MEO, By analyzing the accuracy of observation values ​​of different types of GEO, IGSO and MEO satellites, different weights are set for each type of satellite to optimize the random model in the positioning solution process; the adaptive random model composed of the observation variance factor and variance is: in, represents the random model, where the subscripts Φ and P are phase observations and pseudorange observations respectively; the variance factor R r is the weight ratio of phase observation value to pseudorange observation value, R r =σ b / σ c , σ b , σ c They are respectively the carrier phase observation residual v b and pseudorange observation residual v c The root mean square error value of It can be obtained that the corresponding weight matrix of pseudorange and phase observation values ​​in formula (10) is:

5. The high-precision precise single-point positioning method based on BDS-3 according to claim 4 is characterized in that: The specific method of estimating the parameters to be estimated by using the extended Kalman filter is: After the phase pseudorange observation weight matrix is ​​obtained by formula (19), the extended Kalman filter is introduced in the PPP mode to estimate the parameters to be estimated in formula (10). From formula (10), it can be seen that the vector to be estimated is: The observation data required for the extended Kalman filter estimation are: phase observation Φ IF and pseudo-range observation P IF , these data can provide the relationship between the system state and the observed quantity, thereby determining the expression of the observation vector: In the process of extended Kalman filter estimation, the state estimation is continuously corrected and improved according to the current observation information at each epoch by updating the state vector and covariance matrix to be estimated. k The state vector to be estimated at the epoch time and the corresponding covariance matrix P k Can use the observation vector y at the same time k To estimate, the model is: Where P k (+) is the observed data at t k The covariance matrix after the epoch is updated; P k (-) is the observed data at t k The covariance matrix before the epoch is updated; K k t k The Kalman gain matrix at the epoch time controls the balance between prediction and measurement; h(·) is the observation matrix, which represents the relationship between state variables and observations; H(·) is the partial derivative matrix; R k t k The covariance matrix of the measurement errors at the epoch time; and t k The state vector before and after the epoch time is updated; H(x) is the partial derivative matrix of the observation model to the state vector, which describes the impact of the change of the state vector on the observation value and reflects the linear relationship between the state and the observation. The form of H(x) is: Among them, each term represents the partial derivative of the observed value with respect to the state variable, as follows: In the formula, is the unit vector from the station to the satellite, s = 1, 2, ... m; the superscript T indicates transposition; in, is the projection function of the tropospheric wet delay in the zenith direction.

6. The high-precision precise single-point positioning method based on BDS-3 according to claim 5 is characterized in that: The extended Kalman filtering method updates the state vector and the covariance matrix at each stage, so as to estimate the current state of the system more accurately; specifically, the formula for updating the state vector and the covariance matrix associated therewith is described as: in, From t k epoch time to t k+1 The system noise transfer matrix at the epoch time; t k epoch time to t k+1 System noise covariance matrix at epoch time.

7. The high-precision precise single-point positioning method based on BDS-3 according to claim 6 is characterized in that: In the positioning of the high-precision precise single-point positioning method based on BDS-3, in the static mode, the station position is assumed to be a fixed value for estimation, and the receiver clock error is regarded as white noise and its mean is taken as 0 to reduce the influence of the clock error on the positioning accuracy; in the dynamic mode, the station position is estimated through a random walk process, and the position estimate is updated in real time to meet the positioning needs in a dynamic environment.

Images