An INS-aided GNSS pseudorange gross error detection method

Abstract

The application discloses an INS-assisted GNSS pseudorange gross error detection method; the role of INS navigation information in assisting in resisting errors is fully tapped, and the resistance to gross error of the integrated navigation system is greatly enhanced; firstly, GNSS single-point pseudorange observation values and Doppler observation values are combined with INS, and real-time estimation of satellite receiver clock difference and clock drift is monitored by using a sliding window; thereafter, the real-time updated position of the INS in the integrated system is substituted into satellite observation equations for pre-processing adjustment solution, and residual error checking is performed on the remaining unknowns; in the next step, through the related threshold setting of residual error checking and sliding window monitoring, the faulty satellites with the maximum residual error and exceeding the threshold are cyclically removed through a forward search mode. Finally, the M-LS filtering mode is adopted, and correct positioning results are obtained through GNSS / INS tight combination solution again.

Description

Technical Field

[0001] The present invention belongs to the field of navigation and positioning technology, and in particular relates to an INS-assisted GNSS pseudorange gross error detection method. Background Art

[0002] The rapid development of autonomous driving, drones, robotics, and other fields has led to higher demands for accuracy and robustness. The Global Navigation Satellite System (GNSS) provides users with all-weather, continuous, high-precision position, velocity, and timing information. GNSS boasts significant advantages such as low cost, high accuracy, and bounded errors. Multi-constellation and multi-frequency GNSS positioning offers a more reliable and accurate solution. However, GNSS is an active positioning method, and the environment significantly impacts GNSS signals. Issues such as multipath and non-line-of-sight in complex urban environments, as well as signal interruption in tunnels, pose challenges to continuous, high-precision GNSS positioning. INS (Inertial Navigation System) is a fully independent navigation and positioning method, offering advantages such as short-term high accuracy, high frequency, and immunity to external interference. However, factors such as inertial device drift can cause errors to accumulate over time. GNSS and INS have distinctly complementary characteristics, making GNSS / INS integrated navigation the most widely used and mature integrated navigation method.

[0003] In practical applications of integrated navigation, anomalies are often caused by observational environmental stimuli. This means that GNSS observation data quality is severely disrupted. However, inertial navigation has a strong degree of autonomy, so anomalies often occur in the measurement model. GNSS has a well-established integrity theory for assessing and monitoring GNSS failures. Receiver Autonomous Integrity Monitoring (RAIM) is the most widely used and inexpensive method for receiver-side positioning quality control. It primarily uses a statistical consistency test based on redundant pseudorange observations. While it performs well for single-system single-failure scenarios, it suffers from a high miss-detection rate for multiple-system multiple-failure scenarios. While some improved RAIM algorithms have been proposed to address this problem, these algorithms also place high demands on redundant observations and the distribution of visible satellites. Given these challenges, INS-assisted GNSS pseudorange gross error detection methods can be introduced to improve positioning accuracy and robustness in complex and large-scale scenarios such as urban canyons and tunnels. INS assistance can reduce the requirement for redundant observations and the correlation between observations, potentially addressing the problem of multiple-system multiple-failure scenarios. Summary of the Invention

[0004] In order to solve the above problems, the present invention discloses an INS-assisted GNSS pseudorange gross error detection method to improve the accuracy and robustness of positioning in complex large-scale scenes such as urban canyons and tunnels. The method uses GNSS single-point pseudorange observations and Doppler observations to tightly combine with INS, and performs sliding window monitoring on the real-time estimation of satellite receiver clock error and clock drift. Afterwards, the real-time updated solution position of INS in the combined system is substituted into the satellite observation equation for pre-processing adjustment solution, and the remaining unknowns are subjected to residual test. Next, through the relevant threshold setting of residual test and sliding window monitoring, the faulty satellites with the largest residual and exceeding the threshold are cyclically eliminated by forward search. Finally, the M-LS filtering method is adopted to perform GNSS / INS tight combination solution again to obtain the correct positioning result, which has high engineering application value.

[0005] To achieve the above object, the technical solution of the present invention is as follows:

[0006] An INS-assisted GNSS pseudorange gross error detection method is as follows:

