Beidou robust kalman filter positioning method based on non-gaussian noise model

Abstract

The present application relates to the field of satellite navigation and positioning technology, and particularly relates to a Beidou robust Kalman filtering positioning method based on a non-Gaussian noise model. The present application precisely characterizes non-Gaussian noise generated in Beidou positioning due to shielding and multipath effects by constructing a mixed Gaussian model and a Student's t-distribution model, and updates model parameters online using an expectation maximization algorithm to solve the problem that a traditional Gaussian assumption model cannot adapt to noise statistical characteristics. The non-Gaussian noise model is then fused with a robust filter to construct an adaptive robust weight matrix through a Huber cost function, an IGGIII scheme or a noise probability density, dynamically adjust an observation covariance matrix and suppress the influence of abnormal observation values to overcome the defects of a classical robust filter, such as lack of noise modeling and failure of weight adjustment under continuous strong interference, while avoiding the problem of complex calculation of a particle filter.

Description

TECHNICAL FIELD

[0001] The application relates to the technical field of satellite navigation and positioning, in particular to a Beidou anti-bias Kalman filtering positioning method based on a non-Gaussian noise model. BACKGROUND

[0002] The Beidou satellite navigation system has important application value in the field of health monitoring of major infrastructures, especially for the safety evaluation of key structures such as highway slopes and large bridges. Such infrastructures are affected by environmental loads, geological activities and other factors during long-term operation, and will produce slight deformation, which needs to be monitored at the millimeter level through high-precision positioning technology. Compared with traditional measurement methods, structure monitoring based on satellite navigation has the advantages of all-weather and automation, but also faces the challenge of measurement in special environments.

[0003] In actual engineering applications, there are often complex electromagnetic environments around highway slopes and bridges. The undulating terrain of the slope leads to satellite signal obstruction, and the main body of the bridge structure and the surrounding buildings bring significant multipath effects, so that the observation noise often deviates from the conventional Gaussian distribution, showing sharp peak pulses and heavy tail characteristics. This condition can seriously affect the positioning accuracy of the Kalman filtering algorithm, and even may lead to false deformation warning. The existing filtering method based on the Gaussian assumption cannot distinguish between real displacement and noise interference, and cannot meet the strict requirements of infrastructure safety monitoring on data reliability. SUMMARY

[0004] In view of the deficiencies of the prior art, the Beidou anti-bias Kalman filtering positioning method based on a non-Gaussian noise model is provided, which solves the problems of insufficient positioning accuracy, easy occurrence of false deformation warning and easy divergence of filtering caused by non-Gaussian characteristics of observation noise when compared with the prior art.

[0005] To achieve the above purpose, the following technical scheme is adopted: the Beidou anti-bias Kalman filtering positioning method based on a non-Gaussian noise model comprises the following steps:

[0006] S01: Obtain the original observation data output by the Beidou satellite navigation system, the original observation data comprising at least one of pseudorange observation values, carrier phase observation values and Doppler frequency shift observation values, the original observation data being used for subsequent positioning calculation;

