Advanced timing receiver autonomous integrity monitoring method, device, equipment and storage medium

Abstract

The present invention discloses an autonomous integrity monitoring method, device, equipment and storage medium for an advanced timing receiver, the method comprising: constructing a fault subset, smoothing pseudorange observation values ​​according to the carrier phase of the receiver to obtain a prior residual, using integrity support information and integrity parameters broadcast by the Beidou satellite-based augmentation system to calculate the residual weight matrix of the pseudorange, using a weighted least squares estimator to estimate the timing solution of the full-view subset composed of all visible satellites and the timing solution of all fault subsets; performing fault detection based on solution separation hypothesis detection for each fault-tolerant timing solution and a chi-square test for the full-view timing solution; in the case where the fault detection passes, calculating the timing protection level according to the integrity risk budget, and determining the availability of the timing service. The present invention eliminates the largest error source in the timing process, the ionospheric error; and adopts multi-constellation signals to greatly enhance satellite geometry and observation value redundancy.

Description

Technical Field

[0001] The present invention relates to the field of satellite navigation technology, and in particular to an advanced timing receiver autonomous integrity monitoring method, device, equipment and storage medium. Background Art

[0002] The global satellite navigation system provides users with high-precision and low-cost time and frequency sources, and is widely used in key infrastructure fields such as telecommunications networks, power transmission and the financial industry. Although the accuracy requirements for timing services in most fields are still at the millisecond level, with the rapid promotion and use of digital technologies such as smart grids, 5G and autonomous driving, the accuracy requirements for timing services will be increased to 100 nanoseconds or even higher in the future. At present, global satellite navigation systems, including the United States' Global Positioning System (GPS), China's Beidou Satellite Navigation System (BDS), the European Union's Galileo and Russia's GLONASS global navigation satellite system, have provided global users with timing services with an accuracy better than 40 nanoseconds, which is sufficient to meet the increasingly stringent accuracy requirements of more advanced applications in the future.

[0003] However, due to the inherent fragility of satellite navigation systems, such as low signal power, susceptibility to interference and terrain shielding, although the accuracy of the timing service provided by the satellite navigation system is high enough, the reliability has been difficult to guarantee, which has greatly affected the application of timing technology based on satellite navigation. For example, in the nearly 10 hours from 6:56:10 on January 26, 2016 to 16:18:40 on January 26, 2016 UTC time, due to software failure caused by the removal of retired satellites, the system time tracing data broadcast by 15 GPS satellites caused an error of about 13 microseconds in the US GPS timing, which had a great impact on global users: for example, some areas in the United States disabled the navigation equipment of emergency personnel, and the United Kingdom stopped some radio broadcasts. Therefore, the reliability monitoring of the timing service provided by the global navigation system, that is, integrity monitoring, is of great significance to the promotion and popularization of the global satellite navigation system in the field of timing.

[0004] Receiver Autonomous Integrity Monitoring (T-RAIM) for timing services is a simple and effective technique that can measure the credibility of time information provided by global satellite navigation systems. T-RAIM uses redundant observations to check the consistency of the timing solution, thereby performing fault detection and troubleshooting on the user receiver. Although there are a few studies on T-RAIM technology, these methods still have many shortcomings: 1) Poor availability: relying only on single constellation signals, poor satellite geometry, and insufficient observation redundancy, resulting in poor availability of the T-RAIM algorithm; 2) Poor accuracy: relying only on single-frequency pseudorange signals, it is impossible to fully eliminate ionospheric errors (the largest source of error in satellite navigation positioning and timing), resulting in poor timing solution accuracy; 3) High computational load: The residual-based fault detection method needs to calculate the residuals of all fault subsets at the same time, resulting in high computational complexity and poor real-time performance of the algorithm; 4) Poor scalability: It is impossible to effectively utilize the carrier and pseudorange observations of multiple systems, and the timing error caused by inter-system bias (ISB) is not considered. Summary of the invention

[0005] In order to solve the above problems existing in the prior art, the present invention provides an advanced timing receiver autonomous integrity monitoring method, device, equipment and storage medium to comprehensively enhance the accuracy and reliability of the timing service provided by the Beidou satellite navigation system. Specifically, the advanced timing receiver autonomous integrity monitoring method (T-ARAIM) proposed in the present invention uses dual-frequency signals to eliminate the largest error source in the timing process - ionospheric error; the use of multi-constellation signals greatly enhances the satellite geometry and observation value redundancy, and greatly improves the availability of the integrity monitoring algorithm. Compared with the positioning-oriented advanced receiver autonomous integrity monitoring algorithm, the scheme proposed in the present invention also improves the three key links of risk allocation method, solution estimation and protection level calculation, so that the ARAIM algorithm can focus on the integrity monitoring of the timing solution. In addition, the present invention uses the correction number and corresponding integrity information broadcast by the Beidou Satellite-Based Augmentation System (BDSBAS) to further enhance the accuracy and reliability of the Beidou timing service, which is sufficient to meet the requirements of most applications for timing services.

[0006] According to a first solution of the present invention, a method for autonomous integrity monitoring of an advanced timing receiver is provided, the method comprising:

[0007] For each epoch, based on the number of Beidou satellites observed by the receiver, a list of all possible satellite fault subsets is made, the prior probabilities of all fault subsets are calculated, and the fault subsets are arranged from small to large according to the number of simultaneous faults. For fault subsets with the same number of faults, they are arranged from large to small according to the prior probabilities, and the fault subsets to be monitored are selected in order of priority until the sum of the prior probabilities of the remaining fault subsets in the list is less than the preset threshold, thus constructing the fault subsets;

