A method of estimating the state of a time-varying unknown noise covariance suspension system

Abstract

The present application belongs to the technical field of suspension system state estimation, and relates to a kind of estimation method of time-varying unknown noise covariance suspension system state, including using two degrees of freedom linear suspension model to analyze suspension dynamic response, simulate the complex dynamic characteristics of suspension system, estimate the state of suspension system;Using variational bayes method, time-varying noise covariance can be inferred from inverse wishart distribution, and then the state estimation is optimized through the finite sampling posterior probability distribution function of noise covariance and backward kalman smoothing;Unknown noise covariance is identified by using the road grade classification method based on multi-objective optimization and linear classifier, and the road grade is estimated.The present application combines multi-objective optimization and linear classifier, and proposes a road classification method with suspension dynamic response as input, which can accurately obtain the nominal value of PNCM and MNCM of suspension system;Compared with other filters, the present method has higher state estimation accuracy.

Description

Technical Field

[0001] This invention belongs to the field of suspension system state estimation technology, specifically involving a state estimation method based on the fusion of variational Bayesian adaptive Kalman filter and road classification method, and particularly involving a novel method for estimating the state of a suspension system with time-varying unknown noise covariance. Background Technology

[0002] As a component of a vehicle, the suspension system's dynamic response plays a crucial role in the overall vehicle performance. Therefore, to improve ride comfort, many suspension control methods based on state measurements have been proposed. However, measurement results are affected by sensor uncertainties and unknown road conditions, making accurate measurement of system state variables impossible. The Kalman filter, a recursive state-space model based on optimal estimation, is widely used for dynamic system state estimation, such as localization, tracking, and control.

[0003] In suspension state estimation, filter performance primarily depends on the accuracy of the noise covariance. Inaccurate noise covariance can lead to estimation errors or even filter divergence. Actual vehicle driving conditions often vary over time and are independent of road level, resulting in the process noise covariance matrix (PNCM) and measurement noise covariance matrix (MNCM) having similar characteristics.

[0004] To estimate states with time-varying and unknown noise covariance, many adaptive Kalman filter (AKF) methods have been developed, mainly including the following categories: covariance matching, correlation, maximum likelihood, and Bayesian.

[0005] Among them, the covariance matching method uses a sliding window of the innovation sequence to approximate the true noise covariance. Due to its simplicity of implementation, it is widely used in many publications. However, the estimation accuracy of this method is closely related to the window size and requires prior knowledge of the innovation sequence.

[0006] The basic idea of ​​correlation methods is to correlate the autocorrelation function of the system's output or innovation sequence with unknown system parameters. Existing techniques use correlation methods to estimate the mean constant coefficient curve (MNCM) and identify the normalized non-constant coefficient curve (PNCM) based on road classification methods. However, correlation methods are primarily applicable to systems with constant coefficients, which require the assumption that the system output is in a steady state.

[0007] Maximum likelihood estimation (MLP) defines a probability density function, commonly called the likelihood function, with unknown parameters as independent variables, and solves for its extrema. Some improved MLP methods, for example, use the maximum likelihood principle (MLP) to estimate unknown noise in nonlinear dynamic systems and propose a novel expectation-maximization technique to iteratively optimize the MLP computation. However, this method is still not suitable for estimating unknown noise parameters with time-varying characteristics.

[0008] Bayesian methods use Bayes' theorem to iteratively solve for the posterior PDF of the state and unknown parameters. For some complex probability distributions, variational inference can be combined to simplify the computation. In recent research, estimation errors have been further reduced by combining Kalman smoothers with variational Bayesian methods. However, their outputs ignore the smoothed data sequence and approximate the prior of the covariance matrix with a single data value, leading to the possibility of noisy covariance estimation getting trapped in local optima.

[0009] For existing variational Bayesian adaptive Kalman filter methods, the noise covariance estimation process heavily relies on the parameter selection of the approximate distribution. Convergence of the iterative process can be accelerated by choosing a suitable conjugate prior distribution and assuming that the nominal values ​​of PNCM and MNCM are known. However, for unknown road conditions, the nominal values ​​of the noise covariance corresponding to different road grades are clearly not fixed. Summary of the Invention

