GNSS observation station comprehensive evaluation and optimization method fusing subjective and objective weights

Abstract

The invention relates to the technical field of satellite navigation and positioning, in particular to a GNSS observation station comprehensive evaluation and optimization method fusing subjective and objective weights, which comprises the following steps: acquiring and preprocessing original observation signals and precision products of a GNSS observation station, and resolving by using a multi-frequency non-combined precision single-point positioning model to obtain observation station quality indexes; carrying out space division and clustering on the geographical coordinates of the observation station by utilizing a K-Means + + algorithm; constructing an analytic hierarchy model and an information entropy model based on the observation station quality indexes, respectively calculating subjective weights and objective weights, and carrying out geometric fusion to obtain comprehensive weights; and grouping according to a clustering result, calculating a close degree coefficient of each observation station by using a TOPSIS model in combination with the comprehensive weight, and outputting an optimal observation station set. According to the invention, by fusing multi-frequency UCPPP quality evaluation, spatial clustering and subjective and objective weight decision models, scientization of observation station quality evaluation and equalization of spatial distribution are realized.

Description

Technical Field

[0001] This invention relates to the field of satellite navigation and positioning technology, specifically to a comprehensive evaluation and optimization method for GNSS stations that integrates subjective and objective weights. Background Technology

[0002] The BeiDou-3 Global Navigation Satellite System has fully launched global services and introduced the new B1C / B2a signal system, which offers higher observation accuracy and better anti-interference capabilities. While a large number of IGS MGEX stations are now capable of receiving B1C / B2a signals, high-precision data processing, such as precise satellite clock bias estimation and phase deviation estimation, requires balancing computational efficiency and accuracy. This necessitates selecting high-quality, evenly distributed stations from a vast pool of data. Existing station selection techniques typically employ two strategies: one based on simple geometric distribution principles, such as selecting stations solely based on latitude and longitude grid divisions to ensure global coverage; the other based on basic data integrity, directly eliminating stations with severe data gaps. While some advanced methods attempt to incorporate signal-to-noise ratio (SNR) or multipath error as screening criteria, these indicators are often treated as independent thresholds for blanket filtering. Meanwhile, most existing high-precision GNSS precision products are still based on traditional signal generation to model and correct the signal deviations of the new system. Furthermore, the indicators provided by mainstream station quality analysis software are only applicable to general monitoring and cannot fully reflect the impact of observation quality on the generation of high-precision products. In addition, existing station selection methods rely on single indicators or subjective experience and lack a systematic multi-indicator comprehensive evaluation mechanism, making it difficult to take into account station quality, spatial distribution and data availability. To address this, a comprehensive evaluation and selection method for GNSS stations that integrates subjective and objective weights is proposed. Summary of the Invention

[0003] The purpose of this invention is to provide a comprehensive evaluation and selection method for GNSS stations that integrates subjective and objective weights, thereby achieving comprehensive evaluation and selection of GNSS stations.

[0004] To achieve the above objectives, the present invention provides the following technical solution: A comprehensive evaluation and selection method for GNSS stations that integrates subjective and objective weights includes: The raw observation signals and precision products of the GNSS station are acquired. The raw observation signals are preprocessed, including time synchronization, gross error removal and error correction. A multi-frequency non-combined precision single-point positioning model is constructed to calculate the station's geographic coordinates and pseudorange and phase residual sequences. The pseudorange and phase residual sequences, the number of visible satellites, the continuity of the observation arc, and the integrity of the multi-frequency observations are statistically analyzed to obtain the station's quality indicators. The K-Means++ algorithm is used to spatially divide the geographic coordinates of the stations, calculate and classify the geometric distances between the stations and the regional center, iteratively update the regional center, and output the spatial clustering results of the stations. Based on the station quality indicators, a hierarchical analysis model and an information entropy model are constructed to calculate subjective weights and objective weights, and a comprehensive weight is obtained by synthesizing the weights using a geometric fusion formula. Based on the spatial clustering results of the stations, the stations are grouped. Within each group, a station quality evaluation matrix is ​​constructed by combining the comprehensive weight and the station quality index. The proximity coefficient of each station is calculated using the TOPSIS model. The benchmark station in each group is selected, and the preferred station set is output.

[0005] Preferably, the multi-frequency non-combined precise single-point positioning model includes a linearized pseudorange observation equation and a linearized carrier phase observation equation; The linearized pseudorange observation equation is based on the linear superposition of the geometric distance term, receiver clock error term, tropospheric wet delay term, ionospheric delay term, and receiver inter-frequency bias term, which are the product of the unit direction cosine matrix and the station coordinate increment. The linearized carrier phase observation equation is also composed of the linear superposition of the geometric distance term, the receiver clock bias term, the tropospheric wet delay term, the ionospheric delay term, the ambiguity term based on the product of wavelength and phase ambiguity, and the inter-frequency clock bias term for the satellite system.

[0006] Preferably, the statistical analysis of the pseudorange and phase residual sequence, the number of visible satellites, the continuity of the observation arc, and the integrity of multi-frequency observations yields station quality indicators, which specifically include: Based on the pseudorange and phase residual sequence, the root mean square error of the pseudorange residual and the root mean square error of the carrier phase residual are generated by performing root mean square statistical calculation. The average number of visible satellites is generated by traversing all observation epochs within the observation period, counting the number of valid satellites and calculating the arithmetic mean. Based on the continuous observation arc segments identified by cycle slip detection, the average number of satellite arc segments is generated by counting the total number of arc segments for each satellite and calculating the mean. For different frequency channels, the number of observations at each frequency point is generated by counting the number of valid observation data.