[0008] The pre-verification residual is obtained based on the carrier phase smoothed pseudorange observations of the receiver. The residual weight matrix of the pseudorange is calculated using the integrity support information and the integrity parameters broadcast by the Beidou satellite-based augmentation system. The timing solution of the full-view subset consisting of all visible satellites and the timing solution of all faulty subsets are estimated using a weighted least squares estimator.

[0009] Fault detection is performed based on solution separation hypothesis testing for each fault-tolerant timing solution and chi-square test for the fully visible timing solution;

[0010] When the fault detection passes, the timing protection level is calculated according to the integrity risk budget, and the timing protection level is compared with the set timing alarm threshold to determine the availability of the timing service.

[0011] According to a second technical solution of the present invention, there is provided an autonomous integrity monitoring device for an advanced timing receiver, the device comprising:

[0012] The fault subset determination module lists all possible satellite fault subsets for each epoch according to the number of Beidou system satellites observed by the receiver, calculates the prior probabilities of all fault subsets, and arranges the fault subsets from small to large according to the number of simultaneous faults. For fault subsets with the same number of faults, they are arranged from large to small according to the prior probabilities, and the fault subsets to be monitored are selected in order of priority until the sum of the prior probabilities of the remaining fault subsets in the list is less than the preset threshold, thus constructing the fault subsets;

[0013] The timing solution estimation module is configured to obtain the pre-verification residual by smoothing the pseudorange observation value according to the carrier phase of the receiver, calculate the residual weight matrix of the pseudorange by using the integrity support information and the integrity parameters broadcast by the Beidou satellite-based augmentation system, and estimate the timing solution of the full-view subset composed of all visible satellites and the timing solution of all faulty subsets by using a weighted least squares estimator;

[0014] A fault detection and elimination module is configured to perform fault detection based on a solution separation hypothesis test for each fault-tolerant timing solution and a chi-square test for a full-view timing solution;

[0015] The timing protection level calculation module is configured to calculate the timing protection level according to the integrity risk budget when the fault detection passes, and compare the timing protection level with the set timing alarm threshold to determine the availability of the timing service.

[0016] According to a fourth technical solution of the present invention, a non-transitory computer-readable storage medium storing instructions is provided. When the instructions are executed by a processor, the method described above is executed.

[0017] The method, device, equipment and storage medium for autonomous integrity monitoring of an advanced timing receiver according to each solution of the present invention have at least the following technical effects:

[0018] (1) Enhanced availability: The present invention integrates observation data from multiple constellations (GPS, GLONASS, Galileo, and BDS) to improve the diversity of satellite geometry and the redundancy of observation values, thereby enhancing the availability of the T-RAIM algorithm. It also solves the impact of inter-system bias on timing services when fusing and processing multi-constellation data.

[0019] (2) Improved timing accuracy and reliability: The present invention uses dual-frequency or multi-frequency pseudo-range signals to effectively eliminate ionospheric errors, thereby greatly improving the accuracy of the timing solution. The differential correction numbers and corresponding integrity information broadcast by the Beidou satellite-based augmentation system are used to enhance the accuracy of the broadcast ephemeris, and the weight matrix is ​​optimized using the integrity information, thereby improving the accuracy and reliability of the timing solution.

[0020] (3) Computational efficiency optimization: The present invention adopts a fault detection and identification method with separated multi-hypothesis testing solutions, which greatly reduces the computational load while optimizing performance. It supports the use of parallel computing technology to perform separate calculations on fault subsets, which is suitable for timing integrity monitoring of larger-scale constellations in the future. BRIEF DESCRIPTION OF THE DRAWINGS

[0021] In the drawings, which are not necessarily drawn to scale, the same reference numerals may describe similar parts in different views. The same reference numerals with letter suffixes or different letter suffixes may represent different instances of similar parts. The accompanying drawings generally illustrate various embodiments by way of example and not limitation, and together with the specification and claims, are used to illustrate the embodiments of the invention. When appropriate, the same reference numerals are used throughout the drawings to refer to the same or similar parts. Such embodiments are illustrative and are not intended to be exhaustive or exclusive embodiments of the present apparatus or method.

[0022] Figure 1 A flow chart of an advanced timing receiver autonomous integrity monitoring method according to an embodiment of the present invention is shown.

[0023] Figure 2 A structural diagram of an autonomous integrity monitoring device for an advanced timing receiver according to an embodiment of the present invention is shown. DETAILED DESCRIPTION

[0024] In order to enable those skilled in the art to better understand the technical solution of the present invention, the present invention is described in detail below in conjunction with the accompanying drawings and specific embodiments. The embodiments of the present invention are further described in detail below in conjunction with the accompanying drawings and specific embodiments, but are not intended to limit the present invention. For the various steps described herein, if there is no necessity for a causal relationship between each other, the order in which they are described as examples herein should not be regarded as a limitation, and those skilled in the art should know that they can be adjusted in order, as long as the logic between them is not destroyed, resulting in the inability to implement the entire process.

[0025] The embodiment of the present invention provides an autonomous integrity monitoring method for an advanced timing receiver. Before the method is implemented, integrity risk and continuity risk allocation can be performed first. Since the receiver clock error parameters and the receiver position parameters can be separated during estimation, and in the method proposed in the present invention, only the reliability of the timing solution is concerned without paying attention to the reliability of the positioning solution, all integrity risks and continuity risks in the present invention are allocated to the timing solution. The ISB empirical value can be used to solve the inter-system deviation between the timing solutions calculated between different systems, and is derived from the average value of the ISB of the previous day stored in the receiver. Other design parameters involved in the present invention, including integrity support information (ISM), etc., are adjusted according to different application scenarios and operation modes.

