A method for jointly adjusting observation weights of multi-system navigation satellite orbit determination networks

Abstract

This invention discloses a weighting method for adjustment observations in a multi-system navigation satellite joint orbit determination network, belonging to the field of satellite navigation and precise orbit determination data processing technology. The method includes: S1 observation data preprocessing and quality control; S2 orbit numerical integration and fitting; S3 prior weighted parameter estimation and residual editing; S4 grouping observations by station and navigation system dimensions, and initializing variance factors with equal weights; S5 estimating the variance scaling factor for each group using a simplified Helmert formula omitting the redundancy trace correction term and updating the weight matrix, combined with robust estimation to reduce the weight of outlier residuals, iterating until convergence; S6 re-estimating ambiguity with fixed constraints, outputting precise orbit and clock error products. This invention uses the above method, replacing fixed prior weights with data-driven simplified variance component estimation, to characterize the non-homogeneous characteristics of observations with almost no increase in computational burden, thus improving the accuracy and reliability of multi-system joint orbit determination.

Description

Technical Field

[0001] This invention relates to the field of satellite navigation and precise orbit determination data processing technology, and in particular to a method for weighting adjustment observations of a multi-system navigation satellite joint orbit determination network. Background Technology

[0002] Precise orbit determination for Global Navigation Satellite Systems (GNSS) is fundamental to satellite navigation applications. The accuracy of satellite orbits and clock errors directly determines the quality of user-end positioning, navigation, and timing services. With the comprehensive deployment of multiple systems such as GPS, GLONASS, Galileo, and BeiDou (BDS), the number of navigation satellites in orbit has now exceeded 120, and the global tracking station network has expanded to hundreds. Multi-system joint orbit determination has become the mainstream technical approach for generating precise orbit products by the International GNSS Service (IGS) and various analysis centers. Multi-system joint orbit determination is based on a weighted least squares estimation framework, incorporating observation data from multiple navigation systems into a unified normal equation for joint solution. In this framework, the observation weight matrix determines the contribution ratio of observations from different sources to parameter estimation, and its reasonable setting directly affects the solution accuracy and the statistical optimality of parameter estimation. Therefore, observation weighting is a key technical step in the adjustment of multi-system navigation satellite joint orbit determination networks.

[0003] However, the error characteristics of multi-system station network observation data have significant heterogeneity, including at least the following three levels of differences: (1) Differences in observation quality between systems. Different navigation systems have different orbit types, signal systems, and constellation maturity. For example, BDS includes three orbit types: MEO, GEO, and IGSO. GLONASS uses the FDMA system, while GPS and Galileo use the CDMA system. These differences lead to significant differences in the pseudorange and carrier phase observation noise levels of different systems. (2) Differences in observation quality between stations. Global tracking stations have different receiver models, antenna types, and station environments. Even for the same navigation system, the quality of observation data acquired by different stations varies significantly and changes with time and environmental conditions. (3) Differences in quality between observation types. The noise levels of pseudorange observations and carrier phase observations usually differ by more than two orders of magnitude, and the observation quality also varies at different frequencies.

[0004] To address the aforementioned heterogeneous error characteristics, existing multi-system navigation satellite joint orbit determination network adjustments typically employ a simplified strategy in observation weighting: assigning all stations and navigation systems the same prior variance or a fixed inter-system weighting ratio (e.g., GPS:GLONASS:Galileo:BDS = 1:1:1:1), and combining this with an elevation angle correlation function to provide relative weighting ratios. While this prior weighting strategy is simple to implement and computationally inexpensive, it has the following limitations: First, fixed prior variances or weighting ratios cannot adapt to the actual differences in observation quality under different station environments, making it difficult to achieve optimal estimation of parameters such as satellite orbits; second, inappropriate weighting ratios may negate the accuracy improvement brought about by multi-system fusion; and third, prior weighting cannot reflect the actual fluctuations in observation quality and lacks adaptive weight reduction capabilities for abnormal observation groups.

[0005] Variance Component Estimation (VCE) can estimate the variance components of each observation group by utilizing the statistical properties of posterior residuals, thereby achieving data-driven adaptive weighting. However, applying the complete Helmert variance component estimation to joint orbit determination network adjustment faces significant computational efficiency issues: the redundancy trace correction term in the standard formula requires calculating the trace of the product of the inverse of the normal equation matrix and the sub-normal equation matrices of each observation group. For joint orbit determination problems of multi-system navigation satellites with parameter scales reaching hundreds of thousands of dimensions, the calculation of this term requires storing a large-scale sub-normal equation matrix, resulting in extremely high memory requirements and computation time, which severely restricts the practical application of variance component estimation in operational joint orbit determination.

[0006] In summary, existing multi-system navigation satellite joint orbit determination network adjustments have not yet achieved station-by-station and system-by-system adaptive observation weighting while ensuring computational efficiency meets operational timeliness requirements. Therefore, there is an urgent need for an observation weighting method that balances accuracy and computational efficiency to meet the practical application needs of multi-system navigation satellite joint orbit determination. Summary of the Invention

[0007] The purpose of this invention is to provide a weighting method for adjustment observations in a multi-system navigation satellite joint orbit determination network. This method addresses the problems of existing multi-system navigation satellite joint orbit determination methods, which use fixed prior variances or empirical weight ratios that cannot adapt to the heterogeneous observation quality of different stations and navigation systems. Furthermore, the complete Helmert variance component estimation suffers from excessive computational costs due to the redundancy trace correction term, making it difficult to apply to operational joint orbit determination. This invention achieves data-driven adaptive observation weight estimation without significantly increasing computational burden, through a collaborative design of multi-granularity observation grouping on a station-by-station and system-by-system basis, simplified Helmert variance component estimation omitting the redundancy trace correction term, and an iterative weighted least squares solution framework combined with robust estimation. This improves the accuracy and reliability of multi-system navigation satellite joint orbit determination.