[0010] The purpose of this invention is to provide a state estimation method based on Variational Bayesian Adaptive Kalman Filter (VBAKF) and road classification for suspension systems with time-varying and unknown noise covariance, thereby solving the problem of accurately estimating the suspension state with time-varying and unknown noise covariance. This invention uses the VB method to infer the time-varying noise covariance from the inverse Wishart distribution, and then optimizes the state estimation through the finite-sampled posterior probability distribution function (PDF) of the noise covariance and backward Kalman smoothing. It also proposes a road classification algorithm based on multi-objective optimization and a linear classifier to identify the unknown noise covariance. The objective of this invention is achieved through the following technical solutions:

[0011] A method for estimating the state of a suspension system with time-varying unknown noise covariance, characterized by comprising the following steps:

[0012] A. Analyze the dynamic response of the suspension using a two-degree-of-freedom linear suspension model, simulate the complex dynamic characteristics of the suspension system, and estimate the state of the suspension system.

[0013] B. Approximate the time-varying noise covariance matrix using the inverse Wishart distribution, and then optimize the state estimation using the finite sampling posterior probability distribution function of the noise covariance and backward Kalman smoothing.

[0014] C. Road level estimation is performed using a road level classification method based on multi-objective optimization and linear classifier. When the suspension system is under unknown road conditions, the dynamic response of the suspension system is used to estimate the current road level, and then the nominal value of the noise covariance matrix corresponding to different road levels is matched.

[0015] Further, step A specifically includes the following steps:

[0016] A1. Construct the dynamic equations of the suspension system;

[0017] A2. Define the system state, measurement output, and disturbance vector respectively, and construct the equations for the discrete-time state-space model;

[0018] A3. Solve for the optimal system state in the least squares sense using a Kalman filter, and estimate the state of the suspension system using a state prediction and state correction iterative model.

[0019] Further, step B specifically includes the following steps:

[0020] B1. Use the variational Bayesian method to calculate the joint probability distribution function from the inverse Wishart distribution to approximate the posterior probability distribution function;

[0021] B2. Use the fixed-point iteration method to find the optimal estimates of the time-varying covariance matrices PECM and MNCM, and use the heuristic parameter propagation method to achieve iterative updates of the approximate posterior probability distribution function.

[0022] B3. Solve the filter state smoothing problem based on the maximum likelihood estimation theory and optimize the state estimation results.

[0023] Further, step C specifically includes the following steps:

[0024] C1. Acquire suspension dynamic response signals: Among them, the suspension dynamic response signals include sprung mass acceleration, suspension dynamic travel and tire dynamic deflection;

[0025] C2. Feature set selection: Feature extraction is performed using time-domain feature extraction, frequency-domain feature extraction, and joint time-frequency-domain feature extraction methods.

[0026] C3. Feature Reduction: The NSGA-II algorithm is used to select the optimal feature subset;

[0027] C4. Road Classification: The recognition accuracy of the road classification model is used as the evaluation function of the NSGA-II algorithm. The receiver operating characteristic curve is plotted, and the area under the curve is calculated by testing the recognition results of different classification models. The unknown road surface level is determined by selecting the L-SVM classification model with the largest area under the curve. When the classification algorithm determines the unknown road level, the nominal value of the empirical noise covariance matrix corresponding to different roads can be referenced.

[0028] Further, in step C1, the sampling frequency of the dynamic response of the suspension system is set to 100Hz, which includes not only the natural frequency of the vehicle body mass (1-2Hz) but also the natural frequency of the unsprung mass (10-16Hz).

[0029] Further, in step C2, the time-domain features include signal statistics such as root mean square, variance, and interpeak; the frequency-domain features include the root mean square of the frequency, the standard deviation of the frequency, and the amplitude power spectrum; the time-frequency joint features are the statistical energy of different frequency bands obtained after six-layer wavelet packet decomposition.

[0030] Further, in step C3, the population size is 150, the number of generations is 200, and the crossover and mutation probabilities are 0.9 and 0.1, respectively.

[0031] Further, in step C4, the classification models used for testing include Naive Bayes, narrow-bandwidth neural networks, k-nearest neighbors, linear support vector machines, and decision trees, all computed by the MATLAB Classifier Toolbox.