[0007] Preferably, the specific process of spatially dividing the geographic coordinates of the station using the K-Means++ algorithm includes: The first initial cluster center is randomly selected from the geographic coordinates of the stations; subsequent cluster centers are generated by executing a roulette wheel selection mechanism whose selection probability is proportional to the minimum geometric distance between the remaining geographic coordinates of the stations and the currently selected cluster center set, until the number of cluster centers reaches a preset value. The station's geographic coordinates are calculated and the Euclidean distance between them and all cluster centers are assigned to the nearest category to generate a station geographic partition. The updated cluster centers are generated by performing an arithmetic mean operation on the station geographic coordinates within each station geographic partition until the changes in the cluster centers meet the convergence condition, thus generating the station spatial clustering result.

[0008] Preferably, the specific process of calculating the subjective weight includes: Based on the station quality indicators, a hierarchical analysis model is generated by establishing a hierarchical structure that includes a target layer, a criterion layer, and a sub-criterion layer. Based on the preset importance ranking of indicators, a judgment matrix is ​​generated by performing pairwise comparisons of the indicators in each layer of the hierarchical analysis model using the 1-9 scaling method. Based on the judgment matrix, the validity is verified by calculating the consistency ratio, and the normalized feature vector corresponding to the largest feature value is extracted to generate subjective weights.

[0009] Preferably, the objective weight calculation process is as follows: based on the station quality indicators, a normalized decision matrix is ​​generated by constructing a decision matrix and performing linear normalization processing to distinguish between positive gain indicators and negative loss indicators; based on the normalized decision matrix, an indicator entropy value is generated by applying an information entropy model to quantify the distribution uncertainty of each indicator data; based on the indicator entropy value, an objective weight is generated by calculating the difference coefficient reflecting the contribution of indicator information.

[0010] Preferably, the geometric fusion formula is synthesized as follows: Based on the subjective weight and the objective weight, the geometric feature value of the index is generated by performing a product operation on the subjective weight and the objective weight under the same evaluation index. Based on the geometric characteristic values ​​of the indicators, a comprehensive weight is generated by performing normalization calculations by calculating the ratio of the geometric characteristic value of a single indicator to the sum of the geometric characteristic values ​​of all indicators.

[0011] Preferably, the specific process of calculating the proximity coefficient of each station using the TOPSIS model includes: Based on the weighted normalized decision matrix, positive ideal solution vectors and negative ideal solution vectors are generated by extracting the best and worst performance values ​​from the column vectors of each evaluation index. For each station, calculate the positive ideal solution distance generated by the Euclidean geometric distance between the index vector and the positive ideal solution vector, and calculate the negative ideal solution distance generated by the Euclidean geometric distance between the index vector and the negative ideal solution vector; Based on the positive ideal solution distance and the negative ideal solution distance, the proximity coefficient of each station is generated by calculating the ratio of the negative ideal solution distance to the sum of the positive ideal solution distance and the negative ideal solution distance.

[0012] Compared with the prior art, the beneficial effects of the present invention are as follows: 1. This application constructs a multi-frequency non-combined precise single-point positioning model to directly extract quality indicators that are highly correlated with the performance of precision products, such as pseudorange and phase residual RMS, and satellite arc number, from the solution process. This overcomes the shortcomings of existing quality analysis software, which is only suitable for general monitoring and cannot reflect the actual impact of observations on high-precision estimation. It can more realistically and accurately characterize the performance of the station when performing high-precision tasks such as satellite clock error and phase deviation estimation, and fills the gap in the lack of effective quality assessment tools for new system signals such as Beidou-3 B1C / B2a.

[0013] 2. This application introduces the K-Means++ clustering algorithm to perform intelligent spatial partitioning of global stations based on geographic coordinates, ensuring the representativeness and balance of the preferred station set in terms of geometric structure. Compared with traditional selection strategies based on simple latitude and longitude grid partitioning or manual experience, it effectively avoids the waste of computing power caused by excessive density of local stations and the information loss caused by coverage blind spots. While ensuring global service coverage, it significantly improves the computational efficiency and stability of the network model by eliminating redundant stations.

[0014] 3. This application establishes a hierarchical entropy weight decision model that integrates the Analytic Hierarchy Process (AHP) and the entropy weight method, and combines it with the TOPSIS method for comprehensive evaluation. This achieves a scientific integration of subjective expert experience and objective data patterns, and solves the problems of one-sidedness and instability of results caused by existing station selection methods that rely on single index thresholds or purely subjective experience judgments. It constructs an automated, interpretable and robust station selection system, thus providing a high-quality data foundation for high-precision GNSS data processing. Attached Figure Description

[0015] Figure 1 A flowchart illustrating a comprehensive evaluation and selection method for GNSS stations that integrates subjective and objective weights; Figure 2 A flowchart for extracting observation quality features of GNSS stations driven by multi-frequency non-combined precise single-point positioning; Figure 3 This is a flowchart of GNSS station spatial clustering based on geographic feature constraints; Figure 4 Flowchart for the optimal selection of GNSS stations by integrating subjective and objective weights; Figure 5 This is a schematic diagram of the hierarchical structure of a three-layer hierarchical analysis model. Detailed Implementation