[0008] To achieve the above objectives, this invention provides a method for weighting adjustment observations in a multi-system navigation satellite joint orbit determination network, comprising the following steps: S1. Preprocessing and quality control of observation data: Acquire multi-system observation data and external model data from the global tracking station network, perform cycle slip detection, gross error removal and arc segment division on the multi-system observation data, and apply observation corrections. S2. Numerical integration and fitting of orbits: Numerical integration is performed on multi-system satellites based on the satellite dynamics model to obtain the reference orbit, state transition matrix and parameter sensitivity matrix, and orbit fitting and initial value correction are performed based on the external orbit. S3, Prior Equal Weight Parameter Estimation and Residual Editing: Construct normal equations with prior equal weights and solve them by least squares. Perform residual editing based on the posterior residuals obtained after the solution. Data quality control is achieved through multiple rounds of iteration. S4. Station-by-station and system-by-system observation grouping: All observations are grouped according to the dimensions of the station and the satellite navigation system. Each observation group corresponds to an independent variance component, and the variance factors of each group are initialized with equal weight. S5. Simplified variance component estimation and iterative weighted solution: After each round of parameter estimation, calculate the weighted sum of squares and observation count of the posterior residuals of all observations in the observation group. Use the simplified Helmert variance component estimation formula with the redundancy trace correction term omitted to calculate the variance scaling factor of each group. Update the weight matrix of the corresponding observation group with the variance scaling factor. Combine robust estimation to reduce the weight of abnormal residuals. Repeat the parameter estimation and variance component estimation until the variance scaling factor of each group converges to the preset threshold. S6. Ambiguity Fixing and Constraint Re-estimation: Perform integer fixing on floating-point ambiguities and check the fixing quality. After passing the check, apply ambiguity constraints and re-estimate using the weight matrix updated in step S5. Output the multi-system precise orbit and satellite clock bias.

[0009] Preferably, step S1 specifically includes: S11. Acquire GNSS observation data from multiple systems of the global tracking station network, perform cycle slip detection, gross error removal and arc segmentation on the observation data, and apply antenna phase center correction, tidal correction, phase entanglement correction and relativistic effect correction; S12. Construct an ionospherically-free linear combination of dual-frequency pseudorange and carrier phase observations to eliminate the first-order ionospheric delay. The ionospherically-free combined observation equation is expressed as follows: Pseudo-distance ionosphere-free combination: ; Carrier-free ionosphere combination: ; in, For the measuring station, For satellites, , For two carrier frequencies, The speed of light in a vacuum. , These are dual-frequency pseudorange observations. , These are dual-frequency carrier phase observations. For the station To satellite geometric distance, For station clock error, For satellite clock bias, For tropospheric delay, This indicates the ambiguity parameters in the absence of an ionosphere. , These are the combined ionosphere-free bias terms for pseudorange and carrier wave, respectively. , These represent the ionospheric noise of pseudorange and carrier wave combined observations and the unmodeled residual error, respectively.

[0010] Preferably, step S2 specifically includes: S21. Establish the equations of motion for each satellite, which are expressed as an acceleration model: ; in, For satellites, and Satellites Three-dimensional position and velocity vectors in a geocentric inertial coordinate system. For satellite acceleration vector The acceleration is the perturbation acceleration of Earth's gravitational field, calculated using the Earth's gravitational potential spherical harmonic expansion model. The gravitational acceleration is due to the third body of the Sun and Moon. , These are the position vectors of the Sun and the Moon in the geocentric inertial coordinate system, respectively. This refers to the solar radiation pressure perturbation acceleration. This is a vector of empirical parameters for solar radiation pressure. The perturbation acceleration caused by Earth's solid tides and ocean tides. This refers to the perturbation acceleration caused by general relativistic effects; S22. The satellite motion equations from step S21 are numerically integrated using the Adams-Cowell multistep prediction-correction method, with the integration step size set to... (Unit: seconds), the The predicted and corrected values ​​for each step are expressed as follows: Forecast step: ; Calibration step: ; in, For the first The satellite position vector of the step, For the first The satellite position vector of the step, For the first The predicted position vector of the step, For the first The position vector after step correction, Let the order be the order of Adams' method. Indicates the summation index. These are the difference coefficients of the Adams-Cowell forecast formula. , These are the difference coefficients of the Adams-Cowell correction formula. ,in, Corresponding to the number to be determined Step acceleration, For the first The acceleration vector of the step; in the integration initiation stage, the Runge-Kutta-Fehlberg single-step method is used to provide the preceding acceleration vector for the Adams method. Initial values ​​for each step; S23. While integrating the equations of motion, simultaneously integrate the variational equations to obtain the partial derivatives required for the design matrix. The variational equations are expressed as follows: ; in, , indicating satellite The six-dimensional state vector, For satellite The initial state vector at the beginning of the integration epoch. This is a vector of dynamic parameters, containing empirical parameters of solar radiation pressure and force model parameters to be estimated. For time, For dynamic functions For the state vector The Jacobi matrix, with dimensions of , Defined as , For dynamic functions For dynamic parameter vectors The partial derivative matrix has dimensions of , The number of dynamic parameters, The state transition matrix has dimensions [missing information]. Describes how the initial state perturbation propagates to the current epoch. Let be the parameter sensitivity matrix, with dimension . The state transition matrix and parameter sensitivity matrix are obtained by synchronously integrating with the satellite motion equations, and together they constitute the orbit-related sub-block of the design matrix in step S3.

