GNSS (Global Navigation Satellite System) multi-source precision satellite orbit product synthesis method and equipment

Abstract

The invention discloses a GNSS (Global Navigation Satellite System) multi-source precision satellite orbit product synthesis method and equipment, and belongs to the technical field of satellite navigation and geodetic survey. Precise orbit data of multiple analysis centers are obtained, and orbit pre-alignment is carried out; the satellites are grouped; calculating radial, tangential and normal residual error sequences of each central orbit relative to the reference orbit, and performing abnormal value detection and elimination based on a corrected Z-fraction method; defining a core satellite in each satellite group; the core satellite residual error is used as an observation value, the weight of each analysis center is determined by using a least square variance component estimation method, and satellite laser ranging prior information is introduced to enhance the estimation stability; carrying out weighted average on the track according to the weight to obtain a comprehensive track; and carrying out iterative updating until convergence. According to the method, the precision difference of data of different sources is fully considered, and a comprehensive orbit product with the precision and reliability superior to those of a single analysis center can be automatically and objectively generated.

Description

Technical Field

[0001] This invention relates to the field of satellite navigation and geodesy technology, specifically to a method and equipment for integrating GNSS multi-source precision satellite orbit products. Background Technology

[0002] Precise satellite orbit products are fundamental data for satellite navigation, high-precision positioning, and Earth science research. Currently, organizations such as the International GNSS Service Organization (IGS), the International BeiDou Navigation Satellite System (iGMAS), and the GNSS Monitoring and Evaluation System (IGAS) have aggregated multi-system GNSS precise orbit products from multiple analysis centers worldwide. These products cover data from multiple constellations, including GPS, GLONASS, Galileo, and BeiDou. Due to differences in data processing strategies, software, and models among the analysis centers, their products vary in accuracy and reliability. Directly using products from a single analysis center carries risks, while simple arithmetic averages cannot objectively reflect the true accuracy level of each product. To obtain higher-precision satellite orbit products, comprehensive processing of multi-source GNSS precise satellite orbit products is necessary.

[0003] The prior art, disclosed in CN106680845A, provides a satellite orbit synthesis weighting method. First, it inputs satellite orbit data from each analysis center and uses a segmented method to determine the User Ranging Error (URE) of the satellites. Then, it calculates the weights of the analysis centers and satellites based on the URE, and performs weighted synthesis on the satellite orbits of all analysis centers to obtain a high-precision, high-reliability synthesized orbit. However, this method is only applicable to satellite orbit synthesis for GEO / IGSO / MEO multi-orbit types and does not fully consider the accuracy differences between different data sources, resulting in limited improvement in the accuracy of the synthesized orbit product. Summary of the Invention

[0004] To address the shortcomings of existing technologies, a method and equipment for integrating GNSS multi-source precision satellite orbit products are proposed. This method can systematically process orbit data from multiple analysis centers. Through automated orbit pre-alignment, robust outlier detection, data-driven satellite grouping, optimal weight determination based on variance component estimation, and rigorous iterative convergence control and external verification, it ultimately generates an integrated orbit product with accuracy and reliability superior to any single analysis center product.

[0005] To achieve the above objectives, this invention discloses a method for integrating GNSS multi-source precision satellite orbit products, comprising the following steps: S1. Acquire precise orbit data from multiple analysis centers and preprocess the precise orbit data; S2. Perform preliminary pre-alignment of the center orbit based on the pre-processed precision satellite orbit data; S3. The satellites are grouped using a time series clustering method based on dynamic time warping. S4. Calculate the radial, tangential, and normal residual sequences of each analysis center orbit relative to the reference orbit, and remove outlier data based on the corrected Z-score method; S5. Define core satellites within each satellite group; use the orbital residuals of the core satellites in each analysis center as observations to determine the weights of each analysis center, introduce prior information from satellite laser ranging to enhance estimation stability, and perform a weighted average of the orbits based on the weights of each analysis center to obtain the comprehensive orbits of the core satellites. S6. Repeat S4~S5 iterative updates until convergence, to obtain the final integrated orbital product.

