Method and device for modeling uncertainty of link delay calibration result and medium

Abstract

The application discloses a link time delay calibration result uncertainty modeling evaluation method and device and a medium, and relates to the technical field of satellite navigation high-precision timing. The method comprises the following steps: selecting a time laboratory station and specifying a reference station, constructing a time transfer link, solving the receiver clock error based on an ionosphere-free combination model, and inversely deducing the hardware time delay of a station to be calibrated through clock error difference; the distribution is discriminated through time delay residual statistical characteristics, a model is matched to calculate a type B uncertainty, and an empirical type B uncertainty is obtained through weighted summation; a standard deviation of a calibration time sequence is calculated to obtain a type A standard uncertainty estimation value; a synthetic expanded uncertainty is calculated by fusing the empirical type B uncertainty and the type A standard uncertainty estimation value, the calibration result reliability is evaluated, and uncertainty quantitative analysis of time delay calibration is realized. The application provides a clear and reliable evaluation framework for link calibration error, and significantly improves the reliability and interpretability of the calibration result.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of satellite navigation high-precision time service, more particularly, to a link time delay calibration result uncertainty modeling and evaluation method, device and medium. BACKGROUND

[0002] As the most widely used high-precision space-time reference acquisition means, Global Navigation Satellite System (GNSS) not only plays an important role in positioning and navigation, but also has a key function in high-precision time and frequency transfer. By receiving the high-stability time signal on the satellite, remote time synchronization between stations with a wide geographical distribution can be achieved, which is widely used in national time service system, power dispatching, financial transactions, communication network synchronization and control, scientific observation and other key fields, and has important significance for the safe operation of national infrastructure and stable support for high-tech systems. In a typical GNSS time and frequency transfer system, the receiving link is usually composed of multiple hardware components such as antennas, radio frequency front ends, intermediate frequency processing modules, signal sampling and demodulation boards, etc. Due to the influence of factors such as signal propagation path between the above link modules, circuit processing characteristics, and interface differences, a stable but difficult to directly measure systematic delay will be introduced, which is usually referred to as receiver link hardware delay. This delay is generally in the order of tens of nanoseconds to hundreds of nanoseconds, and if not corrected, it will directly affect the absolute accuracy of GNSS time transfer results, and is one of the key factors restricting the accuracy of nanosecond-level and even sub-nanosecond-level time synchronization.

[0003] In order to improve the accuracy and consistency of GNSS time transfer results, various calibration methods for receiver link hardware delay have been proposed in the prior art, which can be mainly divided into absolute calibration methods and relative calibration methods. Among them, the absolute calibration method usually relies on high-precision GNSS signal simulators, standard delay models and control environments, and calculates the absolute delay of each link module through precise measurement and modeling. Although this method has high accuracy, it requires strict equipment, complex experiments and high cost, and is difficult to be widely applied. In contrast, the relative calibration method uses the same time reference source for the calibrated receiver and the reference receiver, and inversely calculates the relative link delay of the calibrated receiver based on the differential time observation results. This method is simple to implement, has high calculation efficiency and low dependence on equipment, so it is widely used in practical engineering.

[0004] However, the existing uncertainty evaluation of the receiver chain link relative delay calibration method is still insufficient. Most of the existing researches only focus on the delay estimation itself, and lack of systematic modeling and analysis of the sources of uncertainty, error propagation mechanism and statistical characteristics. In the absence of a unified absolute reference source, calibration errors may be caused by the superposition of various factors, such as signal path disturbance, receiving environment change, equipment intrinsic noise, and station data quality. The statistical characteristics are complex and difficult to model. In addition, there is currently a lack of standardized and repeatable modeling and evaluation methods for reasonably distinguishing between A-type and B-type uncertainty sources and establishing a set of uncertainty modeling and evaluation methods suitable for GNSS receiver chain link time delay calibration results. This limits the consistency and reliability of the relative link calibration results in high-precision time and frequency transfer projects to some extent. SUMMARY

[0005] To solve the above technical problems, the present application provides a link time delay calibration result uncertainty modeling and evaluation method, device and medium, which overcomes the deficiencies of the existing GNSS receiver chain link time delay calibration method in uncertainty evaluation. The present application can quantitatively evaluate the reliability of the relative calibration results of the time transfer link, and improve the measurement reliability of the GNSS time transfer link. Specifically, the GNSS receiving chain link hardware delay calibration result is limited by experimental layout, equipment performance and data processing algorithm, and its error sources are complex and diverse. However, the existing technology lacks a systematic method for classifying and modeling these error sources and quantifying uncertainty.