[0026] Exemplarily, this embodiment provides an application scenario of the autonomous integrity monitoring method for an advanced timing receiver. Specifically, a receiver with a known position (30.532°N, 114.357°N) is selected, and the autonomous integrity monitoring method for an advanced timing receiver proposed in the present invention is deployed. A test with a duration of 24 hours and a sampling interval of 30s is performed on March 23, 2023. Satellite ephemeris and satellite clock errors are provided by Beidou broadcast ephemeris products. Beidou satellite-based augmented correction numbers are provided by enhanced telegrams. The operating mode is selected as a single Beidou system with fixed receiver coordinates. As an example only, the deployment method of the autonomous integrity monitoring method for an advanced timing receiver proposed in the present invention in a receiver is to deploy the method proposed in the present invention in a storage medium contained in the receiver in the form of software code. The receiver generally includes a controller or a processor, and the controller or processor is used to extract relevant software code that can implement the method proposed in the present invention from the storage medium to implement the present invention.

[0027] See also Figure 1 , is a flow chart of an advanced timing receiver autonomous integrity monitoring method according to an embodiment of the present invention. The method includes the following steps 1 to 4, which are described in detail as follows.

[0028] Step 1: Fault subset construction: For each epoch, list all possible satellite fault subsets according to the number of Beidou system satellites observed by the receiver, calculate the prior probabilities of all fault subsets, and arrange the fault subsets from small to large according to the number of simultaneous faults. For fault subsets with the same number of faults, arrange them from large to small according to the prior probability, and select the fault subsets to be monitored from front to back until the sum of the prior probabilities of the remaining fault subsets in the list is less than the preset threshold.

[0029] Specifically, step 1 is to determine the fault subsets that need to be monitored and the corresponding prior probabilities. The specific implementation method is to first put all possible fault subsets into the unmonitored list, and then move the fault subsets from the unmonitored list to the monitored list one by one until the total integrity risk of the unmonitored fault subsets is lower than the preset threshold, as shown in formula (1).

[0030]

[0031] Among them, P fault,k is the prior probability of fault subset k, P THRES is a preset threshold, and k not monitored is the sequence number of the fault subset that is not monitored.

[0032] Step 2: Timing solution estimation: The carrier phase smoothed pseudorange observations of the receiver are processed by observation combination, model correction, Beidou satellite-based augmentation system differential correction number correction, etc. to obtain the pre-verification residual. The residual weight matrix of the pseudorange is calculated using the integrity support information and the integrity parameters broadcast by the Beidou satellite-based augmentation system. The weighted least squares estimator is used to estimate the timing solution of the full-view subset consisting of all visible satellites and the timing solution of all faulty subsets.

[0033] Specifically, this embodiment improves the accuracy and reliability of the timing solution by fusing multi-system signals. However, when fusing multi-system signals, ISB must be considered. The multi-system pseudorange observation equation taking ISB into account is as follows:

[0034]

[0035] in, and Represent the pseudorange observation values ​​between satellite s and receiver r in GPS system, BDS system, GLONASS system and Galileo system respectively; and Respectively represent the geometric distance between the satellite and the receiver antenna in the GPS system, BDS system, GLONASS system and Galileo system; c represents the speed of light in vacuum; Represents the receiver clock error in the GPS system including the receiver hardware delay bias; Δt s,G , Δt s ,C , Δt s,R and Δt s,E Respectively represent the satellite clock errors including satellite hardware delay bias in the GPS system, BDS system, GLONASS system and Galileo system; and They represent the slant ionospheric delay of the GPS system, BDS system, GLONASS system and Galileo system respectively; and They represent the slant tropospheric delays of the GPS system, BDS system, GLONASS system, and Galileo system respectively; and They represent the pseudorange measurement noise of GPS system, BDS system, GLONASS system and Galileo system respectively; and They represent the inter-system deviations between the BDS system, GLONASS system and Galileo system and the reference system GPS system respectively.

[0036] The BeiDou Satellite-Based Augmentation System broadcasts differential corrections and corresponding integrity information to ground users, thereby improving the positioning, navigation and timing accuracy of ground users. The embodiment of the present invention decodes the B1C / B2a navigation message broadcast by the BeiDou Satellite-Based Augmentation System at the B2a frequency point to obtain the ephemeris correction number and the satellite clock correction number, as shown in formula (3) and formula (4), respectively.

[0037]

[0038] Among them, δx(t), δy(t) and δz(t) are the ephemeris corrections on the x, y and z axes at time t, respectively, and δx, δy and δz are the reference time t. D The ephemeris correction information, and is the reference time t D The information on the right side of the equation can be obtained from 32 types of telegrams.

[0039]

[0040] Where δΔt sv (t) is the satellite clock correction at time t, δB, is the reference time t D Satellite clock correction information and satellite clock correction rate information.

[0041] The ephemeris correction and satellite clock correction provided by the Beidou Satellite-Based Augmentation System are used for correction, and the ionosphere-free combination of pseudorange observations is used to eliminate the first-order ionosphere error. The tropospheric error is corrected using an empirical model, and the following observation equation can be obtained:

[0042]

[0043] in, and Represent the prior residuals of the ionosphere-free combination of pseudorange observations of the GPS system, BDS system, GLONASS system and Galileo system respectively; and Represents the unit direction vector from the satellite to the receiver in the GPS system, BDS system, GLONASS system and Galileo system, Δx r and Represents the correction compared to the a priori receiver position and GPS system receiver clock error, ΔISB C / G , ΔISB R / G and ΔISB E / G Respectively represent the inter-system bias correction between the BDS system, GLONASS system, Galileo system and GPS system, and They represent the observation noise of the GPS system, BDS system, GLONASS system and Galileo system respectively.

[0044] The T-ARAIM algorithm proposed in the present invention can be divided into four working modes according to whether the receiver coordinates are fixed and whether the ISB is fixed, namely the first working mode, the second working mode, the third working mode and the fourth working mode. The linearization equations in different working modes are as follows:

[0045] The first working mode estimates the ISB and receiver coordinates simultaneously:

[0046]

[0047] The second working mode is to fix the ISB and estimate the receiver coordinates;

[0048]

[0049] The third working mode is to fix the receiver coordinates and estimate the ISB:

[0050]

[0051] The fourth working mode fixes the ISB and receiver coordinates at the same time:

[0052]

[0053] Among them, ΔP a , ΔP b , ΔP c and ΔP d represent the prior residuals in the first working mode, the second working mode, the third working mode and the fourth working mode respectively; G a ,G b ,G c and G d Represent the design matrices under the first working mode, the second working mode, the third working mode and the fourth working mode respectively, and Respectively represent the pre-verification residuals of the GPS system, BDS system, GLONASS system and Galileo system in the first working mode, and Respectively represent the pre-verification residuals of the GPS system, BDS system, GLONASS system and Galileo system in the second working mode, and The sub-tables represent the pre-verification residuals of the GPS system, BDS system, GLONASS system and Galileo system in the third working mode. and The sub-tables represent the pre-verification residuals of the GPS system, BDS system, GLONASS system and Galileo system in the fourth working mode.

[0054] The weighted least squares estimator is used to estimate the solution to the above linearized equation. For ease of presentation, we will Recorded as It represents the increment of GPS receiver clock error in the current epoch compared with the previous epoch. The calculation method is as follows:

[0055] S=(G T WG) -1 G T W (10)

[0056]

[0057] Among them, S represents the conversion matrix, G represents the design matrix, T represents the matrix transpose, W represents the weight matrix, S t represents the row vector corresponding to the GPS receiver clock error in the transformation matrix S, and ΔP represents the a priori residual.

[0058] The weight matrix W is obtained from the pseudo-range covariance matrix as shown in the formula. The pseudo-range covariance matrix has two forms, namely C int (error empirical model for integrity calculation) and C acc(Error empirical model for accuracy and continuity calculations), as follows:

[0059]

[0060] in, and represent the tropospheric error variance and the receiver-related error variance, respectively. represents the integrity parameters broadcast by the BeiDou Satellite-Based Augmentation System, represents the standard error model of broadcast ephemeris, and i represents the satellite number.

[0061] Therefore, according to the GPS receiver clock error at the previous moment The GPS receiver clock error at the current moment can be calculated

[0062]

[0063] The above is the estimation process of the full-view timing solution. For the fault subset, it is necessary to exclude the assumed faulty satellite from participating in the timing solution calculation, which can be quickly achieved by changing the weight matrix:

[0064]

[0065] Among them, idx k represents the satellite number contained in the kth fault subset, W (k) represents the weight matrix of the kth fault subset.

[0066] Full-view timing solution and the timing solution for the fault subset k difference as follows:

[0067]

[0068] S (k) =(G T W (k) G) -1 G T W (k) (18)

[0069] in, The transition matrix S represents the fault subset k k The row vector corresponding to the GPS receiver clock error is, The transition matrix S represents the fault subset k k The row vector corresponding to the GPS receiver clock error is, The transformation matrix S represents the full view subset k The row vector corresponding to the GPS receiver clock error.

[0070] Timing solution for fault subset k The variance of is as follows:

[0071]

[0072] The variance of is as follows:

[0073]

[0074] For the subset of faults k, the maximum possible deviation caused by the standard error model is as follows:

[0075]

[0076] Among them, b nom,i is the deviation of the standard error model of the i-th satellite, N sat Indicates the number of satellites, express The element corresponding to the i-th satellite in .

[0077] Step 3: Fault detection and troubleshooting: Fault detection includes two methods, namely, solution separation hypothesis test (SST) for each fault-tolerant timing solution and chi-square test (CST) for the full-view timing solution. In SST, the false alarm probability budget P FA_TIME Evenly distributed to each fault subset, so the monitoring threshold T of fault subset k k,t It can be calculated by the following formula:

[0078]

[0079] Among them, Q -1 (P) is the (1-p) quantile of the standard Gaussian distribution with mean 0 and variance 1, where p represents probability; N fault subsets is the total number of fault subsets, K fa,t represents the quantile point, express The variance of .

[0080] Next will and T k,t For comparison, the test is considered passed only if all fault subsets meet the conditions shown in formula (24), otherwise the SST test fails and troubleshooting must be performed.

[0081]

[0082] Among them, τ k,t represents the test quantity monitored by the solution separation hypothesis test.

[0083] In CST, construct the chi-square test

[0084] Detection threshold of CST Defined as:

[0085]

[0086] in, represents the chi-square detection threshold, P FA_CHI2 represents the false alarm probability assigned to the chi-square test, It is the cumulative distribution function of the chi-square distribution with degrees of freedom deg. The degrees of freedom for different operation modes are calculated as follows:

[0087]

[0088] Wherein, (a), (b), (c) and (d) represent the degree of freedom calculation formulas in the first working mode, the second working mode, the third working mode and the fourth working mode respectively, N sat Indicates the number of satellites, N const Indicates the constellation number.

[0089] if Greater than Then the CST test fails. When the CST test fails, if the SST of all fault subsets passes, it means that the fault occurs outside the threat model, and the process of this method should be terminated, and a warning that the timing service is unavailable should be issued. Only when both SST and CST pass can the protection level calculation in step 4 be entered, otherwise, troubleshooting must be performed.