[0006] Furthermore, the multiple analysis centers include: the analysis center under the International GNSS Service Organization (IGS) and the analysis center under the Global Continuous Monitoring and Evaluation System (iGMAS); the satellite systems involved include: Global Positioning System (GPS), Global Navigation Satellite System (GLONASS), Galileo, and BeiDou Navigation Satellite System (BDS). Preprocessing of precision orbit data includes integrity checks and gross error removal, among which: Integrity check process: exclude satellites that provide data from only a single analysis center, and ensure that data can be comprehensively compared. That is, exclude satellite epochs with a data integrity rate of less than 80% in a single day to ensure the continuity of the time series. Gross error removal process: Calculate the average of all available analysis center orbit data for each satellite as the initial synthesized orbit for that satellite. Exclude satellite epoch data from the analysis center that differ from the initial integrated orbit by more than 300m; iterate this step 3 times to obtain the integrated orbit result after removing coarsely screened outliers; The initial composite orbit of the satellite is calculated using the following formula. : , in, This is the three-dimensional position vector of the satellite's orbit, where n is the number of analysis centers providing that satellite orbit. For the analysis center index, For the calendar year, Satellite number.

[0007] Furthermore, the rotation parameters of the common synthetic Earth reference frame solution are used for the initial pre-alignment of the analysis center orbits. If the frame rotation parameters are not available, all analysis center orbits are aligned with the initial synthetic orbit coordinates using a Helmert transformation. The process of aligning all analysis center orbits with the initial integrated orbit coordinates using a Helmert transformation is as follows: For each analysis center's orbit, a set of seven parameters relative to the composite orbit is solved. The transformation model of the analysis center relative to the composite orbit is as follows: , In the formula, The initial composite orbital seven parameters to be solved are: These are the three translation components in the coordinate transformation. It consists of three rotational components. For scale parameters; Indicates the satellite's composite orbital coordinates. Represent the center orbit coordinates; establish the error equation based on the indirect adjustment model: , In the formula Represents the residual vector; It is a coefficient matrix; This is a vector of estimated values ​​for the parameter to be estimated; Let N be the vector of observed values; assuming the analysis center has N common coordinate points relative to the composite orbit, and the coefficient matrix... Estimated value vector Observation vector They are represented as follows: , Then, the coordinate transformation parameters are estimated using least squares to minimize the sum of squares of the residuals, i.e.: , Assume that the coordinates of the satellite orbit products in the error equation have the same precision. Treat the coordinates of the satellite orbit products as independent observations with equal precision. Let the weight matrix P be the identity matrix, and the estimated value vector... Represented as: , Using the obtained initial composite orbit seven parameters, the original orbit coordinates of all satellites of the analysis center at all epochs were determined according to the transformation model of the analysis center relative to the composite orbit. Perform the transformation to obtain the orbital coordinates after satellite orbit alignment. : .

[0008] Furthermore, time series clustering is used to group the satellites, as follows: The radial, tangential, and normal residual sequences of the orbits of all satellites and the composite orbits of each satellite system are extracted. Each satellite is treated as a sequence set, forming a multidimensional time series set. If multiple analysis centers are considered, a four-dimensional feature tensor of "analysis center-satellite-direction-epoch" is constructed, and the multidimensional time series is reduced to a two-dimensional time series matrix through vectorization or principal component analysis. The Dynamic Time Warping (DTW) algorithm is used to calculate the similarity distance between any two two-dimensional time series matrices. The satellites are grouped, and the similarity distance matrix is ​​calculated using the following formula. Each element Indicates the first satellite and the first DTW distance of the satellites: , in To align the path, and These are the feature values ​​of the two two-dimensional time series at corresponding time points; and These represent the feature values ​​of the two two-dimensional time series, and These represent specific points in time. The similarity distance matrix calculated by the DTW algorithm is given by the following formula. The K-Means++ time series clustering method was used to group the satellites; the clustering objective was to minimize the squared intra-cluster distance. sum: , Where K is the preset number of clusters, and the optimal value is determined by the silhouette coefficient; This represents the k-th cluster.

