RTK positioning-oriented slow-varying fault detection method based on hybrid control chart

Abstract

The invention relates to an RTK positioning-oriented slowly-varying fault detection method based on a hybrid control chart. The method comprises the following steps of: positioning by using a pseudo-range and a carrier phase observation value received at a GNSS satellite through a user receiver, and obtaining a carrier phase weighted innovation observation value; filtering test statistics generated according to the carrier phase weighted innovation observation value through the hybrid control chart; and comparing the test statistics with a threshold value to judge whether the slow change abnormity exists or not. Weighted information is used as a monitoring object, the windowing expansion multivariate mixing uniform control chart and the parallel multivariate accumulation control chart are connected in series to remarkably improve the response efficiency to the slow-varying fault, the method has the technical advantage of effectively extracting the slow-varying fault in a complex environment, the slow-varying abnormity of an RTK user can be accurately and rapidly monitored, and the method is suitable for being popularized and applied. Therefore, the requirement of future ITS users on satellite navigation service integrity is improved.

Description

Technical Field

[0001] The present invention relates to a technology in the field of satellite navigation, specifically a slow-varying fault detection method based on a hybrid control chart for real-time kinematic (RTK) positioning for intelligent transportation system (ITS) services. Background Art

[0002] The receiver autonomous integrity monitoring (RAIM) technology based on the least squares method is used for fault detection and exclusion (FDE). When the global navigation satellite system (GNSS) is interfered, it may transmit incorrect navigation information. However, when applying RTK services for high-precision positioning, the traditional RTK-RAIM technology mainly targets step faults at a single point, only focuses on the amplitude of the faults, and ignores the trend information of the faults, resulting in insufficient sensitivity to slow-varying faults, making it difficult to detect faults before their amplitudes exceed the threshold, and limiting the performance of FDE in real complex environments. Summary of the Invention

[0003] Aiming at the deficiencies of the existing technology in terms of low positioning accuracy and inability to achieve high-precision applications based on carrier phase, the present invention proposes a slow-varying fault detection method based on a hybrid control chart for RTK positioning. Taking the normalized innovation observation value as the monitoring object, a hybrid control chart is constructed by connecting the windowed extended multivariate mixed uniform control chart (MHWMA) and the parallel multi-variable cumulative sum control chart (MCUSUM) in series, so as to significantly improve the response efficiency to slow-varying faults, have the technical advantage of effectively extracting slow-varying faults in complex environments, and can achieve accurate and rapid monitoring of RTK users suffering from slow-varying anomalies, thereby meeting the requirements of future ITS users for the integrity of satellite navigation services.

[0004] The present invention is realized through the following technical solutions:

[0005] The present invention relates to a slow-varying fault detection method based on a hybrid control chart for RTK positioning, including:

[0006] Step 1, the user receiver uses the observation values received from GNSS satellites for positioning and obtains the normalized innovation observation value, specifically including:

[0007] Step 1-1) Calculate the double-difference observation value: For a reference station with a known precise position and a rover station with an unknown position, synchronously observe two satellites, including calculating the double difference of the GNSS observation values broadcast by dual-frequency GNSS. After setting the satellite with the largest fixed elevation angle as the reference satellite, traverse all satellites of this frequency point under the same constellation and calculate the double-difference value of each satellite.

[0008] The so-called double difference refers to the difference between the single differences of the observation values obtained by the reference station and the rover station when synchronously observing two satellites.

[0009] The single difference mentioned above refers to the difference between the observation values obtained by the reference station and the rover station observing the same satellite synchronously.

[0010] Step 1-2) After initializing the navigation state and covariance, use the extended Kalman filter (EKF) for RTK relative positioning: Combine the linear motion model and the non-linear measurement model, and through continuous iterative state prediction and measurement update steps, estimate the relative position, velocity, and double-difference integer ambiguity of the rover receiver relative to the reference station in real time, so as to realize the estimation and update of the state, where: The linear motion model describes the motion mode of the receiver, and the non-linear measurement model establishes the correlation between the carrier phase observation data and the position of the receiver.