[0090] In SST, the difference between the timing solution of all faulty subsets and the timing solution of the full-view subset is calculated, and the difference is compared with the detection threshold of SST. If all faulty subset tests can pass, the SST test passes, otherwise the test fails and troubleshooting is required. In CST, the chi-square test quantity is calculated using the residuals of all satellites and compared with the threshold of the chi-square distribution. The timing protection level can only be calculated when both CST and SST pass at the same time, otherwise it is necessary to repeat step 3 after eliminating the faulty satellite until the test passes.

[0091] Step 4: Calculation of timing protection level (TPL): The integrity risk of timing is defined as the probability that a user encounters a timing error greater than the alarm limit without triggering an alarm. The timing protection level can be calculated through the integrity risk budget, as shown in formula (29). When the TPL is less than the set timing alarm threshold (TAL), it means that the timing service is safe and reliable at this time. Otherwise, a warning that the timing service is unavailable is issued.

[0092]

[0093] Among them, TPL represents the timing protection level, Nfault subsets is the total number of fault subsets, and k is the serial number of the fault subset. represents the maximum possible deviation caused by the standard error model. represents the timing solution of the k-th fault subset. is the standard deviation of represents the time deviation of the full-view subset. represents the variance of the timing solution of the full-view subset, T k,3 represents the chi-square detection threshold related to time, P fault,k represents the prior fault probability of the k-th fault subset, PHMI TIME represents the integrity risk related to time, P fault,not monitored represents the probability of unmonitored fault subsets.

[0094] Exemplarily, using the allocated integrity risk budget, the standard error model, the difference between the timing solutions of all fault subsets and the timing solution of the full-view subset, and its variance and other parameters, the equation (29) is solved by the method of interval bisection and gradual approximation, so as to obtain the timing protection level TPL. Comparing TPL with TAL, if TPL < TAL, the timing service is reliable at this time, otherwise, the timing service is unavailable.

[0095] In this embodiment, TAL is set to 30 ns. The results of this embodiment are shown in Table 1. The receiver clock offset obtained by precise point positioning estimation is used as the true value, and the difference is made with the receiver clock offset estimated by the T-ARAIM method to obtain the timing error. It can be seen from Table 1 that when the receiver coordinates are fixed, whether the Beidou satellite-based augmentation system is used or not, the availability of T-ARAIM reaches 100%; the assistance of the Beidou satellite-based augmentation system reduces the mean and standard deviation of the timing error by about 0.5 ns and 0.3 ns respectively, and TPL also decreases by nearly 2 ns, and the accuracy and reliability of the Beidou timing service have been greatly improved.

[0096] Table 1 Performance Statistics Table

[0097]

[0098] The embodiment of the present invention also provides an advanced timing receiver autonomous integrity monitoring device. Please refer to Figure 2 , which is the structural diagram of the advanced timing receiver autonomous integrity monitoring device. The device 200 includes:

[0099] The fault subset determination module 201 lists all possible satellite fault subsets for each epoch according to the number of Beidou system satellites observed by the receiver, calculates the prior probabilities of all fault subsets, and arranges the fault subsets from small to large according to the number of simultaneous faults. For fault subsets with the same number of faults, they are arranged from large to small according to the prior probabilities, and the fault subsets to be monitored are selected in order of priority until the sum of the prior probabilities of the remaining fault subsets in the list is less than a preset threshold, thereby constructing a fault subset;

[0100] The timing solution estimation module 202 is configured to obtain the pre-verification residual according to the carrier phase smoothed pseudorange observation value of the receiver, calculate the residual weight matrix of the pseudorange using the integrity support information and the integrity parameters broadcast by the Beidou satellite-based augmentation system, and estimate the timing solution of the full-view subset composed of all visible satellites and the timing solution of all faulty subsets using a weighted least squares estimator;

[0101] A fault detection and elimination module 203 is configured to perform fault detection based on a solution separation hypothesis test for each fault-tolerant timing solution and a chi-square test for the full-view timing solution;

[0102] The timing protection level calculation module 204 is configured to calculate the timing protection level according to the integrity risk budget when the fault detection passes, and compare the timing protection level with the set timing alarm threshold to determine the availability of the timing service.

[0103] The advanced timing receiver autonomous integrity monitoring device mentioned in the embodiment of the present invention belongs to the same technical concept as the method previously described, and the technical effects they achieve are basically the same, which will not be repeated here.

[0104] An embodiment of the present invention further provides an advanced timing receiver autonomous integrity monitoring device, the device comprising:

[0105] Memory for storing computer programs;

[0106] A processor is used to execute the computer program to implement the method for autonomous integrity monitoring of an advanced timing receiver according to any embodiment of the present invention.

[0107] An embodiment of the present invention further provides a non-transitory computer-readable medium storing instructions. When the instructions are executed by a processor, the method for autonomous integrity monitoring of an advanced timing receiver according to any embodiment of the present invention is executed.

[0108] In addition, although exemplary embodiments have been described herein, the scope includes any and all embodiments based on the present invention with equivalent elements, modifications, omissions, combinations (e.g., various embodiments intersecting schemes), adaptations or changes. The elements in the claims will be interpreted broadly based on the language adopted in the claims, and are not limited to the examples described in this specification or during the implementation of this application, and the examples will be interpreted as non-exclusive. Therefore, this specification and examples are intended to be considered as examples only, and the true scope and spirit are indicated by the following claims and the full scope of their equivalents.

[0109] The above description is intended to be illustrative rather than restrictive. For example, the above examples (or one or more of them) can be used in combination with each other. For example, those of ordinary skill in the art can use other embodiments when reading the above description. In addition, in the above-mentioned specific embodiments, various features can be grouped together to simplify the present invention. This should not be interpreted as a feature of an invention that is not claimed for protection being necessary for any claim. On the contrary, the subject matter of the present invention may be less than all the features of the embodiments of a specific invention. Thus, the following claims are incorporated into the specific embodiments as examples or embodiments, wherein each claim is independently used as a separate embodiment, and it is considered that these embodiments can be combined with each other in various combinations or arrangements. The scope of the present invention should be determined with reference to the full scope of the equivalent forms of the attached claims and these claims.

