Factor graph optimization train positioning method for resisting satellite navigation spoofing attack

Abstract

The invention provides a factor graph optimization train positioning method for resisting satellite navigation spoofing attacks. The method comprises the following steps: determining a factor graph model with deception mitigation capability and related parameters; when the information frame of the speed and distance measuring sensor is received, calculating the pre-integration position of the speed and distance measuring sensor; determining a nearest orbit slice in the orbit information, sensing a deception event and identifying a deception satellite in the visual satellites; constructing a self-adaptive pseudo-range factor of a real / deception satellite, a pre-integration factor of a speed and distance measuring sensor and a combined optimization factor graph, and solving a positioning solution of the train; and calculating a pseudo-range residual error according to the positioning solution, compensating the residual error of the deception satellite, updating a residual error pool, fitting residual error distribution of each visual satellite in the residual error pool, calculating an adaptive weight, and continuing iteration to obtain an optimal positioning solution of the train. According to the method, deception attack detection, identification and mitigation are realized based on a factor graph optimization framework, the real operation state of the train before attack is recovered, and the deception defense capability is enhanced.

Description

Technical Field

[0001] The present invention relates to the technical field of rail transit train positioning, and in particular to a factor graph optimization train positioning method for resisting satellite navigation deception attacks. Background Art

[0002] Train positioning technology enabled by the Global Navigation Satellite System (GNSS) has accelerated the transition from traditional train control systems (TCS) to a new generation of systems with higher operating speeds, lower construction and maintenance costs, and greater intelligence. With the introduction of GNSS, its vulnerability has exposed train positioning to unintentional interference caused by environmental obstructions along the railway and radio frequency interference (RFI) caused by malicious attackers. In particular, GNSS spoofing attacks with active induction properties can manipulate the position, velocity, and time (PVT) of train receivers without alerting the victim train, seriously threatening train safety. A single GNSS system is no longer able to meet the high-availability, high-precision, and high-reliability positioning requirements of safety-critical train control systems in complex interference environments. To improve the continuity, accuracy, and robustness of GNSS-based train positioning, the use of multi-source fusion positioning algorithms to fuse measurement information from multiple GNSS-independent sensors is expected to achieve multi-source trusted positioning of trains. Combining multi-source information for deception immunity can enhance train positioning's ability to recover the true PVT solution under deception attacks, further enabling multi-source trusted train positioning in deception attack scenarios. Existing anti-spoofing technologies primarily focus on deception detection, aiming to promptly detect deception attacks and trigger early warning mechanisms. Deception mitigation, a key component of further responding to deception attacks to mitigate their impact on positioning solutions, has been a relatively recent development. In anti-spoofing technologies, state estimation serves as a back-end optimization method to estimate and recover the state of a deceived system. Most existing solutions still rely on least squares estimation and filter fusion frameworks. The emerging factor graph optimization architecture, as a global optimization solver, is gaining increasing attention and has been demonstrated to outperform filter fusion algorithms in positioning accuracy and robustness. Therefore, exploring anti-spoofing technologies that integrate deception detection, mitigation, and state estimation based on factor graph optimization architectures holds promise for not only achieving more robust response strategies against deception attacks but also potentially achieving state recovery capabilities superior to those of filtering methods under attack.

[0003] Under malicious spoofing attacks, the performance of GNSS-based train positioning may degrade or even fail, making it difficult to meet the needs of trusted positioning. In order to cope with the threat of GNSS spoofing attacks, there is an urgent need to combine the existing on-board positioning sensors, advanced anti-spoofing technologies and positioning solution technologies of the train to detect and mitigate spoofing attacks, so as to trigger safety warnings in a timely manner and accurately restore the true state of the train, and realize multi-source trusted positioning of the train under spoofing attacks. The present invention aims to alleviate the impact of spoofing attacks on the satellite positioning of the train, make full use of the spoofing-immune multi-source auxiliary information provided by the existing configuration of the train and the commonly used sensor equipment, integrate spoofing detection, identification and mitigation technologies into a multi-source fusion positioning architecture based on factor graph optimization, and thus simultaneously achieve the mitigation of spoofing attacks and the restoration of the true state of the train.

[0004] Currently, existing anti-spoofing technologies based on measurement domain information processing generally employ two approaches: "identify and eliminate" and "detect and exploit" to mitigate the impact of spoofing attacks on PVT solutions. The "identify and eliminate" countermeasure involves identifying a specific satellite affected by a spoofing attack and then directly eliminating the identified spoofed measurements from the positioning solution to isolate the attack's impact. Fault Detection Exclusion (FDE) and Receiver Autonomous Integrity Monitoring (RAIM) technologies have proven effective in isolating spoofing attacks. However, the "identify and eliminate" approach may result in situations where the remaining true satellite measurements do not meet the minimum required number for positioning, thereby reducing the availability of satellite positioning. Related research has shown that spoofed satellite information can help recover true PVT information. This has led to the development of "detect and exploit" countermeasures. After detecting a spoofing attack, these countermeasures are exploited and mitigated through online weight adjustment or error compensation. Most existing anti-spoofing methods, such as robust estimation and adaptive Kalman filtering, embody this principle. The spoofing envelope is defined as the deviation in satellite measurements caused by spoofing attacks. Previous studies have explored anti-spoofing capabilities by modeling, estimating, and compensating for this spoofing envelope using a filtered positioning framework. Compared to the "identify-and-remove" approach, the "detect-and-exploit" approach retains information about the spoofed satellite, ensuring both availability and continuity while also improving positioning accuracy. Among these, the online weight adjustment strategy, which does not rely on a precise error correction model, is the clear choice for anti-spoofing techniques based on measurement-domain information processing. The error correction model can further restore the true state under spoofing attack conditions, improving positioning accuracy.

[0005] A train positioning method in the prior art for countering satellite navigation spoofing attacks includes: a GNSS / INS tight coupling method based on extended Kalman filtering is used to counter spoofing attacks, which uses Kalman filter information and a spoofing envelope estimated in an iterative process to perform spoofing identification.

[0006] Another existing train positioning method to combat satellite navigation spoofing attacks integrates sparse estimation theory into the extended Kalman filter, transforms the GNSS state estimation problem into the L1 norm, and uses the least absolute convergence and selection algorithm (LASSO) for sparse estimation, thereby alleviating the impact of spoofing pseudorange on positioning solutions while estimating the state.

[0007] The disadvantages of the above-mentioned prior art train positioning method for resisting satellite navigation spoofing attacks include:

[0008] (1) Most of the existing deception mitigation technologies still use filter estimation as a state solver. Although a few scholars have tried to realize the coordinated processing of deception detection and state estimation within the filtering framework, overall there are still problems with insufficient robustness and recovery capabilities. Factor graph optimization, as a state estimation method that has developed rapidly in recent years, has provided a new solution for the development of anti-deception technology with its advantages in information modeling flexibility, multi-source observation fusion and nonlinear optimization. Combining factor graph optimization with deception mitigation strategies is expected to significantly improve the state recovery capability in the face of GNSS deception attacks. However, the current research on deception mitigation technology based on factor graph optimization is still relatively limited, especially in the construction of a unified optimization model that can simultaneously achieve deception mitigation and robust state estimation, which needs further in-depth exploration.

[0009] (2) Most existing spoofing mitigation methods are oriented towards general satellite positioning scenarios and lack specialized designs tailored to the characteristics of train control systems, making them difficult to directly integrate and apply to existing train operation control systems. In train control systems, existing non-GNSS sensors are not easily affected by spoofing attacks and can provide redundant and reliable auxiliary information for positioning. Therefore, there is an urgent need to develop specialized spoofing interference mitigation and state recovery methods for train satellite positioning scenarios, making full use of the train's unique multi-sensor and track prior information to achieve highly robust and reliable spoofing protection capabilities.

[0010] (3) Most existing anti-spoofing technologies use filter estimation as a state solution method. However, the factor graph optimization algorithm that has emerged in recent years has shown many advantages in modeling flexibility, information fusion efficiency and scalability, and is expected to achieve anti-spoofing performance that is superior to traditional filters. Despite this, the anti-spoofing potential of the factor graph optimization framework has not yet been fully explored. In particular, for special application scenarios such as train operation control based on satellite positioning, there is still a lack of dedicated deception interference mitigation methods designed for train positioning characteristics. Summary of the Invention