[0011] Preferably, step S3 specifically includes: S31. Linearize the ionosphere-free combined observation equations from step S12 at the reference orbit to obtain the residual equations: ; in, The OC vector is the difference between the observed and calculated values, with dimension . , For the total number of observations, The parameters to be estimated include initial orbital values, dynamic parameters, clock errors, station coordinates, Earth orientation parameters, tropospheric delay, and ambiguity, with dimensions of [dimension missing]. , The total number of parameters to be estimated. for The design matrix contains the partial derivatives of the observations with respect to each parameter. For the observed noise vector; Under the prior assumption of equal weight, a uniform prior variance is assigned to all stations and the satellite navigation system, and the normal equation is constructed by summing the variances station by station: ; in, For the station Design matrix sub-blocks, For the station The prior weight matrix, The total number of stations; the weight matrix of each station in the a priori equal-weight stage. It relies solely on the elevation angle function, without distinguishing between navigation systems; parameter estimates are obtained by solving the normal equations. Calculate the a posteriori residuals ; S32. Perform residual editing and iterative quality control. Set a progressively decreasing residual threshold sequence. In each iteration, remove or mark observations whose posterior residuals exceed the current threshold. Reconstruct the normal equation and solve it. Iterate repeatedly until the residual sequence is stable to complete data quality control. This provides the processed observation dataset for the observation grouping in step S4 and the variance component estimation in step S5.

[0012] Preferably, in step S4, all observations are divided according to the dimensions of the station and the satellite navigation system. A number of non-overlapping observation groups, among which The total number of stations participating in orbit determination. For the number of navigation systems participating in orbit determination, factor 2 corresponds to two types of observations: pseudorange and carrier phase; for stations... Navigation system Observation type Define observation group The variance factor of the observation group is initialized as follows: ; in, Indicates observation group The variance factor at the beginning of the iteration is set to an initial value of 1, indicating that each observation group starts with equal weight.

[0013] Preferably, the simplified Helmert variance component estimation in step S5 specifically includes the following steps: S51. According to the variance factors of each current observation group Adjust the weight matrix to adjust the observation group The first in The weights of the observations are adjusted as follows: ; in, For observation index, For observation Based on the prior weights given by the elevation angle function, The current weights are adjusted for variance factor. For observation To which observation group The current variance factor; the normal equations are reconstructed with the adjusted weight matrix and the parameters are solved by Cholesky decomposition; S52. Calculate the posterior residuals. For each observation group, calculate the weighted sum of squares of the posterior residuals of all observations in that group. With observation count : ; in, For observation The posterior residuals, For observation The current weight, For the observation group The number of observations in the data; S53. Estimate the variance scaling factor for each group using the simplified Helmert formula that omits the redundant trace correction term. : ; in, For the observation group The variance scaling factor reflects the degree of matching between the current variance factor and the actual observed noise. For the observation group The current variance factor, For the observation group The number of observations; this simplifies to removing the redundancy correction term from the denominator in the standard Helmert formula. Set to zero, based on observation count. Direct replacement of redundancy , Represents the trace operation of a matrix. The coefficient matrix of the normal equation The inverse matrix, For the observation group The corresponding sub-normal equation matrix, For the observation group The corresponding design matrix sub-blocks, For the observation group The corresponding weighted sub-block; S54. Update the variance factors for each observation group: ; in, For the observation group The current variance factor, For the variance scaling factor estimated in step S53, determine the variance scaling factor for all observation groups. Does it meet the convergence condition? , To preset the convergence threshold, When it approaches 1, it indicates that the observation group The variance factor has stabilized. If any observation group does not meet the convergence condition, return to steps S51 to S54.

[0014] Preferably, the robust estimation in step S5 specifically involves: in each round of iterative weighted solution, using the iterative weighted least squares method to estimate the posterior residuals exceeding... Observations with a unit weighted mean error are either downweighted or eliminated. The weighting threshold factor is used, and the weighting process employs an equivalent weight function, with the posterior residual... Equivalent rights Specifically, it is expressed as follows: ; in, For observation Equivalent rights after robust weight reduction For observation The current weights after adjustment by variance component estimation. For observation The posterior residuals, For observation To which observation group The current unit weighted mean square error estimate is specifically expressed as , For the observation group The weighted sum of squares of the posterior residuals, For the observation group The number of observations; the robust estimation and the variance component estimation in step S5 are performed together within the same iterative framework.

[0015] Preferably, step S6 specifically includes: S61. Perform double-difference wide-lane and narrow-lane integer search fixation on floating-point ambiguity; S62. Using the weight matrix updated in step S5, the fixed ambiguity is applied to the normal equation to re-estimate the parameters. S63. Perform a ratio test on the fixed mass. If the test passes, proceed to the next step; otherwise, return to step S61. S64. Update the orbital parameters based on the re-estimation results, and perform the simplified variance component estimation and iterative weighted solution in step S5 again to update the weight matrix after the ambiguity is fixed. S65. After repeating steps S61 to S64 for a preset number of rounds, output multi-system precision orbit and satellite clock error products.