[0007] (1) Tight combination of GNSS single-point pseudorange observations and Doppler observations with INS

[0008] The INS position, velocity, and attitude error equations based on the psi-angle (ψ) error model are:

[0009]

[0010]

[0011]

[0012] In the formula, the superscript n represents the navigation coordinate system, the subscripts i, e, and b represent the inertial coordinate system, the earth coordinate system, and the carrier coordinate system, respectively. n represents the position vector; represents the transfer rate; v n represents the carrier speed; ψ represents the platform misalignment angle; f n Indicates the relative force under the navigation system; represents the Earth's rotation rate; Represents the gravity calculated by the model; represents the transformation matrix from the carrier system to the navigation system; f b Indicates the specific force under the load system; It represents the angular velocity of the navigation system relative to the inertial system; It represents the angular velocity of the carrier system relative to the inertial system.

[0013] The construction of the tightly combined filter state vector is determined by the error state of the inertial device and the error state of the GNSS:

[0014]

[0015] Where, the superscripts n and b represent the navigation coordinate system and the carrier coordinate system respectively. n Indicates the misalignment angle of the three-dimensional platform under the navigation system; δv n represents the three-dimensional velocity error in the navigation system; δr n Represents the three-dimensional position error in the navigation system; ε b Indicates the three-dimensional gyro drift under the load system; Indicates the zero bias of the three-dimensional accelerometer under the load system; dt G Indicates the GNSS receiver clock offset; dt d Indicates the GNSS receiver clock drift.

[0016] The system error dynamics equation is shown below:

[0017]

[0018] Where F represents the state transfer matrix, which can be obtained by equations [1]-[3]; G represents the noise distribution matrix; and w represents the system noise.

[0019] Single-point pseudorange observations and Doppler observations are tightly combined with INS, where the ionospheric delay and tropospheric delay of the GNSS observations are corrected using the Klobuchar and Saastamoinen models, respectively. The measurement model is expressed as follows:

[0020] z k =H k x k +v k [6]

[0021] Where z k represents the measurement vector; subscript k represents the kth epoch; H k represents the measurement matrix; x k represents the state vector; v k represents the measurement noise, which obeys the zero-mean Gaussian distribution.

[0022] The measurement vector is the original GNSS observation value m GNSS and INS predicted value The difference between them is shown below:

[0023]

[0024]

[0025] Where, P s,f 、 They represent the original pseudorange observation value and Doppler observation value respectively; the superscript f represents the frequency; the superscript s represents the satellite system, including GPS, BDS, and Galileo; ρ INS 、 Represents INS predicted pseudorange and predicted pseudorange rate respectively; △dtr P 、 are the sum of the error corrections associated with the pseudorange and Doppler observations from the receiver clock; and are the sum of other error corrections for the pseudorange and Doppler observations.

[0026] The measurement model is expressed as follows:

[0027]

[0028]

[0029]

[0030]

[0031]

[0032] I=[1 1 ... 1] T

[14]

[0033] Where, e s,f represents the direction cosines from the receiver to the satellite; Represents the position error conversion matrix from the navigation coordinate system to the earth coordinate system; Represents the velocity error conversion matrix from the navigation coordinate system to the earth coordinate system.

[0034] The sliding window detection of the real-time estimation of the receiver clock error and clock drift is to calculate the corresponding sample mean and sample mean square error of the real-time estimation value within the sliding window. The calculation formula is as follows:

[0035]

[0036]

[0037] Where, represents the sample mean, and S represents the sample mean square error.

[0038] (2) INS-assisted GNSS residual test

[0039] The position parameters calculated by the inertial navigation update are X = (x, y, z) in the Earth-centered Earth-fixed coordinate system. Substitute it into the following pseudorange observation equation:

[0040]

[0041] Where, the superscript represents the i-th satellite and the subscript r represents the receiver; represents the pseudorange observation value; represents the geometric distance between the satellite and the earth; c represents the speed of light; δt r represents the receiver clock error; δt i represents the satellite clock error; represents the equivalent tropospheric delay; represents the equivalent ionospheric delay; represents the equivalent orbit error; Represents random error. Among them, the error terms such as tropospheric and ionospheric delays and satellite clock errors are corrected using the corresponding models, leaving only one unknown parameter in the equation, the receiver clock error. The above equation is adjusted to obtain δt r The estimates and the corresponding residuals of each equation are as follows:

[0042]

[0043] Where, the superscript represents the i-th satellite; v i represents the residual; l i represents the error of the i-th satellite;

[0044] (3) Forward search cycle to eliminate faulty satellites

[0045] Since sliding window monitoring is adopted, if the sample mean and sample mean square error of the sliding window monitoring exceed the set threshold, a residual test is performed and a forward search cycle is started to eliminate the faulty satellite.

[0046] As mentioned above, it can be seen from

[18] that the residual v of the observation value corresponding to each satellite is i It will be mainly affected by its own error; if there is no gross error for the i-th satellite, but there is a gross error for the j-th satellite, then the corresponding residual v of the j-th satellite is j The impact is The impact on the i-th satellite is only When there are a large number of satellites, there is a strong reason to believe that the observation value corresponding to the satellite with the largest residual that exceeds the threshold is the satellite with a gross error. At this time, this satellite is eliminated, and the observation equation shown in formula

[17] is reconstructed but does not include the jth satellite. This is repeated until the residual meets the requirements. When the number of satellites is small, the sample mean of the sliding window monitoring of the previous output epoch is also substituted into formula

[17] for residual detection, and the satellite with the largest residual is cyclically eliminated until the residual meets the requirements.

[0047] (4) GNSS / INS tight integration solution based on M-LS filtering

[0048] Assuming that the components of the observation vector of the epoch are independent of each other, but may contain abnormal errors and obey the contaminated normal distribution, while the state vector predicted by the dynamic model still obeys the normal distribution, the observation vector is estimated using the robust M method, and the state parameters are still estimated using the least squares (LS) method.

[0049] The M-LS filter extreme value condition is defined using the equivalent weight matrix, and the recursive solution is obtained as follows:

[0050]

[0051] Where K MLS It can still be called the gain matrix, and the expression is:

[0052]

[0053] Where, Represents the equivalent weight matrix of the observation vector, using the IGGIII weight function:

[0054]

[0055] Where, The weight matrix R k The i-th and j-th elements of γ ij for:

[0056]

[0057] in:

[0058]

[0059] Where, represents the normalized innovation, and k0 and k1 are the corresponding thresholds. This can produce a more accurate robust positioning result after eliminating GNSS gross errors.

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

[0061] The present invention discloses an INS-assisted GNSS pseudorange gross error detection method, which utilizes GNSS single-point pseudorange observations and Doppler observations to perform a tight combination with INS, and performs sliding window monitoring on the real-time estimation of satellite receiver clock error and clock drift. Afterwards, the real-time updated solution position of the INS in the combined system is substituted into the satellite observation equation for pre-processing adjustment solution, and a residual test is performed on the remaining unknowns. Next, by setting the relevant thresholds for the residual test and sliding window monitoring, the faulty satellites with the largest residuals and exceeding the threshold are cyclically eliminated through a forward search method. Finally, the M-LS filtering method is adopted to perform the GNSS / INS tight combination solution again to obtain the correct positioning result. The present invention fully exploits the role of INS navigation information in auxiliary anti-error, greatly enhances the anti-gross error capability of the combined navigation system, and has high engineering application value. BRIEF DESCRIPTION OF THE DRAWINGS

[0062] Figure 1 This is the execution flow of the INS-assisted GNSS pseudorange gross error detection method provided by a specific embodiment of the present invention.

[0063] Figure 2 It is the experimental verification platform of the present invention.

[0064] Figure 3-Figure 5 are the experimental verification results, which are the positioning errors in the E, N, and U directions respectively. The dotted line plus the hollow circle marker represents the positioning error of the RAIM algorithm, and the solid line represents the positioning error of this algorithm.

[0065] Figure 6 It is the display of the experimental positioning trajectory on the map.