[0011] An embodiment of the present invention provides a factor graph optimization train positioning method for resisting satellite navigation spoofing attacks, so as to effectively ensure the credibility of train positioning.

[0012] In order to achieve the above-mentioned purpose, the present invention adopts the following technical solutions.

[0013] A factor graph optimization train positioning method to combat satellite navigation spoofing attacks includes:

[0014] determining and initializing a system state node having a deception envelope component, a factor graph model having deception mitigation capabilities, and related parameters;

[0015] When receiving a speed and distance measuring sensor information frame, or when receiving a satellite signal information frame, after aligning the timestamps of the speed and distance measuring sensor information frame using linear interpolation, a pre-integrated position of the speed and distance measuring sensor is calculated;

[0016] Determine the nearest orbital slice from the orbital information, map-match the pre-integrated position, calculate reliable deviations, construct adaptive detection and recognition statistics, sense spoofing events, and identify spoofed satellites among visible satellites;

[0017] Activate the deception envelope state corresponding to the identified spoofed satellite, construct the adaptive pseudorange factors of the real / spoofed satellite, the pre-integration factors of the speed and range sensors, the map matching factors and the prior factors, jointly optimize the factor graph, and preliminarily solve the train positioning solution;

[0018] Based on the preliminary positioning solution, the pseudorange residuals are calculated and the residuals of the deceptive satellites are compensated. The residual pool is updated and the residual distribution of each visible satellite in the residual pool is fitted. The adaptive weights are calculated and fed back into the adaptive deceptive detection, identification, and mitigation steps to continue iterating to obtain the optimal positioning solution for the train.

[0019] Preferably, the determining and initializing of the system state node having the deception envelope component, the factor graph model having the deception mitigation capability, and related parameters includes:

[0020] S1.1. Determine and initialize the system state node to be estimated with the deceptive envelope component;

[0021] In the multi-source tightly coupled positioning, the raw observed pseudorange and Doppler shift of GNSS, the raw observed values of INS three-axis gyroscope and accelerometer, and the raw observed wheel axis pulse count of ODO are processed and sent to the factor graph optimization for joint optimization positioning. If there is a deceptive satellite among the visible satellites, the deceptive envelope state component of the deceptive satellite will be activated. The east-north-sky world coordinate system (w system) with the initial position as the origin is selected, which is also the reference coordinate system of the factor graph. Then, within a time window, the system's multidimensional state node set χ and state node x are k Build as:

[0022]

[0023] Among them, x k is the 18-dimensional state node at the k-th epoch, and is the three-dimensional position and velocity of the train, is the train attitude expressed in quaternion form, and The three-dimensional accelerometer and gyroscope zero bias of the IMU (Inertial Measurement Unit) system, and is the clock error and clock drift of the GNSS receiver, is the scale factor of ODO, s k is the deception envelope, m is the number of identified deception satellites, and n is the sliding window size;

[0024] S1.2. Determine and initialize the prior factor of the global first node of the deception mitigation factor graph model

[0025] The initial state prior information x0 is given by the measurement value or experience value of the measuring device and corresponds one-to-one with the system state node, which can be expressed as:

[0026]

[0027] in, To initialize position, velocity and attitude, and To initialize the IMU accelerometer and gyroscope bias, and To initialize the GNSS receiver clock error and drift, To initialize the scale factor of ODO;

[0028] Assuming that the error w0 of the initialization information obeys the Gaussian distribution N, the measurement equation of the prior factor of the global first node is expressed as:

[0029] x=x0+w0,w0~N(0,Σ0) (3)

[0030] Among them, Σ0 is the noise covariance matrix of the prior factor of the first global node, which represents the uncertainty of the initialization information;

[0031] Get the error function of the prior factor of the global first node

[0032]

[0033] Calibrate the factor graph model based on the error function of the prior factor of the global first node;

[0034] S1.3. Determine and initialize a factor graph model with deception mitigation capabilities;

[0035] The inertial sensor IMU measurements and the mileage measurements of the train wheels provided by the ODO are combined. During the optimization process, IMU / ODO pre-integration is used as the main body for time series state recursion. The IMU / ODO pre-integration factor constructed based on the IMU and ODO measurement models is used to form probabilistic constraints on the states of adjacent epochs. The GNSS pseudorange factor and Doppler velocity factor, based on the observations of visible satellites, construct a measurement model to constrain the relevant state nodes of the current epoch. The map matching factor constructed with the assistance of the electronic track map will constrain the train position to the inherent track line.

[0036] Deception countermeasures are incorporated into the factor graph structure of tightly coupled GNSS / INS / ODO / DTM positioning, including adaptive deception detection, identification, and mitigation. IMU / ODO pre-integration outputs, map matching, and GMM error weights are used to assist in constructing reliable biases. Binary hypothesis testing is performed to determine whether a deception attack has occurred. After a deception attack is detected, the deceptive satellites are identified among the visible satellites by normalizing the reliable biases. An adaptive deception elimination strategy activates the deception envelope nodes of the identified deceptive satellites, directly optimizing and compensating for the deceptive measurements during the optimization process, and determining and initializing a factor graph model with deception mitigation capabilities.

[0037] Preferably, the method of calculating the pre-integrated position of the speed and ranging sensor when receiving the speed and ranging sensor information frame or aligning the timestamps of the speed and ranging sensor information frame using linear interpolation when receiving the satellite signal information frame includes:

[0038] Set up speed and distance measurement sensors including inertial sensor IMU and axle speed and distance measurement sensor ODO;

[0039] The angular velocity of the IMU and acceleration Modeled as:

[0040]

[0041] Among them, n g and n a are the noise of the gyroscope and accelerometer, respectively, is the direction cosine matrix of the transformation from the w system to the IMU carrier system (b system), For its inverse transform, is the projection of the Earth's rotation angular velocity in the Earth-centered Earth-fixed coordinate system (e system) relative to the inertial coordinate system (i system) in the w system, w is the angular velocity caused by the carrier motion and the curvature of the earth, and is the Coriolis acceleration and centripetal acceleration caused by the earth's rotation and the carrier's motion, is the navigation system (n system) earth gravity coordinate transformation matrix The projection of the earth's gravity in the w system is obtained. The gravity at different locations is usually related to the dimension in which it is located. Related to the height h, set the gravity model to:

[0042]

[0043] According to the pulse count of ODO Calculate the one-dimensional mileage increment of the train along the track between two consecutive time stamps

[0044]

[0045] Among them, R w is the radius of the driving wheel;

[0046] The ODO measurement is expressed in the following vector form:

[0047]

[0048] The pre-integration model of IMU / ODO at two adjacent optimization moments is:

[0049]

[0050] in, is the increment of ODO in system b, is the velocity in frame b, ι a , ι g and ι odo is the Gaussian white noise of the IMU accelerometer, IMU gyroscope and ODO scale factor, Ω is a synthetic matrix whose first column is the quaternion The first line is The transpose of the conjugate quaternion of The antisymmetric matrix formed.

[0051] Preferably, the steps of determining the nearest orbital slice in the orbital information, performing map matching on the pre-integrated position, calculating reliable deviations, constructing adaptive detection and recognition statistics, sensing spoofing events, and identifying spoofed satellites among visible satellites include:

[0052] S3.1 performs map matching to determine the on-orbit position of the pre-integrated predicted position of the speed and distance measuring sensor;

[0053] Extract the track segment currently occupied by the train and simplify the track electronic map into a series of points of interest under the e-system The adjacent points of interest are the endpoint coordinates of the track segment, and the positions are predicted by judging the IMU / ODO pre-integration The track segment occupied by the train is determined by its proximity to the track segment;

[0054]

[0055] Where id represents the index of the track segment;

[0056] After confirming the track segment index occupied by the train, the train predicted position Projected onto the track line, train projection position Calculated as:

[0057]

[0058] S3.2 Calculate reliable biases to sense spoofing events and identify spoofed satellites among visible satellites;

[0059] The reference pseudorange of the i-th visible satellite is calculated as:

[0060]

[0061] Reference pseudorange and satellite-corrected observed pseudorange The difference is defined as the reliability deviation and is calculated as:

[0062]

[0063] For real satellites, the reliable biases only contain the unmodeled pseudorange errors ε k,i , and when a deception event occurs, the deception envelope s will be introduced into the deception satellite bias k,i , based on the reliable deviation, a weighted residual sum of squares is constructed as the deception detection statistic;

[0064]

[0065] Where l is the number of visible satellites in the current epoch, W kis an adaptive weighting matrix whose diagonal elements are the inverse of the satellite residual distribution variance derived by the Gaussian mixture model. Under no deception conditions, the WSSE detection index obeys the chi-square distribution with a degree of freedom of l. However, under deception attacks, the WSSE detection index will no longer obey the chi-square distribution. The detection threshold T is calculated based on the pre-set false alarm rate. h , then the binary hypothesis test for adaptive deception detection is:

[0066]

[0067] The carrier noise power spectral density C / N0 of the satellite signal is introduced into the weighting matrix of the detection statistic, and the C / N0 enhanced weight matrix is constructed as follows:

[0068]

[0069] Among them, Υ(ω1≤Υ≤1) is a quadratic function derived from the visible satellite C / N0, about the center Symmetric; ω1 and ω2 are adjustment factors used to control the minimum value and gradient of the function;

[0070] After a spoofing event is detected, adaptive spoofing identification is further performed to locate the spoofed satellite among the visible satellites. The normalized reliable deviation is calculated as:

[0071]

[0072] In the visible satellite set, the satellite with the largest normalized reliable deviation is judged as a deceptive satellite, and its deceptive flag is set to Tag = 1. k and the adaptive weighting matrix W k The component corresponding to the deception signal in , to update α k and W k , and recalculate the WSSE detection index, and perform deception detection again until the updated WSSE index is lower than the detection threshold. The algorithm terminates and it is considered that all deceptive satellites have been identified. The remaining signals are judged to be real signals, and their deception flags are set to Tag = 0.

[0073] Preferably, the activation and identification of the spoofing envelope state corresponding to the spoofed satellite, the construction of the adaptive pseudo-range factor of the real / spoofed satellite, the pre-integration factor of the speed and distance sensor, the map matching factor and the priori factor, the joint optimization of the factor graph, and the preliminary solution of the train positioning solution include:

[0074] S4.1 constructs a GNSS adaptive factor and adds it to the factor graph model;

[0075] The GNSS factor is used to constrain the position, velocity, clock error and clock drift of the system state node in the current epoch. In the non-deception scenario, the pseudorange of the i-th visible satellite is and Doppler shift d k,i Modeled as:

[0076]

[0077] Among them, the superscript A represents the real satellite information. and is the position and velocity of the receiver in the e frame, and is the spatial position, velocity, clock error and clock drift of the visible satellite, λ is the wavelength, is the line-of-sight vector between the receiver and the satellite, and is the equivalent pseudorange error caused by the ionosphere, troposphere and Earth rotation, and is the corresponding equivalent pseudorange rate error, and are the unmodeled pseudorange and pseudorange rate errors;

[0078] Under the deception attack, the pseudo-range observation equation including the deception is modeled as:

[0079]

[0080] Among them, the superscript S represents the information of deceiving satellites;

[0081] The error function of GNSS pseudorange and Doppler velocity factor is calculated as:

[0082]

[0083] in, and is the GNSS measurement corrected by equations (18), (19) and (20), and is the covariance matrix corresponding to the true pseudorange, pseudorange rate and deceptive pseudorange, which is modeled as a zero-mean Gaussian distribution;

[0084] The GNSS adaptive factor with adaptive error modeling is introduced for factor graph optimization. The GNSS adaptive factor is calculated as:

[0085]

[0086] in, and Adaptive variance derived for Gaussian mixture model GMM;

[0087] S4.2 Given the update period ΔT of the IMU / ODO pre-integration factor, calculate the error function of the IMU / ODO pre-integration factor through the state transition of the error term

[0088]

[0089] in, IMU / ODO measurement IMU / ODO The covariance matrix of represents the uncertainty associated with the IMU / ODO pre-integrated equivalent measurement, and The IMU / ODO pre-integration result corrected by the first-order approximation of formula (9);

[0090] S4.3 Map matching factor is calculated as:

[0091]

[0092] in, is the variance associated with the map matching factor;

[0093] S4.4 constructs a priori factors within the sliding window and adds the prior factors to the factor graph model;

[0094] When the number of states in the sliding window exceeds the window capacity, the oldest state will be marginalized, and the sensor measurements corresponding to the marginalized states will be converted into prior factors. The cost function is:

[0095]

[0096] in, is the fixed linearization point of the first system state node in the sliding window, J prior and σ prior is the Jacobian matrix of the prior factors and the residual;

[0097] S4.5 Preliminary estimation of train running status by combining GNSS true pseudorange adaptive factor, GNSS spoofed pseudorange adaptive factor, GNSS Doppler velocity factor, IMU / ODO pre-integration factor, map matching factor and prior factor By minimizing the global error function:

[0098]

[0099] The factor graph solution method used is the Levenberg-Marquardt method in Ceres Solver.

[0100] Preferably, the method of calculating pseudorange residuals based on the preliminary positioning solution, compensating for the residuals of the deceptive satellites, updating the residual pool, fitting the residual distribution of each visible satellite in the residual pool, calculating the adaptive weights and feeding them back to the adaptive deceptive detection, identification and mitigation steps, and continuing to iterate to obtain the optimal positioning solution for the train includes:

[0101] A Gaussian mixture model with up to three Gaussian components is used to model the core and possible tails of the GNSS residual distribution. According to the optimal state of the factor graph joint optimization solution, the pseudorange residual of the i-th satellite is calculated as:

[0102]

[0103] in, and To solve for the optimal position, receiver clock error and spoofing envelope, for the identified spoofed satellite, its residual will be corrected using the estimated spoofing envelope;

[0104] The most recently calculated GNSS residuals are added to the residual pool. When the residual items in the residual pool exceed the capacity, the earliest added residual items will be removed. GMM is used to fit the residual set of each visible satellite. Since the residuals of spoofed satellites are close to 0 after being corrected by the spoofing envelope, the residuals of identified spoofed satellites will not be added to the residual pool.

[0105] The probability density function of GNSS residuals is modeled as:

[0106]

[0107] in, are the parameters of the Gaussian mixture model, μ and Σ are the weight, mean, and variance of the Gaussian components;

[0108] If variational inference is used to solve the parameters of the Gaussian mixture model, the optimal GMM parameters are obtained by iteratively performing the variational E step and variational M step until the likelihood function converges or reaches a predefined maximum iteration threshold. Then the adaptive variance of satellite pseudorange observation is calculated as:

[0109]

[0110] In the iterative positioning solution, the adaptive variance will be fed back to the adaptive deception detection, identification and mitigation, and the optimal estimate of the train operation status will be obtained by continuous iteration (29).

[0111] As can be seen from the technical solutions provided by the aforementioned embodiments of the present invention, this invention establishes a spoofing interference mitigation mechanism tailored to the actual conditions and characteristics of train satellite positioning, ensuring the credibility of train positioning and the security of train control decisions. This invention aims to mitigate the threat of spoofing interference in train satellite positioning by constructing a factor graph optimization framework to detect, identify, and mitigate spoofing attacks, while also restoring the train's true operating state before the attack and enhancing spoofing defense capabilities.

[0112] Additional aspects and advantages of the present invention will be set forth in part in the following description, will become apparent from the following description, or may be learned by practice of the present invention. BRIEF DESCRIPTION OF THE DRAWINGS

[0113] In order to more clearly illustrate the technical solutions of the embodiments of the present invention, the following briefly introduces the drawings required for use in the description of the embodiments. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying any creative work.

[0114] Figure 1 A flow chart of a factor graph optimization train positioning method for counteracting satellite navigation spoofing attacks provided by an embodiment of the present invention;

[0115] Figure 2 A schematic diagram of a possible implementation of a train multi-source fusion factor graph optimization model for countering GNSS spoofing attacks provided by an embodiment of the present invention;

[0116] Figure 3 A schematic diagram of a possible adaptive deception detection result provided by an embodiment of the present invention