[0006] The present application builds a framework combining experimental measurement error modeling, error distribution type determination, statistical fluctuation analysis and system bias evaluation, clearly distinguishes between A-type and B-type uncertainty components, and quantitatively estimates them respectively, and then synthesizes the final expanded uncertainty. This method can provide a quantitative reliability index for the link delay calibration results of the GNSS time transfer system, provide technical support for high-precision time synchronization applications, and is especially suitable for error control and quality evaluation scenarios of GNSS time and frequency transfer links in long-distance, cross-system and multi-station scenarios.

[0007] To achieve the above purpose, the present application adopts the following technical solutions:

[0008] In a first aspect, the present application provides a link time delay calibration result uncertainty modeling and evaluation method, which comprises:

[0009] Select several time laboratory stations and designate one of them as a reference station to build a time transfer link, obtain GNSS observation data of each station, and synchronously obtain precise ephemeris and satellite clock error of the corresponding period; wherein the link hardware time delay of the reference station is a known value;

[0010] The ionosphere-free combination model is used for single point positioning solution of each station to calculate the receiver clock bias sequence of each station, and on the premise that the reference station link hardware delay is known, the clock bias difference is based on the constructed time transfer link to inversely deduce the link receiver hardware delay estimation value of the to-be-calibrated station;

[0011] The receiver hardware delay estimation value of each link is differentiated with the receiver link hardware delay reference value to form an error residual sequence, statistical feature indexes of the error residual sequence are extracted, a distribution type discrimination function is constructed, according to the identified distribution type, a corresponding uncertainty estimation model is automatically matched, the B-type standard uncertainty corresponding to each link is calculated respectively, the uncertainty results of multiple links are weighted and summarized to obtain the empirical B-type uncertainty;

[0012] The time bias sequence is extracted from the time sequence after the link calibration is completed, and the standard deviation thereof is calculated as the A-type standard uncertainty estimation value, the combined expanded uncertainty is calculated according to the empirical B-type uncertainty and the A-type standard uncertainty estimation value, and the reliability of the calibration result is evaluated based on the combined expanded uncertainty.

[0013] Further, the ionosphere-free combination model is represented as:

[0014] (1)

[0015] (2)

[0016] In the formula, 、 are ionosphere-free combination pseudorange and carrier phase observation values of the station to the satellite ; 、 are frequencies of the carrier and carrier phase observation values; 、 are coordinates of the satellite and the station ; is a satellite clock bias, is a troposphere delay error, is a corresponding mapping function; is a dual-frequency ionosphere-free combination ambiguity; 、 are observation noises of the pseudorange and carrier phase ionosphere-free combination respectively; is a light speed.

[0017] Further, the receiver clock bias sequence of each station is calculated by the following formula: ​

[0018] (4)

[0019] wherein, is the receiver clock bias, is the receiver clock, is the satellite precise product time reference, denotes the coordinated universal time at the kth laboratory, denotes the receiver time compared to the hardware time delay bias.

[0020] Further, based on the constructed time transfer link, the clock bias difference is differentiated, and the link receiver hardware time delay estimate of the station to be calibrated is deduced by the following formula (5):

[0021] (5)

[0022] wherein, and denote the receiver clock biases of the reference station and the receiver to be calibrated, respectively, which are calculated by the ionosphere-free combination model; , are the coordinated universal times at the kth laboratory and the lth laboratory, respectively; , are the hardware time delay estimates of the receiver and the receiver to be calibrated, respectively, and the hardware time delay estimate of the receiver to be calibrated is calculated by formula (5) on the premise that the hardware time delay of the reference station link is known.

[0023] Further, the receiver hardware time delay estimate of each link is differentiated with the receiver link hardware time delay reference value to form the calculation process of the error residual sequence:

[0024] (6)

[0025] wherein, is the error residual sequence, is the receiver hardware time delay calibration value at the nth epoch; is the receiver link hardware time delay reference value, is the total number of epochs.

[0026] ​​​​​​​​Furthermore, statistical characteristic indicators of the error residual sequence are extracted, a distribution type discrimination function is constructed, and based on the identified distribution type, the corresponding uncertainty estimation model is automatically matched. For each link, the Type B standard uncertainty is calculated, including:

[0027] The error residual sequence is subjected to index extraction to obtain feature indices; wherein, the feature indices include skewness, kurtosis, normality test statistic, and quantile-quantile fitting plot;

[0028] Based on the aforementioned characteristic indicators, the statistical distribution type of the error residual is determined by setting a threshold; wherein, the statistical distribution type includes normal distribution, triangular distribution, and uniform distribution;

[0029] If the error residuals follow a normal distribution, then the 1st error is calculated using the following formula. j Type B standard uncertainty corresponding to each link :

[0030] (7)

[0031] In the formula, Standard deviation;

[0032] If the error residuals follow a triangular distribution, the Type B standard uncertainty of the link can be calculated using the following formula. :

[0033] (8)