[0066] Figure 7 It is a comparison diagram of experimental positioning trajectories. The dotted line represents the positioning trajectory of the RAIM algorithm, and the solid line represents the positioning trajectory of this algorithm. DETAILED DESCRIPTION

[0067] The present invention will be further explained below with reference to the accompanying drawings and specific embodiments. It should be understood that the following specific embodiments are only used to illustrate the present invention and are not used to limit the scope of the present invention.

[0068] The present invention provides an INS-assisted GNSS pseudorange gross error detection method, specifically as follows:

[0069] (1) Tight combination of GNSS single-point pseudorange observations and Doppler observations with INS

[0070] The INS position, velocity, and attitude error equations based on the psi-angle (ψ) error model are:

[0071]

[0072]

[0073]

[0074] In the formula, the superscript n represents the navigation coordinate system, the subscripts i, e, and b represent the inertial coordinate system, the earth coordinate system, and the carrier coordinate system, respectively. n represents the position vector; represents the transfer rate; v n represents the carrier speed; ψ represents the platform misalignment angle; f n Indicates the relative force under the navigation system; represents the Earth's rotation rate; Represents the gravity calculated by the model; represents the transformation matrix from the carrier system to the navigation system; f b Indicates the specific force under the load system; It represents the angular velocity of the navigation system relative to the inertial system; It represents the angular velocity of the carrier system relative to the inertial system.

[0075] The construction of the tightly combined filter state vector is determined by the error state of the inertial device and the error state of the GNSS:

[0076]

[0077] Where, the superscripts n and b represent the navigation coordinate system and the carrier coordinate system respectively. n Indicates the misalignment angle of the three-dimensional platform under the navigation system; δv n represents the three-dimensional velocity error in the navigation system; δr n Represents the three-dimensional position error in the navigation system; ε b Indicates the three-dimensional gyro drift under the load system; Indicates the zero bias of the three-dimensional accelerometer under the load system; dt G Indicates the GNSS receiver clock offset; dt d Indicates the GNSS receiver clock drift.

[0078] The system error dynamics equation is shown below:

[0079]

[0080] Where F represents the state transfer matrix, which can be obtained by equations [1]-[3]; G represents the noise distribution matrix; and w represents the system noise.

[0081] Single-point pseudorange observations and Doppler observations are tightly combined with INS, where the ionospheric delay and tropospheric delay of the GNSS observations are corrected using the Klobuchar and Saastamoinen models, respectively. The measurement model is expressed as follows:

[0082] z k =H k x k +v k [6]

[0083] Where z k represents the measurement vector; subscript k represents the kth epoch; H k represents the measurement matrix; x k represents the state vector; v k represents the measurement noise, which obeys the zero-mean Gaussian distribution.

[0084] The measurement vector is the original GNSS observation value m GNSS and INS predicted value The difference between them is shown below:

[0085]

[0086]

[0087] Where, P s,f 、 They represent the original pseudorange observation value and Doppler observation value respectively; the superscript f represents the frequency; the superscript s represents the satellite system, including GPS, BDS, and Galileo; ρ INS 、 Represents INS predicted pseudorange and predicted pseudorange rate respectively; △dtr P 、 are the sum of the error corrections associated with the pseudorange and Doppler observations from the receiver clock; and are the sum of other error corrections for the pseudorange and Doppler observations.

[0088] The measurement model is expressed as follows:

[0089]

[0090]

[0091]

[0092]

[0093]

[0094] I=[1 1 ... 1] T

[14]

[0095] Where, e s,f represents the direction cosines from the receiver to the satellite; Represents the position error conversion matrix from the navigation coordinate system to the earth coordinate system; Represents the velocity error conversion matrix from the navigation coordinate system to the earth coordinate system.

[0096] The sliding window detection of the real-time estimation of the receiver clock error and clock drift is to calculate the corresponding sample mean and sample mean square error of the real-time estimation value within the sliding window. The calculation formula is as follows:

[0097]

[0098]

[0099] Where, represents the sample mean, and S represents the sample mean square error.

[0100] (2) INS-assisted GNSS residual test

