A high-precision precise point positioning method based on BDS-3

Abstract

The application provides a high-precision precise point positioning method based on BDS-3, and relates to the technical field of wireless communication. Firstly, a traditional ionosphere-free combination model is introduced to eliminate the influence of ionosphere delay on positioning accuracy, satellite orbit and clock error data provided by IGS are used to eliminate or weaken satellite orbit error and clock error; secondly, a self-adaptive weighted wet delay troposphere correction model is constructed to eliminate troposphere delay error; then, a self-adaptive random model is established to reasonably weight the pseudo-range and carrier phase observation; finally, the estimated parameters are estimated through extended Kalman filtering, so that the positioning solution is obtained. The application has obvious advantages in improving positioning accuracy and can meet the demand of high-precision positioning. Whether in dynamic mode or static mode, the positioning accuracy can reach centimeter level, and in some stations and directions, the positioning accuracy in static mode can even reach millimeter level.

Description

Technical Field

[0001] This invention relates to the field of wireless communication technology, and in particular to a high-precision single-point positioning method based on BDS-3. Background Technology

[0002] The BeiDou Navigation Satellite System (BDS), as one of the world's leading satellite navigation systems, is widely used in transportation, surveying and mapping, agriculture, disaster monitoring, and other fields due to its high-precision positioning capabilities. Precise Point Positioning (PPP) technology is one of the core positioning methods in the BeiDou system, offering advantages such as eliminating the need for a reference station, reducing labor intensity and production costs. Using a single receiver, PPP technology can achieve high-precision static and dynamic positioning globally, solving the problem of traditional relative positioning methods being unusable at long distances or in the absence of reference stations. However, because PPP technology does not rely on reference stations for error correction, its positioning accuracy is affected by various error sources, especially interference from satellite clock errors, orbital errors, and atmospheric delays. To improve positioning accuracy, an effective error correction model must be established to reduce the impact of these errors on the positioning results. Furthermore, PPP positioning calculations rely on two types of observations: pseudorange and carrier phase. Its accuracy and stability are directly affected by the weighting of these observations. Reasonably allocating the weights of the observations and effectively eliminating various errors will directly improve its positioning accuracy and reliability. Therefore, in-depth research on PPP technology in the BeiDou system and optimization of the error correction model are of great theoretical and practical significance 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 single-point positioning method based on BDS-3, which addresses the shortcomings of the prior art.

[0004] To solve the above-mentioned technical problems, the technical solution adopted by the present invention is as follows: A high-precision single-point positioning method based on BDS-3 is proposed. First, a traditional ionospheric-free combined model is introduced to eliminate the impact of ionospheric delay on positioning accuracy. Satellite orbit and clock error data provided by IGS are used to eliminate or reduce satellite orbit and clock error. Second, an adaptive weighted wet delay tropospheric correction model is constructed to eliminate tropospheric delay error. Next, an adaptive stochastic model is established to reasonably assign weights to pseudorange and carrier phase observations. Finally, the parameters to be estimated are estimated by extended Kalman filtering to obtain the positioning solution.

[0005] Furthermore, in the ionospheric-delay-free combination model, the observations of the two frequencies are weighted and combined by constructing a phase combination and a pseudorange combination without ionospheric delay, thereby canceling the interference of ionospheric delay; the phase and pseudorange observations of the ionospheric-delay-free combination model are as follows: (1) (2) in, and These are pseudorange and carrier phase observations without ionospheric assemblies, respectively; superscript Indicates satellite, subscript Indicates the location of the monitoring station; the subscript IF indicates a combination without an ionosphere. This represents the geometric distance between the satellite's position and the station's position. It is the speed of light in a vacuum; It is the receiver clock bias; It is satellite clock bias; It is the tropospheric delay of the signal propagation path; For the station satellite The corresponding noise for the combination of ionospherically unsaturated delayed pseudo-ranges; For the station satellite The corresponding noise for the ionosphere-free delayed phase combination; For phase ambiguity.

[0006] Furthermore, the adaptive weighted wet delay tropospheric correction model is as follows: Tropospheric delay in signal propagation path Represented as: (3) in, and These are the dry tropospheric delay and the wet tropospheric delay, respectively. and The projection functions corresponding to the dry delay component and the wet delay component are given in the following two equations: (4) (5) in, It is the vertical distance between the stations; The elevation angle is the direction of the line connecting the station and the satellite. , , These are constants determined by the NMF model; , , Obtained by interpolation of NMF dryness coefficients; , , Obtained by interpolation of the NMF moisture component coefficient; a northward horizontal gradient incorporating moisture delay is introduced. and horizontal gradient in the east direction Based on this, the weights of the wet delay gradient are further dynamically adjusted, and the projection function of the wet delay is... Further correction to: (6) in, It is the azimuth angle of the line connecting the station and the satellite; The humidity weighting function is expressed as follows: (7) in, and It is an adjustment parameter used to flexibly control the range of change in humidity weight; Indicates humidity; The humidity weighting function In the middle, when the humidity is high, that is When the percentage is >80%, increase the weight of the wet delay gradient and adjust the factor. Set to 0.05~0.15: If the humidity is high but not extremely humid, i.e., around 80%, the adjustment factor should be set to 0.05; if the humidity is higher than 80%, the moisture delay has a greater impact on signal propagation, requiring a stronger correction to compensate for the delay changes caused by humidity, thus using a higher adjustment factor. ; When the humidity is low, that is When the gradient is less than 30%, reduce the weight of the wet delay gradient and adjust the factor. Set to 0.05~0.1; when the humidity is very low, i.e. below 20%, the effect of humidity delay on signal propagation is minimal, so choose a smaller adjustment factor. Set it to 0.05; if the humidity is low, i.e., between 20% and 30%, the effect of the wet delay is small, but it still exists. In this case, appropriately increase the adjustment factor. ; Under moderate humidity conditions, i.e., 30 ≤ ≤80%, keep the weight at the baseline value of 1.