[0009] Furthermore, the orbits of each analysis center are compared with the reference orbit, which is the initial synthesis orbit. For each epoch, the positional differences between the two orbits in the radial, tangential, and normal directions are calculated, i.e., the orbital residuals. , i∈n, where n is the total number of analysis centers; for each satellite, the residual sequences in the radial, tangential, and normal directions are respectively... Calculate the corrected Z-score: Calculate the residual sequence using the following formula the median of : , Calculate the absolute value of the difference between each residual and the median in the residual sequence using the following formula, and then find the median of all absolute values: , The corrected Z-score is calculated using the following formula: , in, It is a constant factor; under the assumption that the data strictly follows a normal distribution, MAD is related to the standard deviation. The relationship is ; Examine the corrected Z-fractions for each satellite at all epochs in the radial, tangential, and normal directions. As long as the Z-fraction is corrected in any direction and at any epoch, its... If the threshold value is exceeded, the satellite is marked as abnormal and its historical data is removed.

[0010] Furthermore, for each group of satellites after the satellites are grouped, all available satellites in the analysis center are selected as core satellites. If a satellite group contains only one satellite, that satellite is automatically regarded as a core satellite. At the same time, satellites whose absolute value of the corrected Z score at any RAC component and at any epoch exceeds the preset discrimination threshold are filtered out. The set of satellites that remain after filtering and are not marked as abnormal is the core satellite set.

[0011] Furthermore, based on the core satellite, in the process of determining the weights of each analysis center using the least squares variance component estimation method, the improved Förstner least squares variance component estimation method is first adopted to estimate the variance components of the orbital products of each analysis center based on the residual information of the orbital products of each analysis center; finally, the reciprocal of the variance component is used as the weight of the core satellite corresponding to each analysis center. Using the orbital residuals of the core satellites of each analysis center as observations, the variance component estimates obtained in the k-th iteration are calculated using the improved Förstner variance component estimation method. : , In the formula, Let r be the vector formed by the orbital residuals of the r-th core satellite in the k-th iteration. This represents the redundancy of the r-th core satellite in the k-th iteration. The weight matrix represents the weight of the r-th core satellite; To stabilize the estimation and incorporate external prior information, the prior information is integrated in the form of an inverse gamma distribution, and the above variance component estimation formula is updated as follows: , in, For the prior variance, The weights are the prior variances. Using the variance component estimates of the analysis center The reciprocal of the normalized weights of the r-th analysis center is calculated. : , Where n is the total number of analysis centers.

[0012] Furthermore, after the orbits of each analysis center are aligned by Helmert transform and outliers are removed, a weighted average is performed using the optimal weights obtained through iteration to generate the comprehensive precise orbit product for this iteration; the comprehensive orbit coordinates of a satellite at a certain epoch. Represented as: , Where n is the number of analysis centers providing the satellite orbit. Let r be the normalized weight corresponding to the r-th analysis center. Let be the satellite coordinates after alignment with the r-th analysis center.

[0013] Furthermore, the method for determining iterative convergence is as follows: a1. After completing the weight determination and orbit combination process of S4~S5, calculate the weight change of the analysis center obtained in this iteration and the previous iteration. a2. Determine whether the convergence condition is met. If the change in weight of each analysis center is less than the preset first threshold of 1 mm, or the number of iterations has reached the preset maximum number of iterations of 5, then it is determined that the iteration has converged. a3. If convergence is determined, output the integrated track product generated in the current iteration as the final integrated track product; if convergence is not achieved, proceed to the next iteration.

[0014] A computer device includes a processor and a memory, the processor being electrically connected to the memory, the memory being used to store instructions and data, and the processor being used to execute a method for integrating GNSS multi-source precision satellite orbit products.

[0015] Beneficial Effects: This invention employs least squares variance component estimation, using the actual residuals between each center and the integrated orbit as observed values, iteratively estimating their variance components, and using their reciprocals as weights. This method is entirely data-driven, automatically and objectively assigning greater weights to analysis center products with higher accuracy and greater stability, thereby mathematically optimizing the integration results and significantly improving the overall accuracy of the integrated orbit. This invention innovatively introduces the Dynamic Time Warping (DTW) algorithm to cluster the satellite residual sequences, grouping satellites with similar error patterns into the same group. Based on this, variance component estimation and weight determination are performed for each group, avoiding the weight estimation bias caused by conflating satellites with vastly different error characteristics. This ensures that the weights better reflect the true processing level of the analysis center on specific types of satellites, improving the scientific rigor and robustness of the integration method. In the variance component estimation formula, this invention innovatively incorporates prior information from satellite laser ranging (SLR) verification. When the internal residual statistics of an analysis center's product have significant uncertainty, the accuracy assessment provided by SLR, a completely independent external observation method, can serve as a strong constraint, preventing unreasonable fluctuations or even divergence in weight estimation. This fusion of internal and external verification information significantly enhances the numerical stability of the iterative optimization process and ensures the reliability of the final comprehensive product accuracy assessment. Attached Figure Description