[0032] Compared with the prior art, the beneficial effects of the present invention are:

[0033] This invention presents a novel state estimation method for a suspension system with time-varying unknown noise covariance. The method is simple, easy to implement, and readily applicable in practice. Specifically, this method has the following advantages:

[0034] 1. This invention uses the inverse Wishart distribution to approximate the time-varying covariance matrices PECM and MNCM, and then optimizes the state estimation through the finite-sampled posterior PDF of the noise covariance and backward Kalman smoothing. This method has higher state estimation accuracy than other filters.

[0035] 2. This invention combines multi-objective optimization and linear classifiers to propose a road classification method with suspension dynamic response as input, which can accurately obtain the nominal values ​​of PNCM and MNCM of the suspension system;

[0036] 3. Compared with existing filters, the new state estimation method has made better improvements in suspension state estimation. Attached Figure Description

[0037] To more clearly illustrate the technical solutions of the embodiments of the present invention, the accompanying drawings used in the embodiments will be briefly introduced below. It should be understood that the following drawings only show some embodiments of the present invention and should not be regarded as a limitation on the scope. For those skilled in the art, other related drawings can be obtained based on these drawings without creative effort.

[0038] Figure 1 Schematic diagram of a two-degree-of-freedom linear suspension model;

[0039] Figure 2 A schematic diagram of six-layer wavelet packet decomposition;

[0040] Figure 3 Schematic diagram of road classification methods;

[0041] Figure 4 Receiver operating characteristic (ROC) curve;

[0042] Figure 5 Flowchart of a state estimation method based on VBAKF and road classification. Detailed Implementation

[0043] The present invention will be further described below with reference to embodiments:

[0044] The present invention will now be described in further detail with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative of the invention and not intended to limit it. Furthermore, it should be noted that, for ease of description, the accompanying drawings show only the parts relevant to the present invention, and not all of the structures.

[0045] It should be noted that similar reference numerals and letters in the following figures indicate similar items; therefore, once an item is defined in one figure, it does not need to be further defined and explained in subsequent figures. Furthermore, in the description of this invention, terms such as "first," "second," etc., are used only to distinguish descriptions and should not be construed as indicating or implying relative importance.

[0046] This invention proposes a state estimation method based on VBAKF and road classification for estimating suspension states with time-varying and unknown noise covariance. In this method, it is assumed that the conjugate prior distribution of the matrix approximates an IW distribution, and the suspension state estimation process is essentially still performed by VBAKF. However, to fully utilize the effective data during the iteration process and improve the accuracy of state estimation, a heuristic propagation method is used to perform finite sampling of the posterior noise variance, and the propagated estimation samples are folded. Then, to overcome the limitation of traditional VBAKF in handling unknown road slopes well, this invention proposes a road classification method combining multi-objective optimization and a linear classifier of suspension dynamic response.

[0047] The present invention provides a method for estimating the state of a suspension system with time-varying unknown noise covariance, comprising the following steps:

[0048] 1. A two-degree-of-freedom linear suspension model is used to analyze the dynamic response of the suspension, simulate the complex dynamic characteristics of the suspension system, and estimate the state of the suspension system.

[0049] Specifically, the following steps are included:

[0050] 11. Construct the dynamic equations of the suspension system.

[0051] 12. Define the system state, measurement output, and disturbance vector respectively, and construct the equations for the discrete-time state-space model.

[0052] 13. Solve for the optimal system state in the least squares sense using a Kalman filter, and estimate the state of the suspension system using a state prediction and state correction iterative model.

[0053] 2. The time-varying noise covariance matrix is ​​approximated by the inverse Wishart distribution, and then the state estimation is optimized by the finite sampling posterior probability distribution function of the noise covariance and the backward Kalman smoothing.

[0054] Specifically, the following steps are included:

[0055] 21. Use the variational Bayesian method to calculate the joint probability distribution function from the inverse Wishart distribution to approximate the posterior probability distribution function.

[0056] 22. Use the fixed-point iteration method to find the optimal estimates of the time-varying covariance matrices PECM and MNCM, and use the heuristic parameter propagation method to achieve iterative updates of the approximate posterior probability distribution function.

[0057] 23. Solve the filter state smoothing problem based on the maximum likelihood estimation theory and optimize the state estimation results.