Claims

1. A method for autonomous integrity monitoring of an advanced timing receiver, characterized in that: The method comprises: For each epoch, based on the number of Beidou satellites observed by the receiver, a list of all possible satellite fault subsets is made, the prior probabilities of all fault subsets are calculated, and the fault subsets are arranged from small to large according to the number of simultaneous faults. For fault subsets with the same number of faults, they are arranged from large to small according to the prior probabilities, and the fault subsets to be monitored are selected in order of priority until the sum of the prior probabilities of the remaining fault subsets in the list is less than the preset threshold, thus constructing the fault subsets; The pre-verification residual is obtained based on the carrier phase smoothed pseudorange observations of the receiver. The residual weight matrix of the pseudorange is calculated using the integrity support information and the integrity parameters broadcast by the Beidou satellite-based augmentation system. The timing solution of the full-view subset consisting of all visible satellites and the timing solution of all faulty subsets are estimated using a weighted least squares estimator. Fault detection is performed based on solution separation hypothesis testing for the timing solutions of each faulty subset and chi-square test for the timing solutions of the full-view subset; If the fault detection passes, the timing protection level is calculated according to the integrity risk budget, and the timing protection level is compared with the set timing alarm threshold to determine the availability of the timing service; The pre-verification residual is obtained based on the carrier phase smoothed pseudorange observations of the receiver. The residual weight matrix of the pseudorange is calculated using the integrity support information and the integrity parameters broadcast by the Beidou satellite-based augmentation system. The timing solution of the full-view subset consisting of all visible satellites and the timing solution of all faulty subsets are estimated using a weighted least squares estimator, including: The multi-system pseudorange observation equation considering the inter-system bias is as follows: in, and Represent the pseudorange observation values ​​between satellite s and receiver r in GPS system, BDS system, GLONASS system and Galileo system respectively; and Respectively represent the geometric distance between the satellite and the receiver antenna in the GPS system, BDS system, GLONASS system and Galileo system; c represents the speed of light in vacuum; Represents the receiver clock error in the GPS system including the receiver hardware delay bias; Δt s,G , Δt s,C , Δt s,R and Δt s,E Respectively represent the satellite clock errors including satellite hardware delay bias in the GPS system, BDS system, GLONASS system and Galileo system; and They represent the slant ionospheric delay of the GPS system, BDS system, GLONASS system and Galileo system respectively; and They represent the slant tropospheric delays of the GPS system, BDS system, GLONASS system, and Galileo system respectively; and They represent the pseudorange measurement noise of GPS system, BDS system, GLONASS system and Galileo system respectively; and They represent the inter-system deviations between the BDS system, GLONASS system, and Galileo system and the reference system GPS system; Decode the navigation message broadcast by the Beidou Satellite-Based Augmentation System at the frequency point to obtain the ephemeris correction number and satellite clock correction number, as shown in formula (3) and formula (4) respectively: Among them, δx(t), δy(t) and δz(t) are the ephemeris corrections on the x, y and z axes at time t, respectively, and δx, δy and δz are the reference time t. D The ephemeris correction information, and is the reference time t D The rate of change of ephemeris corrections; Among them, δΔt sv (t) is the satellite clock correction at time t, δB, is the reference time t D Satellite clock correction information and satellite clock correction rate information; The ephemeris correction and satellite clock correction provided by the Beidou satellite-based augmentation system are used for correction, and the ionosphere-free combination of pseudorange observations is used to eliminate the first-order ionospheric term error. The tropospheric error is corrected using an empirical model, and the following observation equation is obtained: in, and Represent the prior residuals of the ionosphere-free combination of pseudorange observations of the GPS system, BDS system, GLONASS system and Galileo system respectively; and Represents the unit direction vector from the satellite to the receiver in the GPS system, BDS system, GLONASS system and Galileo system, Δx r and Represents the correction compared to the a priori receiver position and GPS system receiver clock error, ΔISB C / G , ΔISB R / G and ΔISB E / G Respectively represent the inter-system bias correction between the BDS system, GLONASS system, Galileo system and GPS system, and Represent the observation noise of GPS system, BDS system, GLONASS system and Galileo system respectively; Four working modes are determined according to the fixing method of the receiver coordinates and the inter-system deviation, namely the first working mode, the second working mode, the third working mode and the fourth working mode, the first working mode is to simultaneously estimate the inter-system deviation and the receiver coordinates, the second working mode is to fix the inter-system deviation and estimate the receiver coordinates, the third working mode is to fix the receiver coordinates and estimate the inter-system deviation, and the fourth working mode is to simultaneously fix the estimated inter-system deviation and the receiver coordinates; The linearized equations in the first working mode, the second working mode, the third working mode and the fourth working mode are respectively shown as follows: Among them, ΔP a , ΔP b , ΔP c and ΔP d represent the prior residuals in the first working mode, the second working mode, the third working mode and the fourth working mode respectively; G a ,G b ,G c and G d Represent the design matrices under the first working mode, the second working mode, the third working mode and the fourth working mode respectively, and Respectively represent the pre-verification residuals of the GPS system, BDS system, GLONASS system and Galileo system in the first working mode, and Respectively represent the pre-verification residuals of the GPS system, BDS system, GLONASS system and Galileo system in the second working mode, and The sub-tables represent the pre-verification residuals of the GPS system, BDS system, GLONASS system and Galileo system in the third working mode. and The sub-tables represent the pre-verification residuals of the GPS system, BDS system, GLONASS system and Galileo system in the fourth working mode; Using a weighted least squares estimator, the solutions to the linearized equations in the four operating modes are estimated by the following formula: S=(G T WG) -1 G T W (10) in, for It is an abbreviation of the change in clock error compared to GPS. S represents the conversion matrix, G represents the design matrix, T represents the matrix transpose, W represents the weight matrix, and S t represents the row vector corresponding to the GPS receiver clock error in the transformation matrix S, and Δp represents the prior residual error; The weight matrix W is calculated by the following formula: in, and represent the tropospheric error variance and the receiver-related error variance, respectively. represents the integrity parameters broadcast by the BeiDou Satellite-Based Augmentation System, represents the standard error model of broadcast ephemeris, C int represents the error empirical model used for integrity calculation, C acc represents the error empirical model used for accuracy and continuity calculation, and I represents the satellite number; According to the GPS receiver clock error at the previous moment The GPS receiver clock error at the current moment can be calculated For the fault subset, the weight matrix is ​​changed by the following formula to exclude the set faulty satellites from participating in the timing solution calculation: Among them, idx k represents the satellite number contained in the kth fault subset, W (k) represents the weight matrix of the kth fault subset; Timing solution for all-view subset and the timing solution for the fault subset k difference as follows: S (k) =(G T W (k) G) -1 G T W (k) (18) in, The transition matrix S represents the fault subset k k The row vector corresponding to the GPS receiver clock error is, The transformation matrix S represents the full view subset k The row vector corresponding to the GPS receiver clock error is, Timing solution for fault subset k Variance as follows: Variance as follows: For a subset of faults k, the maximum possible deviation caused by the standard error model as follows: Among them, b nom,i is the deviation of the standard error model of the i-th satellite, N sat Indicates the number of satellites, express The element corresponding to the i-th satellite in .