[0016] Therefore, the present invention employs the above-mentioned method for weighting adjustment observations of a multi-system navigation satellite joint orbit determination network, which has the following beneficial effects: (1) The present invention adopts a multi-granularity observation grouping strategy for each station and each system, dividing all observations into non-overlapping observation groups according to the dimensions of the station and the satellite navigation system. Each observation group independently estimates the variance components, which can finely characterize the non-homogeneous characteristics of observations caused by differences in receiver model, antenna type and site environment of different stations and differences in orbit type and signal system of different navigation systems. This overcomes the limitations of traditional fixed prior variance or empirical weight ratios that cannot adapt to changes in actual observation conditions. (2) The present invention adopts a simplified Helmert variance component estimation formula that omits the redundancy trace correction term and directly replaces the redundancy correction term that requires calculating the trace of the product of the coefficient matrix of the normal equation and the sub-normal equation matrix with the observation count. This reduces the computational cost of variance component estimation from relying on large-scale matrix inversion and storage to only requiring the calculation of the weighted sum of squares of the posterior residuals and the observation count. This makes variance component estimation not significantly increase the computational burden in the joint orbit determination of multi-system navigation satellites with parameter scales of up to 100,000 dimensions, and meets the requirements of operational orbit determination for computational timeliness. (3) This invention integrates robust estimation and variance component estimation into the same iterative weighted least squares framework for joint execution. After each round of variance component update, equivalent weight reduction processing is applied to the observations whose posterior residuals exceed the threshold, so that the variance component estimation is not affected by a small number of gross errors or abnormal observations, thereby enhancing the robustness of the weighting results in complex observation environments and reducing the impact of abnormal observations on the accuracy of orbit calculation. (4) This invention integrates three technologies: station-by-station and system-by-system adaptive grouping and weighting, simplified variance component estimation, and robust iterative weighted solution. These three technologies form a data-driven closed-loop collaborative relationship, which effectively improves the orbit accuracy and solution reliability of joint orbit determination of multi-system navigation satellites without significantly increasing the computational burden. This provides an efficient and reliable technical solution for the high-precision generation of precision orbit products for global navigation satellite systems.

[0017] The technical solution of the present invention will be further described in detail below with reference to the accompanying drawings and embodiments. Attached Figure Description

[0018] Figure 1 This is a flowchart of a method for weighting adjustment observations in a multi-system navigation satellite joint orbit determination network according to the present invention; Figure 2 This is a three-dimensional RMS statistical result diagram of the orbit DBD for each navigation system and the global average under two weighting schemes in the embodiments of the present invention. Detailed Implementation

[0019] The technical solution of the present invention will be further described below with reference to the accompanying drawings and embodiments.

[0020] Unless otherwise defined, the technical or scientific terms used in this invention shall have the ordinary meaning understood by one of ordinary skill in the art to which this invention pertains. The terms "first," "second," and similar terms used in this invention do not indicate any order, quantity, or importance, but are merely used to distinguish different components. Terms such as "comprising" or "including" mean that the element or object preceding the word encompasses the elements or objects listed following the word and their equivalents, without excluding other elements or objects. Terms such as "connected" or "linked" are not limited to physical or mechanical connections, but can include electrical connections, whether direct or indirect. Terms such as "upper," "lower," "left," and "right" are used only to indicate relative positional relationships; when the absolute position of the described object changes, the relative positional relationship may also change accordingly.

[0021] Example like Figure 1 As shown, a method for weighting adjustment observations in a multi-system navigation satellite joint orbit determination network includes the following steps: S1. Preprocessing and quality control of observation data: Acquire multi-system observation data and external model data from the global tracking station network, perform cycle slip detection, gross error removal and arc segment division on the multi-system observation data, and apply observation corrections. S11. Acquire GNSS observation data from multiple systems of the global tracking station network, perform cycle slip detection, gross error removal and arc segmentation on the observation data, and apply antenna phase center correction, tidal correction, phase entanglement correction and relativistic effect correction; S12. Construct an ionospherically-free linear combination of dual-frequency pseudorange and carrier phase observations to eliminate the first-order ionospheric delay. The ionospherically-free combined observation equation is expressed as follows: Pseudo-distance ionosphere-free combination: ; Carrier-free ionosphere combination: ; in, For the measuring station, For satellites, , For two carrier frequencies, The speed of light in a vacuum. , These are dual-frequency pseudorange observations. , These are dual-frequency carrier phase observations. For the station To satellite geometric distance, For station clock error, For satellite clock bias, For tropospheric delay, This indicates the ambiguity parameters in the absence of an ionosphere. , These are the combined ionosphere-free bias terms for pseudorange and carrier wave, respectively. , These represent the ionospheric noise of pseudorange and carrier wave combined observations and the unmodeled residual error, respectively.

[0022] S2. Numerical integration and fitting of orbits: Numerical integration is performed on multi-system satellites based on the satellite dynamics model to obtain the reference orbit, state transition matrix and parameter sensitivity matrix, and orbit fitting and initial value correction are performed based on the external orbit. S21. Establish the equations of motion for each satellite, which are expressed as an acceleration model: ; in, For satellites, and Satellites Three-dimensional position and velocity vectors in a geocentric inertial coordinate system. For satellite acceleration vector The acceleration is the perturbation acceleration of Earth's gravitational field, calculated using the Earth's gravitational potential spherical harmonic expansion model. The gravitational acceleration is due to the third body of the Sun and Moon. , These are the position vectors of the Sun and the Moon in the geocentric inertial coordinate system, respectively. This refers to the solar radiation pressure perturbation acceleration. This is a vector of empirical parameters for solar radiation pressure. The perturbation acceleration caused by Earth's solid tides and ocean tides. This refers to the perturbation acceleration caused by general relativistic effects; S22. The satellite motion equations from step S21 are numerically integrated using the Adams-Cowell multistep prediction-correction method, with the integration step size set to... (Unit: seconds), the The predicted and corrected values ​​for each step are expressed as follows: Forecast step: ; Calibration step: ; in, For the first The satellite position vector of the step, For the first The satellite position vector of the step, For the first The predicted position vector of the step, For the first The position vector after step correction, Let the order be the order of Adams' method. Indicates the summation index. These are the difference coefficients of the Adams-Cowell forecast formula. , These are the difference coefficients of the Adams-Cowell correction formula. ,in, Corresponding to the number to be determined Step acceleration, For the first The acceleration vector of the step; in the integration initiation stage, the Runge-Kutta-Fehlberg single-step method is used to provide the preceding acceleration vector for the Adams method. Initial values ​​for each step; S23. While integrating the equations of motion, simultaneously integrate the variational equations to obtain the partial derivatives required for the design matrix. The variational equations are expressed as follows: ; in, , indicating satellite The six-dimensional state vector, For satellite The initial state vector at the beginning of the integration epoch. This is a vector of dynamic parameters, containing empirical parameters of solar radiation pressure and force model parameters to be estimated. For time, For dynamic functions For the state vector The Jacobi matrix, with dimensions of , Defined as , For dynamic functions For dynamic parameter vectors The partial derivative matrix has dimensions of , The number of dynamic parameters, The state transition matrix has dimensions [missing information]. Describes how the initial state perturbation propagates to the current epoch. Let be the parameter sensitivity matrix, with dimension . The state transition matrix and parameter sensitivity matrix are obtained by synchronously integrating with the satellite motion equations, and together they constitute the orbit-related sub-block of the design matrix in step S3.