[0016] Figure 1 This is a flowchart illustrating the GNSS multi-source precision satellite orbit product integration method of the present invention. Detailed Implementation

[0017] The embodiments of the present invention will be further described below with reference to the accompanying drawings.

[0018] like Figure 1 As shown, this invention discloses a method for integrating GNSS multi-source precision satellite orbit products, specifically including the following steps: Step 1, Data Preprocessing: Obtain multi-system precision satellite orbit data for 7 days from each analysis center and check the product integrity rate; calculate the average value of orbit data from all available analysis centers as the initial composite orbit; perform gross error removal to exclude data that differs too much from the initial composite orbit; Figure 1 In this context, AC represents the analysis center, and AC1~ACn represent analysis centers with serial numbers 1~n.

[0019] Obtain precise satellite orbit SP3 format files for multiple systems (GPS, GLONASS, Galileo, and BDS) from multiple analysis centers (IGS and iGMAS) over a period of 7 days. For each analysis center and each satellite, check the data integrity. Preprocessing of precision orbit data includes integrity checks and gross error removal, among which: Integrity check process: exclude satellites that provide data from only a single analysis center, and ensure that data can be comprehensively compared. That is, exclude satellite epochs with a data integrity rate of less than 80% in a single day to ensure the continuity of the time series. Gross error removal process: Calculate the average of all available analysis center orbit data for each satellite as the initial synthesized orbit for that satellite. Exclude satellite epoch data from the analysis center that differ from the initial integrated orbit by more than 300m; iterate this step 3 times to obtain the integrated orbit result after removing coarsely screened outliers; The initial composite orbit of the satellite is calculated using the following formula. : , in, This is the three-dimensional position vector of the satellite's orbit, where n is the number of analysis centers providing that satellite orbit. For the analysis center index, For the calendar year, Satellite number.

[0020] Step 2, Orbit Pre-alignment: If each analysis center provides rotation parameters for its solution relative to a common synthesized Earth reference frame solution (such as the IGS Repro3 daily frame synthesis solution), then these parameters are used directly to rotate and align the orbits of each analysis center to the common frame; if no such parameters are available, then all analysis center orbits are aligned with the initial synthesized orbits using the Helmert seven-parameter transformation.

[0021] If the aforementioned fine rotation parameters are unavailable, a Helmert seven-parameter transformation is used for alignment. The process of aligning all analysis center orbits with the initial integrated orbit coordinates using a Helmert transformation is as follows: For each analysis center's orbit, a set of seven parameters relative to the composite orbit is solved. The transformation model of the analysis center relative to the composite orbit is as follows: , In the formula, The initial composite orbital seven parameters to be solved are: These are the three translation components in the coordinate transformation. It consists of three rotational components. For scale parameters; Indicates the satellite's composite orbital coordinates. Represent the center orbit coordinates; establish the error equation based on the indirect adjustment model: , In the formula Represents the residual vector; It is a coefficient matrix; This is a vector of estimated values ​​for the parameter to be estimated; Let N be the vector of observed values; assuming the analysis center has N common coordinate points relative to the composite orbit, and the coefficient matrix... Estimated value vector Observation vector They are represented as follows: , Then, the coordinate transformation parameters are estimated using least squares to minimize the sum of squares of the residuals, i.e.: , Assume that the coordinates of the satellite orbit products in the error equation have the same precision. Treat the coordinates of the satellite orbit products as independent observations with equal precision. Let the weight matrix P be the identity matrix, and the estimated value vector... Represented as: , Using the obtained initial composite orbit seven parameters, the original orbit coordinates of all satellites of the analysis center at all epochs were determined according to the transformation model of the analysis center relative to the composite orbit. Perform the transformation to obtain the orbital coordinates after satellite orbit alignment. : .

