Tool condition monitoring method based on Kalman filter fusion model

Abstract

The present invention proposes a tool state monitoring method based on a fusion model of Kalman filtering. The method collects tool processing signals and the amount of flank wear on the cutting edge, preprocesses the collected signals, uses different methods to remove noise from different preprocessed signals, extracts features from the processed signals, and establishes a fusion prediction model for tool wear state. A physical model containing a recursive relationship between moments is used as the prediction equation. An appropriate data model is selected and substituted into the extracted features to predict tool wear, and the wear values ​​are directly applied to the observation equation. While fully utilizing the characteristics of the sensor, the present invention reduces data dimensionality. By establishing a recursive relationship between the moments before and after wear using a physical model, the method enables real-time online monitoring of tool wear. A simple prediction structure is required to rapidly establish the model and implement real-time monitoring.

Description

Technical Field

[0001] The present invention relates to a real-time monitoring technology for tool status, and in particular to a tool status monitoring method based on a fusion model of Kalman filtering. Background Art

[0002] The wear state of tooling plays a significant role in both workpiece machining quality and machining costs. Uncertainty in tool wear is one of the bottlenecks restricting machining system reliability. Using tool monitoring technology can reduce downtime, increase the proportion of tool usage time relative to the tool's total lifespan, improve effective machining time, and reduce machining costs. Determining tool wear under different machining conditions is complex. Traditional wear determination relies on manual experience, which is highly subjective and generally adopts a conservative approach, requiring tool replacement before reaching maximum tool life. Existing publicly available or known technologies primarily include direct and indirect monitoring methods. Direct monitoring primarily uses high-speed cameras, but is easily affected by chips and cutting fluid. Indirect monitoring involves collecting sensor signals during the machining process and establishing a correlation between these signals and tool wear. However, models constructed solely from data lack generality and interpretability. Therefore, developing a real-time, online monitoring system for tooling status is of great significance. Summary of the Invention

[0003] In view of the above problems, the present invention proposes a tool state monitoring method based on a fusion model of Kalman filtering. By using a fusion method, integrating the advantages of multiple models and improving defects, different denoising methods are selected by analyzing different signal characteristics, and features that are highly correlated with the tool state are selected. The physical and data models are integrated to achieve rapid and time-related prediction of tool wear.

[0004] The present invention adopts the following technical solutions to achieve the above-mentioned purpose. A tool state monitoring method is implemented based on a fusion model of Kalman filtering, and the steps are as follows:

[0005] 1) Collect the signal of tool processing and the wear of cutting edge flank:

[0006] S01: Install a dynamometer on the fixture to measure the force signal F in the XYZ directions x ,F y ,F z ; Install vibration sensors in the tool feed direction, spindle radial direction and axial direction to measure the vibration signal V in the three directions of XYZ x ,V y ,V z ; Install an acoustic emission sensor S on the workpiece to measure the root mean square value of the acoustic emission. Then, seven sets of sensor signals are collected in each experiment, namely: F x ,、F y ,、Fz ,、V x ,、V y ,、V z , and S;

[0007] S02: After each cutting, the flank wear is measured three times and the average value is calculated;

[0008] 2) Due to the existence of external environmental noise and internal factors of the machine tool, the collected signals need to be preprocessed:

[0009] S01: Elimination of signal scenario trend items:

[0010] Since the temperature rise of the machine tool itself and the temperature change of the external environment will cause the scene trend item, t is the sampling time, x is the acquisition signal, (t i ,x i ) represents the time t i Collect signal value x i , assuming that the number of sampling points in t is n, x can be fitted with a k-order polynomial:

[0011]

[0012] Where, α j is the degree of the polynomial with the highest degree k;

[0013] If you remember The coefficients are obtained using the least squares method:

[0014] A=(t T t) -1 t T x (2)

[0015] Where t is the coefficient matrix of sampling time, A is the sampling matrix of polynomial degree, and x is the value of the collected signal;

[0016] S02: Elimination of signal zero drift:

[0017] The initial stage of the signal contains an unprocessed state. Assuming that the first k sampling points are unprocessed state signals, the average value of the first k signals is moved to shift the signal to zero:

[0018]

[0019] Where s i is the preprocessed signal;

[0020] 3) Use wavelet denoising method to remove noise from the preprocessed force signal:

[0021] S01: Performing wavelet decomposition and low-frequency reconstruction on the signal can obtain a signal with less high-frequency noise. The wavelet basis function ψ and the decomposition scale m will affect the denoising effect. Based on the signal-to-noise ratio, root mean square error and smoothness, the entropy method is used to construct a unified evaluation standard for the fusion denoising index.

[0022]

[0023] Where SNR is the signal-to-noise ratio, RMSE is the root mean square error, S is the smoothness, and s real =[s real,1 s real,2 …s real,n ] is an ideal signal, is the denoised signal;

[0024] Each indicator value is normalized to the range of 0.1 to 1:

[0025]

[0026] Where, I represents the normalized denoising index;

[0027] Use the entropy method to determine the weight of each indicator:

[0028]

[0029] Where P is the probability, E is the entropy value, W is the weight, and C is a constant set to the total number of participating ranking indicators;

[0030] Combining equations (5) and (6) to construct the fusion denoising index T under different wavelet parameters:

[0031] T(ψ,m)=W SNR ×I SNR (ψ,m)+W RMSE ×I RMSE (ψ,m)+W S ×I S (ψ,m) (7)

[0032] S02: Since the ideal signal cannot be obtained in the actual process, it is necessary to construct a simulated signal, add noise, and then denoise it to calculate the optimal noise reduction parameters;

[0033] A periodic signal can be viewed as a superposition of different sinusoidal signals. In the frequency domain, the signal frequency corresponding to when the energy exceeds a specified threshold is selected, and the frequency, amplitude, and phase of the sinusoidal signal are calculated, and the simulation signal is constructed by superposition.

[0034] Gaussian white noise is added to the simulation signal to simulate the noise under the actual processing state. Since the noise generated each time is random, a noisy signal of multiple simulation signals is constructed and the denoising effect is measured by the mean of the fusion index.

[0035] 4) Use the smoothing method to remove noise from the preprocessed vibration signal:

[0036] Polynomial approximation is used in the smoothing interval, and the number of signal points in the smoothing interval is used to solve the unknown parameters. The smoothing window is set to 5, and the number of polynomial approximations is set to 3. This is also called the 5-point 3-times mean smoothing method. The calculation formula is as follows:

[0037]

[0038] 5) Extract features from the processed signal;

[0039] S01: Unify the experimental sampling time, assuming that the signal sampling time is t0 and the retention time is t * , the starting sampling point of the retention time is t start , the end sampling point is t end , the sampling frequency is ω, then the starting and ending sampling points after deletion are set to:

[0040]

[0041] S02: Time domain feature analysis method: Time domain analysis is performed on multiple experiments of 7 sets of tool sensor signals; including mean, root mean square value, variance, skewness, kurtosis, peak index, waveform index, pulse index and margin index;

[0042] S03: Frequency Domain Characteristic Analysis Method: The relationship between frequency ω and energy p(ω) is calculated through Fourier transform. Frequency domain analysis is performed on multiple experiments of force and acoustic emission in the seven sets of tool sensor signals. This includes center of gravity frequency, mean square frequency, root mean square frequency, frequency standard deviation, frequency variance, and total frequency energy.

[0043] S04: Time-frequency feature wavelet packet analysis: The vibration signal has concentrated energy at both low and high frequencies. A 4-layer wavelet decomposition was used to decompose the vibration signal into 16 frequency bands. The range and variance of the energy of each frequency band were calculated for all experimental groups. The wavelet packet nodes with the largest range and variance values ​​were selected, and the frequency domain feature analysis of S03 was then performed.

[0044] S05: Feature Dimensionality Reduction: Calculate the autocorrelation between different features of the same sensor, retain the most representative features to achieve feature dimensionality reduction, and use the Pearson correlation coefficient to calculate the correlation analysis as follows:

[0045]

[0046] Where λ i ,λ j are the calculation results of the i-th and j-th feature analysis methods respectively; |rij The closer | is to 1, the i ,λ j The higher the correlation, the greater the significance. After the correlation is calculated, the t-test is used to determine the significance.

[0047] After completing the autocorrelation analysis, the most representative features are retained and the features are subjected to dimensionality reduction. The correlation between the dimensionality reduction data and the tool wear status is then established, and highly correlated features are selected for a second dimensionality reduction. For some sensors, all extracted features may have a low correlation with tool wear. In this case, at least one feature from each sensor is selected as the input of the prediction model.

[0048] 6) Establish a fusion prediction model for tool wear status: Establish a physical model containing a time-recursive relationship as the prediction equation, select a suitable data model to substitute the extracted features to predict tool wear, and directly apply the wear value to the observation equation, completing the fusion of multiple models to reduce errors;

[0049] S01: Classification of tool wear stages: Tool wear is divided into three different wear stages, and a different prediction model is established for each stage. Tool wear conditions are converted into first-order derivatives that vary with time to determine how to classify tool wear states. A decision tree is selected as the classification prediction model, with the final features extracted in step 5) as input and the tool wear classification stages as output, enabling rapid and accurate classification of tool wear stages.

[0050] S02: Establish a physical model for tool wear: The initial, mid-term, and late-term tool wear satisfy different prediction equations, and wear failure follows a Weibull distribution. The demarcation point of the three wear states is determined by the derivative image of the tool wear time domain image. It is assumed that the initial wear change rate decreases with time; the mid-term wear change rate remains constant and non-zero with time; and the late-term wear change rate increases with time. The Weibull distribution of initial and late-term wear is assumed to be uniformly accelerated linear motion according to the change law:

[0051]

[0052] The mid-term wear Weibull distribution is assumed to be uniform linear motion according to the variation law:

[0053]

[0054] It is only necessary to measure the wear amount, wear velocity and initial value of wear acceleration in each wear period to establish the recursive relationship expression of the physical model of the entire cycle;

[0055] S03: Establish a data model for tool wear: Establish a connection between the extracted optimal features and tool wear; select multiple data prediction models with relatively simple structures to quickly predict the general trend; for each data model, calculate the error between the predicted wear and the actual wear;

[0056] S04: Establish a fusion prediction model based on Kalman filtering: The change of tool wear over time can be regarded as a continuous random process. The result of S02 in step 6) is used as the prediction equation, and the result of S03 in step 6) is used as the observation equation. The numerical values ​​of the two are combined to accurately predict the tool wear;

[0057] Define the predicted quantity of the tool wear physical model at time t as Y0, and the predicted quantity is y0; select a sensor to extract the feature data model and define the observed quantity as Y1, and the observed quantity is y1; define the probability density function as f; then the posterior probability of tool wear is expressed as:

[0058]

[0059] Where, It is called the likelihood function;

[0060] Select the normal distribution model to simplify the calculation of the likelihood function:

[0061] The time t is extended to the continuous time process, and the prediction and observation equations are constructed. The prediction equation needs to find the recursive relationship between the physical model before and after the time, while the output value of the data model can be directly applied to the observation equation. Assuming that the recursive relationship between the time before and after the prediction equation is F, the result of a data model is used as the judgment value G(Y) of the observation equation at time t. 0,t ), the model noise of prediction and observation is expressed as Q t and R t , the prediction model and observation model in the continuous process are expressed as:

[0062]

[0063] Set the initial value to f0(y0) and make a prediction f1 - (y0) and update f1 + Probability density function calculation of (y0):

[0064]

[0065] Where, Subsequent f2 - ,f2 + ,f3 - ,f3 + ,… are calculated recursively through formula (15);

[0066] The specific updated value after each time t is directly solved through the probability density function:

[0067]

[0068] On this basis, two assumptions are made to solve the infinite integral; the prediction and observation equations are assumed to be linear, and the prediction and observation noises follow a normal distribution:

[0069]

[0070] Where, F c ,G c is a constant, and the assumption is substituted into the recursive process to obtain the KF calculation formula:

[0071]

[0072] Where, is the Kalman gain, are the predicted value and updated value of the physical model at time t, are the errors of the predicted value and the updated value at time t, y 1,t Output value of the data model at time t;

[0073] The prediction equation contains information about tool wear, wear velocity and wear acceleration, which can be rewritten in vector form: The error is rewritten as the covariance matrix σ t →∑ t ; The selected data model is extended from a single to multiple for tool wear prediction and rewritten as Y t =[y 1,t y 2,t y 3,t …]; F c ,G c ,Q t ,R t is a matrix constant, and the KF in matrix form is as follows:

[0074]

[0075] Where, is the Kalman gain, M I is the identity matrix;

[0076] According to the prediction error percentage at different wear stages, the confidence level of prediction and observation is set; the final tool wear curve is obtained.

[0077] The beneficial effects or advantages achieved by the present invention are: selecting different denoising methods for different sensors; extracting features from a time-frequency perspective, and reducing the dimensionality of features through autocorrelation to ensure that each type of sensor contains at least one feature, thereby fully utilizing the characteristics of the sensor while reducing the data dimension; the extracted features are connected to tool wear through a variety of relatively simple data models, which requires less time, and a recursive relationship between the wear amount before and after the physical model is established to achieve real-time online monitoring of tool wear. The errors at each stage of the different sub-models are used to set the error of the fusion model. The advantage of the fusion model is that the errors in the wear amount predicted by different methods are continuously reduced, and only a simple prediction structure is required to achieve rapid model establishment and real-time monitoring. BRIEF DESCRIPTION OF THE DRAWINGS

[0078] Figure 1 is a graph showing the wear measurement value of tool C1 in the present invention;

[0079] Figure 2 is a graph of the early signal processing steps in the present invention;

[0080] Figure 3a is a frequency domain characteristic diagram of the force signal sensor in the present invention;

[0081] Figure 3b It is a frequency domain characteristic diagram of the vibration signal sensor in the present invention;

[0082] Figure 3c It is a frequency domain characteristic diagram of the acoustic emission signal sensor in the present invention;

[0083] Figure 4a It is a histogram of the range of 16 groups of wavelet packet nodes in the present invention;

[0084] Figure 4b It is a histogram of the variance of 16 groups of wavelet packet nodes in the present invention;

[0085] Figure 5a It is the time domain signal correlation heat map of the X-direction force signal in the present invention;

[0086] Figure 5b It is the time domain signal correlation heat map of the X-direction vibration signal in the present invention;

[0087] Figure 5c It is a heat map of the time domain signal correlation of the acoustic emission signal in the present invention;

[0088] Figure 6a It is the frequency domain signal correlation heat map of the X-direction force signal in the present invention;

[0089] Figure 6b It is a heat map of the frequency domain signal correlation of the acoustic emission signal in the present invention;

[0090] Figure 6c It is the frequency domain signal correlation heat map of the (4,5) node of the X-direction vibration signal in the present invention;

[0091] Figure 6d It is the heat map of the frequency domain signal correlation of the (4,13) node of the X-direction vibration signal in the present invention;

[0092] Figure 7a is a suitable time domain feature histogram in the dimension reduction matrix of the present invention;

[0093] Figure 7b is a suitable frequency domain feature histogram in the dimension reduction matrix of the present invention;

[0094] Figure 8a It is a curve diagram for dividing tool wear stages in the present invention;

[0095] Figure 8b It is a dividing curve diagram of the first-order derivative distribution of the tool wear stage in the present invention;

[0096] Figure 9 It is the decision tree classification diagram of the present invention;

[0097] Figure 10a It is a curve diagram of the predicted value of tool C1 before the fusion of the physical model in the present invention;

[0098] Figure 10b This is a curve diagram of the predicted value of tool C1 before the first data model fusion in the present invention;

[0099] Figure 10c This is a curve diagram of the predicted value of tool C1 before the second data model is fused in the present invention;

[0100] Figure 10d This is a graph showing the predicted value of tool C1 before the third data model is fused in the present invention;

[0101] Figure 11a is the error histogram before the physical model fusion in the present invention;

[0102] Figure 11b is the error histogram before the first data model fusion in the present invention;

[0103] Figure 11c is the error histogram before the fusion of the second data model in the present invention;

[0104] Figure 11d is the error histogram before the third data model fusion in the present invention;

[0105] Figure 12 This is a curve diagram showing the predicted wear of tool C1 in the present invention;

[0106] Figure 13a is a graph showing the wear measurement values ​​of tool C4 in the present invention;

[0107] Figure 13b is a graph showing the wear measurement values ​​of tool C6 in the present invention;

[0108] Figure 14a It is a curve diagram for dividing the wear stages of C4 tools in the present invention;

[0109] Figure 14b It is a dividing curve diagram of the first-order derivative distribution of the C4 tool wear stage in the present invention;

[0110] Figure 15a It is a curve diagram for dividing the wear stages of C6 tool in the present invention;

[0111] Figure 15b It is a dividing curve diagram of the first-order derivative distribution of the C6 tool wear stage in the present invention;

[0112] Figure 16a It is a curve diagram of the predicted value of tool C4 before the fusion of the physical model in the present invention;

[0113] Figure 16b This is a curve diagram of the predicted value of tool C4 before the first data model fusion in the present invention;

[0114] Figure 16c This is a curve diagram of the predicted value of tool C4 before the second data model is fused in the present invention;

[0115] Figure 16d This is a graph showing the predicted value of tool C4 before the third data model is fused in the present invention;

[0116] Figure 17a It is a curve diagram of the predicted value of tool C6 before the fusion of the physical model in the present invention;

[0117] Figure 17b This is a curve diagram of the predicted value of tool C6 before the first data model fusion in the present invention;

[0118] Figure 17c This is a curve diagram of the predicted value of tool C6 before the second data model is fused in the present invention;

[0119] Figure 17d This is a graph showing the predicted value of tool C6 before the third data model is integrated in the present invention;

[0120] Figure 18a This is a curve diagram showing the predicted results of the tool C4 fusion model wear in the present invention;

[0121] Figure 18bThis is a curve chart showing the predicted wear results of the tool C6 fusion model in the present invention. DETAILED DESCRIPTION

[0122] The present invention will be further described below with reference to the accompanying drawings and embodiments. Figures 1 to 18b ,A tool condition monitoring method based on the fusion model of Kalman filter is implemented with the following steps:

[0123] Step 1: Collect tool signals and cutting edge flank wear during machining. Specific implementation: Using the open data from the PHM2010 High-Speed ​​CNC Machine Tool Tool Health Prediction Competition, published by the New York Society for Predictive and Health Management in 2010, 315 sets of cutting edge flank wear were measured using C1 tools. The signal sampling frequency was set to 50 kHz. The machining parameters for the milling experiment are shown in Table 1:

[0124] Table 1 Processing parameters of PHM2010 milling experiment

[0125] Spindle speed (rpm) Feed speed (mm / min) Cutting depth (mm) Cutting width (mm) 10400 1555 0.2 0.125

[0126] S01: A dynamometer is installed on the fixture to measure force signals in three directions. Vibration sensors are installed in the tool feed direction, spindle radial direction, and spindle axial direction to measure vibration signals in three directions. An acoustic emission sensor is installed on the workpiece to measure the RMS value of the acoustic emission. A total of seven sets of sensor signals are measured.

[0127] S02: After each cutting, measure the flank wear three times and calculate the average value (such as Figure 1 As shown in the figure, the dotted line, dotted line and dash-dot line represent the measurement results of the tool flank wear of the 1st, 2nd and 3rd experiments respectively, the solid line represents the average of the 3 measurements, the horizontal axis is the number of experiments changing with time, and the vertical axis is the measured tool wear). The above steps are performed for each cutting experiment, and a total of 315 cutting experiments are carried out.

[0128] Step 2: Due to the existence of external environmental noise and internal factors of the machine tool, the collected signal needs to be preprocessed. The preprocessing content includes the elimination of signal scene trend items and zero drift (such as Figure 2 As shown, assuming that the collected sensor signal has a trend term (represented by the circle in the figure), first eliminate the trend term of the signal (represented by the square in the figure), and then eliminate the zero drift phenomenon so that the initial point is located at the 0 scale line (represented by the diamond in the figure):

[0129] S01: Elimination of signal scenario trend items:

[0130] Since the temperature rise of the machine tool itself and the temperature change of the external environment will cause the scene trend item, t is the sampling time, x is the acquisition signal, (t i ,xi ) represents the time t i Collect signal value x i , assuming that the number of sampling points in t is n, x can be fitted with a k-order polynomial:

[0131]

[0132] Where, α j is the degree of the polynomial with the highest degree k;

[0133] If you remember The coefficients are obtained using the least squares method:

[0134] A=(t T t) -1 t T x (2)

[0135] S02: Elimination of signal zero drift:

[0136] The initial stage of the signal contains an unprocessed state. Assuming that the first k sampling points are unprocessed state signals, the average value of the first k signals is moved to shift the signal to zero:

[0137]

[0138] Where s i is the preprocessed signal.

[0139] Step 3: Convert the time domain signal into the frequency domain signal through Fourier transform. The frequency domain signals of different sensors are different (such as Figure 3a Frequency domain image showing part of the force signal, where the main energy is concentrated in the low frequencies; Figure 3b Frequency domain image showing part of the vibration signal, where energy is present in both low and high frequencies; Figure 3c The frequency domain image represents part of the acoustic emission signal, with energy in the low frequency and almost no noise. The appropriate denoising method is selected according to the different frequency domain signals. The useful energy of the force signal is mainly concentrated in the low frequency. Therefore, the wavelet denoising method is used to remove the noise of the preprocessed force signal:

[0140] S01: Performing wavelet decomposition and low-frequency reconstruction on the signal can obtain a signal with less high-frequency noise. The wavelet basis function ψ and the decomposition scale m will affect the denoising effect. Based on the signal-to-noise ratio, root mean square error and smoothness, the entropy method is used to construct a unified evaluation standard for the fusion denoising index.

[0141]

[0142] Where SNR is the signal-to-noise ratio, RMSE is the root mean square error, S is the smoothness, and s real =[s real,1s real,2 …s real,n ] is an ideal signal, is the denoised signal.

[0143] Each indicator value is normalized to the range of 0.1 to 1:

[0144]

[0145] Where I represents the normalized denoising index.

[0146] Use the entropy method to determine the weight of each indicator:

[0147]

[0148] Where P is the probability, E is the entropy value, W is the weight, and C is a constant set to the total number of participating ranking indicators.

[0149] Combining equations (5) and (6) to construct the fusion denoising index T under different wavelet parameters:

[0150] T(ψ,m)=W SNR ×I SNR (ψ,m)+W RMSE ×I RMSE (ψ,m)+W S ×I S (ψ,m) (7)

[0151] S02: Since ideal signals cannot be obtained in the actual process, it is necessary to construct a simulated signal, add noise, and then denoise it to calculate the optimal noise reduction parameters. The periodic signal can be viewed as a superposition of different sinusoidal signals. In the frequency domain, the signal frequency corresponding to the energy exceeding the specified threshold is selected, and the frequency, amplitude, and phase of the sinusoidal signal are calculated. The simulated signal is then superimposed to construct the simulated signal. Taking the tool C1X direction as an example, the different sinusoidal signal parameters included in the calculated simulation signal are shown in Table 2:

[0152] Table 2 Different sinusoidal signal parameters contained in the simulation signal

[0153] Serial number Amplitude frequency Phase 1 8.397036 0 0 2 0.727680627 100 1.660768302 3 2.440561399 150 1.620169753 4 2.822006941 200 -1.85121348 5 0.632422867 250 -0.693428881 6 1.251043333 350 -1.811953406 7 0.52901464 450 -0.632435039 8 1.591572735 500 -0.634078298 9 0.720878262 550 2.537788547 10 1.564030271 700 -3.016074022 11 0.866284518 1000 2.737025144 12 0.791319018 1050 -0.845327962 13 1.398901765 1550 0.578470043 14 1.090421148 2050 -0.501961713 15 1.116381724 2100 2.643372299 16 0.62776108 3650 2.358726536

[0154] Gaussian white noise was added to the simulation signal to simulate the noise in the actual machining state. Since the noise generated each time was random, a noisy signal of multiple simulation signals was constructed and the denoising effect was measured using the mean of the fusion index. The fusion index corresponding to different denoising parameters is shown in Table 3. Ultimately, the wavelet basis function db8 and the number of decomposition layers 2 were determined as the wavelet denoising parameters for the force signal.

[0155] Table 3 Fusion indicators corresponding to different denoising parameters

[0156] Wavelet basis Number of layers SNR RMSE R T 8 2 26.054 0.281 1.086 0.915 6 2 26.015 0.282 1.097 0.913 7 2 25.964 0.284 1.093 0.913 9 2 25.956 0.284 1.091 0.913 5 2 25.864 0.287 1.106 0.909 4 2 25.752 0.291 1.127 0.905 3 2 25.275 0.308 1.194 0.889 2 2 24.173 0.349 1.449 0.843 8 3 21.520 0.473 0.837 0.827 9 3 21.482 0.476 0.827 0.827 … … … … … …

[0157] Step 4: The useful energy of the vibration signal is mainly concentrated in the low and high frequencies, so the pre-processed vibration signal is smoothed to remove noise:

[0158] Polynomial approximation is used in the smoothing interval, and the number of signal points in the smoothing interval is used to solve for the unknown parameters. Taking the smoothing window as 5 and the polynomial approximation order as 3, it is also called the 5-point 3-times mean smoothing method. Its calculation formula is as follows:

[0159]

[0160] Step 5: Extract features from the processed signal.

[0161] S01: Unify the experimental sampling time, assuming that the signal sampling time is t0 and the retention time is t * , the starting sampling point of the retention time is t start , the end sampling point is t end , the sampling frequency is ω, then the starting and ending sampling points after deletion are set to:

[0162]

[0163] Since the acquisition time of most signals is about 4s, the final signal of each group of experiments mainly retains t * =4s.

[0164] S02: Time domain feature analysis method: Time domain analysis is performed on multiple experiments of 7 sets of sensor signals of the tool.

[0165] Table 4 Main methods of time domain analysis

[0166]

[0167] S03: Frequency domain feature analysis method: The relationship between frequency ω and energy p(ω) is calculated through Fourier transform, and frequency domain analysis is performed on multiple experiments of force and acoustic emission in 7 sets of sensor signals of the tool.

[0168] Table 5 Main methods of frequency domain analysis

[0169]

[0170]

[0171] S04: Wavelet packet analysis of time-frequency characteristics. The vibration signal has energy concentration at both low and high frequencies. The vibration signal of tool C1 is decomposed into 16 frequency bands using a 4-layer wavelet decomposition. The range and variance of the energy of each frequency band in 315 groups of experiments are calculated (e.g. Figure 4aIt represents the range of 16 wavelet packet node coefficients in 315 groups of experiments; Figure 4b Indicates the variance of 16 wavelet packet node coefficients in 315 groups of experiments), select the wavelet packet nodes with larger range and variance values. The greater the degree of change in range or variance, the more it can reflect the degree of change in tool wear state at the node. After comprehensive consideration, the wavelet packet nodes are determined to be (4,5) and (4,13), and then the frequency domain feature analysis of S03 is used.

[0172] S05: Feature Dimensionality Reduction. Calculate the autocorrelation between different features of the same sensor and retain the most representative features to achieve feature dimensionality reduction. Correlation analysis uses the Pearson correlation coefficient as follows:

[0173]

[0174] Where λ i ,λ j are the calculation results of the i-th and j-th feature analysis methods respectively. ij The closer | is to 1, the i ,λ j After the correlation is calculated, the t-test is used to determine the significance.

[0175] Take the X direction as an example to show the time domain of force, vibration and acoustic emission signals (such as Figure 5a is the time domain correlation heat map of the force signal; Figure 5b It is the time domain correlation heat map of the vibration signal; Figure 5c is the time domain correlation heat map of acoustic emission signals) and frequency domain (such as Figure 6a is the frequency domain correlation heat map of the force signal; Figure 6b It is the frequency domain correlation heat map of the vibration signal (4,5) node; Figure 6c It is the frequency domain correlation heat map of the vibration signal (4,13) node; Figure 6dIt is the frequency domain correlation heat map of the acoustic emission signal. After completing the autocorrelation analysis, the most representative features are retained and the features are subjected to dimensionality reduction. The dimensionality reduction method is shown in Table 6. Table 6 Time domain features of the X-direction force signal. The 9 time domain features can be clustered into 2 categories. If only the representative features in each category are extracted, the dimensionality reduction of the data from 9 dimensions to 2 dimensions can be achieved. First, the dimensionality reduction of the first category features is performed, and the first 3 groups of features of the entire matrix are taken to form a 3×3 matrix. The sum of the values ​​in the first row or the first column is the same, and the calculated value is 2.89; similarly, the sum of the values ​​in the second row is 2.95; and the sum of the values ​​in the third row is 2.88. Therefore, feature 2 is selected as the representative feature of the time domain features 1, 2, and 3 of the X-direction force signal. The same method is used to reduce the dimension of the second type of features. The summation results of rows 1 to 6 are 5.06, 4.94, 5.34, 5.33, 5.49, and 5.51, respectively. Feature 9 is selected as the representative feature of time domain features 4 to 9 of the X-direction force signal.

[0176] Table 6 Time domain characteristics of X-direction force signal

[0177] feature 1 2 3 4 5 6 7 8 9 1 1 0.98 0.91 0.26 0.59 0.59 0.38 0.56 0.54 2 0.98 1 0.97 0.23 0.54 0.55 0.36 0.52 0.50 3 0.91 0.97 1 0.19 0.42 0.50 0.32 0.47 0.45 4 0.26 0.23 0.19 1 0.81 0.75 0.89 0.80 0.81 5 0.59 0.54 0.42 0.81 1 0.77 0.78 0.80 0.79 6 0.59 0.55 0.50 0.75 0.77 1 0.84 0.99 0.98 7 0.38 0.36 0.32 0.89 0.78 0.84 1 0.90 0.92 8 0.56 0.52 0.47 0.80 0.80 0.99 0.90 1 1.00 9 0.54 0.50 0.45 0.81 0.79 0.98 0.92 1.00 1

[0178] The same method is applied to the extraction of all features. Finally, the features after dimensionality reduction of all original features are shown in Table 7.

[0179] Table 7 Signal characteristics after dimensionality reduction

[0180]

[0181] Then establish the correlation between the dimensionality reduction data and the tool wear status, select highly correlated features for the second dimensionality reduction, and select features with a correlation coefficient greater than 90% (such as Figure 7a The horizontal axis represents the 23 features of the time domain. Figure 7b The horizontal axis represents the 23 features in the frequency domain, and the vertical axis in both figures represents the magnitude of the correlation. The dashed line indicates a correlation coefficient of 90%. The feature number corresponding to the portion above the dashed line is selected. For some sensors, all extracted features may have a low correlation with tool wear. In this case, at least one feature from each sensor is selected as the input to the prediction model.

[0182] Table 8 Signal characteristics after dimensionality reduction

[0183]

[0184] Step 6: Establish a fusion prediction model for tool wear. A physical model containing a recursive time relationship is established as the prediction equation. An appropriate data model is selected to extract features to predict tool wear. The wear values ​​are directly applied to the observation equation, completing the fusion of multiple models to reduce errors.

[0185] S01: Classification of tool wear stages. Divide tool wear into three different wear stages, and establish different prediction models for each stage. Convert tool wear conditions into first-order derivatives that change over time to determine how to divide tool wear states. You can intuitively judge and divide tool wear states (such as Figure 8a is the wear of tool C1; Figure 8b is the first-order derivative of tool C1 wear versus time, and the dotted line indicates that the tool wear is divided into three stages by the first-order derivative. In the 315 experiments, 1 to 90 are divided into the initial wear stage, 91 to 250 are divided into the middle wear stage, and 251 to 315 are divided into the late wear stage.

[0186] The decision tree is selected as the classification prediction model. The input is the final feature extracted in step 5, and the output is the stage of tool wear classification. The classification of tool wear stages can be achieved quickly and accurately. The decision tree structure and parameters (such as Figure 9 As shown in the figure, the root node selects the 15th feature to complete the judgment of the early wear stage, and the child node selects the 4th feature to complete the judgment of the middle wear stage and the late wear stage). The accuracy of wear stage judgment using this decision tree reaches 100%.

[0187] S02: Establish a physical model for tool wear. Initial, intermediate, and late wear of a tool satisfy different prediction equations, and wear failure follows a Weibull distribution. The demarcation point between the three wear states is determined by the derivative image of the tool wear time domain image. It is assumed that the rate of change of initial wear decreases over time; the rate of change of intermediate wear remains constant and non-zero over time; and the rate of change of late wear increases over time. The Weibull distributions for initial and late wear are assumed to be uniformly accelerated linear motion based on the variation pattern:

[0188]

[0189] The mid-term wear Weibull distribution is assumed to be uniform linear motion according to the variation law:

[0190]

[0191] Simply measuring the wear volume, wear velocity, and initial wear acceleration at each wear stage allows for the development of a recursive physical model for the entire cycle. For tool C1, the physical model predictions for the three classification stages are performed with initial wear values ​​of w = 40, v = 1.4, and a = -0.02; mid-term wear values ​​of w = 90, v = 0.3; and late-term wear values ​​of w = 135, v = 0.4, and a = 0.04.

[0192] S03: Build a data model for tool wear. Use the extracted optimal features to establish a correlation with tool wear. Select multiple, relatively simple data prediction models to quickly predict general trends. For each data model, calculate the error between the predicted wear and the actual wear.

[0193] Construct three data models for tool wear. Select 20 decision trees to build the first random forest data model; use BP neural network with a structure of 18-30-10-1 to build the second data model; select support vector machine with Gaussian kernel function to build the third data model. Obtain tool wear amount of different prediction models (such as Figure 10a Wear prediction for physical models; Figure 10b The wear amount prediction for the first data model is random forest; Figure 10c Wear prediction for the second data model, namely the neural network; Figure 10d The wear prediction of the third data model, namely the support vector machine) and the errors of different prediction models (such as Figure 11a is the prediction error histogram of the physical model; Figure 11b is the prediction error histogram of the first data model, namely random forest; Figure 11c is the prediction error histogram of the second data model, namely the neural network; Figure 11d is the prediction error histogram of the third data model, namely the support vector machine).

[0194] S04: Establish a fusion prediction model based on Kalman filtering. The change in tool wear over time can be considered a continuous random process. The result of S02 in step 6 is used as the prediction equation, and the result of S03 in step 6 is used as the observation equation. The combined values ​​of the two are used to accurately predict tool wear.

[0195] Define the predicted quantity of the tool wear physical model at time t as Y0, and the predicted quantity is y0; select one sensor to extract the feature data model and define the observed quantity as Y1, and the observed quantity is y1; define the probability density function as f. The posterior probability form of tool wear is expressed as:

[0196]

[0197] Where, It is called the likelihood function.

[0198] Select the normal distribution model to simplify the calculation of the likelihood function.

[0199] The time t is extended to the continuous time process to construct the prediction and observation equations. The prediction equation needs to find the recursive relationship between the physical model before and after the time, while the output value of the data model can be directly applied to the observation equation. Assuming that the recursive relationship between the time before and after the prediction equation is F, the result of a data model is used as the judgment value G(Y) of the observation equation at time t. 0,t ), the model noise of prediction and observation is expressed as Q t and R t , the prediction model and observation model in the continuous process are expressed as:

[0200]

[0201] Set the initial value to f0(y0) and make a prediction f1 - (y0) and update f1 + Probability density function calculation of (y0):

[0202]

[0203] Where, Subsequent f2 - ,f2 + ,f3 - ,f3 + ,… are calculated recursively through formula (15).

[0204] The specific updated value after each time t is directly solved through the probability density function:

[0205]

[0206] Based on this, we make two assumptions to solve the infinite integral. Assume that the prediction and observation equations are linear and the prediction and observation noises follow a normal distribution:

[0207]

[0208] Where, F c ,G c is a constant, and the assumption is substituted into the recursive process to obtain the KF calculation formula:

[0209]

[0210] Where, is the Kalman gain, are the predicted value and updated value of the physical model at time t, are the errors of the predicted value and the updated value at time t, y 1,t Output value of the data model at time t.

[0211] The prediction equation contains information about tool wear, wear velocity and wear acceleration, which can be rewritten in vector form: The error is rewritten as the covariance matrix σ t →∑ t The selected data model is extended from a single to multiple for tool wear prediction and rewritten as Y t =[y 1,t y 2,t y 3,t …]. F c ,G c ,Q t ,R t is a matrix constant, and the KF in matrix form is as follows:

[0212]

[0213] Where, is the Kalman gain, M I is the identity matrix.

[0214] According to the percentage of prediction error at different wear stages, the confidence level of prediction and observation is set. The final tool wear curve (such as Figure 12 is the fusion prediction result of tool C1. The area is the average of three flank surface measurements of tool C1. The solid line is the prediction result after KF fusion, and its smoothness is better than all the curves in Figure 10).

[0215] Example:

[0216] Step 1: Apply the analysis method of tool C1 to C4 and C6 to obtain the actual wear measurement value (such as Figure 13a is the measured value of the actual wear of tool C4; Figure 13b is the measured value of the actual wear of tool C6).

[0217] Step 2: Comparative analysis of tools C1, C4, and C6 revealed similar acquisition times and amplitudes in the time domain images, while the primary frequency bands corresponding to each sensor in the frequency domain images were roughly the same. The time duration selected for each experiment remained roughly constant at 4 seconds. The nine time domain features and six frequency domain features of the three tools exhibited similar trends. Wavelet packet analysis of the vibration signals also used nodes 5 and 13 based on range and variance. The heatmap correlation analysis showed similar results, indicating that the extracted features were appropriate.

[0218] Step 3: According to the first-order derivative, for tool C4 (such as Figure 14a It is divided into the C4 wear stage of the tool; Figure 14b is the first-order derivative of tool C4, and the dotted line indicates that the tool wear is divided into three stages by the first-order derivative) and C6 ( Figure 15aIt is divided into the C6 wear stage of the tool; Figure 15b is the first-order derivative of tool C6, and the dotted line indicates that the tool wear is divided into three stages by the first-order derivative. The wear stage is divided using the physical model and the three data models. Figure 16a C4 wear prediction for the physical model; Figure 16b The first data model is the C4 wear prediction of random forest; Figure 16c The second data model is the C4 wear prediction of the neural network; Figure 16d For the third data model, namely the support vector machine C4 wear prediction) and C6 (such as Figure 17a The wear amount prediction of C6 for the physical model; Figure 17b The first data model is the C6 wear prediction of random forest; Figure 17c The second data model is the C6 wear prediction of the neural network; Figure 17d The third data model is the C6 wear prediction of the support vector machine) wear prediction.

[0219] Step 4: Fusion of different models by KF to obtain the prediction results of tool C4 and C6 wear (e.g. Figure 18a is the fusion prediction result of tool C4. The area is the average of three flank surface measurements of tool C4. The solid line is the prediction result after KF fusion, and its smoothness is better than all the curves in Figure 16. Figure 18b is the fusion prediction result of tool C6. The area is the average of three flank surface measurements of tool C6. The solid line is the prediction result after KF fusion, and its smoothness is better than all the curves in Figure 17).

[0220] Step 5: Comparative analysis is performed on the predictions of tools C1, C4, and C6, using the average error percentage of the entire wear stage as the evaluation index. The comparison results are shown in Table 9.

[0221] Table 9 Average error percentage of different tools in the full wear stage

[0222] RF BPNN Support Vector Machine KF C1 0.005497 0.004302 0.036134 0.003995 C4 0.009674 0.010875 0.049578 0.007128 C6 0.010904 0.010178 0.048123 0.009453

[0223] Table 9 shows that, using the average percentage error as the evaluation criterion, RF performs best among the single prediction models, and RF and BPNN perform better than SVM. The proposed fusion method also achieves better prediction errors than the single prediction models. This comparison demonstrates the advantages of the present invention in measuring wear.

Claims

1. A tool condition monitoring method based on a fusion model of Kalman filtering is characterized by: The steps are as follows: 1) Collect the signal of tool processing and the wear of cutting edge flank: S01: Install a dynamometer on the fixture to measure the force signal F in the XYZ directions x ,F y ,F z ; Install vibration sensors in the tool feed direction, spindle radial direction and axial direction to measure the vibration signal V in the three directions of XYZ x ,V y ,V z ; Install an acoustic emission sensor S on the workpiece to measure the root mean square value of the acoustic emission. Then, seven sets of sensor signals are collected in each experiment, namely: F x ,、F y ,、F z ,、V x ,、V y ,、V z , and S; S02: After each cutting, the flank wear is measured three times and the average value is calculated; 2) Due to the existence of external environmental noise and internal factors of the machine tool, the collected signals need to be preprocessed: S01: Elimination of signal scenario trend items: Since the temperature rise of the machine tool itself and the temperature change of the external environment will cause the scene trend item, t is the sampling time, x is the acquisition signal, (t i ,x i ) represents the time t i Collect signal value x i , assuming that the number of sampling points in t is n, x can be fitted with a k-order polynomial: Where, α j is the degree of the polynomial with the highest degree k; If you remember The coefficients are obtained using the least squares method: Α=(t T t) -1 t T x (2) Where t is the coefficient matrix of sampling time, A is the sampling matrix of polynomial degree, and x is the value of the collected signal; S02: Elimination of signal zero drift: The initial stage of the signal contains an unprocessed state. Assuming that the first k sampling points are unprocessed state signals, the average value of the first k signals is moved to shift the signal to zero: Where s i is the preprocessed signal; 3) Use wavelet denoising method to remove noise from the preprocessed force signal: S01: Performing wavelet decomposition and low-frequency reconstruction on the signal can obtain a signal with less high-frequency noise. The wavelet basis function ψ and the decomposition scale m will affect the denoising effect. Based on the signal-to-noise ratio, root mean square error and smoothness, the entropy method is used to construct a unified evaluation standard for the fusion denoising index. Where SNR is the signal-to-noise ratio, RMSE is the root mean square error, S is the smoothness, and s real =[s real,1 s real,2 …s real,n ] is an ideal signal, is the denoised signal; Each indicator value is normalized to the range of 0.1 to 1: Where, I represents the normalized denoising index; Use the entropy method to determine the weight of each indicator: Where P is the probability, E is the entropy value, W is the weight, and C is a constant set to the total number of participating ranking indicators; Combining equations (5) and (6) to construct the fusion denoising index T under different wavelet parameters: T(ψ,m)=W SNR ×I SNR (ψ,m)+W RMSE ×I RMSE (ψ,m)+W S ×I S (ψ,m) (7) S02: Since the ideal signal cannot be obtained in the actual process, it is necessary to construct a simulated signal, add noise, and then denoise it to calculate the optimal noise reduction parameters; A periodic signal can be viewed as a superposition of different sinusoidal signals. In the frequency domain, the signal frequency corresponding to when the energy exceeds a specified threshold is selected, and the frequency, amplitude, and phase of the sinusoidal signal are calculated, and the simulation signal is constructed by superposition. Gaussian white noise is added to the simulation signal to simulate the noise under the actual processing state. Since the noise generated each time is random, a noisy signal of multiple simulation signals is constructed and the denoising effect is measured by the mean of the fusion index. 4) Use the smoothing method to remove noise from the preprocessed vibration signal: Polynomial approximation is used in the smoothing interval, and the number of signal points in the smoothing interval is used to solve the unknown parameters. The smoothing window is set to 5, and the number of polynomial approximations is set to 3. This is also called the 5-point 3-times mean smoothing method. The calculation formula is as follows: 5) Extract features from the processed signal; S01: Unify the experimental sampling time, assuming that the signal sampling time is t0 and the retention time is t * , the starting sampling point of the retention time is t start , the end sampling point is t end , the sampling frequency is ω, then the starting and ending sampling points after deletion are set to: S02: Time domain feature analysis method: Time domain analysis is performed on multiple experiments of 7 sets of tool sensor signals; including mean, root mean square value, variance, skewness, kurtosis, peak index, waveform index, pulse index and margin index; S03: Frequency Domain Characteristic Analysis Method: The relationship between frequency ω and energy p(ω) is calculated through Fourier transform. Frequency domain analysis is performed on multiple experiments of force and acoustic emission in the seven sets of tool sensor signals. This includes center of gravity frequency, mean square frequency, root mean square frequency, frequency standard deviation, frequency variance, and total frequency energy. S04: Time-frequency feature wavelet packet analysis: The vibration signal has concentrated energy at both low and high frequencies. A 4-layer wavelet decomposition was used to decompose the vibration signal into 16 frequency bands. The range and variance of the energy of each frequency band were calculated for all experimental groups. The wavelet packet nodes with the largest range and variance values ​​were selected, and the frequency domain feature analysis of S03 was then performed. S05: Feature Dimensionality Reduction: Calculate the autocorrelation between different features of the same sensor, retain the most representative features to achieve feature dimensionality reduction, and use the Pearson correlation coefficient to calculate the correlation analysis as follows: Where λ i ,λ j are the calculation results of the i-th and j-th feature analysis methods respectively; |r ij The closer | is to 1, the i ,λ j The higher the degree of relevance; After the correlation calculation, the t test was used to determine the significance; After completing the autocorrelation analysis, the most representative features are retained and the features are subjected to dimensionality reduction. The correlation between the dimensionality reduction data and the tool wear status is then established, and highly correlated features are selected for a second dimensionality reduction. For some sensors, all extracted features may have a low correlation with tool wear. In this case, at least one feature from each sensor is selected as the input of the prediction model. 6) Establish a fusion prediction model for tool wear status: Establish a physical model containing a time-recursive relationship as the prediction equation, select a suitable data model to substitute the extracted features to predict tool wear, and directly apply the wear value to the observation equation, completing the fusion of multiple models to reduce errors; S01: Classification of tool wear stages: Tool wear is divided into three different wear stages, and a different prediction model is established for each stage. Tool wear conditions are converted into first-order derivatives that vary with time to determine how to classify tool wear states. A decision tree is selected as the classification prediction model, with the final features extracted in step 5) as input and the tool wear classification stages as output, enabling rapid and accurate classification of tool wear stages. S02: Establish a physical model for tool wear: The initial, mid-term, and late-term tool wear satisfy different prediction equations, and wear failure follows a Weibull distribution. The demarcation point of the three wear states is determined by the derivative image of the tool wear time domain image. It is assumed that the initial wear change rate decreases with time; the mid-term wear change rate remains constant and non-zero with time; and the late-term wear change rate increases with time. The Weibull distribution of initial and late-term wear is assumed to be uniformly accelerated linear motion according to the change law: The mid-term wear Weibull distribution is assumed to be uniform linear motion according to the law of change: It is only necessary to measure the wear amount, wear velocity and initial value of wear acceleration in each wear period to establish the recursive relationship expression of the physical model of the entire cycle; S03: Establishing a data model for tool wear: Establishing a connection with tool wear through the extracted optimal features; Select multiple data prediction models with relatively simple structures to quickly predict the general trend; for each data model, calculate the error between the predicted wear amount and the actual wear amount; S04: Establish a fusion prediction model based on Kalman filtering: The change of tool wear over time can be regarded as a continuous random process. The result of S02 in step 6) is used as the prediction equation, and the result of S03 in step 6) is used as the observation equation. The numerical values ​​of the two are combined to accurately predict the tool wear; Define the predicted quantity of the tool wear physical model at time t as Y0, and the predicted quantity is y0; select a sensor to extract the feature data model and define the observed quantity as Y1, and the observed quantity is y1; define the probability density function as f; then the posterior probability of tool wear is expressed as: Where, It is called the likelihood function; Select the normal distribution model to simplify the calculation of the likelihood function: The time t is extended to the continuous time process, and the prediction and observation equations are constructed. The prediction equation needs to find the recursive relationship between the physical model before and after the time, while the output value of the data model can be directly applied to the observation equation. Assuming that the recursive relationship between the time before and after the prediction equation is F, the result of a data model is used as the judgment value G(Y) of the observation equation at time t. 0,t ), the model noise of prediction and observation is expressed as Q t and R t , the prediction model and observation model in the continuous process are expressed as: Set the initial value to f0(y0) and make a prediction f1 - (y0) and update f1 + Probability density function calculation of (y0): Where, Follow-up The calculation of is recursively deduced through formula (15); The specific updated value after each time t is directly solved through the probability density function: On this basis, two assumptions are made to solve the infinite integral; the prediction and observation equations are assumed to be linear, and the prediction and observation noises follow a normal distribution: Where, F c ,G c is a constant, and the assumption is substituted into the recursive process to obtain the KF calculation formula: Where, is the Kalman gain, are the predicted value and updated value of the physical model at time t, are the errors of the predicted value and the updated value at time t, y 1,t Output value of the data model at time t; The prediction equation contains information about tool wear, wear velocity and wear acceleration, which can be rewritten in vector form: The error is rewritten as the covariance matrix σ t →∑ t ; The selected data model is extended from a single to multiple for tool wear prediction and rewritten as Y t =[y 1,t y 2,t y 3,t …]; F c ,G c ,Q t ,R t is a matrix constant, and the KF in matrix form is as follows: Where, is the Kalman gain, M I is the identity matrix; Setting confidence levels for predictions and observations based on the percentage of prediction error at different wear stages; The final tool wear curve is obtained.

Images