2. The method according to claim 1, characterized in that The fault subsets to be monitored are selected in order until the sum of the prior probabilities of the remaining fault subsets in the list is less than the preset threshold. The calculation formula is as follows: ∑ k not monitored P fault,k ≤P THRES (1) Among them, P fault,k is the prior probability of fault subset k, P THRES is a preset threshold, and k not monitored is the sequence number of the fault subset that is not monitored.

3. The method according to claim 1, characterized in that Fault detection is performed based on the separation hypothesis detection of the timing solution for each fault subset, specifically including: The monitoring threshold T of the fault subset k is calculated by the following formula: k,t : Among them, Q -1 (P) is the 1-p quantile of the standard Gaussian distribution with mean 0 and variance 1, where p represents probability; N faultsubsets is the total number of fault subsets, P FA_TIME represents the false alarm probability budget, K fa,t represents the quantile point, express The variance of Will and T k,t For comparison, when all fault subsets satisfy the conditions shown in formula (24), the solution separation hypothesis detection of the timing solution for each fault subset passes: Among them, τ k,t represents the test quantity monitored by the solution separation hypothesis test.

4. The method according to claim 1, characterized in that: The timing protection level is calculated using the following formula: Among them, TPL represents the timing protection level, N fault subsets is the total number of fault subsets, k is the sequence number of the fault subset, represents the maximum possible deviation caused by the standard error model, represents the timing solution for fault subset k The standard deviation of represents the time deviation of the full-view subset, represents the variance of the timing solution for the full-view subset, T k,3 represents the time-dependent chi-square detection threshold, P fault,k represents the prior failure probability of the kth failure subset, PHMI TIME represents the time-related integrity risk, P fault,not monitored Indicates the probability of not monitoring the failure subset.

5. The method according to claim 4, characterized in that Comparing the timing protection level with the set timing alarm threshold to determine the availability of the timing service, specifically including: When the timing protection level is less than the set timing alarm threshold, it indicates that the timing service is available; When the timing protection level is not less than the set timing alarm threshold, it indicates that the timing service is unavailable.

