Real-time GIM product integrity monitoring method based on non-difference and non-combination

Abstract

The invention relates to a real-time GIM product integrity monitoring method based on non-difference and non-combination. The method comprises the following steps: estimating a first vertical ionosphere delay at an ionosphere puncture point by using a non-difference and non-combination model; receiving the real-time GIM product and obtaining a second vertical ionospheric delay at the ionospheric grid point; selecting at least three first vertical ionospheric delays near ionospheric grid points, and calculating to obtain third vertical ionospheric delays; subtracting the second vertical ionosphere delay from the third vertical ionosphere delay to construct a detection statistic; calculating a detection threshold according to the false alarm rate, comparing the detection statistics of the grid points with the detection threshold to carry out fault detection, calculating the minimum detectable deviation of the ranging domain of the grid points passing the fault detection, and comparing the minimum detectable deviation of the grid points with a preset alarm limit. The grid points with the minimum detectable deviation smaller than the alarm limit are normal grid points, that is, the grid points pass the integrity monitoring of the GIM product.

Description

Technical Field

[0001] The present invention relates to the field of satellite navigation technology, and in particular to a real-time GIM product integrity monitoring method based on non-differential and non-combination. Background Art

[0002] The Global Navigation Satellite System (GNSS) has powerful functions such as navigation, positioning, and timing. It can provide all-round, all-weather, and all-time navigation and positioning services, and can be widely used in various military and civilian vehicles such as land, sea, and air. The ionosphere is a key risk source affecting the atmospheric segment of GNSS space, directly affecting the propagation direction and speed of radio waves, thereby generating ionospheric delay errors. The total electron content of the ionosphere (TEC) is an important parameter in ionospheric research. In order to establish a complete ionospheric model and provide a reliable global ionospheric TEC grid (Global Ionospheric Maps, GIM), IGS established the Ionosphere Working Group in 1998 to approve and implement the ionospheric total electron content data exchange format file IONEX (Ionosphere Exchange Format), with a product accuracy of 2 to 8 TECU.

[0003] GIM products use data from long-term reference stations distributed around the world, derive TEC of each grid point through spherical harmonic functions and algorithms, and generate real-time, fast, final and forecast products, providing rich data resources for global ionosphere research and applications. GIM products are not only basic data products for inverting the laws of ionosphere changes and revealing the laws of ionosphere evolution, but also one of the important ways to enhance the positioning performance of Precise Points Positioning (PPP).

[0004] Factors such as time, season, solar activity and geomagnetic latitude will affect the electron density of the ionosphere, making the ionosphere concentration show strong regularity accompanied by random inhomogeneity, affecting the ionosphere observation information; at the same time, the GNSS observation stations used to extract ionosphere parameters are unevenly distributed around the world, resulting in a lack of effective observation data in large areas, so the accuracy of the model in this area is significantly reduced. Although the quality of GIM products has gradually improved, the reliability of GIM products is difficult to guarantee due to the above factors. Therefore, the ionospheric delay error corrected by unreliable GIM products will affect the accuracy of GNSS services and even threaten its integrity, limiting the application of GIM products in life safety-related fields.

[0005] The traditional GIM product integrity monitoring method relies on differentially acquired ionospheric information, and analyzes whether the RMS map (Root mean square map) can reliably characterize the corresponding TEC error characteristics from the perspective of the rationality of ionospheric model accuracy. However, the existing GIM product integrity monitoring method does not strictly control the false alarms and missed detection errors of GIM products, which makes it impossible to effectively ensure the integrity of GIM products. Summary of the invention

[0006] In order to solve or improve the above problems, the present invention proposes a real-time GIM product integrity monitoring method based on non-differential and non-combination, comprising:

[0007] Step S1) collecting raw observation data from global integrity monitoring stations and receiving real-time satellite clock and orbit correction products, and estimating the first vertical ionospheric delay at the ionospheric puncture point between the monitoring station and the satellite using a non-differenced non-combined model;

[0008] Step S2) receiving the real-time GIM product and obtaining the second vertical ionospheric delay at the ionospheric grid point;