[0034] In the formula, a, b, and c are the minimum, maximum, and mode of the error, respectively;

[0035] If the error residuals follow a uniform distribution, the Type B standard uncertainty of the link can be calculated using the following formula. :

[0036] (9).

[0037] Furthermore, the uncertainty results from multiple links are weighted and summed using the following formula to obtain the empirical Type B uncertainty. :

[0038] (10)

[0039] In the formula, Indicates the first Link, =1,2,…,M; M is the total number of links; This is the weighting factor.

[0040] Further, the combined expanded uncertainty is calculated by the following formula based on the empirical B-type uncertainty and the A-type standard uncertainty estimation value :

[0041] (11)

[0042] wherein, is the empirical B-type uncertainty, is the A-type standard uncertainty estimation value, is the confidence factor.

[0043] In a second aspect, the present application provides a device for modeling and evaluating the uncertainty of link delay calibration results, and the device comprises:

[0044] A link construction module is configured to select a plurality of time laboratory stations and designate one of the stations as a reference station, construct a time transfer link, obtain GNSS observation data of each station, and synchronously obtain precise ephemeris and satellite clock error corresponding to a time period; wherein the link hardware delay of the reference station is a known value.

[0045] A hardware delay estimation module is configured to perform single-point positioning solution on each station by using an ionosphere-free combination model, calculate the receiver clock error sequence of each station, and based on the known link hardware delay of the reference station, perform clock error difference based on the constructed time transfer link to inversely deduce the link receiver hardware delay estimation value of the station to be calibrated.

[0046] A B-type uncertainty calculation module is configured to perform difference based on the receiver hardware delay estimation value and the receiver link hardware delay reference value of each link to form an error residual sequence, extract statistical feature indexes of the error residual sequence, construct a distribution type discrimination function, automatically match a corresponding uncertainty estimation model according to the identified distribution type, calculate the B-type standard uncertainty of each link respectively, weight and aggregate the uncertainty results of a plurality of links to obtain the empirical B-type uncertainty.

[0047] A combined expanded uncertainty calculation module is configured to extract a time bias sequence in a time sequence in which link calibration has been completed, calculate the standard deviation thereof as an A-type standard uncertainty estimation value, calculate the combined expanded uncertainty based on the empirical B-type uncertainty and the A-type standard uncertainty estimation value, and evaluate the reliability of the calibration result based on the combined expanded uncertainty.

[0048] In a third aspect, the present application provides a readable storage medium storing one or more programs, which can be executed by one or more processors to implement the method as described above.

[0049] The present application has at least the following beneficial effects:

[0050] (1) effectively eliminate the link hardware time delay error, improve the time transfer accuracy and consistency: the hardware time delay error in the receiver link in the prior art often becomes the main bottleneck limiting the time transfer accuracy, especially in the scene of multi-system cooperative operation, the inconsistency of hardware delay will lead to systematic time deviation. The present application realizes the accurate modeling and elimination of the link time delay error of each station by constructing a relative link hardware time delay calibration model, without relying on absolute time reference source, which significantly reduces the system error and improves the consistency and accuracy of the multi-link time comparison result.

[0051] (2) build A / B type uncertainty evaluation system, enhance the reliability of calibration result: the present application systematically establishes independent evaluation process of A type and B type uncertainty, wherein the A type uncertainty is obtained by statistical fluctuation analysis of time transfer data, and the B type uncertainty is evaluated by error modeling and residual distribution identification method, and finally the complete quantitative result is given through the synthesis uncertainty and the extended uncertainty model (99.7% confidence level). The process provides a clear and reliable evaluation framework for the link calibration error, which significantly improves the reliability and interpretability of the calibration result.

[0052] (3) with adaptive modeling capability and residual distribution identification function: for different link residual data characteristics, the present application introduces an automatic distribution identification mechanism, which classifies and identifies the probability distribution type of residual error sequence through statistical characteristic indexes such as skewness, kurtosis, quantile-quantile fitting graph and normality test, and supports the adaptation of various error models such as normal, triangular and uniform. The mechanism enables the present application to have strong data adaptive ability, which can automatically select the optimal modeling method according to the observation characteristics, and improve the robustness and applicability of the overall algorithm.

[0053] (4) adapt to multi-system environment: the present application supports multi-band observation data under mainstream GNSS systems (including GPS, BDS, Galileo), is compatible with different types of GNSS receiver platforms, and has good system adaptability and expansion capability. The method can be widely used in national time center, regional / wide area time and frequency experiment network time comparison engineering, and provides stable and reliable technical support for GNSS link time delay calibration and uncertainty control, especially for high-precision time and frequency service system with high requirements for link time delay consistency and time transfer accuracy. BRIEF DESCRIPTION OF DRAWINGS

