Hierarchical bayesian estimation method and device for measuring noise covariance and storage medium

Abstract

The application discloses a hierarchical Bayesian measurement noise covariance estimation method and device and a storage medium, and belongs to the field of state estimation and filtering. In view of the problems of insufficient anisotropic noise modeling capability, confusion deviation of process noise and measurement noise, and lack of unbiased repair strategy for negative definite result in the prior art, a global inverse Wishart and a dimension-by-dimension gamma double-layer hierarchical Bayesian prior is constructed, an unbiased empirical estimate is obtained by subtracting a prediction uncertainty term from an innovation outer product through Cholesky factor triangular solution, and a positive semi-definite repair is projected on an eigenvalue, a forgetting factor is used to drive two-layer prior and posterior parameter recursive update, and a measurement noise covariance matrix estimation is obtained by weighted fusion of two-layer posterior expectation adaptive weights. The application is suitable for target tracking, navigation positioning, sensor fusion and other scenes requiring recursive state estimation and unknown or time-varying measurement noise, and can serve the fields of industrial internet of things information sensing and sensing technology service and industrial information and data processing service.

Description

Technical Field

[0001] This invention belongs to the field of state estimation and filtering technology, and relates to industrial Internet of Things (IoT) information sensing, sensing technology services, and industrial information and data processing services. More specifically, this invention relates to a method, apparatus, and storage medium for estimating the hierarchical Bayesian measurement noise covariance. Background Technology

[0002] Kalman filters are a fundamental tool for state estimation in target tracking, navigation and positioning, and multi-sensor fusion. When the system is linear and the noise follows a Gaussian distribution, the Kalman filter provides the optimal state estimate in the sense of minimum mean square error. However, this optimality depends strictly on the precise knowledge of the process noise covariance matrix and the measurement noise covariance matrix. In practical engineering applications, the measurement noise covariance matrix is ​​often unknown or varies over time. The ranging accuracy of sensors decreases with increasing target distance, the noise level of infrared sensors drifts with changes in ambient temperature, and the noise characteristics of different sensor channels in a multi-sensor fusion system can differ by several orders of magnitude. When the measurement noise covariance matrix is ​​inaccurate, the estimation accuracy of the Kalman filter will decrease significantly, and may even lead to filter divergence.

[0003] To address the problem of unknown or time-varying measurement noise covariance matrices, various adaptive estimation methods have been developed in existing technologies. The variational Bayesian R-adaptive method models the measurement noise covariance matrix as an inverse Wissaud distribution and uses variational Bayesian inference to recursively update its posterior distribution from the innovation sequence. This method uses a single degree-of-freedom parameter to control the estimation of all measurement dimensions, assuming that the noise from each measurement channel follows the same prior distribution. This assumption often fails in engineering practice. For example, in sensor target tracking, the variances of distance measurement noise and angle measurement noise may differ by six orders of magnitude. A uniform degree-of-freedom parameter cannot simultaneously match these two vastly different noise characteristics, causing the estimation accuracy of smaller noise dimensions to be hampered by the larger noise dimension.

[0004] The Sage-Husa subtraction estimation method estimates the measurement noise covariance matrix by using the difference between the innovation covariance and the theoretical innovation covariance. This method directly uses the innovation outer product as an empirical estimate of the measurement noise covariance matrix. Innovation is the difference between the actual and predicted measurements, encompassing contributions from both prediction errors caused by process noise and measurement noise itself. When the innovation outer product is directly used to estimate measurement noise, the contribution of prediction errors is inevitably included in the measurement noise estimation, leading to a systematic overestimation of the measurement noise covariance matrix. This problem is further exacerbated when the process noise covariance matrix itself is uncertain, creating a positive feedback loop between the measurement noise and process noise estimates, making it difficult to maintain filter consistency.

[0005] Some methods assume that the measurement noise covariance matrix is ​​a diagonal matrix and estimate the noise variance independently for each channel on a dimension-by-dimensional basis. These methods abandon the correlation information between measurement channels, and this assumption introduces model error when off-diagonal covariances exist between measurement channels. Furthermore, these methods are highly sensitive to anomalous information; anomalies in a single frame can cause drastic jumps in the noise variance estimate for a particular dimension.

[0006] Recently, methods have emerged that combine the inverse Wieshard prior with the gamma prior into a hierarchical Bayesian framework, but their technical purpose is mainly focused on modeling heavy-tailed distributions or detecting outliers. The hierarchical priors of these methods serve to model the shape of the noise distribution, rather than to handle the magnitude differences between different dimensions in the measurement noise covariance matrix.

[0007] Furthermore, all the aforementioned methods lack reasonable repair strategies when negative eigenvalues ​​occur in the subtraction estimation. When the filter's prediction covariance is large, the prediction uncertainty term may numerically exceed the innovation outer product, leading to negative eigenvalues ​​in the subtraction result, which is physically unreasonable. Existing methods mainly handle this situation by falling back to the original innovation outer product or using a half-decay mixing strategy, both of which introduce uncontrollable systematic biases into the estimation.

[0008] In summary, the existing technology lacks an adaptive estimation method for measurement noise covariance that can simultaneously handle the differences in noise levels across measurement dimensions, eliminate the confusion between process noise and measurement noise, and possess a reasonable negative definite repair strategy. Summary of the Invention

[0009] The purpose of this invention is to provide a method, device, and storage medium for estimating hierarchical Bayesian measurement noise covariance. By constructing a two-layer hierarchical Bayesian prior structure consisting of a global inverse Wissaud distribution layer and a dimension-wise gamma distribution layer, independent adaptive modeling of each measurement dimension is achieved, solving the problem of noise magnitude differences between measurement dimensions. Unbiased innovation decomposition explicitly subtracts the prediction uncertainty term from the outer product of the innovation, eliminating the confusion between process noise and measurement noise. Eigenvalue projection positive semi-definite repair preserves all effective information in the positive semi-definite directions while discarding only the non-physical solutions in the negative definite directions, ensuring asymptotic unbiasedness of the estimation and solving the problem of lacking unbiased repair strategies for negative definite results. Hierarchical fusion achieves an automatic balance between global structural constraints and dimension-wise independent adaptation, eliminating the need for manual setting of independent learning rates for each dimension. This results in an unbiased and anisotropic online adaptive estimation of the measurement noise covariance matrix, suitable for applications requiring recursive state estimation and where the statistical characteristics of measurement noise are unknown or time-varying, such as target tracking, navigation and positioning, and sensor fusion. It can be applied to industrial IoT information sensing and sensing technology services and industrial information and data processing services.

[0010] To address the aforementioned problems and achieve the objectives and other advantages of this invention, a method for estimating the covariance of hierarchical Bayesian measurement noise is provided, comprising the following steps:

[0011] S1, Constructing the hierarchical Bayesian prior structure: This is the measurement noise covariance matrix. Construct a two-level hierarchical Bayesian prior, where m is the dimension of the measurement vector;

[0012] The two-level hierarchical Bayesian prior includes:

[0013] Global Inverse Wissaud Distribution Layer: Measurement Noise Covariance Matrix Obeying the distribution of the reverse Saudi Arabia ,in The global scale matrix, For global degrees of freedom;

[0014] Dimensionally Gamma Distribution Layer: For each measurement dimension i, the measurement noise covariance matrix diagonal elements corresponding precision Follows a gamma distribution ,in For shape parameters, For rate parameters, Let molecule be the shape parameter of the gamma distribution in the i-th measurement dimension, satisfying , Let be the rate parameter molecule of the gamma distribution in the i-th measurement dimension;

[0015] S2, Unbiased innovation decomposition: At each recursive time k, obtain the innovation vector of the current frame. ,in For the current measurement vector, For measurement prediction vectors;

[0016] According to the definition of new information covariance:

[0017] ;

[0018] The superscript T indicates matrix transpose. To measure the Jacobian matrix, For predicting the state covariance matrix;

[0019] The unbiased empirical estimation matrix is ​​calculated using the following sub-steps:

[0020] S21, Obtain the prediction covariance Cholesky factor for the current frame. ,satisfy and the measurement-state cross-covariance matrix ;

[0021] S22, by Solve the trigonometric equations by substituting forward. The intermediate variable matrix is ​​obtained. ,but:

[0022] ;

[0023] S23, Calculate the unbiased empirical estimation matrix:

[0024] ;

[0025] S3, Eigenvalue projection positive semidefinite repair: Based on the result obtained in S23... The matrix is ​​symmetricized, and then eigenvalue decomposition is performed on the symmetricized matrix. The negative eigenvalues ​​of the decomposed matrix are projected to zero to reconstruct the empirically estimated matrix that satisfies the positive semidefinite requirement. ;

[0026] S4, Global Inverse Wissaud Distribution Layer Posterior Recursive Update: The hyperparameters of the global inverse Wissaud prior are updated exponentially by using the forgetting factor ρ to drive the exponential moving average recursive update.

[0027] , ;

[0028] in and These represent the global degrees of freedom and the global scale matrix at the current time k, respectively. and These are the global degrees of freedom and the global scale matrix at the previous time step k-1, respectively.

[0029] like Then, based on the mean formula of the inverse Wissaud distribution, the full matrix estimate of the global inverse Wissaud distribution layer posterior expectation of the current frame is calculated:

[0030] ;

[0031] like Then the global inverse Wissaud distribution layer posterior expectation full matrix estimate of the current frame Keep the value from the previous moment unchanged;

[0032] Output the full matrix estimate of the global inverse Wissaud distribution posterior expectation of the current frame. ;

[0033] S5, Posterior recursive update of the dimensional gamma distribution layer: For each measurement dimension i, use the innovation vector The i-th component The square value Perform recursive updates of the underlying prior parameters:

[0034] , ;

[0035] in and Let be the shape parameter molecules of the i-th measurement dimension at the current time k and the previous time k-1, respectively; and Let be the numerators of the rate parameters for the i-th measurement dimension at the current time k and the previous time k-1, respectively.

[0036] like Then, calculate the posterior expectation estimate for that dimension:

[0037] ;

[0038] S6, Hierarchical Fusion: S6 receives the global inverse Wissaud distribution layer posterior expectation full matrix estimation from the output of S4. Construct the measurement noise covariance matrix of the final output of the current frame. ,in Represents the measurement noise covariance matrix The i-th diagonal element, Represents the measurement noise covariance matrix The non-diagonal element in the i-th row and j-th column; This represents the full matrix estimation of the posterior expectation of the global inverse Wissaud distribution layer. The i-th diagonal element; This represents the full matrix estimation of the posterior expectation of the global inverse Wissaud distribution layer. The non-diagonal element in the i-th row and j-th column;

[0039] For each measurement dimension i, determine The i-th diagonal element:

[0040] like ,but:

[0041] ;

[0042] Where the fusion weight coefficient ω∈[0.5,0.8];

[0043] like ,but:

[0044] ;

[0045] For all off-diagonal elements, i.e., i ≠ j and i,j∈{1,…,m}, the full matrix estimate of the posterior expectation of the global inverse Wissaud distribution layer is provided:

[0046] ;

[0047] At this point, the measurement noise covariance matrix... All elements have been determined, and then the entire measurement noise covariance matrix is... Perform symmetry transformation and restrict each diagonal element to a preset range;

[0048] S7, Output measurement noise covariance matrix This is used for the measurement update step of the filter in the next frame.

[0049] Preferably, in the hierarchical Bayesian method for estimating the measurement noise covariance, step S6 involves estimating the entire measurement noise covariance matrix. Perform symmetry transformation, specifically by transforming the measurement noise covariance matrix. Assigned value .

[0050] Preferably, in the method for estimating the hierarchical Bayesian measurement noise covariance, the initial values ​​of each hyperparameter in S1 are set as follows:

[0051] Initial values ​​of global degrees of freedom ;

[0052] Initial value of the global scale matrix ,in The initial measurement noise covariance matrix;

[0053] Initial shape parameters of gamma distribution in each measurement dimension ;

[0054] Initial rate parameters of gamma distribution in each measurement dimension ,in The initial measurement noise covariance matrix The diagonal elements.

[0055] Preferably, in the hierarchical Bayesian method for estimating the noise covariance, the eigenvalue projection positive semidefinite repair in S3 specifically includes:

[0056] right Perform eigenvalue decomposition ,in The eigenvector matrix, For matrix The i-th eigenvalue;

[0057] Projecting negative eigenvalues ​​to zero: ,in These are the projected eigenvalues;

[0058] Reconstruct the positive semidefinite matrix:

[0059] .

[0060] Preferably, in the method for estimating the hierarchical Bayesian measurement noise covariance, the forgetting factor ρ takes a value range of [0.9, 0.99].

[0061] Preferably, in the hierarchical Bayesian method for estimating the measurement noise covariance, the exponential forgetting factor is recursively calculated. The steady-state value is It approaches 1 / (1-ρ).

[0062] Preferably, in the method for estimating the hierarchical Bayesian measurement noise covariance, the exponential forgetting factor of the dimensionally gamma distribution layer is recursively calculated. The steady-state value is It approaches 1 / (1-ρ).