[0009] Step S3) selecting at least three first vertical ionospheric delays obtained in step S1 near the ionospheric grid point, and calculating the third vertical ionospheric delay at the ionospheric grid point by inverse distance weighted interpolation method;

[0010] Step S4) subtracting the second vertical ionospheric delay from the third vertical ionospheric delay, thereby constructing a detection statistic for real-time GIM product integrity monitoring;

[0011] Step S5) Calculate the detection threshold according to the false alarm rate required for abnormal integrity monitoring of real-time GIM products, traverse each ionospheric grid point and compare its detection statistic and the detection threshold to perform fault detection, if the detection statistic is less than or equal to the detection threshold, execute step S6; otherwise, if the detection statistic is greater than the detection threshold, mark the corresponding ionospheric grid point as abnormal;

[0012] Step S6) calculates the minimum detectable deviation of the ranging domain of the ionospheric grid point that has passed the fault detection, traverses the minimum detectable deviation at each ionospheric grid point and compares it with the preset alarm limit. If the minimum detectable deviation is less than the preset alarm limit, the corresponding ionospheric grid point is normal; otherwise, if the minimum detectable deviation is greater than or equal to the preset alarm limit, the corresponding ionospheric grid point fails the detection, the ionospheric grid point is marked as abnormal, and an alarm is issued to the user.

[0013] The beneficial effects of the present invention are:

[0014] The present invention provides a method for monitoring the integrity of real-time GIM products based on non-differential and non-combined methods. The method takes the vertical ionospheric delay at the grid points estimated by non-differential and non-combined methods as a reference, constructs a detection statistic by subtracting the vertical ionospheric delay at the grid points corresponding to the real-time GIM products, constructs a detection threshold for the required false alarm rate, and performs fault detection. Finally, according to the missed detection rate control requirements, the minimum detectable deviation in the ranging domain of the GIM product grid points is calculated, and compared with the alarm limit to achieve missed detection error control. Through the dual control of false alarm errors and missed detection errors, abnormal real-time GIM products can be effectively detected, which effectively ensures the integrity of real-time GIM products, and is of great significance in the application of navigation and positioning in fields related to life safety. BRIEF DESCRIPTION OF THE DRAWINGS

[0015] Figure 1 The present invention provides a flow chart of a method for real-time GIM product integrity monitoring based on non-differential and non-combination. DETAILED DESCRIPTION

[0016] The technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the drawings in the embodiments of the present invention. All other embodiments obtained by ordinary technicians in the field without making creative work based on the embodiments of the present invention shall fall within the scope of protection of the present invention.

[0017] The non-differential and non-combination-based real-time GIM product integrity monitoring method proposed in the present invention can achieve dual control of false alarm errors and missed detection errors of GIM products, effectively monitor abnormal real-time GIM products, and fully ensure the integrity of GIM products. The technical solution of the present invention is further specifically described below in conjunction with the accompanying drawings.

[0018] like Figure 1 As shown, the present invention provides a real-time GIM product integrity monitoring method based on non-differential and non-combination, comprising the following steps:

[0019] Step S1) collects raw observation data from global integrity monitoring stations and receives real-time satellite clock and orbit correction products, and estimates the first vertical ionospheric delay at the ionospheric puncture point (IPP) between the monitoring station and the satellite using a non-differenced non-combined model.

[0020] Collect raw observation data from global integrity monitoring stations and receive real-time satellite clock and orbit correction products. In the non-differenced, non-combined PPP model, the oblique ionospheric delay is retained in the observation equation as a parameter to be estimated. When performing PPP solution, it is necessary to first use the precise satellite clock error product to correct the clock error. The satellite clock error in the precise satellite clock error product absorbs the satellite hardware delay bias, while the receiver clock error absorbs the receiver hardware delay bias. In the non-differenced, non-combined PPP model, the satellite hardware delay bias and the receiver hardware delay bias cannot be offset in the pseudorange observation value, but can be absorbed by the ionospheric parameters. The satellite hardware delay bias and the receiver hardware delay bias can be absorbed by the ionospheric parameters and the periodic ambiguity parameters in the carrier phase observation value. The specific expression is as follows:

[0021]