[0023] S3, Prior Equal Weight Parameter Estimation and Residual Editing: Construct normal equations with prior equal weights and solve them by least squares. Perform residual editing based on the posterior residuals obtained after the solution. Data quality control is achieved through multiple rounds of iteration. S31. Linearize the ionosphere-free combined observation equations from step S12 at the reference orbit to obtain the residual equations: ; in, The OC vector is the difference between the observed and calculated values, with dimension . , For the total number of observations, The parameters to be estimated include initial orbital values, dynamic parameters, clock errors, station coordinates, Earth orientation parameters, tropospheric delay, and ambiguity, with dimensions of [dimension missing]. , The total number of parameters to be estimated. for The design matrix contains the partial derivatives of the observations with respect to each parameter. For the observed noise vector; Under the prior assumption of equal weight, a uniform prior variance is assigned to all stations and the satellite navigation system, and the normal equation is constructed by summing the variances station by station: ; in, For the station Design matrix sub-blocks, For the station The prior weight matrix, The total number of stations; the weight matrix of each station in the a priori equal-weight stage. It relies solely on the elevation angle function, without distinguishing between navigation systems; parameter estimates are obtained by solving the normal equations. Calculate the a posteriori residuals ; S32. Perform residual editing and iterative quality control. Set a progressively decreasing residual threshold sequence. In each iteration, remove or mark observations whose posterior residuals exceed the current threshold. Reconstruct the normal equation and solve it. Iterate repeatedly until the residual sequence is stable to complete data quality control. This provides the processed observation dataset for the observation grouping in step S4 and the variance component estimation in step S5.

[0024] S4. Station-by-Station, System-by-System Observation Grouping: Group all observations according to the station and satellite navigation system dimensions. Each observation group corresponds to an independent variance component, and the variance factors of each group are initialized with equal weights. Divide all observations according to the station and satellite navigation system dimensions into... A number of non-overlapping observation groups, among which The total number of stations participating in orbit determination. For the number of navigation systems participating in orbit determination, factor 2 corresponds to two types of observations: pseudorange and carrier phase; for stations... Navigation system Observation type Define observation group The variance factor of the observation group is initialized as follows: ; in, Indicates observation group The variance factor at the beginning of the iteration is set to an initial value of 1, indicating that each observation group starts with equal weight.

[0025] S5. Simplified variance component estimation and iterative weighted solution: After each round of parameter estimation, calculate the weighted sum of squares and observation count of the posterior residuals of all observations in the observation group. Use the simplified Helmert variance component estimation formula with the redundancy trace correction term omitted to calculate the variance scaling factor of each group. Update the weight matrix of the corresponding observation group with the variance scaling factor. Combine robust estimation to reduce the weight of abnormal residuals. Repeat the parameter estimation and variance component estimation until the variance scaling factor of each group converges to the preset threshold. Simplified Helmert variance component estimation specifically includes the following steps: S51. According to the variance factors of each current observation group Adjust the weight matrix to adjust the observation group The first in The weights of the observations are adjusted as follows: ; in, For observation index, For observation Based on the prior weights given by the elevation angle function, The current weights are adjusted for variance factor. For observation To which observation group The current variance factor; the normal equations are reconstructed with the adjusted weight matrix and the parameters are solved by Cholesky decomposition; S52. Calculate the posterior residuals. For each observation group, calculate the weighted sum of squares of the posterior residuals of all observations in that group. With observation count : ; in, For observation The posterior residuals, For observation The current weight, For the observation group The number of observations in the data; S53. Estimate the variance scaling factor for each group using the simplified Helmert formula that omits the redundant trace correction term. : ; in, For the observation group The variance scaling factor reflects the degree of matching between the current variance factor and the actual observed noise. For the observation group The current variance factor, For the observation group The number of observations; this simplifies to removing the redundancy correction term from the denominator in the standard Helmert formula. Set to zero, based on observation count. Direct replacement of redundancy , Represents the trace operation of a matrix. The coefficient matrix of the normal equation The inverse matrix, For the observation group The corresponding sub-normal equation matrix, For the observation group The corresponding design matrix sub-blocks, For the observation group The corresponding weighted sub-block; S54. Update the variance factors for each observation group: ; in, For the observation group The current variance factor, For the variance scaling factor estimated in step S53, determine the variance scaling factor for all observation groups. Does it meet the convergence condition? , To preset the convergence threshold, When it approaches 1, it indicates that the observation group The variance factor has stabilized. If any observation group does not meet the convergence condition, return to steps S51 to S54.

[0026] The robust estimation in step S5 specifically involves: in each round of iterative weighted solution, using the iterative weighted least squares method, for posterior residuals exceeding... Observations with a unit weighted mean error are either downweighted or eliminated. The weighting threshold factor is used, and the weighting process employs an equivalent weight function, with the posterior residual... Equivalent rights Specifically, it is expressed as follows: ; in, For observation Equivalent rights after robust weight reduction For observation The current weights after adjustment by variance component estimation. For observation The posterior residuals, For observation To which observation group The current unit weighted mean square error estimate is specifically expressed as , For the observation group The weighted sum of squares of the posterior residuals, For the observation group The number of observations; the robust estimation and the variance component estimation in step S5 are performed together within the same iterative framework.

