Factor graph optimization train positioning method against spoofing attacks of satellite navigation

Abstract

The application provides a factor graph optimization train positioning method against satellite navigation spoofing attack. The method comprises the following steps: determining a factor graph model with spoofing mitigation capability and related parameters; when a speed and distance measurement sensor information frame is received, calculating a speed and distance measurement sensor pre-integration position; determining a nearest track piece in track information, perceiving a spoofing event and identifying a spoofing satellite in a visible satellite; constructing an adaptive pseudo-range factor of a real / spoofing satellite, a speed and distance measurement sensor pre-integration factor, jointly optimizing a factor graph, and solving a positioning solution of the train; calculating a pseudo-range residual according to the positioning solution, compensating for the residual of the spoofing satellite, updating a residual pool, fitting a residual distribution of each visible satellite in the residual pool, calculating an adaptive weight, and continuously iterating to obtain an optimal positioning solution of the train. The application realizes spoofing attack detection, identification and mitigation based on a factor graph optimization framework, recovers a real running state of the train before the attack, and enhances spoofing defense capability.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of rail transit train positioning technology, and in particular to a factor graph optimization train positioning method against satellite navigation spoofing attacks. BACKGROUND

[0002] Global Navigation Satellite System (GNSS) enabled train positioning technology has accelerated the transformation of the traditional train control system (hereinafter referred to as "train control system") to a new generation of train control system with higher running speed, lower construction and maintenance cost and more intelligence. With the introduction of GNSS, its vulnerability exposes train positioning to unintentional interference caused by environmental shielding along the railway and radio frequency interference caused by malicious attackers. Especially in radio frequency interference, GNSS spoofing attacks with active induction properties can manipulate the Position, Velocity and Time (PVT) solution of the train-mounted receiver without arousing the vigilance of the victim train, seriously threatening the safety of train operation. A single GNSS system has been difficult to meet the requirements of the train control system for high-availability, high-precision and high-reliability positioning technology in a complex interference environment. In order to improve the continuity, accuracy and robustness of GNSS-based train positioning, the use of multi-source fusion positioning algorithm to fuse the measurement information of multiple independent GNSS sensors is expected to achieve multi-source trusted positioning of the train. Joint spoofing immune multi-source information can enhance the ability of train positioning to recover the true PVT solution under spoofing attacks, and further achieve multi-source trusted positioning of the train in a spoofing attack scenario. Existing anti-spoofing technologies mainly focus on spoofing detection technology, which aims to realize timely perception of spoofing attacks and trigger the warning mechanism. Spoofing mitigation technology, as a key link to further respond to the influence of spoofing attacks on positioning solution after the perception of spoofing attacks, has started relatively late. In anti-spoofing technology, state estimation is used as a back-end optimization to estimate and recover the spoofed system state. Most existing schemes still mainly use the least square estimation and filtering fusion framework. The emerging factor graph optimization architecture, as a global optimization solver, has gradually attracted the attention of researchers and has been proven to be superior to the positioning accuracy and robustness of filtering fusion algorithms. Therefore, exploring anti-spoofing technology based on the factor graph optimization architecture to integrate spoofing detection, mitigation and state estimation is not only expected to achieve a more robust response strategy to spoofing attacks, but also may achieve a state recovery capability superior to filtering methods in attack scenarios.

[0003] Under malicious spoofing attacks, GNSS-based train positioning performance may degrade or even fail, making it difficult to meet the demand for trusted positioning. To deal with the threat of GNSS spoofing attacks, it is urgent to combine existing train on-board positioning sensors, advanced anti-spoofing technologies, and positioning solution technologies to detect and mitigate spoofing attacks, trigger safety warnings in a timely manner, and accurately restore the true state of the train, achieving multi-source trusted positioning of the train under spoofing attacks. The present invention aims to mitigate the impact of spoofing attacks on train satellite positioning, fully utilize the spoofing-immune multi-source auxiliary information provided by the existing configuration and commonly used sensor equipment of the train, integrate spoofing detection, identification, and mitigation techniques into a multi-source fusion positioning architecture based on factor graph optimization, and thus simultaneously achieve spoofing attack mitigation and train true state restoration.

[0004] Currently, anti-spoofing technologies based on measurement domain information processing in the prior art generally use two measures, "identification-exclusion" and "detection-utilization", to suppress the impact of spoofing attacks on PVT solutions. The "identification-exclusion" countermeasure is to directly exclude the identified spoofed measurements from the positioning solution after confirming that a specific satellite is affected by spoofing attacks, to isolate the impact of spoofing attacks. Fault Detection Exclusion (FDE) and Receiver Autonomous Integrity Monitoring (RAIM) technologies have been proven to be able to isolate spoofing attacks. However, the "identification-exclusion" technique may result in a situation where the remaining real satellite measurements do not meet the minimum number required for positioning solution, thereby reducing the availability of satellite positioning. Related research shows that spoofed satellite information is helpful for restoring true PVT information, therefore, the "detection-utilization" countermeasure is derived, which uses weight online adjustment or error compensation to utilize and mitigate the impact of spoofed measurements on positioning solutions after detecting spoofing attacks. Most existing anti-spoofing methods are based on this idea, such as robust estimation and adaptive Kalman filtering. Spoofing envelope is defined as the deviation caused by spoofing attacks in satellite measurements, and existing research has explored the use of filtering positioning framework to model, estimate, and compensate for spoofing envelope for anti-spoofing capability. Compared to the "identification-exclusion" measure, the "detection-utilization" measure retains the information of spoofed satellites, which not only ensures the availability and continuity of positioning, but also has certain improvement in positioning accuracy. Among them, the weight online adjustment strategy does not need to rely on an accurate error correction model, and obviously becomes the best choice for anti-spoofing technologies based on measurement domain information processing. The error correction model may further restore the true state under spoofing attacks and improve positioning accuracy.

[0005] A train positioning method against spoofing attacks of satellite navigation in the prior art includes: a GNSS / INS tight coupling method based on extended Kalman filtering for resisting spoofing attacks, which uses Kalman filter innovation and estimated spoofing envelope in the iteration process to perform spoofing identification.

[0006] Another train positioning method against spoofing attacks of satellite navigation in the prior art integrates sparse estimation theory into extended Kalman filtering, converts the GNSS state estimation problem into L1 norm, and adopts a least absolute shrinkage and selection algorithm (LASSO) for sparse estimation, so as to estimate the state while mitigating the influence of spoofed pseudoranges on positioning solution.

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

[0008] (1) Most spoofing mitigation techniques in the prior art still use filter estimation as a state solver, although a few scholars have tried to realize the cooperative processing of spoofing detection and state estimation in the filter framework, but overall there are still problems of insufficient robustness and recovery ability. Factor graph optimization, as a state estimation method that has developed rapidly in recent years, provides a new solution for the development of anti-spoofing technology due to its advantages in information modeling flexibility, multi-source observation fusion and nonlinear optimization. Combining factor graph optimization with spoofing mitigation strategies is expected to significantly improve the state recovery ability when facing GNSS spoofing attacks. However, the current research on spoofing mitigation techniques based on factor graph optimization is still limited, especially in constructing a unified optimization model that can realize spoofing mitigation and robust state estimation at the same time, which needs further exploration.

[0009] (2) Most spoofing mitigation methods in the prior art are designed for general satellite positioning scenarios, lack of special design tailored for the characteristics of train control systems, and are difficult to be directly integrated and applied in existing train operation control systems. In the train control system, existing non-GNSS sensors are not easily affected by spoofing attacks and can provide redundant and reliable auxiliary information for positioning. Therefore, it is urgent to develop spoofing interference mitigation and state recovery methods for train satellite positioning scenarios, make full use of the unique multi-sensor and track prior information of trains, and realize high robustness and high reliability of spoofing protection capability.