[0007] After combining the above components, the final tropospheric delay along the signal propagation path is... Calculated using the following formula: (8) in, Indicates tropospheric delay in the zenith direction; Obtained from the projection function model NMF, The dry component delay in the zenith direction The tropospheric delay in the zenith direction was calculated using the Saastamoinen model based on standard meteorological parameters. and the northward horizontal gradient of the wet delay and horizontal gradient in the east direction They are treated 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 expressions of which are: (9) in, and The vector represents the difference between the observed and calculated phase and pseudorange values; A is the design matrix. The parameters to be estimated include the location of the station. Distance error corresponding to receiver clock bias Stratotropic delay at the zenith and horizontal gradient , ; For phase ambiguity; It is the identity matrix; and These are the residuals of phase and pseudorange observations; After solving the observation equation (9), the parameters to be estimated for joint positioning using pseudorange and phase observations are obtained as follows: (10) in, ; ; , and Let X and Y represent the estimated values, respectively. This is the corresponding weight matrix for pseudorange and phase observations; and Determined based on the PPP stochastic model.

[0008] Furthermore, the establishment of the adaptive stochastic model first involves evaluating the observation accuracy of GEO, IGSO, and MEO satellites using residual values, and then constructing an adaptive stochastic model by combining the satellite elevation angle model; the specific method is as follows: In the ionosphere-free combined model, subtracting the ionosphere-free phase and pseudorange combination yields: (11) By combining the standard deviations of phase and pseudorange observations, the delayed standard deviation of the combined difference between ionospheric phase and pseudorange is calculated to quantify the magnitude and distribution of residual errors in the difference; given the known phase standard deviation... and pseudorange standard deviation In this case, by using the error propagation rule, the errors of phase and pseudorange are combined, and the root mean square error corresponding to their difference is calculated as follows: (12) in, For the station satellite The root mean square error corresponding to the combined difference of phase and pseudorange without ionosphere; exist When no cycle slip occurs in any of the observation epochs, the ionospheric-free combination can effectively reduce the impact of ionospheric errors. In this case, 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 this difference is described by the following expression: (13) in, This is the residual corresponding to the combined difference between the ionospheric delay phase and pseudorange; The overall accuracy of satellite observation data is quantified by calculating the root mean square error (RMSE). For residuals related to ionospheric combination differences, the RMSE is used to assess the magnitude of errors in phase and pseudorange differences. For ionospheric combination differences, the RMSE of the residuals is calculated using the following expression: (14) The mean square root error of satellite GEO, IGSO, and MEO was calculated separately, and the accuracy of each satellite data type was measured based on these errors, as shown in the following formula: (15) in, This indicates the type of BeiDou satellite, namely GEO, IGSO, or MEO; This represents the average root mean square error value for one of the three types of BeiDou satellites. Representing the The number of observations of satellite-like objects; According to equation (15), if the root mean square error of MEO is the largest, that is, the observation accuracy of MEO satellite is the lowest, then the variance of the stochastic model is expressed as: (16) in, , ; The solution is based on the elevation angle, as shown in the following formula: (17) in, and All are set to 0.003m; Indicates the elevation angle; If, according to equation (15), the observation accuracy of the IGSO satellite is the lowest, then in the above formula (16), =IGSO ; =GEO, MEO , , ; If, according to equation (15), the observation accuracy of GEO satellites is the lowest, then in the above formula (16), When =GEO, ; =IGSO, MEO , , ; By analyzing the accuracy of observations from different types of satellites (GEO, IGSO, and MEO), different weights are assigned to each type of satellite to optimize the stochastic model in the positioning calculation process. The adaptive stochastic model, composed of the observation variance factor and variance, is as follows: (18) in, To indicate a stochastic model, its subscript , These are the phase observations and pseudorange observations, respectively; variance factor. This is the weighting ratio of phase observations to pseudorange observations. = , , These are the residuals of the carrier phase observations. and pseudorange observation residuals The root mean square error value; Therefore, the corresponding weight matrices for pseudorange and phase observations in equation (10) are: (19).