[0117] Figure 4 A schematic diagram of a possible adaptive deception identification result provided by an embodiment of the present invention;

[0118] Figure 5 A schematic diagram of a possible adaptive deception mitigation result provided by an embodiment of the present invention. DETAILED DESCRIPTION

[0119] The embodiments of the present invention are described in detail below, examples of which are shown in the accompanying drawings, wherein the same or similar reference numerals throughout represent the same or similar elements or elements having the same or similar functions. The embodiments described below with reference to the accompanying drawings are exemplary and are only used to explain the present invention, and are not to be construed as limiting the present invention.

[0120] It will be understood by those skilled in the art that, unless expressly stated otherwise, the singular forms "a", "an", "said" and "the" used herein may also include the plural forms. It should be further understood that the term "comprising" used in the description of the present invention refers to the presence of the features, integers, steps, operations, elements and / or components, but does not exclude the presence or addition of one or more other features, integers, steps, operations, elements, components and / or groups thereof. It should be understood that when we refer to an element as being "connected" or "coupled" to another element, it may be directly connected or coupled to the other element, or there may be intermediate elements. In addition, "connected" or "coupled" as used herein may include wireless connections or couplings. The term "and / or" used herein includes any unit and all combinations of one or more associated listed items.

[0121] It will be understood by those skilled in the art that, unless otherwise defined, all terms (including technical and scientific terms) used herein have the same meaning as commonly understood by those skilled in the art in the art to which the present invention pertains. It should also be understood that terms such as those defined in common dictionaries should be understood to have meanings consistent with their meanings in the context of the prior art and, unless defined as such herein, will not be interpreted in an idealized or overly formal sense.

[0122] To facilitate understanding of the embodiments of the present invention, several specific embodiments will be further explained below with reference to the accompanying drawings, and each embodiment does not constitute a limitation on the embodiments of the present invention.

[0123] This invention provides a method for implementing spoofing resistance in the GNSS measurement domain. This method aims to mitigate the impact of spoofing attacks on train satellite positioning. By fully leveraging the spoofing-immune multi-source auxiliary information provided by existing train configurations and commonly used sensor devices, the invention integrates spoofing detection, identification, and mitigation technologies into a multi-source fusion positioning architecture based on factor graph optimization, thereby simultaneously mitigating spoofing attacks and restoring the train's true state.

[0124] This paper adopts a "detection-exploitation approach" based on deception error compensation. Unlike existing methods, this paper further integrates deception detection, identification, and mitigation into a factor graph-based train multi-source fusion architecture to achieve coordinated processing of deception mitigation and true train state recovery.

[0125] This embodiment provides an application scenario for a factor graph optimized train positioning method for countering satellite navigation spoofing attacks. Upon receiving satellite signal information frames, linear interpolation is performed on typical speed and ranging sensor inertial measurement unit (IMU) / axle speed and ranging sensor (ODO) information frames to align timestamps, and IMU / ODO pre-integration is performed. Reliable biases for adaptive spoofing detection and identification are constructed using the IMU / ODO pre-integrated positions and track information provided by the electronic track map (DTM). Hypothesis testing is then performed to detect the occurrence of spoofing events and identify spoofed satellites among visible satellites. If a spoofing event occurs, the spoofing envelope state components corresponding to the identified spoofed satellites in the factor graph system state nodes are activated. Factors are constructed and added to the factor graph model. The factor graph automatically estimates and compensates for the spoofing envelope of the spoofed measurements during joint optimization, thereby restoring the train's true state. Pseudorange residuals are calculated based on the estimated optimal position and updated to a residual pool. A Gaussian mixture model (GMM) is used to fit the pseudorange residual distribution, thereby assisting in adaptive spoofing detection, identification, and suppression during the system cycle.

[0126] The present invention aims to propose a factor graph optimization train positioning method to resist satellite navigation spoofing attacks, so as to timely perceive spoofing attacks and accurately restore the true state of the train, and realize multi-source trusted positioning of the train under spoofing attack scenarios. The flowchart of the factor graph optimization train positioning method to resist satellite navigation spoofing attacks provided by the embodiment of the present invention is as follows: Figure 1 As shown, the processing steps include the following:

[0127] Step S1: Determine and initialize a system state node with a deception envelope component, a factor graph model with deception mitigation capabilities, and related parameters.

[0128] The factor graph model is primarily used in step S4 to achieve spoofing mitigation and optimal estimation of the train's operating status. Based on anti-spoofing technology requirements and sensor fusion positioning methods, the system state nodes and factor graph model under the spoofing attack state are determined. The prior factors and related parameters of the state nodes to be estimated, the global first node, and related parameters are initialized. The resulting factor graph model is then used for subsequent spoofing mitigation and optimal estimation of the train's operating status.

[0129] S1.1. Determine and initialize the system state node to be estimated with the deceptive envelope component.

[0130] A sliding window factor graph model is used to balance the computational cost and real-time performance of graph optimization. The fused sensors include the GNSS (Global Navigation Satellite System), INS (Inertial Navigation System), Odometer (Odometer), and DTM (Digital Track Map), using a tightly coupled approach for high accuracy. In this multi-source tightly coupled positioning, the raw observed pseudoranges and Doppler shifts from the GNSS, the raw observations from the INS triaxial gyroscope and accelerometer, and the raw observed wheel axle pulse counts from the ODO are processed and fed into a factor graph optimization for joint positioning optimization. Typical states to be estimated include position, velocity, attitude, and the GNSS receiver's clock error and drift. The estimated GNSS receiver clock error and drift are used to correct for errors in pseudoranges and pseudorange rates (calculated from Doppler shift). Taking into account the error accumulation characteristics of the inertial measurement unit (IMU), the zero bias of the IMU gyroscope and accelerometer also needs to be considered and added to the state to be estimated to provide error correction. Similarly, due to problems such as actual train wheel radius loss and slippage, deviations are also introduced in the original ODO measurement, and its scale parameters are also added to the state to be estimated. If there is a deceptive satellite among the visible satellites, the deceptive envelope state component of the deceptive satellite will be activated to compensate for the pseudorange of the deceptive satellite. To facilitate multi-sensor fusion and coordinate transformation, the east-north-sky with the initial position as the origin is selected as the world coordinate system (w system), which is also the reference coordinate system of the factor graph. Then, within a time window, the multidimensional state node set χ of the system can be constructed as:

[0131]

[0132] Among them, x k is the 18-dimensional state node at the k-th epoch, and is the three-dimensional position and velocity of the train, is the train attitude expressed in quaternion form, and is the three-dimensional accelerometer and gyroscope zero bias of the IMU carrier system (system b), and is the clock error and clock drift of the GNSS receiver, is the scale factor of ODO, s k is the deception envelope, m is the number of identified deception satellites, and n is the sliding window size.

[0133] S1.2. Determine and initialize the prior factor of the global first node of the deception mitigation factor graph model.

[0134] The prior factor of the first global node determines the initial value of the entire system and is used to align the initial state of the system, directly affecting the convergence and convergence speed of subsequent optimization. A reliable and accurate initial state prior information x0 can be given by the measurement value or empirical value of high-precision measurement equipment, and has a one-to-one correspondence with the system state node, which can be expressed as:

[0135]

[0136] in, To initialize position, velocity and attitude, and To initialize the IMU accelerometer and gyroscope bias, and To initialize the GNSS receiver clock error and drift, The scale factor for initializing ODO.

[0137] Assuming that the error w0 of the initialization information obeys the Gaussian distribution N, the measurement equation of the prior factor of the global first node can be expressed as

[0138] x=x0+w0,w0~N(0,Σ0) (3)

[0139] Among them, Σ0 is the noise covariance matrix of the prior factor of the global first node, which represents the uncertainty of the initialization information. This value can be given by the measurement standard or empirical value of high-precision measurement equipment.

[0140] Then we can get the error function of the prior factor of the global first node

[0141]

[0142] The factor graph model is calibrated based on the error function of the prior factor of the global first node.

[0143] S1.3. Determine and initialize a factor graph model with deception mitigation capabilities.

[0144] Figure 2 A schematic diagram of a possible implementation of a factor graph optimization train positioning method for countering satellite navigation spoofing attacks in this embodiment is shown.