[0011] Step 1-3) Obtain the normalized innovation observation value, specifically: For dual-frequency observables, obtain the carrier phase normalized innovation vector where: z k is the innovation vector output during the Kalman filter iteration at time k; the inverse of the covariance matrix of the z vector is the Choleskey decomposition of the W matrix at the current time, H is the linearized design matrix at the current time, Φ k is the state prediction transfer matrix at the current time, is the state estimate value after covariance correction at time k-1, Q is the state noise matrix, and R is the observation noise matrix. When using GNSS observations of one frequency point, q is a (M-1)×1 dimensional vector, and M is the number of visible satellites.

[0012] Step 2, after initializing the parameters, filter and smooth the carrier phase normalized innovation observation value q through the window-expanded MHWMA control chart, and input the smoothed observation value into the parallel MCUSUM filter to extract the growth trend of the slow-varying fault, and use the output of the hybrid control chart as the test statistic, specifically including:

[0013] Step 2-1) Initialize the parameters, including the sliding window length m and the prior fault amplitude (a, b), where: a is the upper bound of the prior fault amplitude, and b is the lower bound;

[0014] Step 2-2) Use the window-expanded MHWMA control chart to filter and smooth q k For the normalized innovation observation value at time k, calculate the filtered and smoothed result where: ω1∈(0,1] is the weight coefficient 1, ω2∈(0,ω1] is the weight coefficient 2, is the average value of q within the window before time k. When k≤m, that is, when the control chart is in the initial operation stage, When k>m, that is, when the control chart is in the stable initial stage,

[0015] The essence of the windowed extended MHWMA control chart is to assign a higher weight to the current observed quantity q k After that, through windowing processing, only the influence of the previous m moments within a finite epoch is considered, and the weights of earlier moments are all assigned a value of 0.

[0016] When , the windowed extended MHWMA control chart degenerates into a multivariate Shewhart control chart. According to the formula, the mean expectation of WEH k and the covariance expectation Σ can be obtained. WEH,k .

[0017] Step 2 - 3) Input the smoothed observed value WEH k into the parallel MCUSUM filter to extract the growth trend of the slow - varying fault, and calculate the test statistic of the hybrid control chart at time k, specifically: where: μ0 is the expectation of q, assumed to be a vector of all 0s; A 1,k is the sum of WEH from time k - n 1,k +1 to time k, A 2,k is the sum of WEH from time k - n 2,k +1 to time k, n 1,k and n 2,k are counting ordinals, and their values are related to B 1,k-1 and B 2,k-1 respectively. When B 1,k-1 ≤0, n 1,k =1; when B 1,k-1 >0, n 1,k =1 + n 1,k-1 . When B 2,k-1 ≤0, n 2,k =1; when B 2,k-1 >0, n 2,k =1 + n 2,k-1 . (θ1, θ2) is the reference value.

[0018] The reference value (θ1, θ2) is related to the amplitude of the fault to be detected. When it is desired to test for a fault δ with an amplitude between a and b, k1=(3a + b) / 8 and k2=(a + 3b) / 8.

[0019] Step 3, use the test statistic to compare with the threshold to determine whether there is a slow - varying anomaly, specifically including:

[0020] Step 3 - 1) Initialize the false - alarm probability P fa ;

[0021] Step 3-2) Calculate the thresholds h1 and h2 offline using the Monte Carlo method, and the thresholds simultaneously satisfy the constraints ①k1h1=k2h2; ②k1+h1>k2+h2;

[0022] The Monte Carlo method is to obtain the statistical characteristics of q based on a large amount of historical data, including the mean μ0 and the covariance Σ0. Offline generation Random generation 20 / P fa A multidimensional Gaussian sample with mean μ0 and covariance Σ0 is obtained by combining these Gaussian samples through a joint control chart and obtaining 20 / P fa B 1,k With B 2,k Set the initial value of h1 to 0.01, and h2 is always equal to h1k1 / k2, and gradually increase h1 and h2. If B is satisfied 1,k >h1 or B 2,k >h1's B 1,k With B 2,k When the number of samples reaches 100 and k1+h1>k2+h2, h1 and h2 are output as the thresholds. Otherwise, h1 and h2 are gradually increased until the above conditions are met. Repeat the experiment 100 times and use h1 and h2 of all experiments as the final threshold.