[0063] Preferably, in the hierarchical Bayesian method for estimating the noise covariance, in step S6, each diagonal element is restricted to a preset interval. Inside, among which and These are the lower and upper bounds of the noise variance for the i-th measurement dimension, respectively.

[0064] The preset interval For the sensor ranging dimension, take the lower bound of the noise variance. upper limit For the sensor's angular measurement dimension, take the lower bound of the noise variance. upper limit .

[0065] Preferably, in the hierarchical Bayesian method for estimating measurement noise covariance, step S6 involves calculating the measurement noise covariance matrix. The non-diagonal element in the i-th row and j-th column Then, constrain its value to an interval. Internally, this prevents the cumulative numerical explosion of off-diagonal elements during the recursion of the global scale matrix.

[0066] Preferably, in the hierarchical Bayesian method for estimating the measurement noise covariance, in S2, when numerical anomalies occur in the Cholesky triangulation, leading to an unbiased empirical estimation matrix... When computation fails, the empirically estimated matrix is ​​used from the most recently successfully computed matrix. Used as the empirical estimation matrix for the current frame; if the current frame is the first frame and no historical frames are available, then the outer product matrix of the current frame's innovation is used. As the unbiased empirical estimation matrix for the current frame.

[0067] Preferably, in the method for estimating the hierarchical Bayesian measurement noise covariance, the fusion weight coefficient ω in S6 is adaptively determined in the following manner:

[0068]

[0069] in The normalized index (NIS) is the exponential moving average of the innovation normalization statistic, and τ is the expected threshold of the NIS. When the filter is functioning correctly... ω approaches 0.5; when the process noise covariance matrix mismatch leads to When deviating from τ, ω is adaptively increased to 0.8, which increases the unbiased anchoring weight of the global inverse Wieshard distribution layer and suppresses the accumulation of bias introduced by process noise pollution in the dimension-by-dimensional gamma distribution layer.

[0070] Preferably, in the hierarchical Bayesian method for estimating measurement noise covariance, in the first N frames, N≥2, steps S1 to S6 are executed normally and the estimated value of the measurement noise covariance matrix is ​​obtained. Subsequently, the operation in S7 used for measurement updates in the next frame is temporarily suspended, and the measurement update step of the filter in the next frame uses the initial measurement noise covariance matrix. Perform filtering updates, and after the number of recursive frames exceeds N frames, execute step S7 to... The measurement update step is passed to the Kalman filter in the next frame.

[0071] A hierarchical Bayesian device for estimating the noise covariance of a measurement includes:

[0072] The prior construction module is used to construct the measurement noise covariance matrix during the Kalman filter initialization phase. Construct a two-layer hierarchical Bayesian prior structure, where the global layer is an inverse Wissaud distribution. The dimension-wise gamma distribution layer is a gamma distribution that is independent of each measurement dimension. ;

[0073] The unbiased innovation decomposition module is used to obtain the innovation vector at each recursive time k. Obtain the prediction covariance Cholesky factor for the current frame. ,satisfy and the measurement-state cross-covariance matrix Solve the trigonometric equations by substituting forward. Obtain the intermediate variable matrix ,but And explicitly subtract the forecast uncertainty term from the new information extrapolation product. The unbiased empirical estimation matrix is ​​obtained. ;

[0074] The eigenvalue projection repair module is used to repair the unbiased empirical estimation matrix. Perform symmetry transformation, eigenvalue decomposition, project negative eigenvalues ​​to zero, and reconstruct the positive semidefinite matrix. ;

[0075] The global inverse Wissaud update module is used to recursively update the posterior parameters of the global inverse Wissaud distribution layer using an exponential moving average recursive method driven by a forgetting factor ρ. The recursive formula is as follows: and ;when At that time, the full matrix estimate of the posterior expectation of the global inverse Wissaud distribution layer is calculated based on the mean formula of the inverse Wissaud distribution. ;when At that time, the global inverse Wissaud distribution layer posterior expectation full matrix estimation Keeping the values ​​from the previous time step unchanged, and outputting the full matrix estimate of the posterior expectation of the global inverse Wissaud distribution layer for the current frame. ;

[0076] The dimension-wise gamma update module is used to apply the innovation vector for each measurement dimension i. The i-th component The square value The recursive formula for updating the gamma posterior parameters is as follows: and ;like Then calculate the posterior expectation of each dimension. ;

[0077] The hierarchical fusion module is used to receive the full matrix estimate of the posterior expectation of the global inverse Wissaud distribution layer. Construct the measurement noise covariance matrix of the final output of the current frame. ,in For the measurement noise covariance matrix The i-th diagonal element, For the measurement noise covariance matrix The non-diagonal element in the i-th row and j-th column; Full matrix estimation of the posterior expectation of the global inverse Wissaud distribution layer The i-th diagonal element, Full matrix estimation of the posterior expectation of the global inverse Wissaud distribution layer The non-diagonal element in the i-th row and j-th column; for each measurement dimension i, when At that time, the diagonal elements of the full matrix estimate of the posterior expectation of the global inverse Wissaud distribution layer are... With dimensional posterior expectation By performing weighted fusion, we obtain ;when At that time, the corresponding value of the full matrix estimate of the posterior expectation of the global inverse Wissaud distribution layer is directly taken. For all off-diagonal elements, i.e., i ≠ j and i,j∈{1,…,m}, the full matrix estimate of the posterior expectation of the global inverse Wissaud distribution layer is provided. ; and the measurement noise covariance matrix Perform symmetry transformation and restrict each diagonal element to a preset range.

[0078] A computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the method for estimating the hierarchical Bayesian measurement noise covariance.

[0079] The present invention has at least the following beneficial effects:

[0080] This invention achieves independent adaptive modeling of noise in each measurement dimension by constructing a two-layer hierarchical Bayesian prior structure consisting of a global inverse Wissaud distribution layer and a dimension-by-dimensional gamma distribution layer. The global inverse Wissaud distribution layer captures off-diagonal correlation information in the full matrix structure, providing overall structural constraints; the dimension-by-dimensional gamma distribution layer independently models each measurement channel, with each dimension having independent shape and rate parameters, allowing each dimension to independently adapt at different rates. The adaptive weight fusion strategy (ω defaults to 0.5, i.e., equal weight fusion; when the NIS deviates, ω is linearly and continuously adjusted to a maximum of 0.8 within the range of [0.5, 0.8] according to the degree of deviation) fundamentally solves the problem that the estimation of small noise dimensions is overwhelmed by large noise dimensions when the variance of noise in each dimension varies across orders of magnitude in a single inverse Wissaud prior.

[0081] This invention eliminates the ambiguity between process noise and measurement noise through unbiased innovation decomposition. It utilizes the definition of innovation covariance. Explicitly subtract the prediction uncertainty term from the new information extrapolation product. This enables empirical estimation of measurement noise. In the expected sense, it is equal to the true measurement noise covariance matrix. The Cholesky factor trigonometric method enables efficient calculation of the prediction uncertainty term, avoiding the direct inversion of the prediction covariance matrix. This improves computational efficiency and numerical stability while eliminating the systematic bias of process noise contribution being incorrectly attributed to measurement noise.

[0082] This invention provides an unbiased negative definite repair strategy through eigenvalue projection positive semidefinite repair. When subtraction produces negative eigenvalues, eigenvalue decomposition is performed on the unbiased empirical estimation matrix. Projecting negative eigenvalues ​​to zero The reconstructed matrix is ​​a positive semidefinite matrix. It discards only non-physical solutions in the negative definite direction while retaining all effective information in the positive semi-definite direction, and does not introduce system bias in the positive semi-stator space compared to traditional full-matrix backoff or half-decay mixing strategies. When the filter enters steady state, the predicted state covariance matrix... As the value approaches a stable level, the probability of a significant negative eigenvalue approaches zero, and the truncation bias introduced by the repair gradually disappears, ensuring the asymptotic unbiasedness of the estimation.

[0083] This invention achieves an automatic balance between global structural constraints and dimensionally independent adaptive fusion through hierarchical fusion. An adaptive weight fusion strategy (default 0.5:0.5, tilted to 0.5:0.8 when NIS deviates) is employed to fuse the posterior expectations of the global inverse Wissaud distribution layer and the posterior expectations of the dimensionally independent gamma distribution layer. The measurement dimensions undergo adaptive weight fusion, with the fusion weight coefficient ω defaulting to 0.5 (i.e., equal weight fusion), ensuring that the unbiased anchoring of the global inverse Wissaud distribution layer and the independent adaptive capability of the dimension-by-dimensional gamma distribution layer each contribute half. When the exponential moving average of the innovation normalization statistic deviates from the expected threshold, ω adaptively increases to a maximum of 0.8, increasing the unbiased anchoring weight of the global inverse Wissaud distribution layer and effectively suppressing the accumulation of bias in the dimension-by-dimensional gamma distribution layer when the process noise covariance matrix is ​​mismatched. The measurement dimensions are directly taken from the corresponding posterior expectation values ​​of the global inverse Wissaud distribution layer, avoiding bias introduced by unreliable dimension-by-dimensional estimation. All off-diagonal elements are provided by the global inverse Wissaud distribution layer. When the noise levels of each dimension are similar, the estimates of the two layers are close, and the adaptive weight fusion degenerates into a simple average. When the noise in a certain dimension increases abnormally, the dimension-by-dimensional gamma posterior expectation deviates significantly from the diagonal of the posterior expectation of the global inverse Wissaud distribution layer, and the equal-weight fusion automatically tilts towards the deviation direction. This fusion mechanism does not require manual setting of the independent learning rate for each dimension, simplifying the engineering deployment process.

[0084] This invention ensures that the method functions correctly from the first frame by setting initial values ​​for each hyperparameter. Initial values ​​for global degrees of freedom. Strictly greater than m+1, ensuring that the first moment of the posterior expectation of the global inverse Wissaud distribution layer exists from the first frame; the initial shape parameters of the gamma distribution in each measurement dimension. The condition of being greater than 2 ensures that the first moment of the inverse gamma distribution exists before the maintenance correction begins. Initial value of the global scaling matrix. and initial rate parameters The setting makes the posterior expectation estimate of the global inverse Wissaud distribution layer in the first frame equal to The posterior expectation estimate of the dimensional gamma distribution layer in the first frame is equal to That is, the estimated values ​​of both levels in the first frame are equal to the initial measurement noise covariance matrix preset by the engineer according to the sensor specifications. The subsequent two layers evolve independently according to different update mechanisms, without the need for a warm-up phase or additional parameter adjustments.

[0085] This invention achieves an optimal balance between response speed and estimation stability by setting the forgetting factor ρ = 0.95. Global degrees of freedom The steady-state value is 20, and the equivalent memory window is approximately 20 frames, which is long enough to filter out random fluctuations in the innovation sequence, yet short enough to track the time-varying characteristics of the measurement noise covariance matrix. For a typical configuration with a measurement dimension of m=3, the steady-state value of the global degrees of freedom, 20, is strictly greater than 4, and the first moment of the posterior expectation of the global inverse Wissaud distribution layer is always present in steady state; the steady-state value of the dimension-by-dimensional shape parameter, 20, is strictly greater than 2, and the first moment of the posterior expectation of the dimension-by-dimensional gamma distribution layer is always present in steady state. The global inverse Wissaud distribution layer does not freeze but enters the steady-state working region, achieving anchored and non-rigid long-term operational robustness.

[0086] This invention eliminates the asymmetric error introduced by numerical rounding during calculation through a symmetry operation. Symmetry is explicitly performed after hierarchical fusion is completed. This ensures that the final output measurement noise covariance matrix estimate strictly meets the symmetry requirement and conforms to the basic definition of the covariance matrix. This guarantees the stability and reliability of numerical operations such as matrix inversion and Cholesky decomposition involved in subsequent Kalman filter measurement updates, and improves the numerical robustness of the filter during long-term operation.

[0087] This invention prevents abnormal amplification of the Kalman gain by limiting each diagonal element to a preset interval independently set for each measurement dimension. When the diagonal elements of the measurement noise covariance matrix are estimated too small, the inversion operation will cause an abnormal increase in the Kalman gain. The filter over-relys on the current measurement and ignores the prediction information of the dynamic model, and the state estimation may exhibit violent oscillations or even divergence. (Lower limit of sensor ranging dimension) The corresponding standard deviation is 10m, and the upper limit is... The corresponding standard deviation is 100m, and the lower limit of the angular dimension is... The corresponding standard deviation is 0.001 rad, and the upper limit is... The corresponding standard deviation is approximately 0.032 rad. The lower limit ensures numerical stability without being too high to suppress normal noise variance estimation. The corresponding upper limit prevents the filter from completely ignoring measurement information due to overestimation, thus ensuring the filter's basic response capability to the measurement.