[0054] Figure 1 A flowchart of a link time delay calibration result uncertainty modeling and evaluation method according to an embodiment of the present application is shown;

[0055] Figure 2A structural diagram of a link system according to an embodiment of the present application is shown.

[0056] Figure 3 A flow chart of calculation of empirical type B uncertainty according to an embodiment of the present application is shown.

[0057] Figure 4 A structural diagram of an uncertainty modeling and evaluation device for a link delay calibration result according to an embodiment of the present application is shown. DETAILED DESCRIPTION

[0058] In order for those skilled in the art to better understand the technical solutions of the present application, the present application will be described in detail below in combination with the drawings and specific embodiments. The embodiments of the present application will be further described in detail below in combination with the drawings and specific embodiments, but not as a limitation on the present application. The order in which each step is described herein as an example should not be considered as a limitation, and those skilled in the art should know that the order can be adjusted, as long as the logic between them is not destroyed and the whole process cannot be realized.

[0059] Figure 1 A flow chart of a method for modeling and evaluating uncertainty of a link delay calibration result according to an embodiment of the present application is shown. The embodiment of the present application provides a method for modeling and evaluating uncertainty of a link delay calibration result, as shown in Figure 1 The method comprises the following steps S100-S400.

[0060] S100: Select a number of time laboratory stations, and designate one of the stations as a reference station, construct a time transfer link, obtain GNSS observation data of each station, and synchronously obtain precise ephemeris and satellite clock error of the corresponding period; wherein the link hardware time delay of the reference station is a known value.

[0061] In this embodiment, the GNSS observation data acquisition and link construction of the station can be realized by step S100. A number of time laboratory stations with high stability local time reference are selected, and one of the stations is designated as a master node to construct a time transfer link. Each station needs to continuously record multi-system, multi-frequency GNSS observation data, including pseudorange and carrier phase observation values, and the data is saved in RINEX format. Precise satellite orbit, clock error products, BIPMT series announcements and reference hardware time delay parameters are synchronously obtained to support subsequent link model construction and delay calibration.

[0062] For example, receivers deployed in several time laboratories with UTC(k) time reference are selected, and one of the receivers is selected as a reference station. Dual-frequency GNSS observation data of each station is collected, and the data format is standard RINEX format, the sampling interval is 30 seconds, and the observation time is not less than 5 consecutive days. At the same time, high-precision orbit products such as precise ephemeris and satellite clock error are obtained in the corresponding period. According to the spatial layout of the stations, a plurality of GNSS time transfer relative links are constructed, and the link hardware delay of the reference station is a known value. A constructed link system is shown in Figure 2 The link system includes a plurality of calibration links and a time transfer link, each calibration link is connected to an error modeling module and a feature judgment module, the B-type uncertainty of each calibration link is calculated based on the error modeling module and the feature judgment module, and the B-type uncertainty of each calibration link is weighted to obtain an empirical B-type uncertainty. For the time transfer link, the calibrated time transfer link can be obtained based on the standard estimation of the hardware delay, and the A-type uncertainty is calculated. Based on the empirical B-type uncertainty and the A-type uncertainty, the combined expanded uncertainty can be calculated, and the calibration and evaluation application can be performed using the combined expanded uncertainty. It should be noted that the specific calculation method as mentioned above will be described in detail in the subsequent steps.

[0063] S200: Single point positioning solution is performed on each station using the ionosphere-free combination model, and the receiver clock error sequence of each station is calculated. On the premise that the link hardware delay of the reference station is known, the clock error difference is calculated based on the constructed time transfer link, and the link receiver hardware delay estimation of the calibrated station is inversely calculated.

[0064] The relative link hardware delay calibration and error estimation can be realized by step S200.

[0065] In some embodiments, the ionosphere-free combination model is represented as:

[0066] (1)

[0067] (2)

[0068] In the formula, , are ionosphere-free combination pseudorange and carrier phase observation values of the station to the satellite ; , are the frequencies of the carrier and carrier phase observation values, respectively; , are the coordinates of the satellite and the station , respectively; , For satellite clock bias, For tropospheric delay error, For the corresponding mapping function; For dual-frequency, ionosphere-free combined ambiguity; , These are observation noises consisting of pseudorange and carrier phase without ionosphere, respectively. It is the speed of light.

[0069] In some embodiments, for a receiver connected to a time-frequency source in a time laboratory, based on PPP time-frequency transfer technology, its receiver clock bias can be expressed as:

[0070] (3)

[0071] In the formula, For receiver clock bias, For the receiver clock, This serves as the time reference for satellite precision products. When considering the hardware delay of the Timekeeping Laboratory UTC(k) receiver, it can be expressed as:

[0072] (4)

