GNSS positioning slowly varying fault detection method based on residual-svr regression

Abstract

The GNSS positioning slowly-varying fault detection method based on residual-SVR regression has the steps that: the observation innovation sequence is obtained based on a Kalman filter, and a covariance matrix thereof is calculated; a chi-square statistic is constructed by accumulating multi-step innovations through a sliding window; a training set is constructed based on fault-free data, with the innovation sequence as input and the chi-square statistic as output, and an innovation-statistic mapping function is generated; a normal slope threshold is fitted based on an SVR predicted value, a least square fitting curve of the observation statistic is realized in a real-time sliding window, and it is judged whether to start a slowly-varying fault alarm or not.The present application directly takes the residual sequence as the model input, replaces the traditional two-classification label with the dynamic chi-square statistic, reduces the detection delay, utilizes the SVR detection model based on the grid search and cross-validation collaborative optimization, selects the optimal parameter combination minimizing the mean square error (MSE) through the logarithmic uniform sampling, interval linear sampling and five-fold cross-validation, and reduces the slowly-varying fault detection average absolute error.

Description

Technical Field

[0001] This invention relates to a method for detecting slowly varying faults in GNSS positioning based on residual-SVR regression, belonging to the field of satellite positioning technology. Background Technology

[0002] Currently, the detection of slowly varying faults in Global Navigation Satellite Systems (GNSS) (such as satellite clock drift and multipath cumulative error) has the following limitations: the RAIM algorithm based on geometric distribution is insensitive to continuous small changes; Kalman filtering relies on linear assumptions and precise state equations, making it difficult to model nonlinear slowly varying processes; while deep learning (such as LSTM) can capture temporal features, it requires massive amounts of fault samples and has high computational complexity, making it difficult to meet the real-time requirements of scenarios such as aviation. In contrast, Support Vector Regression (SVR) directly models nonlinear features through kernel function techniques (such as RBF), requiring only a small amount of normal data to build a high-precision prediction model. Its ε-insensitive band design can adaptively filter out GNSS environmental noise, and it achieves dynamic threshold adjustment through sliding window residual statistics (such as mean ± 3σ).

[0003] In fault detection of integrated navigation systems, traditional methods often face the challenge of balancing detection efficiency and accuracy, specifically manifested in the following technical bottlenecks: (1) Redundancy of input features and insufficient adaptability of output statistics: Existing methods mostly use raw data from multiple sensors (such as accelerometer and gyroscope outputs) or high-order fusion features (such as attitude angle error and velocity residual) as input to the fault detection model, resulting in high feature dimensions and large computational load. In addition, most schemes directly use binary classification labels (normal / fault) as model outputs, lacking the ability to quantitatively evaluate the fault intensity, which makes it difficult to meet the incremental decision-making requirements for slowly changing fault detection. For example, although neural network-based classifiers can handle complex features, they have problems such as high inference latency (>50ms) and inability to output fault confidence in real time; (2) Coarse optimization of model parameters and limited generalization performance: Traditional support vector regression (SVR) methods rely heavily on empirical parameters (such as penalty factor C and kernel function bandwidth γ) or local grid search in fault detection, resulting in significant prediction bias of the model in dynamic noise environments. Especially when the sensor noise distribution changes abruptly due to environmental interference, the fixed parameter model is prone to the phenomenon of threshold being too tight (high false alarm rate) or too loose (high false alarm rate); (3) Insufficient real-time detection in dynamic scenarios: The existing residual analysis method needs to calculate multi-dimensional state estimation errors (such as position, velocity, and attitude residuals) simultaneously, and generate statistics through complex hypothesis testing (such as generalized likelihood ratio test), resulting in a single detection taking more than 100ms, which is difficult to meet the real-time requirements of highly dynamic scenarios (such as drone obstacle avoidance). In addition, the traditional chi-square detection is only used as a post-test verification method and is not deeply coupled with the machine learning model, so it cannot achieve end-to-end adaptive threshold adjustment. Summary of the Invention

[0004] This invention provides a GNSS positioning gradually varying fault detection method based on residual-SVR regression. This method has low input feature dimension and low computational load, and can achieve continuous quantization output of fault intensity, reducing detection latency. It also has strong generalization performance for the characteristics of gradually varying faults, meeting the real-time detection requirements of high dynamic scenarios.