[0088] This invention offers controllable computational complexity and strong engineering practicality. All steps employ analytical computation methods, with the computational cost per frame primarily stemming from a single Cholesky triangulation and an eigenvalue decomposition. Under a typical measurement dimension of 3, the incremental computational cost is approximately 5%-10% of that of a basis filter, significantly lower than variational Bayesian methods and Monte Carlo or particle filter methods requiring iterative optimization. Its modular design facilitates implementation on various hardware platforms, including embedded platforms, digital signal processors, graphics processors, and field-programmable gate arrays. It can be used in conjunction with various nonlinear filters such as volumetric Kalman filters, unscented Kalman filters, or extended Kalman filters. Adaptive estimation of measurement noise covariance can be achieved without modifying the core structure of the basis filter, demonstrating excellent integrability and scalability. This meets the stringent requirements of real-time systems such as flight control, autonomous driving, and robot navigation for deterministic processing time boundaries.

[0089] Other advantages, objectives and features of the present invention will become apparent in part from the following description, and in part from those skilled in the art through study and practice of the invention. Attached Figure Description

[0090] Figure 1 This is a flowchart illustrating the hierarchical Bayesian measurement noise covariance adaptive estimation method.

[0091] Figure 2 yes Matrix Trace Compare the bar charts;

[0092] Figure 3 It is dimension-wise Estimation accuracy comparison chart;

[0093] Figure 4 These are comparison images before and after eigenvalue projection repair;

[0094] Figure 5 This is a comparison chart of three repair strategies. Detailed Implementation

[0095] The present invention will now be described in further detail so that those skilled in the art can implement it based on the description.

[0096] The hierarchical Bayesian measurement noise covariance estimation method of this invention is not limited to a specific type of basis filter, and can be used in conjunction with various nonlinear filters such as capillary Kalman filters, unscented Kalman filters, or extended Kalman filters. The following detailed description uses a capillary Kalman filter as the basis filter.

[0097] In the field of Kalman filtering, adaptive estimation of the measurement noise covariance matrix has been a long-standing technical challenge. The optimality of a Kalman filter strictly depends on the precise knowledge of the measurement noise covariance matrix. However, in practical engineering, the measurement noise covariance matrix is ​​often unknown or changes over time. For example, the ranging accuracy of a sensor varies with the target distance, the noise level of an infrared sensor is affected by temperature, and the noise characteristics of each channel in multi-sensor fusion can differ by several orders of magnitude.

[0098] Existing methods suffer from three fundamental difficulties. First, they lack the ability to model anisotropic noise. Variational Bayesian methods use a single inverse Wissaud prior, with all measurement dimensions sharing a uniform degree-of-freedom parameter. When the noise variances of different dimensions vary by orders of magnitude, the estimation of smaller noise dimensions is overwhelmed by larger noise dimensions. Second, there is the confounding bias between process noise and measurement noise. The Sage-Husa method directly uses the innovation outer product as an empirical estimate of measurement noise, but the innovation simultaneously includes prediction errors driven by process noise and measurement noise, leading to a systematic overestimation of measurement noise. Third, there is a lack of unbiased repair strategies for negative definite results. When the prediction covariance is large, subtraction estimation may produce negative eigenvalues. Existing methods, such as simple back-off or half-decay mixing, introduce uncontrollable systematic biases. Although methods using hierarchical Bayesian priors have emerged recently, their hierarchical priors are all used for heavy-tailed distribution modeling or outlier detection, rather than for anisotropic estimation of the noise covariance of each measurement dimension, and none of them include the two core technical features of unbiased innovation decomposition and eigenvalue projection positive semidefinite repair.

[0099] In summary, existing technologies lack an adaptive estimation method for measurement noise covariance that can independently model the noise level of each measurement dimension using hierarchical priors, separate measurement noise from prediction uncertainty through unbiased innovation decomposition, and possess an unbiased negative definite repair strategy.

[0100] This invention addresses the aforementioned technical challenges by providing a hierarchical Bayesian method for estimating the measurement noise covariance. The overall approach involves modeling anisotropic noise through a two-layer hierarchical Bayesian prior structure, eliminating the ambiguity between process noise and measurement noise through unbiased innovation decomposition, providing an unbiased negative definite repair strategy through eigenvalue projection positive semidefinite repair, and achieving an automatic balance between global structure estimation and dimensionally independent estimation through hierarchical fusion. Ultimately, this yields an unbiased and anisotropic online adaptive estimate of the measurement noise covariance matrix. The entire process forms a closed-loop adaptive flow encompassing prior construction, innovation observation, unbiased decomposition, hierarchical fusion, and measurement updates.

[0101] The specific solution of the present invention is as follows:

[0102] Step S1, construct the hierarchical Bayesian prior structure: for the measurement noise covariance matrix A two-layer hierarchical Bayesian prior is constructed, where m is the dimension of the measurement vector. This two-layer hierarchical Bayesian prior includes a global inverse Wissaud distribution layer and a dimension-wise gamma distribution layer. In the global inverse Wissaud distribution layer, the measurement noise covariance matrix... Obeying the distribution of the reverse Saudi Arabia ,in The global scaling matrix is ​​a symmetric positive definite matrix. For global degrees of freedom.

[0103] Dimensionally Gamma Distribution Layer: For each measurement dimension i, the measurement noise covariance matrix diagonal elements corresponding precision Follows a gamma distribution ,in For shape parameters, For rate parameters, For the shape parameter molecule of the gamma distribution in the i-th measurement dimension, it must satisfy... , Let be the rate parameter molecule of the gamma distribution in the i-th measurement dimension.

[0104] The global inverse Wissaud distribution layer is used to capture the full matrix structure of the measurement noise covariance, including off-diagonal correlation information between each measurement dimension. The dimension-wise gamma distribution layer is used to independently model each measurement dimension, with each dimension having independent shape and rate parameters, allowing each measurement dimension to adapt independently at different rates.

[0105] The purpose of step S1 is to address the anisotropic noise modeling problem. A two-layer structure, consisting of a global inverse Wissaud distribution layer and a dimensional gamma distribution layer, is employed. The global inverse Wissaud distribution layer provides overall structural constraints, while the dimensional gamma distribution layer allows each measurement dimension to independently adapt at different rates, thus handling noise differences spanning orders of magnitude, such as distance and angle noise, within the same framework. During the initialization phase, initial values ​​for each hyperparameter and the global degrees of freedom need to be set. The initial value of the global scale matrix Initial shape parameters of gamma distribution in each measurement dimension Initial rate parameters ,in The initial measurement noise covariance matrix is ​​preset by the engineers according to the sensor specifications. Shape parameters The condition must be greater than 2 to guarantee the existence of the first moment of the inverse gamma distribution. The two-layer posterior expectation is integrated using a weighted fusion strategy with fusion weight coefficients ω, the specific value of which is determined by the formula. The adaptive determination is that its value changes linearly and continuously within the range of [0.5, 0.8] as the deviation ratio increases from zero to 1. When the deviation ratio reaches or exceeds 1, ω increases linearly from 0.5 to 0.8. When the deviation ratio reaches or exceeds 1, ω saturates at 0.8 and tilts towards the unbiased anchoring of the global inverse Wieshard distribution layer.

[0106] Step S2, Unbiased innovation decomposition: At each recursive time k, obtain the innovation vector of the current frame. ,in For the current measurement vector, For measurement prediction vectors; measurement prediction vectors The Kalman filter determines the predicted state. and measurement function Calculated. According to the definition of information covariance: ;

[0107] The superscript T indicates matrix transpose. To measure the Jacobian matrix, For predicting the state covariance matrix;

[0108] The purpose of step S2 is to eliminate the ambiguity between process noise and measurement noise. The expectation of the innovation extraproduct equals the prediction uncertainty. Add the actual measurement noise covariance matrix If we directly use the new-innovation outer product as... The estimation of process noise contribution will be incorrectly attributed to measurement noise. Therefore, it is necessary to subtract the prediction uncertainty term from the innovation outer product to obtain an unbiased estimate of the measurement noise. In the global inverse Wissaud distribution layer, the empirical estimate of the measurement noise covariance is separated from the process noise contribution by explicitly subtracting the prediction uncertainty term, thus solving the fundamental problem of systematic overestimation of the measurement noise covariance matrix in the traditional Sage-Husa method.

[0109] The unbiased empirical estimation matrix is ​​calculated using the following sub-steps:

[0110] Step S21: Obtain the prediction covariance Cholesky factor for the current frame. ,satisfy and the measurement-state cross-covariance matrix .

[0111] The purpose of step S21 is to obtain the input data required for subsequent trigonometric solving. (Cholesky factor) It is the result of predicting the square root decomposition of the covariance matrix. This is a lower triangular matrix, and this decomposition guarantees numerical stability. Measurement-state cross-covariance matrix. It describes the degree of correlation between state prediction error and measurement prediction error.

[0112] Step S22, by Solve the trigonometric equations by substituting forward. The intermediate variable matrix is ​​obtained. .

[0113] The derivation process of step S22 is as follows: From ,Will Substitution Solve the trigonometric equations by forward substitution. ,get .then Therefore, the prediction of uncertainties is... This step avoids modifying the predicted covariance matrix. Direct inversion improves computational efficiency and numerical stability.

[0114] Step S23, calculate the unbiased empirical estimation matrix: ;

[0115] The purpose of step S23 is to obtain an unbiased empirical estimate of the measurement noise covariance matrix in the current frame. Subtracting the prediction uncertainty term calculated in step S22 from the innovation outer product eliminates the confusion between process noise and measurement noise, thus obtaining an unbiased empirical estimate of the measurement noise covariance matrix. The hat symbol ^ indicates that the quantity is an estimate, and the subscript emp indicates that the quantity is an empirical estimate. This unbiased separation is crucial for the statistical consistency of the filter because when the measurement noise covariance matrix is ​​overestimated, the filter over-relyes on the dynamic model, leading to biased state estimates.

[0116] Step S3, Eigenvalue projection positive semidefinite repair: The eigenvalue obtained in step S23... The matrix is ​​symmetricized, and then eigenvalue decomposition is performed on the symmetricized matrix. The negative eigenvalues ​​of the decomposed matrix are projected to zero to reconstruct the empirically estimated matrix that satisfies the positive semidefinite requirement. This repair strategy is triggered only when the unbiased empirical estimation matrix exhibits negative eigenvalues, specifically when the prediction covariance is large and the prediction uncertainty term numerically exceeds the innovation outer product. Negative eigenvalues ​​typically occur during the initial convergence phase of the filter, at which point the prediction covariance matrix... The difference is relatively large. Once the filter reaches steady state, the predicted state covariance matrix tends to a stable value. The difference between the innovation outer product and the prediction uncertainty term is always positive definite in the expected sense, and the probability of significant negative eigenvalues ​​approaches zero. Therefore, the correction amount of the eigenvalue projection repair operation on the estimation result approaches zero, essentially degenerating into an identity transformation. Thus, this repair strategy will not be frequently triggered during the long-term operation of the filter and will not affect the real-time computational efficiency of the method.

[0117] Step S3 aims to address the potential for negative eigenvalues ​​in subtraction estimation. This is particularly relevant when the prediction covariance is large, such as during the initial convergence phase of the filter, leading to prediction uncertainty. It may exceed the new interest rate. ,lead to The occurrence of negative eigenvalues ​​is physically impossible. The specific fix is ​​as follows: first, the unbiased empirical estimation matrix is ​​symmetricized to eliminate potential asymmetric errors introduced during computation; then, eigenvalue decomposition is performed. Projecting negative eigenvalues ​​to zero Negative eigenvalues ​​are set to zero while positive eigenvalues ​​remain unchanged, and finally, a positive semidefinite matrix is ​​reconstructed. This strategy discards only non-physical solutions in the negative definite direction, retaining all valid information in the positive semi-definite direction, unlike the traditional full-matrix backoff strategy. Compared to not introducing system bias in the positive semi-stator space, the predicted state covariance matrix is ​​improved after the filter reaches steady state. As the system approaches a stable value, the probability of significant negative eigenvalues ​​approaches zero, and the truncation bias introduced by the repair gradually disappears, achieving asymptotic unbiasedness. Unbiasedness in this context has two meanings: the first is the unbiasedness of subtraction, meaning the subtraction result is, in the expected sense, an unbiased estimate of the measurement noise covariance matrix; the second is the unbiasedness of repair, meaning the eigenvalue projection repair operation does not change the projected value of the original matrix in the positive semi-definite direction, only discarding information in the negative definite direction while retaining the original estimate in the positive semi-definite direction, contrasting with the full matrix backoff strategy, which also introduces bias in the positive semi-definite direction.

[0118] Step S4, Global Inverse Wissaud Distribution Layer Posterior Recursive Update: The hyperparameters of the global inverse Wissaud prior are recursively updated using an exponential moving average driven by the forgetting factor ρ∈(0.9,1.0).

[0119] , ;

[0120] in and These represent the global degrees of freedom and the global scale matrix at the current time k, respectively. and These are the global degrees of freedom and the global scale matrix at the previous time step k-1, respectively.

[0121] like Then, based on the mean formula of the inverse Wissaud distribution, the full matrix estimate of the global inverse Wissaud distribution layer posterior expectation of the current frame is calculated: ;

[0122] like Then the global inverse Wissaud distribution layer posterior expectation full matrix estimate of the current frame Keep the value from the previous moment unchanged;