[0145] In traditional combined positioning systems, IMU pre-integration is used to improve computational efficiency. Before performing nonlinear optimization, IMU observations between adjacent nodes are pre-integrated to generate relative position, velocity, and attitude increments that are independent of the starting position of the integration. During the optimization process, the IMU pre-integration factor combines multiple observations between adjacent nodes and provides a single relative motion constraint. This embodiment further combines IMU measurements with odometry measurements of the train's wheel travel direction provided by the ODO to avoid errors introduced when converting ODO odometry to velocity, while also satisfying the integral form of the nonlinear optimization problem. During the optimization process, IMU / ODO pre-integration serves as the primary mechanism for time-series state recursion. The IMU / ODO pre-integration factor, constructed based on the IMU and ODO measurement models, is used to form probabilistic constraints on the states of adjacent epochs. GNSS pseudorange factors and Doppler velocity factors are used to construct a measurement model based on observations from visible satellites, thereby constraining the relevant state nodes of the current epoch. Furthermore, a map matching factor, constructed with the assistance of the electronic track map (DTM), constrains the train position to the inherent track alignment. A marginalization strategy is used to remove old nodes within the window and to add edge constraints to minimize the loss of historical observation information.

[0146] To achieve multi-source trusted train positioning under spoofing attacks, spoofing countermeasures, including adaptive spoofing detection, identification, and mitigation, are incorporated into the factor graph structure of the tightly coupled GNSS / INS / ODO / DTM positioning. The adaptive spoofing detection and identification mechanism, implemented prior to joint optimization, utilizes IMU / ODO pre-integration outputs, map matching, and GMM error weights to assist in constructing reliable biases and perform binary hypothesis testing to determine whether a spoofing attack has occurred. After a spoofing attack is detected, the spoofed satellites are identified among the visible satellites by normalizing the reliable biases. The adaptive spoofing mitigation strategy activates the spoofing envelope nodes of the identified spoofed satellites, directly optimizing and compensating for the spoofed measurements during the optimization process. The optimal estimated position is used to calculate the GNSS residuals, which are then sampled using a GMM model to fit the residual distribution, supporting adaptive spoofing detection and identification, as well as the construction of GNSS adaptive factors.

[0147] Step S2: When receiving a speed and distance measuring sensor information frame, or when receiving a satellite signal information frame, aligning the timestamps of the speed and distance measuring sensor information frame using linear interpolation, calculate the pre-integration position of the speed and distance measuring sensor.

[0148] A sophisticated IMU pre-integration is implemented to compensate for the Earth's rotation, shape, and gravity changes. The IMU's angular velocity and acceleration It can be modeled as:

[0149]

[0150] Among them, ng and n a are the noise of the gyroscope and accelerometer, respectively, is the direction cosine matrix for transforming from system w to system b, For its inverse transform, It is the projection of the Earth's rotation angular velocity in the Earth-centered, Earth-fixed coordinate system (e system) relative to the inertial coordinate system (i system) in the w system. w is the angular velocity caused by the carrier motion and the curvature of the earth, and are the Coriolis acceleration and centripetal acceleration caused by the earth's rotation and the carrier's motion. is the navigation system (n system) earth gravity coordinate transformation matrix The projection of the earth's gravity in the w system is obtained. The gravity at different locations is usually related to the dimension in which it is located. Related to the height h, set the gravity model to

[0151]

[0152] According to the pulse count of ODO The one-dimensional mileage increment of the train running along the track between two consecutive time stamps can be calculated

[0153]

[0154] Among them, R w is the radius of the driving wheel. Ignoring the lateral and vertical movement of the contact point between the train wheel and the ground, and taking into account the effects of wheel diameter loss and slippage, the distance increment measured by ODO needs to be corrected using a proportional factor. The ODO measurement can be expressed as the following vector form

[0155]

[0156] The pre-integration model of IMU / ODO at two adjacent optimization moments is:

[0157]

[0158] in, is the increment of ODO in system b, is the velocity in frame b, ι a , ι g and ι odo is the Gaussian white noise of the IMU accelerometer, IMU gyroscope and ODO scale factor, Ω is a synthetic matrix whose first column is the quaternion The first line is The transpose of the conjugate quaternion of The antisymmetric matrix formed.

[0159] Step S3: Determine the nearest orbital slice in the orbital information, perform map matching on the pre-integrated position, calculate reliable deviations, construct adaptive detection and recognition statistics, sense deception events, and identify deceptive satellites among visible satellites.

[0160] Before optimizing state estimation, an adaptive spoofing detection and identification method is implemented to confirm spoofing events and identify the specific satellite measurements affected by the spoofing attack. First, the orbit segment closest to the predicted position derived from the IMU / ODO pre-integration at the current epoch is determined based on the range association model. The predicted position is projected onto the orbital slice using the perpendicular projection principle to determine the on-orbit position. Then, the reference pseudorange and reliable deviation are calculated based on the on-orbit position and the spatial coordinates of the visible satellites. The detection statistic and threshold are then calculated, and a binary hypothesis test is performed to detect whether a spoofing event has occurred. After confirming the occurrence of a spoofing event, the normalized reliable deviation is used to identify the satellites affected by the spoofing among the visible satellites.

[0161] S3.1 performs map matching to determine the on-orbit position of the pre-integrated predicted position of the speed and distance measuring sensor.

[0162] A distance association model is constructed to extract the track segment currently occupied by the train. For ease of representation, the track electronic map is simplified to a series of points of interest under the e system. The adjacent points of interest are the coordinates of the endpoints of the track segment. The distance association model predicts the position by judging the IMU / ODO pre-integration The track section occupied by the train is determined by its proximity to the track section.

[0163]

[0164] Where id represents the index of the track segment.

[0165] After confirming the track segment index occupied by the train, the vertical projection model is used to predict the position Projected onto the track line to correct vertical track errors. Projected position It can be calculated as

[0166]

[0167] S3.2 Calculate reliable biases, sense spoofing events, and identify spoofed satellites among visible satellites.

[0168] The equivalent observation quantity of the above speed and distance measuring sensor is the reference pseudorange referred to by formula (12). For the i-th visible satellite, its reference pseudorange can be calculated as:

[0169]

[0170] Reference pseudorange and satellite-corrected observed pseudorange The difference between the two is defined as the reliability deviation and can be calculated as

[0171]

[0172] For real satellites, the reliable biases only contain the unmodeled pseudorange errors ε k,i , and when a deception event occurs, the deception envelope s will be introduced into the deception satellite bias k,i , based on the reliable deviation, the weighted sum of squared residual errors (WSSE) can be constructed as a deception detection statistic.

[0173]

[0174] Where l is the number of visible satellites in the current epoch, W k is an adaptive weighting matrix whose diagonal elements are the inverse of the variance of the satellite residual distribution derived from the Gaussian mixture model. Under non-spoofing conditions, the WSSE detection index follows a chi-square distribution with 1 degree of freedom. Under a spoofing attack, the WSSE detection index no longer follows a chi-square distribution. The detection threshold T is calculated based on the pre-set false alarm rate. h , then the binary hypothesis test for adaptive deception detection is

[0175]

[0176] Under a deception attack, a deception envelope is introduced into the deviation vector so that the detection statistic has a distribution that is different from the normal state and is more likely to exceed the detection threshold, thereby quickly responding to deception events. For a more covert deception attack, in order to avoid alerting the target receiver during the intrusion phase, a pseudorange offset that is small enough or even 0 is often set at the beginning of the intrusion, resulting in a weak deception envelope of the deceptive satellite, which is difficult to support a rapid response and may even cause the deception detection to fail. A successful deception attack will inevitably cause changes in the amplitude of the signal received by the target receiver in the early stage, which is directly reflected in the fluctuation of the carrier noise power spectrum density C / N0 of the satellite signal. Therefore, in order to enhance the detection capability of adaptive deception detection for weak deception attacks, C / N0 is further introduced into the weighting matrix of the detection statistic to compensate for the detection defects caused by the weak deception envelope. The weight matrix for C / N0 enhancement can be constructed as follows:

[0177]

[0178] Among them, Υ(ω1≤Υ≤1) is a quadratic function derived from the visible satellite C / N0, about the center Symmetric; ω1 and ω2 are adjustment factors used to control the minimum value and gradient of the function.