[0073] In the formula, This represents the Coordinated Universal Time (UTC) of laboratory k. Indicates the receiver time relative to Hardware delay deviation. Using receiver hardware delay calibration data published in the BIPM T Bulletin, a time laboratory was selected. receiver Using known hardware latency calibration values ​​as a reference, the time laboratory will be located... receiver to be calibrated With reference receiver A hardware delay calibration link is constructed, and PPP time transfer technology is used to obtain the time comparison results. The expression is as follows:

[0074] (5)

[0075] in and Calculated using PPP , They are respectively Laboratory and The laboratory's Coordinated World Time (UTC(k)) is published by the BIPM T Bulletin, using a receiver. With the known hardware delay calibration value as a reference, the receiver to be calibrated can be determined. Hardware latency The estimated value is compared with the corresponding reference value published by the BIPM, and the error of the current to-be-calibrated receiver hardware time delay calibration is obtained, and the calibration error and uncertainty level of the method are further evaluated.

[0076] In some embodiments, the process of implementing relative link hardware time delay calibration and error estimation can be: performing precise point positioning solution for each station using an ionosphere-free combination model, and the model parameters include station coordinates, receiver clock error, tropospheric delay, etc. The receiver clock error sequence of each station is calculated using a self-developed PPP solution module, and the calculation method is shown in formula (4). On the premise that the reference station link hardware time delay is known, the clock error difference is performed based on the constructed relative link, and the link receiver hardware time delay estimation value of the to-be-calibrated station is inversely deduced according to formula (5).

[0077] S300: Difference is performed between the receiver hardware time delay estimation value of each link and the reference value of the receiver link hardware time delay, an error residual sequence is formed, statistical feature indexes of the error residual sequence are extracted, a distribution type discrimination function is constructed, and according to the identified distribution type, a corresponding uncertainty estimation model is automatically matched, the B-type standard uncertainty corresponding to each link is calculated respectively, the uncertainty results of multiple links are weighted and summarized, and the empirical B-type uncertainty is obtained.

[0078] In this embodiment, the B-type uncertainty evaluation based on error modeling is realized by step S300. The B-type uncertainty mainly comes from non-repeated observation, empirical knowledge or external reference information, and the modeling evaluation of systematic errors or long-term change factors is very important. In the GNSS receiver relative link time delay calibration, the B-type uncertainty may involve the uncertainty of the reference time source UTC(k) itself, the station coordinate error, the multipath effect and other long-term or systematic error factors. Since the above uncertainty sources are complex and difficult to model quantitatively, the embodiment of the present application performs statistical analysis on the long-term multi-link hardware time delay calibration residual sequence, constructs an uncertainty estimation model matched with the actual error distribution type, and thus obtains the empirical B-type uncertainty evaluation result reflecting the overall stability and bias level of the method.

[0079] In some embodiments, as shown in FIG. 3, step S300 can be implemented by steps S310-S340. Figure 3

[0080] S310: Calibration error residual sequence construction.

[0081] Based on the to-be-calibrated receiver hardware time delay calibration value obtained in S200, the calibration value is compared with the reference hardware delay value published by the BIPM, and the residual sequence of the calibration error is obtained from formula (6):

[0082] ​ (6)

[0083] wherein, is the receiver hardware time delay calibration value for the th epoch; is the receiver chain hardware time delay reference value. This sequence reflects the statistical fluctuation of receiver calibration error under actual measurement conditions, providing basic data for subsequent uncertainty analysis.

[0084] S320: Error distribution type identification.

[0085] Based on the error residual sequence obtained in S310 Step S320 aims to classify the distribution pattern of the sequence to support the automatic matching and scheduling of the subsequent standard B uncertainty calculation model. To achieve the above purpose, statistical feature indicators of the residual sequence are extracted, and a distribution type discrimination function is constructed, including the following processes:

[0086] S321, statistical feature index extraction: for the residual sequence The following feature indicators are extracted: skewness (a statistical quantity that measures the symmetry of the distribution); kurtosis (measures the steepness of the tail of the distribution); normality test statistic (used to determine whether the distribution is normal); quantile-quantile fitting graph (judges the difference between the sample distribution and the target distribution). The above features are automatically extracted by the distribution identification module, supporting batch residual input and automated evaluation without human intervention.

[0087] S322, set error distribution classification rules: according to the extracted statistical indicators, the system determines the statistical distribution type that the residual may conform to by setting thresholds:

[0088] I. Normal distribution: skewness is approximately 0; kurtosis is close to 3; normality test statistic passes; quantile-quantile fitting graph is a diagonal straight line.

[0089] II. Triangular distribution: kurtosis is less than 3, boundary is obvious, distribution is symmetric "peak" shape; skewness is close to 0; normality test deviates.

[0090] III. Uniform distribution: skewness is about 0; kurtosis is less than 2.5; quantile-quantile fitting graph is a flat straight line distribution.