[0123] Output the full matrix estimate of the global inverse Wissaud distribution posterior expectation of the current frame. .

[0124] Step S4 aims to update the posterior parameters of the global inverse Wissaud distribution layer using the repaired empirical estimate, thereby obtaining a full matrix estimate of the posterior expectation of the global inverse Wissaud distribution layer. A forgetting factor controls the weight distribution between historical information and current observations, and an exponential moving average is used to recursively update the global degrees of freedom and global scale matrix. The forgetting factor ρ∈(0.9,1.0), with a recommended value of 0.95, at which point the global degrees of freedom... It approaches 1 / (1-0.95)=20, and always satisfies this condition. Under the given conditions, the global inverse Wissaud distribution layer remains actively updated at any frame rate. The posterior expectation is calculated using the mean formula of the inverse Wissaud distribution. When the global degrees of freedom are insufficient, the estimate from the previous time step is kept unchanged to avoid unreliable posterior expectation calculations due to insufficient degrees of freedom. Global degrees of freedom The global scale matrix controls the concentration of the inverse Wieshart distribution. It controls the scale of the distribution. The exponential forgetting factor is recursively calculated. The steady-state value is When ρ=0.95 , always maintain Under these conditions, the global inverse Wissaud distribution layer remains actively updated at any frame rate. To provide additional stability for long-term operation, when After reaching steady state, approximately 1 / (1-ρ) frames are passed, and the process continues at... The equivalent sample size is exponentially smoothed, and the global inverse Wissaud distribution layer does not freeze but enters the steady-state working region, achieving long-term robustness of anchoring without rigidity. The global inverse Wissaud distribution layer provides unbiased anchoring by subtracting the prediction uncertainty term, while the dimension-by-dimensional gamma distribution layer provides independent adaptation for each dimension. The two layers run in parallel and are continuously updated. The complementary design achieves long-term stability.

[0125] Step S5, posterior recursive update of the dimensional gamma distribution layer: For each measurement dimension i, use the innovation vector The i-th component The square value Perform recursive updates of the underlying prior parameters:

[0126] , ;

[0127] in and Let be the shape parameter molecules of the i-th measurement dimension at the current time k and the previous time k-1, respectively; and Let be the numerators of the rate parameters for the i-th measurement dimension at the current time k and the previous time k-1, respectively.

[0128] like Then, calculate the posterior expectation estimate for that dimension:

[0129] ;

[0130] like If the first moment of the inverse gamma distribution does not exist, then the posterior expectation estimate of that dimension is... Not calculated.

[0131] The purpose of step S5 is to update the posterior parameters of the dimensional gamma distribution layer using the square of the innovation component, thereby obtaining independent noise variance estimates for each measurement dimension. The gamma posterior parameters are recursively updated independently for each measurement dimension using the square of the innovation component of that dimension. The denominator of the dimensional posterior expectation is the shape parameter minus 2, rather than the shape parameter minus 1, because in the gamma conjugate system, the precision parameter follows a gamma distribution. noise variance It follows an inverse gamma distribution, and its posterior expectation is... To ensure the existence of the first moment, α > 2 must be satisfied. The forgetting factor ρ is also recommended at 0.95. At this point, the dimension-wise shape parameter... Approaching 20, it always satisfies the condition of being greater than 2. The dimensional gamma distribution layer intentionally avoids subtracting the prediction uncertainty term because the forgetting factor provides an effective memory window of approximately 20 frames, allowing the contribution of transient prediction errors to decay over time averaging. Furthermore, the adaptive weight fusion strategy mixes the dimensional gamma estimate with the global inverse Wissat estimate corrected for the prediction uncertainty term subtraction. The unbiased anchoring of the global inverse Wissat distribution layer effectively suppresses the accumulation of bias in the dimensional gamma layer. When the process noise covariance matrix is ​​systematically overestimated, the prediction covariance expands, and the prediction uncertainty term becomes non-negligible. In this case, omitting the subtraction in the dimensional gamma layer would further overestimate the noise variance in that dimension. However, in the two-layer architecture of this invention, the unbiased subtraction of the global inverse Wissat distribution layer provides anchoring suppression for such biases through adaptive weight fusion. During the transient period, i.e., the initial convergence phase of the filter, when the prediction covariance is large, the unbiased subtraction of the global inverse Wissat distribution layer provides anchoring correction for the dimensional gamma layer through adaptive weight fusion, suppressing the accumulation of transient biases.

[0132] Step S6, Layer Fusion: S6 receives the full matrix estimate of the posterior expectation of the global inverse Wissaud distribution layer output from S4. Construct the measurement noise covariance matrix of the final output of the current frame. ,in Represents the measurement noise covariance matrix The i-th diagonal element, Represents the measurement noise covariance matrix The non-diagonal element in the i-th row and j-th column; This represents the full matrix estimation of the posterior expectation of the global inverse Wissaud distribution layer. The i-th diagonal element; This represents the full matrix estimation of the posterior expectation of the global inverse Wissaud distribution layer. The non-diagonal element in the i-th row and j-th column;

[0133] Step S6 receives the output from step S4. Regardless of whether the value is a new value after the update or a value from the previous time step that remains unchanged, it participates in hierarchical fusion.

[0134] For each measurement dimension i, determine The i-th diagonal element: if ,but:

[0135] ;

[0136] Where the fusion weight coefficient ω∈[0.5,0.8];

[0137] like ,but:

[0138] ;

[0139] For all off-diagonal elements, i.e., i ≠ j and i,j∈{1,…,m}, the full matrix estimate of the posterior expectation of the global inverse Wissaud distribution layer is provided: ;

[0140] At this point, the measurement noise covariance matrix... All elements have been determined, and then the entire measurement noise covariance matrix is... Perform symmetry transformation and restrict each diagonal element to a preset range;

[0141] S7, the output measurement noise covariance matrix Measurement update steps for the Kalman filter in the next frame.

[0142] The purpose of step S6 is to fuse the posterior expectation of the global inverse Wissaud distribution layer with the dimensionally progressive posterior expectation to obtain a measurement noise covariance matrix estimate that has both unbiased and anisotropic resolution. For The measurement dimensions undergo adaptive weight fusion (the fusion weight coefficient ω is determined adaptively, specifically ω = 0.5 + 0.3·min(|NIS)). EMA / τ-1|,1), where |NIS EMA / τ-1| represents the deviation ratio of the exponential moving average of the innovation normalization statistic from the expected threshold. When the deviation ratio is zero, ω is 0.5, i.e., equal-weighted fusion; as the deviation ratio increases from zero to 1, ω linearly increases from 0.5 to 0.8; when the deviation ratio reaches or exceeds 1, ω saturates at 0.8. This linear saturation mechanism ensures a smooth transition of fusion weights when there is a slight mismatch in the process noise covariance matrix, avoiding estimation jumps caused by abrupt weight changes. This allows the unbiased anchoring of the global inverse Wissaud distribution layer and the independent adaptive capability of the dimension-by-dimensional gamma distribution layer to each contribute half; for The measurement dimensions are directly taken from the corresponding values ​​of the global inverse Wissaud distribution layer, avoiding bias introduced by unreliable dimension-wise estimation. All off-diagonal elements are provided by the global inverse Wissaud distribution layer, because the dimension-wise gamma distribution layer only estimates diagonal elements and does not provide estimates for off-diagonal elements. This fusion mechanism utilizes the inherent differences in the hierarchical prior structure to achieve automatic balancing, eliminating the need for manually setting independent learning rates for each dimension. When the noise levels of each dimension are similar, the global inverse Wissaud distribution layer estimation and the dimension-wise gamma distribution layer estimation are close, and the fusion degenerates into a simple average; when the noise in a certain dimension is abnormal, the dimension-wise gamma distribution layer estimation automatically deviates from the global inverse Wissaud distribution layer, and the fusion result automatically tilts in the direction of deviation.

[0143] Lower limit of the preset interval and upper limit Each measurement dimension should be set independently. A lower limit is set to prevent abnormal amplification of the Kalman gain due to extremely small variance, while an upper limit is set to match engineering realities where noise levels can differ by several orders of magnitude across dimensions. For example, the lower limit of noise variance for the sensor ranging dimension can be set to 100m. 2 This corresponds to a standard deviation of approximately 10m, with an upper limit that can be set to 10,000m. 2 The corresponding standard deviation is 100m; the angular dimension is taken as the lower limit of the noise variance, which can be set to 10. -6 rad 2 This corresponds to a standard deviation of approximately 0.001 rad, with an upper limit that can be set to 10. -3 rad 2 The corresponding standard deviation is approximately 0.032 rad.

[0144] The forgetting factor ρ=0.95 introduces an inherent scale bias in steady state, affecting the posterior expectation of the global inverse Wissaud distribution layer. Approaching of The posterior expectation of a multi-dimensional gamma distribution layer Approaching of The diagonal elements after fusion with equal weights of ω=0.5 tend to be ,in Let be the true value of the noise variance for the i-th measurement dimension. This bias is conservative, slightly overestimating the measurement noise covariance matrix. To make the filter more cautious in its measurement, in order to improve the process noise covariance matrix Paying a small price to achieve accuracy It provides additional stability protection in the event of mismatch. The engineering rationale for this conservative bias lies in its ability to achieve stability in the event of mismatch. A small trade-off in accuracy. Additional stability guarantees in case of mismatch. Measurement noise covariance matrix. The effect of an 18% amplification on the Kalman gain is that the increased gain denominator reduces the gain by approximately 8%-15%, depending on the specific gain. and The relative magnitude. Reduced gain means the filter trusts dynamic model predictions more than measurements, in the process noise covariance matrix. When correct, this will slightly increase the estimation error, but... The underestimation of gain provided an additional protection margin, partially mitigating the impact. The risk of underestimating the error and causing the filter to over-rely on the incorrect model predictions.

[0145] Compared with existing variational Bayesian R adaptive methods, this invention has the following significant differences and advantages:

[0146] Regarding noise modeling capabilities, existing variational Bayesian R-adaptive methods use a single inverse Wissaud prior, with all measurement dimensions sharing a unified degree-of-freedom parameter. When the noise variances of different dimensions vary by orders of magnitude, the estimation of smaller noise dimensions is overwhelmed by larger noise dimensions. This invention employs a two-layer hierarchical Bayesian prior structure. The global inverse Wissaud distribution layer captures the off-diagonal correlations within the entire matrix structure, while the dimension-by-dimensional gamma distribution layer independently adapts to each measurement channel, enabling each dimension to converge independently at different rates, fundamentally solving the anisotropic noise modeling problem.

[0147] Regarding the accuracy of measurement noise estimation, existing variational Bayesian R-adaptive methods directly use the innovation outer product as an empirical estimate of measurement noise, incorrectly attributing the contribution of process noise to measurement noise. This leads to a systematic overestimation of the measurement noise covariance matrix, causing the filter to over-rely on the dynamic model and resulting in biased state estimation. This invention explicitly subtracts the prediction uncertainty term from the innovation outer product through unbiased innovation decomposition, ensuring that the empirical estimate of measurement noise equals the true measurement noise covariance matrix in an expected sense, thus eliminating the confusion between process noise and measurement noise.

[0148] Regarding numerical stability, existing variational Bayesian R-adaptive methods lack reasonable repair strategies for negative eigenvalues ​​generated by subtraction, and simple backoff introduces uncontrollable system biases. This invention employs eigenvalue projection positive semi-definite repair, discarding only non-physical solutions in the negative definite direction while retaining all effective information in the positive semi-definite direction. It does not introduce biases in the positive semi-stator space and achieves asymptotic unbiasedness after the filter enters steady state.

[0149] Regarding the hierarchical fusion mechanism, while existing hierarchical Bayesian methods have adopted a hierarchical prior structure with inverse Wissat and gamma, their technical purpose is to model heavy-tailed noise distributions or detect outliers, rather than to estimate the anisotropy of the noise covariance matrix. The hierarchical prior of this invention serves the anisotropy estimation of the noise covariance matrix, implicitly coupling the two prior layers through adaptive weight fusion.

[0150] In terms of computational efficiency, the computational overhead per frame of this invention mainly comes from one Cholesky triangulation and one eigenvalue decomposition. Under a typical measurement dimension of 3, the incremental computational cost is about 5%-10% of that of the basis filter, which is far lower than Monte Carlo or particle filter methods.

[0151] In another approach, in the hierarchical Bayesian estimation method for measurement noise covariance, step S6 describes the estimation of the entire measurement noise covariance matrix. Perform symmetry transformation, specifically by transforming the measurement noise covariance matrix. Assigned value .

[0152] In the unbiased innovation decomposition in step S2, the innovation outer product Theoretically, it is a symmetric matrix, predicting uncertainty terms. It must also be a symmetric matrix, and the unbiased empirical estimation matrix obtained by subtracting the two is... Mathematically, symmetry should be maintained. However, in actual floating-point operations, due to the accumulation of numerical rounding errors, This may introduce a small asymmetric bias. If a matrix with asymmetric bias is directly used for subsequent filtering calculations, it may lead to numerical instability in the Kalman filter, or even filter divergence.