[0009] Furthermore, the specific method for estimating the parameters to be estimated using extended Kalman filtering is as follows: After obtaining the phase pseudorange observation weight matrix through equation (19), in PPP mode, an extended Kalman filter is introduced to estimate the parameters to be estimated in equation (10). From equation (10), the vector to be estimated is: (20) The extended Kalman filter estimation requires the following observation data as input: phase observations. and pseudo-distance observation These data can provide the relationship between the system state and the observations, thus determining the expression for the observation vector as follows: (twenty one) (twenty two) (twenty three) In the extended Kalman filter estimation process, the state estimate is continuously corrected and improved based on the current observation information at each epoch by updating the state vector and covariance matrix to be estimated; The state vector to be estimated at each epoch and the corresponding covariance matrix Able to utilize observation vectors at the same time The estimation is performed, and the model is as follows: (twenty four) (25) (26) In the formula, For observation data in The covariance matrix updated at each epoch; For observation data in The covariance matrix before the epoch update; for The Kalman gain matrix at each epoch controls the balance between prediction and measurement; This is the observation matrix, representing the relationship between state variables and observed values; It is a partial derivative matrix; for The covariance matrix of measurement errors at epoch times; and They are respectively The state vector before and after the update at each epoch; It is the partial derivative matrix of the observation model with respect to the state vector, which describes the impact of changes in the state vector on the observed values ​​and reflects the linearized relationship between the state and the observation. The form is: (27) Each term represents the partial derivative of the observed value with respect to the state variable, as detailed below: (28) In the formula, The unit vector from the station to the satellite. The superscript T indicates transpose. (29) (30) (31) in, It is the projection function of the tropospheric wet delay in the zenith direction.

[0010] Furthermore, the extended Kalman filter method updates the state vector and covariance matrix at each stage, thereby estimating the current state of the system more accurately; specifically, the formulas for updating the state vector and its associated covariance matrix are described as follows: (32) (33) in, From The epochal time has arrived The system noise transfer matrix at each epoch; for The epochal time has arrived The system noise covariance matrix at each epoch.

[0011] Furthermore, in the BDS-3-based high-precision single-point positioning method, the station position is assumed to be a fixed value for estimation in static mode, and the receiver clock error is treated as white noise with a mean of 0 to reduce the impact of clock error on positioning accuracy; while in dynamic mode, the station position is estimated through a random walk process and the position estimate is updated in real time to adapt to the positioning needs in dynamic environments.

[0012] The beneficial effects of adopting the above technical solution are as follows: The high-precision single-point positioning method based on BDS-3 provided by the present invention has obvious advantages in improving positioning accuracy and can meet the needs of high-precision positioning. It can achieve centimeter-level positioning accuracy in both dynamic and static modes, and the positioning accuracy in static mode can even reach millimeter level in some stations and directions. Attached Figure Description

[0013] Figure 1 A schematic diagram of the high-precision single-point positioning method based on BDS-3 provided in this embodiment of the invention; Figure 2 The positioning error statistics of three precise single-point positioning methods for stations URUM, SGOC, CAS1 and WUH2 in the E direction in static mode provided for embodiments of the present invention; Figure 3The positioning error statistics of three precise single-point positioning methods for stations URUM, SGOC, CAS1 and WUH2 in the N direction in static mode provided for embodiments of the present invention; Figure 4 The positioning error statistics of three precise single-point positioning methods for stations URUM, SGOC, CAS1 and WUH2 in the U direction in static mode provided for embodiments of the present invention; Figure 5 Error comparison diagram of three dynamic precise single-point positioning methods for URUM stations provided in the embodiments of the present invention; Figure 6 Error comparison diagram of three dynamic precise single-point positioning methods for SGOC stations provided in the embodiments of the present invention; Figure 7 Error comparison diagram of three dynamic precise single-point positioning methods for WUH2 station provided in the embodiments of the present invention; Figure 8 Error comparison chart of three dynamic precise single-point positioning methods for CAS1 station provided in the embodiments of the present invention. Detailed Implementation

[0014] The specific embodiments of the present invention will be described in further detail below with reference to the accompanying drawings and examples. The following examples are for illustrative purposes only and are not intended to limit the scope of the invention.

[0015] This embodiment provides a high-precision precise point positioning method based on BDS-3. Addressing the issues of uneven tropospheric wet delay distribution and humidity differences under varying conditions, this method constructs an adaptive weighted wet delay tropospheric correction model. Based on the ZTD parameter estimation model, this model introduces northward and eastward horizontal gradients of the wet delay to correct the wet delay projection function, and constructs a humidity weight function to further dynamically adjust the wet delay projection function. Furthermore, considering that traditional stochastic models fail to fully reflect the differences in observation accuracy among different types of BDS-3 satellites, this method establishes an adaptive stochastic model 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. Figure 1 As shown, the method of this embodiment is described below.

[0016] Step 1: The high-precision single-point positioning method based on BDS-3 introduces a traditional ionospheric-delayed combination model. This model can effectively eliminate the impact of ionospheric delay on positioning accuracy during signal propagation. In this model, by constructing ionospheric-delay-free phase combination and pseudorange combination, the observations of the two frequencies are weighted and combined to cancel out the interference of ionospheric delay. The observed values ​​of phase and pseudorange in the ionospheric-delayed combination model are: (1) (2) In the formula, and These are pseudorange and carrier phase observations without ionospheric assemblies, respectively; superscript Indicates satellite, subscript This indicates the coordinates of the station; the subscript IF indicates a combination without an ionosphere. This represents the geometric distance between the satellite's position and the station's position. It is the speed of light in a vacuum; It is the receiver clock bias; It is satellite clock bias; It is the tropospheric delay of the signal propagation path; For the station satellite The corresponding noise for the combination of ionospherically unsaturated delayed pseudo-ranges; For the station satellite The corresponding noise for the ionosphere-free delayed phase combination; For phase ambiguity.