[0022] Step 3: Satellite Grouping: Based on the statistical characteristics of satellites in each analysis center, mainly the mean and standard deviation of the weighted time series, a time series clustering method based on dynamic time warping is used to group the satellites to identify satellite groups with similar behavioral patterns.

[0023] The radial, tangential, and normal residual sequences of each satellite after alignment across all analysis centers are extracted to form a multidimensional time series set. If multiple analysis centers are considered, a four-dimensional feature tensor of "analysis center-satellite-direction-epoch" is constructed, and the multidimensional time series is reduced to a two-dimensional time series matrix through vectorization or principal component analysis. The Dynamic Time Warping (DTW) algorithm is used to calculate the similarity distance between any two two-dimensional time series matrices. The satellites are grouped, and the similarity distance matrix is ​​calculated using the following formula. Each element Indicates the first satellite and the first DTW distance of the satellites: , in To align the path, and These are the feature values ​​of the two two-dimensional time series at corresponding time points; and These represent the feature values ​​of the two two-dimensional time series, and These represent specific points in time. The similarity distance matrix calculated by the DTW algorithm is given by the following formula. The K-Means++ time series clustering method was used to group the satellites; the clustering objective was to minimize the squared intra-cluster distance. sum: , Where K is the preset number of clusters, and the optimal value is determined by the silhouette coefficient; This represents the k-th cluster.

[0024] Profile coefficient calculation: 1) For a single data point : calculate :point The average distance to all other points within its cluster (cohesion); calculate :point The average distance to all points in the nearest neighbor cluster (separation); 2.) Point contour coefficient for: , The range is between [-1, 1], with the closer to 1 indicating better clustering. All points The mean value is the overall contour score of the clustering scheme; by plotting the contour scores corresponding to different numbers of clusters, the number of clusters with the highest score is selected as the optimal solution.

[0025] Step 4, Outlier Removal: Calculate the positional differences of each analysis center's orbit relative to the current composite orbit in the radial, tangential, and normal directions, i.e., the orbital residuals. , i∈n, where n is the total number of analysis centers; for each core satellite, the residual sequences in the radial, tangential, and normal directions are respectively... Calculate the corrected Z-score: Use the corrected Z-score method for outlier detection and remove data that is identified as anomalous.

[0026] Using the formula: Calculate the residual sequence the median of : Calculate the absolute value of the difference between each residual and the median in the residual sequence using the following formula, and then find the median of all absolute values: , Using the formula: Calculate the corrected Z-score; where, It is a constant factor, under the assumption that the data strictly follows a normal distribution. with standard deviation The relationship is ; Examine the corrected Z-fractions for each satellite at all epochs in the radial, tangential, and normal directions. As long as the Z-fraction is corrected in any direction and at any epoch, its... If the threshold value is exceeded, the satellite is marked as abnormal and its historical data is removed.

[0027] Step 5, Core Satellite Definition: Within each satellite group, select all satellites with complete data from all analysis centers that have not been marked as abnormal in step S4 as core satellites; if a group contains only one satellite, it will be automatically regarded as a core satellite.

[0028] For each satellite group after the satellites are grouped, satellites with complete data from all analysis centers and not marked as anomalies are selected as core satellites. If a satellite group contains only one satellite, that satellite is automatically considered a core satellite. At the same time, satellites whose absolute value of the corrected Z-score exceeds the preset discrimination threshold at any RAC component and at any epoch are removed. The set of satellites remaining after the screening and not marked as anomalies is the core satellite set. The core satellite set constitutes the stable observation basis for variance component estimation.

[0029] Step 6, Weight Determination: Using the orbital residuals of the core satellites of each analysis center as observed values, the least squares variance component estimation method is used to estimate the variance components of the orbital products of each analysis center; the inverse of the variance component estimate is normalized and used as the weight of each analysis center in the synthesis.