[0022] In formula (1), the superscript s represents the satellite and the subscript r represents the receiver. It represents the ionospheric parameters absorbing the hardware delay, where the hardware delay includes the satellite hardware delay and the receiver hardware delay; Represents the pure ionospheric parameters; is the ionosphere-free combination coefficient, γ j (j=1,2) represent the ionospheric mapping factor at frequency j, which can be expressed as DCB r Indicates the differential code deviation at the receiver, DCB r =b r,1 -b r,2 , b r,j (j=1,2) represents the hardware delay of pseudorange at the receiver end; DCB s Indicates the satellite differential code deviation, represents the hardware delay of the pseudorange at the satellite end; j represents the wavelength corresponding to frequency j; represents the ambiguity parameter absorbing the hardware delay, where the hardware delay includes the satellite hardware delay and the receiver hardware delay; represents the pure ambiguity parameter corresponding to frequency j; B r,j Indicates the hardware delay of the carrier phase at the receiver; Indicates the hardware delay of the satellite carrier phase; t r,IF represents the receiver clock error of the ionosphere-free combination, Represents the satellite clock error of the ionosphere-free combination.

[0023] Since the receiver hardware delay bias values ​​of each base station are the same, they can be eliminated through inter-satellite differentials without introducing redundant model error terms. Converted to the first vertical ionospheric delay Iest .

[0024]

[0025] In formula (2), I est is the first vertical ionospheric delay, Re is the radius of the earth, α is the satellite altitude angle, h is the distance between the ionospheric shell model and the earth's surface, It represents the ionospheric parameters absorbing the hardware delay, which includes the satellite hardware delay and the receiver hardware delay.

[0026] Step S2) Receive the real-time GIM product and obtain the second vertical ionospheric delay at the ionospheric grid point (IGP).

[0027] The TEC (total electron content of the ionosphere) at the ionospheric grid point corresponding to the real-time GIM product is calculated by the constant 40.3 / f 2 Convert to second vertical ionospheric delay:

[0028]

[0029] In formula (3), I GIM_grid Represents the second vertical ionospheric delay at the ionospheric grid point; TEC grid represents the total electron content at the ionospheric grid point, and f represents the frequency.

[0030] Step S3) m (m≥3) first vertical ionospheric delays obtained in step S1 are selected near the ionospheric grid point, and the third vertical ionospheric delay at the ionospheric grid point is calculated by inverse distance weighted interpolation method.

[0031] If there is no ionospheric piercing point between the monitoring station and the satellite in the vicinity of some individual grid points, the status of the grid point is marked as not monitored.

[0032] In the satellite navigation system, a large number of global ground observation stations will be established. Among them, a large number of global ground observation stations are selected as monitoring stations to collect original observation data and use them to solve the non-difference and non-combined ionospheric delay. Among them, a few monitoring stations that are not used to calculate the ionospheric delay are used as stations to be solved. The ionospheric delay of the stations to be solved is obtained through the ionospheric delay of the monitoring stations near the corresponding grid points.

[0033] Select m monitoring stations around the station to be solved, m ≥ 3, and record the first vertical ionospheric delay obtained in step S1 as I est_i , calculate the distance d from each of the m monitoring stations to the station to be solved i , the third vertical ionospheric delay at the ionospheric grid point is obtained by the inverse proportional weighted interpolation method:

[0034]

[0035] In formula (4), I est_grid is the third vertical ionospheric delay at the ionospheric grid point, I est_i is the first vertical ionospheric delay calculated by the i-th monitoring station, d i is the distance from the ith monitoring station to the station to be solved, m is the number of monitoring stations selected around the station to be solved, and m≥3.

[0036] Step S4) Subtract the second vertical ionospheric delay obtained in step S2 from the third vertical ionospheric delay obtained in step S3, thereby constructing a detection statistic for real-time GIM product integrity monitoring.

[0037] Based on the vertical ionospheric delay obtained by non-differenced and non-combined estimation, a detection statistic res is constructed for real-time GIM product integrity monitoring. grid , which is expressed as:

[0038] res grid =I GIM_grid -I est_gr id (5)