[0005] To achieve the above objectives, this invention provides a method for detecting slowly varying GNSS positioning faults based on residual-SVR regression, comprising the following steps:

[0006] S1. Obtain the observation innovation sequence based on the Kalman filter (EKF) and calculate its covariance matrix; accumulate multiple steps of innovation through a sliding window, construct the chi-square statistic, and quantify the fault distribution characteristics.

[0007] S2. Based on fault-free data, a training set is constructed with the innovation sequence as input and the chi-square statistic as output. The radial basis function is used, and the SVR parameters are optimized by grid search combined with five-fold cross-validation to generate the innovation-statistic mapping function.

[0008] S3. Fit the normal slope threshold based on the SVR predicted value, and fit the observed statistical curve using the real-time sliding window least squares. If the slope exceeds the threshold range multiple times in a row, a slow change fault alarm will be triggered.

[0009] Furthermore, the specific process of S1 is as follows:

[0010] S1.1 Define the state vector of the GNSS positioning system. The state vector includes position error components δp corresponding to the three orthogonal directions: East (E), North (N), and U. e ,δp n ,δp u and velocity error component δv e ,δv n ,δv u and receiver clock deviation δt clk and its drift rate δt drift ;

[0011] The dynamic evolution model of clock skew characterizes the cumulative error characteristics of the receiver clock, and its expression is:

[0012] Where k represents the discrete-time index of the current moment. Δ represents the clock drift at the previous time k-1. t For the filter period, w t For clock noise that follows a Gaussian distribution, This represents the noise variance, reflecting the stability of the clock crystal oscillator.

[0013] S1.2, Define the state transition matrix The position and velocity block matrices used to predict the state variables at the next time step are as follows: Δt∈[10ms,100ms];

[0014] S1.3 Calculate the actual pseudorange observation Z based on the Extended Kalman Filter (EKF) framework. k Compared with the predicted value H k X k / k-1 residual r k , Z k This is the actual pseudorange observation vector. ρ M This is the pseudorange measurement value of the Mth visible satellite; This is the state prediction vector. x r y r z r For receiver position, v x v y v z For receiver speed, δt clk δt represents the receiver clock offset. drift H represents the receiver drift rate. k Here is the observation matrix; observation matrix H k Each row corresponds to the geometric relationship of a visible satellite, and its elements are determined by the geometric partial derivatives of the satellite's position with respect to the receiver's position. The expression is as follows: in, This represents the partial derivative of the pseudorange of the i-th satellite with respect to the receiver's three-dimensional position, with 1 and 0 corresponding to the observation coefficients of clock bias and drift rate, respectively.

[0015] S1.4 Update the residual covariance matrix, the residual covariance matrix V k From the predicted covariance matrix P k / k-1 With the observation noise matrix R k The joint calculation uses the following update formula: Among them, P k / k-1 =φP k-1 Φ T +Q is the state prediction covariance matrix, P k-1 This is the state estimation covariance matrix at time k-1, representing the confidence level of the system state at the previous time step. In the calculation of the state prediction covariance: P k / k-1 =ΦP k-1 Φ T +Q, where Φ represents the state transition matrix and Q represents the process noise covariance;

[0016] The process noise covariance matrix is... and These represent the process noise variances for position and velocity, respectively, reflecting the uncertainty of the system model; The process noise variance of the receiver clock skew. The process noise variance is the clock drift rate;

[0017] Observation noise matrix It is a diagonal matrix, and its elements are... The pseudorange measurement noise variance of the i-th satellite is calculated in real time by the satellite signal quality assessment module.

[0018] S1.5. Constructing a sliding window chi-square statistic: To enhance sensitivity to slowly changing faults, a sliding window is used to accumulate multi-step residual information, constructing a dynamic chi-square statistic as follows: r i V represents the residual vector at time step i. i Represents the residual r i The covariance matrix, window length m = [T] resp / Δt]∈[5,15], response time T resp The statistic ∈ [0.5s, 1.5s] follows a chi-square distribution with m degrees of freedom.