[0016] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0017] Example 1: Please see Figure 1 This invention provides a comprehensive evaluation and optimization method for GNSS stations that integrates subjective and objective weights. The technical solution is as follows: The raw observation signals and precision products of the GNSS station are acquired. The raw observation signals are preprocessed, including time synchronization, gross error removal and error correction. A multi-frequency non-combined precision single-point positioning model is constructed to calculate the station's geographic coordinates and pseudorange and phase residual sequences. Statistical analysis is performed on the pseudorange and phase residual sequences, the number of visible satellites, the continuity of the observation arc, and the completeness of multi-frequency observations to obtain the station's quality indicators. The K-Means++ algorithm is used to spatially divide the geographic coordinates of the stations, calculate and classify the geometric distances between the stations and the regional center, iteratively update the regional center, and output the spatial clustering results of the stations. Based on the station quality indicators, a hierarchical analysis model and an information entropy model are constructed to calculate subjective weights and objective weights, and a comprehensive weight is obtained by synthesizing the weights using a geometric fusion formula. Based on the spatial clustering results of the stations, the stations are grouped. Within each group, a station quality evaluation matrix is ​​constructed by combining the comprehensive weight and the station quality index. The proximity coefficient of each station is calculated using the TOPSIS model. Benchmark stations within each group are selected, and a set of preferred stations is output.

[0018] See Figure 2 The process involves acquiring observation files from all GNSS stations capable of multi-frequency and multi-system observations for the selected day, along with corresponding precise ephemeris and precise clock bias products. The raw observation data is then synchronized with the time data and gross errors are removed to ensure the integrity and consistency of the observations. Satellites with missing dual-frequency pseudorange and carrier phase observations are deleted, and the pseudorange observations and signal-to-noise ratio (SNR) at each observation epoch are checked to ensure they are within a reasonable range. If abnormal pseudorange observations or excessively low SNRs are detected, the corresponding observations are automatically removed to avoid interference from low-quality data in subsequent results. To comprehensively characterize the observation characteristics and error distribution patterns at various frequencies across different systems, multiple combinations of dual-frequency observations are defined for each satellite, primarily including pseudo range difference (PI), carrier phase difference (LI), and the Melbourne-Wübbena combination (MW): ; In the formula, , , These represent combinations of PI, LI, and MW, respectively. , , , These represent the observed values ​​of the pseudorange at the first and second frequencies, and the phase at the first and second frequencies, respectively. and These represent the first frequency point and the second frequency point, respectively.

[0019] Cycle slip detection was performed using a combination of MW and LI observations. The MW combination, by fusing pseudorange and phase observations, can sensitively detect discontinuous changes in integer ambiguity, while the LI combination, by analyzing phase difference changes between different frequency points, supplements the detection of weak cycle slips and anomalous ionospheric disturbances. Simultaneously, the continuity of observation arcs for each satellite was determined and tracked, recording the number of continuous observation arcs for each satellite throughout the entire observation epoch.

[0020] Error corrections are applied to the retained high-quality observations to ensure the consistency of accuracy in subsequent multi-frequency UCPPP calculations. These corrections include code bias, BDS-2 satellite multipath error, tropospheric delay, gravity delay error, antenna phase entanglement, antenna phase center deviation, and tidal corrections, ensuring consistency between observations from different systems and frequencies.

[0021] The multi-frequency non-combined precise single-point positioning model includes a linearized pseudorange observation equation and a linearized carrier phase observation equation. The linearized pseudorange observation equation is based on the linear superposition of the geometric distance term, receiver clock error term, tropospheric delay term, ionospheric delay term, and receiver inter-frequency bias term, which are the product of the unit direction cosine matrix and the station coordinate increment. The linearized carrier phase observation equation is also composed of the linear superposition of the geometric distance term, the receiver clock bias term, the tropospheric delay term, the ionospheric delay term, the ambiguity term based on the product of wavelength and phase ambiguity, and the inter-frequency clock bias term for a specific satellite system.

[0022] Before constructing the equations, it is necessary to determine the set of parameters to be estimated, including station coordinate increments, receiver clock errors, tropospheric delay, ionospheric parameters, inter-frequency bias, and ambiguity parameters. By establishing a linear relationship between the observed values ​​and these parameters, the extended Kalman filter is used for state estimation. Specifically, the corresponding linearized pseudo-range With phase The observation equation is: ; In the formula, Represents the cosine matrix in the unit direction; Indicates the increment of station coordinates; This represents the tropospheric delay mapping function; Represents the ionospheric coefficient matrix; Indicates inter-frequency deviation (IFB); and These represent the wavelength and phase ambiguity at different frequencies, respectively. This indicates inter-frequency clock bias, which exists only in GPS, GLONASS, and BDS-2. and These represent the noise in the code and phase observations, respectively.