[0039] In formula (5), res grid represents the detection statistic used for real-time GIM product integrity monitoring, I GIM_grid represents the second vertical ionospheric delay at the ionospheric grid point, I est_grid Represents the third vertical ionospheric delay at the ionospheric grid point.

[0040] The RMS map in the GIM product represents the standard deviation of the corresponding TEC map. In theory, the ionospheric delay calculated from the TEC map should have a mean of zero and a standard deviation of σ. GIM_grid The Gaussian distribution of the ionospheric delay standard deviation σ GIM_grid It is determined by the ionospheric delay accuracy value derived from the RMS map.

[0041] The standard deviation of the ionospheric delay estimated by non-difference and non-combination using the known techniques is expressed as σ est_gird , the standard deviation of the ionospheric delay σ est_gird It can be obtained through the covariance matrix. Therefore, the ionospheric residual standard deviation can be linearly convolved as:

[0042]

[0043] In formula (6), σ resgrid represents the standard deviation of the ionospheric residual, σ est_gird represents the standard deviation of the ionospheric delay estimated by undifference and non-combination, σ GIM_gridRepresents the standard deviation of ionospheric delay calculated from the TEC map.

[0044] For the ionospheric grid point residual, the detection statistic res is assumed to be fault-free under the condition H0 grid It has zero mean and standard deviation The unbiased Gaussian distribution of Under the condition H1, the detection statistic res is assumed to be faulty. grid Subject to mean μ and standard deviation σ resgrid has a biased Gaussian distribution, that is, The mean μ represents the deviation caused by abnormal conditions of real-time GIM products.

[0045] The present invention will carry out false alarm error control and missed detection error control based on this assumption.

[0046] Step S5) Calculate the detection threshold according to the false alarm rate required for abnormal integrity monitoring of real-time GIM products, traverse each ionospheric grid point and compare its detection statistic and the detection threshold to perform fault detection. If the detection statistic is less than or equal to the detection threshold, execute step S6; otherwise, if the detection statistic is greater than the detection threshold, mark the corresponding ionospheric grid point as abnormal.

[0047] Under the assumption of no fault condition H0, the probability that the detection statistic is greater than the detection threshold must be less than or equal to the required false alarm rate P fa , the false alarm rate can be expressed as:

[0048]

[0049] In formula (7), P fa is the false alarm rate required for ionospheric grid point monitoring, which is a preset value; x represents the random variable in the cumulative probability density function, The meaning of is shown in formula (6), T grid is the detection threshold at the ionospheric grid point, T grid It can be expressed as:

[0050]

[0051] In formula (8), K fa is the required false alarm rate P fa The corresponding false alarm factor, the value of which can be obtained through the false alarm rate, is the quantile of the Gaussian inverse cumulative distribution function of the corresponding false alarm rate, represents the standard deviation of the ionospheric residual, σ est_gird represents the standard deviation of the ionospheric delay estimated by undifference and non-combination, σ GIM_grid Represents the standard deviation of ionospheric delay calculated from the TEC map.

[0052] By comparing the test statistic res grid and the detection threshold T grid Carry out detection of ionospheric grid point residual abnormalities. When the detection statistic is less than or equal to the detection threshold, it means that the ionospheric grid point residuals in the current epoch have passed the fault detection, and the current ionospheric grid point is normal, and execute step S6; if the detection statistic is greater than the detection threshold, it means that the ionospheric grid point residuals in the current epoch have not passed the fault detection, and the corresponding ionospheric grid point needs to be marked as abnormal.

[0053] Step S6) If the fault detection is passed in step S5, that is, the current ionospheric grid point is normal, then the minimum detectable deviation of the ranging domain of the ionospheric grid point that has passed the fault detection is calculated, and the minimum detectable deviation at each ionospheric grid point is traversed and compared with the preset alarm limit. If the minimum detectable deviation is less than the preset alarm limit, the corresponding ionospheric grid point is normal; otherwise, if the minimum detectable deviation is greater than or equal to the preset alarm limit, the corresponding ionospheric grid point fails the detection, the ionospheric grid point is marked as abnormal, and an alarm is issued to the user.