[0019] Furthermore, the specific process of S2 is as follows:

[0020] S2.1 Construct a training set and collect historical residual sequences under fault-free conditions. As input, the corresponding chi-square statistic As output, this constitutes the regression training set, where d = dim(r k ) represents the residual dimension (equal to the number of visible satellites), and n≥1000 represents the sample size, ensuring the model's generalization ability;

[0021] S2.2 Define the model structure, using the Support Vector Regression (SVR) model, whose mapping function is: Among them, kernel function The radial basis function (RBF) is used, where γ > 0 is the bandwidth parameter, controlling the decay rate of similarity between samples; α i , Let be the Lagrange multipliers and b be the bias term, which are solved using an optimization algorithm;

[0022] S2.3 Joint parameter optimization, traversing the candidate space of bandwidth parameter γ and penalty factor C through grid search: γ∈{10 -3 10 -2 ,0.1,1,10,10 2 10 3}, C∈{0.1,1,10,100}, this search range covers feature space mappings from sparse to dense, ensuring the robustness of the model in complex noisy environments;

[0023] S2.4 Cross-validation mechanism: Five-fold time series cross-validation is adopted, and the objective function is the weighted sum of mean squared error (MSE) and detection delay time: Score = 0.6 × MSE + 0.4 × T detect The formula for calculating MSE is: n val This indicates the number of samples in the validation set.

[0024] Furthermore, the specific process of S3 is as follows:

[0025] S3.1, Statistical series for SVR prediction Perform 5th-order polynomial fitting: Where a0 is the constant term of the polynomial, representing the value of the fitted curve at t=0; a1 is the coefficient of the first-order term, representing the linear rate of change of the fitted curve; a2 is the coefficient of the second-order term, representing the quadratic curvature of the fitted curve; a3 is the coefficient of the third-order term, representing the cubic curvature of the fitted curve; a4 is the coefficient of the fourth-order term, representing the quartic curvature of the fitted curve; and a5 is the coefficient of the fifth-order term, representing the quintic curvature of the fitted curve. The fifth-order polynomial can effectively approximate the nonlinear trend of the statistic while avoiding overfitting caused by higher-order terms.

[0026] S3.2 Calculate the instantaneous slope and differentiate the fitted curve to obtain the instantaneous slope function: This function quantifies the rate of change of statistics over time, providing a dynamic threshold benchmark for fault detection.

[0027] S3.3 Determine the confidence interval and calculate the mean slope μ based on historical normal data. s With standard deviation σ s Determine the 99% confidence interval: [s min ,s max ] = [μ s -2.576σ s ,μ s +2.576σ s ], where 2.576 is the quantile of the standard normal distribution at the 99% confidence level;

[0028] S3.4 Outlier Removal: During the preprocessing stage, outliers satisfying the condition |s i -μ s |>3σ s Abnormal samples are identified to avoid interference from extreme values ​​in threshold calculation;

[0029] S3.5. Within the sliding window [k-L+1,k], where L≤m, fit the slope of the current statistic using the weighted least squares method: weight w i With statistic variance Inversely proportional, ensuring that high-confidence data dominates the fitting; variance The standard deviation is the square of the standard deviation. Where, n val That is the number of samples in the validation set. These are the observations for each data point. It is the mean of the observed values;

[0030] S3.6, Define a binary event sequence {E} k}: This sequence is used to mark whether the slope at each time point exceeds the normal range; exceeding the normal range is marked as 1, and being within the normal range is marked as 0.

[0031] S3.7 Continuous Exceedance Decision: An alarm is triggered when three consecutive exceedance events are detected, i.e.: This decision condition effectively distinguishes between slowly changing faults and transient disturbances by utilizing the continuous deviation characteristics of the cumulative statistic slope.

[0032] Continuous exceedance determination distinguishes between slowly varying faults and transient disturbances by analyzing the sustained deviation of the cumulative statistic slope. Specifically, the trigger condition for continuous exceedance determination is three consecutive detections of a slope exceeding the normal range, i.e., an exceedance event. This method effectively distinguishes between slowly varying faults and transient disturbances for the following reasons:

[0033] Characteristics of slowly varying faults: Slowly varying faults are usually characterized by a slow and continuous deviation of the slope of the statistical quantity from the normal range. Since this deviation is gradually accumulated, the slope is likely to continue to exceed the normal range in multiple consecutive tests.

[0034] Characteristics of transient disturbances: Transient disturbances usually manifest as a sudden deviation in the slope of a statistic, but this deviation is usually brief and does not last.

[0035] Therefore, in multiple consecutive tests, the slope of the transient interference may exceed the normal range only once or twice, but it is unlikely to exceed it three times in a row.

[0036] By utilizing the continuous deviation characteristics of slowly varying faults and the brief deviation characteristics of instantaneous disturbances, the two can be effectively distinguished. Specifically: Slowly varying faults: due to their continuous deviation in slope, they are likely to exceed the normal range three times in a row, thus triggering an alarm; Instantaneous disturbances: because their slope deviation is usually brief, they are unlikely to exceed the normal range three times in a row, thus reducing the possibility of false alarms.

[0037] S3.8 False Alarm Rate Control, the false alarm rate is calculated as: P false =α 3 =10 -6 α = 0.01, assuming the false alarm probability of a single exceedance is 1% and the joint probability of three consecutive independent exceedances is 10. -6 Then the security requirements are met.

[0038] Under fault-free conditions, the real-time slope s obs Follows a normal distribution: Where, μ s σ is the mean of the historical normal slope. s The standard deviation is given. The 99% confidence interval for the normal slope is set as follows:

[0039] [s min ,s max ] = [μ s -z 0.995 σ s ,μ s +z 0.995 σ s ];

[0040] Among them, z 0.995 =2.576 is the 99th percentile of the standard normal distribution. The one-sided out-of-limit probability is: α = P(s) obs >s max )=P(s obs <s min ) = 0.5%;

[0041] The total probability of exceeding the limit is: α total =1% (i.e., 0.5% + 0.5%);

[0042] According to DO-178C Level A standard (aviation-grade reliability requirements), the maximum permissible false alarm rate for critical systems is:

[0043] P false ≤10 -9 / hour, via α and P false The relationship determines the minimum number of consecutive out-of-limit times N min :

[0044] If α = 1%, P false =10 -6 →N min =3;