[0153] By explicitly performing a symmetry operation after layer blending is complete, i.e. It can eliminate the asymmetric error introduced by numerical rounding during the calculation process, and ensure the estimation of the measurement noise covariance matrix in the final output. Strict symmetry requirements are met. This operation ensures that the measurement noise covariance matrix mathematically conforms to the basic definition of a covariance matrix, thereby guaranteeing the stability and reliability of numerical operations such as matrix inversion and Cholesky decomposition involved in subsequent Kalman filter measurement updates, and improving the numerical robustness of the filter during long-term operation.

[0154] In another approach, the initial values ​​of each hyperparameter in S1 are set in the hierarchical Bayesian method for estimating the noise covariance as follows:

[0155] Initial values ​​of global degrees of freedom ;

[0156] Initial value of the global scale matrix ,in The initial measurement noise covariance matrix;

[0157] Initial shape parameters of gamma distribution in each measurement dimension ;

[0158] Initial rate parameters of gamma distribution in each measurement dimension ,in The initial measurement noise covariance matrix The diagonal elements.

[0159] pass The initial value of the global degrees of freedom is strictly greater than m+1, ensuring that the first moment of the posterior expectation of the global inverse Wissaud distribution layer exists from the first frame. This allows the layer to reliably output the global measurement noise covariance matrix estimate during the initialization phase, without waiting for accumulation over multiple frames. Simultaneously, the initial value of the global scale matrix... This ensures that the initial global posterior expectation, in the absence of observation information, is exactly equal to the initial measurement noise covariance matrix pre-set by the engineer according to the sensor specifications. That is, the prior expectation of the inverse Wissaud distribution is set as This prior setting incorporates engineering experience, ensuring that the filter operates within a reasonable noise level range from startup, thus avoiding filter instability or divergence caused by severe estimation errors in the initial stage.

[0160] The initial shape parameters of the gamma distribution in each measurement dimension The settings satisfy The condition ensures that the first moment of the inverse gamma distribution exists before the maintenance correction begins. Initial rate parameter The setting makes the initial dimension-wise posterior expectation equal to the diagonal elements of the corresponding dimension of the initial measurement noise covariance matrix. This ensures that the initial estimates of the dimensionally progressive gamma distribution layer and the global inverse Wissaud distribution layer remain consistent even without observational information. Thus, the estimates of both layers in the first frame are equal to the initial measurement noise covariance matrix pre-set by the engineers according to the sensor specifications. Subsequently, as the number of recursive frames increases, the two layers evolve independently according to different update mechanisms. The global inverse Wissaud distribution layer updates based on the outer product of the innovations minus the prediction uncertainty term, providing unbiased overall structural constraints. The dimension-by-dimensional gamma distribution layer updates independently based on the squares of the innovation components in each dimension, providing anisotropic resolution. The two layers gradually differentiate based on the same initial values, each performing its own function.

[0161] With the initial shape parameter set to 3 and the forgetting factor ρ taking the recommended value of 0.95, the shape parameter of the dimensional gamma distribution layer is recursively calculated. The steady-state value is 20, which always satisfies the condition of being greater than 2, ensuring the constant existence of the first moment of the dimensionally a posteriori expectation throughout the entire operating cycle. The above initial value setting allows the method of this invention to work normally after startup, without the need for a warm-up stage or additional parameter debugging, simplifying the engineering deployment process and enhancing the practicality of the method.

[0162] In another approach, the eigenvalue projection positive semidefinite repair in S3 of the hierarchical Bayesian measurement noise covariance estimation method specifically includes:

[0163] right Perform eigenvalue decomposition ,in The eigenvector matrix, For matrix The i-th eigenvalue;

[0164] Projecting negative eigenvalues ​​to zero: ,in These are the projected eigenvalues;

[0165] Reconstruct the positive semidefinite matrix: .

[0166] When predicting the covariance matrix When the value is large, such as during the initial convergence phase of the filter, the prediction uncertainty term... It may exceed the new interest rate. This leads to the unbiased empirical estimation matrix. Negative eigenvalues ​​are present. Physically, negative eigenvalues ​​correspond to negative variance or covariance, which contradicts the fundamental definition of a covariance matrix and cannot be directly used in subsequent Kalman filter measurement update steps. The eigenvalue projection positive semi-definite repair method of this invention accurately identifies and sets negative eigenvalues ​​to zero through eigenvalue decomposition, while retaining valid information in the positive semi-definite directions corresponding to all positive eigenvalues. The reconstructed matrix... It meets the positive semidefinite requirement and can be safely used for subsequent global inverse Wissaud distribution layer a posteriori recursive updates.

[0167] Compared with existing technologies, this invention has significant advantages in unbiased repair. Existing methods typically employ a full matrix backoff strategy when encountering negative eigenvalues, such as directly weighting the unbiased empirical estimation matrix with the original innovation outer product, or directly backoffing to the original innovation outer product. These methods also alter the projection values ​​of the original matrix in the positive semi-definite direction, introducing systematic bias in the positive semi-definite subspace. In contrast, the eigenvalue projection strategy of this invention only discards non-physical solutions in the negative definite direction, maintaining the eigenvalues ​​in the positive semi-definite direction. It does not change the projection value of the original matrix, and therefore does not introduce additional systematic bias.

[0168] Once the filter reaches steady state, the predicted state covariance matrix... As the system approaches a stable value, the difference between the innovation outer product and the prediction uncertainty term remains positive definite in the expected sense, and the probability of a significant negative eigenvalue approaches zero. At this point, the correction amount of the eigenvalue projection repair operation to the estimation result approaches zero, essentially degenerating into an identity transformation. The truncation bias introduced by the repair gradually disappears with the increase of the number of recursive frames; therefore, this repair strategy possesses asymptotic unbiasedness. This enables the estimation of this invention to converge to an unbiased estimate of the true measurement noise covariance matrix during long-term operation, ensuring the consistency and stability of the filter in long-term operation.

[0169] In another embodiment, in the method for estimating the hierarchical Bayesian measurement noise covariance, the forgetting factor ρ has a value range of [0.9, 0.99], preferably ρ = 0.95.

[0170] The forgetting factor ρ = 0.95 achieves the optimal balance between response speed and estimation stability. This is based on the exponential moving average recursive formula. When ρ=0.95, the global degrees of freedom The steady-state value is 1 / (1-0.95)=20, and its equivalent memory window is approximately 20 frames. This timescale is long enough to filter out random fluctuations in the innovation sequence and ensure the statistical stability of the estimation, while being short enough to track the time-varying characteristics of the measurement noise covariance matrix. When the forgetting factor is below 0.9, the equivalent memory window is less than 10 frames, and the estimation result is too sensitive to noise fluctuations, easily causing drastic jumps due to anomalous innovations in a single frame. When the forgetting factor is above 0.99, the equivalent memory window exceeds 100 frames, resulting in a slow response speed that cannot effectively track the rapid time-varying nature of the measurement noise covariance, such as scenarios where the sensor ranging accuracy changes drastically with the target distance or noise level drift is caused by sudden changes in sensor temperature.

[0171] This value also possesses important mathematical guarantees. For the typical configuration with a measurement dimension of m=3, the steady-state value of the global degrees of freedom, 20, is strictly greater than m+1=4, and the first moment of the posterior expectation of the global inverse Wieshard distribution layer remains constant in steady state. The shape parameters of the dimensionally progressive gamma distribution layer are recursively derived. The steady-state value is also 20, which is strictly greater than 2, ensuring that the first moment of the dimensionally successive posterior expectation always exists in the steady state. This makes the posterior expectations of the two priors always computable after entering the steady-state working region, without the need for special processing.

[0172] In practical engineering applications, ρ=0.95 is the default configuration that requires no tuning. This value has been validated across various target motion modes and multi-dimensional anisotropic noise scenarios, and can be directly used as the recommended parameter for engineering deployment, simplifying the parameter tuning work for engineers and enhancing the practicality and scalability of the method. The intentional omission of prediction uncertainty term subtraction in the dimensional gamma distribution layer is based on the engineering fact that the forgetting factor ρ=0.95 provides an effective memory window of approximately 20 frames, effectively attenuating the transient prediction error contribution in the innovation component during time averaging.

[0173] In another approach, the method for estimating the hierarchical Bayesian measurement noise covariance involves recursion of the exponential forgetting factor. The steady-state value is It approaches 1 / (1-ρ).

[0174] Recurrence formula The steady-state value 1 / (1-ρ) has a clear physical meaning, namely, the equivalent sample size. When the forgetting factor ρ=0.95, the global degrees of freedom The steady-state value is 20, meaning that in the steady-state operating region, the posterior estimate of the global inverse Wissaud distribution layer is equivalent to statistical inference based on approximately 20 frames of valid samples. This equivalent sample size provides sufficient statistical significance to suppress the impact of random fluctuations on the estimate while maintaining a moderate degree of adaptability to track the time-varying statistical properties of noise.

[0175] This steady-state value also provides important mathematical guarantees. Global degrees of freedom in steady state. When ρ = 0.95 For a typical measurement dimension configuration of m=3, 20 is strictly greater than m+1=4, and the first moment of the posterior expectation of the global inverse Wieshard distribution layer always exists in steady state. This mathematical property makes the statement in step S4, "If...", correct. The conditional judgment of "then calculate, otherwise keep the previous time step unchanged" is always true in steady state. The global inverse Wissaud distribution layer enters the continuously active and updated steady-state working area without freezing, which ensures that the estimation of the global measurement noise covariance matrix is ​​continuously optimized in long-term operation.

[0176] Global degrees of freedom The steady-state value also serves to guarantee the existence of moments in the stepwise gamma distribution layer. Since the forgetting factors of both layers take the same ρ value, the exponential forgetting factor of the stepwise gamma distribution layer is recursively calculated. The steady-state value is 1 / (1-ρ), which is 20, strictly greater than 2, ensuring that the first moment of the posterior expectation of the dimension-wise gamma distribution layer always exists in steady state. This consistency between the steady-state value of the global degrees of freedom and the steady-state value of the dimension-wise shape parameter reflects the inherent coordination of the hierarchical Bayesian framework, providing a dual mathematical guarantee for the long-term stable integration of the global layer and the dimension-wise layer.

[0177] This steady-state value ensures that the global inverse Wissaud distribution layer remains actively updated rather than frozen at any frame rate, achieving robust long-term operation with non-rigid anchoring. The global inverse Wissaud distribution layer provides unbiased anchoring through the subtraction of prediction uncertainty terms, while the dimensional gamma distribution layer provides independent adaptation for each dimension. The two layers operate in parallel and are continuously updated. Their complementary design relies on the existence of moments under steady-state conditions to ensure long-term maintenance.

[0178] In another approach, the method for estimating the hierarchical Bayesian measurement noise covariance involves recursively applying the exponential forgetting factor of the dimensional gamma distribution layer. The steady-state value is It approaches 1 / (1-ρ).

[0179] Recurrence formula The steady-state value 1 / (1-ρ) has a clear physical meaning, namely, the equivalent sample size of the dimension-wise gamma distribution layer. When the forgetting factor ρ=0.95, the dimension-wise shape parameter The steady-state value is 20, meaning that in the steady-state operating region, the posterior estimation of the dimension-wise gamma distribution layer is equivalent to statistical inference based on approximately 20 effective samples. This equivalent sample size provides sufficient statistical significance to suppress the impact of random fluctuations in the innovation components of each dimension on the estimation, while maintaining a moderate degree of adaptability to track the independent time-varying noise statistical characteristics of each measurement channel.

[0180] This steady-state value also provides a mathematical guarantee for the existence of the first moment of the posterior expectation in each dimension. The posterior shape parameter in steady state... When ρ=0.95 20 is strictly greater than 2, and the first moment of the posterior expectation of the layer of the dimensional gamma distribution always exists in steady state. This mathematical property makes step S5... The conditional judgment that calculates if the condition is true and not calculated otherwise is always true in steady state. The dimension-wise gamma distribution layer enters the continuously active and updated steady-state working region without interruption, ensuring that the dimension-wise measurement noise variance estimation is continuously optimized in long-term operation.

[0181] The consistency between this steady-state value and the steady-state value of the global degrees of freedom reflects the inherent coordination of the hierarchical Bayesian framework. Since the forgetting factors in both layers take the same ρ value, the global degrees of freedom... The steady-state values ​​are both 1 / (1-ρ), i.e., 20. The steady-state value of the global degrees of freedom, 20, is strictly greater than m+1, and the steady-state value of the dimension-by-dimensional shape parameter, 20, is strictly greater than 2. The first-order moments of the posterior expectations of both layers of priors are always present in steady state. This dual mathematical guarantee provides a solid foundation for the long-term stable integration of the global layer and the dimension-by-dimensional layer, ensuring that the adaptive weight fusion strategy remains effective in steady state. The global inverse Wissaud distribution layer provides unbiased anchoring, and the dimension-by-dimensional gamma distribution layer provides anisotropic resolution. The two layers operate continuously in parallel, complementing and cooperating with each other.