[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, preset a noise model that the observation noise and the process noise conform to a non-Gaussian distribution, to accurately represent the statistical characteristics thereof, the non-Gaussian noise model being capable of online adaptive updating of parameters;

[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 a posteriori probability of a mixture component, the expectation maximization algorithm being used to implement adaptive parameter optimization of the Gaussian mixture model.

[0030] Further, in step S02, the non-Gaussian noise model further comprises a heavy-tailed distribution correction term, the heavy-tailed distribution correction term being corrected by introducing a Cauchy distribution, the corrected observation noise probability density function is defined as:

[0031] ;

[0032] wherein, is a Gaussian body part, is a Cauchy correction part;

[0033] The Cauchy distribution density function is expressed as:

[0034] ;

[0035] wherein, is the noise dimension, is a scale matrix, is a mixing proportion coefficient, the heavy-tailed distribution correction term aiming to enhance the representation ability of the non-Gaussian noise model to impulse noise.

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

[0037] ;

[0038] wherein, is a innovation, is a tuning parameter, the Huber cost function classifying the innovation, performing quadratic penalty on errors smaller than the tuning parameter, and performing linear penalty on errors larger than the tuning parameter, so as to reduce the influence of abnormal observations.

[0039] Further, in step S032, the adaptive robust weight matrix is constructed using an IGGIII scheme to dynamically adjust the observation weight, a weight function of the IGGIII scheme is defined as:

[0040] ;

[0041] wherein, is a standardized residual, is a 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] The application solves the problem that the traditional Gaussian assumption model cannot adapt to the statistical characteristics of noise by constructing a mixed Gaussian model and a student t distribution model to accurately represent the non-Gaussian noise caused by occlusion and multipath effect in Beidou positioning, and updating the model parameters online by using the expectation maximization algorithm; the non-Gaussian noise model is fused with the robust filtering, an adaptive robust weight matrix is constructed by using the Huber cost function, the IGGIII scheme or the noise probability density, the observation covariance matrix is dynamically adjusted, the influence of abnormal observation values is suppressed, the defects of the classic robust filtering, such as lack of noise modeling and failure of weight adjustment under continuous strong interference, are overcome, and the problem of complex calculation of particle filtering is avoided; finally, the problem of insufficient Beidou positioning accuracy and easy divergence of filtering in a complex environment is effectively solved by iteratively performing robust filtering, the positioning robustness and real-time performance are improved, and the demand for positioning reliability of highway slope, bridge and other infrastructure monitoring is met. BRIEF DESCRIPTION OF DRAWINGS

[0049] Figure 1 The method flowchart of the application. DETAILED DESCRIPTION

[0050] The technical solutions in the embodiments of the application will be clearly and completely described below with reference to the drawings in the embodiments of the application. Obviously, the described embodiments are only part of the embodiments of the application, rather than all the embodiments of the application. Based on the embodiments in the application, all other embodiments obtained by those skilled in the art without creative labor fall within the protection scope of the application.

[0051] Please refer to Figure 1 The application provides a Beidou robust Kalman filtering positioning method based on a non-Gaussian noise model, which comprises the following steps:

[0052] S01: Obtain original observation data output by a Beidou satellite navigation system, wherein the original observation data comprises at least one of pseudo-range observation values, carrier phase observation values and Doppler frequency shift observation values, and the original observation data is used for subsequent positioning calculation;

[0053] S02: Establish a non-Gaussian noise model, preset a noise model that the observation noise in the original observation data and the process noise in the Beidou positioning system are subject to a non-Gaussian distribution, so as to accurately represent the statistical characteristics, and the non-Gaussian noise model can update parameters online and adaptively;

[0054] S03: Fuse the non-Gaussian noise model with robust filtering, and the fusion comprises:

[0055] S031: Design a state equation and an observation equation, wherein the state equation is represented as , and the observation equation is represented 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] The high-precision receiver supporting Beidou multi-frequency observation (for example, a certain type of Beidou dual-frequency receiver) is used to set up the receiver at the highway slope monitoring point, the sampling interval is set to 1Hz (adapted to the slow deformation characteristics of the slope), and the original observation data output by the Beidou satellite navigation system are collected, including the pseudo-range observation values, carrier phase observation values and Doppler frequency shift observation values of three frequency bands B1I, B2I and B3I of Beidou. Among them, the pseudo-range observation value reflects the distance observation information from the satellite to the receiver, the carrier phase observation value provides higher precision distance change information, and the Doppler frequency shift observation value assists in representing the dynamic motion state of the receiver, providing multi-dimensional data support for subsequent positioning calculation.

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

[0068] For the observation noise in the original observation data and the process noise of the Beidou positioning system, the observation noise and the process noise are preset to follow a non-Gaussian noise model of mixed Gaussian distribution (the mixed Gaussian model can accurately represent the noise multi-peak and heavy tail characteristics caused by the multi-path effect), and the model parameters support online adaptive update. The probability density function of the mixed Gaussian model is:

[0069] ;

[0070] Among them, is the observation noise vector at the epoch , with the dimension of meter; is the number of mixed components, for example, taking 3 (corresponding to three main sources of noise under the multi-path effect); is the weight coefficient of the th mixed component, satisfying , dimensionless; is the mean vector of the th Gaussian component, with the same dimension as , and the initial value can be obtained by statistical calculation of historical clean observation data; is the covariance matrix of the th Gaussian component, with the dimension of square meter, and the initial value is also calculated based on historical data; is the probability density function of a single Gaussian component, and the expression is:

[0071] ;

[0072] Among them, is the dimension of the observation noise, for example, when B1I and B2I pseudo-range observations are used at the same time, , dimensionless; is the determinant of the covariance matrix , with the same dimension as ; 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 The element of the matrix is 1. The observation noise vector is consistent with the observation noise in step S02.

[0084] Step S032: Construct adaptive robust weight matrix

[0085] Based on the Gaussian mixture model in step S02, construct an adaptive robust weight matrix for dynamic adjustment of the observation covariance matrix. The weight matrix is a diagonal matrix, and the expression is:

[0086] ;

[0087] Wherein, The number of observations, for example (corresponding to B1I, B2I pseudo-range); The weight of the first observation, the calculation formula is:

[0088] ;

[0089] The probability density value of the first observation noise under the Gaussian mixture model, dimensionless; The constant is prevented from being zero, dimensionless. When the observation value is abnormal, small, followed by a decrease, achieving dynamic weight reduction of abnormal observation values.

[0090] Step S04: Iterative execution of robust filtering

[0091] Step S041: Time update

[0092] According to the state equation, using the optimal state estimation and error covariance matrix of the last epoch , the state vector prior estimate value and the predicted error covariance matrix of the current epoch , the formulas are respectively:

[0093] ;

[0094] ;

[0095] Wherein, The error covariance matrix of the last epoch state estimation, the dimension is the square of the corresponding state parameter; is the process noise covariance matrix, dimension of which is consistent with and its value is determined based on the parameters of the Gaussian Mixture Model; is the transpose of the state transition matrix. reflects the uncertainty of the prior estimate of the state and provides the basis for the subsequent measurement update.

[0096] Step S042: Measurement Update

[0097] o (a) Calculate the innovation sequence , whose formula is:

[0098] ;

[0099] wherein, is the observation vector, is the observation prediction value, both of which are in meters and have consistent dimensions, and thus do not need to be normalized. In combination with the Gaussian Mixture Model of step S02, the probability density of each element of the innovation is calculated. If is less than a preset probability threshold (for example, 0.05), it is determined that the corresponding observation value is of low reliability. (b) Based on the adaptive robust weight matrix of step S032 , the observation value of low reliability is processed with reduced weight. For example, when the observation value corresponds to

[0100] , it is considered to be an abnormal observation value, and its contribution in the observation covariance matrix is reduced through the weight matrix, so as to reduce its interference on the positioning result. (c) Update the posterior estimate value of the state vector and the error covariance matrix

[0101] , whose formulas are respectively: ;

[0102] ;

[0103] ;

[0104] ;

[0105] wherein, is the Kalman gain matrix, which is dimensionless; is the original observation noise covariance matrix, which is in meters²; is the unit matrix, which is dimensionless; is the inverse matrix of the weight matrix, which adjusts the effective contribution of the observation covariance matrix through the weight.

