Beidou robust Kalman filtering positioning method based on non-Gaussian noise model

Abstract

The invention relates to the technical field of satellite navigation positioning, in particular to a Beidou robust Kalman filtering positioning method based on a non-Gaussian noise model. According to the method, non-Gaussian noise generated by shielding and a multipath effect in Beidou positioning is accurately represented by constructing a Gaussian mixture model and a student t distribution model, model parameters are updated online by an expectation maximization algorithm, and the problem that a traditional Gaussian hypothesis model cannot adapt to noise statistical characteristics is solved; the non-Gaussian noise model is fused with robust filtering, a self-adaptive robust weight matrix is constructed through a Huber cost function, an IGGIII scheme or noise probability density, an observation covariance matrix is dynamically adjusted, the influence of abnormal observation values is suppressed, the defects that classical robust filtering lacks noise modeling and weight adjustment fails under continuous strong interference are overcome, and the robustness of the robust filtering is improved. And meanwhile, the problem of complex particle filtering calculation is avoided.

Description

Technical Field

[0001] This invention relates to the field of satellite navigation and positioning technology, specifically to a robust Kalman filter positioning method for BeiDou based on a non-Gaussian noise model. Background Technology

[0002] The BeiDou Navigation Satellite System has significant application value in the health monitoring of major infrastructure, especially for the safety assessment of critical structures such as highway slopes and large bridges. These types of infrastructure undergo minute deformations during long-term operation due to environmental loads and geological activities, requiring millimeter-level displacement monitoring through high-precision positioning technology. While satellite navigation-based structural monitoring offers advantages over traditional measurement methods, such as all-weather operation and automation, it also faces measurement challenges in special environments.

[0003] In practical engineering applications, highway slopes and the areas surrounding bridges often present complex electromagnetic environments. The undulating terrain of slopes obstructs satellite signals, and the main bridge structure and surrounding buildings introduce significant multipath effects, causing observation noise to deviate from the conventional Gaussian distribution, exhibiting spikes, impulses, and heavy tails. This situation severely impacts the positioning accuracy of Kalman filtering algorithms and may even lead to erroneous deformation warnings. Existing filtering methods based on the Gaussian assumption struggle to distinguish between actual displacement and noise interference, failing to meet the stringent data reliability requirements of infrastructure safety monitoring. Summary of the Invention

[0004] To address the shortcomings of existing technologies, this invention provides a robust Kalman filter positioning method for BeiDou based on a non-Gaussian noise model. This method solves the problems of insufficient positioning accuracy, erroneous deformation warnings, and easy filter divergence caused by the non-Gaussian nature of observation noise in existing technologies.

[0005] To achieve the above objectives, the present invention provides the following technical solution: a robust Kalman filter positioning method for BeiDou based on a non-Gaussian noise model, comprising the following steps:

[0006] S01: Obtain the raw observation data output by the BeiDou Navigation Satellite System. The raw observation data includes at least one of pseudorange observations, carrier phase observations, and Doppler frequency shift observations. The raw observation data is used for subsequent positioning calculations.

[0007] S02: Establish a non-Gaussian noise model. For the observation noise in the original observation data and the process noise in the BeiDou positioning system, a noise model is preset to show that the observation noise and the process noise follow a non-Gaussian distribution in order to accurately characterize their statistical characteristics. The non-Gaussian noise model can adaptively update its parameters online.

[0008] S03: The non-Gaussian noise model and robust filtering are fused together, the fusion including:

[0009] S031: Design the state equations and observation equations, wherein the state equations are expressed as follows: The observation equation is expressed as ,in, For the calendar The observation vector, Let the state vector be... Here is the state transition matrix. For the observation matrix, Let be the observed noise vector. The process noise vector is the state vector used to describe the dynamic parameters of the BeiDou receiver;

[0010] S032: Based on the non-Gaussian noise model established in step S02, construct an adaptive robust weight matrix. The adaptive robust weight matrix is ​​used to dynamically adjust the covariance matrix of the observed values ​​during the filtering process to suppress the influence of abnormal observed values.

[0011] S04: Iteratively execute the robust filtering, the iterative execution including:

[0012] S041: Time update: Based on the state equation, use the optimal state estimate of the previous epoch to predict the prior estimate of the state vector of the current epoch, and calculate the corresponding prediction error covariance matrix. The prediction error covariance matrix reflects the uncertainty of the state estimate.

[0013] S042: Measurement update, including:

[0014] (a) Calculate the innovation sequence, which is the residual between the observation vector and the prior estimate of the state vector, and evaluate the reliability of the observation based on the non-Gaussian noise model.

[0015] (b) Based on the adaptive robust weight matrix, abnormal observations in the innovation sequence that exceed a preset threshold are weighted to reduce their impact on the positioning results;

[0016] (c) Update the posterior estimate of the state vector and the error covariance matrix, wherein the posterior estimate is the state estimate after correction of the observations;

[0017] (d) Output the posterior estimate of the state vector to obtain a high-precision positioning result.

[0018] Further, in step S02, the non-Gaussian noise model is a Gaussian mixture model, and the Gaussian mixture model affects the observed noise. probability density function Defined as:

[0019] ;

[0020] in, The number of mixture components in the Gaussian mixture model. For the first The weighting coefficients of each mixed component, and The first The mean vector and covariance matrix of each Gaussian component;

[0021] The single Gaussian component The expression is:

[0022] ;

[0023] in, Let be the dimension of the observed noise.

[0024] Further, in step S02, the non-Gaussian noise model is a Student's t-distribution model, and the Student's t-distribution model affects the observed noise. probability density function Defined as:

[0025] ;

[0026] in, Let be the dimension of the observed noise. Let be the degrees of freedom of the student t-distribution. Let be the mean vector of the student's t-distribution. Let be the scaling matrix of the student t-distribution. The student t-distribution model is used to describe observation noise with heavy-tailed characteristics, and is a gamma function.