[0054] Only performing fault detection on whether the ionospheric grid point residuals are abnormal cannot effectively guarantee the integrity of the GIM product. In addition, under fault condition H1, the ionospheric grid point residual distribution model based on the zero-mean Gaussian assumption may fail, and missed detection errors are prone to occur. Combined with the missed detection rate requirement, the minimum detectable deviation for the integrity of the ionospheric grid point detection that passes the fault detection in step S5 is calculated. The missed detection rate can be expressed as:

[0055]

[0056] In formula (9), P MD is the missed detection rate, which is the preset value, T grid is the detection threshold at the ionospheric grid point, and μ is the unknown bias caused by the ionospheric grid point residual anomaly in the detection statistics.

[0057] Calculate the minimum detectable deviation (MDE) used to detect the integrity of the ionospheric grid points grid The method is:

[0058]

[0059] In formula (10), MDE grid represents the minimum detectable deviation in the ranging domain of the ionospheric grid point that has passed the fault detection, T grid is the detection threshold at the ionospheric grid point, K MD Expressed with the required missed detection rate P MD The corresponding Gaussian distribution quantile, represents the standard deviation of the ionospheric residual, σest_gird represents the standard deviation of the ionospheric delay estimated by undifference and non-combination, σ GIM_grid Represents the standard deviation of ionospheric delay calculated from the TEC map.

[0060] By using formulas (6) and (8), we can calculate and T grid , based on the preset missed detection rate P MD You can get K MD (This operation can be realized using existing technology), so that the minimum detectable deviation MDE can be obtained grid .

[0061] Minimum detectable deviation (MDE) grid Compared with the preset alarm limit, if the minimum detectable deviation MDE grid If the constraint of the alarm limit is exceeded, that is, the minimum detectable deviation is greater than or equal to the preset alarm limit, the corresponding ionospheric grid point fails the detection, the current ionospheric grid point is marked as abnormal and an alarm is issued to the user; conversely, if the minimum detectable deviation is less than the preset alarm limit, it indicates that the status of the corresponding ionospheric grid point is normal and the ionospheric grid point has passed the integrity monitoring of the GIM product.

[0062] The above content is a further detailed description of the present invention in combination with specific implementation methods. It cannot be determined that the specific implementation of the present invention is limited to these descriptions. For ordinary technicians in the technical field to which the present invention belongs, some simple deductions or substitutions can be made without departing from the concept of the present invention, which should be regarded as belonging to the scope of protection determined by the claims submitted for the present invention.

Claims

1. A real-time GIM product integrity monitoring method based on non-differential and non-combination, comprising: Step S1) collecting raw observation data from global integrity monitoring stations and receiving real-time satellite clock and orbit correction products, and estimating the first vertical ionospheric delay at the ionospheric puncture point between the monitoring station and the satellite using a non-differenced non-combined model; Step S2) receiving the real-time GIM product and obtaining the second vertical ionospheric delay at the ionospheric grid point; Step S3) selecting at least three first vertical ionospheric delays obtained in step S1 near the ionospheric grid point, and calculating the third vertical ionospheric delay at the ionospheric grid point by inverse distance weighted interpolation method; Step S4) subtracting the second vertical ionospheric delay from the third vertical ionospheric delay, thereby constructing a detection statistic for real-time GIM product integrity monitoring; Step S5) Calculate the detection threshold according to the false alarm rate required for abnormal integrity monitoring of real-time GIM products, traverse each ionospheric grid point and compare its detection statistic and the detection threshold to perform fault detection, if the detection statistic is less than or equal to the detection threshold, execute step S6; otherwise, if the detection statistic is greater than the detection threshold, mark the corresponding ionospheric grid point as abnormal; Step S6) calculates the minimum detectable deviation of the ranging domain of the ionospheric grid point that has passed the fault detection, traverses the minimum detectable deviation at each ionospheric grid point and compares it with the preset alarm limit. If the minimum detectable deviation is less than the preset alarm limit, the corresponding ionospheric grid point is normal; otherwise, if the minimum detectable deviation is greater than or equal to the preset alarm limit, the corresponding ionospheric grid point fails the detection, the ionospheric grid point is marked as abnormal, and an alarm is issued to the user.