[0023] The parameter estimation using extended Kalman filtering results in a state vector comprising the following five components: station three-dimensional coordinate increment, receiver clock error, tropospheric wet delay, ionospheric delay, inter-frequency deviation of each system, and uncombined ambiguity. The process noise covariance matrix in the filtering process is set according to the physical characteristics of the parameters: the static station coordinates are considered constant parameters, and their process noise is set to 0; the receiver clock error is not modeled for time correlation, but is regarded as an independent white noise parameter between epochs, and is estimated by setting a large initial variance; the tropospheric delay is treated as a slowly varying parameter, and the process noise is set to 3 × 10⁻⁶ per second. -8 Square meters; the ionospheric delay parameter is explicitly estimated in UCPPP, and its time variation is relatively rapid, with the process noise set to 9 × 10⁻⁶ per second. -4 Square meters; inter-frequency deviation is considered as a random walk process, and the process noise is set to 1×10⁻⁶ per second. -8 square meters; The measurement noise covariance matrix is ​​constructed based on an elevation angle weighting model. For pseudorange observations, the standard deviation is set to 0.3 meters multiplied by a weighting factor equal to the reciprocal of the elevation angle sine plus one; for carrier phase observations, the standard deviation is set to 0.003 meters multiplied by a weighting factor equal to the reciprocal of the elevation angle sine plus one. Through the above parameter configuration, the convergence and solution accuracy of the filter when processing multi-frequency data are ensured. The multi-frequency non-combination model avoids the noise amplification effect introduced by traditional linear combination. By refining the modeling of inter-frequency deviation and inter-frequency clock deviation, the hardware delay and system error of the station are effectively separated, ensuring that the pseudorange and phase residual sequences can truly and purely reflect the original signal quality of the station, thus providing a high-confidence data foundation for the comprehensive evaluation of the station.

[0024] The statistical analysis of the pseudorange and phase residual sequences, the number of visible satellites, the continuity of the observation arc, and the completeness of multi-frequency observations yields station quality indicators, which specifically include: Based on the pseudorange and phase residual sequence, the root mean square error of the pseudorange residual and the root mean square error of the carrier phase residual are generated by performing root mean square statistical calculation. The average number of visible satellites is generated by traversing all observation epochs within the observation period, counting the number of valid satellites and calculating the arithmetic mean. Based on the continuous observation arc segments identified by cycle slip detection, the average number of satellite arc segments is generated by counting the total number of arc segments for each satellite and calculating the mean. For different frequency channels, the number of observations at each frequency point is generated by counting the number of valid observation data.

[0025] Specifically, after obtaining the parameters through Kalman filtering, the observed residuals of pseudorange and phase at each frequency are calculated. The formula is: ; In the formula, Represents the original observation values; These are the theoretical observations calculated after parameter estimation.

[0026] The entire epoch and satellites were traversed, and the observation quality indicators of each station were statistically analyzed. Specifically, this included: calculating the number of consecutive observation arcs for each satellite and its average value; obtaining the number of observations for each satellite at each frequency and their average value; and further calculating the average number of visible satellites for each station throughout the entire observation period. Simultaneously, the pseudorange and phase residual RMS were statistically analyzed to obtain the average RMS.

[0027] The root mean square error of the pseudorange residual and the root mean square error of the carrier phase residual are obtained by acquiring the pseudorange residual sequence and phase residual sequence of each station at each frequency point throughout the entire observation period. By averaging the sum of squared residuals over all valid epochs, the pseudorange residual and the carrier phase residual are calculated respectively, which directly reflect the degree to which the station's observation values ​​are affected by multipath effects, observation noise, and unmodeled errors. The average number of visible satellites is obtained by identifying and counting the number of effective satellites participating in the calculation at each observation epoch during the observation period, summing the number of effective satellites in all epochs, and dividing by the total number of epochs to obtain the average number of visible satellites at the station, which reflects the visibility conditions and obstruction situation of the station in its geographical location. The average number of satellite arc segments is based on the cycle slip detection results in the preprocessing stage. It identifies the continuous observation arc segments of each satellite during the observation period. First, the total number of continuous arc segments of each satellite is counted. Then, the arithmetic mean of the number of arc segments of all observed satellites is calculated to generate the average number of satellite arc segments, which is used to characterize the continuity and stability of the observation data. The fewer the number of arc segments (meaning the more continuous the data), the higher the quality of the station. The number of observations at each frequency point is calculated for different frequency channels. After preprocessing, the actual number of valid observation data points at that frequency point is counted. The number of observations for each satellite at each frequency point and its average value are further calculated, reflecting the tracking capability of the station receiver for different signals and the data integrity rate.

[0028] The above indicators are systematically integrated to generate station quality assessment documents, providing quantitative input for subsequent comprehensive evaluation. The format refers to the SINEX document standard published by IGS, and systematically records the multi-dimensional observation quality indicators of each station, including the average number of visible satellites, observation residual RMS, arc continuity, and the number of observations at each frequency point. Through the standardized output quality assessment documents, the traceability of station performance and cross-system comparison can be achieved, providing highly consistent and quantifiable basic data support for subsequent station selection.

[0029] See Figure 3 The specific process of spatially dividing the geographic coordinates of the station using the K-Means++ algorithm includes: The first initial cluster center is randomly selected from the geographic coordinates of the stations; subsequent cluster centers are generated by executing a roulette wheel selection mechanism whose selection probability is proportional to the minimum geometric distance between the remaining geographic coordinates of the stations and the currently selected cluster center set, until the number of cluster centers reaches a preset value. The station's geographic coordinates are calculated and the Euclidean distance between them and all cluster centers are assigned to the nearest category to generate a station geographic partition. The updated cluster centers are generated by performing an arithmetic mean operation on the station geographic coordinates within each station geographic partition until the changes in the cluster centers meet the convergence condition, thus generating the station spatial clustering result.

[0030] Specifically, let the set of geographic coordinates of all stations be: ; in, Indicates the first The three-dimensional coordinates of each station, Let represent the desired number of stations to be selected. First, randomly select one station from the set of stations as the initial cluster center. For each of the remaining stations, calculate the minimum Euclidean distance between it and the set of currently selected cluster centers, reflecting the station's independence from existing cluster centers in geographic space.