[0010] (3) Existing anti-deception technologies mostly use filtering estimation as the state solution method. However, the factor graph optimization algorithm, which has emerged in recent years, exhibits many advantages in modeling flexibility, information fusion efficiency and scalability, and is expected to achieve better anti-deception performance than traditional filters. Nevertheless, the anti-deception potential based on the factor graph optimization framework has not been fully tapped. Especially for special application scenarios such as satellite positioning-based train operation control, there is still a lack of special deception interference mitigation methods designed for train positioning characteristics. SUMMARY

[0011] Embodiments of the present application provide a factor graph optimization train positioning method against satellite navigation deception attack to effectively guarantee the credibility of train positioning.

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

[0013] A factor graph optimization train positioning method against satellite navigation deception attack, comprising:

[0014] determining and initializing system state nodes with deception envelope components, factor graph models with deception mitigation capabilities, and related parameters;

[0015] calculating speed and distance sensor pre-integrated positions after linear interpolation timestamp alignment of speed and distance sensor information frames upon receiving the speed and distance sensor information frames, or upon receiving satellite signal information frames;

[0016] determining the nearest track segment in the track information, map matching the pre-integrated positions, calculating reliable biases, constructing adaptive detection and identification statistics, perceiving deception events and identifying deception satellites among visible satellites;

[0017] activating deception envelope states corresponding to the identified deception satellites, constructing adaptive pseudo-range factors of real / deception satellites, speed and distance sensor pre-integrated factors, map matching factors and prior factors, jointly optimizing the factor graph, and preliminarily solving the positioning solution of the train;

[0018] calculating pseudo-range residuals according to the preliminary positioning solution, compensating for the residuals of the deception satellites, updating the residual pool, fitting the residual distribution of each visible satellite in the residual pool, calculating adaptive weights and feeding back to the adaptive deception detection, identification and mitigation steps, and continuing to iterate to obtain the optimal positioning solution of the train.

[0019] Preferably, the determination and initialization of system state nodes with deception envelope components, factor graph models with deception mitigation capabilities and related parameters comprise:

[0020] S1.1, determining and initializing system state nodes with deception envelope components;

[0021] In multi-source tightly coupled positioning, the raw observations of pseudo-range and Doppler shift of GNSS, the raw observations of three-axis gyroscope and accelerometer of INS, and the raw observations of wheel axle pulse count of ODO are processed and sent into factor graph optimization for joint optimization positioning. If there is a spoofing satellite in the visible satellite, the spoofing envelope state component of the spoofing satellite will be activated. 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. In a time window, the multi-dimensional state node set χ and the state node x of the system k are constructed as:

[0022]

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

[0024] S1.2, determine and initialize the prior factor of the global first node of the spoofing mitigation factor graph model

[0025] The initial state prior information x0 is given by the measurement value or empirical value of the measurement device, and is one-to-one corresponding to the system state node, and is expressed as:

[0026]

[0027] wherein, is the initialized position, velocity and attitude, and are the initialized IMU accelerometer and gyroscope bias, and are the initialized GNSS receiver clock error and clock drift, is the initialized scale factor of ODO;

[0028] Assuming that the error w0 of the initialization information is subject to 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] where∑0is the noise covariance matrix of the prior factor of the global first node, indicating the uncertainty of the initialization information;

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

[0032]

[0033] calibrate the factor graph model according to the error function of the prior factor of the global first node;

[0034] S1.3, determine and initialize the factor graph model with spoofing mitigation capability;

[0035] Combine the inertial sensor IMU measurement and the mileage measurement of the train wheel forward direction provided by ODO. In the optimization process, the IMU / ODO pre-integration is used as the main body to perform time series state recursion. The IMU / ODO pre-integration factor constructed based on the IMU and ODO measurement model is used to form a probability constraint on the adjacent epoch state. The GNSS pseudorange factor and the Doppler velocity factor construct a measurement model based on the observation of the visible satellite to constrain the related state node of the current epoch. The map matching factor constructed with the aid of the electronic track map will constrain the train position to the inherent track line;

[0036] Spoofing countermeasures are added to the tightly coupled GNSS / INS / ODO / DTM positioning factor graph structure, including adaptive spoofing detection, identification and mitigation. The IMU / ODO pre-integration output, map matching and GMM error weight are used to assist in constructing reliable bias. Binary hypothesis testing is performed to determine whether there is a spoofing attack. After detecting a spoofing attack, the spoofing satellite in the visible satellite is determined by the standardized reliable bias. The adaptive spoofing elimination strategy activates the spoofing envelope node of the identified spoofing satellite, directly optimizes and solves the spoofing measurement in the optimization process, and determines and initializes the factor graph model with spoofing mitigation capability.

[0037] Preferably, after the linear interpolation timestamp alignment of the speed and distance sensor information frame is performed when the speed and distance sensor information frame is received, or when the satellite signal information frame is received, the speed and distance sensor pre-integration position is calculated, comprising:

[0038] The speed and distance sensor includes an inertial sensor IMU and an axle speed and distance sensor ODO;

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

[0040]

[0041] where n g and n a are the noise of the gyroscope and accelerometer, respectively, is the direction cosine matrix of w-frame to IMU carrier frame (b-frame), is its inverse, is the projection of the earth rotation angular velocity of the earth-fixed frame (e-frame) to the w-frame, is the dependent angular velocity of the w-frame due to the carrier motion and earth curvature, and are the Coriolis acceleration and centripetal acceleration due to the earth rotation and carrier motion, is the coordinate transformation matrix of the navigation frame (n-frame) to the w-frame The projection of the earth gravity in the w-frame is obtained by The earth gravity is usually related to the latitude

[0042]

[0043] According to the pulse count of ODO The one-dimensional mileage increment of the train along the track within two consecutive time stamps is calculated

[0044]

[0045] where R w is the radius of the driving wheel;

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

[0047]

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

[0049]

[0050] where, is the increment of ODO in b-frame, is the velocity in b-frame, ι a , ι g and ι odo are the Gaussian white noises of IMU accelerometer, IMU gyroscope and ODO scale factor, Ω is a synthetic matrix whose first column is the quaternion the first row is the transpose of the conjugate quaternion of and the rest is The antisymmetric matrix is constructed.

[0051] Preferably, the determination of the nearest track segment in the track information, the map matching of the pre-integration position, the calculation of the reliable bias, the construction of the adaptive detection and identification statistics, the sensing of the spoofing event and the identification of the spoofing satellite in the visible satellite, comprise:

[0052] S3.1 Perform map matching to determine the in-orbit position of the pre-integration prediction position of the speed and distance measuring sensor;

[0053] Extract the track segment currently occupied by the train, and simplify the track electronic map to a series of points of interest in the e system The adjacent points of interest are the end point coordinates of the track segment, and the track segment occupied by the train is determined by judging the proximity of the IMU / ODO pre-integration prediction position 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 prediction position is projected onto the track line, and the train projection position is calculated as:

[0057]

[0058] S3.2 Calculate the reliable bias, sense the spoofing event and identify the spoofing satellite in the visible satellite;

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

[0060]

[0061] The difference between the reference pseudo-range and the satellite-corrected observed pseudo-range is defined as the reliable bias, which is calculated as:

[0062]