[0182] This steady-state value also ensures that the dimensional gamma distribution layer remains actively updated rather than interrupted at any number of frames, achieving continuous adaptive estimation of the dimensional measurement channel noise variance. The dimensional gamma distribution layer intentionally avoids subtracting the prediction uncertainty term, precisely because the approximately 20-frame effective memory window provided by the steady-state value of the dimensional shape parameter 20 allows the transient prediction error contribution in the innovation component to be effectively attenuated in time averaging, while simultaneously reducing the process noise covariance matrix... When systematically overestimated, the unbiased subtraction method of the global inverse Wissaud distribution layer provides anchoring suppression for such biases through adaptive weight fusion. This steady-state value ensures that the dimension-by-dimensional gamma distribution layer converges independently in each measurement dimension. Combined with the overall structural constraints provided by the global layer, each dimension adapts independently at different rates, fundamentally realizing online real-time estimation of anisotropic noise.

[0183] In another approach, in the hierarchical Bayesian method for estimating the noise covariance, step S6 restricts each diagonal element to a preset interval. Inside, among which and These are the lower and upper bounds of the noise variance for the i-th measurement dimension, respectively; the lower bound of the preset interval. Take 0.01-0.1 times the true value of the noise variance in this measurement dimension, with the upper limit being... Take 2-10 times the true value of the noise variance of this measurement dimension;

[0184] Lower limit of the preset interval and upper limit Each measurement dimension can be set independently. A lower limit is set to prevent abnormal amplification of the Kalman gain due to extremely small variance, while an upper limit is set to match engineering realities where noise levels can differ by several orders of magnitude across dimensions. For example, the lower limit for the sensor's ranging dimension can be set to 100m. 2 This corresponds to a standard deviation of approximately 10m, with an upper limit that can be set to 10,000m. 2 The corresponding standard deviation is 100m; the lower limit of the angular dimension can be set to 10. -6 rad 2 This corresponds to a standard deviation of approximately 0.001 rad, with an upper limit that can be set to 10. -3 rad 2The corresponding standard deviation is approximately 0.032 rad.

[0185] The beneficial effects of the above settings are reflected in the following aspects:

[0186] First, prevent abnormal amplification of the Kalman gain. In the measurement and update steps of the Kalman filter, the formula for calculating the Kalman gain is: When the measurement noise covariance matrix When the diagonal elements are estimated to be too small, the matrix... The inverse operation will cause an abnormal increase in Kalman gain, leading the filter to over-rely on the current measurement and ignore the prediction information of the dynamic model. In extreme cases, an abnormally large gain may cause the filter to become overly sensitive to measurement noise, resulting in violent oscillations or even divergence in state estimation. By setting a lower bound on the noise variance, the measurement noise covariance matrix is ​​ensured... All diagonal elements are greater than or equal to This fundamentally avoids the problem of abnormal amplification of Kalman gain caused by excessively small variance estimation.

[0187] Second, each measurement dimension has its own independently set intervals to match actual engineering conditions. The noise levels of different measurement channels may differ by several orders of magnitude. Taking sensor target tracking as an example, the typical noise variance for the distance dimension is 2500m. 2 The corresponding standard deviation is 50m, while the typical noise variance for the angular dimension is 0.0007rad. 2 The corresponding standard deviation is approximately 0.026 rad. If a uniform lower bound is used for all dimensions, the lower bound will either be too large, causing the estimation sensitivity of small noise dimensions to be lost, or too small, failing to effectively prevent abnormal gain amplification in large noise dimensions. This is addressed by setting the lower bound independently for each measurement dimension. and This allows the noise variance range of each measurement channel to be adapted to the actual noise level of that channel, thus improving the engineering applicability of the method.

[0188] Third, the engineering rationality of the specific values. The lower limit of the sensor's ranging dimension can be set to 100m. 2 This corresponds to a standard deviation of 10m, which is typically higher than the limit of accuracy for practical ranging sensors, effectively preventing abnormal gain amplification. The lower limit of the angular dimension can be set to 10. -6 rad 2 This corresponds to a standard deviation of 0.001 rad, or approximately 0.057°, which is also higher than the limit accuracy of most practical angle sensors. These lower limits ensure numerical stability without being so high as to suppress normal noise variance estimation. The upper limit of the sensor's ranging dimension can be set to 10000m. 2 This corresponds to a standard deviation of 100m. The upper limit of the angular dimension can be set to 10. -3 rad 2The corresponding standard deviation is approximately 0.032 rad. Upper limit. The corresponding settings prevent the filter from completely ignoring the measurement information due to overestimation, thus ensuring the filter's basic response capability to the measurement.

[0189] In another approach, in the hierarchical Bayesian estimation method for measurement noise covariance, step S6 involves calculating the measurement noise covariance matrix. The non-diagonal element in the i-th row and j-th column Then, constrain its value to an interval. Internally, this prevents the cumulative numerical explosion of off-diagonal elements during the recursion of the global scale matrix.

[0190] This invention uses interval constraints on the off-diagonal elements of S6 to define the off-diagonal elements of the measurement noise covariance matrix. Restricted to diagonal elements and Within the defined range of the Cauchy-Schwarz inequality, that is, guaranteed This constraint ensures that the covariance matrix strictly satisfies the positive definiteness requirement mathematically, avoiding filter instability or even divergence caused by the accumulation of numerical errors during the global scaling matrix recursion, which leads to off-diagonal elements exceeding physically reasonable ranges. Furthermore, this interval constraint only pulls off-diagonal elements back to the boundary when they exceed the limit, without altering the overall structure of the covariance matrix, thus guaranteeing the consistency of the estimation results and improving the algorithm's numerical robustness over long periods.

[0191] In another approach, in the hierarchical Bayesian method for estimating the noise covariance, in S2, when numerical anomalies occur in the Cholesky triangulation, leading to an unbiased empirical estimation matrix... When computation fails, the empirically estimated matrix is ​​used from the most recently successfully computed matrix. Used as the empirical estimation matrix for the current frame; if the current frame is the first frame and no historical frames are available, then the outer product matrix of the current frame's innovation is used. As the unbiased empirical estimation matrix for the current frame.

[0192] In another approach, in the method for estimating the hierarchical Bayesian measurement noise covariance, the fusion weight coefficient ω in S6 is adaptively determined in the following manner:

[0193]

[0194] in The normalized index (NIS) is the exponential moving average of the innovation normalization statistic, and τ is the expected threshold of the NIS. When the filter is functioning correctly... ω approaches 0.5; when the process noise covariance matrix mismatch leads to When deviating from τ, ω is adaptively increased to 0.8, which increases the unbiased anchoring weight of the global inverse Wieshard distribution layer and suppresses the accumulation of bias introduced by process noise pollution in the dimension-by-dimensional gamma distribution layer.

[0195] This invention employs an adaptive determination scheme for the fusion weight coefficient ω in S6, enabling the weight allocation of hierarchical fusion to be dynamically adjusted according to the actual operating state of the filter. When a mismatch occurs in the process noise covariance matrix, the exponential moving average of the innovation normalized statistic is used. It will deviate from its expected threshold τ. At this point, ω adaptively increases, and the unbiased anchoring weight of the global inverse Wissaud distribution layer increases accordingly, thereby effectively suppressing the accumulation of bias in the dimension-by-dimensional gamma distribution layer caused by process noise contamination. When the filter is working normally, When ω approaches 0.5, the fusion degenerates into equal-weight fusion, maintaining a balanced contribution between the global inverse Wissaud distribution layer and the dimensional gamma distribution layer. This adaptive mechanism requires no manual intervention and automatically tilts towards the unbiased anchoring direction of the global inverse Wissaud distribution layer when the process noise covariance matrix mismatches, enhancing the filter's robustness and adaptability in complex environments.

[0196] In another approach, in the hierarchical Bayesian method for estimating measurement noise covariance, steps S1 to S6 are executed normally in the first N frames (N≥2) and the estimated value of the measurement noise covariance matrix is ​​obtained. Subsequently, the operation in S7 used for measurement updates in the next frame is temporarily suspended, and the measurement update step of the filter in the next frame uses the initial measurement noise covariance matrix. Perform filtering updates, and after the number of recursive frames exceeds N frames, execute step S7 to... The measurement update step is passed to the Kalman filter in the next frame.

[0197] This invention employs a preheating phase design for the first N frames to temporarily delay the adaptive estimation of the measurement noise covariance matrix during the initial convergence of the filter. Measurement updates are still performed using the initial measurement noise covariance matrix. The filter is updated, and steps S1 to S6 are executed in their entirety to accumulate the exponential moving average states of each hyperparameter. This scheme solves the problem of predicting the covariance matrix during the initial convergence phase of the filter. When the measurement noise covariance matrix generated by adaptive estimation is large and the innovation statistics are not yet stable, The problem of insufficient reliability is addressed by updating each hyperparameter recursively to approach the true value during the warm-up period. The adaptive estimation results are then officially activated only after the filter enters the stable operating region. This avoids the initial transient noise misleading the adaptive estimation of the measurement noise covariance, and prevents the estimation bias in the initial stage of the filter from forming a positive feedback amplification through measurement updates. This improves the stability and estimation accuracy of the filter during the startup phase, eliminating the need for complex manual initialization and debugging in engineering deployments, and demonstrating good engineering practicality and robustness.

[0198] A hierarchical Bayesian device for estimating the noise covariance of a measurement includes:

[0199] The prior construction module is used to construct the measurement noise covariance matrix during the Kalman filter initialization phase. Construct a two-layer hierarchical Bayesian prior structure, where the global layer is an inverse Wissaud distribution. The dimension-wise gamma distribution layer is a gamma distribution that is independent of each measurement dimension. ;

[0200] The unbiased innovation decomposition module is used to obtain the innovation vector at each recursive time step. Obtain the prediction covariance Cholesky factor for the current frame. ,satisfy and the measurement-state cross-covariance matrix Solve the trigonometric equations by substituting forward. Obtain the intermediate variable matrix ,but And explicitly subtract the forecast uncertainty term from the new information extrapolation product. The unbiased empirical estimation matrix is ​​obtained. ;

[0201] The eigenvalue projection repair module is used to repair the unbiased empirical estimation matrix. Perform symmetry transformation, eigenvalue decomposition, project negative eigenvalues ​​to zero, and reconstruct the positive semidefinite matrix. ;

[0202] The global inverse Wissaud update module is used to recursively update the posterior parameters of the global inverse Wissaud distribution layer using an exponential moving average recursive method driven by a forgetting factor ρ. The recursive formula is as follows: and ;when At that time, the full matrix estimate of the posterior expectation of the global inverse Wissaud distribution layer is calculated based on the mean formula of the inverse Wissaud distribution. ;when At that time, the global inverse Wissaud distribution layer posterior expectation full matrix estimation Keeping the values ​​from the previous time step unchanged, and outputting the full matrix estimate of the posterior expectation of the global inverse Wissaud distribution layer for the current frame. ;

[0203] The dimension-wise gamma update module is used to apply the innovation vector to each measurement dimension i. The i-th component The square value The recursive formula for updating the gamma posterior parameters is as follows: and ;like Then calculate the posterior expectation of each dimension. ;

[0204] The hierarchical fusion module is used to receive the full matrix estimate of the posterior expectation of the global inverse Wissaud distribution layer. Construct the measurement noise covariance matrix of the final output of the current frame. ,in For the measurement noise covariance matrix The i-th diagonal element, For the measurement noise covariance matrix The non-diagonal element in the i-th row and j-th column; Full matrix estimation of the posterior expectation of the global inverse Wissaud distribution layer The i-th diagonal element, Full matrix estimation of the posterior expectation of the global inverse Wissaud distribution layer The non-diagonal element in the i-th row and j-th column; for each measurement dimension i, when At that time, the diagonal elements of the full matrix estimate of the posterior expectation of the global inverse Wissaud distribution layer are... With dimensional posterior expectation Perform adaptive weight fusion to obtain ;when At that time, the corresponding value of the full matrix estimate of the posterior expectation of the global inverse Wissaud distribution layer is directly taken. For all off-diagonal elements, i.e., i ≠ j and i, j ∈ {1, …, m}, the full matrix estimate of the posterior expectation of the global inverse Wissaud distribution layer is provided. ; and the measurement noise covariance matrix Perform symmetry transformation and restrict each diagonal element to a preset range.

[0205] The hierarchical Bayesian measurement noise covariance estimation device described in this invention implements the hierarchical Bayesian measurement noise covariance estimation method in the form of functional modules, including a priori construction module, an unbiased innovation decomposition module, an eigenvalue projection repair module, a global inverse Wissaud update module, a dimension-by-dimensional gamma update module, and a hierarchical fusion module. The overall beneficial effect of this device lies in modularizing the complete process of online adaptive estimation of the measurement noise covariance matrix. Each module undertakes its corresponding function and is sequentially connected through data flow, realizing a complete closed-loop processing flow from innovation input to the output of the estimated value of the measurement noise covariance matrix. This allows the invention to be deployed in various real-time systems such as navigation and positioning, target tracking, and multi-sensor fusion in the form of software, hardware, or firmware, serving the fields of industrial artificial intelligence network system services and industrial Internet of Things information perception and sensing technology services.