[0031] The probability of a station being selected as a new center is determined by the squared distance between the station and the existing center: ; in, Indicates the first The probability that a test station is selected as the next cluster center Indicates the first The minimum Euclidean distance between each station and the currently selected set of cluster centers; This indicates the total number of stations. This represents the sum of the squares of the minimum Euclidean distances between all stations; A roulette wheel selection mechanism is introduced to adaptively select new cluster centers based on the above probabilities, ensuring that the geographical distribution of the initial centers has good representativeness and diversity. This process is repeated until a new cluster center is selected. The initial cluster centers.

[0032] The distance between each station and all cluster centers is calculated, and each station is assigned to the cluster corresponding to the nearest center, thus achieving initial geographical partitioning. Subsequently, the arithmetic mean of the coordinates of all stations within each cluster is calculated to obtain new cluster centers. : ; In the formula, Indicates the first The number of cluster internal testing stations.

[0033] Repeat the calculation of cluster centers until the coordinate change of all cluster centers in two consecutive iterations is less than a preset threshold of 1×10. -6 The iterations can continue until the maximum number of iterations, 100, is reached. This method optimizes the spatial distribution of GNSS stations globally, avoiding regional information redundancy caused by station concentration, thereby improving the stability and representativeness of subsequent parameter estimation and multi-system joint solutions.

[0034] The convergence condition is specifically defined as follows: Calculate the change in Euclidean distance between the spatial locations of all cluster centers in two adjacent iterations; when the maximum change is less than 1 × 10⁻⁶, the convergence condition is satisfied. -6 The algorithm is considered converged when the number of iterations reaches a preset limit (e.g., 100 times).

[0035] By introducing the K-Means++ algorithm and establishing a direct proportional relationship between the probability of a station being selected and the square of the minimum geometric distance, the initial cluster centers are ensured to have maximum discrete distribution and diversity in geographic space. This effectively overcomes the defect of traditional random initialization being prone to getting trapped in local optima. This intelligent partitioning mechanism based on geographic feature constraints can achieve a balanced layout of stations globally, effectively avoiding regional information redundancy and computing power waste caused by excessive station density, while also preventing the generation of coverage blind spots. This significantly improves the geometric strength and computational stability of the subsequent precise positioning network solution model.

[0036] See Figure 4 The specific process for calculating the subjective weight includes: Based on the station quality indicators, a hierarchical analysis structure model is generated by establishing a hierarchical structure including a target layer, a criterion layer, and a sub-criterion layer. A judgment matrix is ​​generated by performing pairwise comparisons of the indicators at each level of the hierarchical analysis structure model using the 1-9 scaling method. Based on the judgment matrix, validity verification is performed by calculating the consistency ratio, and the normalized eigenvector corresponding to the largest eigenvalue is extracted to generate subjective weights.

[0037] Specifically, the described model is a three-layer hierarchical analysis structure, such as... Figure 5 As shown, the target layer selects the optimal GNSS station for clock error estimation; the criterion layer contains multiple station quality indicators, including the average number of satellites, pseudorange residuals, phase residuals, number of satellite arcs, and the number of observations at each frequency point; the sub-criterion layer further refines to specific systems and frequencies, such as the number of BDS-3 satellites, the RMS of pseudorange and phase residuals of C1 / C5 / L1 / L5, the number of BDS-3 satellite arcs, and the number of observations for B1C and B2a; The 1-9 scale method compares indicators pairwise, using numbers from 1 to 9 and their reciprocals to represent relative importance. Here, 1 indicates that the two indicators are equally important; 3 indicates that the former is slightly more important than the latter; 5 indicates that the former is significantly more important than the latter; 7 indicates that the former is strongly more important than the latter; 9 indicates that the former is extremely more important than the latter; and 2, 4, 6, and 8 represent the median values ​​of these adjacent judgments.

[0038] See Figure 5A criterion-level judgment matrix was constructed using the analytic hierarchy process (AHP). A 1-9 scale was used for pairwise comparisons of the indicators. Taking clock error estimation as an example, considering that it mainly relies on the ionospherically unaffected combination of pseudorange and carrier phase observations, and that carrier phase has higher accuracy and a more significant impact on clock error stability, the phase residual RMS was designated as the most important indicator. Although pseudorange residuals have lower accuracy, they also play a crucial role in clock error estimation, thus their importance is secondary. The number of satellite arcs reflects signal continuity, and the average number of satellites represents spatial geometry and visibility; these two indicators have similar impacts and are assigned equal weights. In contrast, the number of observations at each frequency point has a more indirect impact on clock error estimation, therefore its importance is the lowest. The resulting order of indicator importance is as follows: Phase residual RMS > pseudorange residual RMS > number of satellite averages = number of satellite arcs > number of observations at each frequency. The corresponding comparison matrix is: ; in, This represents the importance of each factor. .

[0039] Perform a consistency check on the following formulas and calculate the consistency ratio to verify the reasonableness of the judgment: ; In the formula, Indicators of consistency Represents the largest eigenvalue. This indicates the number of indicators in the criteria layer. , Indicates the consistency ratio. This represents the average random consistency index.

[0040] When the consistency ratio is less than 0.1, the judgment matrix is ​​considered to meet the consistency requirements. When the test passes, the normalized eigenvalue corresponding to the largest eigenvalue is taken as the result of the subjective weight of the analytic hierarchy process. By employing the analytic hierarchy process (AHP) to construct a hierarchical structure model and using the 1-9 scaling method to quantitatively describe expert experience and actual observation needs, this invention effectively solves the problem of excessive reliance on subjective arbitrariness in the allocation of indicator weights in traditional methods. Particularly in clock error estimation applications, this invention assigns higher weights to high-precision indicators such as phase residuals based on physical characteristics, and ensures the rigor and scientific nature of the judgment matrix logic through consistency checks. This approach allows the generated subjective weights to accurately reflect the actual contribution of each quality indicator to the high-precision solution, significantly improving the relevance and reliability of the station optimization results.