[0023] Step 3-3) Use the binary hypothesis test theory to determine whether the user-side RTK system has a slope failure. ① Null hypothesis H0: No failure, B 1,k ≤h1 and B 2,k ≤h2;②Alternative hypothesis H1: There is a fault, B 1,k >h1 or B 2,k ≤h2. When the alternative hypothesis is true, the test statistic is greater than the threshold, and it is judged that a fault exists.

[0024] The present invention relates to an RTK-oriented sequential fault monitoring system for realizing the above-mentioned method, comprising: a GNSS observation receiving and data preprocessing unit, a filtering unit based on a hybrid control diagram, and a fault detection unit based on a binary hypothesis theory, wherein: the GNSS observation receiving and data preprocessing unit performs signal acquisition and Kalman filtering processing according to GNSS observation information to obtain normalized new information observations, the filtering unit based on the hybrid control diagram performs filtering according to the normalized new information observation information using a windowed extended MHWMA control diagram in parallel with an MCUSUM filter, the filtering result obtained by the hybrid control diagram is used as a test statistic, and the fault detection unit based on the binary hypothesis theory calculates a threshold value according to the above-mentioned test statistic and by a Monte Carlo method, compares the two, and obtains a judgment result of whether the current RTK system has a fault. Technical Effects

[0025] In view of the lack of means for detecting RTK client ramp faults in existing research / inventions, the present invention proposes a new hybrid control chart for RTK high-precision positioning services and applies it to RTK fault detection. The hybrid control chart first includes a new windowed extended MHWMA. Compared with the existing ordinary MHWMA, the windowed extended MHWMA uses an extended form while only smoothing within the window, and improves the degree of freedom by introducing new parameters; the hybrid control chart also includes a parallel MCUSUM control chart. Compared with the known ordinary MCUSUM control chart that requires a priori fault amplitude, it can detect faults when only the a priori fault range is known, with higher flexibility. Therefore, the fault detection method proposed by the present invention has a smaller detection delay when a ramp fault occurs. BRIEF DESCRIPTION OF THE DRAWINGS

[0026] Figure 1 is a flow chart of the present invention;

[0027] Figure 2 is the normalized innovation diagram calculated in the embodiment without faults;

[0028] FIG. 3 is the detection response time of the embodiment when a fault is injected;

[0029] FIG. 3(a) shows the test statistic and threshold calculated by the present method, and FIG. 3(b) shows the test statistic and threshold obtained by the prior art. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0030] As Figure 1 shown, this embodiment relates to a slow-varying fault detection method based on a hybrid control chart for RTK positioning, including:

[0031] Step 1: The receiver receives dual-frequency pseudorange and carrier phase observations from GNSS satellites, calculates the double difference and uses EKF for positioning, and outputs the normalized innovation observable. Specifically, it includes:

[0032] Step 1-1) After the receiver receives the pseudorange and carrier phase observations from GNSS satellites, it calculates the double difference. Specifically, assuming there are M visible satellites in the line of sight, for single-frequency observables, there are (M - 1) groups of double-difference observables. Taking the satellite with the largest fixed elevation as the reference satellite (numbered 1), it traverses all satellites of this frequency point under the same constellation and calculates the double-difference value of each satellite.

[0033] The double difference is calculated in the following way: Taking the pseudorange as an example, the pseudorange observations of the rover r with unknown position received from satellite No. 1 and satellite No. 2 are respectively where: f is the observation frequency. The pseudorange observations of the reference station r with known position received from satellite No. 1 and satellite No. 2 are respectively In meters, the single differences in pseudorange between the rover and the reference station for satellites 1 and 2 are The double difference in pseudorange is Where: is the double difference geometric distance; is the double difference pseudorange noise. Similarly, the double difference pseudorange observations for other satellites and the reference satellite can be obtained. Similar to the double difference pseudorange observations, the double difference carrier observations are Where: is the double difference carrier noise between the reference satellite 1 and satellite j, λ f is the wavelength, is the double difference integer ambiguity, is the double difference carrier phase noise.