[0045] Therefore, for GNSS / INS integrated navigation systems, the combination of α=1% and N=3 can balance real-time performance (detection delay ≤1.5s) and safety (P false =10 -6It meets the requirements of DO-178C Level B.

[0046] This invention avoids the computational burden of high-order features by directly using the residual sequence as model input, and replaces traditional binary classification labels with dynamic chi-square statistics, achieving continuous quantification of fault intensity output and significantly reducing detection latency. Utilizing an SVR detection model based on grid search and cross-validation co-optimization, it designs a parameter search space tailored to the characteristics of slowly varying faults. Through logarithmic uniform sampling, it matches the low-frequency characteristics of slowly varying signals; through interval linear sampling, it balances the model's sensitivity to small residual fluctuations and significant faults; and through five-fold cross-validation, it selects the optimal parameter combination that minimizes the mean square error (MSE), solving the prediction bias problem caused by traditional empirical parameters and reducing the mean absolute error of slowly varying fault detection. This invention exhibits strong generalization performance, meeting the real-time detection requirements of highly dynamic scenarios. Attached Figure Description

[0047] Figure 1 This is a schematic diagram of the workflow of the method of the present invention;

[0048] Figure 2 (a) and (b) in the examples are comparison charts of RMSE of the LSTM model training set and test set, respectively.

[0049] Figure 3 (a) and (b) in the examples are comparison charts of RMSE of the LSSVM model training set and test set, respectively.

[0050] Figure 4 In the examples, (a) and (b) are comparison charts of RMSE of the training set and test set of the model on which the present invention is based in the embodiments.

[0051] Figure 5 In the figure, (a) is the chi-square statistic and threshold plot of the AIME algorithm in the embodiment, and (b) is the normal slope range predicted by the model on which the present invention is based and the actual slope plot. Detailed Implementation

[0052] The invention will now be further described with reference to the accompanying drawings.

[0053] like Figure 1 As shown, a method for detecting slowly varying faults in GNSS positioning based on residual-SVR regression includes the following steps:

[0054] S1. Obtain the observation innovation sequence based on the Kalman filter (EKF) and calculate its covariance matrix; accumulate multiple steps of innovation through a sliding window, construct the chi-square statistic, and quantify the fault distribution characteristics.

[0055] S2. Based on fault-free data, a training set is constructed with the innovation sequence as input and the chi-square statistic as output. The radial basis function is used, and the SVR parameters are optimized by grid search combined with five-fold cross-validation to generate the innovation-statistic mapping function.

[0056] S3. Fit the normal slope threshold based on the SVR predicted value, and fit the observed statistical curve using the real-time sliding window least squares. If the slope exceeds the threshold range for 3 consecutive times, a slow change fault alarm is triggered.

[0057] Example: To verify the superiority of the GNSS positioning gradual change fault detection method based on grid search optimized support vector regression (GS-SVR) proposed in this invention, two sets of systematic experiments were conducted to comprehensively evaluate its model generalization ability and detection timeliness, and a comparative analysis was performed with mainstream algorithms.

[0058] The first set of experiments was used to compare the predictive performance of different regression models. The experiments used historical residual sequences collected by a GNSS receiver under fault-free conditions as input to train LSTM, LSSVM, and the GS-SVR model proposed in this invention, respectively, with the root mean square error (RMSE) as the core evaluation index.

[0059] The LSTM model exhibited a low RMSE (1.8341 m) on the training set, but its error on the test set increased significantly to 3.2687 m, with an overfit rate as high as 78.1%. This indicates that its ability to model the non-stationary temporal characteristics of slowly varying GNSS faults is limited, making it difficult to adapt to dynamic noise interference. Figure 2 As shown in (a) and (b), the RMSE of the LSTM model on the training and test sets differs significantly, revealing its insufficient generalization ability.

[0060] Although the LSSVM model achieved a low error of 0.8339 m on the training set, the RMSE on the test set rose to 3.5809 m, with an overfit rate exceeding 300%, further exposing the inherent flaw of fixed-parameter models in their sensitivity to data distribution. Figure 3 As shown in (a) and (b), the training set error of LSSVM is extremely low, while the test set error fluctuates wildly, verifying that it is heavily dependent on the distribution of the training data.

[0061] This invention, through grid search optimization and five-fold cross-validation, achieved RMSEs of 1.8632m and 1.3544m on the training and test sets respectively, with an overfitting rate of only 27.3%, demonstrating excellent generalization ability. Figure 4 As shown in (a) and (b), the errors of the training and test sets of this invention are both within a reasonable range and the difference is minimal, proving its anti-overfitting properties. This advantage stems from the global optimization of the regularization parameter C and the bandwidth parameter γ by grid search (C∈{0.1, 1, 10, 100}, γ∈{10...-3 10 -2 The model is robust in noisy environments by combining cross-validation with 0.1 and 1 to correct for data bias.

[0062] The second set of experiments focused on the timeliness of fault detection, comparing the response speed of the present invention with that of the traditional Adaptive Iterative Least Squares (AIME) algorithm in slowly varying fault scenarios. The experiments injected slowly varying faults (initial deviation 0.1 m, slope 0.1 m / s, duration 100 seconds) with linearly increasing amplitudes into GNSS pseudorange observations and recorded the first alarm time for both methods.

[0063] The AIME algorithm relies on a threshold judgment mechanism based on the chi-square test statistic, requiring a sufficient accumulation of residual deviations to trigger an alarm. This resulted in it only detecting the fault after 82.6 seconds, by which time the fault amplitude had already accumulated to 5.26 meters, significantly lagging behind the actual fault development process. For example... Figure 5 As shown in (a), the chi-square statistic of AIME (blue curve) requires a long accumulation period to exceed the threshold (gray dashed line).

[0064] This invention presents a GNSS positioning gradually varying fault detection method based on residual-SVR regression. By calculating the slope of the chi-square test statistic curve after the injected fault in real time and comparing it with the slope range predicted based on historical normal data, it detects the slope exceeding the normal range within 62 seconds. At this point, the fault amplitude is only 3.2 meters, 24.9% earlier than AIME. Figure 5 As shown in (b), the normal slope range is defined by the 99% confidence interval of historical data ([s min ,s max ]=μ s ±2.576σ s The actual slope after the injected fault deviated significantly from this range after 62 seconds, verifying the invention's ability to keenly capture gradually changing trends. This difference stems from the real-time sliding window fitting strategy: using a weighted least squares method with a window length of L=5, the instantaneous slope of the current residual sequence is calculated, and early fault identification is achieved by comparing it with the normal range of the prediction model.

[0065] The combined results of the two sets of experiments demonstrate that this invention achieves a comprehensive improvement in GNSS gradually varying fault detection performance by integrating parameter optimization, dynamic thresholding, and real-time slope monitoring. In terms of model generalization, the test set RMSE is reduced by 58.6% and 62.2% compared to LSTM and LSSVM, respectively, effectively overcoming the overfitting problem. Regarding detection timeliness, the alarm response speed is improved by 24.9% compared to AIME. These characteristics make it particularly suitable for highly dynamic GNSS / INS integrated navigation systems, providing accurate and reliable fault perception capabilities for scenarios such as UAVs and autonomous driving. Figures 2 to 4The comparison results of the algorithms are presented intuitively, further demonstrating the effectiveness of the present invention. Experimental data show that the fault detection accuracy and response speed of the present invention are significantly better than traditional methods in complex noise environments, providing a high reliability guarantee for GNSS positioning systems.

[0066] Furthermore, this invention can be used not only for single systems such as GPS, GLONASS, BDS, and GALILEO, but also for combined positioning of different satellite systems and inertial systems, making it widely applicable.

Claims

1. A method for detecting slowly varying faults in GNSS positioning based on residual-SVR regression, characterized in that, Includes the following steps: S1. Obtain the observed information sequence based on the Kalman filter EKF and calculate its covariance matrix; By accumulating multiple steps of information through a sliding window, a chi-square statistic is constructed to quantify the characteristics of fault distribution. S2. Based on fault-free data, a training set is constructed with the innovation sequence as input and the chi-square statistic as output. Radial basis function is used, and the SVR parameters are optimized by grid search combined with five-fold cross-validation to generate the innovation-chi-square statistic mapping function. S3. Fit the normal slope threshold based on the SVR predicted value, and fit the chi-square statistic curve of the observed value using the real-time sliding window least squares. If the slope exceeds the threshold range multiple times in a row, a slow-change fault alarm will be triggered. The specific process of S3 is as follows: S3.1, Chi-square statistic series for SVR prediction Perform 5th-order polynomial fitting: ;in, It is the constant term of the polynomial, representing the value of the fitted curve at t=0. It is the coefficient of the first-order term, representing the linear rate of change of the fitted curve. These are the coefficients of the quadratic term, representing the quadratic curvature of the fitted curve. These are the coefficients of the cubic term, representing the cubic curvature of the fitted curve. These are the coefficients of the quartic term, representing the quartic curvature of the fitted curve. These are the coefficients of the fifth-order term, representing the fifth-order curvature of the fitted curve; S3.2 Calculate the instantaneous slope and differentiate the fitted curve to obtain the chi-square statistic instantaneous slope function of the SVR prediction: ; S3.3 Determine the confidence interval and calculate the mean slope based on historical normal data. with standard deviation Determine the confidence interval: Where 2.576 is the quantile of the standard normal distribution at the confidence level; S3.4 Outlier Removal: During the preprocessing stage, outliers that meet the following conditions are removed. Abnormal samples are identified to avoid interference from extreme values ​​in threshold calculation; S3.5, in the sliding window We use the weighted least squares method to fit the slope of the chi-square statistic of the current observations: Weight Variance of chi-square statistic Inversely proportional, ensuring that high-confidence data dominates the fitting; variance The standard deviation is the square of the standard deviation. ,in, That is the number of samples in the validation set. It is the chi-square statistic of the observed values ​​for each data point. It is the mean of the chi-square statistic of the observed values; S3.6, Define a binary event sequence : This sequence is used to mark whether the slope of the chi-square statistic of the observation at each time point exceeds the normal range; exceeding the normal range is marked as 1, and being within the normal range is marked as 0. S3.7 Continuous Exceedance Decision: An alarm is triggered when three consecutive exceedance events are detected, i.e.: The decision to trigger an alarm after multiple consecutive over-limits effectively distinguishes between slow-changing faults and instantaneous disturbances by utilizing the continuous deviation characteristics of the cumulative chi-square statistic slope. S3.8 False Alarm Rate Control, the false alarm rate is calculated as follows: Assuming the false alarm probability of a single exceedance is 1%, the joint probability of three consecutive independent exceedances is... Then the security requirements are met.

2. The GNSS positioning slowly varying fault detection method based on residual-SVR regression according to claim 1, characterized in that, The specific process of S1 is as follows: S1.1 Define the state vector of the GNSS positioning system. The state vector contains position error components corresponding to the three orthogonal directions of East, North, and Sky. and velocity error components and receiver clock deviation and its drift rate ; The dynamic evolution model of clock skew characterizes the cumulative error characteristics of the receiver clock, and its expression is: Where k represents the discrete-time index of the current moment. This represents the clock drift at the previous time k-1. For the filter period, For clock noise that follows a Gaussian distribution, This represents the noise variance, reflecting the stability of the clock crystal oscillator. S1.2, Define the state transition matrix The position and velocity block matrices are used to predict the state vector at the next moment. ; S1.3 Calculate the actual pseudorange observations based on the Extended Kalman Filter (EKF) framework. Compared with the predicted value residual , , This is the actual pseudorange observation vector. , This is the pseudorange measurement value of the Mth visible satellite; To predict the state vector, , , , For receiver location, , , For receiver speed, For receiver clock deviation, For receiver drift rate; The observation matrix; the observation matrix Each row corresponds to the geometric relationship of a visible satellite, and its elements are determined by the geometric partial derivatives of the satellite's position with respect to the receiver's position. The expression is as follows: ,in, This represents the partial derivative of the pseudorange of the i-th satellite with respect to the receiver's three-dimensional position. 1 and 0 correspond to the observation coefficients of clock bias and drift rate, respectively. S1.4 Update the residual covariance matrix, the residual covariance matrix From the predicted covariance matrix With observation noise matrix The joint calculation uses the following update formula: , among which, among which The state prediction covariance matrix, This is the state estimation covariance matrix at time k-1, representing the confidence level of the system state at the previous time step. In the calculation of the state prediction covariance: Φ represents the state transition matrix. Represents the process noise covariance; The process noise covariance matrix is... and These represent the process noise variances for position and velocity, respectively, reflecting the uncertainty of the system model; The process noise variance of the receiver clock skew. The process noise variance is the clock drift rate; Observation noise matrix It is a diagonal matrix, and its elements are... Indicates the first The pseudorange measurement noise variance of each satellite is calculated in real time by the satellite signal quality assessment module; S1.

5. Constructing a sliding window chi-square statistic: To enhance sensitivity to slowly changing faults, a sliding window is used to accumulate multi-step residual information, constructing a dynamic chi-square statistic as follows: , This represents the residual vector at time step i. Residual The covariance matrix, window length Response time The chi-square statistic follows a chi-square distribution with m degrees of freedom.

3. The GNSS positioning slowly varying fault detection method based on residual-SVR regression according to claim 1, characterized in that, The specific process of S2 is as follows: S2.1 Construct a training set and collect historical residual sequences under fault-free conditions. As input, the corresponding chi-square statistic As output, this constitutes the regression training set, where... The residual dimension is equal to the number of visible satellites; For sample size; S2.2 Define the model structure, using the Support Vector Regression (SVR) model, whose mapping function is: Among them, kernel function Radial basis function (RBF) This is a bandwidth parameter that controls the decay rate of similarity between samples. , For Lagrange multipliers, The bias term is solved using an optimization algorithm; S2.3 Joint parameter optimization, traversing bandwidth parameters through grid search. With penalty factor Candidate space: The candidate space covers a feature space mapping from sparse to dense. S2.4 Cross-validation mechanism: Five-fold time series cross-validation is adopted, and the objective function is the weighted sum of mean squared error (MSE) and detection delay time. The formula for calculating MSE is: , This indicates the number of samples in the validation set.

Images