[0027] Further, in step S02, the parameters of the Gaussian mixture model are updated online using an expectation-maximization algorithm, the update equation of which includes:

[0028] ;

[0029] in, To adjust the sliding window size, Represents the epoch Observations Belongs to the The posterior probabilities of the mixture components are used to implement adaptive parameter optimization of the Gaussian mixture model using the expectation-maximization algorithm.

[0030] Furthermore, in step S02, the non-Gaussian noise model further includes a heavy-tailed distribution correction term, which is corrected by introducing a Cauchy distribution. The corrected observation noise probability density function... Defined as:

[0031] ;

[0032] in, This is the main body of Gaussian. For Cauchy's corrections;

[0033] The Cauchy distribution density function The expression is:

[0034] ;

[0035] in, For the noise dimension, The scaling matrix, As a mixing scaling factor, the heavy-tailed distribution correction term is intended to enhance the non-Gaussian noise model's ability to characterize impulse noise.

[0036] Further, in step S032, the adaptive robust weight matrix is ​​constructed using the Huber cost function. Defined as:

[0037] ;

[0038] in, For the new interest, To adjust the parameters, the Huber cost function classifies the information, applies a secondary penalty to errors smaller than the adjustment parameters, and applies a linear penalty to errors larger than the adjustment parameters, in order to reduce the impact of outlier observations.

[0039] Further, in step S032, the adaptive robust weight matrix is ​​constructed using the IGGIII scheme to dynamically adjust the observation weights, and the weight function of the IGGIII scheme is... Defined as:

[0040] ;

[0041] in, To standardize the residuals, The first threshold, The second threshold is used as the weighting function to reduce the weight of the observed values ​​in segments according to the magnitude of the standardized residual. When the standardized residual exceeds the second threshold, the weight is zero, thereby realizing flexible handling of outliers of different degrees.

[0042] Further, in step S032, the adaptive robust weight matrix The construction method is as follows:

[0043] ;

[0044] in, The number of observations, For the first The weight of each observation, For the first The probability density values ​​of each observation in the non-Gaussian noise model. To prevent the removal of zero constants, the weight matrix evaluates the reliability of the observations in real time using the probability density of the non-Gaussian noise model, and adaptively weights the observation data accordingly to enhance the robustness of the filter.

[0045] Further, in step S031, the state vector The state vector, comprising at least four of the following: three-dimensional position, three-dimensional velocity, receiver clock bias, and receiver clock drift, is used to comprehensively describe the dynamic motion state of the BeiDou receiver in space, as well as its clock bias and drift; the raw observation data The observations include pseudorange and carrier phase observations from multiple frequency bands such as BeiDou B1I, B2I, and B3I. These multi-frequency band observations are used for high-precision state updates through the observation equations.

[0046] Furthermore, in step S042(b), when the amplitude of the new information is detected to exceed the preset threshold, the robustness mechanism is activated. The preset threshold is dynamically determined based on the satellite elevation angle, carrier-to-noise ratio, and historical noise statistical characteristics to adapt to the complex and ever-changing actual positioning environment. In addition, the weighting process uses linear weighting for observations with slight deviations and completely removes observations with significant outliers. The weighting process aims to retain effective information while blocking the error propagation path to ensure the reliability of the positioning results.

[0047] Compared with the prior art, the beneficial effects of the present invention are as follows:

[0048] This invention accurately characterizes non-Gaussian noise generated by occlusion and multipath effects in BeiDou positioning by constructing a Gaussian mixture model and a Student's t-distribution model. It then uses an expectation-maximization algorithm to update model parameters online, addressing the problem that traditional Gaussian assumption models cannot adapt to the statistical characteristics of noise. Furthermore, this non-Gaussian noise model is fused with robust filtering. An adaptive robust weight matrix is ​​constructed using the Huber cost function, the IGGIII scheme, or noise probability density to dynamically adjust the observation covariance matrix and suppress the influence of outlier observations. This overcomes the shortcomings of classical robust filtering, such as lack of noise modeling and failure of weight adjustment under continuous strong interference, while avoiding the computational complexity of particle filtering. Finally, through iterative execution of robust filtering, the invention effectively solves the problems of insufficient BeiDou positioning accuracy and easy filter divergence in complex environments, improving positioning robustness and real-time performance, and meeting the positioning reliability requirements of infrastructure monitoring such as highway slopes and bridges. Attached Figure Description

[0049] Figure 1 This is a flowchart of the method of the present invention. Detailed Implementation

[0050] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0051] Please see Figure 1 This invention provides a robust Kalman filter positioning method for BeiDou based on a non-Gaussian noise model, comprising the following steps:

[0052] S01: Obtain the raw observation data output by the BeiDou Navigation Satellite System. The raw observation data includes at least one of pseudorange observations, carrier phase observations, and Doppler frequency shift observations. The raw observation data is used for subsequent positioning calculations.

[0053] S02: Establish a non-Gaussian noise model. For the observation noise in the original observation data and the process noise in the BeiDou positioning system, a noise model is preset to show that the observation noise and the process noise follow a non-Gaussian distribution in order to accurately characterize their statistical characteristics. The non-Gaussian noise model can adaptively update its parameters online.

[0054] S03: The non-Gaussian noise model and robust filtering are fused together, the fusion including:

[0055] S031: Design the state equations and observation equations, wherein the state equations are expressed as follows: The observation equation is expressed as ,in, For the calendar The observation vector, Let the state vector be... Here is the state transition matrix. For the observation matrix, Let be the observed noise vector. The process noise vector is the state vector used to describe the dynamic parameters of the BeiDou receiver;

[0056] S032: Based on the non-Gaussian noise model established in step S02, construct an adaptive robust weight matrix. The adaptive robust weight matrix is ​​used to dynamically adjust the covariance matrix of the observed values ​​during the filtering process to suppress the influence of abnormal observed values.