6. An advanced timing receiver autonomous integrity monitoring device, characterized in that: The device comprises: The fault subset determination module lists all possible satellite fault subsets for each epoch according to the number of Beidou system satellites observed by the receiver, calculates the prior probabilities of all fault subsets, and arranges the fault subsets from small to large according to the number of simultaneous faults. For fault subsets with the same number of faults, they are arranged from large to small according to the prior probabilities, and the fault subsets to be monitored are selected in order of priority until the sum of the prior probabilities of the remaining fault subsets in the list is less than the preset threshold, thus constructing the fault subsets; The timing solution estimation module is configured to obtain the pre-verification residual by smoothing the pseudorange observation value according to the carrier phase of the receiver, calculate the residual weight matrix of the pseudorange by using the integrity support information and the integrity parameters broadcast by the Beidou satellite-based augmentation system, and estimate the timing solution of the full-view subset composed of all visible satellites and the timing solution of all faulty subsets by using a weighted least squares estimator; A fault detection and elimination module is configured to perform fault detection based on a solution separation hypothesis test for the timing solution of each all-fault subset and a chi-square test for the timing solution of the all-view subset; A timing protection level calculation module is configured to calculate a timing protection level according to an integrity risk budget when the fault detection passes, and compare the timing protection level with a set timing alarm threshold to determine the availability of the timing service; The pre-verification residual is obtained based on the carrier phase smoothed pseudorange observations of the receiver. The residual weight matrix of the pseudorange is calculated using the integrity support information and the integrity parameters broadcast by the Beidou satellite-based augmentation system. The timing solution of the full-view subset consisting of all visible satellites and the timing solution of all faulty subsets are estimated using a weighted least squares estimator, including: The multi-system pseudorange observation equation considering the inter-system bias is as follows: in, and Represent the pseudorange observation values ​​between satellite s and receiver r in GPS system, BDS system, GLONASS system and Galileo system respectively; and Respectively represent the geometric distance between the satellite and the receiver antenna in the GPS system, BDS system, GLONASS system and Galileo system; c represents the speed of light in vacuum; Represents the receiver clock error in the GPS system including the receiver hardware delay bias; Δt s,G , Δt s,C , Δt s,R and Δt s,E Respectively represent the satellite clock errors including satellite hardware delay bias in the GPS system, BDS system, GLONASS system and Galileo system; and They represent the slant ionospheric delay of the GPS system, BDS system, GLONASS system and Galileo system respectively; and They represent the slant tropospheric delays of the GPS system, BDS system, GLONASS system, and Galileo system respectively; and They represent the pseudorange measurement noise of GPS system, BDS system, GLONASS system and Galileo system respectively; and They represent the inter-system deviations between the BDS system, GLONASS system, and Galileo system and the reference system GPS system; Decode the navigation message broadcast by the Beidou Satellite-Based Augmentation System at the frequency point to obtain the ephemeris correction number and satellite clock correction number, as shown in formula (3) and formula (4) respectively: Among them, δx(t), δy(t) and δz(t) are the ephemeris corrections on the x, y and z axes at time t, respectively, and δx, δy and δz are the reference time t. D The ephemeris correction information, and is the reference time t D The rate of change of ephemeris corrections; Among them, δΔt sv (t) is the satellite clock correction at time t, δB, is the reference time t D Satellite clock correction information and satellite clock correction rate information; The ephemeris correction and satellite clock correction provided by the Beidou satellite-based augmentation system are used for correction, and the ionosphere-free combination of pseudorange observations is used to eliminate the first-order ionospheric term error. The tropospheric error is corrected using an empirical model, and the following observation equation is obtained: in, and Represent the prior residuals of the ionosphere-free combination of pseudorange observations of the GPS system, BDS system, GLONASS system and Galileo system respectively; and Represents the unit direction vector from the satellite to the receiver in the GPS system, BDS system, GLONASS system and Galileo system, Δx r and Represents the correction compared to the a priori receiver position and GPS system receiver clock error, ΔISB C / G , ΔISB R / G and ΔISB E / G Respectively represent the inter-system bias correction between the BDS system, GLONASS system, Galileo system and GPS system, and Represent the observation noise of GPS system, BDS system, GLONASS system and Galileo system respectively; Four working modes are determined according to the fixing method of the receiver coordinates and the inter-system deviation, namely the first working mode, the second working mode, the third working mode and the fourth working mode, the first working mode is to simultaneously estimate the inter-system deviation and the receiver coordinates, the second working mode is to fix the inter-system deviation and estimate the receiver coordinates, the third working mode is to fix the receiver coordinates and estimate the inter-system deviation, and the fourth working mode is to simultaneously fix the estimated inter-system deviation and the receiver coordinates; The linearized equations in the first working mode, the second working mode, the third working mode and the fourth working mode are respectively shown as follows: Among them, ΔP a , ΔP b , ΔP c and ΔP d represent the prior residuals in the first working mode, the second working mode, the third working mode and the fourth working mode respectively; G a ,G b ,G c and G d Represent the design matrices under the first working mode, the second working mode, the third working mode and the fourth working mode respectively, and Respectively represent the pre-verification residuals of the GPS system, BDS system, GLONASS system and Galileo system in the first working mode, and Respectively represent the pre-verification residuals of the GPS system, BDS system, GLONASS system and Galileo system in the second working mode, and The sub-tables represent the pre-verification residuals of the GPS system, BDS system, GLONASS system and Galileo system in the third working mode. and The sub-tables represent the pre-verification residuals of the GPS system, BDS system, GLONASS system and Galileo system in the fourth working mode; Using a weighted least squares estimator, the solutions to the linearized equations in the four operating modes are estimated by the following formula: S=(G T WG) -1 G T W (10) in, for It is an abbreviation of the change in clock error compared to GPS. S represents the conversion matrix, G represents the design matrix, T represents the matrix transpose, W represents the weight matrix, and S t represents the row vector corresponding to the GPS receiver clock error in the transformation matrix S, and ΔP represents the a priori residual; The weight matrix W is calculated by the following formula: in, and represent the tropospheric error variance and the receiver-related error variance, respectively. represents the integrity parameters broadcast by the BeiDou Satellite-Based Augmentation System, represents the standard error model of broadcast ephemeris, C int represents the error empirical model used for integrity calculation, C acc represents the error empirical model used for accuracy and continuity calculation, i represents the satellite number; According to the GPS receiver clock error at the previous moment The GPS receiver clock error at the current moment can be calculated For the fault subset, the weight matrix is ​​changed by the following formula to exclude the set faulty satellites from participating in the timing solution calculation: Among them, idx k represents the satellite number contained in the kth fault subset, W (k) represents the weight matrix of the kth fault subset; Timing solution for all-view subset and the timing solution for the fault subset k difference as follows: S (k) =(G T W (k) G) -1 G T W (k) (18) in, The transition matrix S represents the fault subset k k The row vector corresponding to the GPS receiver clock error is, The transformation matrix S represents the full view subset k The row vector corresponding to the GPS receiver clock error is, Timing solution for fault subset k Variance as follows: Variance as follows: For a subset of faults k, the maximum possible deviation caused by the standard error model as follows: Among them, b nom,i is the deviation of the standard error model of the i-th satellite, N sat Indicates the number of satellites, express The element corresponding to the i-th satellite in .

7. An advanced timing receiver autonomous integrity monitoring device, characterized in that: The device comprises: Memory for storing computer programs; A processor, configured to execute the computer program to implement the method according to any one of claims 1 to 5. 8 . A non-transitory computer-readable storage medium storing instructions, which, when executed by a processor, executes the method according to claim 1 .

Images