[0034] Step 1-2) After initializing the navigation state and covariance, use the Extended Kalman Filter (EKF) for RTK relative positioning. Specifically: The state to be solved at time k is Where: is the relative position of the reference station and the rover in the x / y / z directions, is the velocity of the rover in the x / y / z directions.

[0035] To estimate X k , first construct the state prediction equation Where: is the predicted state, and Φ k,k-1 is the dynamics transition matrix used to determine the motion model; at the same time, construct the covariance prediction matrix Where: is the predicted covariance, P k-1 is the covariance matrix at time k-1, which is also the initial covariance matrix at time k, and Q is the noise matrix; construct the state update equation as follows: Where: z k is the innovation, y k is the true double difference value calculated using the observations at time k, is the double difference value predicted using the state prediction equation, where: h is the design matrix; at the same time, construct the inverse matrix of the covariance matrix of z Where: H k is the linearized design matrix, and R is the measurement noise matrix; calculate the Kalman gain Finally, complete the update of the state and covariance and Where: is the estimated value of the state at time k, is the estimated value of the covariance at time k. Therefore, the relative position and velocity of the rover can be obtained.

[0036] Step 1-3) Output the normalized innovation observable, specifically: for the single-frequency observable, where: q is the normalized innovation observation value of the carrier phase, and is an (M - 1)×1 dimensional vector. is the Choleskey decomposition of the W k matrix.

[0037] Step 2: After initializing the parameters, filter and smooth the normalized innovation observation value q of the carrier phase through the window-expanded MHWMA control chart, and input the smoothed observation value into the parallel MCUSUM filter to extract the growth trend of the slow-varying fault. Use the output of the hybrid control chart as the test statistic, specifically including:

[0038] Step 2-1) Initialize the parameters, including the sliding window length m and the prior fault amplitudes (a, b).

[0039] Step 2-2) Use the window-expanded MHWMA control chart to filter and smooth q: for the normalized innovation observable at time k, calculate the filtered and smoothed result where: ω1 ∈ (0, 1] is the weight coefficient 1, ω2 ∈ (0, ω1] is the weight coefficient 2, is the average value of q within the window before time k, which is related to the current time and the window length. When k ≤ m, that is, when the control chart is in the initial running stage, When k > m, that is, when the control chart is in the stable initial stage,

[0040] The essence of the window-expanded MHWMA control chart is to assign a higher weight to the observable q at the current time, and through window processing, only consider the influence of the m previous moments within a finite epoch, and assign all weights of earlier moments to 0. k When

[0041] When , the window-expanded MHWMA control chart degenerates into a multivariate Shewhart control chart. At the same time The value of cannot be too small, otherwise false alarms are likely to occur. Therefore, in this embodiment

[0042] After obtaining WEH k , it is necessary to calculate its mean and covariance under the fault-free condition. WEH k The expected mean is easily proven to be the expected mean μ0 of q k , which is 0 under normal conditions; WEH k The covariance matrix can be obtained through mathematical derivation. If the original covariance of q is Σ0, then there is

[0043] Step 2-3) Input the smoothed observed value WEH k into the parallel MCU SUM filters to extract the growth trend of the slow-varying fault, and calculate the test statistic of the hybrid control chart at time k, specifically: where: μ0 is the mean expectation of q, which is a vector of all 0s under nominal conditions; A 1,k is the sum of WEH from time k-n 1,k +1 to time k, and A 2,k is the sum of WEH from time k-n 2,k +1 to time k, n 1,k and n 2,k are counting ordinals, and their values are related to B 1,k-1 and B 2,k-1 respectively. When B 1,k-1 ≤0, n 1,k =1; when B 1,k-1 >0, n 1,k =1 + n 1,k-1 . When B 2,k-1 ≤0, n 2,k =1; when B 2,k-1 >0, n 2,k =1 + b 2,k-1 . (θ1, θ2) are reference values.

[0044] The reference values (θ1, θ2) are related to the amplitude of the fault to be detected. When it is desired to test for a fault δ with an amplitude between a and b, k1 = (3a + b) / 8 and k2 = (a + 3b) / 8.

[0045] Step 3: Use the test statistic to compare with the threshold to determine whether there is a slow-varying anomaly, specifically including:

[0046] Step 3-1) Initialize the false alarm probability P fa ;

[0047] Step 3-2) Use the Monte Carlo method to calculate the thresholds h1 and h2 offline, and the thresholds simultaneously satisfy the constraints ① k1h1 = k2h2; ② k1 + h1 > k2 + h2, specifically including:

[0048] Step 3-2-1) Initialize the parameters, let h1 = 0.01, h2 = h1(a + 3b) / (3a + b), ordinal l is 1; obtain the statistical characteristics of q based on a large amount of historical data, that is, its mean μ0 and covariance Σ0;

[0049] Step 3-2-2) Perform the following operations: ① Generate multi-dimensional Gaussian random numbers with a mean of μ0 and a covariance of Σ0, and the number of random numbers is 20 / P fa ; ② Expand the random numbers through windowed MHWMA and multi-source parallel CUSUM according to Step 2 to obtain 20 / P fa h1's and h2's, and sort h1 and h2 by size; ③ If both p1 + h1 > p2 + h2 and sum(T1 < Thre l,1 ) + sum(T2 < Thre l,2 ) - sum(T1 < Thre l,1 ∩T2 < Thre l,2 ) < 20, where: sum(*) represents the sum of the number of events satisfying (), then jump out of the loop and store the h l,1 and h l,2 at this time. Otherwise, increase h l,1 by 0.01 each time, assign h l,2 to h l,1 (a + 3)(3a + b), and return to Step ③;

[0050] Step 3-2-3) Execute Step 3-2-2) 200 times, each time increasing the value of l by 1, and store all the h l,1 and h l,1 ;

[0051] Step 3-2-4) Assign h1 to and assign h2 to where: and are the means of all h l,1 and h l,1 , and output h1 and h2 as the thresholds for Monte Carlo calculation;

[0052] Step 3-3) Use the binary hypothesis testing theory to determine whether the user-side RTK system has a ramp fault. ① The null hypothesis H0: no fault, B 1,k ≤h1 and B 2,k ≤h2; ② The alternative hypothesis H1: there is a fault, B 1,k >h1 or B 2,k ≤h2. When the null hypothesis holds, if the test statistic is less than the threshold, it is judged that no fault occurs, increment k by one epoch, and feedback to Step 1. When the alternative hypothesis holds, if the test statistic is greater than the threshold, it is judged that a fault exists, and an alarm message is broadcast.

[0053] Through specific experiments, an example of RTK positioning with the dual-frequency observation signals (L1 and L2) of the GPS navigation system is analyzed. The observation data is sourced from the Hong Kong Satellite Positioning Service Reference Station, with a data sampling rate of 1HZ and a satellite cutoff elevation angle of 15 degrees. In this embodiment, the HKPC station is used as the reference station with a known position, and the KYC1 station is used as the roving station with an unknown position. The baseline distance between the two is less than 5 km, belonging to a short baseline. There are a total of 8 visible GPS satellites, and G19 is used as the reference satellite. Figure 2 It is calculated based on the observations of the observation file obtained on January 2, 2022. It is the normalized innovation diagram of different satellites and different frequency points without faults between 17:35 and 17:40 on that day.

[0054] Since the fault occurrence frequency is relatively low, in this embodiment, a slowly varying ramp fault is artificially injected into the L1 carrier phase observation value of the GPS G02 satellite in the fault-free observation file to simulate a ramp fault with a change rate of 0.003 cycles / s in the observed quantity, lasting for 300 s. To illustrate the superiority of the method based on the hybrid control chart in terms of detection time, the detection method proposed in this embodiment is compared with the existing method, and the inspection time is shown in Figure 3. For the control chart parameters, m = 40, and it is assumed that δ ∈ (5, 7), so (k1, k2) = (2.75, 3.25).

[0055] As shown in Figure 3(a), it is the test statistic and threshold obtained by the proposed hybrid control chart. As shown in Figure 3(b), it is the test statistic and threshold obtained by the traditional fault detection method. Obviously, the method based on the hybrid control chart disclosed in this embodiment detected the fault at 230 seconds, while the test statistic of the traditional method was always less than the threshold, resulting in a missed detection event. This proves the superiority of the proposed method.