[0027] S6. Ambiguity Fixing and Constraint Re-estimation: Perform integer fixing on floating-point ambiguities and check the fixing quality. After passing the check, apply ambiguity constraints and re-estimate using the weight matrix updated in step S5. Output the multi-system precise orbit and satellite clock bias.

[0028] S61. Perform double-difference wide-lane and narrow-lane integer search fixation on floating-point ambiguity; S62. Using the weight matrix updated in step S5, the fixed ambiguity is applied to the normal equation to re-estimate the parameters. S63. Perform a ratio test on the fixed mass. If the test passes, proceed to the next step; otherwise, return to step S61. S64. Update the orbital parameters based on the re-estimation results, and perform the simplified variance component estimation and iterative weighted solution in step S5 again to update the weight matrix after the ambiguity is fixed. S65. After repeating steps S61 to S64 for a preset number of rounds, output multi-system precision orbit and satellite clock error products.

[0029] To verify the technical effectiveness of the method of the present invention, the following description is based on specific experimental data.

[0030] Example 1 This embodiment addresses the weighting problem of adjustment observations in a multi-system navigation satellite joint orbit determination network, conducting a comparative study between the proposed method (simplified VCE adaptive weighting) and the traditional equal-weighting scheme. The experiment uses observation data from approximately 120 IGS tracking stations globally, performing joint precise orbit determination calculations on approximately 100 satellites from four navigation systems: GPS, GLONASS, Galileo, and BDS. The orbit determination arc is set to 24 hours, the data sampling interval is set to 300 seconds, and the experimental period covers 7 days from day 68 to day 74 of 2025 (corresponding to GPS week 2357). Orbit accuracy is evaluated using DayBoundary Discontinuity (DBD) as the accuracy index, calculating the root mean square (RMS) of the three-dimensional position discontinuities of adjacent day orbits at the boundary.

[0031] Under the same experimental conditions, the two weighting schemes are as follows: (1) The simplified VCE adaptive weighting scheme proposed in this invention estimates the variance components group by group according to the dimensions of the station and the navigation system and iteratively updates the weight matrix; (2) The traditional equal weight scheme assigns the same prior weight ratio of 1:1:1:1 to the four systems GPS, GLONASS, Galileo and BDS.

[0032] Apart from the difference in the observation weighting strategy, the two schemes are completely identical in terms of preprocessing procedures, dynamic models, parameter settings, and ambiguity fixing strategies.

[0033] The orbital boundary jump reflects the positional continuity of the orbit determination solutions for adjacent days at the boundary, and is an important indicator for measuring the orbit determination accuracy of navigation satellites. For example... Figure 2 As shown in Table 1, the three-dimensional RMS statistics of orbit DBD for each navigation system and the global average are listed under the two weighting schemes.

[0034] Table 1. Three-dimensional RMS statistics of orbital date jumps for two weighting schemes.

[0035] It can be seen that the VCE scheme outperforms the equal-weighted scheme in all four navigation systems. The Galileo system shows the most significant improvement, with its DBD 3D RMS decreasing from 47.1mm to 38.6mm, an improvement of 18.0%; the BDS system from 66.6mm to 58.3mm, an improvement of 12.4%; the GPS system from 41.2mm to 37.9mm, an improvement of 8.1%; and the GLONASS system by 1.6%. Overall, the VCE scheme's DBD 3D RMS is 55.7mm, a reduction of 4.8mm compared to the equal-weighted scheme's 60.5mm, representing an improvement of 8.0%.

[0036] The experimental results above demonstrate that the simplified VCE adaptive weighting method proposed in this invention achieved stable improvements in orbit accuracy during the 7-day experiment. The global orbit DBD 3D RMS was reduced by 8.0% compared to the traditional equal-weighting scheme, with the Galileo system showing the highest improvement of up to 18.0%. Therefore, by employing the above method, this invention can replace the traditional fixed prior weight ratio with data-driven adaptive variance component estimation without significantly increasing the computational burden, effectively improving the overall orbit accuracy and solution reliability of joint orbit determination for multi-system navigation satellites.

[0037] Therefore, the present invention employs the above-mentioned method for weighting adjustment observations of a multi-system navigation satellite joint orbit determination network, which has the following beneficial effects: (1) The present invention adopts a multi-granularity observation grouping strategy for each station and each system, dividing all observations into non-overlapping observation groups according to the dimensions of the station and the satellite navigation system. Each observation group independently estimates the variance components, which can finely characterize the non-homogeneous characteristics of observations caused by differences in receiver model, antenna type and site environment of different stations and differences in orbit type and signal system of different navigation systems. This overcomes the limitations of traditional fixed prior variance or empirical weight ratios that cannot adapt to changes in actual observation conditions. (2) The present invention adopts a simplified Helmert variance component estimation formula that omits the redundancy trace correction term and directly replaces the redundancy correction term that requires calculating the trace of the product of the coefficient matrix of the normal equation and the sub-normal equation matrix with the observation count. This reduces the computational cost of variance component estimation from relying on large-scale matrix inversion and storage to only requiring the calculation of the weighted sum of squares of the posterior residuals and the observation count. This makes variance component estimation not significantly increase the computational burden in the joint orbit determination of multi-system navigation satellites with parameter scales of up to 100,000 dimensions, and meets the requirements of operational orbit determination for computational timeliness. (3) This invention integrates robust estimation and variance component estimation into the same iterative weighted least squares framework for joint execution. After each round of variance component update, equivalent weight reduction processing is applied to the observations whose posterior residuals exceed the threshold, so that the variance component estimation is not affected by a small number of gross errors or abnormal observations, thereby enhancing the robustness of the weighting results in complex observation environments and reducing the impact of abnormal observations on the accuracy of orbit calculation. (4) This invention integrates three technologies: station-by-station and system-by-system adaptive grouping and weighting, simplified variance component estimation, and robust iterative weighted solution. These three technologies form a data-driven closed-loop collaborative relationship, which effectively improves the orbit accuracy and solution reliability of joint orbit determination of multi-system navigation satellites without significantly increasing the computational burden. This provides an efficient and reliable technical solution for the high-precision generation of precision orbit products for global navigation satellite systems.