[0057] S04: Iteratively execute the robust filtering, the iterative execution including:

[0058] S041: Time update: Based on the state equation, use the optimal state estimate of the previous epoch to predict the prior estimate of the state vector of the current epoch, and calculate the corresponding prediction error covariance matrix. The prediction error covariance matrix reflects the uncertainty of the state estimate.

[0059] S042: Measurement update, including:

[0060] (a) Calculate the innovation sequence, which is the residual between the observation vector and the prior estimate of the state vector, and evaluate the reliability of the observation based on the non-Gaussian noise model.

[0061] (b) Based on the adaptive robust weight matrix, abnormal observations in the innovation sequence that exceed a preset threshold are weighted to reduce their impact on the positioning results;

[0062] (c) Update the posterior estimate of the state vector and the error covariance matrix, wherein the posterior estimate is the state estimate after correction of the observations;

[0063] (d) Output the posterior estimate of the state vector to obtain a high-precision positioning result.

[0064] Specifically, this implementation method is applied to the scenario of highway slope deformation monitoring to solve the problem of insufficient positioning accuracy caused by the obstruction of BeiDou positioning and the non-Gaussian observation noise due to the multipath effect.

[0065] Step S01: Obtain raw observation data

[0066] A high-precision receiver supporting BeiDou multi-frequency observation (such as a certain model of BeiDou dual-frequency receiver) was deployed at a highway slope monitoring point. The sampling interval was set to 1Hz (to adapt to the slow deformation characteristics of the slope). Raw observation data output by the BeiDou satellite navigation system was collected, including pseudorange observations, carrier phase observations, and Doppler shift observations from the BeiDou B1I, B2I, and B3I frequency bands. The pseudorange observations reflect the distance observation information from the satellite to the receiver, the carrier phase observations provide more precise distance change information, and the Doppler shift observations help characterize the receiver's dynamic motion state, providing multi-dimensional data support for subsequent positioning calculations.

[0067] Step S02: Establish a non-Gaussian noise model

[0068] To address the observation noise in the raw observation data and the process noise from the BeiDou positioning system, a non-Gaussian noise model is pre-defined, where both observation and process noise follow a Gaussian mixture distribution (the Gaussian mixture model can accurately characterize the multi-peak and heavy-tailed characteristics of noise caused by multipath effects). Furthermore, the model parameters support online adaptive updates. The probability density function of the Gaussian mixture model is:

[0069] ;

[0070] in, For the calendar The observed noise vector, with dimensions in meters; This represents the number of mixed components, for example, taking 3 (corresponding to the three main sources of noise under multipath effects). For the first The weighting coefficients of the mixed components satisfy the following conditions: Dimensionless; For the first The mean vector of Gaussian components, with dimensions equal to... Consistent, the initial value can be obtained through statistical analysis of historical cleanliness observation data; For the first The covariance matrix of Gaussian components, with the dimension of square meters, is calculated based on historical data. The probability density function of a single Gaussian component is expressed as:

[0071] ;

[0072] in, To determine the dimension of observation noise, for example, when using B1I and B2I pseudorange observations simultaneously. Dimensionless; Covariance matrix Determinant, dimensions and of The power-law consistency is achieved. This model can accurately characterize the non-Gaussian statistical properties of observation noise, laying the foundation for subsequent robust filtering.

[0073] Step S03: Fusion of non-Gaussian noise model and robust filtering

[0074] Step S031: Design state equations and observation equations