[0058] 3. A road level classification method based on multi-objective optimization and linear classifier is adopted to estimate the road level. When the suspension system is under unknown road conditions, the dynamic response of the suspension system is used to estimate the current road level, and then the nominal value of the noise covariance matrix corresponding to different road levels is matched.

[0059] Specifically, the following steps are included:

[0060] 31. Acquire suspension dynamic response signals: The suspension dynamic response signals include sprung mass acceleration, suspension dynamic travel, and tire dynamic deflection. The acquisition frequency of the suspension system dynamic response is set to 100Hz, which includes not only the natural frequency of the vehicle body mass (1-2Hz) but also the natural frequency of the unsprung mass (10-16Hz).

[0061] 32. Feature Set Selection: Feature extraction is performed using time-domain feature extraction, frequency-domain feature extraction, and joint time-frequency-domain feature extraction. The time-domain features include signal statistics such as root mean square, variance, and interpeak frequency; the frequency-domain features include the root mean square of frequency, standard deviation of frequency, and amplitude power spectrum; the joint time-frequency features are the statistical energy of different frequency bands obtained after six-level wavelet packet decomposition.

[0062] 33. Feature Reduction: The NSGA-II algorithm is used to select the optimal feature subset. Population size = 150, number of generations = 200, crossover and mutation probabilities are 0.9 and 0.1, respectively.

[0063] 34. Road Classification: The recognition accuracy of the road classification model is used as the evaluation function of the NSGA-II algorithm. A receiver operating characteristic (ROC) curve is plotted, and the area under the curve (AUC) is calculated by testing the recognition results of different classification models. The L-SVM classification model with the largest AUC is selected to determine the unknown road surface level. When the classification algorithm determines the unknown road level, the nominal value of the empirical noise covariance matrix corresponding to different roads can be referenced. The classification models used for testing include Naive Bayes, narrow-bandwidth neural networks, k-nearest neighbors, linear support vector machines, and decision trees, all computed using the MATLAB Classifier Toolbox.

[0064] Example 1

[0065] A method for estimating the state of a suspension system with time-varying unknown noise covariance includes the following steps:

[0066] Step 1: To simulate the complex dynamic characteristics of the suspension system, this invention uses a two-degree-of-freedom (2-DOF) linear suspension model to analyze the suspension dynamic response. The structure of the suspension model is as follows: Figure 1 As shown.

[0067] According to Newton's second law, the dynamic equations of the suspension system are (1) and (2):

[0068]

[0069]

[0070] When the system state, measurement output, and disturbance vector are defined as follows:

[0071]

[0072] The equations for the discrete-time state-space model can be written as:

[0073] x k =Ax k-1 +Bw k-1 #(3)

[0074] y k =Cx k +v k #(4)

[0075] Where k is the instantaneous discrete sampling time, and x is the value with d x A system state vector of dimension y is a system state vector with d dimensions. y The system measurement vector is dimensional, and A and C are the state transition matrix and observation matrix, respectively. The process disturbance vector w and measurement noise vector v are defined as w ~ N(0,Q) and v ~ N(0,R), respectively.

[0076] The basic idea of ​​the standard Kalman filter is to solve for the optimal system state in the least squares sense and estimate the state of the suspension system through the iterative model of state prediction (5) and state correction (6).

[0077] Prediction process:

[0078]

[0079] Update process:

[0080]

[0081] For suspension systems, imprecise PNCM(Q) and MNCM(R) often exhibit time-varying and unknown characteristics. The time-varying characteristic is that the covariance fluctuates slowly over time, while the unknown characteristic is that the magnitude of the covariance is unknown for different road levels.

[0082] Step 2: Variational Bayes