[0206] The functional division and data flow of each module in this device are as follows: The prior construction module is responsible for the measurement noise covariance matrix during the Kalman filter initialization phase. A two-level hierarchical Bayesian prior structure is constructed to provide initial values ​​for the hyperparameters of subsequent modules. The unbiased innovation decomposition module obtains the innovation vector at each recursive time step. Obtain the predicted covariance Cholesky factor Measurement-state cross-covariance matrix Solve the trigonometric equations by substituting forward. Obtain the intermediate variable matrix Calculate the prediction uncertainty term The prediction uncertainty term is then explicitly subtracted from the outer product of the new information to obtain the unbiased empirical estimation matrix. This module uses an analytical calculation method, which eliminates the need for iterative optimization and results in high computational efficiency.

[0207] The eigenvalue projection repair module performs symmetry processing on the unbiased empirical estimation matrix, performs eigenvalue decomposition, projects negative eigenvalues ​​to zero, and reconstructs the positive semi-definite matrix. This ensures the positive semi-qualitative nature of the input data for subsequent modules. The global inverse Wissaud update module uses the forgetting factor ρ to drive the posterior parameters of the global inverse Wissaud distribution layer to be updated recursively using an exponential moving average. The recursive formula is as follows: and And based on the current global degrees of freedom Whether it is greater than m+1 determines the calculation of the new global layer posterior expectation full matrix estimate. Keeping the values ​​from the previous time step unchanged, output the full matrix estimate of the global posterior expectation for the current frame. The dimension-wise gamma update module uses an innovation vector for each measurement dimension i. The i-th component The square value The recursive formula for updating the gamma posterior parameters is as follows: and ,when Time-wise computation of dimensional posterior expectation The hierarchical fusion module receives the full matrix estimate of the global layer's posterior expectation. and dimensional posterior expectation For each measurement dimension i, when When and Perform adaptive weight fusion when Take directly off-diagonal elements Provided by the global layer, and performing symmetry and diagonal element pruning, the final measurement noise covariance matrix estimate is output. .

[0208] The modular architecture features clear data dependencies and low coupling between modules, facilitating implementation and independent maintenance on different hardware platforms. Each module employs analytical computation, resulting in short and deterministic single-frame processing times, meeting the stringent requirements of real-time systems for deterministic processing time boundaries.

[0209] A computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the method for estimating the hierarchical Bayesian measurement noise covariance.

[0210] The following uses a specific discrete-time nonlinear state-space model as an example to illustrate the implementation of the present invention.

[0211] Consider the following discrete-time nonlinear state-space model:

[0212] ,

[0213] ,

[0214] in Let n be the state vector (where n is the dimension of the state vector). is the measurement vector (m is the dimension of the measurement vector). This is the state transition function. For measurement functions, and Let Q and R be the independent process noise and measurement noise, respectively, and Q and R be the corresponding covariance matrices. This invention focuses on problems where R is unknown or time-varying, while Q can be known or handled by an independent Q-adaptive method.

[0215] Taking target tracking in three-dimensional space as an example: State vector n=6, where px ,p y ,p z v represents the position component of the target in three-dimensional space. x ,v y ,v z Let n be the velocity component of the target in three-dimensional space, and n be the dimension of the state vector.

[0216] The measurement vector is in sensor spherical coordinates: m=3.

[0217] Where r is the distance and θ is the azimuth angle. Let m be the pitch angle and m be the dimension of the measurement vector.

[0218] Measurement function Nonlinear coordinate transformation: ;

[0219] This invention uses a capacitive Kalman filter as the basis filter. The capacitive Kalman filter generates 2n equally weighted capacitive points using a third-order spherical-radial capacitive rule: , , i=1,…,n;

[0220] in Let represent the i-th column of the n-dimensional identity matrix. The weights of each volume point are: .

[0221] At each recursive time k, the square root capacitive Kalman filter performs the following steps:

[0222] (a) Time update: based on the Cholesky factor from the previous time step. Generate volume points through the state transition function. Propagation, calculation of predicted state and the predicted covariance Cholesky factor (Through QR decomposition).

[0223] (b) Measurement update: from predicted state and the predicted covariance Cholesky factor Generate volume points through measurement functions. Propagation, computational measurement and prediction Cholesky factor for new information covariance Covariance Then calculate the Kalman gain. And update the state estimate Cholesky factor for sum and covariance .

[0224] Initialization phase:

[0225] Global Reversal Saudi Priority:

[0226] The initial values ​​for the global degrees of freedom are set as follows: ;

[0227] The initial value of the global scaling matrix is ​​set as follows: ;

[0228] in The initial measurement noise covariance matrix (set by the engineer according to the sensor specifications).

[0229] Dimensionally successive Gamma priors:

[0230] The initial shape parameter molecule of the gamma distribution for each measurement dimension is set as follows: , i=1,…,m;

[0231] The initial rate parameter numerator for the gamma distribution in each measurement dimension is set as follows:

[0232] ;

[0233] The forgetting factor was set to ρ = 0.95;

[0234] The lower limit of noise variance is set separately for each measurement dimension: for example, the lower limit of noise variance for the sensor ranging dimension is... The lower bound of the noise variance for the angular dimension is The upper limit of noise variance is set accordingly; for example, the upper limit of noise variance in the sensor ranging dimension is... The upper limit of the noise variance in the angular dimension is This is to match engineering realities where the noise levels in various dimensions can differ by several orders of magnitude.

[0235] Automatic warm-up phase: The first N=5 frames of the recursion are the automatic warm-up period. During this period, the two prior layers continue to perform hyperparameter recursive updates to accumulate the exponential moving average state, but the output of step S6 is not updated. Applied to measurement updates, the measurement updates still use the initial measurement noise covariance matrix. After the warm-up period, it will automatically switch to normal working mode.

[0236] Execution steps per frame:

[0237] Step 1: Obtain the latest information

[0238] New information obtained after measurement update: ;

[0239] Predicting the Cholesky factor for covariance: ;

[0240] Measurement-state cross-covariance matrix: .

[0241] Step 2: Unbiased Information Decomposition

[0242] Calculate the new outside product: ;

[0243] Solve the trigonometric equations by substituting forward. The intermediate variable matrix is ​​obtained. Use a lower triangular solver to solve the problem;

[0244] The uncertainty term for prediction is: ;

[0245] Calculate the unbiased empirical estimation matrix: .

[0246] Step 3: Eigenvalue projection positive semidefinite repair

[0247] right After symmetry transformation, eigenvalue decomposition is performed. If all eigenvalues ​​are greater than or equal to zero, then the eigenvalues ​​are directly used. Otherwise, prune the negative eigenvalues ​​to zero and reconstruct the positive semidefinite matrix. .

[0248] Step 4: Global Reverse Saudi Arabia Distribution Layer Update

[0249] Global degrees of freedom recursively updated: ;

[0250] Steady-state value: ;

[0251] Always satisfied conditions.

[0252] Global scale matrix recursive update: ;

[0253] Global layer posterior expectation full matrix estimation: .

[0254] Step 5: Update the dimensional gamma distribution layer

[0255] For each measurement dimension i:

[0256] Shape parameter molecular recursive update: ;

[0257] The steady-state value is 20, and it always satisfies the condition of being greater than 2.

[0258] Rate parameter molecular recursive update: ;

[0259] Dimensionally successive posterior expectation estimation: ;

[0260] And cut to the preset area .

[0261] Step 6: Hierarchical Fusion. The diagonal elements are obtained by weighted fusion of the global posterior expectation and the dimension-wise posterior expectation:

[0262] ;

[0263] Off-diagonal elements are provided by the global layer posterior expectation full matrix estimate:

[0264] i ≠ j;

[0265] Symmetry is performed after fusion:

[0266] ;

[0267] Crop each diagonal element to the preset area .

[0268] Step 7: Output the estimated value of the measurement noise covariance matrix: ;

[0269] Will Measurement update steps for the Kalman filter in the next frame.

[0270] This invention focuses on realizing the measurement noise covariance matrix using hierarchical Bayesian prior structures and unbiased innovation decomposition. The core innovation of adaptive estimation lies in:

[0271] (i) The two-layer hierarchical prior consists of a global inverse Wissaud distribution layer and a dimension-by-dimensional gamma distribution layer to achieve anisotropic noise modeling;

[0272] (ii) through Subtraction achieves unbiased decomposition of new information;

[0273] (iii) Eigenvalue projection provides an unbiased guarantee for negative definite repair.

[0274] Figure 1 This is a flowchart illustrating the hierarchical Bayesian measurement noise covariance adaptive estimation method; it supports the entire process from S1 to S6.

[0275] First, S1 constructs a two-layer hierarchical Bayesian prior: the global inverse Wieshard distribution layer is... Capture the correlation of the entire matrix structure, and perform a dimension-wise gamma distribution layer. Each measurement dimension is modeled independently.

[0276] Then, the measurement update step is performed to obtain the new information vector: ;

[0277] S2 performs unbiased innovation decomposition: subtracts the prediction uncertainty term from the outer product of the innovations. The unbiased empirical estimation matrix is ​​obtained as follows: ;

[0278] S3 performs eigenvalue projection positive semidefinite repair: after symmetric eigenvalue decomposition, negative eigenvalues ​​are projected to zero.

[0279] ;

[0280] Reconstructing the positive semidefinite matrix .

[0281] The next step involves a parallel, two-layer posterior recursion:

[0282] The S4 global inverse Wittner layer performs a forgetting factor-driven hyperparameter moving average update, and calculates the posterior expectation after determining the degrees of freedom. ;

[0283] The S5 dimensional gamma layer independently performs hyperparameter updates driven by the square of the innovation component in each dimension, and calculates the dimensional posterior expectation after determining the shape parameter conditions. .

[0284] The two-layer posterior results are fused at the S6 level: diagonal elements are weighted adaptively. The two-layer expectation is integrated, and the off-diagonal elements are provided by the global layer. After symmetry, they are pruned dimension by dimension. Final output Next frame measurement update step.

[0285] Figure 2 yes Matrix Trace Comparing the histograms, the results support the S2 unbiased innovation decomposition and the S3 eigenvalue projection positive semi-definite repair. The truth benchmark is: ;

[0286] The green dashed line represents the truth baseline. The trace estimation results of the three methods are as follows:

[0287] direct The trace convergence value of the method is 3575.4 (orange bar, relative error 43%). The bias is systematically overestimated due to process noise confounding.

[0288] Unbiased subtraction method of the present invention The convergence value of the trace is 2550.5 (blue bar, relative error of only 2%, close to the true value).

[0289] Truth value Trace 2502.0 (green column).

[0290] The figure includes annotations indicating an improvement in accuracy of over 20 times. The figure uses a Nature colorblind-friendly color palette, with blue representing the method of this invention, orange representing the direct comparison method, and green representing the truth benchmark.

[0291] Figure 3 It is dimension-wise Estimation accuracy comparison chart; technical effect of supporting S1-level Bayesian prior processing of anisotropic noise.