[0101] The position parameters calculated by the inertial navigation update are X = (x, y, z) in the Earth-centered Earth-fixed coordinate system. Substitute it into the following pseudorange observation equation:

[0102]

[0103] Where, the superscript represents the i-th satellite and the subscript r represents the receiver; represents the pseudorange observation value; represents the geometric distance between the satellite and the earth; c represents the speed of light; δt r represents the receiver clock error; δt i represents the satellite clock error; represents the equivalent tropospheric delay; represents the equivalent ionospheric delay; represents the equivalent orbit error; Represents random error. Among them, the error terms such as tropospheric and ionospheric delays and satellite clock errors are corrected using the corresponding models, leaving only one unknown parameter in the equation, the receiver clock error. The above equation is adjusted to obtain δt r The estimates and the corresponding residuals of each equation are as follows:

[0104]

[0105] Where, the superscript represents the i-th satellite; v i represents the residual; l i represents the error of the i-th satellite;

[0106] (3) Forward search cycle to eliminate faulty satellites

[0107] Since sliding window monitoring is adopted, if the sample mean and sample mean square error of the sliding window monitoring exceed the set threshold, a residual test is performed and a forward search cycle is started to eliminate the faulty satellite.

[0108] As mentioned above, it can be seen from formula

[18] that the residual v of the observation value corresponding to each satellite is i It will be mainly affected by its own error; if there is no gross error for the i-th satellite, but there is a gross error for the j-th satellite, then the corresponding residual v of the j-th satellite is j The impact is The impact on the i-th satellite is only When there are a large number of satellites, there is a strong reason to believe that the observation value corresponding to the satellite with the largest residual that exceeds the threshold is the satellite with a gross error. At this time, this satellite is eliminated, and the observation equation shown in formula

[17] is reconstructed but does not include the jth satellite. This is repeated until the residual meets the requirements. When the number of satellites is small, the sample mean of the sliding window monitoring of the previous output epoch is also substituted into formula

[17] for residual detection, and the satellite with the largest residual is cyclically eliminated until the residual meets the requirements.

[0109] (4) GNSS / INS tight integration solution based on M-LS filtering

[0110] Assuming that the components of the observation vector of the epoch are independent of each other, but may contain abnormal errors and obey the contaminated normal distribution, while the state vector predicted by the dynamic model still obeys the normal distribution, the observation vector is estimated using the robust M method, and the state parameters are still estimated using the least squares (LS) method.

[0111] The M-LS filter extreme value condition is defined using the equivalent weight matrix, and the recursive solution is obtained as follows:

[0112]

[0113] Where K MLS It can still be called the gain matrix, and the expression is:

[0114]

[0115] Where, Represents the equivalent weight matrix of the observation vector, using the IGGIII weight function:

[0116]

[0117] Where, The weight matrix R k The i-th and j-th elements of γ ij for:

[0118]

[0119] in:

[0120]

[0121] Where, represents the normalized innovation, and k0 and k1 are the corresponding thresholds. This can produce a more accurate robust positioning result after eliminating GNSS gross errors.

[0122] The following experiments are conducted in a real environment to verify the effect and accuracy of the technical solution of the present invention, and compare the trajectory accuracy obtained by the improved algorithm with that of the conventional algorithm. The experimental site is the Changsha Internet of Vehicles Test Center, and the ADIS16488IMU and UM482GNSS board are used to collect data for IE post-calculation, which is used as the true value to evaluate the accuracy. The test results are as follows Figure 3-Figure 5 , as shown in Table 1.

[0123] Table 1 Comparison of RMSE (m) in E, N, and U directions between the RAIM algorithm and this algorithm

[0124]

[0125] Depend on Figure 3-Figure 5 It can be seen that in actual test experiments, the RAIM algorithm still has large positioning errors in the positioning results of some epochs. This algorithm performs better in the case of GNSS pseudorange gross errors. As shown in Table 1, the RMSE of this algorithm in the E and U directions is significantly improved compared to the RAIM algorithm, by approximately 68.1% and 21.5%, respectively. Through testing, it can be seen that the INS-assisted GNSS pseudorange gross error detection method proposed in this invention has significant advantages in positioning accuracy, continuity, and robustness in complex outdoor environments that may cause GNSS pseudorange gross errors.