[0083] According to Bayesian optimal filtering theory, in order to obtain the state estimate x k ,p(x k ,P k|k-1 ,R k |y 1:k The joint posterior probability distribution function (PDF) of the given condition needs to be solved first. Formally, this PDF only has approximate solutions, and a well-known solution is the VB approximation. Using the VB method, this posterior PDF can be extended as follows:

[0084] p(x k ,P k|k-1 ,R k |y 1:k )≈q(x k )q(P k|k-1 )q(R k )#(7)

[0085] In the formula, the optimal approximation q(·) can be solved by minimizing the Kullback–Leibler (KL) divergence between the approximate distribution and the true joint posterior:

[0086]

[0087] Using the variational Bayesian computation method, the solution to (8) can be simplified to the following form:

[0088]

[0089] In the formula, / x k x represents k Other elements and with x k Complementary Indicates P k∣k-1 and R k The expected value. Next, we only need to calculate the joint probability distribution function p(x). k ,P k∣k-1 ,R k ,y 1:k To approximate the posterior probability distribution function q(x) k ),q(P k∣k-1 ) and q(R) k ).

[0090] In the recursive model of KF(5)-(6), the filtered PDFp(x) k |y 1:k-1 ) and p(y k |x k It is approximately a Gaussian distribution.

[0091] p(y k |x k )≈N(y k ;Cx k ,R k )#(10)

[0092]

[0093] Within the theoretical framework of VB, P k∣k-1 ,R k The conjugate prior distribution of can be assumed to be the inverse Wishart distribution because its prior distribution has the same functional form as the multivariate normal distribution, i.e.,

[0094]

[0095]

[0096] Considering that the inverse Wishart distribution and the Gaussian distribution have many constant terms that do not affect the shape of the PDF, it is safe to delete the PDF. Therefore, the simplified density function can be written as:

[0097]

[0098] in, Indicates having d Ω ×d Ω The noise covariance of the matrix dimension, IW(Ω; λ, ψ), represents the IW distribution with DOF parameter λ and inverse scaling matrix ψ. The mean Ω is...

[0099]

[0100] System status x k and PECMPk∣k-1 The measurement noise v is only relevant to data at time k-1. k These are random, independent vectors. Therefore, using the concept of conditional probability, the joint probability distribution function p(x) k ,P k∣k-1 ,R k ,y 1:k It can be broken down into:

[0101]

[0102] Using (14) and (16), logp(x k ,P k∣k-1 ,R k ,y 1:k This can be represented as:

[0103]

[0104] Using (17) and (9), we can obtain the following:

[0105]

[0106]

[0107]

[0108] in,

[0109]

[0110]

[0111] From equations (18) and (19), it can be seen that the posterior probability distribution function (PDF) of time k conforms to the functional form of the IW distribution. Furthermore, according to equations (12) and (13), ... it has the following density:

[0112]

[0113]

[0114] By matching the elements of (18)-(23), the posterior probability distribution function parameters of (22)(23) can be updated as follows:

[0115]

[0116] Since the one-step approximation method cannot accurately estimate (24)PECM(P) k∣k-1 ) and MNCN(R kThe distribution parameters of q(P) are such that a fixed-point iteration method is used to find the optimal estimate. Furthermore, a heuristic parameter propagation method is used to achieve an approximate posterior PDF q(P). k∣k-1 ),q(R k The iterative update of ( ). As follows:

[0117]

[0118] Where ρ∈(0,1] is the forgetting factor value, which reflects the time-varying fluctuation of the noise covariance. When ρ=1, it indicates that the noise covariance is stable.

[0119] According to (15), the i-th iteration estimate of PECM and MNCM can be expressed as the expected value of the IW distribution.

[0120]

[0121]

[0122] Using the delta function, the posterior probability distribution function of the MNCM can be obtained by averaging over finite approximate solutions. Therefore, the posterior PDF p(R∣y) 1:k It can be expressed in the following form

[0123]

[0124] Using Bayes' theorem The calculation can be performed as follows:

[0125]

[0126] in, The random sampling is independent and identically distributed, so

[0127] The random conditional probability distribution function has a mean.

[0128]

[0129] covariance

[0130]

[0131] Therefore, there is

[0132]

[0133] To minimize The mean square error of the state estimate must be equal to the conditional mean. The calculation is as follows:

[0134]

[0135] When updating PDF R through N steps of fixed-point iteration k When a finite approximate estimation set is obtained, the integral operation in (33) is replaced by the sum of this finite set, and the state estimate can be written as

[0136]

[0137] Furthermore, the state estimation error covariance can be derived from... Reconstructing a finite set into different The weighted sum is as follows:

[0138]

[0139] Output time series through time series measurement N-step solution x k The optimal estimation problem is called fixed interval smoothing. Using maximum likelihood estimation theory

[26] , the solution to the smoothing problem is expressed as:

[0140]

[0141] The above process of the adaptive Kalman filtering scheme proposed in this invention is summarized in the variational Bayes algorithm.

[0142] INPUT:A,C, u0,t0,ρ,N

[0143] OUTPUT:

[0144] FORWARD PROCESS:

[0145] for k=1to T do

[0146]

[0147]

[0148] Calculate t k|k-1 ,T k|k-1 ,u k|k-1 U k|k-1 ,employing(25)

[0149] Initialization:

[0150] for i=1to N do

[0151]

[0152]

[0153]

[0154] Update:

[0155]

[0156]

[0157] Update:

[0158]

[0159]

[0160]

[0161] end for

[0162]

[0163]

[0164]

[0165] end for

[0166] Backward Process:

[0167] Initialization:

[0168] for k=T-1to 0do

[0169] Smoothing:

[0170]

[0171]

[0172]

[0173] end for

[0174] Step 3: Road Classification Algorithm

[0175] When the suspension system is under unknown road conditions, the dynamic response of the suspension system is used to estimate the current road level. Then, by matching the nominal values ​​corresponding to different road levels, the nominal PNCM can be selected for Algorithm 1. and MNCN

[0176] To improve the accuracy of road grade estimation, a road grade classification method based on multi-objective optimization and a linear classifier is proposed, as shown below:

[0177] 1. Suspension dynamic response signal acquisition:

[0178] The suspension dynamic response signal includes sprung mass acceleration, suspension dynamic travel, and tire dynamic deflection. Considering that the road spectrum corresponding to speeds of 10-40 km / h is only 0.44-28.3 Hz, the sampling frequency of the suspension system response can be set to 100 Hz. This includes not only the natural frequency of the vehicle body mass (1-2 Hz) but also the natural frequency of the unsprung mass (10-16 Hz).

[0179] 2. Feature set selection:

[0180] Commonly used feature extraction methods include time-domain feature extraction, frequency-domain feature extraction, and joint time-frequency domain feature extraction.

[0181] Time-domain features include common signal statistics such as root mean square (RMS), variance, and peak-to-peak values. Frequency-domain features include the RMS of frequency, standard deviation of frequency, and amplitude power spectrum. Joint time-frequency features are the statistical energy of different frequency bands obtained after six-level wavelet packet decomposition, and their structure is as follows: Figure 2 As shown.

[0182] 3. Feature reduction:

[0183] The feature set generated in step 2 contains a large amount of irrelevant and redundant information, which will reduce the performance of the road classification model. In order to quickly and accurately obtain the optimal feature combination, this invention transforms the feature dimensionality reduction problem into a multi-objective optimization problem and uses the NSGA-II algorithm to select the optimal features, with the following parameters: population size = 150, number of generations = 200, and crossover and mutation probabilities of 0.9 and 0.1, respectively.

[0184] 4. Classification:

[0185] The recognition accuracy of the road classification model needs to be used as an evaluation function for the NSGA-II algorithm. Therefore, in order to fully learn the dynamic response characteristics of the suspension, Figure 4Receiver Operating Characteristic (ROC) curves were plotted, and the Area Under the ROC Curve (AUC) value was calculated by testing the recognition results of different classification models. The classification models used for testing included Naive Bayes (NB), Narrow Bandwidth Neural Network (N-NN), k-Nearest Neighbor (KNN), Linear Support Vector Machine (L-SVM), and Decision Tree (DT), all computed using the MATLAB Classifier Toolbox. By selecting the L-SVM classification model with the highest AUC value, the unknown road surface level could be determined effectively.

[0186] Finally, when the classification algorithm determines the level of an unknown road, it can refer to the empirical nominal values ​​corresponding to different roads.

[0187] The steps of the road classification scheme proposed above are as follows: Figure 3 As shown.

[0188] Note that the above description is merely a preferred embodiment of the present invention and the technical principles employed. Those skilled in the art will understand that the present invention is not limited to the specific embodiments described herein, and various obvious changes, readjustments, and substitutions can be made without departing from the scope of protection of the present invention. Therefore, although the present invention has been described in detail through the above embodiments, the present invention is not limited to the above embodiments, and may include many other equivalent embodiments without departing from the concept of the present invention, the scope of which is determined by the scope of the appended claims.

Claims

1. A method for estimating the state of a suspension system with time-varying unknown noise covariance, characterized in that, Includes the following steps: A. Analyze the dynamic response of the suspension using a two-degree-of-freedom linear suspension model, simulate the complex dynamic characteristics of the suspension system, and estimate the state of the suspension system. B. Approximate the time-varying noise covariance matrix using the inverse Wishart distribution, and then optimize the state estimation using the finite sampling posterior probability distribution function of the noise covariance and backward Kalman smoothing. C. Road level estimation is performed using a road level classification method based on multi-objective optimization and linear classifier. When the suspension system is under unknown road conditions, the dynamic response of the suspension system is used to estimate the current road level, and then the nominal value of the noise covariance matrix corresponding to different road levels is matched.

2. The method for estimating the state of a time-varying unknown noise covariance suspension system according to claim 1, characterized in that: Step A specifically includes the following steps: A1. Construct the dynamic equations of the suspension system; A2. Define the system state, measurement output, and disturbance vector respectively, and construct the equations for the discrete-time state-space model; A3. Solve for the optimal system state in the least squares sense using a Kalman filter, and estimate the state of the suspension system using a state prediction and state correction iterative model.

3. The method for estimating the state of a time-varying unknown noise covariance suspension system according to claim 1, characterized in that, Step B specifically includes the following steps: B1. Use the variational Bayesian method to calculate the joint probability distribution function from the inverse Wishart distribution to approximate the posterior probability distribution function; B2. Use the fixed-point iteration method to find the optimal estimates of the time-varying covariance matrices PECM and MNCM, and use the heuristic parameter propagation method to achieve iterative updates of the approximate posterior probability distribution function. B3. Solve the filter state smoothing problem based on the maximum likelihood estimation theory and optimize the state estimation results.

4. The method for estimating the state of a time-varying unknown noise covariance suspension system according to claim 1, characterized in that, Step C specifically includes the following steps: C1. Acquire suspension dynamic response signals: Among them, the suspension dynamic response signals include sprung mass acceleration, suspension dynamic travel and tire dynamic deflection; C2. Feature set selection: Feature extraction is performed using time-domain feature extraction, frequency-domain feature extraction, and joint time-frequency-domain feature extraction methods. C3. Feature Reduction: The NSGA-II algorithm is used to select the optimal feature subset; C4. Road Classification: The recognition accuracy of the road classification model is used as the evaluation function of the NSGA-II algorithm. The receiver operating characteristic curve is plotted, and the area under the curve is calculated by testing the recognition results of different classification models. The unknown road surface level is determined by selecting the L-SVM classification model with the largest area under the curve. When the classification algorithm determines the unknown road level, the nominal value of the empirical noise covariance matrix corresponding to different roads can be referenced.

5. The method for estimating the state of a time-varying unknown noise covariance suspension system according to claim 4, characterized in that: Step C1: The sampling frequency of the suspension system dynamic response is set to 100 Hz, which includes not only the natural frequency of the vehicle body mass (1-2 Hz) but also the natural frequency of the unsprung mass (10-16 Hz).

6. The method for estimating the state of a time-varying unknown noise covariance suspension system according to claim 4, characterized in that: Step C2: Time-domain features include root mean square, variance, and interpeak signal statistics; frequency-domain features include root mean square of frequency, standard deviation of frequency, and amplitude power spectrum; joint time-frequency features are the statistical energy of different frequency bands obtained after six-level wavelet packet decomposition.

7. The method for estimating the state of a time-varying unknown noise covariance suspension system according to claim 4, characterized in that: Step C3: Population size = 150, number of generations = 200, crossover and mutation probabilities = 0.9 and 0.1, respectively.

8. The method for estimating the state of a time-varying unknown noise covariance suspension system according to claim 1, characterized in that: Step C4 uses classification models for testing, including Naive Bayes, narrow-bandwidth neural networks, k-nearest neighbors, linear support vector machines, and decision trees, all computed by the MATLAB Classifier Toolbox.

Images