[0030] Let the residual vector of the r-th analysis center be... Its cofactor array is Redundancy is Iterative calculations were performed using an improved Förstner variance component estimation formula: , in, It is the variance component estimate of the (k+1)th iteration; Let r be the vector formed by the orbital residuals of the r-th core satellite in the k-th iteration. This represents the redundancy of the r-th core satellite in the k-th iteration. The weight matrix represents the weight of the r-th core satellite; and These are the prior variances and their weights, which can be used to introduce external accuracy information, i.e., the SLR verification results from satellite laser ranging, to increase the stability of the estimation; if there is no reliable prior information, they can be set as follows: If , then the above expression degenerates into the standard form.

[0031] No. Weight of each analysis center It is determined by normalizing the inverse of the estimated variance components: , Where n is the total number of analysis centers; the analysis center with higher accuracy will receive a greater weight.

[0032] Step 7, Orbit Synthesis: After Helmert alignment and outlier removal, the orbits of each analysis center are weighted and averaged using the optimal weights obtained through iteration to generate the comprehensive precise orbit product for this iteration; the comprehensive orbit coordinates of a satellite at a specific epoch. Represented as: , Where n is the number of analysis centers providing the satellite orbit. Let be the normalized weight corresponding to the i-th analysis center. The coordinates of the satellite after alignment with the i-th analysis center.

[0033] Step 8, Iterative convergence judgment: Combine the orbit As a reference track for the next iteration, steps 4 to 7 are repeated to form an iterative optimization process; After each iteration, the change in the 3D root mean square error (3D RMS) of each analysis center orbital obtained in this iteration compared to the previous iteration is calculated. A convergence threshold and a maximum number of iterations are preset. If the conditions are met... If the iteration has been completed 5 times and the 6th iteration is about to be performed, the iteration is considered to have converged, and the current integrated track product is output as the final result; otherwise, the next iteration continues. Based on actual needs, non-core satellites are processed using weighted analysis by the available analysis center (AC).

[0034] Step 9, Product Accuracy Assessment and Verification: Using satellite laser ranging (SLR) data, the final integrated orbital product undergoes independent external verification to assess its accuracy and reliability. Specific steps include: b1. Data preparation: Obtain SLR observations from the same period as the integrated orbital products, as well as the precise coordinates of the SLR station under the International Earth Reference Frame (ITRF). b2. Residual calculation: Interpolate the integrated orbit product to each SLR observation time, calculate the theoretical distance from the satellite to the station, and compare it with the SLR measured distance after various corrections to obtain the SLR observation residual; b3. Statistical evaluation: Calculate the root mean square error, mean and standard deviation of the SLR residual sequence as statistical measures to evaluate the three-dimensional accuracy and systematic deviation of the integrated track product; b4. Application of Results: By analyzing the statistical characteristics of SLR residuals under different satellite types and different solar altitude angles, the accuracy, reliability and systematic error characteristics of integrated orbit products are comprehensively evaluated.

[0035] This method significantly enhances the numerical stability of the iterative optimization process by integrating internal and external verification information, thus ensuring the reliability of the final comprehensive product accuracy assessment.

[0036] The above description is merely one embodiment of the present invention and is not intended to limit the present invention. Any minor modifications, equivalent substitutions, and improvements made to the above embodiment based on the technical essence of the present invention should be included within the protection scope of the present invention.

Claims

1. A method for GNSS multi-source precise satellite orbit product synthesis, characterized in that, Includes the following steps: S1. Acquire precise orbit data from multiple analysis centers and preprocess the precise orbit data; S2. Perform preliminary pre-alignment of the center orbit based on the pre-processed precision satellite orbit data; S3. The satellites are grouped using a time series clustering method based on dynamic time warping. S4. Calculate the radial, tangential, and normal residual sequences of each analysis center orbit relative to the reference orbit, and remove outlier data based on the corrected Z-score method; S5. Define core satellites within each satellite group; use the orbital residuals of the core satellites in each analysis center as observations to determine the weights of each analysis center, introduce prior information from satellite laser ranging to enhance estimation stability, and perform a weighted average of the orbits based on the weights of each analysis center to obtain the comprehensive orbits of the core satellites. S6. Repeat S4~S5 iterative updates until convergence, to obtain the final integrated orbital product.