[0091] S323, discrimination process and output results.

[0092] System output includes: distribution type label, summary table of each statistical indicator, distribution goodness-of-fit score, and distribution comparison graph for auxiliary judgment, for subsequent automatic selection or review of estimation model.

[0093] S330: Distribution type mapping to uncertainty estimation model.

[0094] According to the distribution type identified in S320, the system automatically matches the corresponding B-type standard uncertainty estimation model. The mapping relationship is as follows:

[0095] I. Normal distribution: if the measurement value obeys the normal distribution, the expected value is , and the standard deviation is , then the B-type standard uncertainty is directly equal to the sample standard deviation, as shown in formula (7).

[0096] (7)

[0097] II. Triangular distribution: if the error obeys the triangular distribution in the interval [a, b], the distribution is determined by the minimum value a, the maximum value b, and the mode c (i.e. the position where the error is most likely to occur, satisfying a≤c≤b), and the B-type standard uncertainty is as shown in formula (8).

[0098] (8)

[0099] III. Uniform distribution: if the error obeys the uniform distribution in the interval [a, b], the B-type standard uncertainty is as shown in formula (9).

[0100] (9)

[0101] S340: Final output and report generation.

[0102] After completing the identification of the residual error distribution type and the mapping of the standard uncertainty model in S330, the system automatically performs this step to generate the final evaluation result. For each calibrated link, the system calculates the B-type standard uncertainty of the link according to the characteristics of the residual error sequence of the corresponding receiver hardware time delay error, according to the matched error distribution model , where represents the th link, =1,2,…,M. To obtain the overall empirical B-type uncertainty of the calibration method , the weighted average of the B-type uncertainty of each link is calculated, and the weight factor can be set according to the degrees of freedom of the residual error sequence, the goodness of fit, the confidence of the error distribution type, or the residual error variance. Formula (10) is as follows:

[0103] (10)

[0104] where is the weight factor of link .

[0105] S400: Based on the time deviation sequence extracted from the time series of completed link calibration, the standard deviation is calculated as the A-type standard uncertainty estimate, the combined expanded uncertainty is calculated according to the experience B-type uncertainty and the A-type standard uncertainty estimate, and the reliability of the calibration result is evaluated based on the combined expanded uncertainty.

[0106] According to the definition of A-type uncertainty by BIPM, A-type uncertainty refers to the measurement uncertainty obtained by statistical analysis of repeated observation data. In the present application, the A-type uncertainty evaluation is based on the time transfer result sequence after the completion of link calibration, reflecting the statistical fluctuation of the result. In some embodiments, S400 has the following specific steps:

[0107] S410: A-type uncertainty calculation: Extract the time deviation sequence from the time series of completed link calibration, calculate the standard deviation, and take it as the A-type standard uncertainty estimate of the link .

[0108] S420: Combined uncertainty calculation: In order to comprehensively consider the joint influence of A-type and B-type uncertainties, the internationally accepted combined uncertainty calculation method is adopted to combine the two types of uncertainties by square and root, and the standard combined uncertainty of the receiver hardware time delay calibration of the link is obtained. In order to improve the credibility level of the final uncertainty estimate result and enhance the coverage ability of the extreme observation error scenario, the expression form of the expanded uncertainty is introduced. The expanded uncertainty is weighted and enlarged by introducing the confidence factor to cover the specified confidence interval. In this embodiment, the confidence factor = 3, corresponding to a confidence level of about 99.7%, and the expanded uncertainty calculation formula is as follows:

[0109] (11)

[0110] The final combined expanded uncertainty is represented by : Under standard measurement conditions, the probability that the receiver link hardware time delay calibration error will be within this range is about 99.7%, meeting the needs of high-precision time comparison systems for measurement credibility and reliability.

[0111] Through the above implementation steps, the uncertainty modeling and quantitative calculation of the receiver hardware time delay calibration result in any GNSS relative time transfer link can be realized, with good degree of automation, system consistency and statistical robustness. This method is suitable for high-precision time and frequency application scenarios such as multi-GNSS time comparison, time transfer link uncertainty evaluation and inter-laboratory metrological calibration of international time laboratories.

[0112] The embodiment of the present application also provides a link time delay calibration result uncertainty modeling evaluation device, as shown in the figure, the device comprises: Figure 4

[0113] A link construction module 401 is configured to select a plurality of time laboratory stations and designate one of the stations as a reference station, construct a time transfer link, obtain GNSS observation data of each station, and synchronously obtain precise ephemeris and satellite clock error of a corresponding period; wherein the link hardware time delay of the reference station is a known value;

[0114] A hardware time delay estimation module 402 is configured to perform single point positioning solution on each station by using an ionosphere-free combination model, calculate a receiver clock error sequence of each station, and on the premise that the link hardware time delay of the reference station is known, perform clock error difference based on the constructed time transfer link, and inversely deduce a link receiver hardware time delay estimation value of a to-be-calibrated station;