[0126] It should be noted that the above content merely illustrates the technical idea of ​​the present invention and cannot be used to limit the scope of protection of the present invention. For ordinary technicians in this technical field, several improvements and modifications can be made without departing from the principles of the present invention. These improvements and modifications all fall within the scope of protection of the claims of the present invention.

Claims

1. An INS-assisted GNSS pseudorange gross error detection method, characterized by: The following steps are involved: (1) Tight combination of GNSS single-point pseudorange observations and Doppler observations with INS The INS position, velocity, and attitude error equations based on the psi-angle (ψ) error model are: Where, the superscript n represents the navigation coordinate system, and the subscripts i, e, and b represent the inertial coordinate system, the earth coordinate system, and the carrier coordinate system, respectively; r n represents the position vector; represents the transfer rate; v n represents the carrier speed; ψ represents the platform misalignment angle; f n Indicates the relative force under the navigation system; represents the Earth's rotation rate; Represents the gravity calculated by the model; represents the transformation matrix from the carrier system to the navigation system; f b Indicates the specific force under the load system; It represents the angular velocity of the navigation system relative to the inertial system; It represents the angular velocity of the carrier system relative to the inertial system; The construction of the tightly combined filter state vector is determined by the error state of the inertial device and the error state of the GNSS: Where, the superscripts n and b represent the navigation coordinate system and the carrier coordinate system respectively; ψ n Indicates the misalignment angle of the three-dimensional platform under the navigation system; δv n represents the three-dimensional velocity error in the navigation system; δr n Represents the three-dimensional position error in the navigation system; ε b Indicates the three-dimensional gyro drift under the load system; Indicates the zero bias of the three-dimensional accelerometer under the load system; dt G Indicates the GNSS receiver clock offset; dt d Indicates the GNSS receiver clock drift; The system error dynamics equation is shown below: Where, F represents the state transfer matrix, which is obtained by equations [1]-[3]; G represents the noise distribution matrix; w represents the system noise; Single-point pseudorange observations and Doppler observations are tightly combined with INS. The ionospheric delay and tropospheric delay of the GNSS observations are corrected using the Klobuchar and Saastamoinen models, respectively. The measurement model is expressed as follows: z k =H k x k +v k [6] Where z k represents the measurement vector; subscript k represents the kth epoch; H k represents the measurement matrix; x k represents the state vector; v k represents the measurement noise, which obeys the zero-mean Gaussian distribution; The measurement vector is the original GNSS observation value m GNSS and INS predicted value The difference between ; as shown below: Where, P s,f 、 They represent the original pseudorange observation value and Doppler observation value respectively; the superscript f represents the frequency; the superscript s represents the satellite system, including GPS, BDS, and Galileo; ρ INS 、 They represent the INS predicted pseudorange and predicted pseudorange rate respectively; Δdtr P 、 are the sum of the error corrections associated with the pseudorange and Doppler observations from the receiver clock; Δδ P 、 denote the sum of other error corrections of pseudorange and Doppler observations respectively; The measurement model is expressed as follows: I＝[11...1] T [14] Where, e s,f represents the direction cosines from the receiver to the satellite; Represents the position error conversion matrix from the navigation coordinate system to the earth coordinate system; Represents the velocity error conversion matrix from the navigation coordinate system to the earth coordinate system; The sliding window detection of the real-time estimation of the receiver clock error and clock drift is to calculate the corresponding sample mean and sample mean square error of the real-time estimation value within the sliding window. The calculation formula is as follows: Where, represents the sample mean, S represents the sample mean square error; (2) INS-assisted GNSS residual test The position parameters calculated by the inertial navigation update are X = (x, y, z) in the Earth-centered Earth-fixed coordinate system. Substitute it into the following pseudorange observation equation: Where, the superscript represents the i-th satellite and the subscript r represents the receiver; represents the pseudorange observation value; represents the geometric distance between the satellite and the earth; c represents the speed of light; δt r represents the receiver clock error; δt i represents the satellite clock error; represents the equivalent tropospheric delay; represents the equivalent ionospheric delay; represents the equivalent orbit error; represents random error; among them, the tropospheric and ionospheric delays and satellite clock error errors are corrected using the corresponding models, leaving only one unknown parameter in the equation, the receiver clock error. The above equation is adjusted to obtain δt r The estimates and the corresponding residuals of each equation are as follows: Where, the superscript represents the i-th satellite; v i represents the residual; l i represents the error of the i-th satellite; (3) Forward search cycle to eliminate faulty satellites Since sliding window monitoring is adopted, if the sample mean and sample mean square error of the sliding window monitoring exceed the set threshold, a residual test is performed and a forward search cycle is started to eliminate the faulty satellite; As mentioned above, it can be seen from formula [18] that the residual v of the observation value corresponding to each satellite is i It will be mainly affected by its own error; if there is no gross error for the i-th satellite, but there is a gross error for the j-th satellite, then the corresponding residual v of the j-th satellite is j The impact is The impact on the i-th satellite is only When there are a large number of satellites, there is a strong reason to believe that the observation value corresponding to the satellite with the largest residual that exceeds the threshold is the satellite with a gross error. At this time, this satellite is eliminated, and the observation equation shown in formula [17] is reconstructed but does not include the jth satellite. This is repeated until the residual meets the requirements. When the number of satellites is small, the sample mean of the sliding window monitoring of the previous output epoch is also substituted into formula [17] for residual detection, and the satellite with the largest residual is cyclically eliminated until the residual meets the requirements. (4) GNSS / INS tight integration solution based on M-LS filtering Assuming that the components of the observation vector of the epoch are independent of each other, but may contain abnormal errors and obey the contaminated normal distribution, while the state vector predicted by the dynamic model still obeys the normal distribution, the observation vector is estimated using the robust M method, and the state parameters are still estimated using the least squares (LS) method. The M-LS filter extreme value condition is defined using the equivalent weight matrix, and the recursive solution is obtained as follows: Where K MLS Still called the gain matrix, the expression is: Where, Represents the equivalent weight matrix of the observation vector, using the IGGIII weight function: Where, is the weight matrix R k The i-th and j-th elements of γ ij for: in: Where, represents the standardized new information, k0 and k1 are the corresponding thresholds set; thus, a more accurate robust positioning result is obtained after the GNSS gross error is eliminated.