2. The method for integrating GNSS multi-source precision satellite orbit products according to claim 1, characterized in that, The multiple analysis centers include: the analysis center under the International GNSS Service Organization (IGS) and the analysis center under the Global Continuous Monitoring and Evaluation System (iGMAS); the satellite systems involved include: Global Positioning System (GPS), Global Navigation Satellite System (GLONASS), Galileo, and BeiDou Navigation Satellite System (BDS). Preprocessing of precision orbit data includes integrity checks and gross error removal, among which: Integrity check process: exclude satellites that provide data from only a single analysis center, and ensure that data can be comprehensively compared. That is, exclude satellite epochs with a data integrity rate of less than 80% in a single day to ensure the continuity of the time series. Gross error removal process: Calculate the average of all available analysis center orbit data for each satellite as the initial synthesized orbit for that satellite. Exclude satellite epoch data from the analysis center that differ from the initial integrated orbit by more than 300m; iterate this step 3 times to obtain the integrated orbit result after removing coarsely screened outliers; The initial composite orbit of the satellite is calculated using the following formula. : ， in, This is the three-dimensional position vector of the satellite's orbit, where n is the number of analysis centers providing that satellite orbit. For the analysis center index, For the calendar year, Satellite number.

3. The method for integrating GNSS multi-source precision satellite orbit products according to claim 2, characterized in that, Use the rotation parameters of the common synthetic Earth reference frame solution to perform the initial pre-alignment of the analysis center orbits. If the frame rotation parameters are not available, perform a Helmert transformation to align all analysis center orbits with the initial synthetic orbit coordinates. The process of aligning all analysis center orbits with the initial integrated orbit coordinates using a Helmert transformation is as follows: For each analysis center's orbit, a set of seven parameters relative to the composite orbit is solved. The transformation model of the analysis center relative to the composite orbit is as follows: ， In the formula, The initial composite orbital seven parameters to be solved are: These are the three translation components in the coordinate transformation. It consists of three rotational components. For scale parameters; Indicates the satellite's composite orbital coordinates. Indicates the coordinates of the center orbit; Establish the error equation based on the indirect adjustment model: ， In the formula Represents the residual vector; It is a coefficient matrix; This is a vector of estimated values ​​for the parameter to be estimated; Let N be the vector of observed values; assuming the analysis center has N common coordinate points relative to the composite orbit, and the coefficient matrix... Estimated value vector Observation vector They are represented as follows: ， Then, the coordinate transformation parameters are estimated using least squares to minimize the sum of squares of the residuals, i.e.: ， Assume that the coordinates of the satellite orbit products in the error equation have the same precision. Treat the coordinates of the satellite orbit products as independent observations with equal precision. Let the weight matrix P be the identity matrix, and the estimated value vector... Represented as: ， Using the obtained initial composite orbit seven parameters, the original orbit coordinates of all satellites of the analysis center at all epochs were determined according to the transformation model of the analysis center relative to the composite orbit. Perform the transformation to obtain the orbital coordinates after satellite orbit alignment. : 。 4. The method for integrating GNSS multi-source precision satellite orbit products according to claim 3, characterized in that, The satellites are grouped using time series clustering, as follows: The radial, tangential, and normal residual sequences of the orbits and composite orbits of all satellites in each satellite system are extracted. Each satellite is treated as a sequence set, forming a multidimensional time series set. If multiple analysis centers are considered, a four-dimensional feature tensor of "analysis center-satellite-direction-epoch" is constructed, and the multidimensional time series is reduced to a two-dimensional time series matrix through vectorization or principal component analysis. The Dynamic Time Warping (DTW) algorithm is used to calculate the similarity distance between any two two-dimensional time series matrices. The satellites are grouped, and the similarity distance matrix is ​​calculated using the following formula. Each element Indicates the first satellite and the first DTW distance of the satellites: ， in To align the path, and These are the feature values ​​of the two two-dimensional time series at corresponding time points; and These represent the feature values ​​of the two two-dimensional time series, and These represent specific points in time. The similarity distance matrix calculated by the DTW algorithm is given by the following formula. The K-Means++ time series clustering method was used to group the satellites; the clustering objective was to minimize the squared intra-cluster distance. sum: ， Where K is the preset number of clusters, and the optimal value is determined by the silhouette coefficient; This represents the k-th cluster.