[0115] A type B uncertainty calculation module 403 is configured to perform difference between a receiver hardware time delay estimation value and a receiver link hardware time delay reference value of each link, form an error residual sequence, extract a statistical feature index of the error residual sequence, construct a distribution type discrimination function, automatically match a corresponding uncertainty estimation model according to the identified distribution type, calculate a type B standard uncertainty corresponding to each link respectively, weight and aggregate uncertainty results of a plurality of links, and obtain an empirical type B uncertainty;

[0116] A synthetic extended uncertainty calculation module 404 is configured to extract a time bias sequence in a time sequence in which link calibration has been completed, calculate a standard deviation as a type A standard uncertainty estimation value, calculate a synthetic extended uncertainty according to the empirical type B uncertainty and the type A standard uncertainty estimation value, and evaluate the reliability of the calibration result based on the synthetic extended uncertainty.

[0117] It should be noted that the structure of the link time delay calibration result uncertainty modeling evaluation device described in the embodiment and the previously described link time delay calibration result uncertainty modeling evaluation method belong to the same technical concept, and the same beneficial effects are achieved through the same principles, which will not be described here.

[0118] The embodiment of the present application also provides a readable storage medium, the readable storage medium stores one or more programs, the one or more programs can be executed by one or more processors to implement the method described in any of the above embodiments.

[0119] ​The above embodiments are only used for illustrating the present application, and are not intended to limit the present application. Those skilled in the art can make various changes and modifications without departing from the spirit and scope of the present application. Therefore, all equivalent technical solutions belong to the scope of the present application, and the patent protection scope of the present application should be defined by the claims.

Claims

1. A method for modeling and evaluating uncertainty of link latency calibration results, the method comprising: determining a plurality of link latency calibration results; determining a plurality of link latency calibration result uncertainties; and determining a plurality of link latency calibration result uncertainty correlations. The method comprises: selecting a plurality of time laboratory stations, and specifying one of the stations as a reference station to construct a time transfer link, obtain GNSS observation data of each station, and synchronously obtain precise ephemeris and satellite clock error of a corresponding period; wherein the link hardware time delay of the reference station is a known value; single point positioning is solved for each station by using an ionosphere-free combination model, receiver clock error sequences of each station are calculated, and on the premise that the link hardware time delay of the reference station is known, clock error difference is performed based on the constructed time transfer link to inversely deduce the link receiver hardware time delay estimate value of the station to be calibrated; based on the receiver hardware time delay estimate value of each link and the receiver link hardware time delay reference value, error residual sequences are formed by difference, statistical feature indexes of the error residual sequences are extracted, a distribution type discrimination function is constructed, corresponding uncertainty estimation models are automatically matched according to the identified distribution type, B-class standard uncertainty corresponding to each link is calculated respectively, the uncertainty results of a plurality of links are weighted and summarized to obtain an empirical B-class uncertainty; a time bias sequence is extracted from a time sequence in which link calibration has been completed, the standard deviation thereof is calculated as an A-class standard uncertainty estimate value, a combined expanded uncertainty is calculated according to the empirical B-class uncertainty and the A-class standard uncertainty estimate value, and the reliability of the calibration result is evaluated based on the combined expanded uncertainty; statistical feature indexes of the error residual sequences are extracted, a distribution type discrimination function is constructed, corresponding uncertainty estimation models are automatically matched according to the identified distribution type, B-class standard uncertainty corresponding to each link is calculated respectively, including: index extraction is performed on the error residual sequences to obtain feature indexes; wherein the feature indexes include skewness, kurtosis, normality test statistics, and quantile-quantile fitting graphs; based on the feature indexes, a threshold is set to judge the statistical distribution type to which the error residual conforms; wherein the statistical distribution type includes normal distribution, triangular distribution and uniform distribution; If the error residual is subject to a normal distribution, the first j B-type standard uncertainty corresponding to the link : (7) In the formula, is the standard deviation; If the error residuals follow a triangular distribution, the type B standard uncertainty for the link is calculated by the following formula : (8) in the formula, a, b and c are respectively the minimum value, maximum value and mode of the error; If the error residual obeys uniform distribution, the B-type standard uncertainty corresponding to the link is calculated by the following formula : (9)。 2. The method of link latency calibration result uncertainty modeling evaluation according to claim 1, characterized in that, the ionosphere-free combination model is represented as: (1) (2) where, , are the ionosphere-free combined pseudorange and carrier phase observations from the satellite to the receiver; , are the frequencies of the carrier and carrier phase observations, respectively; , are the coordinates of the satellite and the receiver , respectively; is the satellite clock error, is the tropospheric delay error, is the corresponding mapping function; is the dual-frequency ionosphere-free combined ambiguity; , are the observation noises of the pseudorange and carrier phase ionosphere-free combinations, respectively; is the speed of light; is the receiver clock error.​​ 3. The method of link latency calibration result uncertainty modeling evaluation of claim 2, wherein, receiver clock error sequences of each station are calculated by the following formula: (4) wherein is the receiver clock, is the satellite precise product time reference, denotes the coordinated universal time at k laboratory, denotes the receiver time compared to the hardware time delay bias.