[0056] Compared with the prior art, the windowed extended MHWMA of the present invention uses an extended form while only using the smoothing inside the window, and improves the degree of freedom by introducing new parameters; the hybrid control chart also includes a parallel MCUSUM control chart. Compared with the known ordinary MCUSUM control chart that depends on the prior fault amplitude, it can perform detection only knowing the prior fault range, with higher flexibility. Compared with the existing classical fault detection technologies, the present invention first performs anomaly detection on the slow and slowly varying faults that affect RTK positioning and service security without the need for complex numerical calculations. The system has a significant improvement in response speed for slow-varying anomaly detection, and has important theoretical and application values for the high-precision positioning service and safety guarantee of ITS users.

[0057] Those skilled in the art can make partial adjustments to the above specific embodiments in different ways without departing from the principles and purposes of the present invention. The protection scope of the present invention is subject to the claims and is not limited by the above specific embodiments. All implementation solutions within its scope are subject to the present invention.

Claims

1. A slow-varying fault detection method based on a hybrid control chart for RTK positioning, characterized in that Including: Step 1: The user receiver uses the observations received from GNSS satellites for positioning and obtains normalized innovation observations, specifically including: Step 1-1) Calculate double-difference observations: For a reference station with known position and a rover station with unknown position, synchronously observe two satellites, including calculating the double-difference of GNSS observables broadcast by dual-frequency GNSS; after setting the satellite with the largest fixed elevation angle as the reference satellite, traverse all satellites of this frequency point under the same constellation, and calculate the double-difference value of each satellite; Step 1-2) After initializing the navigation state and covariance, use the extended Kalman filter for RTK relative positioning: Combine the linear motion model and the nonlinear measurement model, and through continuous iterative state prediction and measurement update steps, estimate the relative position, velocity, and double-difference integer ambiguity of the rover receiver relative to the reference station in real time, realizing the estimation and update of the state, where: The linear motion model describes the motion mode of the receiver, and the nonlinear measurement model establishes the correlation between the carrier phase observation data and the position of the receiver; Step 1-3) Obtain the normalized innovation observation value, specifically: for the dual-frequency observation quantity, obtain the carrier phase normalized innovation vector where: z k is the innovation vector output during the Kalman filter iteration at time k; the inverse of the covariance matrix of the z vector is the Choleskey decomposition of the W matrix at the current time, H is the linearized design matrix at the current time, Φ k is the state prediction transfer matrix at the current time, is the state estimate value after covariance correction at time k-1, Q is the state noise matrix, R is the observation noise matrix; when using GNSS observations of one frequency point, q is an (M-1)×1 dimensional vector, and M is the number of visible satellites; Step 2: After initializing the parameters, filter and smooth the carrier phase normalized innovation observation value q through a windowed extended MHWMA control chart, and input the smoothed observation value into a parallel MCUSUM filter to extract the growth trend of slow-varying faults, using the output of the hybrid control chart as the test statistic, specifically including: Step 2-1) Initialize the parameters, including the sliding window length m and the prior fault amplitude (a, b), where: a is the upper bound of the prior fault amplitude, and b is the lower bound; Step 2-2) Use the windowed extended MHWMA control chart to filter and smooth q k For the normalized innovation observation value at time k, calculate the result after filtering and smoothing where: ω1 ∈ (0, 1] is the weight coefficient 1, ω2 ∈ (0, ω1] is the weight coefficient 2, is the average value of q within the window before time k. When k ≤ m, that is, when the control chart is in the initial operation stage, When k > m, that is, when the control chart is in the initial stable stage, Step 2-3) Input the smoothed observed value WEH k into the parallel MCU SUM filters to extract the growth trend of the slow-varying fault, and calculate the test statistic of the hybrid control chart at time k, specifically as follows: where: μ0 is the expectation of q, set as a vector of all zeros; A 1，k is the sum of WEH from time k - n 1，k +1 to time k, A 2，k is the sum of WEH from time k - n 2，k +1 to time k, n 1，k and n 2，k are counting ordinals, and their values are related to B 1，k-1 and B 2，k-1 respectively; when B 1，k-1 ≤0, n 1，k =1; when B 1，k-1 >0, n 1，k =1 + n 1，k-1 ; when B 2，k-1 ≤0, n 2，k =1; when B 2，k-1 >0, n 2，k =1 + n 2，k-1 ; (θ1, θ2) are reference values; Step 3: Use the test statistic to compare with the threshold to determine whether there is a slow-varying anomaly, specifically including: Step 3-1) Initialize the false alarm probability P fa ; Step 3-2) Use the Monte Carlo method to calculate the thresholds h1 and h2 offline, and the thresholds simultaneously satisfy the constraints ① k1h1 = k2h2; ② k1 + h1 > k2 + h2; Step 3-3) Use the binary hypothesis testing theory to determine whether a ramp fault has occurred in the user-side RTK system. ① Null hypothesis H0: No fault, B 1，k ≤h1 and B 2，k ≤h2; ② Alternative hypothesis H1: Fault exists, B 1，k >h1 or B 2，k ≤h2; When the alternative hypothesis holds, the test statistic is greater than the threshold, and it is determined that a fault exists.