2. The method for real-time GIM product integrity monitoring based on non-differential and non-combined according to claim 1, wherein in step S1, the method for estimating the first vertical ionospheric delay at the ionospheric penetration point between the monitoring station and the satellite using the non-differential and non-combined model is: in, I est is the first vertical ionospheric delay, Re is the radius of the earth, α is the satellite altitude angle, h is the distance between the ionospheric shell model and the earth's surface, represents the ionospheric parameters that absorb the hardware delay, further, The superscript s represents the satellite, and the subscript r represents the receiver. Represents the ionospheric parameters absorbing the hardware delay; Represents the pure ionospheric parameters; is the ionosphere-free combination coefficient, γ j (j=1,2) represent the ionospheric mapping factor at frequency j, which can be expressed as DCB r Indicates the differential code deviation at the receiver, DCB r =b r,1 -b r,2 , b r,j (j=1,2) represents the hardware delay of pseudorange at the receiver end; DCB s Indicates the satellite differential code deviation, represents the hardware delay of the pseudorange at the satellite end; j represents the wavelength corresponding to frequency j; Indicates the blur parameter that absorbs the hardware delay; represents the pure ambiguity parameter corresponding to frequency j; B r,j Indicates the hardware delay of the carrier phase at the receiver; Indicates the hardware delay of the satellite carrier phase; t r,IF represents the receiver clock error of the ionosphere-free combination, Represents the satellite clock error of the ionosphere-free combination.

3. The method for real-time GIM product integrity monitoring based on non-differential and non-combined according to claim 1, wherein in step S2, the method for obtaining the second vertical ionospheric delay at the ionospheric grid point is: Among them I GIM_grid represents the second vertical ionospheric delay at the ionospheric grid point, TEC grid represents the total electron content at the ionospheric grid point, and f represents the frequency.

4. The method for real-time GIM product integrity monitoring based on non-differential and non-combination according to claim 1, wherein in step S3, the method for obtaining the third vertical ionospheric delay at the ionospheric grid point is: Among them I est_grid is the third vertical ionospheric delay at the ionospheric grid point, I est_i is the first vertical ionospheric delay calculated by the i-th monitoring station, d i is the distance from the ith monitoring station to the station to be solved, m is the number of monitoring stations selected around the station to be solved, and m≥3.

5. The method for real-time GIM product integrity monitoring based on non-differential and non-combined according to claim 1, wherein in step S3, if there is a lack of ionospheric puncture points between the monitoring station and the satellite in the vicinity of the ionospheric grid point, the status of the ionospheric grid point is marked as not monitored.

6. The method for real-time GIM product integrity monitoring based on non-differential and non-combination according to claim 1, wherein in step S4, the method for constructing the detection statistic for real-time GIM product integrity monitoring is: res grid =I GIM_grid -I est_grid where res grid represents the detection statistic used for real-time GIM product integrity monitoring, I GIM_grid represents the second vertical ionospheric delay at the ionospheric grid point, I est_grid Represents the third vertical ionospheric delay at the ionospheric grid point.

7. The method for real-time GIM product integrity monitoring based on non-differential and non-combination according to claim 1, wherein in step S5, the method for calculating the detection threshold is: Where T grid is the detection threshold at the ionospheric grid point, K fa is the false alarm factor corresponding to the required false alarm rate, σ est_gird represents the standard deviation of the ionospheric delay estimated by undifference and non-combination, σ GIM_grid Represents the standard deviation of ionospheric delay calculated from the TEC map.

8. The method for real-time GIM product integrity monitoring based on non-differential non-combination according to claim 1, wherein in step S6, the method for calculating the minimum detectable deviation of the ranging domain of the ionospheric grid point that has passed the fault detection is: Among them, MDE grid represents the minimum detectable deviation in the ranging domain of the ionospheric grid point that has passed the fault detection, T grid is the detection threshold at the ionospheric grid point, K MD represents the Gaussian distribution quantile corresponding to the required missed detection rate, σ est_gird represents the standard deviation of the ionospheric delay estimated by undifference and non-combination, σ GIM_grid Represents the standard deviation of ionospheric delay calculated from the TEC map.

Images