4. The method of link latency calibration result uncertainty modeling evaluation of claim 3, wherein, on the premise that the link hardware time delay of the reference station is known, clock error difference is performed based on the constructed time transfer link to inversely deduce the link receiver hardware time delay estimate value of the station to be calibrated by formula (5): (5) In the formula, and respectively represent the receiver clock difference of the receiver of the reference station and the receiver to be calibrated, which is calculated by the ionosphere-free combination model; , respectively represent the coordinated universal time of the laboratory and the laboratory ; , respectively represent the hardware delay estimation of the receiver and the receiver to be calibrated , and the hardware delay estimation of the receiver to be calibrated is calculated by formula (5) on the premise that the hardware delay of the reference station link is known.​ 5. The method of link latency calibration result uncertainty modeling evaluation of claim 1, wherein, the calculation process of the error residual sequences based on the difference between the receiver hardware time delay estimate value of each link and the receiver link hardware time delay reference value is: (6) wherein is the error residual sequence, is the receiver hardware time delay calibration value for the th epoch; is the receiver link hardware time delay reference value, is the total number of epochs.

6. The method of link latency calibration result uncertainty modeling evaluation of claim 1, wherein, The uncertainty results of multiple links are weighted and summarized by the following formula to obtain the empirical B-type uncertainty : (10) In the formula, represents the first link, =1,2,…,M; M is the total number of links; is a weight factor.

7. The method of link latency calibration result uncertainty modeling evaluation of claim 1, wherein, The combined expanded uncertainty is calculated from the empirical class B uncertainty and the class A standard uncertainty estimate by the formula : (11) wherein is the empirical class B uncertainty, is the class A standard uncertainty estimate, is the confidence factor.

8. A device for modeling and evaluating the uncertainty of link delay calibration results, characterized in that, The device comprises: a link construction module configured to select a plurality of time laboratory stations, and specify one of the stations as a reference station to construct a time transfer link, obtain GNSS observation data of each station, and synchronously obtain precise ephemeris and satellite clock error of a corresponding period; wherein the link hardware time delay of the reference station is a known value; The hardware delay estimation module is configured to perform single point positioning calculation on each station by using the ionosphere-free combination model, to calculate the receiver clock difference sequence of each station, and to perform clock difference difference based on the constructed time transfer link under the premise that the reference station link hardware delay is known, and to inversely deduce the link receiver hardware delay estimation value of the to-be-calibrated station; The B-type uncertainty calculation module is configured to perform difference between the receiver hardware delay estimation value and the receiver link hardware delay reference value of each link, to form an error residual sequence, to extract statistical feature indexes of the error residual sequence, to construct a distribution type discrimination function, to automatically match a corresponding uncertainty estimation model according to the identified distribution type, to calculate the B-type standard uncertainty corresponding to each link for each link, to weight and aggregate the uncertainty results of multiple links, and to obtain an empirical B-type uncertainty; The synthetic extended uncertainty calculation module is configured to extract a time bias sequence in a time sequence in which link calibration has been completed, to calculate a standard deviation as an A-type standard uncertainty estimation value, to calculate a synthetic extended uncertainty according to the empirical B-type uncertainty and the A-type standard uncertainty estimation value, and to evaluate the reliability of the calibration result based on the synthetic extended uncertainty; The method comprises the following steps: index extraction is performed on the error residual sequence to obtain feature indexes; wherein the feature indexes comprise skewness, kurtosis, normality test statistics, and quantile-quantile fitting graph; based on the feature indexes, a threshold is set to judge the statistical distribution type to which the error residual conforms; wherein the statistical distribution type comprises normal distribution, triangular distribution and uniform distribution; If the error residual is subject to a normal distribution, the first j The B-type standard uncertainty corresponding to the link : (7) In the formula, is the standard deviation; If the error residuals follow a triangular distribution, the type B standard uncertainty for the link is calculated by the following formula : (8) in the formula, a, b and c are respectively the minimum value, the maximum value and the mode of the error; If the error residual obeys uniform distribution, the B-type standard uncertainty corresponding to the link is calculated by the following formula : (9)。 9.A non-transitory computer-readable storage medium storing instructions, wherein, When the instructions are executed by the processor, the method according to any one of claims 1 to 7 is executed.

Images