[0041] The objective weight calculation process is as follows: based on the station quality indicators, a normalized decision matrix is ​​generated by constructing a decision matrix and distinguishing between positive gain indicators and negative loss indicators through linear normalization; based on the normalized decision matrix, an indicator entropy value is generated by applying an information entropy model to quantify the distribution uncertainty of each indicator data; based on the indicator entropy value, an objective weight is generated by calculating the difference coefficient reflecting the contribution of indicator information.

[0042] Specifically, a decision matrix is ​​constructed based on the station quality indicators and the candidate stations obtained under each cluster: ; In the formula, Indicates the first The first monitoring station was at the Performance values ​​under each quality indicator The number of candidate stations, The quantity of quality indicators; To eliminate the differences in the dimensions of different indicators, linear normalization is performed on each quality indicator. For positive indicators, namely the number of satellites and the number of observations, gain normalization is performed, while for negative indicators, namely the observation residuals and the number of arc segments, inverse normalization is performed, so that the values ​​of each indicator are mapped to the interval [0, 1], ensuring comparability between different latitudes.

[0043] The information entropy model is chosen to reflect the uncertainty of each indicator and its contribution to the overall variability. For the first... The entropy value of each quality indicator is calculated as follows: in ; In the formula, The normalized index value, This is the adjustment coefficient.

[0044] Based on the information entropy of the calculated quality indicators, the weight of each quality indicator is calculated. The smaller the entropy value, the more significant the distribution difference of the indicator, the greater the amount of information, and the higher the weight accordingly.

[0045] By introducing an information entropy model to calculate objective weights, the potential for human bias and empirical limitations in subjective weighting is effectively overcome. Differentiated linear normalization of positive and negative indicators successfully eliminates dimensional differences between multidimensional quality indicators (such as the number of satellites and the RMS residual of observations), ensuring data comparability. Simultaneously, by quantifying the dispersion and uncertainty of each indicator data using entropy values, higher-discriminating indicators can be automatically identified and assigned greater weight based on the data's inherent distribution characteristics. This objective weighting mechanism fully leverages the practical value inherent in the observation data, ensuring the objectivity and scientific rigor of the station optimization results and their accurate reflection of the actual data state.

[0046] The geometric fusion formula is synthesized as follows: Based on the subjective weight and the objective weight, the geometric feature value of the index is generated by performing a product operation on the subjective weight and the objective weight under the same evaluation index. Based on the geometric characteristic values ​​of the indicators, a comprehensive weight is generated by performing normalization calculations by calculating the ratio of the geometric characteristic value of a single indicator to the sum of the geometric characteristic values ​​of all indicators.

[0047] Based on the aforementioned objective weights, subjective and objective information are integrated through geometric fusion. The geometric fusion formula is as follows: ; in, Indicates the first The overall weight after integrating the quality indicators For the first Subjective weighting of each quality indicator For objective weighting, This represents the total number of quality indicators involved in the evaluation, with the denominator being the normalization factor. A geometric fusion model is adopted to scientifically integrate subjective and objective weights, effectively avoiding the one-sidedness of single weighting methods. By constructing a product-form fusion mechanism, this method not only takes into account the guiding role of expert experience in key physical indicators and the information entropy characteristics of the observation data itself, but also strengthens the weights of indicators with high consistency between subjective and objective judgments through the multiplicative effect, effectively suppressing extreme biases from a single perspective. This dynamic balancing strategy achieves self-correction and optimization of weight allocation, significantly improving the adaptability, robustness, and reliability of the comprehensive evaluation system in complex multi-frequency, multi-system GNSS environments.

[0048] The specific process of calculating the proximity coefficient of each station using the TOPSIS model includes: Based on the weighted normalized decision matrix, positive ideal solution vectors and negative ideal solution vectors are generated by extracting the best and worst performance values ​​from the column vectors of each evaluation index. For each station, the positive ideal solution distance is generated by calculating the Euclidean geometric distance between the index vector and the positive ideal solution vector, and the negative ideal solution distance is generated by calculating the Euclidean geometric distance between the index vector and the negative ideal solution vector. Based on the positive ideal solution distance and the negative ideal solution distance, the proximity coefficient of each station is generated by calculating the proportion of the negative ideal solution distance to the sum of the positive ideal solution distance and the negative ideal solution distance. Specifically, a weighted normalized decision matrix is ​​constructed. Based on the normalized station quality index values ​​obtained in the preceding steps and the calculated comprehensive weights, the normalized index values ​​of each station are multiplied by their corresponding comprehensive weights to generate weighted index data. In the weighted normalized decision matrix, the optimal performance value under each column, i.e. each evaluation index, is extracted to form a positive ideal solution vector, which represents the theoretically optimal station state; at the same time, the worst performance value of each column is extracted to form a negative ideal solution vector, which represents the theoretically worst station state. For each GNSS station to be evaluated, the Euclidean distance between its weighted index vector and the positive ideal solution vector is calculated, denoted as the positive ideal solution distance; simultaneously, the Euclidean distance between it and the negative ideal solution vector is calculated, denoted as the negative ideal solution distance. This distance is obtained by calculating the arithmetic square root of the sum of the squares of the corresponding index differences; For each station, calculate the ratio of its negative ideal solution distance to the sum of the positive and negative ideal solution distances. This ratio is the proximity coefficient, i.e.: ; In the formula, For the positive ideal solution distance, The distance is the negative ideal solution. This is the proximity coefficient; a larger proximity coefficient value indicates that the overall performance of the station is closer to the ideal state. Through analysis of... By sorting the values, an optimal set of stations can be obtained, providing a high-quality, spatiotemporally balanced observation benchmark for satellite clock bias estimation.