[0075] State vector The dynamic parameters used to describe the BeiDou receiver include three-dimensional position ( (dimensions in meters), three-dimensional velocity ( The unit is meters per second), receiver clock error ( The unit of measurement is meters, which has been converted to distance via the speed of light; receiver clock drift ( (The unit is meters per second, also converted from the speed of light), totaling 8 dimensions, i.e. .

[0076] The state equation is constructed based on the receiver's uniform motion model, and its expression is:

[0077] ;

[0078] in, for The dimensionless state transition matrix has the following form:

[0079] ;

[0080] It is a 3-order identity matrix. It is a 3rd order zero matrix. The sampling interval is 1 second (i.e., 1 second corresponds to 1 Hz). It is a 1-dimensional zero vector; The process noise vector follows a Gaussian mixture model established in step S02, with dimensions equal to... Consistency is used to characterize receiver motion uncertainty and satellite orbital errors, etc.

[0081] The observation equation is constructed based on pseudorange observations, and its expression is:

[0082] ;

[0083] in, For the calendar The observation vector, for example, takes the pseudorange observations of B1I and B2I, with the dimension of meters; for A dimensionless observation matrix whose elements are determined by the geometric relationship between pseudorange observations and state vectors, for example, corresponding to... The elements are for satellites in The unit vector component of the direction, corresponding to The element is 1; The observation noise vector is consistent with the observation noise in step S02.

[0084] Step S032: Construct the adaptive robust weight matrix

[0085] Based on the Gaussian mixture model in step S02, an adaptive robust weight matrix is ​​constructed. This is used to dynamically adjust the covariance matrix of the observations. The weight matrix is ​​a diagonal matrix, expressed as:

[0086] ;

[0087] in, For example, the number of observations. (Corresponding to B1I and B2I pseudoranges); For the first The weight of each observation is calculated using the following formula:

[0088] ;

[0089] For the first The probability density value of the observed noise under the Gaussian mixture model is dimensionless. To prevent the division by zero constant, it is dimensionless. When the observed value is abnormal, Smaller This reduces the weighting of outlier observations, thus achieving dynamic deweighting.

[0090] Step S04: Iteratively perform robust filtering

[0091] Step S041: Time Update

[0092] Based on the state equation, using the previous epoch ( Optimal state estimation And error covariance matrix Predict the current epoch ( The prior estimate of the state vector and prediction error covariance matrix The formulas are as follows:

[0093] ;

[0094] ;

[0095] in, Let be the error covariance matrix of the state estimate in the previous epoch, with dimensions equal to the squares of the corresponding state parameters; Let be the process noise covariance matrix, with dimensions and Consistent, its value is determined based on the parameters of the Gaussian mixture model; This is the transpose of the state transition matrix. This reflects the uncertainty of the prior state estimate, providing a basis for subsequent measurement updates.

[0096] Step S042: Measurement Update

[0097] o(a) Calculate the new information sequence The formula is:

[0098] ;

[0099] in, For the observation vector, For the observed and predicted values, both are measured in meters, so they are dimensionless and do not require normalization. Combining the Gaussian mixture model from step S02, each innovation is calculated. ( The The probability density corresponding to (number of elements) ,like If the probability is less than a preset probability threshold (e.g., 0.05), the corresponding observation is considered to have low reliability.

[0100] (b) Adaptive robust weight matrix based on step S032 Observations with low reliability are downweighted. For example, when an observation corresponds to... If an observation is considered an anomaly, its contribution to the observation covariance matrix is ​​reduced using a weighting matrix to decrease its interference with the positioning results.

[0101] (c) Update the posterior estimate of the state vector With error covariance matrix The formulas are as follows:

[0102] ;

[0103] ;

[0104] ;

[0105] in, The Kalman gain matrix is ​​dimensionless. The original observation noise covariance matrix has the dimension of meters². It is an identity matrix, dimensionless; It is the inverse of the weight matrix, which adjusts the effective contribution of the observation covariance matrix through weight adjustment.

[0106] (d) Posterior estimate of the output state vector The three-dimensional location information is extracted from this data to obtain the high-precision positioning results of the highway slope monitoring points. This result can effectively suppress the influence of non-Gaussian noise, improve the accuracy of slope deformation monitoring, and avoid erroneous deformation warnings caused by noise interference.

[0107] In this embodiment, in step S02, the non-Gaussian noise model is a Gaussian mixture model, and the Gaussian mixture model affects the observed noise. probability density function Defined as:

[0108] ;

[0109] in, The number of mixture components in the Gaussian mixture model. For the first The weighting coefficients of each mixed component, and The first The mean vector and covariance matrix of each Gaussian component;

[0110] The single Gaussian component The expression is:

[0111] ;

[0112] in, Let be the dimension of the observed noise.

[0113] Specifically, in this embodiment, the non-Gaussian noise model adopts a Gaussian mixture model to accurately characterize the multi-peak and heavy-tailed characteristics of observation noise caused by the multipath effect in BeiDou positioning.

[0114] Gaussian mixture model for observation noise The probability density function is defined as:

[0115] ;

[0116] The parameters are explained as follows: The number of mixed components is determined based on the noise complexity of the actual monitoring scenario. For example, in highway slope monitoring, since noise mainly originates from three paths: direct satellite radiation, single reflection, and multiple reflections, the number of mixed components is taken as... ; For the first The weighting coefficients of the mixed components satisfy the following conditions: For example, initially set to , , These correspond to the noise contribution percentages of direct, single-reflection, and multiple-reflection paths, respectively. For the first The mean vector of the Gaussian components, with dimensions equal to the observation noise. Consistent (meters), for example, the mean noise level along the direct path. Meters, mean path noise of a single reflection Meters, mean of multiple reflection path noise rice; For the first The covariance matrix of Gaussian components, with dimensions in meters², for example... , , ( (This is the identity matrix), reflecting the degree of dispersion of noise from different paths.

[0117] Single Gaussian component The expression is:

[0118] ;

[0119] in, To observe the dimensions of noise, for example when using pseudorange observations in the BeiDou B1I and B2I frequency bands, ; Covariance matrix The determinant of the matrix, with dimensions in meters. ; for The transpose of both A dimensional vector, with units of meters, and (Dimension: meter) The product is dimensionless, ensuring that the exponent is dimensionless and that the formula has consistent dimensions.

[0120] This Gaussian mixture model can accurately fit the distribution characteristics of non-Gaussian observation noise, providing an accurate statistical basis for subsequent robust filtering and thus improving the stability of the positioning results.

[0121] In this embodiment, in step S02, the non-Gaussian noise model is a Student's t-distribution model, and the Student's t-distribution model affects the observed noise. probability density function Defined as:

[0122] ;

[0123] in, Let be the dimension of the observed noise. Let be the degrees of freedom of the student t-distribution. Let be the mean vector of the student's t-distribution. Let be the scaling matrix of the student t-distribution. The student t-distribution model is used to describe observation noise with heavy-tailed characteristics, and is a gamma function.

[0124] Specifically, in this embodiment, the non-Gaussian noise model adopts the Student's t-distribution model, which focuses on dealing with the impulse noise (heavy tail characteristic) caused by transient electromagnetic interference in BeiDou positioning, such as the impact of electromagnetic pulses generated by high-voltage lines around highways on the observation data.

[0125] Student t-distribution model for observation noise The probability density function is defined as:

[0126] ;

[0127] The parameters are explained as follows: To measure the dimension of noise, for example, when using only BeiDou B1I pseudorange observations, ; Let be the degrees of freedom of the student's t-distribution, whose value determines the degree of heavy tails in the noise distribution. For example, when frequent electromagnetic pulses exist, take . (The fewer the degrees of freedom, the more pronounced the heavy tail, making it more suitable for impulse noise.) Let be the mean vector of the student's t-distribution, with dimensions ∝ ... Consistency (meters), usually set to 0, represents the baseline noise level; Let be the scale matrix of the student t-distribution, with dimensions in meters², for example, calculated from historical observation data during periods of no disturbance. ( (Time is a scalar) The gamma function is used to normalize the probability density, for example... , .

[0128] Formula dimension check: The unit of measurement is meters. The unit of measurement is meter The product of the two is dimensionless, and (Dimensionless) together form the dimensionless term within the parentheses. The exponent is dimensionless, and the final dimension of the entire formula is meters. (The dimensions of probability density) are consistent.

[0129] The Student's t-distribution model can effectively characterize the heavy-tailed characteristics of impulse noise, avoid the underestimation of impulse noise by the traditional Gaussian model, thereby reducing the state estimation bias caused by impulse noise during the filtering process and improving the anti-interference capability of the positioning.

[0130] In this embodiment, in step S02, the parameters of the Gaussian mixture model are updated online using the expectation-maximization algorithm. The update equation of the expectation-maximization algorithm includes:

[0131] ;

[0132] in, To adjust the sliding window size, Represents the epoch Observations Belongs to the The posterior probabilities of the mixture components are used to implement adaptive parameter optimization of the Gaussian mixture model using the expectation-maximization algorithm.

[0133] Specifically, in this embodiment, the parameters of the Gaussian mixture model are updated online using the expectation-maximization (EM) algorithm to adapt to the dynamic changes in the noise characteristics of BeiDou positioning, such as the changes in multipath noise characteristics caused by changes in vegetation cover in different seasons during highway slope monitoring.

[0134] The update equation for the EM algorithm is as follows:

[0135] ;

[0136] The parameters are explained as follows: The sliding window size is determined based on the rate of noise change, for example, in a vegetation seasonal variation scenario. Era (corresponding to 50 seconds); For the calendar observation noise Belongs to the The posterior probability of each mixed component, dimensionless; For the updated number The weighting coefficients of the mixed components are dimensionless and satisfy the following conditions: ; For the updated number The mean vector of Gaussian components, with dimensions of... Consistent (meters); For the updated number A Gaussian component covariance matrix with dimensions in meters²; The number of mixed components, compared with the number in step S02 Consistency, for example .

[0137] The update process is as follows: First, in each epoch... Calculate based on current model parameters (E-step); then based on the sliding window... Individual metadata is updated sequentially. , , (M-step). For example, as vegetation cover increases, the proportion of noise contribution from multiple reflection paths rises. It will be gradually adjusted from 0.1 to 0.2. This will also increase accordingly to ensure that the model always adapts to the current noise characteristics.

[0138] By updating the Gaussian mixture model online using the EM algorithm, the model can dynamically track noise changes, avoiding model mismatch caused by parameter fixation and further improving the robustness of filtering.

[0139] In this embodiment, in step S02, the non-Gaussian noise model further includes a heavy-tailed distribution correction term, which is corrected by introducing a Cauchy distribution. The corrected observation noise probability density function... Defined as:

[0140] ;

[0141] in, This is the main body of Gaussian. For Cauchy's corrections;

[0142] The Cauchy distribution density function The expression is:

[0143] ;

[0144] in, For the noise dimension, The scaling matrix, As a mixing scaling factor, the heavy-tailed distribution correction term is intended to enhance the non-Gaussian noise model's ability to characterize impulse noise.

[0145] Specifically, in this embodiment, the non-Gaussian noise model introduces Cauchy distribution as a heavy-tailed distribution correction term to enhance its ability to characterize extreme impulse noise, such as strong impulse noise generated by sudden construction interference near highway slopes.

[0146] The corrected observation noise probability density function is:

[0147] ;

[0148] The parameters are explained as follows: This is a mixing proportionality coefficient used to adjust the contribution of the Cauchy correction, for example, during periods of frequent construction disturbances. Take during periods of no interference Dimensionless; This is the main body of Gaussian. The expression is consistent with step S02. Let Gaussian component covariance matrix have units of meters², for example... ( hour); This is the Cauchy correction part, used to characterize extreme impulse noise.

[0149] Cauchy distribution density function The expression is:

[0150] ;

[0151] in, For example, noise dimension (B1I, B2I pseudorange); Let be the scaling matrix of the Cauchy distribution, with dimensions . Consistent (meters²), for example The larger the scale, the wider the tail of the Cauchy distribution, and the better it can adapt to strong impulse noise; scale matrix The determinant of the matrix, with dimensions in meters. .

[0152] Formula dimension check: The unit of measurement is meters. The unit of measurement is meter , Dimensionless, the terms within parentheses are dimensionless, and the entire Cauchy distribution density function has a dimension of meters. The dimensions are consistent with those of the main Gaussian part, ensuring Dimensions are consistent.

[0153] By introducing a Cauchy distribution correction term, the non-Gaussian noise model can effectively capture extreme impulse noise, avoid filtering divergence caused by such noise, and improve the reliability of positioning results under sudden interference.

[0154] In this embodiment, in step S032, the adaptive robust weight matrix is ​​constructed using the Huber cost function. Defined as:

[0155] ;

[0156] in, For the new interest, To adjust the parameters, the Huber cost function classifies the information, applies a secondary penalty to errors smaller than the adjustment parameters, and applies a linear penalty to errors larger than the adjustment parameters, in order to reduce the impact of outlier observations.

[0157] Specifically, in this embodiment, the adaptive robust weight matrix is ​​constructed using the Huber cost function. By classifying the innovation, differentiated suppression of abnormal observations of different degrees is achieved. For example, in bridge monitoring, it can be used to deal with short-term observation anomalies caused by bridge deck vibration.

[0158] The Huber cost function is defined as:

[0159] ;

[0160] The parameters are explained as follows: The new information is the residual between the observed value and the prior state estimate, and its dimension is meter; The parameters are adjusted to classify error types, for example, based on bridge vibration amplitude statistics. Meters (when the information is less than this value, it is considered normal error; when it is greater than this value, it is considered abnormal error).

[0161] Weight calculation process: First, calculate the information corresponding to each observation. ( (The observation index is used); then the cost weight is calculated based on the Huber cost function. For example, when Rice time, (Normal weight), when Rice time, (Demotion); Finally, An adaptive robust weight matrix is ​​incorporated to achieve dynamic weight reduction of anomalous information.

[0162] Formula dimension check: and Both units are in meters. After addition, subtraction, multiplication and division, they become dimensionless, ensuring that the cost function calculation result is dimensionless and the weight values ​​are reasonable.

[0163] By using the Huber cost function, the impact of abnormal errors can be suppressed while retaining normal error information, avoiding the loss of accuracy caused by excessive weighting, and balancing positioning accuracy and robustness.

[0164] In this embodiment, in step S032, the construction of the adaptive robust weight matrix uses the IGGIII scheme to dynamically adjust the observation weights, and the weight function of the IGGIII scheme is... Defined as:

[0165] ;

[0166] in, To standardize the residuals, The first threshold, The second threshold is used as the weighting function to reduce the weight of the observed values ​​in segments according to the magnitude of the standardized residual. When the standardized residual exceeds the second threshold, the weight is zero, thereby realizing flexible handling of outliers of different degrees.

[0167] Specifically, in this embodiment, the IGGIII scheme is used to dynamically adjust the observation weights in order to flexibly handle abnormal observations of different degrees. For example, in highway slope monitoring, it distinguishes between observation anomalies caused by slight multipath interference and severe occlusion.

[0168] The weight function of the IGGIII scheme is defined as follows:

[0169]

[0170] The parameters are explained as follows: To standardize the residuals, through calculate( For the new interest, To observe the standard deviation of noise, all units are in meters; after standardization, they are dimensionless. The first threshold is used to determine minor anomalies, for example, taking... ; The second threshold is used to determine severe anomalies; for example, taking... .

[0171] Weight adjustment example: when (less than) )hour, (Normal weight); when (between) and (between) (Significantly reduced weight); when (greater than) )hour, (Completely removed).

[0172] Formula dimension check: For standardized residuals (dimensionless). , All are dimensionless thresholds to ensure that the weight calculation results are dimensionless and take values ​​between 0 and 1, which conforms to the weight definition.

[0173] The IGGIII scheme can progressively reduce the weight of slightly anomalous observations and completely remove severely anomalous observations, thus avoiding interference from outliers while preserving the most effective observation information and improving the filter's adaptability to complex environments.

[0174] In this embodiment, in step S032, the adaptive robust weight matrix The construction method is as follows:

[0175] ;

[0176] in, The number of observations, For the first The weight of each observation, For the first The probability density values ​​of each observation in the non-Gaussian noise model. To prevent the removal of zero constants, the weight matrix evaluates the reliability of the observations in real time using the probability density of the non-Gaussian noise model, and adaptively weights the observation data accordingly to enhance the robustness of the filter.

[0177] Specifically, in this embodiment, the adaptive robust weight matrix evaluates the reliability of the observations in real time through the probability density of the non-Gaussian noise model, thereby achieving more accurate weighting adjustments, such as in the monitoring scenario at the bottom of a bridge with complex multipath effects.

[0178] Adaptive robust weight matrix The construction method is as follows:

[0179]

[0180] The parameters are explained as follows: This represents the number of observations, for example, 3 (corresponding to BeiDou B1I, B2I, and B3I pseudoranges). For the first The weights of each observation are dimensionless. For the first The probability density value of an observation in a non-Gaussian noise model (such as a Gaussian mixture model) is dimensionless. To prevent division by zero constant, for example, take To avoid due to Too small a value will cause the weight calculation to be abnormal.

[0181] Example of weight calculation: In bridge bottom monitoring, a certain B1I pseudorange observation information probability density (High reliability) New information from a certain B2I pseudorange observation probability density (Low reliability) (Reduced weight).

[0182] Formula dimension check: and All are dimensionless, weights Dimensionless, but the dimensions are consistent throughout the calculation process.

[0183] This method directly correlates the reliability of observations under non-Gaussian noise models with weight adjustment, achieving a precise fit of "high reliability, high weight; low reliability, low weight," thereby further enhancing the anti-interference capability of filtering.

[0184] In this embodiment, in step S031, the state vector The state vector, comprising at least four of the following: three-dimensional position, three-dimensional velocity, receiver clock bias, and receiver clock drift, is used to comprehensively describe the dynamic motion state of the BeiDou receiver in space, as well as its clock bias and drift; the raw observation data The observations include pseudorange and carrier phase observations from multiple frequency bands such as BeiDou B1I, B2I, and B3I. These multi-frequency band observations are used for high-precision state updates through the observation equations.

[0185] Specifically, in this embodiment, the selection of state vectors and raw observation data aims to comprehensively describe receiver dynamics and support high-precision positioning, making it suitable for millimeter-level deformation monitoring of infrastructure such as highway slopes and large bridges.

[0186] State vector Includes three-dimensional position ( ), three-dimensional velocity ( ), receiver clock bias ( ), receiver clock drift ( There are a total of 8 parameters. Among them, the three-dimensional position is used to directly output the coordinates of the monitoring point, with the dimension of meters; the three-dimensional velocity is used to characterize the movement trend of the monitoring point, with the dimension of meters per second; and the receiver clock error is used to correct the deviation between the receiver clock and the satellite clock, with the dimension of meters (measured by the speed of light). Convert time deviation to distance deviation, i.e. , The time deviation is measured in seconds; receiver clock drift is used to correct the rate of change of clock deviation, measured in meters per second. , The change rate of time deviation is expressed in seconds per second. For example, in bridge deformation monitoring, the change in three-dimensional position directly reflects the settlement or displacement of the bridge, while the three-dimensional velocity can predict the deformation trend in advance.

[0187] The raw observation data includes pseudorange and carrier phase observations from multiple frequency bands of BeiDou B1I, B2I, and B3I. Pseudorange observations provide rapid coarse positioning results and are measured in meters; carrier phase observations offer higher precision (millimeter level) and are also measured in meters (based on carrier wavelength). Phase difference Converted to distance, i.e. For example, by using carrier phase observations of B1I (wavelength approximately 0.19 meters) and B2I (wavelength approximately 0.16 meters), ionospheric errors can be eliminated and positioning accuracy improved through dual-frequency combination.

[0188] By selecting the aforementioned state vectors and raw observation data, the dynamic parameters and high-precision observation information required for infrastructure monitoring can be fully covered, providing data support for millimeter-level positioning.

[0189] In this embodiment, in step S042(b), when the amplitude of the new information is detected to exceed the preset threshold, the robustness mechanism is activated. The preset threshold is dynamically determined based on the satellite elevation angle, carrier-to-noise ratio, and historical noise statistical characteristics to adapt to the complex and ever-changing actual positioning environment. Furthermore, the weighting process uses linear weighting for observations with slight deviations and completely removes significantly outlier observations. The weighting process aims to retain effective information while blocking error propagation paths, ensuring the reliability of the positioning results.

[0190] Specifically, in the implementation method, the activation and weighting of the robustness mechanism are based on dynamic thresholds and hierarchical strategies to adapt to complex and ever-changing positioning environments, such as in highway slope monitoring, to cope with changes in satellite elevation angle and carrier-to-noise ratio at different times (such as day and night).

[0191] The preset threshold is determined dynamically based on satellite elevation angle, carrier-to-noise ratio (CNR), and historical noise statistics. For example, the lower the satellite elevation angle, the higher the probability of obstruction, so the threshold is set to a smaller value (e.g., threshold = 0.8 meters when elevation angle < 10°); the lower the CNR, the worse the observation quality, so the threshold is set to a smaller value (e.g., threshold = 0.9 meters when CNR < 35 dB-Hz); combined with historical noise statistics (e.g., historical nighttime noise standard deviation = 0.7 meters), the final dynamic threshold is the minimum value calculated from the above factors.

[0192] Robustness mechanism activation: When the detected innovation amplitude exceeds a preset threshold, the robustness mechanism is activated. For example, for a satellite at a certain epoch with an elevation angle of 8°, a carrier-to-noise ratio of 32dB-Hz, a historical noise standard deviation of 0.7 meters, and a dynamic threshold of 0.7 meters, if the innovation amplitude is 0.9 meters (exceeding the threshold), the robustness mechanism is activated.

[0193] Weight reduction: A tiered strategy is adopted. Observations with slight bias (news exceeding the threshold but less than 1.2 times the threshold) are linearly weighted. For example, when the news = 0.77m (1.1 times the threshold), the weight = 1 - 0.1 × (0.77 - 0.7) / (1.2 × 0.7 - 0.7) = 0.9. Observations with significant outliers (news ≥ 1.2 times the threshold) are completely removed. For example, when the news = 0.84m (1.2 times the threshold), the weight = 0.

[0194] By using dynamic thresholds and hierarchical weighting, the robustness mechanism can be activated at the appropriate time, effectively blocking error propagation and avoiding the loss of effective information due to overprocessing, thereby improving the reliability of positioning results.

[0195] In summary, this invention accurately characterizes non-Gaussian noise (heavy-tailed, impulse characteristics) caused by occlusion and multipath effects in BeiDou positioning by constructing a Gaussian mixture model and a Student's t-distribution model (including a Cauchy distribution heavy-tailed correction term). It then uses an expectation-maximization algorithm to update model parameters online, addressing the problem that traditional Gaussian assumption models cannot adapt to the statistical characteristics of noise. Furthermore, this non-Gaussian noise model is fused with robust filtering. An adaptive robust weight matrix is ​​constructed using the Huber cost function, the IGGIII scheme, or noise probability density to dynamically adjust the observation covariance matrix and suppress the influence of outlier observations. This overcomes the shortcomings of classical robust filtering, such as lack of noise modeling and failure of weight adjustment under continuous strong interference, while avoiding the computational complexity of particle filtering. Finally, through iterative execution of robust filtering (time-based updates of predicted state and covariance, measurement updates to assess observation reliability, graded weight reduction or removal of outliers and state correction), it effectively solves the problems of insufficient BeiDou positioning accuracy and easy filtering divergence in complex environments, improving positioning robustness and real-time performance, and meeting the positioning reliability requirements of infrastructure monitoring such as highway slopes and bridges.

[0196] It should be noted that, in this document, relational terms such as "first" and "second" are used only to distinguish one entity or operation from another, and do not necessarily require or imply any such actual relationship or order between these entities or operations. Furthermore, the terms "comprising," "including," or any other variations thereof are intended to cover non-exclusive inclusion, such that a process, method, article, or apparatus that comprises a list of elements includes not only those elements but also other elements not expressly listed, or elements inherent to such process, method, article, or apparatus.

[0197] Although embodiments of the invention have been shown and described, it will be understood by those skilled in the art that various changes, modifications, substitutions and alterations can be made to these embodiments without departing from the principles and spirit of the invention, the scope of which is defined by the appended claims and their equivalents.

Claims

1. A robust Kalman filter positioning method for BeiDou based on a non-Gaussian noise model, characterized in that, Includes the following steps: S01: Obtain the raw observation data output by the BeiDou Navigation Satellite System. The raw observation data includes at least one of pseudorange observations, carrier phase observations, and Doppler frequency shift observations. The raw observation data is used for subsequent positioning calculations. S02: Establish a non-Gaussian noise model. For the observation noise in the original observation data and the process noise in the BeiDou positioning system, a noise model is preset to show that the observation noise and the process noise follow a non-Gaussian distribution in order to accurately characterize their statistical characteristics. The non-Gaussian noise model can adaptively update parameters online. S03: The non-Gaussian noise model and robust filtering are fused together, the fusion including: S031: Design the state equations and observation equations, wherein the state equations are expressed as follows: The observation equation is expressed as ,in, For the calendar The observation vector, Let the state vector be... Here is the state transition matrix. For the observation matrix, Let be the observed noise vector. The process noise vector is the state vector used to describe the dynamic parameters of the BeiDou receiver; S032: Based on the non-Gaussian noise model established in step S02, construct an adaptive robust weight matrix. The adaptive robust weight matrix is ​​used to dynamically adjust the covariance matrix of the observed values ​​during the filtering process to suppress the influence of abnormal observed values. S04: Iteratively execute the robust filtering, the iterative execution including: S041: Time update: Based on the state equation, use the optimal state estimate of the previous epoch to predict the prior estimate of the state vector of the current epoch, and calculate the corresponding prediction error covariance matrix. The prediction error covariance matrix reflects the uncertainty of the state estimate. S042: Measurement update, including: (a) Calculate the innovation sequence, which is the residual between the observation vector and the prior estimate of the state vector, and evaluate the reliability of the observation based on the non-Gaussian noise model. (b) Based on the adaptive robust weight matrix, abnormal observations in the innovation sequence that exceed a preset threshold are weighted to reduce their impact on the positioning results; (c) Update the posterior estimate of the state vector and the error covariance matrix, wherein the posterior estimate is the state estimate after correction of the observations; (d) Output the posterior estimate of the state vector to obtain a high-precision positioning result.

2. The BeiDou robust Kalman filter positioning method based on a non-Gaussian noise model according to claim 1, characterized in that, In step S02, the non-Gaussian noise model is a Gaussian mixture model, and the Gaussian mixture model affects the observed noise. probability density function Defined as: ； in, The number of mixture components in the Gaussian mixture model. For the first The weighting coefficients of each mixed component, and The first The mean vector and covariance matrix of each Gaussian component; The single Gaussian component The expression is: ； in, Let be the dimension of the observed noise.

3. The BeiDou robust Kalman filter positioning method based on a non-Gaussian noise model according to claim 1, characterized in that, In step S02, the non-Gaussian noise model is a Student's t-distribution model, and the Student's t-distribution model affects the observed noise. probability density function Defined as: ； in, Let be the dimension of the observed noise. Let be the degrees of freedom of the student t-distribution. Let be the mean vector of the student's t-distribution. Let be the scaling matrix of the student t-distribution. The student t-distribution model is used to describe observation noise with heavy-tailed characteristics, and is a gamma function.

4. The BeiDou robust Kalman filter positioning method based on a non-Gaussian noise model according to claim 2, characterized in that, In step S02, the parameters of the Gaussian mixture model are updated online using the expectation-maximization algorithm. The update equation of the expectation-maximization algorithm includes: ； in, To adjust the sliding window size, Represents the epoch Observations Belongs to the The posterior probabilities of the mixture components are used to implement adaptive parameter optimization of the Gaussian mixture model using the expectation-maximization algorithm.

5. The BeiDou robust Kalman filter positioning method based on a non-Gaussian noise model according to claim 1, characterized in that, In step S02, the non-Gaussian noise model further includes a heavy-tailed distribution correction term, which is corrected by introducing a Cauchy distribution. The corrected observation noise probability density function... Defined as: ； in, This is the main body of Gaussian. For Cauchy's corrections; The Cauchy distribution density function The expression is: ； in, For the noise dimension, The scaling matrix, As a mixing scaling factor, the heavy-tailed distribution correction term is intended to enhance the non-Gaussian noise model's ability to characterize impulse noise.

6. The BeiDou robust Kalman filter positioning method based on a non-Gaussian noise model according to claim 1, characterized in that, In step S032, the adaptive robust weight matrix is ​​constructed using the Huber cost function. Defined as: ； in, For the new interest, To adjust the parameters, the Huber cost function classifies the information, applies a secondary penalty to errors smaller than the adjustment parameters, and applies a linear penalty to errors larger than the adjustment parameters, in order to reduce the impact of outlier observations.

7. The BeiDou robust Kalman filter positioning method based on a non-Gaussian noise model according to claim 1, characterized in that, In step S032, the adaptive robust weight matrix is ​​constructed using the IGGIII scheme to dynamically adjust the observation weights. The weight function of the IGGIII scheme is... Defined as: ； in, To standardize the residuals, The first threshold, The second threshold is used as the weighting function to reduce the weight of the observed values ​​in segments according to the magnitude of the standardized residual. When the standardized residual exceeds the second threshold, the weight is zero, thereby realizing flexible handling of outliers of different degrees.

8. The BeiDou robust Kalman filter positioning method based on a non-Gaussian noise model according to claim 1, characterized in that, In step S032, the adaptive robust weight matrix The construction method is as follows: ； in, The number of observations, For the first The weight of each observation, For the first The probability density values ​​of each observation in the non-Gaussian noise model. To prevent the removal of zero constants, the weight matrix evaluates the reliability of the observations in real time using the probability density of the non-Gaussian noise model, and adaptively weights the observation data accordingly to enhance the robustness of the filter.

9. The BeiDou robust Kalman filter positioning method based on a non-Gaussian noise model according to claim 1, characterized in that, In step S031, the state vector The state vector, comprising at least four of the following: three-dimensional position, three-dimensional velocity, receiver clock bias, and receiver clock drift, is used to comprehensively describe the dynamic motion state of the BeiDou receiver in space, as well as its clock bias and drift; the raw observation data The observations include pseudorange and carrier phase observations from multiple frequency bands such as BeiDou B1I, B2I, and B3I. These multi-frequency band observations are used for high-precision state updates through the observation equations.

10. The BeiDou robust Kalman filter positioning method based on a non-Gaussian noise model according to claim 1, characterized in that, In step S042(b), when the amplitude of the new information is detected to exceed the preset threshold, the robustness mechanism is activated. The preset threshold is dynamically determined based on the satellite elevation angle, carrier-to-noise ratio, and historical noise statistical characteristics to adapt to the complex and ever-changing actual positioning environment. Furthermore, the weighting process uses linear weighting for observations with slight deviations and completely removes observations with significant outliers. The weighting process aims to retain effective information while blocking error propagation paths to ensure the reliability of the positioning results.

Images