[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 mixture of Gaussian models 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. is a scale matrix of the Student's t-distribution, is a gamma function, the Student's t-distribution model is used to describe the observation noise with heavy-tailed characteristics.

[0124] Specifically, in this embodiment, the non-Gaussian noise model adopts the Student's t-distribution model, which focuses on the pulse noise (heavy-tailed characteristics) caused by instantaneous electromagnetic interference in Beidou positioning, for example, the influence of electromagnetic pulses generated by high-voltage lines around the highway on observation data.

[0125] The Student's t-distribution model defines the probability density function of the observation noise as follows:

[0126]

[0127] The parameters are explained as follows: is the dimension of the observation noise, for example, when only Beidou B1I pseudorange observations are used, ; is the degree of freedom of the Student's t-distribution, the value of which determines the degree of heavy tail of the noise distribution, for example, when there are frequent electromagnetic pulses, take (the smaller the degree of freedom, the more obvious the heavy tail, and the more suitable for pulse noise); is the mean vector of the Student's t-distribution, which is consistent in dimension with (meters), and is usually set to 0, representing the baseline level of the noise; is the scale matrix of the Student's t-distribution, which is in the dimension of meters², for example, according to the observation data of the history of the interference-free period ( is a scalar); is a gamma function, used to normalize the probability density, for example , .

[0128] Formula dimension check: the dimension is meters, the dimension is meters , and after multiplication, it is dimensionless, together with (dimensionless), to form a dimensionless item in the parentheses, the exponential part is dimensionless, and the final dimension of the entire formula is meters (the dimension of the probability density), which is consistent in dimension.

[0129] Through the Student's t-distribution model, the heavy-tailed characteristics of pulse noise can be effectively characterized, avoiding the underestimation of pulse noise by traditional Gaussian models, thereby reducing the state estimation bias caused by pulse noise in the filtering process and improving the anti-interference ability of positioning.

[0130] ​In the embodiment, in step S02, the parameters of the Gaussian mixture model are updated online by an expectation maximization algorithm, and an update equation of the expectation maximization algorithm comprises:

[0131] ;

[0132] wherein, is a size of a sliding window, denotes an observation value at an epoch , and is a posterior probability of the observation value belonging to a mixed component, the expectation maximization algorithm is used to realize adaptive parameter optimization of the Gaussian mixture model.

[0133] Specifically, in the embodiment, the parameters of the Gaussian mixture model are updated online by an expectation maximization (EM) algorithm to adapt to dynamic changes of Beidou positioning noise characteristics, for example, multipath noise characteristics change caused by vegetation coverage change in different seasons in highway slope monitoring.

[0134] An update equation of the EM algorithm is as follows:

[0135] ;

[0136] The parameters are explained as follows: is a size of a sliding window, determined according to a noise change rate, for example, taken as in a vegetation seasonal change scenario; is an observation noise at an epoch ; is a posterior probability of the observation value belonging to a mixed component, ; is an updated mixed component weight coefficient, , dimensionless, satisfying ; is an updated mean vector of a Gaussian component, , the dimension is consistent with (meters); is an updated covariance matrix of a Gaussian component, , the dimension is meters2; is a number of mixed components, consistent with in step S02, for example, .

[0137] The update process is as follows: first, at each epoch , the (E step) is calculated according to the current model parameters; then, based on data of epochs in the sliding window, the , , (M step). For example, when the vegetation coverage increases, the noise contribution of multiple reflection paths rises, will be adjusted from 0.1 to 0.2, will also increase accordingly, ensuring that the model always adapts to the current noise characteristics.

[0138] By updating online through the EM algorithm, the Gaussian mixture model can dynamically track noise changes and avoid model mismatch caused by parameter solidification, 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 modified by introducing a Cauchy distribution, and the modified observation noise probability density function is defined as:

[0140] ;

[0141] wherein, is the Gaussian main part, is the Cauchy correction part;

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

[0143] ;

[0144] wherein, is the noise dimension, is the scale matrix, is the mixing proportion coefficient, and the heavy-tailed distribution correction term aims to enhance the representation ability of the non-Gaussian noise model to impulse noise.

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

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

[0147] ;

[0148] The parameters are explained as follows: is the mixing proportion coefficient, used to adjust the contribution of the Cauchy correction part, for example, taking during periods of frequent construction interference, and dimensionless; is the Gaussian main part, the expression of which is consistent with step S02, is the covariance matrix of the Gaussian component, dimension of meter is the Cauchy correction part, used to characterize extreme impulse noise.

[0149] The expression of the Cauchy distribution density function is as follows:

[0150]

[0151] wherein, is the dimension of the noise, for example (B1I, B2I pseudo-range); is the scale matrix of the Cauchy distribution, dimension consistent with (meter ), for example , the larger the scale, the wider the tail of the Cauchy distribution, and the better the adaptation to strong impulse noise; is the determinant of the scale matrix .

[0152] Formula dimension check: dimension of meter, dimension of meter , dimensionless, the items in the brackets are dimensionless, and the entire Cauchy distribution density function has a dimension of meter , consistent with the dimension of the Gaussian main part, ensuring dimensional uniformity.

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

[0154] In this embodiment, in step S032, the construction of the adaptive robust weight matrix adopts a Huber cost function, and the Huber cost function is defined as:

[0155]

[0156] wherein, is the innovation, is an adjustment parameter, and the Huber cost function classifies the innovation, performs quadratic punishment on errors smaller than the adjustment parameter, and performs linear punishment on errors larger than the adjustment parameter, so as to reduce the influence of abnormal observation values.

[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 to determine a serious abnormality, and the weight function is taken as 0.5.

[0167] Specifically, in the embodiment, the IGGIII scheme is adopted to dynamically adjust the observation weight, so as to flexibly process abnormal observation values of different degrees. For example, in the monitoring of the slope of the highway, the slight multipath interference and the observation abnormality caused by serious shielding are distinguished.

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

[0169]

[0170] The parameters are explained as follows: The standardized residual is calculated by The innovation is The observation noise standard deviation is , and the dimension is meter, and the dimension is dimensionless after standardization. The first threshold is used to determine a slight abnormality, and the weight function is taken as 0.5. ; The second threshold is used to determine a serious abnormality, and the weight function is taken as 0.5. .

[0171] The weight adjustment example is as follows: when (smaller than ), (normal weight); when (between and ), (significantly reduced weight); and when (larger than ), (completely removed).

[0172] Dimension check of the formula: The standardized residual is dimensionless, , are dimensionless thresholds, which ensure that the weight calculation result is dimensionless, and the value is between 0 and 1, which meets the weight definition.

[0173] Through the IGGIII scheme, the slight abnormal observation value can be gradually reduced in weight, and the serious abnormal observation value can be completely removed, so as to avoid the interference of the abnormal value and maximize the retention of effective observation information, and the adaptability of the filter to the complex environment is improved.

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

[0175] ;

[0176] wherein, is the number of observations, is the weight of the th observation, is the probability density value of the th observation innovation under the non-Gaussian noise model, is a constant to prevent zero, the weight matrix evaluates the reliability of the observation in real time by using the probability density of the non-Gaussian noise model, and the observation data is adaptively weighted accordingly to enhance the robustness of the filter.

[0177] Specifically, in the embodiment, the adaptive robust weight matrix evaluates the observation reliability in real time by the probability density of the non-Gaussian noise model, realizes more accurate weighting adjustment, and for example, in a bridge bottom monitoring scene with complex multipath effect.

[0178] The construction method of the adaptive robust weight matrix is as follows:

[0179]

[0180] The parameters are explained as follows: is the number of observations, for example, 3 (corresponding to Beidou B1I, B2I and B3I pseudo ranges); is the weight of the th observation, dimensionless; is the probability density value of the th observation innovation under the non-Gaussian noise model (such as a mixed Gaussian model), dimensionless; is a constant to prevent zero, for example, 0.0001, to avoid abnormal weight calculation due to too small .

[0181] Weight calculation example: in bridge bottom monitoring, the probability density of a certain B1I pseudo range observation innovation is high (high reliability), and ; the probability density of a certain B2I pseudo range observation innovation is low (low reliability), and .

[0182] Formula dimension check: and are dimensionless, the weight is dimensionless, and the dimensions are consistent in 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 It includes pseudorange observations 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, the carrier phase observations of B1I (wavelength of about 0.19 meters) and B2I (wavelength of about 0.16 meters) can be combined to eliminate ionospheric errors and improve positioning accuracy.

[0188] By selecting the state vector and the original observation data, the dynamic parameters and high-precision observation information required for infrastructure monitoring can be comprehensively covered, and data support can be provided for millimeter-level positioning.

[0189] In the embodiment, in step S042(b), the robustness mechanism is activated when it is detected that the innovation amplitude exceeds the preset threshold, and the preset threshold is dynamically determined based on the satellite elevation angle, the carrier-to-noise ratio, and the historical noise statistical characteristics to adapt to complex and variable actual positioning environments; and the weight reduction processing adopts linear weight reduction processing for slightly deviated observation values and completely eliminates significantly outlying observation values, and the weight reduction processing aims to retain effective information while blocking error propagation paths to ensure the reliability of the positioning results.

[0190] Specifically, in the embodiment, the activation of the robustness mechanism and the weight reduction processing are based on dynamic thresholds and hierarchical strategies to adapt to complex and variable positioning environments, for example, to cope with satellite elevation angle changes and carrier-to-noise ratio fluctuations in different time periods (such as daytime and nighttime) in highway slope monitoring.

[0191] Determination of the preset threshold: dynamically calculated based on the satellite elevation angle, the carrier-to-noise ratio, and the historical noise statistical characteristics. For example, the lower the satellite elevation angle, the higher the probability of obstruction, and the threshold is set to a smaller value (such as when the elevation angle < 10°, the threshold = 0.8 meters); the lower the carrier-to-noise ratio, the worse the observation quality, and the threshold is set to a smaller value (such as when the carrier-to-noise ratio < 35 dB-Hz, the threshold = 0.9 meters); combined with the historical noise statistical value (such as the historical nighttime noise standard deviation = 0.7 meters), the final dynamic threshold is the minimum value of the calculation results of the above factors.

[0192] Activation of the robustness mechanism: when it is detected that the innovation amplitude exceeds the preset threshold, the robustness mechanism is activated. For example, in a certain epoch, the satellite elevation angle = 8°, the carrier-to-noise ratio = 32 dB-Hz, the historical noise standard deviation = 0.7 meters, the dynamic threshold = 0.7 meters, and if the innovation = 0.9 meters (exceeds the threshold), the robustness mechanism is activated.

[0193] Weight reduction processing: hierarchical strategy is adopted, linear weight reduction is adopted for slightly deviated observation values (innovation exceeds the threshold but is less than 1.2 times the threshold), for example, when the innovation = 0.77 meters (1.1 times the threshold), the weight = 1 - 0.1x(0.77 - 0.7) / (1.2x0.7 - 0.7) = 0.9; and completely eliminated for significantly outlying observation values (innovation ≥ 1.2 times the threshold), for example, when the innovation = 0.84 meters (1.2 times the threshold), the weight = 0.

[0194] By dynamic threshold and hierarchical weight reduction, the robustness mechanism can be activated at the right time, effectively blocking error propagation and avoiding effective information loss caused by excessive processing, and improving the reliability of the positioning result.

[0195] In summary, the application accurately characterizes the non-Gaussian noise (heavy-tailed, impulse characteristics) generated by occlusion and multipath effects in Beidou positioning by constructing a mixed Gaussian model and a student t distribution model (including a Cauchy distribution heavy-tail correction term), and updates the model parameters online using the expectation maximization algorithm, solving the problem that the traditional Gaussian assumption model cannot adapt to the statistical characteristics of the noise; then the non-Gaussian noise model is fused with the robust filtering, the adaptive robust weight matrix is constructed by Huber cost function, IGGIII scheme or noise probability density, the observation covariance matrix is dynamically adjusted, the influence of abnormal observation values is suppressed, the defects of the classical robust filtering such as lack of noise modeling and weight adjustment failure under continuous strong interference are overcome, while the problem of complex calculation of particle filtering is avoided; finally, the robust filtering is iteratively executed (time update to predict the state and covariance, measurement update to evaluate the reliability of the observation, hierarchical weight reduction or rejection of abnormal values and correction of the state), effectively solving the problems of insufficient Beidou positioning accuracy and filter divergence in complex environments, improving the positioning robustness and real-time performance, and meeting the demand for positioning reliability in monitoring of highway slopes, bridges and other infrastructure.

[0196] It should be noted that the relational terms herein such as first and second and the like are used solely to distinguish one entity or action from another entity or action without necessarily requiring or implying any such actual relationship or order between such entities or actions. Moreover, the terms "comprises", "comprising", or any other variations thereof, are intended to cover a non-exclusive inclusion, such that a process, method, article, or apparatus that comprises a list of elements does not include only those elements but can include other elements not expressly listed or inherent to such process, method, article, or apparatus.

[0197] Although embodiments of the present application have been shown and described, it is to be understood that various modifications, substitutions, replacements and changes can be made to these embodiments without departing from the principles and spirit of the present application, and the scope of the present application 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, For state vectors, 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; 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 observations in segments according to the magnitude of the standardized residuals. When the standardized residuals exceed the second threshold, the weights are zero, thus enabling 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 It includes pseudorange observations and carrier phase observations from multiple frequency bands of 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