[0038] 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 preferred embodiments, those skilled in the art should understand that modifications or equivalent substitutions can still be made to the technical solutions of the present invention, and these modifications or equivalent substitutions cannot cause the modified technical solutions to deviate from the spirit and scope of the technical solutions of the present invention.

Claims

1. A method for weighting adjustment observations of a multi-system navigation satellite joint orbit determination network, characterized in that, Includes the following steps: S1. Preprocessing and quality control of observation data: Acquire multi-system observation data and external model data from the global tracking station network, perform cycle slip detection, gross error removal and arc segment division on the multi-system observation data, and apply observation corrections. S2. Numerical integration and fitting of orbits: Numerical integration is performed on multi-system satellites based on the satellite dynamics model to obtain the reference orbit, state transition matrix and parameter sensitivity matrix, and orbit fitting and initial value correction are performed based on the external orbit. S3, Prior Equal Weight Parameter Estimation and Residual Editing: Construct normal equations with prior equal weights and solve them by least squares. Perform residual editing based on the posterior residuals obtained after the solution. Data quality control is achieved through multiple rounds of iteration. S4. Station-by-station and system-by-system observation grouping: All observations are grouped according to the dimensions of the station and the satellite navigation system. Each observation group corresponds to an independent variance component, and the variance factors of each group are initialized with equal weight. S5. Simplified variance component estimation and iterative weighted solution: After each round of parameter estimation, calculate the weighted sum of squares and observation count of the posterior residuals of all observations in the observation group. Use the simplified Helmert variance component estimation formula with the redundancy trace correction term omitted to calculate the variance scaling factor of each group. Update the weight matrix of the corresponding observation group with the variance scaling factor. Combine robust estimation to reduce the weight of abnormal residuals. Repeat the parameter estimation and variance component estimation until the variance scaling factor of each group converges to the preset threshold. S6. Ambiguity Fixing and Constraint Re-estimation: Perform integer fixing on floating-point ambiguities and check the fixing quality. After passing the check, apply ambiguity constraints and re-estimate using the weight matrix updated in step S5. Output the multi-system precise orbit and satellite clock bias.

2. The method for weighting adjustment observations of a multi-system navigation satellite joint orbit determination network according to claim 1, characterized in that, Step S1 specifically includes: S11. Acquire GNSS observation data from multiple systems of the global tracking station network, perform cycle slip detection, gross error removal and arc segmentation on the observation data, and apply antenna phase center correction, tidal correction, phase entanglement correction and relativistic effect correction; S12. Construct an ionospherically-free linear combination of dual-frequency pseudorange and carrier phase observations to eliminate the first-order ionospheric delay. The ionospherically-free combined observation equation is expressed as follows: Pseudo-distance ionosphere-free combination: ； Carrier-free ionosphere combination: ； in, For the measuring station, For satellites, , For two carrier frequencies, The speed of light in a vacuum. , These are dual-frequency pseudorange observations. , These are dual-frequency carrier phase observations. For the station To satellite geometric distance, For station clock error, For satellite clock bias, For tropospheric delay, This indicates the ambiguity parameters in the absence of an ionosphere. , These are the combined ionosphere-free bias terms for pseudorange and carrier wave, respectively. , These represent the ionospheric noise of pseudorange and carrier wave combined observations and the unmodeled residual error, respectively.

3. The method for weighting adjustment observations of a multi-system navigation satellite joint orbit determination network according to claim 1, characterized in that, Step S2 specifically includes: S21. Establish the equations of motion for each satellite, which are expressed as an acceleration model: ； in, For satellites, and Satellites Three-dimensional position and velocity vectors in a geocentric inertial coordinate system. For satellite acceleration vector The acceleration due to the perturbation of Earth's gravitational field. The gravitational acceleration is due to the third body of the Sun and Moon. , These are the position vectors of the Sun and the Moon in the geocentric inertial coordinate system, respectively. This refers to the solar radiation pressure perturbation acceleration. This is a vector of empirical parameters for solar radiation pressure. The perturbation acceleration caused by Earth's solid tides and ocean tides. This refers to the perturbation acceleration caused by general relativistic effects; S22. The satellite motion equations from step S21 are numerically integrated using the Adams-Cowell multistep prediction-correction method, with the integration step size set to... , No. The predicted and corrected values ​​for each step are expressed as follows: Forecast step: ； Calibration step: ； in, For the first The satellite position vector of the step, For the first The satellite position vector of the step, For the first The predicted position vector of the step, For the first The position vector after step correction, Let the order be the order of Adams' method. Indicates the summation index. These are the difference coefficients of the Adams-Cowell forecast formula. , These are the difference coefficients of the Adams-Cowell correction formula. ,in, Corresponding to the number to be determined Step acceleration, For the first The acceleration vector of the step; in the integration initiation stage, the Runge-Kutta-Fehlberg single-step method is used to provide the preceding acceleration vector for the Adams method. Initial values ​​for each step; S23. While integrating the equations of motion, simultaneously integrate the variational equations to obtain the partial derivatives required for the design matrix. The variational equations are expressed as follows: ； in, , indicating satellite The six-dimensional state vector, For satellite The initial state vector at the beginning of the integration epoch. For dynamic parameter vectors, For time, For dynamic functions For the state vector The Jacobi matrix, with dimensions of , Defined as , For dynamic functions For dynamic parameter vectors The partial derivative matrix has dimensions of , The number of dynamic parameters, The state transition matrix has dimensions [missing information]. , Let be the parameter sensitivity matrix, with dimension . The state transition matrix and the parameter sensitivity matrix are obtained by synchronous integration with the satellite motion equation, and together they constitute the orbit-related sub-block of the design matrix in step S3.