[0292] The horizontal axis is divided into distance dimensions according to the measurement dimensions ( ), azimuth dimension ( ) and pitch angle dimension ( The noise variance spans six orders of magnitude across the three dimensions, with the vertical axis representing each dimension. The relative error percentage is estimated. The figure uses a Nature color-blind friendly palette; blue represents the method of this invention, and orange represents the variational Bayesian R-adaptive method (single inverse Wissaud prior). The variational Bayesian R-adaptive method suffers errors as high as 324% and 298% in the angle dimension (the single prior collapses in the small noise dimension), while the method of this invention has corresponding errors of only 23% and 26%, representing an accuracy improvement of over 10 times (14 times and 11 times), validating the independent adaptive modeling capability of the two-layer hierarchical prior structure for anisotropic noise across orders of magnitude. In the distance dimension, the error of the method of this invention (12%) is also superior to the 18% of the variational Bayesian R-adaptive method.

[0293] Figure 4 These are before-and-after comparison images of eigenvalue projection restoration; used to measure dimensionality. For example, it supports positive semi-definite repair of S3 eigenvalue projection.

[0294] The image uses the Nature colorblind-friendly color palette, with blue representing positive eigenvalue directions, red representing negative eigenvalue anomalies, and gray representing discard directions projected to zero. Left panel. Figure 4 (a) Before restoration The eigenvalue decomposition results.

[0295] The eigenvalues ​​of the positive semi-definite direction (blue) are: , ;

[0296] Negative eigenvalues ​​(marked in red) are: ;

[0297] Labeling the physical impossibility of negative eigenvalues ​​and their causes: Predicting covariance An excessively large value will result in a negative semidefinite result from subtraction.

[0298] Right panel Figure 4 (b) Display the repair results after the eigenvalue projection operation. Projection operation: ;

[0299] The result after repair (in blue) is: , , .

[0300] positive semidefinite direction and Strictly maintain the original value (blue). Projected to (Gray markings indicate that only negative definite directions of non-physical solutions are discarded).

[0301] The key conclusion is: positive semi-stator spatial orientation. and Strictly maintain that no systematic deviation is introduced within the positive semi-stator space.

[0302] Figure 5 This is a comparison chart of three repair strategies; supporting S3 eigenvalue projection positive semi-definite repair. The chart compares the three strategies in... , , Compare the results on the three eigenvalues:

[0303] Original value (orange, (Including non-physical understanding) is:

[0304] , , ;

[0305] The projection method of this invention (blue, (Strictly maintaining the positive semi-definite direction) is as follows:

[0306] , , ;

[0307] Traditional rollback method (gray, (where all directions are arbitrarily changed, introducing systematic bias) is:

[0308] , , .

[0309] The diagram is marked This strictly maintains the advantage of the projection method of this invention being unbiased in the positive semi-stator space, as well as the drawback of the backtracking method introducing arbitrary deviations in all directions.

[0310] The core conclusion is that the key advantage of the eigenvalue projection method lies in discarding only the negative definite components of non-physical solutions, while maintaining all positive semi-definite directions strictly unchanged, without introducing additional systematic bias. The figure uses a Nature colorblind-friendly palette; blue represents the projection method of this invention, orange represents the original result containing negative eigenvalues, and gray represents the traditional backtracking method.

[0311] Although the embodiments of the present invention have been disclosed above, they are not limited to the applications listed in the specification and embodiments. They can be applied to various fields suitable for the present invention. For those skilled in the art, other modifications can be easily made. Therefore, without departing from the general concept defined by the claims and their equivalents, the present invention is not limited to the specific details.

Claims

1. A method for estimating the covariance of hierarchical Bayesian measurement noise, characterized in that, Includes the following steps: S1, Constructing the hierarchical Bayesian prior structure: This is the measurement noise covariance matrix. Construct a two-level hierarchical Bayesian prior, where m is the dimension of the measurement vector; The two-level hierarchical Bayesian prior includes: Global Inverse Wissaud Distribution Layer: Measurement Noise Covariance Matrix Obeying the distribution of the reverse Saudi Arabia ,in The global scale matrix, For global degrees of freedom; Dimensionally Gamma Distribution Layer: For each measurement dimension i, the measurement noise covariance matrix diagonal elements corresponding precision Follows a gamma distribution ,in For shape parameters, For rate parameters, Let molecule be the shape parameter of the gamma distribution in the i-th measurement dimension, satisfying , Let be the rate parameter molecule of the gamma distribution in the i-th measurement dimension; S2, Unbiased innovation decomposition: At each recursive time k, obtain the innovation vector of the current frame. ,in For the current measurement vector, For measurement prediction vectors; According to the definition of new information covariance: ； The superscript T indicates matrix transpose. To measure the Jacobian matrix, For predicting the state covariance matrix; The unbiased empirical estimation matrix is ​​calculated using the following sub-steps: S21, Obtain the prediction covariance Cholesky factor for the current frame. ,satisfy and the measurement-state cross-covariance matrix ; S22, by Solve the trigonometric equations by substituting forward. The intermediate variable matrix is ​​obtained. ,but: ； S23, Calculate the unbiased empirical estimation matrix: ； S3, Eigenvalue projection positive semidefinite repair: Based on the result obtained in S23... The matrix is ​​symmetricized, and then eigenvalue decomposition is performed on the symmetricized matrix. The negative eigenvalues ​​of the decomposed matrix are projected to zero to reconstruct the empirically estimated matrix that satisfies the positive semidefinite requirement. ; S4, Global Inverse Wissaud Distribution Layer Posterior Recursive Update: The hyperparameters of the global inverse Wissaud prior are updated exponentially by using the forgetting factor ρ to drive the exponential moving average recursive update. ， ； in and These represent the global degrees of freedom and the global scale matrix at the current time k, respectively. and These are the global degrees of freedom and the global scale matrix at the previous time step k-1, respectively. like Then, based on the mean formula of the inverse Wissaud distribution, the full matrix estimate of the global inverse Wissaud distribution layer posterior expectation of the current frame is calculated: ； like Then the global inverse Wissaud distribution layer posterior expectation full matrix estimate of the current frame Keep the value from the previous moment unchanged; Output the full matrix estimate of the global inverse Wissaud distribution posterior expectation of the current frame. ; S5, Posterior recursive update of the dimensional gamma distribution layer: For each measurement dimension i, use the innovation vector The i-th component The square value Perform recursive updates of the underlying prior parameters: ， ； in and Let be the shape parameter molecules of the i-th measurement dimension at the current time k and the previous time k-1, respectively; and Let be the numerators of the rate parameters for the i-th measurement dimension at the current time k and the previous time k-1, respectively. like Then, calculate the posterior expectation estimate for that dimension: ； S6, Hierarchical Fusion: S6 receives the global inverse Wissaud distribution layer posterior expectation full matrix estimation from the output of S4. Construct the measurement noise covariance matrix of the final output of the current frame. ,in Represents the measurement noise covariance matrix The i-th diagonal element, Represents the measurement noise covariance matrix The non-diagonal element in the i-th row and j-th column; This represents the full matrix estimation of the posterior expectation of the global inverse Wissaud distribution layer. The i-th diagonal element; This represents the full matrix estimation of the posterior expectation of the global inverse Wissaud distribution layer. The non-diagonal element in the i-th row and j-th column; For each measurement dimension i, determine The i-th diagonal element: like ,but: ； Where the fusion weight coefficient ω∈[0.5,0.8]; like ,but: ； For all off-diagonal elements, i.e., i ≠ j and i,j∈{1,…,m}, the full matrix estimate of the posterior expectation of the global inverse Wissaud distribution layer is provided: ； At this point, the measurement noise covariance matrix... All elements have been determined, and then the entire measurement noise covariance matrix is... Perform symmetry transformation and restrict each diagonal element to a preset range; S7, Output measurement noise covariance matrix This is used for the measurement update step of the filter in the next frame.

2. The method for estimating the hierarchical Bayesian measurement noise covariance as described in claim 1, characterized in that, The entire measurement noise covariance matrix described in S6 Perform symmetry transformation, specifically by transforming the measurement noise covariance matrix. Assigned value .

3. The method for estimating the hierarchical Bayesian measurement noise covariance as described in claim 1, characterized in that, The initial values ​​of each hyperparameter described in S1 are set as follows: Initial values ​​of global degrees of freedom ; Initial value of the global scale matrix ,in The initial measurement noise covariance matrix; Initial shape parameters of gamma distribution in each measurement dimension ; Initial rate parameters of gamma distribution in each measurement dimension ,in The initial measurement noise covariance matrix The diagonal elements.

4. The method for estimating the hierarchical Bayesian measurement noise covariance as described in claim 1, characterized in that, The eigenvalue projection positive semidefinite repair described in S3 specifically includes: right Perform eigenvalue decomposition ,in The eigenvector matrix, For matrix The i-th eigenvalue; Projecting negative eigenvalues ​​to zero: ,in These are the projected eigenvalues; Reconstruct the positive semidefinite matrix: 。 5. The method for estimating the covariance of hierarchical Bayesian measurement noise as described in claim 1, characterized in that, The forgetting factor ρ has a value range of [0.9, 0.99].

6. The method for estimating the covariance of hierarchical Bayesian measurement noise as described in claim 1, characterized in that, Exponential forgetting factor recursion The steady-state value is It approaches 1 / (1-ρ).

7. The method for estimating the hierarchical Bayesian measurement noise covariance as described in claim 1, characterized in that, Recursive calculation of the exponential forgetting factor of the dimensional gamma distribution layer The steady-state value is It approaches 1 / (1-ρ).

8. The method for estimating the covariance of hierarchical Bayesian measurement noise as described in claim 1, characterized in that, In S6, each diagonal element is restricted to a preset range. Inside, among which and These are the lower and upper bounds of the noise variance for the i-th measurement dimension, respectively. The preset interval For the sensor ranging dimension, take the lower bound of the noise variance. upper limit For the sensor's angular measurement dimension, take the lower bound of the noise variance. upper limit .

9. The method for estimating the hierarchical Bayesian measurement noise covariance as described in claim 1, characterized in that, In S6, the measurement noise covariance matrix is ​​calculated. The non-diagonal element in the i-th row and j-th column Then, constrain its value to an interval. Internally, this prevents the cumulative numerical explosion of off-diagonal elements during the recursion of the global scale matrix.

10. The method for estimating the covariance of hierarchical Bayesian measurement noise as described in claim 1, characterized in that, In S2, when numerical anomalies occur in the Cholesky triangulation, the unbiased empirical estimation matrix is ​​lost. When computation fails, the empirically estimated matrix is ​​used from the most recently successfully computed matrix. As the empirical estimation matrix for the current frame; If the current frame is the first frame and no historical frames are available, then the outer product matrix of the current frame's information is used. As the unbiased empirical estimation matrix for the current frame.

11. The method for estimating the hierarchical Bayesian measurement noise covariance as described in claim 1, characterized in that, The fusion weighting coefficient ω described in S6 is adaptively determined in the following manner:

12. Among them The normalized index (NIS) is the exponential moving average of the innovation normalization statistic, and τ is the expected threshold of the NIS. When the filter is functioning correctly... ω approaches 0.5; when the process noise covariance matrix mismatch leads to When deviating from τ, ω is adaptively increased to 0.8, which increases the unbiased anchoring weight of the global inverse Wieshard distribution layer and suppresses the accumulation of bias introduced by process noise pollution in the dimension-by-dimensional gamma distribution layer.

13. The method for estimating the hierarchical Bayesian measurement noise covariance as described in claim 1, characterized in that, In the first N frames, N≥2, steps S1 to S6 are executed normally and the estimated value of the measurement noise covariance matrix is ​​obtained. Subsequently, the operation in S7 used for measurement updates in the next frame is temporarily suspended, and the measurement update step of the filter in the next frame uses the initial measurement noise covariance matrix. Perform filtering updates, and after the number of recursive frames exceeds N frames, execute step S7 to... The measurement update step is passed to the Kalman filter in the next frame.

14. A hierarchical Bayesian method for estimating the noise covariance, characterized in that, include: The prior construction module is used to construct the measurement noise covariance matrix during the Kalman filter initialization phase. Construct a two-layer hierarchical Bayesian prior structure, where the global layer is an inverse Wissaud distribution. The dimension-wise gamma distribution layer is a gamma distribution that is independent of each measurement dimension. ; The unbiased innovation decomposition module is used to obtain the innovation vector at each recursive time k. Obtain the prediction covariance Cholesky factor for the current frame. ,satisfy and the measurement-state cross-covariance matrix Solve the trigonometric equations by substituting forward. Obtain the intermediate variable matrix ,but And explicitly subtract the forecast uncertainty term from the new information extrapolation product. The unbiased empirical estimation matrix is ​​obtained. ; The eigenvalue projection repair module is used to repair the unbiased empirical estimation matrix. Perform symmetry transformation, eigenvalue decomposition, project negative eigenvalues ​​to zero, and reconstruct the positive semidefinite matrix. ; The global inverse Wissaud update module is used to recursively update the posterior parameters of the global inverse Wissaud distribution layer using an exponential moving average recursive method driven by a forgetting factor ρ. The recursive formula is as follows: and ;when At that time, the full matrix estimate of the posterior expectation of the global inverse Wissaud distribution layer is calculated based on the mean formula of the inverse Wissaud distribution. ;when At that time, the global inverse Wissaud distribution layer posterior expectation full matrix estimation Keeping the values ​​from the previous time step unchanged, and outputting the full matrix estimate of the posterior expectation of the global inverse Wissaud distribution layer for the current frame. ; The dimension-wise gamma update module is used to apply the innovation vector for each measurement dimension i. The i-th component The square value The recursive formula for updating the gamma posterior parameters is as follows: and ;like Then calculate the posterior expectation of each dimension. ; The hierarchical fusion module is used to receive the full matrix estimate of the posterior expectation of the global inverse Wissaud distribution layer. Construct the measurement noise covariance matrix of the final output of the current frame. ,in For the measurement noise covariance matrix The i-th diagonal element, For the measurement noise covariance matrix The non-diagonal element in the i-th row and j-th column; Full matrix estimation of the posterior expectation of the global inverse Wissaud distribution layer The i-th diagonal element, Full matrix estimation of the posterior expectation of the global inverse Wissaud distribution layer The non-diagonal element in the i-th row and j-th column; for each measurement dimension i, when At that time, the diagonal elements of the full matrix estimate of the posterior expectation of the global inverse Wissaud distribution layer are... With dimensional posterior expectation By performing weighted fusion, we obtain ;when At that time, the corresponding value of the full matrix estimate of the posterior expectation of the global inverse Wissaud distribution layer is directly taken. For all off-diagonal elements, i.e., i ≠ j and i,j∈{1,…,m}, the full matrix estimate of the posterior expectation of the global inverse Wissaud distribution layer is provided. ; and the measurement noise covariance matrix Perform symmetry transformation and restrict each diagonal element to a preset range.

15. A computer-readable storage medium having a computer program stored thereon, characterized in that, When the computer program is executed by the processor, it implements the method for estimating the hierarchical Bayesian measurement noise covariance as described in any one of claims 1-12.

Images