5. The method for integrating GNSS multi-source precision satellite orbit products according to claim 4, characterized in that, The orbits of each analysis center are compared with the reference orbit, which is the initial composite orbit. For each epoch, the positional differences between the two orbits in the radial, tangential, and normal directions are calculated, i.e., the orbital residuals. , i∈n, where n is the total number of analysis centers; for each satellite, the residual sequences in the radial, tangential, and normal directions are respectively... Calculate the corrected Z-score: Calculate the residual sequence using the following formula the median of : ， Calculate the absolute value of the difference between each residual and the median in the residual sequence using the following formula, and then find the median of all absolute values: ， The corrected Z-score is calculated using the following formula: ， in, It is a constant factor; under the assumption that the data strictly follows a normal distribution, MAD is related to the standard deviation. The relationship is ; Examine the corrected Z-fractions for each satellite at all epochs in the radial, tangential, and normal directions. As long as the Z-fraction is corrected in any direction and at any epoch, its... If the threshold value is exceeded, the satellite is marked as abnormal and its historical data is removed.

6. The method for integrating GNSS multi-source precision satellite orbit products according to claim 5, characterized in that, After the satellites are grouped, select the available satellites from all analysis centers as core satellites for each group. If a satellite group contains only one satellite, that satellite is automatically considered a core satellite. At the same time, satellites whose absolute value of the corrected Z-score exceeds the preset discrimination threshold at any RAC component and at any epoch are filtered out. The set of satellites that remain after filtering and are not marked as abnormal is the core satellite set.

7. The method for integrating GNSS multi-source precision satellite orbit products according to claim 6, characterized in that, Based on the core satellite, in determining the weights of each analysis center using the least squares variance component estimation method, the improved Förstner least squares variance component estimation method is first used to estimate the variance components of the orbital products of each analysis center based on the residual information of the orbital products of each analysis center; finally, the reciprocal of the variance component is used as the weight of the core satellite corresponding to each analysis center. Using the orbital residuals of the core satellites of each analysis center as observations, the variance component estimates obtained in the k-th iteration are calculated using the improved Förstner variance component estimation method. : ， In the formula, Let r be the vector formed by the orbital residuals of the r-th core satellite in the k-th iteration. This represents the redundancy of the r-th core satellite in the k-th iteration. The weight matrix represents the weight of the r-th core satellite; To stabilize the estimation and incorporate external prior information, the prior information is integrated in the form of an inverse gamma distribution, and the above variance component estimation formula is updated as follows: ， in, For the prior variance, The weights are the prior variances. Using the variance component estimates of the analysis center The reciprocal of the normalized weights of the r-th analysis center is calculated. : ， Where n is the total number of analysis centers.

8. The method for integrating GNSS multi-source precision satellite orbit products according to claim 7, characterized in that, After alignment by Helmert transform, the center orbits of each analysis point are weighted and averaged using the optimal weights obtained through iteration to generate the comprehensive precise orbit product for that iteration; the comprehensive orbit coordinates of a satellite at a specific epoch. Represented as: ， Where n is the number of analysis centers providing the satellite orbit. Let r be the normalized weight corresponding to the r-th analysis center. Let be the satellite coordinates after alignment with the r-th analysis center.

9. The method for integrating GNSS multi-source precision satellite orbit products according to claim 8, characterized in that, The method for determining iterative convergence is as follows: a1. After completing the weight determination and orbit combination process of S4~S5, calculate the weight change of the analysis center obtained in this iteration and the previous iteration. a2. Determine whether the convergence condition is met. If the change in weight of each analysis center is less than the preset first threshold of 1 mm, or the number of iterations has reached the preset maximum number of iterations of 5, then it is determined that the iteration has converged. a3. After determining convergence, output the integrated track product generated in the current iteration as the final integrated track product; If convergence is not achieved, proceed to the next iteration.

10. A computer device, characterized in that, It includes a processor and a memory, the processor being electrically connected to the memory, the memory being used to store instructions and data, and the processor being used to execute the GNSS multi-source precision satellite orbit product integration method according to any one of claims 1-9.

Images