[0179] After a spoofing event is detected, adaptive spoofing identification is further performed to locate the spoofed satellites among the visible satellites. The normalized reliable deviation can be calculated as

[0180]

[0181] In the visible satellite set, the satellite with the largest normalized reliable deviation is determined to be a deceptive satellite, and its deceptive flag is set to Tag = 1. k and the adaptive weighting matrix W k The component corresponding to the deception signal in , to update α k and W k , and recalculate the WSSE detection index, and perform deception detection again until the updated WSSE index is lower than the detection threshold. The algorithm terminates and it is considered that all deceptive satellites have been identified. The remaining signals are judged to be real signals, and their deception flags are set to Tag = 0.

[0182] Figure 3 shows a possible WSSE metric statistic for adaptive deception detection. Figure 4 A possible standardized reliable deviation statistic is shown for adaptive deception identification.

[0183] Step S4: Activate the deception envelope state corresponding to the identified deception satellite, construct the adaptive pseudorange factors of the real / deception satellite, the pre-integration factors of the speed and distance measurement sensors, the map matching factors and the prior factors, jointly optimize the factor graph, and preliminarily solve the train positioning solution.

[0184] S4.1 constructs the GNSS adaptive factor and adds it to the factor graph model.

[0185] The GNSS factor is used to constrain the position, velocity, clock error, and clock drift of the system state node at the current epoch. In a non-deception scenario, the pseudorange of the i-th visible satellite is and Doppler shift d k,i Can be modeled as

[0186]

[0187] Among them, the superscript A represents the real satellite information. and is the position and velocity of the receiver in the e frame, and is the spatial position, velocity, clock error and clock drift of the visible satellite, λ is the wavelength, is the line-of-sight vector between the receiver and the satellite, and is the equivalent pseudorange error caused by the ionosphere, troposphere and Earth rotation, and is the corresponding equivalent pseudorange rate error, and are the unmodeled pseudorange and pseudorange rate errors.

[0188] Under deception attack, deception will be introduced into the pseudo-range observation equation, which can be modeled as

[0189]

[0190] The superscript S represents the information of the spoofed satellite. Here, we only consider the pseudorange spoofing attack. A similar modeling process can be easily extended to the pseudorange rate measurement.

[0191] Then the error function of GNSS pseudorange and Doppler velocity factor can be calculated as:

[0192]

[0193] in, and is the GNSS measurement corrected by equations (18), (19) and (20), and is the covariance matrix corresponding to the true pseudorange, pseudorange rate and deceptive pseudorange, which is generally modeled as a zero-mean Gaussian distribution.

[0194] The covariance matrix of the GNSS factor characterizes the uncertainty of the measurement and is crucial to the state estimation accuracy of the factor graph. A measurement error distribution that does not match the model will increase the difficulty of finding the optimal solution for the factor graph. There are restricted environments such as stations, road cuts, and urban canyons along the railway, which may cause the GNSS measurement error distribution to exhibit non-Gaussian distribution phenomena such as multi-peaks and long tails. And under the influence of deception attacks, the non-Gaussian distribution characteristics of the measurement error will be further enhanced. The single Gaussian distribution assumption will degrade the state estimation performance and lead to a decrease in positioning accuracy. In order to solve the non-Gaussian and time-varying error distribution characteristics, a GNSS adaptive factor with adaptive error modeling is further introduced for factor graph optimization. The GNSS adaptive factor can be calculated as:

[0195]

[0196] in, and Adaptive variance derived for Gaussian mixture models GMM.

[0197] S4.2 constructs the pre-integration factor of the speed and distance measurement sensor and adds it to the factor graph model.

[0198] Given the update period ΔT of the IMU / ODO pre-integration factor, the error function of the IMU / ODO pre-integration factor can be calculated through the state transition of the error term:

[0199]

[0200] in, IMU / ODO measurement IMU / ODO The covariance matrix of represents the uncertainty associated with the IMU / ODO pre-integrated equivalent measurement. and IMU / ODO pre-integration results corrected by the first-order approximation of Equation (9).

[0201] S4.3 Construct a map matching factor and add it to the factor graph model.

[0202] In step S3.1, the track index and on-orbit position closest to the current state have been determined. The vertical track error is the distance between the positioning position and the on-orbit position, which describes the degree to which the positioning position deviates from the track. A smaller vertical track error is more likely to indicate a more accurate positioning solution. Therefore, a map matching factor is constructed based on the vertical track error and added to the optimization process of the factor graph to directly constrain the position. The map matching factor can be calculated as

[0203]

[0204] in, is the variance associated with the map matching factor. It is usually assumed that the vertical track error follows a normal distribution.

[0205] S4.4 constructs the prior factors within the sliding window and adds them to the factor graph model.

[0206] When the number of states in the sliding window exceeds the window capacity, the oldest state will be marginalized, and the sensor measurements corresponding to the marginalized states will be converted into prior factors. The cost function is

[0207]

[0208] in, is the fixed linearization point of the first system state node in the sliding window, J prior and σ prior is the Jacobian matrix of the prior factors and the residuals.

[0209] S4.5 joint optimization to preliminarily solve the train's operating status.

[0210] Finally, by combining the GNSS true pseudorange adaptive factor, GNSS spoofed pseudorange adaptive factor, GNSS Doppler velocity factor, IMU / ODO pre-integration factor, map matching factor and prior factor, the preliminary estimate of the train running state is obtained. It can be obtained by minimizing the global error function

[0211]

[0212] The factor graph solution method used is the Levenberg-Marquardt method in Ceres Solver.

[0213] Step S5: Calculate the pseudorange residuals based on the preliminary positioning solution, compensate for the residuals of the deceptive satellites, update the residual pool, fit the residual distribution of each visible satellite in the residual pool, calculate the adaptive weights, and feed them back into the adaptive deception detection, identification, and mitigation steps to continue iterating to obtain the optimal positioning solution for the train.

[0214] Figure 5 A possible optimal state estimation result is provided for adaptive spoofing mitigation. The Gaussian mixture model (GMM) is an extension of the single Gaussian model. It smoothly fits a density distribution of arbitrary shape by combining multiple Gaussian components with different weights. To balance fitting accuracy and computational timeliness, a Gaussian mixture model with up to three Gaussian components is used to model the core and possible tail of the GNSS residual distribution. Based on the optimal state solved by the factor graph joint optimization, the pseudorange residual of the i-th satellite can be calculated as:

[0215]

[0216] in, and The optimal position, receiver clock error and spoofing envelope are solved. For the identified spoofed satellites, their residuals will be corrected using the estimated spoofing envelope.

[0217] The most recently calculated GNSS residuals are added to the residual pool. When the residual pool exceeds the capacity, the oldest residuals are removed. Subsequently, a GMM is used to fit the residual set for each visible satellite. Because the residuals of spoofed satellites are close to zero after being corrected by the spoofing envelope, their distribution differs significantly from that of real satellites. Using GMM fitting would over-amplify the measurement weights of spoofed satellites in the positioning solution, significantly weakening the contribution of real satellites to the positioning solution. To avoid this, the residuals of identified spoofed satellites are not added to the residual pool.

[0218] The probability density function (PDF) of GNSS residuals can be modeled as

[0219]

[0220] in, are the parameters of the Gaussian mixture model, μ and Σ are the weight, mean, and variance of the Gaussian components.

[0221] If variational inference is used to solve the parameters of the Gaussian mixture model, the optimal GMM parameters are obtained by iteratively executing the variational E step (Expectation) and variational M step (Maximization) until the likelihood function converges or reaches a predefined maximum iteration threshold. Then the adaptive variance of satellite pseudorange observation is calculated as

[0222]

[0223] The covariance matrix formed by it will be fed back into the adaptive deception detection and identification, as well as the construction of the GNSS pseudorange adaptive factor, and the global error function (29) of the factor graph will be continuously iterated to obtain the optimal solution for the train operation status.