4. The method for weighting adjustment observations of a multi-system navigation satellite joint orbit determination network according to claim 2, characterized in that, Step S3 specifically includes: S31. Linearize the ionosphere-free combined observation equations from step S12 at the reference orbit to obtain the residual equations: ； in, The OC vector is the difference between the observed and calculated values, with dimension . , For the total number of observations, The parameters to be estimated include initial orbital values, dynamic parameters, clock errors, station coordinates, Earth orientation parameters, tropospheric delay, and ambiguity, with dimensions of [dimension missing]. , The total number of parameters to be estimated. for The design matrix contains the partial derivatives of the observations with respect to each parameter. For the observed noise vector; Under the prior assumption of equal weight, a uniform prior variance is assigned to all stations and the satellite navigation system, and the normal equation is constructed by summing the variances station by station: ； in, For the station Design matrix sub-blocks, For the station The prior weight matrix, The total number of stations; the weight matrix of each station in the a priori equal-weight stage. It relies solely on the elevation angle function, without distinguishing between navigation systems; parameter estimates are obtained by solving the normal equations. Calculate the a posteriori residuals ; S32. Perform residual editing and iterative quality control. Set a progressively decreasing residual threshold sequence. In each iteration, remove or mark observations whose posterior residuals exceed the current threshold. Reconstruct the normal equation and solve it. Iterate repeatedly until the residual sequence is stable to complete data quality control. This provides the processed observation dataset for the observation grouping in step S4 and the variance component estimation in step S5.

5. The method for weighting adjustment observations of a multi-system navigation satellite joint orbit determination network according to claim 1, characterized in that, In step S4, all observations are divided according to the dimensions of the station and the satellite navigation system. A number of non-overlapping observation groups, among which The total number of stations participating in orbit determination. For the number of navigation systems participating in orbit determination, factor 2 corresponds to two types of observations: pseudorange and carrier phase; for stations... Navigation system Observation type Define observation group The variance factor of the observation group is initialized as follows: ； in, Indicates observation group The variance factor at the beginning of the iteration is set to an initial value of 1, indicating that each observation group starts with equal weight.

6. The method for weighting adjustment observations of a multi-system navigation satellite joint orbit determination network according to claim 1, characterized in that, Step S5, which simplifies the Helmert variance component estimation, specifically includes the following steps: S51. According to the variance factors of each current observation group Adjust the weight matrix to adjust the observation group The first in The weights of the observations are adjusted as follows: ； in, For observation index, For observation Based on the prior weights given by the elevation angle function, The current weights are adjusted for variance factor. For observation To which observation group The current variance factor; the normal equations are reconstructed with the adjusted weight matrix and the parameters are solved by Cholesky decomposition; S52. Calculate the posterior residuals. For each observation group, calculate the weighted sum of squares of the posterior residuals of all observations in that group. With observation count : ； in, For observation The posterior residuals, For observation The current weight, For the observation group The number of observations in the data; S53. Estimate the variance scaling factor for each group using the simplified Helmert formula that omits the redundant trace correction term. : ； in, For the observation group The variance scaling factor, For the observation group The current variance factor, For the observation group The number of observations; this simplifies to removing the redundancy correction term from the denominator in the standard Helmert formula. Set to zero, based on observation count. Direct replacement of redundancy , Represents the trace operation of a matrix. The coefficient matrix of the normal equation The inverse matrix, For the observation group The corresponding sub-normal equation matrix, For the observation group The corresponding design matrix sub-blocks, For the observation group The corresponding weighted sub-block; S54. Update the variance factors for each observation group: ； in, For the observation group The current variance factor, For the variance scaling factor estimated in step S53, determine the variance scaling factor for all observation groups. Does it meet the convergence condition? , To preset the convergence threshold, When it approaches 1, it indicates that the observation group The variance factor has stabilized. If any observation group does not meet the convergence condition, return to steps S51 to S54.

7. The method for weighting adjustment observations of a multi-system navigation satellite joint orbit determination network according to claim 1, characterized in that, The robust estimation in step S5 specifically involves: in each round of iterative weighted solution, using the iterative weighted least squares method, for posterior residuals exceeding... Observations with a unit weighted mean error are either downweighted or eliminated. The weighting threshold factor is used, and the weighting process employs an equivalent weight function, with the posterior residual... Equivalent rights Specifically, it is expressed as follows: ； in, For observation Equivalent rights after robust weight reduction For observation The current weights after adjustment by variance component estimation. For observation The posterior residuals, For observation To which observation group The current unit weighted mean square error estimate is specifically expressed as , For the observation group The weighted sum of squares of the posterior residuals, For the observation group The number of observations; the robust estimation and the variance component estimation in step S5 are performed together within the same iterative framework.

8. The method for weighting adjustment observations of a multi-system navigation satellite joint orbit determination network according to claim 1, characterized in that, Step S6 specifically includes: S61. Perform double-difference wide-lane and narrow-lane integer search fixation on floating-point ambiguity; S62. Using the weight matrix updated in step S5, the fixed ambiguity is applied to the normal equation to re-estimate the parameters. S63. Perform a ratio test on the fixed mass. If the test passes, proceed to the next step; otherwise, return to step S61. S64. Update the orbital parameters based on the re-estimation results, and perform the simplified variance component estimation and iterative weighted solution in step S5 again to update the weight matrix after the ambiguity is fixed. S65. After repeating steps S61 to S64 for a preset number of rounds, output multi-system precision orbit and satellite clock error products.

Images