[0017] As can be seen from equations (1) and (2), although the traditional ionosphere-free model eliminates the influence of the ionosphere, satellite clock bias still exists. Receiver clock bias Observation noise and and tropospheric delay error To mitigate the impact of tropospheric delay errors, satellite orbit and clock bias data provided by IGS are used in the data processing for precise point positioning. The impact of this has become a key issue.

[0018] Step 2: Considering tropospheric delay error To investigate the uneven distribution of tropospheric wet delay and the impact of humidity differences under different climatic conditions, an adaptive weighted tropospheric correction model for wet delay is constructed based on the ZTD parameter estimation model.

[0019] Tropospheric delay in signal propagation path Represented as: (3) In the formula, and These are the dry tropospheric delay and the wet tropospheric delay, respectively. and These are the projection functions corresponding to the dry and wet delay components, as shown in the following two equations: (4) (5) in, It is the vertical distance between the stations; The elevation angle is the direction of the line connecting the station and the satellite. , , These are constants determined by the NMF model; , , Obtained by interpolation of NMF dryness coefficients; , , Obtained by interpolation of the NMF wet component coefficient.

[0020] To correct the uneven distribution of tropospheric wet delay in different directions and its distribution under different climatic conditions, a northward horizontal gradient of the wet delay is introduced. and horizontal gradient in the east direction Based on this, the weights of the wet delay gradient are further dynamically adjusted. For this, the projection function of the wet delay... Further correction to: (6) in, It is the azimuth angle of the line connecting the station and the satellite; The humidity weighting function is expressed as follows: (7) In the formula, and This is an adjustment parameter used to flexibly control the magnitude of changes in humidity weighting. When humidity... higher ( When the gradient is >80%, increase the weight of the wet delay gradient and adjust the factor. Set to 0.05~0.15; if the humidity is high but not extremely humid (around 80%), adjust the factor accordingly. Setting it to 0.05 is necessary when humidity exceeds 80%, as humidity delay significantly impacts signal propagation. Therefore, a stronger correction is needed to compensate for the delay changes caused by humidity, requiring a higher adjustment factor. When humidity Lower ( When the percentage is less than 30%, reduce the weight of the wet delay gradient and adjust the factor. Set to 0.05~0.1; when humidity is very low (e.g., below 20%), the effect of humidity delay on signal propagation is minimal, and over-correction may introduce unnecessary errors. Therefore, a smaller adjustment factor should be chosen. Set it to 0.05; if the humidity is low (between 20% and 30%), the effect of the wet delay is small, but it still exists. In this case, the adjustment factor can be increased appropriately. To better correct for errors caused by humidity. Under moderate humidity conditions (30 ≤ (≤80%), keep the weight at the baseline value of 1.

[0021] After combining the above components, the final tropospheric delay along the signal propagation path is... It can be calculated using the following formula: (8) in, This indicates tropospheric delay in the zenith direction.

[0022] In the process of high-precision single-point positioning based on BDS-3 Obtained from the projection function model NMF, The dry component delay in the zenith direction Calculated by the Saastamoinen model based on standard meteorological parameters. Tropospheric delay in the zenith direction. and the northward horizontal gradient of the wet delay and horizontal gradient in the east direction It is estimated as an unknown parameter in order to more accurately correct the impact of tropospheric delay on the positioning results.

[0023] After processing the tropospheric delay error, equations (1) and (2) are linearized to simplify the solution process using the observation equations. The simplified observation equations are obtained through linearization, and their specific expressions are as follows: (9) In the formula, and The vector represents the difference between the observed and calculated phase and pseudorange values; A is the design matrix. Parameters to be estimated: including station location Distance error corresponding to receiver clock bias Stratotropic delay at the zenith and horizontal gradient , ; For phase ambiguity; It is the identity matrix; and It is the residual of phase and pseudorange observations.

[0024] After solving the observation equation (9), the parameters to be estimated for joint positioning using pseudorange and phase observations can be obtained: (10) In the formula, ; ; , and Let X and Y represent the estimated values, respectively. This is the weight matrix for the pseudorange and phase observations. and Determined based on the PPP stochastic model.

[0025] Step 3: As can be seen from equation (10), the key to solving the parameter to be estimated lies in accurately calculating the corresponding weights of the pseudorange observations. The corresponding weights of phase observations Therefore, establishing a reasonable stochastic model is crucial. Considering the limitations of traditional stochastic models in failing to fully reflect the differences in observation accuracy among different satellites, an adaptive stochastic model is proposed based on the accuracy evaluation results of BDS-3 observation data to weight pseudorange and carrier phase observations. This method first evaluates the observation accuracy of GEO, IGSO, and MEO satellites using residual values, and then constructs an adaptive stochastic model by combining the satellite elevation angle model.

[0026] In the ionospheric-free combined model, the influence of ionospheric delay can be effectively eliminated by subtracting the phase and pseudorange combinations, thus obtaining a relatively pure expression for the observation difference, which helps to more accurately assess various observation errors. Subtracting the ionospheric-free phase and pseudorange combinations yields: (11) Furthermore, by combining the standard deviations of phase and pseudorange observations, the delayed standard deviation of the combined difference between ionospheric phase and pseudorange can be calculated to quantify the magnitude and distribution of residual errors in the difference. Given the phase standard deviation... and pseudorange standard deviation In this case, by using the error propagation rule, the errors of phase and pseudorange are combined, and the root mean square error corresponding to their difference can be calculated as follows: (12) exist In the absence of cycle slips in any observation epoch, the ionospheric-free combination can effectively reduce the impact of ionospheric errors. In this case, the difference between the phase and pseudorange observations (i.e., the residual) is typically caused by system noise and other error sources. This residual, related to the difference, can be described by the following expression: (13) in, This represents the residual corresponding to the combined difference between the ionospheric delay phase and pseudorange.

[0027] The overall accuracy of satellite observation data is quantified by calculating the root mean square error (RMSE). For residuals related to ionospheric combination differences, the RMS error can be used to assess the magnitude of errors in phase and pseudorange differences. For ionospheric combination differences, the RMS error of the residuals can be calculated using the following expression: (14) Next, the satellites were calculated separately. , and The average root mean square error, and based on these errors, the accuracy of each type of satellite data is measured: (15) In the formula, Indicates the type of BeiDou satellite, i.e. , or ; This represents the average root mean square error value for one of the three types of BeiDou satellites. Representing the The number of observations by satellite-like objects.

[0028] According to equation (15), if we obtain The root mean square error value is the largest, that is... Satellites have the lowest observation accuracy, therefore the variance of the stochastic model can be expressed as: (16) in, ; The solution is based on the elevation angle, as shown in the following formula: (17) in, and All are set to 0.003m; Indicates the altitude angle.

[0029] If, according to equation (15), the observation accuracy of the IGSO satellite is the lowest, then in the above formula (16), =IGSO ; =GEO, MEO , , .

[0030] If, according to equation (15), the observation accuracy of GEO satellites is the lowest, then in the above formula (16), When =GEO, ; =IGSO, MEO , , .

[0031] Through the , and The analysis of the accuracy of observations from different types of satellites assigns different weights to each type to optimize the stochastic model in the positioning calculation process. The adaptive stochastic model, combining the observation variance factor and variance, is as follows: (18) in, To indicate a stochastic model, its subscript , These are the phase observations and pseudorange observations, respectively; variance factor. This is the weighting ratio of phase observations to pseudorange observations. = , , These are the residuals of the carrier phase observations. and pseudorange observation residuals The root mean square error value.

[0032] Therefore, the corresponding weight matrices for pseudorange and phase observations in equation (10) are: (19).

[0033] Step 4: After obtaining the phase pseudorange observation weight matrix through equation (19), in PPP mode, an extended Kalman filter (EKF) is introduced to estimate the parameters to be estimated in equation (10). From equation (10), the vector to be estimated is: (20) The EKF filter estimation requires the following observation data as input: phase observations. and pseudo-distance observation These data can provide the relationship between the system state and the observations, thus determining the expression for the observation vector as follows: (twenty one) (twenty two) (twenty three) In the EKF filtering estimation process, by updating the state vector and covariance matrix to be estimated, the state estimate can be continuously corrected and improved based on the current observation information at each epoch. The state vector to be estimated at each epoch and the corresponding covariance matrix Able to utilize observation vectors at the same time The estimation is performed, and the model is as follows: (twenty four) (25) (26) In the formula, For observation data in The covariance matrix updated at each epoch; For observation data in The covariance matrix before the epoch update; for The Kalman gain matrix at each epoch controls the balance between prediction and measurement; This is the observation matrix, representing the relationship between state variables and observed values; It is a partial derivative matrix; for The covariance matrix of measurement errors at epoch times; and They are respectively The state vector before and after the update at each epoch.

[0034] It is the partial derivative of the observation model with respect to the state vector. It describes the impact of changes in the state vector on the observed values ​​and reflects the linearized relationship between the state and the observation. Specifically, The form is: (27) Partial derivative matrix In the formula, each term represents the partial derivative of the observed value with respect to the state variable, as follows: (28) In the formula, ( ) is the unit vector in the direction of the satellite at the station.

[0035] (29) (30) (31) In the formula, It is the projection function of the tropospheric wet delay in the zenith direction, and the meanings of the other symbols are the same as in formula (6).

[0036] The calculation and determination of each term in H(x) depend on the state vector to be estimated. The parameters in the M. For example, M T Elements in the matrix The projection function of the zenith-direction tropospheric wet delay is related to parameters such as the zenith-direction tropospheric delay and the tropospheric horizontal gradient. This is because the calculation of the tropospheric delay and the determination of the projection function involve these tropospheric parameters, which are contained within the state vector to be estimated. middle.

[0037] The EKF filtering method updates the state vector and covariance matrix at each stage, thereby estimating the current state of the system more accurately. Specifically, the formulas for updating the state vector and its associated covariance matrix are described as follows: (32) (33) In the formula, From The epochal time has arrived The system noise transfer matrix at time t; for Time's up The system noise covariance matrix at time t.

[0038] In high-precision point positioning based on BDS-3, in static mode, the station position is typically estimated as a fixed value, and the receiver clock error is treated as white noise with a mean of 0 to reduce the impact of clock error on positioning accuracy. In dynamic mode, the station position is estimated through a random walk process, with the position estimate updated in real time to adapt to positioning requirements in dynamic environments.

[0039] To verify the effectiveness and superiority of the high-precision point positioning method based on BDS-3 in this embodiment, a series of positioning simulation experiments were conducted to evaluate the technology's ability to accurately calculate station position coordinates. Multiple stations in the MGEX network were selected for the experiments, with an observation period of 24 hours from December 1, 2023, and an epoch interval of 30 seconds, containing a total of 2880 epochs of data. To achieve high-precision PPP positioning, this experiment used 30-second precise ephemeris and satellite clock bias products provided by a university's IGS data center and utilized the BDS-3 system for positioning calculation. The specific configuration was as follows: the satellite elevation angle was set to 15°, the BDS-3 B1 / B3 ionospheric elimination combination observations were selected, and the receiver sampling interval was set to 30 seconds.

[0040] To evaluate the positioning performance of the high-precision point positioning method based on BDS-3, it is compared and analyzed with two widely used positioning methods: one is the WH-ER method, which uses a tropospheric delay model based on parameter estimation and weighting according to the experience ratio (WH-ER); the other is the WH-RPO method, which uses a tropospheric delay model based on parameter estimation and weighting according to the residual ratio of phase and pseudo-range observations (WH-RPO). Simulation experiments in both dynamic and static modes were conducted on the BDS-3-based high-precision point positioning method, the WH-ER method, and the WH-RPO method. By comparing the positioning accuracy of different methods, the positioning performance of each method is evaluated.

[0041] The following section analyzes the positioning accuracy of the BDS-3-based high-precision positioning method in static mode.

[0042] Figure 2 , Figure 3 and Figure 4 The diagram shows the positioning error statistics in the E, N, and U directions for three precise single-point positioning methods at the URUM, SGOC, CAS1, and WUH2 stations in static mode. Figure 2 As can be clearly seen, the high-precision point positioning method based on BDS-3 achieves centimeter-level positioning errors in the E, N, and U directions, meeting the required positioning standards and verifying the method's effectiveness. Furthermore, compared to the WH-ER and WH-RPO precise point positioning methods, this method exhibits a significant advantage in positioning errors across all directions. Although there is a slight increase in error at some stations and in specific directions (such as the E direction at the URUM station), exceeding that of the WH-RPO method, overall, the high-precision point positioning method based on BDS-3 maintains the lowest positioning error level in all directions. Therefore, the BDS-3-based precise point positioning method has a significant advantage in positioning accuracy compared to traditional methods, providing higher-precision positioning results and meeting the requirements of high-precision positioning.

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

[0044] Table 1. Positioning accuracy improvement of the high-precision single-point positioning method compared to the precision single-point positioning WH-ER method.

[0045] As shown in Table 1, the high-precision point positioning method based on BDS-3 demonstrates improved positioning accuracy in the E, N, and U directions compared to the WH-ER method for precise point positioning. Particularly at the SGOC station, the improvement in positioning accuracy in the E, N, and U directions is significant, with an average improvement of 3.40 cm. The average improvements at other stations are 3.17 cm (URUM), 2.46 cm (CAS1), and 2.82 cm (WUH2), respectively. These results further validate the advantage of the high-precision point positioning method over the WH-ER method in improving positioning accuracy. Table 2 shows the average improvement in positioning accuracy of the high-precision point positioning method compared to the WH-RPO method.

[0046] As shown in Table 2, the high-precision point positioning method based on BDS-3 improves the positioning accuracy to varying degrees compared to the WH-RPO method at each station. 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, with an overall average improvement of 2.01 cm. Other stations, such as SGOC, CAS1, and WUH2, also show varying degrees of improvement, with SGOC showing the largest average improvement of 2.38 cm. These results indicate that the high-precision point positioning method in this embodiment has a relatively consistent accuracy advantage compared to the WH-RPO method.

[0047] The following section analyzes the positioning accuracy of the BDS-3-based high-precision positioning method in dynamic mode.

[0048] Figure 5 , Figure 6 , Figure 7 as well as Figure 8This paper presents a comparison of errors for three dynamic precise point positioning methods at the URUM, SGOC, WUH2, and CAS1 stations. The figures clearly show that the positioning error fluctuation of the high-precision point positioning method based on BDS-3 is significantly smaller than that of the other two methods. Particularly at the initial epochs of the URUM and SGOC stations, the initial epoch of the WUH2 station, and around epoch 1600, the positioning error fluctuation of the high-precision point positioning method based on BDS-3 is smaller, exhibiting more stable characteristics. Furthermore, at the CAS1 station, the positioning error fluctuation range of the high-precision point positioning method based on BDS-3 is significantly lower than that of the precise point positioning methods WH-ER and WH-RPO, further validating the superiority of this method in dynamic positioning.

[0049] Table 3 shows the positioning accuracy of three dynamic high-precision point positioning methods in different directions. Statistical analysis leads to the following conclusions: At the URUM station, the dynamic high-precision point positioning method improved the average three-dimensional positioning error by 3.70 cm and 2.42 cm compared to the WH-RPO and WH-ER methods, respectively, demonstrating the significant advantage of the dynamic high-precision method in positioning accuracy. For the SGOC station, the improvement was even more significant, at 3.68 cm and 2.31 cm, further validating the method's superiority in improving positioning accuracy. Meanwhile, at the WUH2 station, the dynamic high-precision point positioning method improved the accuracy by 4.99 cm and 1.58 cm compared to the WH-RPO and WH-ER methods, again showcasing its outstanding performance in improving positioning accuracy. Based on the above results, it can be concluded that the dynamic high-precision single-point positioning method has higher positioning accuracy than the WH-RPO and WH-ER methods, demonstrating its significant advantage in accuracy improvement.

[0050] Table 3 Positioning accuracy of dynamic precision single-point positioning

[0051] In summary, the high-precision point positioning method based on BDS-3 achieves centimeter-level positioning accuracy in both dynamic and static modes, and even reaches millimeter-level accuracy in static mode for certain stations and directions. This result fully verifies the effectiveness of the method. Furthermore, compared to the WH-RPO and WH-ER precise point positioning methods, the high-precision point positioning method exhibits superior 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 point positioning method based on BDS-3 has significant advantages in improving positioning accuracy and can meet the needs of high-precision positioning.

[0052] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, and not to limit them; although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some or all of the technical features therein; and these modifications or substitutions 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 single-point positioning method based on BDS-3, characterized in that: First, a traditional ionospheric-free combined model is introduced to eliminate the impact of ionospheric delay on positioning accuracy. Satellite orbit and clock error data provided by IGS are used to eliminate or reduce satellite orbit and clock error. Second, an adaptive weighted wet delay tropospheric correction model is constructed to eliminate tropospheric delay error. Next, an adaptive stochastic model is established to reasonably assign weights to pseudorange and carrier phase observations. Finally, the parameters to be estimated are estimated by extended Kalman filtering to obtain the positioning solution. In the ionospheric-delay-free combined model, the observations of the two frequencies are weighted and combined by constructing a phase combination and a pseudorange combination without ionospheric delay, thereby canceling the interference of ionospheric delay. The phase and pseudorange observations of the ionospheric-delay-free combined model are as follows: （1） （2） in, and These are pseudorange and carrier phase observations without ionospheric assemblies, respectively; superscript Indicates satellite, subscript Indicates the location of the monitoring station; the subscript IF indicates a combination without an ionosphere. This represents the geometric distance between the satellite's position and the station's position. It is the speed of light in a vacuum; It is the receiver clock bias; It is satellite clock bias; It is the tropospheric delay of the signal propagation path; For the station satellite The corresponding noise for the combination of ionospherically unsaturated delayed pseudo-ranges; For the station satellite The corresponding noise for the ionosphere-free delayed phase combination; For phase ambiguity; The adaptive weighted wet delay tropospheric correction model is as follows: Tropospheric delay in signal propagation path Represented as: （3） in, and These are the dry tropospheric delay and the wet tropospheric delay, respectively. and The projection functions corresponding to the dry delay component and the wet delay component are given in the following two equations: （4） （5） in, It is the vertical distance between the stations; The elevation angle is the direction of the line connecting the station and the satellite. , , These are constants determined by the NMF model; , , Obtained by interpolation of NMF dryness coefficients; , , Obtained by interpolation of the NMF moisture component coefficient; a northward horizontal gradient incorporating moisture delay is introduced. and horizontal gradient in the east direction Based on this, the weights of the wet delay gradient are further dynamically adjusted, and the projection function of the wet delay is... Further correction to: （6） in, It is the azimuth angle of the line connecting the station and the satellite; The humidity weighting function is expressed as follows: （7） in, and It is an adjustment parameter used to flexibly control the range of change in humidity weight; Indicates humidity; The humidity weighting function In the middle, when the humidity is high, that is When the percentage is >80%, increase the weight of the wet delay gradient and adjust the factor. Set to 0.05~0.15: If the humidity is high but not extremely humid, i.e., around 80%, the adjustment factor should be set to 0.05; if the humidity is higher than 80%, the moisture delay has a greater impact on signal propagation, requiring a stronger correction to compensate for the delay changes caused by humidity, thus using a higher adjustment factor. ; When the humidity is low, that is When the gradient is less than 30%, reduce the weight of the wet delay gradient and adjust the factor. Set to 0.05~0.1; when the humidity is very low, i.e. below 20%, the effect of humidity delay on signal propagation is minimal, so choose a smaller adjustment factor. Set it to 0.05; if the humidity is low, i.e., between 20% and 30%, the effect of the wet delay is small, but it still exists. In this case, appropriately increase the adjustment factor. ; Under moderate humidity conditions, i.e., 30 ≤ ≤80%, keep the weight at the baseline value of 1; After combining all the above components, the final tropospheric delay along the signal propagation path is... Calculated using the following formula: （8） in, Indicates tropospheric delay in the zenith direction; Obtained from the projection function model NMF, The dry component delay in the zenith direction The tropospheric delay in the zenith direction was calculated using the Saastamoinen model based on standard meteorological parameters. and the northward horizontal gradient of the wet delay and horizontal gradient in the east direction They are treated 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 expressions of which are: （9） in, and The vector represents the difference between the observed and calculated phase and pseudorange values; A is the design matrix. The parameters to be estimated include the location of the station. Distance error corresponding to receiver clock bias Stratotropic delay at the zenith and horizontal gradient , ; For phase ambiguity; It is the identity matrix; and These are the residuals of phase and pseudorange observations; After solving the observation equation (9), the parameters to be estimated for joint positioning using pseudorange and phase observations are obtained as follows: （10） in, ; ; , and Let X and Y represent the estimated values, respectively. This is the corresponding weight matrix for pseudorange and phase observations; and Determined based on the PPP stochastic model; The adaptive stochastic model is established by first evaluating the observation accuracy of GEO, IGSO, and MEO satellites using residual values, and then constructing an adaptive stochastic model by combining the satellite elevation angle model; the specific method is as follows: In the ionosphere-free combined model, subtracting the ionosphere-free phase and pseudorange combination yields: （11） By combining the standard deviations of phase and pseudorange observations, the delayed standard deviation of the combined difference between ionospheric phase and pseudorange is calculated to quantify the magnitude and distribution of residual errors in the difference; given the known phase standard deviation... and pseudorange standard deviation In this case, by using the error propagation rule, the errors of phase and pseudorange are combined, and the root mean square error corresponding to their difference is calculated as follows: （12） in, For the station satellite The root mean square error corresponding to the combined difference of phase and pseudorange without ionosphere; exist When no cycle slip occurs in any of the observation epochs, the ionospheric-free combination can effectively reduce the impact of ionospheric errors. In this case, 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 this difference is described by the following expression: （13） in, This is the residual corresponding to the combined difference between the ionospheric delay phase and pseudorange; The overall accuracy of satellite observation data is quantified by calculating the root mean square error (RMSE). For residuals related to ionospheric combination differences, the RMSE is used to assess the magnitude of errors in phase and pseudorange differences. For ionospheric combination differences, the RMSE of the residuals is calculated using the following expression: （14） The mean square root error of satellite GEO, IGSO, and MEO was calculated separately, and the accuracy of each satellite data type was measured based on these errors, as shown in the following formula: （15） in, This indicates the type of BeiDou satellite, namely GEO, IGSO, or MEO; This represents the average root mean square error value for one of the three types of BeiDou satellites. Representing the The number of observations of satellite-like objects; According to equation (15), if the root mean square error of MEO is the largest, that is, the observation accuracy of MEO satellite is the lowest, then the variance of the stochastic model is expressed as: （16） in, , ; The solution is based on the elevation angle, as shown in the following formula: （17） in, and All are set to 0.003m; Indicates the elevation angle; If, according to equation (15), the observation accuracy of the IGSO satellite is the lowest, then in the above formula (16), =IGSO ; =GEO, MEO , , ; If, according to equation (15), the observation accuracy of GEO satellites is the lowest, then in the above formula (16), When =GEO, ; =IGSO, MEO , , ; By analyzing the accuracy of observations from different types of satellites (GEO, IGSO, and MEO), different weights are assigned to each type of satellite to optimize the stochastic model in the positioning calculation process. The adaptive stochastic model, composed of the observation variance factor and variance, is as follows: （18） in, To indicate a stochastic model, its subscript , These are the phase observations and pseudorange observations, respectively; variance factor. This is the weighting ratio of phase observations to pseudorange observations. = , , These are the residuals of the carrier phase observations. and pseudorange observation residuals The root mean square error value; Therefore, the corresponding weight matrices for pseudorange and phase observations in equation (10) are: （19）。 2. The high-precision single-point positioning method based on BDS-3 according to claim 1, characterized in that: The specific method for estimating the parameters to be estimated using extended Kalman filtering is as follows: After obtaining the phase pseudorange observation weight matrix through equation (19), in PPP mode, an extended Kalman filter is introduced to estimate the parameters to be estimated in equation (10). From equation (10), the vector to be estimated is: （20） The extended Kalman filter estimation requires the following observation data as input: phase observations. and pseudo-distance observation These data can provide the relationship between the system state and the observations, thus determining the expression for the observation vector as follows: （21） （22） （23） In the extended Kalman filter estimation process, the state estimate is continuously corrected and improved based on the current observation information at each epoch by updating the state vector and covariance matrix to be estimated; The state vector to be estimated at each epoch and the corresponding covariance matrix Able to utilize observation vectors at the same time The estimation is performed, and the model is as follows: （24） （25） （26） In the formula, For observation data in The covariance matrix updated at each epoch; For observation data in The covariance matrix before the epoch update; for The Kalman gain matrix at each epoch controls the balance between prediction and measurement; This is the observation matrix, representing the relationship between state variables and observed values; It is a partial derivative matrix; for The covariance matrix of measurement errors at epoch times; and They are respectively The state vector before and after the update at each epoch; It is the partial derivative matrix of the observation model with respect to the state vector, which describes the impact of changes in the state vector on the observed values ​​and reflects the linearized relationship between the state and the observation. The form is: （27） Each term represents the partial derivative of the observed value with respect to the state variable, as detailed below: （28） In the formula, The unit vector from the station to the satellite. The superscript T indicates transpose. （29） （30） （31） in, It is the projection function of the tropospheric wet delay in the zenith direction.

3. The high-precision single-point positioning method based on BDS-3 according to claim 2, characterized in that: The extended Kalman filter method updates the state vector and covariance matrix at each stage to more accurately estimate the current state of the system; specifically, the formulas for updating the state vector and its associated covariance matrix are described as follows: （32） （33） in, From The epochal time has arrived The system noise transfer matrix at each epoch; for The epochal time has arrived The system noise covariance matrix at each epoch.

4. The high-precision single-point positioning method based on BDS-3 according to claim 3, characterized in that: In the BDS-3-based high-precision single-point positioning method, the station position is estimated as a fixed value in static mode, and the receiver clock error is treated as white noise with a mean of 0 to reduce the impact of clock error on positioning accuracy. In dynamic mode, the station position is estimated through a random walk process, and the position estimate is updated in real time to adapt to the positioning needs in dynamic environments.

Images