[0224] In summary, in the embodiment (1) of the present invention, existing deception mitigation methods generally use least squares optimization and filter estimators as state solvers, and anti-deception technology for fusion positioning with factor graph optimization needs further research. The present invention implements a factor graph optimization method to combat satellite navigation deception attacks, integrating deception detection, identification, and mitigation mechanisms into a multi-source tightly coupled factor graph optimization architecture. During the system state optimization solution process, the method not only estimates and compensates the deception envelope of the deception signal in real time, but also accurately restores the actual operating state of the train.

[0225] (2) At present, the research on anti-spoofing methods for satellite positioning in specific railway scenarios is still in the initial exploratory stage. The present invention specifically studies the problem of satellite navigation deception threats in the train operation environment, makes full use of the train speed and distance sensors and prior track information, and constructs auxiliary constraints through optimization and integration strategies to achieve rapid recovery of state estimation under deception attacks. Without significantly increasing the system hardware cost, the anti-spoofing positioning capability of the train is effectively improved. The anti-spoofing method provided by the present invention is highly consistent with the railway application scenario, and can be widely used in various applications driven by satellite navigation-based train positioning, with good engineering feasibility and promotion prospects.

[0226] Those skilled in the art will appreciate that the accompanying drawings are merely schematic diagrams of an embodiment, and the modules or processes in the accompanying drawings are not necessarily required to implement the present invention.

[0227] From the above description of the embodiments, it can be seen that those skilled in the art can clearly understand that the present invention can be implemented by means of software plus the necessary general-purpose hardware platform. Based on this understanding, the technical solution of the present invention, or the portion that contributes to the prior art, can be embodied in the form of a software product. This computer software product can be stored in a storage medium such as ROM / RAM, a magnetic disk, or an optical disk, and includes a number of instructions for enabling a computer device (which can be a personal computer, a server, or a network device, etc.) to execute the methods described in various embodiments of the present invention or certain parts of the embodiments.

[0228] Each embodiment in this specification is described in a progressive manner. The same or similar parts between the embodiments can be referred to each other. Each embodiment focuses on the differences from other embodiments. In particular, for the device or system embodiments, since they are basically similar to the method embodiments, the description is relatively simple. For the relevant parts, refer to the partial description of the method embodiments. The device and system embodiments described above are merely schematic, wherein the units described as separate components may or may not be physically separated, and the components displayed as units may or may not be physical units, that is, they may be located in one place, or they may be distributed on multiple network units. Some or all of the modules can be selected according to actual needs to achieve the purpose of the scheme of this embodiment. A person of ordinary skill in the art can understand and implement it without making any creative efforts.

[0229] The above description is merely a preferred embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any changes or substitutions that can be easily conceived by a person skilled in the art within the technical scope disclosed in the present invention should be included in the scope of protection of the present invention. Therefore, the scope of protection of the present invention should be based on the scope of protection of the claims.

Claims

1. A factor graph optimization train positioning method to counter satellite navigation spoofing attacks, characterized by: include: determining and initializing a system state node having a deception envelope component, a factor graph model having deception mitigation capabilities, and related parameters; When receiving a speed and distance measuring sensor information frame, or when receiving a satellite signal information frame, after aligning the timestamps of the speed and distance measuring sensor information frame using linear interpolation, a pre-integrated position of the speed and distance measuring sensor is calculated; Determine the nearest orbital slice from the orbital information, map-match the pre-integrated position, calculate reliable deviations, construct adaptive detection and recognition statistics, sense spoofing events, and identify spoofed satellites among visible satellites; Activate the deception envelope state corresponding to the identified spoofed satellite, construct the adaptive pseudorange factors of the real / spoofed satellite, the pre-integration factors of the speed and range sensors, the map matching factors and the prior factors, jointly optimize the factor graph, and preliminarily solve the train positioning solution; Based on the preliminary positioning solution, the pseudorange residuals are calculated and the residuals of the deceptive satellites are compensated. The residual pool is updated and the residual distribution of each visible satellite in the residual pool is fitted. The adaptive weights are calculated and fed back into the adaptive deceptive detection, identification, and mitigation steps to continue iterating to obtain the optimal positioning solution for the train.

2. The method according to claim 1, characterized in that The determining and initializing of a system state node having a deception envelope component, a factor graph model having deception mitigation capabilities, and related parameters includes: S1.

1. Determine and initialize the system state node to be estimated with the deceptive envelope component; In the multi-source tightly coupled positioning, the raw observed pseudorange and Doppler shift of the global navigation satellite system (GNSS), the raw observed values of the three-axis gyroscope and accelerometer of the inertial navigation system (INS), and the raw observed wheel axle pulse count of the wheel axle speed sensor (ODO) are processed and sent to the factor graph optimization for joint optimization positioning. If there is a deceptive satellite among the visible satellites, the deceptive envelope state component of the deceptive satellite will be activated. The world coordinate system (W system) with the initial position as the origin is selected, which is also the reference coordinate system of the factor graph. Then, within a time window, the system's multidimensional state node set χ and state node x are k Build as: Among them, x k is the 18-dimensional state node at the k-th epoch, and is the three-dimensional position and velocity of the train, is the train attitude expressed in quaternion form, and The three-dimensional accelerometer and gyroscope zero bias of the IMU system, and is the clock error and clock drift of the GNSS receiver, is the scale factor of ODO, s k is the deception envelope, m is the number of identified deception satellites, and n is the sliding window size; S1.

2. Determine and initialize the prior factor of the global first node of the deception mitigation factor graph model The initial state prior information x0 is given by the measurement value or experience value of the measuring device and corresponds one-to-one with the system state node, which can be expressed as: in, To initialize position, velocity and attitude, and To initialize the IMU accelerometer and gyroscope bias, and To initialize the GNSS receiver clock error and drift, To initialize the scale factor of ODO; Assuming that the error w0 of the initialization information obeys the Gaussian distribution N, the measurement equation of the prior factor of the global first node is expressed as: x=x0+w0,w0~N(0,Σ0) (3) Among them, Σ0 is the noise covariance matrix of the prior factor of the first global node, which represents the uncertainty of the initialization information; Get the error function of the prior factor of the global first node Calibrate the factor graph model based on the error function of the prior factor of the global first node; S1.

3. Determine and initialize a factor graph model with deception mitigation capabilities; The inertial measurement unit (IMU) measurements and the mileage measurements of the train wheels provided by the ODO (Operational Detection and Tracking) are combined. During the optimization process, IMU / ODO pre-integration is used as the main body for time series state recursion. The IMU / ODO pre-integration factor, constructed based on the IMU and ODO measurement models, is used to form probabilistic constraints on the states of adjacent epochs. The GNSS pseudorange factor and Doppler velocity factor, constructed based on visible satellite observations, constrain the relevant state nodes of the current epoch. The map matching factor, constructed with the assistance of the electronic track map (DTM), constrains the train position to the inherent track line. Deception countermeasures are incorporated into the factor graph structure of tightly coupled GNSS / INS / ODO / DTM positioning, including adaptive deception detection, identification, and mitigation. IMU / ODO pre-integration outputs, map matching, and GMM error weights are used to assist in constructing reliable biases. Binary hypothesis testing is performed to determine whether a deception attack has occurred. After a deception attack is detected, the deceptive satellites are identified among the visible satellites by normalizing the reliable biases. An adaptive deception elimination strategy activates the deception envelope nodes of the identified deceptive satellites, directly optimizing and compensating for the deceptive measurements during the optimization process, and determining and initializing a factor graph model with deception mitigation capabilities.

3. The method according to claim 2, characterized in that The method of calculating the pre-integrated position of the speed and ranging sensor when receiving the speed and ranging sensor information frame or aligning the timestamp of the speed and ranging sensor information frame using linear interpolation when receiving the satellite signal information frame includes: Set up speed and distance measurement sensors including inertial sensor IMU and axle speed and distance measurement sensor ODO; The angular velocity of the IMU and acceleration Modeled as: Among them, n g and n a are the noise of the gyroscope and accelerometer, respectively, is the direction cosine matrix of the transformation from the w system to the IMU carrier system (b system), For its inverse transform, is the projection of the Earth's rotation angular velocity in the Earth-centered Earth-fixed coordinate system (e system) relative to the inertial coordinate system (i system) in the w system, w is the angular velocity caused by the carrier motion and the curvature of the earth, and is the Coriolis acceleration and centripetal acceleration caused by the earth's rotation and the carrier's motion, is the navigation system (n system) earth gravity coordinate transformation matrix The projection of the earth's gravity in the w system is obtained. The gravity at different locations is usually related to the dimension in which it is located. Related to the height h, set the gravity model to: According to the pulse count of ODO Calculate the one-dimensional mileage increment of the train along the track between two consecutive time stamps Among them, R w is the radius of the driving wheel; The ODO measurement is expressed in the following vector form: The pre-integration model of IMU / ODO at two adjacent optimization moments is: in, is the increment of ODO in system b, is the velocity in frame b, ι a , ι g and ι odo is the Gaussian white noise of the IMU accelerometer, IMU gyroscope and ODO scale factor, Ω is a synthetic matrix whose first column is the quaternion The first line is The transpose of the conjugate quaternion of The antisymmetric matrix formed.