2. The slow-varying fault detection method based on a hybrid control chart for RTK positioning according to claim 1, wherein The so-called double-difference refers to: the difference between the single-differences of the observations obtained by the reference station and the rover station synchronously observing two satellites; The so-called single-difference refers to: the difference between the observations obtained by the reference station and the rover station synchronously observing the same satellite.

3. The slow-varying fault detection method based on a hybrid control chart for RTK positioning according to claim 1, characterized in that, The essence of the windowed extended MHWMA control chart is to assign a higher weight to the observed quantity q at the current moment. k After windowing processing, only the influence of the m moments before the current moment is considered, and the weights of earlier moments are all assigned a value of 0.

4. The slow-varying fault detection method based on a hybrid control chart for RTK positioning according to claim 1, characterized in that, When the windowed extended MHWMA control chart degenerates into a multivariate Shewhart control chart, and the mean expectation k of WEH and the covariance expectation ∑ WEH，k are calculated.

5. The slow-varying fault detection method based on a hybrid control chart for RTK positioning according to claim 1, wherein The reference values (θ1, θ2) are related to the amplitude of the fault to be detected. When it is desired to test a fault δ with an amplitude between a and b, k1 = (3a + b) / 8 and k2 = (a + 3b) / 8.

6. The slow-varying fault detection method based on a hybrid control chart for RTK positioning according to claim 1, characterized in that, The so-called Monte Carlo method refers to: i) Obtain the mean μ0 and covariance ∑0 of q from historical data; ii) Generate 20 / P random numbers offline fa multidimensional Gaussian samples that follow a mean of μ0 and a covariance of ∑0. By using the combined control chart, obtain 20 / P fa samples of B 1，k and samples of B 2，k ; iii) Set the initial value of h1 to 0.01, and h2 is always equal to h1k1 / k2, and gradually increase h1 and h2; iv) When B is satisfied 1，k > h1 or B 2，k > h1 of B 1，k and B 2，k The number of samples reaches 100, and at this time k1 + h1 > k2 + h2, then output h1 and h2 as thresholds, otherwise continue to gradually increase h1 and h2 until the above conditions are met; v) Repeat steps i) to iv) 100 times, and use the mean values of h1 and h2 obtained as the final thresholds.

7. A sequential fault monitoring system for RTK that implements the above method of the claims, comprising: GNSS observation receiving and data preprocessing unit, a filtering unit based on a hybrid control chart, and a fault detection unit based on the binary hypothesis theory, where: The GNSS observation receiving and data preprocessing unit performs signal acquisition and Kalman filtering processing according to the GNSS observation information, and obtains the normalized innovation observation. The filtering unit based on the hybrid control chart filters according to the normalized innovation observation information by using a windowed extended MHWMA control chart in parallel with an MCUSUM filter. The filtering result obtained by the hybrid control chart is used as the test statistic. The fault detection unit based on the binary hypothesis theory calculates the threshold through the Monte Carlo method according to the obtained test statistic above, and compares the two to obtain the judgment result of whether the current RTK system has a fault.

Images