2. The INS-assisted GNSS pseudorange gross error detection method according to claim 1, characterized in that: In step (1), the GNSS single-point pseudorange observation value and Doppler observation value are tightly combined with the INS. The pseudorange observation value and Doppler observation value are used, and the three systems of GPS, BDS and Galieo are used to model and estimate the receiver clock error and receiver clock drift respectively. The sliding window monitoring method is used to estimate the receiver clock error and receiver clock drift in real time to solve the sample variance and sample mean square error.

3. The INS-assisted GNSS pseudorange gross error detection method according to claim 1, characterized in that: In the INS-assisted GNSS residual test described in step (2), the carrier position calculated by the INS update is first substituted into the GNSS observation equation, and then the corresponding error terms are corrected using the model. Finally, the remaining unknown parameters are adjusted and the residual test is performed.

4. The INS-assisted GNSS pseudorange gross error detection method according to claim 1, characterized in that: In step (3), the forward search cycle is used to eliminate faulty satellites, and residual tests are performed on epochs that exceed the sliding window monitoring threshold and a forward search cycle is started to eliminate faulty satellites. When the number of satellites is less than 8, the sample mean of the sliding window monitoring is substituted into the GNSS observation equation and a forward search cycle is performed to eliminate faulty satellites.

5. The INS-assisted GNSS pseudorange gross error detection method according to claim 1, characterized in that: The GNSS / INS tight combination solution based on M-LS filtering described in step (4) uses robust M estimation for the observation vector after removing the faulty satellite in step (3), and still uses least squares (LS) estimation for the state parameters; by setting the corresponding threshold, the noise matrix is ​​replaced by the equivalent weight matrix of the observation vector to obtain the robust Kalman filter gain, and finally obtain a more accurate positioning result.

Images