[0049] In summary, this invention employs multi-frequency non-combined precise single-point positioning to extract station quality indicators, and uses the K-Means++ clustering algorithm and hierarchical entropy weight decision model to achieve intelligent optimization of GNSS stations. This not only improves the robustness and objectivity of the station optimization results, but also enhances the overall observation quality and spatial balance of the station set in the network solution model, such as satellite clock error and phase deviation estimation. This provides reliable technical support for high-precision GNSS time and frequency transmission and precise positioning.

[0050] Example 2: This embodiment takes the precise estimation of satellite clock bias for new signals of the BeiDou-3 global system as an example. It selects a preferred set of stations with a target number from a massive number of candidate stations worldwide. The specific implementation process is as follows: The process involves data acquisition and index extraction based on non-combined precise point positioning. Multi-frequency observation data and precision products from candidate stations are collected daily. After preprocessing the raw data, a multi-frequency non-combined precise point positioning model is constructed and solved. During the solution process, quality indicators for each station are directly generated statistically. Taking one candidate station as an example, its calculated average number of visible satellites throughout the day is 5, the RMS statistical value of B1C pseudorange residual is 0.27 meters, the RMS statistical value of B1C carrier phase residual is 0.003 meters, and the number of valid observations at a specific new signal frequency B1C is 612. Through this step, evaluation files containing detailed quality data for all candidate stations are generated in batches, realizing the digital quantification of the station's physical performance. The K-Means++ algorithm is used to perform spatial clustering based on geographic features. The geographic coordinates of all candidate stations are used as input data, and the number of cluster centers is set according to a preset selection criteria. In the initial stage of the algorithm, the first cluster center is randomly selected. Subsequently, based on the principle that the greater the geometric distance between other stations and the current center, the higher the probability of selection, the remaining initial centers are selected sequentially to ensure the global discreteness of the initial distribution. After multiple rounds of iteration, the Euclidean distance from each station to each center is calculated and the center position is updated until the algorithm converges. Finally, the global candidate stations are divided into multiple non-overlapping geographic clusters, ensuring that each geographic region contains several candidate stations and avoiding overly dense or sparse stations in some areas. The comprehensive weighting of subjective experience and objective data was calculated. Based on the analytic hierarchy process (AHP), considering that the accuracy of carrier phase data plays a decisive role in the final product quality in clock error estimation scenarios, the importance of the carrier phase residual RMS index was assigned to the highest level, and the subjective weight of carrier phase residual RMS was calculated to be 0.25. Simultaneously, the uncertainty of the distribution of actual observation data was quantified using an information entropy model to calculate the objective weight reflecting the inherent differences in the data. Subsequently, the two were integrated using a geometric fusion model. After normalization, the comprehensive weight of the B1C carrier phase residual RMS index was adjusted to 0.19, and the comprehensive weight of the B2a carrier phase residual RMS index was adjusted to 0.18, while the weights of less influential auxiliary indicators were correspondingly reduced, thus generating a scoring standard that conforms to both the physical mechanism and the current data situation. The TOPSIS model is used to optimize station selection within each cluster. Taking a geographical group containing multiple candidate stations as an example, a weighted normalized decision matrix is ​​constructed by combining comprehensive weights, and the positive ideal solution vector (representing the theoretical optimum) and negative ideal solution vector (representing the theoretical worst value) are extracted within the group. Calculations show that the third station in the group has the closest Euclidean distance to the positive ideal solution, with a proximity coefficient of 0.83, while the lowest proximity coefficient in the group is 0.16. Based on the principle of maximizing the proximity coefficient, the first station is selected as the representative station for this geographical region. The above steps are performed for all geographical groups, ultimately outputting a uniformly distributed and optimal set of stations.

[0051] Although embodiments of the invention have been shown and described, it will be understood by those skilled in the art that various changes, modifications, substitutions and alterations can be made to these embodiments without departing from the principles and spirit of the invention, the scope of which is defined by the appended claims and their equivalents.

Claims

1. A comprehensive evaluation and optimization method for GNSS stations that integrates subjective and objective weights, characterized in that, include: The system acquires raw observation signals and precise products from GNSS stations. It performs time synchronization, gross error removal, cycle slip detection, and error correction preprocessing on the raw observation signals. Through multi-frequency non-combined precise single-point positioning calculations, it obtains the station's geographic coordinates, pseudorange and phase residual sequences, number of visible satellites, number of consecutive observation arcs for each satellite, and number of multi-frequency observations. A comprehensive statistical analysis is then performed on the pseudorange and phase residual sequences, the number of visible satellites, the number of consecutive observation arcs for each satellite, and the number of multi-frequency observations to construct station quality indicators, including: root mean square error of pseudorange and phase residuals, average number of visible satellites for each satellite system, average number of satellite arcs, and number of observations for each frequency point. The K-Means++ algorithm is used to spatially divide the geographic coordinates of the stations, calculate and classify the geometric distances between the stations and the regional center, iteratively update the regional center, and output the spatial clustering results of the stations. Based on the station quality indicators, a hierarchical analysis model and an information entropy model are constructed to calculate subjective weights and objective weights, and a comprehensive weight is obtained by synthesizing the weights using a geometric fusion formula. Based on the spatial clustering results of the stations, the stations are grouped. Within each group, a station quality evaluation matrix is ​​constructed by combining the comprehensive weight and the station quality index. The proximity coefficient of each station is calculated using the TOPSIS model. The benchmark station in each group is selected, and the preferred station set is output.