[0063] For a real satellite, the reliable bias only contains the unmodeled pseudo-range error ε k,i , and when a spoofing event occurs, the spoofing envelope s k,i will be introduced into the bias of the spoofing satellite. Based on the reliable bias, a weighted residual sum of squares is constructed as a spoofing detection statistic;

[0064]

[0065] Where l is the number of visible satellites at the current epoch, W kThe adaptive weighting matrix has its diagonal elements being 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 degrees of freedom. However, under a spoofing attack, the WSSE detection index will no longer follow a chi-square distribution. The detection threshold T is calculated based on the pre-set false alarm rate. h The binary hypothesis test for adaptive deception detection is:

[0066]

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

[0068]

[0069] Where Υ(ω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 detecting a spoofing event, adaptive spoofing identification is further performed to locate the spoofing satellite among the visible satellites. The standardized reliability bias is calculated as follows:

[0071]

[0072] In the visible satellite set, the satellite with the largest standardized reliability deviation is identified as a spoofing satellite, and its spoofing flag is set to Tag=1. This is achieved by removing the satellite with the highest reliability deviation α. k and the adaptive weighting matrix W k The components corresponding to the deception signal are used to update α. k and W k The algorithm then recalculates the WSSE detection index and performs spoofing detection again until the updated WSSE index is lower than the detection threshold. At this point, the algorithm terminates, deeming all spoofing satellites identified and the remaining signals as genuine signals, with their spoofing flag set to Tag=0.

[0073] Preferably, the activation identification of the deception satellite's corresponding deception envelope state, constructing an adaptive pseudorange factor for the real / deception satellite, a pre-integration factor for the velocity and ranging sensor, a map matching factor, and a priori factors, jointly optimizing the factor graph, and initially solving the train's positioning solution, includes:

[0074] S4.1 Construct GNSS adaptive factors and add them to the factor graph model;

[0075] The GNSS factors are used to constrain the position, velocity, clock bias and clock drift of the system state nodes at the current epoch. In the non-spoofing scenario, the pseudorange and Doppler shift d and Doppler shift d k,i are modeled as:

[0076]

[0077] where the superscript A denotes the true satellite information, and are the position and velocity of the receiver in the e-frame, and are the spatial position, velocity, clock bias and clock drift of the visible satellites, and λ is the wavelength, is the line-of-sight vector between the receiver and the satellites, and are the equivalent pseudorange errors caused by ionosphere, troposphere and earth rotation, and are the corresponding equivalent pseudorange rate errors, and are the unmodeled pseudorange and pseudorange rate errors;

[0078] In the spoofing attack, the pseudorange observation equation including the spoofing is modeled as:

[0079]

[0080] where the superscript S denotes the spoofing satellite information;

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

[0082]

[0083] where, and are the GNSS measurements corrected by equations (18), (19) and (20), and are the covariance matrices of the true pseudorange, pseudorange rate and spoofed pseudorange, which are modeled as zero-mean Gaussian distributions;

[0084] The GNSS adaptive factors with adaptive error modeling are introduced for factor graph optimization, which are calculated as:

[0085]

[0086] where, and are the adaptive variances derived from the Gaussian mixture model (GMM).

[0087] S4.2 Given the update period AT of IMU / ODO pre-integration factor, the error function of IMU / ODO pre-integration factor is calculated by the state transition of error term

[0088]

[0089] where, is the covariance matrix of IMU / ODO measurement z IMU / ODO , which represents the uncertainty of IMU / ODO pre-integration equivalent measurement correlation, and the IMU / ODO pre-integration result corrected by the first-order approximation of equation (9);

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

[0091]

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

[0093] S4.4 The prior factor within the sliding window is constructed, and the prior factor is added to the factor graph model;

[0094] When the number of states within the sliding window exceeds the window capacity, the oldest state will be marginalized, and the sensor measurement corresponding to the marginalized state will be converted into a prior factor, and its cost function is:

[0095]

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

[0097] S4.5 The preliminary estimation of train running state is obtained by combining GNSS real pseudorange adaptive factor, GNSS spoofing pseudorange adaptive factor, GNSS Doppler velocity factor, IMU / ODO pre-integration factor, map matching factor and prior factor The global error function is minimized to obtain:

[0098]

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

[0100] Preferably, the pseudo-range residuals are calculated from the preliminary positioning solution, and the residuals of spoofed satellites are compensated, the residual pool is updated, the residual distribution of each visible satellite in the residual pool is fitted, the adaptive weights are calculated and fed back to the adaptive spoofing detection, identification and mitigation steps, and the optimal positioning solution of the train is obtained by continuous iteration, including:

[0101] A Gaussian Mixture Model (GMM) with at most three Gaussian components is used to model the core and possible tail of the GNSS residual distribution, and the optimal state is solved by factor graph joint optimization. The pseudo-range residual of the i-th satellite is calculated as:

[0102]

[0103] wherein, and are the optimal position, receiver clock bias and spoofing envelope solved, and for the identified spoofed satellite, its residual will be corrected using the estimated spoofing envelope;

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

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

[0106]

[0107] wherein, is the parameter of the Gaussian Mixture Model (GMM), μ and Σ are the weight, mean and variance of the Gaussian component;

[0108] If variational inference is used to solve the parameters of the Gaussian Mixture Model (GMM), the variational E-step and the variational M-step are performed by iteration until the likelihood function converges or reaches a predefined maximum iteration threshold, and the optimal GMM parameters are obtained and the adaptive variance of the satellite pseudo-range observation is calculated as:

[0109]

[0110] In the iterative positioning solution, the adaptive variance will be fed back to the adaptive spoofing detection, identification and mitigation, and the optimal estimation of the train running state is obtained by continuous iteration of formula (29).

[0111] From the technical solutions provided by the above-mentioned embodiments of the present application, it can be seen that the present application constructs a spoofing interference mitigation mechanism suitable for the actual situation and characteristics of train satellite positioning, so as to guarantee the credibility of train positioning and the safety of train control decision-making. The present application aims to mitigate the threat of spoofing interference of train satellite positioning, and constructs a set of factor graph optimization framework to realize spoofing attack detection, identification and mitigation, while recovering the real running state of the train before the attack, and enhancing the spoofing defense capability.

[0112] Additional aspects and advantages of the present application will be described in the following description and will be apparent from the following description and the organization of the present application. BRIEF DESCRIPTION OF DRAWINGS

[0113] In order to more clearly illustrate the technical solutions of the embodiments of the present application, the drawings needed in the embodiment description will be briefly introduced as follows. Obviously, the drawings in the following description are only some embodiments of the present application, and other drawings can be obtained by those skilled in the art without creative labor on the basis of these drawings.

[0114] Figure 1 A factor graph optimization train positioning method for resisting satellite navigation spoofing attack provided by the embodiment of the present application is shown in the flowchart.

[0115] Figure 2 An implementation schematic diagram of a possible train multi-source fusion factor graph optimization model for resisting GNSS spoofing attack provided by the embodiment of the present application is shown in the flowchart.

[0116] Figure 3 A possible adaptive spoofing detection result schematic diagram provided by the embodiment of the present application is shown in the flowchart.

[0117] Figure 4 A possible adaptive spoofing identification result schematic diagram provided by the embodiment of the present application is shown in the flowchart.

[0118] Figure 5 A possible adaptive spoofing mitigation result schematic diagram provided by the embodiment of the present application is shown in the flowchart. DETAILED DESCRIPTION

[0119] The embodiments of the present application are described in detail below, and examples of the embodiments are shown in the drawings, wherein the same or similar reference numerals represent the same or similar elements or elements having the same or similar functions throughout. The embodiments described below by referring to the drawings are exemplary and are only used to explain the present application, and cannot be interpreted as a limitation on the present application.

[0120] As will be understood by one skilled in the art, the singular forms "a," "an," and "the" include plural referents unless the context clearly dictates otherwise. It should be further understood that the terms "comprise," "comprises," "comprising," "include," "includes," and "including" when used in this specification and in the following claims are taken to specify the presence of stated features, integers, steps, operations, elements, and / or components, but do not preclude 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 one element being "connected" or "coupled" to another element, it can be directly connected or coupled to the other element or intervening elements can also be present. In addition, the use of "connection" or "coupling" herein also includes wireless connection or coupling. As used herein, the term "and / or" includes any and all combinations of one or more of the associated listed items.

[0121] As will be understood by one skilled in the art, unless otherwise defined, all terms (including technical and scientific terms) used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this application belongs. It will be further understood that terms, such as those defined in commonly used dictionaries, should be interpreted as having a meaning that is consistent with their meaning in the context of the relevant art and will not be interpreted in an idealized or overly formal sense unless expressly so defined herein.

[0122] For the purpose of promoting an understanding of the principles of the present application, reference will now be made to the embodiments illustrated in the drawings and specific language will be used to describe the same. It will, nevertheless, be understood that no limitation of the scope of the application is thereby intended, such alterations and further modifications in the illustrated device being contemplated as falling within the scope of the application.

[0123] The present application provides a method for implementing spoofing resistance in the GNSS measurement domain. The present application aims to alleviate the impact of spoofing attacks on train satellite positioning, fully utilizes the spoofing-immune multi-source auxiliary information provided by the existing configuration and commonly used sensor equipment of the train, integrates spoofing detection, identification and mitigation techniques into a multi-source fusion positioning architecture based on factor graph optimization, and thus simultaneously realizes spoofing attack mitigation and real train state recovery.

[0124] The present application adopts a "detection-utilization approach" based on spoofing error compensation. Unlike existing methods, the present application further integrates spoofing detection, identification and mitigation into a train multi-source fusion architecture based on factor graph, to realize collaborative processing of spoofing mitigation and real train state recovery.

[0125] The application scenario of the factor graph optimization train positioning method against satellite navigation spoofing attack provided by the embodiment is as follows: when receiving a satellite signal information frame, linear interpolation is performed on a typical speed and distance sensor inertial measurement unit (IMU) / wheel axle speed and distance sensor (ODO) information frame to align the time stamp, and IMU / ODO pre-integration is performed; the IMU / ODO pre-integrated position and track information provided by a track digital terrain model (DTM) are used to construct a reliable bias for adaptive spoofing detection and identification, on the basis of which, hypothesis testing is implemented, a spoofing event is perceived to occur, and a spoofing satellite in a visible satellite is identified; if the spoofing event occurs, a spoofing envelope state component corresponding to the spoofing satellite identified in the factor graph system state node is activated, each factor is constructed and added to the factor graph model, and the factor graph automatically estimates and compensates for the spoofing envelope of the spoofing measurement in joint optimization, so that the real state of the train is restored; pseudo-range residuals are calculated according to the estimated optimal position, and are updated to a residual pool, a Gaussian mixture model (GMM) is used to fit the pseudo-range residual distribution, so that in the system cycle process, adaptive spoofing detection, identification and suppression are assisted to be realized.

[0126] The present application aims to provide a factor graph optimization train positioning method against satellite navigation spoofing attack, and timely perceive spoofing attack and accurately restore the real state of the train, so as to realize multi-source trusted positioning of the train in a spoofing attack scenario. The flowchart of the factor graph optimization train positioning method against satellite navigation spoofing attack provided by the embodiment of the present application is as shown in Figure 1 The present application aims to provide a factor graph optimization train positioning method against satellite navigation spoofing attack, and timely perceive spoofing attack and accurately restore the real state of the train, so as to realize multi-source trusted positioning of the train in a spoofing attack scenario. The flowchart of the factor graph optimization train positioning method against satellite navigation spoofing attack provided by the embodiment of the present application is as shown in

[0127] Step S1, determine and initialize a system state node with a spoofing envelope component, a factor graph model with spoofing mitigation capability and related parameters.

[0128] The factor graph model is mainly used in step S4 to realize spoofing mitigation and optimal estimation of the train running state. According to the anti-spoofing technology requirement and the sensor fusion positioning mode, the system state node and the factor graph model in the spoofing attack state are determined, and the prior factor of the first global node and the related parameters to be estimated are initialized, and the obtained factor graph model is used for subsequent spoofing mitigation and optimal estimation of the train running state.

[0129] S1.1, determine and initialize a system to-be-estimated state node with a spoofing envelope component.

[0130] The sliding window factor graph model is adopted to balance the computational cost and real-time performance of graph optimization. The sensors include GNSS (Global Navigation Satellite System), INS (Inertial Navigation System), ODO (Odometer) and DTM (Digital Track Map). The fusion mode is a tightly coupled mode with high precision. In the multi-source tightly coupled positioning, the original observation pseudo-range and Doppler shift of GNSS, the original observation values of three-axis gyroscope and accelerometer of INS, and the original observation wheel shaft pulse count of ODO are processed and sent to the factor graph optimization for joint optimization positioning. The typical states to be estimated include position, velocity, attitude, clock error and clock drift of the GNSS receiver, wherein the estimation of the clock error and clock drift of the GNSS receiver is used to correct the errors of the pseudo-range and pseudo-range rate (calculated from the Doppler shift). Considering the error accumulation characteristics of the inertial measurement unit (IMU), the IMU gyroscope and accelerometer bias also need to be considered and added to the estimated state to provide error correction. Similarly, due to the actual train wheel diameter loss and slipping, the bias is also introduced in the original measurement of ODO, and the scale parameter is also added to the estimated state. If there is a spoofing satellite in the visible satellite, the spoofing envelope state component of the spoofing satellite will be activated to compensate for the pseudo-range of the spoofing satellite. In order to facilitate multi-sensor fusion and coordinate conversion, the east-north-sky is selected as the world coordinate system (w system) with the initial position as the origin, which is also the reference coordinate system of the factor graph. Then in a time window, the multi-dimensional state node set χ of the system can be constructed as:

[0131]

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

[0133] S1.2, determine and initialize the prior factor of the global first node of the spoofing mitigation factor graph model.

[0134] The prior factor of the first global node determines the initial value of the whole system, is used for the initial state alignment of the system, and directly affects the convergence and convergence speed of the subsequent optimization. A reliable and accurate initial state prior information x0 can be given by the measurement value of a high-precision measurement device or an empirical value, and is one-to-one corresponding to the system state node, which can be expressed as:

[0135]

[0136] wherein, is the initialized position, velocity and attitude, and is the initialized IMU accelerometer and gyroscope bias, and is the initialized GNSS receiver clock error and clock drift, is the scale factor of the initialized 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 first global node can be expressed as

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

[0139] wherein, Σ0 is the noise covariance matrix of the prior factor of the first global node, representing the uncertainty of the initialization information, which can be given by the measurement standard of a high-precision measurement device or an empirical value.

[0140] Further, the error function of the prior factor of the first global node can be obtained

[0141]

[0142] According to the error function of the prior factor of the first global node, the factor graph model is calibrated.

[0143] S1.3, determine and initialize the factor graph model with spoofing mitigation capability.

[0144] Figure 2 A possible implementation schematic diagram of the factor graph optimization train positioning method for resisting satellite navigation spoofing attack in the embodiment is shown.

[0145] In traditional bundle adjustment systems, IMU pre-integration is used to improve computational efficiency. Before nonlinear optimization, IMU observations between adjacent nodes are pre-integrated to obtain relative position, velocity and attitude increments independent of the initial pose. In the optimization process, the IMU pre-integration factor combines several observations between adjacent nodes and provides a single relative motion constraint. The embodiment further combines the IMU measurement and the ODO mileage measurement of the train wheel forward direction to avoid the error introduced when the ODO mileage measurement is converted into speed, while satisfying the integral form of the nonlinear optimization problem. In the optimization process, the IMU / ODO pre-integration is the main body of the time sequence state recursion. The IMU / ODO pre-integration factor constructed based on the IMU and ODO measurement model is used to form a probability constraint on the adjacent epoch state. The GNSS pseudorange factor and the Doppler velocity factor construct a measurement model based on the observation of the visible satellite, and then constrain the related state node of the current epoch. In addition, the map matching factor constructed with the aid of the electronic track map DTM will constrain the train position to the inherent track line. The marginalization strategy is used to remove old nodes within the window, increase edge constraints and minimize the loss of historical observation information as much as possible.

[0146] In order to realize the multi-source reliable positioning of the train under the spoofing attack, spoofing countermeasures are added to the factor graph structure of the tightly coupled GNSS / INS / ODO / DTM positioning, including adaptive spoofing detection, identification and mitigation. The adaptive spoofing detection and identification mechanism acts before joint optimization, and uses the IMU / ODO pre-integration output, map matching and GMM error weight to assist in constructing reliable bias, performing binary hypothesis testing, and deciding whether there is a spoofing attack. After detecting the spoofing attack, the spoofing satellite in the visible satellite is determined by the standardized reliable bias. The adaptive spoofing elimination strategy activates the spoofing envelope node of the identified spoofing satellite, so as to directly optimize and compensate the spoofing measurement in the optimization process. The optimal estimated position is used to calculate the GNSS residual, and the GMM model is further sampled to fit the residual distribution, to support adaptive spoofing detection and identification, and GNSS adaptive factor construction.

[0147] Step S2, after receiving the speed and distance measurement sensor information frame, or after aligning the time stamp of the speed and distance measurement sensor information frame with linear interpolation, the speed and distance measurement sensor pre-integrated position is calculated.

[0148] A fine IMU pre-integration is implemented by compensating for the Earth's rotation, shape and gravity changes. The angular velocity and acceleration of the IMU can be modeled as:

[0149]

[0150] where ng and n a are the noises of the gyroscope and the accelerometer, respectively, is the direction cosine matrix from w-frame to b-frame, is its inverse, is the projection of the earth rotation angular velocity of the earth-centered earth-fixed frame (e-frame) in the w-frame with respect to the inertial frame (i-frame). is the dependent angular velocity of the w-frame due to the carrier motion and the earth curvature, and are the Coriolis acceleration and the centripetal acceleration due to the earth rotation and the carrier motion. is the coordinate transformation matrix from the navigation frame (n-frame) to the earth gravity in the w-frame is the projection of the earth gravity in the w-frame, which is usually related to the latitude and the altitude h of the location, and the gravity model is set as

[0151]

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

[0153]

[0154] where R w is the radius of the driving wheel. Ignoring the lateral and vertical motion of the contact point between the train wheel with the ODO installed and the ground, considering the influence of the wheel diameter loss, skidding, etc., the distance increment measured by the ODO needs to be corrected using a scale factor. Then the ODO measurement can be expressed in the following vector form

[0155]

[0156] The pre-integrated model of the IMU / ODO at two adjacent optimization instants is

[0157]

[0158] where, is the increment of the ODO in the b-frame, is the velocity in the b-frame, a , ι g and ι odo are the Gaussian white noises of the IMU accelerometer, the IMU gyroscope and the ODO scale factor, and Ω is a composite matrix whose first column is the quaternion the first row is the transpose of the conjugate quaternion of and the rest is the skew-symmetric matrix composed of .

[0159] Step S3, determine the nearest orbit segment in the orbit information, map match the pre-integrated position, calculate the reliable bias, construct the adaptive detection and identification statistics, perceive the spoofing event and identify the spoofed satellite in the visible satellites.

[0160] Before the optimization state estimation, the adaptive spoofing detection and identification method is implemented to confirm the spoofing event and identify the satellite measurements that are spoofed. First, the nearest orbit segment to the IMU / ODO pre-integrated predicted position is determined according to the distance association model, the predicted position is projected onto the orbit segment using the vertical projection principle to determine the on-track position; then, the reference pseudorange and the reliable bias are calculated according to the on-track position and the spatial coordinates of the visible satellites, and then the detection statistics and threshold are calculated, the binary hypothesis test is implemented to perceive whether the spoofing event occurs; after confirming that the spoofing event occurs, the standardized reliable bias is used to identify the satellites in the visible satellites that are affected by spoofing.

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

[0162] The distance association model is constructed to extract the track segment currently occupied by the train. For ease of expression, the electronic track map is simplified to a series of interest points in the e system The adjacent interest points are the end point coordinates of the track segment. The distance association model determines the track segment occupied by the train by judging the proximity of the IMU / ODO pre-integrated predicted position to the track segment.

[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 project the predicted position onto the track line to correct the vertical track error. The projected position can be calculated as

[0166]

[0167] S3.2 Calculate the reliable bias, perceive the spoofing event and identify the spoofed satellite in the visible satellites.

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

[0169]

[0170] The difference between the reference pseudo-range and the satellite corrected observed pseudo-range is defined as the reliable bias, which can be calculated as

[0171]

[0172] For a real satellite, the reliable bias only contains the un-modeled pseudo-range error ε k,i , while when a spoofing event occurs, a spoofing envelope s k,i will be introduced into the bias of the spoofing satellite. Based on the reliable bias, a Weighted Sum of the Squared Errors (WSSE) can be constructed as a spoofing detection statistic.

[0173]

[0174] where l is the number of visible satellites at the current epoch, W k is an adaptive weight matrix, whose diagonal elements are the inverse of the satellite residual distribution variance derived from a Gaussian Mixture Model. Under the condition of no spoofing, the WSSE detection statistic follows a chi-square distribution with degree of freedom l. While under the spoofing attack, the WSSE detection statistic will not follow the chi-square distribution. According to the pre-set false alarm rate, the detection threshold T h is calculated, then the binary hypothesis test for adaptive spoofing detection is

[0175]

[0176] Under the spoofing attack, the spoofing envelope is introduced into the bias vector, which can make the detection statistic deviate from the normal state distribution and more easily exceed the detection threshold, so as to quickly respond to the spoofing event. For a more covert spoofing attack, in order to avoid arousing the vigilance of the target receiver in the invasion stage, a pseudo-range offset small enough or even 0 is often set in the early invasion stage, which leads to a weak spoofing envelope of the spoofing satellite, making it difficult to support fast response and even leading to the failure of spoofing detection. A successful spoofing attack will inevitably cause changes in the signal amplitude received by the target receiver, which is directly reflected in the fluctuation of the carrier-to-noise power spectral density C / N0 of the satellite signal. Therefore, in order to enhance the detection ability of adaptive spoofing detection to weak spoofing attacks, C / N0 is further introduced into the weight matrix of the detection statistic to make up for the detection defects caused by the weak spoofing envelope. The C / N0 enhanced weight matrix can be constructed as

[0177]

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

[0179] After the spoofing event is perceived, adaptive spoofing identification is further performed to locate the spoofing satellite in the visible satellite set. The normalized reliability bias can be calculated as

[0180]

[0181] In the visible satellite set, the satellite with the largest normalized reliability bias is determined as the spoofing satellite, and its spoofing flag is set as Tag = 1. Subsequently, by deleting the component corresponding to the spoofing signal in the reliability bias a k and the adaptive weighting matrix W k , a k and W k are updated, the WSSE detection index is recalculated, and spoofing detection is performed again until the updated WSSE index is lower than the detection threshold. The algorithm terminates, considering that all spoofing satellites have been identified, and the remaining signals are determined as real signals, with their spoofing flags set as Tag = 0.

[0182] Figure 3 A possible WSSE index statistical result for adaptive spoofing detection is shown, Figure 4 A possible normalized reliability bias statistical result for adaptive spoofing identification is shown.

[0183] Step S4, activate the spoofing envelope state corresponding to the identified spoofing satellite, construct adaptive pseudorange factors of real / spoofing satellites, pre-integration factors of velocity and range sensors, map matching factors, and prior factors, jointly optimize the factor graph, and preliminarily solve the positioning solution of the train.

[0184] S4.1 Construct the GNSS adaptive factor and add it to the factor graph model.

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

[0186]

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

[0188] Under spoofing attack, the spoofing included in the pseudorange observation equation can be modeled as

[0189]

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

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

[0192]

[0193] where, and are the GNSS measurements corrected by equations (18), (19) and (20), and are the covariance matrices of the true pseudorange, pseudorange rate and spoofed pseudorange, which are generally modeled as zero-mean Gaussian distribution.

[0194] The covariance matrix of GNSS factor represents the uncertainty of the measurement, which is crucial to the state estimation accuracy of the factor graph. The mismatched measurement error distribution will increase the difficulty of solving the optimal solution of the factor graph. There are stations, road cuts, urban canyons and other restricted environments along the railway, which may cause the GNSS measurement error distribution to present non-Gaussian distribution phenomena such as multi-peak and long tail. And under the influence of spoofing attack, the non-Gaussian distribution characteristics of measurement error will be further enhanced. The single Gaussian distribution assumption will degrade the state estimation performance, resulting in 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] where, and are the adaptive variances derived by Gaussian mixture model (GMM).

[0197] S4.2 Build a velocity and distance sensor pre-integration factor and add it to the factor graph model.

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

[0199]

[0200] where, is the covariance matrix of IMU / ODO measurement z IMU / ODO , which represents the uncertainty of IMU / ODO pre-integration equivalent measurement correlation. and IMU / ODO pre-integration result corrected by first order approximation of equation (9).

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

[0202] The track index and in-track position closest to the current state have been determined in step S3.1. The vertical track error is the distance between the positioning position and the in-track position, which describes the degree of deviation of the positioning position from the track. Smaller vertical track error is more likely to represent more accurate positioning solution. Therefore, the map matching factor is constructed with the vertical track error, which is added to the optimization process of the factor graph to directly constrain the position. The map matching factor can be calculated as

[0203]

[0204] where, 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 Construct the prior factor within the sliding window and add it to the factor graph model.

[0206] When the number of states within the sliding window exceeds the window capacity, the oldest state will be marginalized, and the sensor measurement corresponding to the marginalized state will be converted into a prior factor, and its cost function is

[0207]

[0208] where, is the fixed linearization point of the first system state node within the sliding window, J prior and σ prior are the Jacobian matrix and residual of the prior factor.

[0209] S4.5 Joint optimization, preliminary solution of the running state of the train.

[0210] Finally, the preliminary estimation of train running state is obtained by jointly optimizing GNSS true pseudorange adaptive factor, GNSS spoofing pseudorange adaptive factor, GNSS Doppler velocity factor, IMU / ODO pre-integration factor, map matching factor and priori factor The global error function can be minimized to obtain

[0211]

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

[0213] Step S5, calculate the pseudorange residuals according to the preliminary positioning solution, compensate the residuals of the spoofing satellites, update the residual pool, fit the residual distribution of each visible satellite in the residual pool, calculate the adaptive weight, and feed back to the adaptive spoofing detection, identification and mitigation step to continue iteration to obtain the optimal positioning solution of the train.

[0214] Figure 5 A possible optimal state estimation result is provided for adaptive spoofing mitigation. Gaussian Mixture Model (GMM) is an extension of single Gaussian model, which smoothes the fitting of density distribution of arbitrary shape by combining multiple Gaussian components with different weights. In order to balance the fitting accuracy and the calculation timeliness, Gaussian Mixture Model with a maximum of three Gaussian components is used to model the core of GNSS residual distribution and the possible tail. According to the optimal state obtained by joint optimization of the factor graph, the pseudorange residual of the i-th satellite can be calculated as:

[0215]

[0216] wherein, and are the optimal position, receiver clock bias and spoofing envelope obtained by solving. For the identified spoofing satellite, its residual will be corrected using the estimated spoofing envelope.

[0217] The latest calculated GNSS residual is added to the residual pool, and when the residual items in the residual pool exceed the capacity, the earliest added residual item will be removed. Then, GMM is used to fit the residual set of each visible satellite. Since the residual of the spoofing satellite is close to 0 value after being corrected by the spoofing envelope, its distribution will be significantly different from that of the real satellite. If GMM fitting is used, the measurement weight of the spoofing satellite will be excessively enlarged in the positioning calculation, which greatly weakens the contribution of the real satellite to the positioning calculation. To avoid such a situation, the residual of the identified spoofing satellite will not be added to the residual pool.

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

[0219]

[0220] wherein, are parameters of the Gaussian Mixture Model, μ and Σ are the weight, mean and variance of the Gaussian component.

[0221] If variational inference is used to solve the parameters of the Gaussian Mixture Model, the variational E step (Expectation) and the variational M step (Maximization) are performed iteratively until the likelihood function converges or a predefined maximum iteration threshold is reached, thereby obtaining the optimal GMM parameters Further, the adaptive variance of the satellite pseudorange observation is calculated as

[0222]

[0223] The covariance matrix formed by the same will be fed back to the adaptive spoofing detection and identification, and the construction of the GNSS pseudorange adaptive factor, and the global error function formula (29) of the factor graph is continuously iteratively optimized to obtain the optimal solution of the train running state.

[0224] In summary, the existing spoofing mitigation method of the embodiment (1) generally uses least squares optimization and filter estimator as a state solver, and the anti-spoofing technology of the fusion positioning based on factor graph optimization needs to be further studied. The present application realizes a factor graph optimization method against satellite navigation spoofing attacks, integrates the spoofing detection, identification and mitigation mechanism into the multi-source tightly coupled factor graph optimization architecture, and in the process of system state optimization solving, the method not only estimates and compensates the spoofing envelope of the spoofing signal in real time, but also realizes the accurate recovery of the real running state of the train.

[0225] (2) The current research on satellite positioning anti-spoofing methods for railway specific scenarios is still in the initial exploration stage, and the present application is specially developed for the satellite navigation spoofing threat problem in the train running environment, fully utilizes the train speed and distance measurement sensors and prior track information, and realizes the rapid recovery of state estimation under spoofing attack through optimization integration strategy to construct auxiliary constraints. On the basis of not significantly increasing the hardware cost of the system, the anti-spoofing positioning ability of the train is effectively improved. The anti-spoofing method provided by the present application is highly compatible with the railway application scenario, and can be widely applied to various applications based on satellite navigation driven train positioning, and has good engineering feasibility and promotion prospect.

[0226] Those skilled in the art can understand that the drawings are only schematic diagrams of an embodiment, and the modules or processes in the drawings are not necessarily necessary for implementing the present application.

[0227] Those skilled in the art can clearly understand the present application by the description of the above embodiments. Based on such an understanding, the technical solutions of the present application can be embodied in the form of a software product, and the computer software product can be stored in a storage medium, such as a ROM / RAM, a magnetic disk, an optical disk, and the like, and includes a plurality of instructions to cause a computer device (which can be a personal computer, a server, or a network device, and the like) to execute the methods described in the various embodiments or some parts of the embodiments.

[0228] Each of the embodiments in the specification is described in a progressive manner, and the same or similar parts between the embodiments can be referred to each other. Each of the embodiments mainly describes the difference 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, and the related parts can be referred to the part of the method embodiments. The above-described device and system embodiments are merely illustrative, and the units described as separate components can be or can not be physically separated, and the components displayed as units can be or can not be physical units, that is, they can be located in one place, or can be distributed on multiple network units. Part or all of the modules can be selected to achieve the purpose of the embodiments according to the actual needs. Those skilled in the art can understand and implement it without creative labor.

[0229] The above is only the preferred embodiment of the present application, but the protection scope of the present application is not limited to this. Any skilled person in the art can easily think of changes or replacements within the technical range disclosed by the present application, which should be covered in the protection scope of the present application. Therefore, the protection scope of the present application should be subject to the protection scope of the claims.

Claims

1. A factor graph optimization train positioning method against spoofing attacks of satellite navigation, characterized in that, The method comprises the following steps: Determine and initialize the system state node with the spoofing envelope component, the factor graph model with the spoofing mitigation capability, and the related parameters; Upon receiving the speed and distance measurement sensor information frame, or upon receiving the satellite signal information frame, the speed and distance measurement sensor pre-integrated position is calculated by linear interpolation alignment timestamp; Determine the nearest track segment in the track information, map match the pre-integrated position, calculate the reliable bias, construct the adaptive detection and identification statistics, perceive the spoofing event and identify the spoofing satellite in the visible satellite; Activate the spoofing envelope state corresponding to the identified spoofing satellite, construct the adaptive pseudo-range factor of the real / spoofing satellite, the speed and distance measurement sensor pre-integrated factor, the map matching factor, and the prior factor, jointly optimize the factor graph, and preliminarily solve the positioning solution of the train; According to the preliminary positioning solution, the pseudo-range residual is calculated, the residual of the spoofing satellite is compensated, the residual pool is updated, the residual distribution of each visible satellite in the residual pool is fitted, the adaptive weight is calculated and fed back to the adaptive spoofing detection, identification and mitigation step, and the optimal positioning solution of the train is continuously obtained by iteration.

2. The method of claim 1, wherein, The method comprises the following steps: S1.1, determine and initialize the system state node with the spoofing envelope component; In multi-source tightly coupled positioning, the raw observations of pseudorange and Doppler shift of global navigation satellite system (GNSS), the raw observations of three-axis gyroscope and accelerometer of inertial navigation system (INS), and the raw observations of wheel axle pulse count of wheel axle speed sensor (ODO) are processed and sent into factor graph optimization for joint optimization positioning. If there is a spoofing satellite in the visible satellites, the spoofing envelope state component of the spoofing satellite will be activated. The east-north-sky with the initial position as the origin is selected as the world coordinate system, i.e., the w system, which is also the reference coordinate system of the factor graph. In a time window, the multi-dimensional state node set χ and the state node x of the system are constructed as: k ​ where x k denotes the 18-dimensional state node, and is the three-dimensional position and velocity of the train, is the train pose in quaternion form, and are the three-dimensional accelerometer and gyroscope biases of the IMU carrier, and are the clock error and clock drift of the GNSS receiver, is the scale factor of the ODO, s k is the spoofing envelope, m is the number of spoofed satellites identified, and n is the size of the sliding window. S1.2, determine and initialize the prior factor of the global first node of the spoofing mitigation factor graph model The initial state prior information x0 is given by the measurement value or the empirical value of the measurement device, and is one-to-one corresponding to the system state node, and is expressed as: wherein, are the initialized position, velocity, and attitude, and are the initialized IMU accelerometer and gyroscope biases, and are the initialized GNSS receiver clock bias and clock drift, is the initialized scale factor for 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: Where ∑0 is the noise covariance matrix of the prior factor of the global first node, indicating the uncertainty of the initialization information; The error function of the prior factor of the global first node is obtained The factor graph model is calibrated according to the error function of the prior factor of the global first node; S1.3, determine and initialize the factor graph model with the spoofing mitigation capability; The inertial measurement unit IMU measurement and the mileage measurement of the train wheel forward direction provided by the ODO are combined, in the optimization process, the IMU / ODO pre-integration is used as the main body to perform time sequence state recursion, the IMU / ODO pre-integration factor constructed based on the IMU and ODO measurement model is used to form a probability constraint on the adjacent epoch state, the GNSS pseudo-range factor and the Doppler velocity factor are used to construct a measurement model based on the observation of the visible satellite to constrain the related state node of the current epoch, and the map matching factor constructed with the aid of the electronic track map DTM will constrain the train position to the inherent track line. In the factor graph structure of the tightly coupled GNSS / INS / ODO / DTM positioning, spoofing countermeasures are added, including adaptive spoofing detection, identification and mitigation, the IMU / ODO pre-integrated output, map matching and GMM error weight are used to assist in building reliable bias, binary hypothesis test is carried out to decide whether there is a spoofing attack, after detecting the spoofing attack, the spoofing satellite in the visible satellite is determined by the standardized reliable bias, the adaptive spoofing elimination strategy activates the spoofing envelope node of the identified spoofing satellite, directly optimizes and solves in the optimization process and compensates the spoofing measurement, and determines and initializes the factor graph model with spoofing mitigation capability.

3. The method of claim 2, wherein, The linear interpolation timestamp alignment is performed on the speed and distance sensor information frame when the speed and distance sensor information frame is received, or when the satellite signal information frame is received, and the speed and distance sensor pre-integration position is calculated, comprising: The speed and distance sensor includes an inertial sensor IMU and an axle speed and distance sensor ODO; The angular velocity of the IMU and the acceleration are modeled as: where n g and n a are the noise of the gyroscope and accelerometer respectively, is the direction cosine matrix of w-frame to IMU-frame, is its inverse, is the dependent angular velocity of w-frame due to the carrier motion and the earth curvature, is the projection of the earth rotation angular velocity of e-frame to i-frame in w-frame, and are the Coriolis acceleration and centripetal acceleration due to the earth rotation and the carrier motion, is the earth gravity in n-frame, is the projection of the earth gravity in w-frame, the earth gravity at different position is usually related to the latitude and the altitude h, the model of the earth gravity is set as: Pulse count according to ODO Calculating one-dimensional mileage increments of the train along the track within two consecutive time stamps wherein R w is the radius of the drive wheel and π is the circle constant. The ODO measurement is expressed in the following vector form: At two adjacent optimization periods [t i-1 , t i ], the pre-integrated model of IMU / ODO is: where, is the velocity under b, i a , g and i odo are IMU accelerometer, IMU gyroscope, and ODO scale factor Gaussian white noise, Ω is a composite matrix whose first column is the quaternion The first row is the transpose of the conjugate quaternion of , and the rest is the skew-symmetric matrix composed of .

4. The method of claim 3, wherein, The reliable bias is calculated by determining the nearest track in the track information and performing map matching on the pre-integrated position, and the adaptive detection and identification statistics are constructed to perceive the spoofing event and identify the spoofing satellite in the visible satellite, comprising: S3.1 performing map matching to determine the in-orbit position of the speed and distance sensor pre-integrated predicted position; Extract the track segment currently occupied by the train, simplify the track electronic map into a series of points of interest under e The adjacent point of interest is the end point coordinate of the track segment, and the IMU / ODO pre-integrated prediction position is judged The proximity to the track segment is used to determine the track segment occupied by the train Wherein, id represents the index of the point of interest; After confirming the train occupies the track segment index, the train predicted position Project to the track line, train projection position The calculation is: S3.2 calculating the reliable bias, perceiving the spoofing event and identifying the spoofing satellite in the visible satellite; The reference pseudorange of the i th visible satellite is calculated as: The difference between the reference pseudo-range and the satellite corrected observed pseudo-range is defined as the reliable bias a , calculated as: k，i , calculated as: For real satellites, the reliable bias only contains un-modeled pseudorange errors ε k，i While a spoofing event occurs, a spoofing envelope s k，i A weighted residual sum of squares is constructed based on the reliable bias as a spoofing detection statistic. where l is the number of visible satellites at the current epoch, W k is the adaptive weighting matrix, whose diagonal elements are the inverse of the variance of the satellite residual distribution derived from the Gaussian mixture model. Under the condition of no spoofing, the WSSE detection index obeys the chi-square distribution with degree of freedom l, while under the spoofing attack, the WSSE detection index will not obey the chi-square distribution. According to the pre-set false alarm rate, the detection threshold T is calculated h Then the binary hypothesis test for adaptive spoofing detection is: The carrier noise power spectral density C / N0 of the satellite signal is introduced into the weight matrix of the detection statistics, and the C / N0 enhanced weight matrix W is constructed as: where γ (ω1≤ γ ≤ 1) is a quadratic function derived from the visible satellite C / N0 about the center symmetric; ω1and ω2are adjustment factors to control the minimum value and gradient of the function; After perceiving the spoofing event, adaptive spoofing identification is further performed to locate the spoofing satellite in the visible satellite, and the standardized reliable bias is calculated as: In the visual satellite set, the satellite with the largest normalized reliability bias is determined as the spoofing satellite, and its spoofing flag is set as Tag = 1, and the reliability bias a k and the adaptive weighting matrix W k corresponding to the spoofing signal in the adaptive weighting matrix W k and W k , and the WSSE detection index is recalculated, and the spoofing detection is performed again until the updated WSSE index is lower than the detection threshold, the algorithm is terminated, and it is considered that all spoofing satellites are identified, and the remaining signals are determined as real signals, and the spoofing flag is set as Tag = 0.

5. The method of claim 4, wherein, The spoofing envelope state corresponding to the identified spoofing satellite is activated, the adaptive pseudorange factor of the real / spoofing satellite, the speed and distance sensor pre-integration factor, the map matching factor and the prior factor are constructed, the factor graph is optimized, and the positioning solution of the train is solved, comprising: S4.1 constructing the GNSS adaptive factor and adding it to the factor graph model; The position, velocity, clock bias and clock drift of the system state node at the current epoch are constrained by using GNSS factors. In the non-spoofing scenario, the pseudorange and Doppler frequency shift d k，i are modeled as: where the upper index A denotes the true satellite information, and is the spatial position, velocity, clock bias and clock drift of the visible satellites, and λ is the wavelength, is the line-of-sight vector between the receiver and the satellite, and is the ionosphere, troposphere and earth rotation induced equivalent pseudorange error, and is the corresponding equivalent pseudorange rate error, and are unmodeled pseudorange and pseudorange rate errors; h Aρ () and h d () are the observation models for the true pseudorange and Doppler frequency shift, respectively; Under spoofing attack, the pseudorange observation equation with spoofing envelope is modeled as: where the upper index S denotes the information of the spoofing satellite; h Sρ () is the observation model of the spoofed pseudorange; s k，i is the spoofing envelope of the spoofed pseudorange; Then the error function of the GNSS pseudorange and Doppler velocity factor is calculated as: wherein, and are GNSS measurements corrected through equations (18), (19) and (20), and are the real pseudorange, pseudorange rate and spoofed pseudorange corresponding covariance matrices, modeled as zero-mean Gaussian distributions; The GNSS adaptive factor with adaptive error modeling is introduced for factor graph optimization, and the GNSS adaptive factor is calculated as: wherein, and are adaptive variances of real and spoofed satellite pseudoranges derived from a Gaussian Mixture Model, GMM. S4.2 Given the update period ΔT of the IMU / ODO pre-integration factor, the error function of the IMU / ODO pre-integration factor at k to k-1 is calculated through the state transition of the error term wherein is the IMU / ODO measurement z IMU / ODO is the covariance matrix of the IMU / ODO pre-integrated equivalent measurement, representing the uncertainty related to the IMU / ODO pre-integrated equivalent measurement, and is the IMU / ODO pre-integrated result corrected by the first order approximation of equation (9); and is defined in equation (1) for the state node components at time k, representing the position, velocity, attitude, accelerometer bias, gyroscope bias and the scale factor of ODO, respectively. S4.3 The map matching factor is calculated as: wherein, is the variance associated with the map matching factor; S4.4 The prior factor in the sliding window is constructed, and the prior factor is added 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 measurement corresponding to the marginalized state will be converted into a prior factor, and the cost function is: wherein is a fixed linearization point of the first system state node within the sliding window, J prior and σ prior are Jacobian matrix of the prior factor and the residual; S4.5 The preliminary estimation of train running state is obtained by combining GNSS real pseudorange adaptive factor, GNSS spoofing pseudorange adaptive factor, GNSS Doppler velocity factor, IMU / ODO pre-integration factor, map matching factor and priori factor By minimizing the global error function, we get: The factor graph solving method used is the Levenberg-Marquardt method in Ceres Solver.

6. The method of claim 5, wherein, The pseudo-range residuals are calculated according to the preliminary positioning solution, and the residuals of the spoofed satellites are compensated, the residual pool is updated, the residual distribution of each visible satellite in the residual pool is fitted, the adaptive weight is calculated and fed back to the adaptive spoofing detection, identification and mitigation step, and the optimal positioning solution of the train is obtained by continuous iteration, including: The Gaussian Mixture Model with at most three Gaussian components is used to model the core and possible tail of the GNSS residual distribution, and the optimal state is solved by joint optimization according to the factor graph, the pseudo-range residual of the ith satellite is calculated as: wherein, and For the optimal position, receiver clock bias and spoofing envelope to be solved, the residual of the identified spoofing satellite will be corrected using the estimated spoofing envelope; The latest calculated GNSS residual is added to the residual pool, when the residual items in the residual pool exceed the capacity, the earliest added residual item will be removed, the residual set of each visible satellite is fitted by GMM, since the residual of the spoofed satellite is close to 0 value after being corrected by the spoofing envelope, the residual of the identified spoofed satellite will not be added to the residual pool; The probability density function of the Gaussian Mixture Model of the ith satellite pseudo-range residual at time k is modeled as: wherein are parameters of a Gaussian mixture model, are weights of the Gaussian components, μ and ∑ are mean and variance; If variational inference is used to solve the parameters of the Gaussian mixture model, the variational E-step and the variational M-step are performed iteratively until the likelihood function converges or a predefined maximum iteration threshold is reached, and the optimal GMM parameters are obtained Further, the adaptive variance of the satellite pseudo-range observation is calculated as In the iterative positioning solution, the adaptive variance will be fed back to the adaptive spoofing detection, identification and mitigation, and the optimal estimation of the train running state is obtained by continuous iteration of formula (29).

Images