4. The method according to claim 3, characterized in that The steps of determining the nearest orbital slice in the orbital information, performing map matching on the pre-integrated position, calculating reliable deviations, constructing adaptive detection and recognition statistics, sensing deception events, and identifying deceptive satellites among visible satellites include: S3.1 performs map matching to determine the on-orbit position of the pre-integrated predicted position of the speed and distance measuring sensor; Extract the track segment currently occupied by the train and simplify the track electronic map into a series of points of interest under the e-system The adjacent points of interest are the endpoint coordinates of the track segment, and the positions are predicted by judging the IMU / ODO pre-integration The track segment occupied by the train is determined by its proximity to the track segment; Where id represents the index of the track segment; After confirming the track segment index occupied by the train, the train predicted position Projected onto the track line, train projection position Calculated as: S3.2 Calculate reliable biases to sense spoofing events and identify spoofed satellites among visible satellites; The reference pseudorange of the i-th visible satellite is calculated as: Reference pseudorange and satellite-corrected observed pseudorange The difference is defined as the reliability deviation and is calculated as: For real satellites, the reliable biases only contain the unmodeled pseudorange errors ε k,i , and when a deception event occurs, the deception envelope s will be introduced into the deception satellite bias k,i , based on the reliable deviation, a weighted residual sum of squares is constructed as the deception detection statistic; Where l is the number of visible satellites in the current epoch, W k is an adaptive weighting matrix whose diagonal elements are the inverse of the satellite residual distribution variance derived by the Gaussian mixture model. Under no deception conditions, the WSSE detection index obeys the chi-square distribution with a degree of freedom of l. However, under deception attacks, the WSSE detection index will no longer obey the chi-square distribution. The detection threshold T is calculated based on the pre-set false alarm rate. h , then the binary hypothesis test for adaptive deception detection is: The carrier noise power spectral density C / N0 of the satellite signal is introduced into the weighting matrix of the detection statistic, and the C / N0 enhanced weight matrix is constructed as follows: Among them, Υ(ω1≤Υ≤1) is a quadratic function derived from the visible satellite C / N0, about the center Symmetric; ω1 and ω2 are adjustment factors used to control the minimum value and gradient of the function; After a spoofing event is detected, adaptive spoofing identification is further performed to locate the spoofed satellite among the visible satellites. The normalized reliable deviation is calculated as: In the visible satellite set, the satellite with the largest normalized reliable deviation is judged as a deceptive satellite, and its deceptive flag is set to Tag = 1. k and the adaptive weighting matrix W k The component corresponding to the deception signal in , to update α k and W k , and recalculate the WSSE detection index, and perform deception detection again until the updated WSSE index is lower than the detection threshold. The algorithm terminates and it is considered that all deceptive satellites have been identified. The remaining signals are judged to be real signals, and their deception flags are set to Tag = 0.

5. The method according to claim 4, characterized in that The deception envelope state corresponding to the activated and identified deception satellites is constructed, and adaptive pseudorange factors of real / deception satellites, pre-integration factors of speed and distance sensors, map matching factors, and prior factors are constructed. The factor graph is jointly optimized to preliminarily solve the train positioning solution, including: S4.1 constructs a GNSS adaptive factor and adds it to the factor graph model; The GNSS factor is used to constrain the position, velocity, clock error and clock drift of the system state node in the current epoch. In the non-deception scenario, the pseudorange of the i-th visible satellite is and Doppler shift d k,i Modeled as: Among them, the superscript A represents the real satellite information. and is the position and velocity of the receiver in the e frame, and is the spatial position, velocity, clock error and clock drift of the visible satellite, λ is the wavelength, is the line-of-sight vector between the receiver and the satellite, and is the equivalent pseudorange error caused by the ionosphere, troposphere and Earth rotation, and is the corresponding equivalent pseudorange rate error, and are the unmodeled pseudorange and pseudorange rate errors; Under the deception attack, the pseudo-range observation equation including the deception is modeled as: Among them, the superscript S represents the information of deceiving satellites; The error function of GNSS pseudorange and Doppler velocity factor is calculated as: in, and is the GNSS measurement corrected by equations (18), (19) and (20), and is the covariance matrix corresponding to the true pseudorange, pseudorange rate and deceptive pseudorange, which is modeled as a zero-mean Gaussian distribution; The GNSS adaptive factor with adaptive error modeling is introduced for factor graph optimization. The GNSS adaptive factor is calculated as: in, and Adaptive variance derived for Gaussian mixture model GMM; S4.2 Given the update period ΔT of the IMU / ODO pre-integration factor, calculate the error function of the IMU / ODO pre-integration factor through the state transition of the error term in, IMU / ODO measurement IMU / ODO The covariance matrix of represents the uncertainty associated with the IMU / ODO pre-integrated equivalent measurement, and The IMU / ODO pre-integration result corrected by the first-order approximation of formula (9); S4.3 Map matching factor is calculated as: in, is the variance associated with the map matching factor; S4.4 constructs a priori factors within the sliding window and adds the prior factors to the factor graph model; When the number of states in the sliding window exceeds the window capacity, the oldest state will be marginalized, and the sensor measurements corresponding to the marginalized states will be converted into prior factors. The cost function is: in, is the fixed linearization point of the first system state node in the sliding window, J prior and σ prior is the Jacobian matrix of the prior factors and the residual; S4.5 Preliminary estimation of train running status by combining GNSS true pseudorange adaptive factor, GNSS spoofed pseudorange adaptive factor, GNSS Doppler velocity factor, IMU / ODO pre-integration factor, map matching factor and prior factor By minimizing the global error function: The factor graph solution method used is the Levenberg-Marquardt method in Ceres Solver.

6. The method according to claim 5, characterized in that The aforementioned steps of calculating pseudorange residuals based on the preliminary positioning solution, compensating for the residuals of the deceptive satellites, updating the residual pool, fitting the residual distribution of each visible satellite in the residual pool, calculating the adaptive weights and feeding them back into the adaptive deceptive detection, identification, and mitigation steps, and continuing to iterate to obtain the optimal positioning solution for the train include: A Gaussian mixture model with up to three Gaussian components is used to model the core and possible tails of the GNSS residual distribution. According to the optimal state of the factor graph joint optimization solution, the pseudorange residual of the i-th satellite is calculated as: in, and To solve for the optimal position, receiver clock error and spoofing envelope, for the identified spoofed satellite, its residual will be corrected using the estimated spoofing envelope; The most recently calculated GNSS residuals are added to the residual pool. When the residual items in the residual pool exceed the capacity, the earliest added residual items will be removed. GMM is used to fit the residual set of each visible satellite. Since the residuals of spoofed satellites are close to 0 after being corrected by the spoofing envelope, the residuals of identified spoofed satellites will not be added to the residual pool. The probability density function of GNSS residuals is modeled as: in, are the parameters of the Gaussian mixture model, μ and Σ are the weight, mean, and variance of the Gaussian components; If variational inference is used to solve the parameters of the Gaussian mixture model, the optimal GMM parameters are obtained by iteratively performing the variational E step and variational M step until the likelihood function converges or reaches a predefined maximum iteration threshold. Then the adaptive variance of satellite pseudorange observation is calculated as: In the iterative positioning solution, the adaptive variance will be fed back to the adaptive deception detection, identification and mitigation, and the optimal estimate of the train operation status will be obtained by continuous iteration (29).

Images