2. The method for comprehensive evaluation and selection of GNSS stations integrating subjective and objective weights as described in claim 1, characterized in that, The multi-frequency non-combined precise single-point positioning model includes a linearized pseudorange observation equation and a linearized carrier phase observation equation. The linearized pseudorange observation equation is based on the linear superposition of the geometric distance term, receiver clock error term, tropospheric wet delay term, ionospheric delay term, and receiver inter-frequency bias term, which are the product of the unit direction cosine matrix and the station coordinate increment. The linearized carrier phase observation equation is also composed of the linear superposition of the geometric distance term, the receiver clock bias term, the tropospheric wet delay term, the ionospheric delay term, the ambiguity term based on the product of wavelength and phase ambiguity, and the inter-frequency clock bias term for the satellite system.

3. The method for comprehensive evaluation and optimization of GNSS stations integrating subjective and objective weights as described in claim 1, characterized in that, The statistical analysis of the pseudorange and phase residual sequences, the number of visible satellites, the continuity of the observation arc, and the completeness of multi-frequency observations yields station quality indicators, which specifically include: Based on the pseudorange and phase residual sequence, the root mean square error of the pseudorange residual and the root mean square error of the carrier phase residual are generated by performing root mean square statistical calculation. The average number of visible satellites is generated by traversing all observation epochs within the observation period, counting the number of valid satellites and calculating the arithmetic mean. Based on the continuous observation arc segments identified by cycle slip detection, the average number of satellite arc segments is generated by counting the total number of arc segments for each satellite and calculating the mean. For different frequency channels, the number of observations at each frequency point is generated by counting the number of valid observation data.

4. The GNSS station comprehensive evaluation and optimization method integrating subjective and objective weights as described in claim 1, characterized in that, The specific process of spatially dividing the geographic coordinates of the station using the K-Means++ algorithm includes: The first initial cluster center is randomly selected from the geographic coordinates of the stations; subsequent cluster centers are generated by executing a roulette wheel selection mechanism whose selection probability is proportional to the minimum geometric distance between the remaining geographic coordinates of the stations and the currently selected cluster center set, until the number of cluster centers reaches a preset value. The station's geographic coordinates are calculated and the Euclidean distance between them and all cluster centers are assigned to the nearest category to generate a station geographic partition. The updated cluster centers are generated by performing an arithmetic mean operation on the station geographic coordinates within each station geographic partition until the changes in the cluster centers meet the convergence condition, thus generating the station spatial clustering result.

5. The method for comprehensive evaluation and optimization of GNSS stations integrating subjective and objective weights as described in claim 1, characterized in that, The specific process for calculating the subjective weight includes: Based on the station quality indicators, a hierarchical analysis model is generated by establishing a hierarchical structure that includes a target layer, a criterion layer, and a sub-criterion layer. The judgment matrix is ​​generated by performing pairwise comparisons of the indicators of each layer in the hierarchical analysis model using the 1-9 scaling method; based on the judgment matrix, the validity is verified by calculating the consistency ratio, and the normalized feature vector corresponding to the largest eigenvalue is extracted to generate subjective weights.

6. The method for comprehensive evaluation and optimization of GNSS stations integrating subjective and objective weights as described in claim 1, characterized in that, The objective weight calculation process is based on the station quality index, and generates a normalized decision matrix by constructing a decision matrix and performing linear normalization processing to distinguish between positive gain index and negative loss index. Based on the normalized decision matrix, the index entropy value is generated by quantifying the distribution uncertainty of each index data using the information entropy model. Based on the index entropy value, the objective weight is generated by calculating the difference coefficient reflecting the contribution of the index information.

7. The GNSS station comprehensive evaluation and optimization method integrating subjective and objective weights as described in claim 1, characterized in that, The geometric fusion formula is synthesized as follows: Based on the subjective weight and the objective weight, the geometric feature value of the index is generated by performing a product operation on the subjective weight and the objective weight under the same evaluation index. Based on the geometric characteristic values ​​of the indicators, a comprehensive weight is generated by performing normalization calculations by calculating the ratio of the geometric characteristic value of a single indicator to the sum of the geometric characteristic values ​​of all indicators.

8. The method for comprehensive evaluation and optimization of GNSS stations integrating subjective and objective weights as described in claim 1, characterized in that, The specific process of calculating the proximity coefficient of each station using the TOPSIS model includes: Based on the weighted normalized decision matrix, positive ideal solution vectors and negative ideal solution vectors are generated by extracting the best and worst performance values ​​from the column vectors of each evaluation index. For each station, calculate the positive ideal solution distance generated by the Euclidean geometric distance between the index vector and the positive ideal solution vector, and calculate the negative ideal solution distance generated by the Euclidean geometric distance between the index vector and the negative ideal solution vector; Based on the positive ideal solution distance and the negative ideal solution distance, the proximity coefficient of each station is generated by calculating the ratio of the negative ideal solution distance to the sum of the positive ideal solution distance and the negative ideal solution distance.

Images