A monitoring and fault prediction system for CNC machining centers based on multi-sensor fusion

Abstract

This invention relates to the field of CNC machining equipment monitoring technology, and discloses a CNC machining center monitoring and fault prediction system based on multi-sensor fusion. The system includes: a data acquisition module for acquiring data from lubricating oil circuit dielectric constant sensors, spindle temperature sensor array data, and tool tip three-dimensional position deviation data; a data synchronization module for synchronizing multi-source sensor data to a unified time reference; a feature extraction module for extracting high-frequency pulse disturbance features of the dielectric constant and position-temperature correlation parameters; a causal analysis module for constructing a multi-physical quantity causal topology graph and identifying the dominant causal path; a risk assessment module for calculating a lubrication contamination risk score; and an early warning decision module for generating targeted maintenance instructions. This invention solves the technical problem of existing technologies being unable to distinguish between lubrication contamination and anomalies caused by external heat sources.

Description

Technical Field

[0001] This invention relates to the field of CNC machining equipment monitoring technology, and more specifically, to a CNC machining center monitoring and fault prediction system based on multi-sensor fusion. Background Technology

[0002] In the field of precision CNC machining, the spindle tip position accuracy is a key factor affecting machining quality. During the operation of a CNC machining center, the axial position deviation of the spindle tip is affected by both lubrication and thermal conditions. When lubricating oil becomes contaminated with trace amounts of moisture due to condensation, the moisture catalyzes the oxidation of the lubricating oil, leading to increased bearing friction, increased heat generation, and spindle thermal expansion, resulting in an increased axial position deviation of the tool tip. Simultaneously, external heat sources (such as ambient temperature fluctuations and cutting heat accumulation) also cause the spindle temperature to rise and thermal expansion, similarly manifesting as an increased axial position deviation of the tool tip. While these two abnormal situations appear similar in terms of temperature and macroscopic manifestations of tool tip axial deviation, their physical mechanisms and maintenance requirements are completely different: the former requires lubricating oil dehydration or replacement, while the latter only requires adjustment of thermal compensation parameters.

[0003] Existing monitoring technologies rely on single-parameter monitoring methods and fixed-threshold judgment mechanisms for early warning based on temperature or position deviation thresholds. These technologies suffer from the following drawbacks: First, trace amounts of moisture (water content <0.1%) have less impact on the average dielectric constant of lubricating oil than on measurement accuracy; however, water droplets in oil flow cause high-frequency pulse disturbances in the dielectric constant, which are masked during the averaging process. Second, the three-dimensional position deviation of the cutting edge is a spatial coupling effect of thermal expansion and mechanical runout; undecoupled monitoring cannot distinguish between thermal and mechanical effects. Third, single-parameter threshold early warning systems cannot establish a causal relationship between lubrication status, temperature status, and position deviation, and therefore cannot distinguish the root cause of anomalies.

[0004] Misclassifying frictional heat caused by lubrication contamination as a normal operating heat source will delay lubrication maintenance and lead to bearing damage. Conversely, misclassifying thermal expansion caused by external heat sources as lubrication contamination will result in unnecessary lubrication oil changes. Misjudgments in maintenance strategies increase costs and fail to address the real problem. Therefore, a monitoring method for CNC machining centers that can integrate multi-sensor data, decouple multi-physical domain effects, and identify the root causes of anomalies is needed. Summary of the Invention

[0005] This invention provides a monitoring and fault prediction system for CNC machining centers based on multi-sensor fusion, which solves the technical problems in related technologies such as the inability to distinguish between abnormal spindle tip position deviation caused by lubrication contamination and external heat sources, increased costs due to misjudgment of maintenance strategies, and unresolved problems.

[0006] This invention provides a monitoring method for CNC machining centers based on multi-sensor fusion, comprising: Acquire high-sampling-rate raw data from the dielectric constant sensor of the lubrication circuit of the CNC machining center, data from the distributed temperature sensor array of the spindle, and three-dimensional position deviation data of the tool tip. Synchronize the multi-source sensor data to a unified time reference based on the timestamp alignment algorithm. High-pass filtering is performed on the dielectric constant time series to generate a high-frequency pulse perturbation sequence of dielectric constant. The pulse event frequency and average pulse amplitude within the sliding window are calculated. The skewness and kurtosis parameters of the pulse amplitude distribution are calculated using a probability density estimation algorithm to generate a lubrication contamination feature vector. The three-dimensional position deviation of the tool tip is decomposed into axial and radial components. The average temperature sequence of the spindle is calculated. Correlation analysis is used to calculate the linear correlation coefficient between axial deviation and temperature. The thermal elongation coefficient is calculated based on the linear regression relationship between axial deviation and temperature. A four-variable causal topological directed graph is constructed, consisting of lubrication contamination feature vector, temperature, axial deviation, and radial deviation. The Granger causality test is used to calculate the causal relationship strength of the three paths from lubrication contamination feature vector to temperature, temperature to axial deviation, and lubrication contamination feature vector to axial deviation. Based on the causal relationship strength, the dominant causal path is identified and the lubrication contamination weight and heat source weight are calculated. An anomaly detection boundary is established on healthy lubricating oil data using a single-class support vector machine, and the distance from the current lubricating contamination feature vector to the boundary is calculated to generate a lubricating contamination risk score. Determine whether the product of the lubrication contamination weight and the lubrication contamination risk score exceeds the lubrication maintenance threshold, or whether the heat source weight exceeds the heat compensation threshold, and generate a targeted early warning signal. The Granger causality test is performed as follows: For two standardized time series variables X and Y, an autoregressive model including X is established. The historical values ​​of Y and X are used as independent variables to predict the current value of Y. The sum of squared residuals of the complete model including X and the constrained model not including X is calculated. The F-statistic is calculated based on the difference between the sums of squared residuals of the two models. The value of the F-statistic is used as the strength of the causal relationship between X and Y.

[0007] Furthermore, a Butterworth high-pass filter is used for the high-pass filtering, with a cutoff frequency set to 100Hz and a filter order of 4. The method for identifying the pulse event is as follows: set an amplitude threshold, and identify a pulse event when the absolute value of the high-frequency disturbance sequence exceeds the amplitude threshold. Count the number of pulse events within the time window. The pulse event frequency is equal to the number of pulse events divided by the time window length, and the average pulse amplitude is equal to the arithmetic mean of the amplitudes of all pulse events within the window. The amplitude threshold is set to 3 times the standard deviation of the high-frequency perturbation sequence under healthy lubricating oil conditions, and the time window length is set to 60 seconds. The standard deviation under healthy lubricating oil conditions is calculated by continuously collecting high-frequency disturbance sequence data for 1 hour under stable operating conditions with new lubricating oil on the spindle, and using the standard deviation of the high-frequency disturbance sequence data as the baseline standard deviation of the healthy state.

[0008] Furthermore, the probability density estimation adopts the kernel density estimation method, with the input being a sequence of pulse amplitude samples, and the probability density function is estimated based on the Gaussian kernel function; Based on the estimated probability density function, the ratio of the third central moment to the cube of the standard deviation is calculated as the skewness parameter, and the ratio of the fourth central moment to the fourth power of the standard deviation is calculated and then 3 is subtracted as the kurtosis parameter. The lubrication contamination feature vector is a four-dimensional vector, including pulse event frequency, average pulse amplitude, skewness parameter, and kurtosis parameter; Z-score standardization is performed on the lubrication contamination feature vector to convert each component into a standardized feature with zero mean and unit variance.

[0009] Furthermore, the decomposition method for the three-dimensional positional deviation is as follows: Let the components of the three-dimensional position deviation of the tool tip measured by the laser rangefinder in the Cartesian coordinate system be the X-axis component, the Y-axis component, and the Z-axis component. The Z-axis is parallel to the main axis. The axial component is directly taken as the Z-axis component. The radial component is calculated as the square root of the sum of the squares of the X-axis component and the squares of the Y-axis component. The calculation method for the spindle average temperature sequence is as follows: the arithmetic average of the measurement values ​​of multiple temperature sensors deployed on the spindle at the same time is performed. The linear correlation coefficient between the axial deviation and temperature is the Pearson correlation coefficient, calculated by dividing the covariance of the axial deviation sequence and the average temperature sequence of the main shaft by the square root of the product of the variance of the axial deviation and the variance of the temperature.

[0010] Furthermore, the method for calculating the coefficient of thermal elongation is as follows: A linear regression model is established, in which the axial deviation is expressed as the thermal elongation coefficient multiplied by the average temperature of the main shaft, plus a constant term and a residual term. The thermal elongation coefficient and the constant term are estimated using the least squares method. The estimated value of the thermal elongation coefficient is equal to the covariance of the axial deviation and the temperature divided by the variance of the temperature. The temperature-independent axial runout residual sequence is equal to the measured axial deviation minus the estimated value of the thermal elongation coefficient multiplied by the average temperature of the spindle and then minus the estimated value of the constant term; The training samples were collected during the complete thermal equilibrium process of the spindle under normal operating conditions, with a data collection time of at least 2 hours to cover the temperature change range of the spindle from cold start to hot stability.

[0011] Furthermore, the Granger causality test standardizes the time series of each variable before execution to eliminate dimensional differences; The method for calculating the F-statistic is as follows: establish a constrained model without independent variable X and a complete model with X, calculate the sum of squared residuals of the two models respectively, and the F-statistic is equal to the difference between the sum of squared residuals of the constrained model and the complete model divided by the number of constraints, then divided by the sum of squared residuals of the complete model divided by the degrees of freedom. The lag order of the autoregressive model is determined by the Akaike Information Criterion. Different combinations of lag orders are traversed within the candidate range, and the combination of orders that minimizes the value of the Akaike Information Criterion is selected as the optimal lag order.

[0012] Furthermore, the rules for identifying the dominant causal path are as follows: Let the causal strength from the lubrication contamination feature vector to temperature be the first causal strength, the causal strength from temperature to axial deviation be the second causal strength, and the causal strength from the lubrication contamination feature vector to axial deviation be the third causal strength. If the first causal strength is greater than the causal significance threshold and the third causal strength is greater than the causal significance threshold, then it is determined to be the lubrication contamination dominant mode. The lubrication contamination weight is equal to the product of the first causal strength and the third causal strength divided by the product of the first causal strength and the third causal strength plus the second causal strength. If the first causal strength is less than the causal significance threshold and the second causal strength is greater than the causal significance threshold, it is determined to be an external heat source dominant mode. The heat source weight is equal to the second causal strength divided by the product of the first causal strength and the third causal strength plus the second causal strength. The causal significance threshold is set as the critical value of the F distribution at a significance level of 0.05.

[0013] Furthermore, the training process for a single-class support vector machine is as follows: Collect a sample set of lubrication contamination feature vectors of the spindle under healthy lubrication oil conditions, and use a single-class support vector machine algorithm to learn the decision boundary of the healthy state. The optimization objective is to minimize the square of the norm of the weight vector divided by 2 minus the regularization term of the bias parameter plus the relaxation variable. The constraint is that the inner product of the feature vectors after the weight vector is mapped to the kernel is greater than or equal to the bias parameter minus the relaxation variable and the relaxation variable is non-negative. The single-class support vector machine uses a radial basis function kernel. The kernel function is equal to an exponential function multiplied by the square of the Euclidean distance between two feature vectors with a negative bandwidth parameter. The bandwidth parameter is set to the reciprocal of the feature vector dimension. The method for calculating the lubrication contamination risk score is as follows: input the current lubrication contamination feature vector into the decision function to calculate the decision function value. The lubrication contamination risk score is equal to the larger of 0 and the negative value of the decision function. When the decision function value is less than 0, it indicates that the current feature vector is outside the health boundary.

[0014] Furthermore, the method for generating targeted early warning signals is as follows: The lubrication maintenance joint index equals the lubrication contamination weight multiplied by the lubrication contamination risk score; If the combined lubrication maintenance index is greater than the lubrication maintenance threshold, a lubricating oil dehydration treatment instruction or a replacement instruction will be generated. If the lubrication contamination risk score is less than the replacement threshold, a dehydration treatment instruction will be output. If the lubrication contamination risk score is greater than or equal to the replacement threshold, a lubricating oil replacement instruction will be output. If the heat source weight is greater than the heat compensation threshold, the heat expansion is calculated based on the heat expansion coefficient and the current spindle temperature. The heat expansion is equal to the heat expansion coefficient multiplied by the difference between the current spindle temperature and the reference temperature. The heat compensation parameter output is the negative value of the heat expansion as the axial position correction amount. The lubrication maintenance threshold and heat compensation threshold adopt an adaptive adjustment mechanism. Based on the quantile of the joint index of historical data and the weight of heat source within the most recent time window, the lubrication maintenance threshold is set to the 95th quantile of historical data, and the heat compensation threshold is set to the 90th quantile of historical data.

[0015] This invention provides a fault prediction system for CNC machining centers based on multi-sensor fusion, comprising: The data acquisition module is used to acquire data from the dielectric constant sensor of the lubricating oil circuit, data from the spindle temperature sensor array, and data from the three-dimensional position deviation of the tool tip; The data synchronization module is used to synchronize multi-source sensor data to a unified time reference; The feature extraction module is used to extract high-frequency pulse perturbation features of dielectric constant and position-temperature correlation parameters; The causal analysis module is used to construct a causal topology graph of multiple physical quantities and identify the dominant causal path; The risk assessment module is used to calculate the lubrication contamination risk score; The early warning and decision-making module is used to generate targeted maintenance instructions.

[0016] The beneficial effects of this invention are as follows: This invention establishes a collaborative sensing mechanism across three physical domains—lubrication, thermal, and mechanical—by integrating data from three types of multi-source sensors: dielectric constant sensors, temperature sensor arrays, and position deviation sensors. This overcomes the technical limitation of single sensors failing to comprehensively reflect equipment status. High-pass filtering extracts high-frequency pulse disturbance characteristics of the dielectric constant, constructing a four-dimensional feature vector containing pulse event frequency, average pulse amplitude, skewness parameter, and kurtosis parameter. This overcomes the technical limitation of average value calculation masking pulse disturbances caused by trace moisture, enabling early detection at the trace moisture (water content <0.1%) stage. Axial-radial decomposition of the blade tip position deviation, temperature correlation analysis, and linear regression are used to calculate the thermal elongation coefficient. This overcomes the technical limitation of the inability to distinguish the physical mechanism of three-dimensional position deviation due to the coupling effect of thermal elongation and mechanical runout, achieving decoupling of thermal and mechanical effects.

[0017] This invention utilizes the Granger causality test to construct a four-variable causal topology graph of lubrication contamination feature vector, temperature, axial deviation, and radial deviation. It quantifies the strength of two causal paths: lubrication contamination → temperature → axial deviation and external heat source → temperature → axial deviation. It identifies the dominant causal path and calculates the weights of lubrication contamination and heat sources, overcoming the technical limitations of single-parameter threshold warnings that cannot establish causal relationships between physical quantities or distinguish abnormal root causes. By using a single-class support vector machine to establish anomaly detection boundaries for healthy lubricating oil conditions and calculating a lubrication contamination risk score, it combines the lubrication contamination weight and risk score to form a joint index. This enables separate decision-making for lubrication maintenance and thermal compensation, avoiding misjudging lubrication contamination as an external heat source or vice versa, leading to misjudgments in maintenance strategies.

[0018] This invention solves the technical problems in the prior art, such as the detection failure caused by the influence of trace moisture on the average dielectric constant being less than that of measurement accuracy, the inability to distinguish the mechanism caused by the coupling effect of tool tip position deviation, and the inability to identify the root cause of abnormalities caused by single threshold warning. It achieves the technical effects of realizing early detection at the trace moisture stage, accurately distinguishing lubrication contamination from external heat sources, and generating targeted maintenance instructions to avoid the cost of misjudgment. Attached Figure Description

[0019] Figure 1 This is a flowchart of a CNC machining center monitoring method based on multi-sensor fusion according to the present invention; Figure 2 This is a multi-Y-axis piecewise linear plot of the time series of multi-sensor synchronized data from the present invention; Figure 3 This is a mixed plot of line and scatter plots showing the high-frequency pulse perturbation feature extraction of dielectric constant in this invention; Figure 4 This is a mixed plot of scatter points and regression lines from the linear regression analysis of axial deviation and temperature of the present invention; Figure 5 This invention provides a hybrid graph of line graph and threshold bar graph for time-series monitoring and threshold judgment of early warning indicators. Detailed Implementation

[0020] The subject matter described herein will now be discussed with reference to exemplary embodiments. It should be understood that these embodiments are discussed only to enable those skilled in the art to better understand and implement the subject matter described herein, and changes may be made to the function and arrangement of the elements discussed without departing from the scope of this specification. Various processes or components may be omitted, substituted, or added as needed in the examples. Furthermore, some features described in the examples may be combined in other examples.

[0021] At least one embodiment of the present invention discloses a monitoring method for CNC machining centers based on multi-sensor fusion, such as... Figure 1 As shown, it includes the following steps: Step 100: Acquire multi-sensor synchronization data and generate a multi-physical quantity fusion time series dataset; Acquire high-sampling-rate raw data (sampling rate ≥ 1kHz) from the dielectric constant sensor in the lubrication circuit of the CNC machining center, data from the distributed temperature sensor array on the spindle, and three-dimensional position deviation data of the tool tip. Note that the sampling rates of different sensors are typically inconsistent; the dielectric constant sensor has a sampling rate above 1kHz, the temperature sensor typically has a sampling rate of 10Hz, and the position sensor typically has a sampling rate of 100Hz. Based on a timestamp alignment algorithm, the multi-source sensor data with different sampling rates are synchronized to a unified time base, generating a synchronized time-series dataset that fuses multiple physical quantities.

[0022] It should be noted that the timestamp alignment algorithm described above uses a linear interpolation method. The input consists of multi-source sensor data with different sampling rates and their timestamp sequences, and the output is a synchronized dataset under a unified time reference. The specific execution steps are as follows: the timestamp sequence of the sensor with the lowest sampling rate is selected as the reference time series, and linear interpolation is performed on the other high-sampling-rate sensor data at the reference timestamp to generate the values ​​of each sensor under a unified timestamp.

[0023] Furthermore, the calculation method for the above linear interpolation is as follows: for the reference timestamp Find timestamps in high sampling rate sensor data that satisfy The two adjacent data points have corresponding measured values. and Then in the timestamp The interpolation result at point is calculated as follows: .

[0024] When the CNC machining center in a precision parts processing workshop started operating at 10:15 AM, the monitoring system simultaneously collected data from three types of sensors. A dielectric constant sensor continuously recorded the dielectric constant of the lubrication oil circuit at a sampling rate of 1000 Hz. Four temperature sensors deployed on the spindle measured the temperatures of the front bearing, rear bearing, motor end, and spindle housing at a sampling rate of 10 Hz. A laser rangefinder measured the positional deviation of the tool tip in the X, Y, and Z directions at a sampling rate of 100 Hz. The spindle was currently rotating at 3000 rpm, machining aluminum alloy. Table 1 shows the raw sampling data segments of each sensor from 10:15:00.0 to 10:15:00.2.

[0025] Table 1. Raw sampling data from multiple sensors (partial): To achieve time synchronization of multi-sensor data, the temperature sensor with the lowest sampling rate (10Hz, sampling interval 0.1s) was selected as the reference time series. Linear interpolation was performed on the dielectric constant data and position data at the reference timestamp. Taking s as an example, the measured value of the dielectric constant at that moment is 2.154, so no interpolation is needed; position data exist The value at point s is -0.5 μm. The value at point s is -0.8 μm, and the interpolation calculation is as follows: μm. Table 2 shows the fused dataset after synchronization and alignment.

[0026] Table 2. Multi-physical quantity fusion dataset after time synchronization (partial): Figure 2 The results of synchronous data acquisition from multiple sensors are shown.

[0027] Step 200: Extract the high-frequency pulse perturbation features of dielectric constant to generate a lubrication contamination feature vector; High-pass filtering is performed on the dielectric constant time series obtained in step 100 to remove DC components and low-frequency drift, generating a high-frequency pulse perturbation sequence of the dielectric constant. The pulse event frequency and average pulse amplitude within the sliding window are calculated for the high-frequency perturbation sequence. The skewness and kurtosis parameters of the pulse amplitude distribution are calculated using a probability density estimation algorithm to generate a lubrication contamination feature vector.

[0028] It should be noted that the high-pass filter described above uses a Butterworth high-pass filter with a cutoff frequency set to 100Hz and a filter order of 4. This filter is used to remove the DC component of the dielectric constant and low-frequency drift below 100Hz. The filtered signal is then a high-frequency pulse perturbation sequence of the dielectric constant.

[0029] It should be noted that the method for identifying the aforementioned pulse events is as follows: setting an amplitude threshold. When the absolute value of the high-frequency perturbation sequence exceeds This is identified as a single pulse event within the time window. Internal statistical pulse event count The pulse event frequency is calculated as follows: The average pulse amplitude is calculated as the arithmetic mean of the amplitudes of all pulse events within the window. ,in For the summation index, For the number of pulse events, For the first The amplitude of the secondary pulse event.

[0030] Furthermore, the aforementioned amplitude threshold The time window length is set to 3 times the standard deviation of the high-frequency perturbation sequence under healthy lubricating oil conditions. A 60-second time window ensures that the time window includes enough pulse event samples to guarantee the stability of statistical characteristics while also responding promptly to changes in lubrication status. The standard deviation under healthy lubrication conditions is calculated by continuously collecting high-frequency disturbance sequence data for one hour under stable operating conditions with new lubricating oil on the spindle. This standard deviation is used as the baseline standard deviation for the healthy state.

[0031] It should be noted that the above probability density estimation uses the kernel density estimation method, and the input is a pulse amplitude sample sequence. The output consists of skewness and kurtosis parameters, where The amplitude of the first pulse event. The amplitude of the second pulse event. The amplitude of the Nth pulse event. Let be the total number of impulse events. The probability density function is estimated as follows: ,in For the summation index, For Gaussian kernel function, The bandwidth parameter is used. The skewness parameter is calculated based on the estimated probability density function. and kurtosis parameter ,in The mean, Standard deviation This is the expectation operator. The lubrication contamination feature vector is a four-dimensional vector, represented as... ,in This indicates transpose.

[0032] Furthermore, the aforementioned expectation operator The calculation is estimated through sample statistics, mean. Standard deviation Third central moment Fourth-order central moments ,in For the summation index, the sample sequence is: .

[0033] Furthermore, the aforementioned bandwidth parameters The calculation is performed using Silverman's empirical rule, and the formula is as follows: ,in The standard deviation of the sample is 1. The interquartile range of the sample. This represents the number of samples.

[0034] Furthermore, the aforementioned interquartile range The calculation method is as follows: [The sample sequence is used for calculation.] Sort in ascending order and calculate the 25th percentile. and the 75th percentile Interquartile range .

[0035] It should be noted that the components of the above-mentioned lubrication contamination feature vector have different dimensions, and the pulse event frequency... The dimension of the average pulse amplitude is frequency. The dimension of is dielectric constant, while the skewness and kurtosis parameters are dimensionless. To eliminate the influence of dimensions on subsequent distance calculations, Z-score standardization is performed on the lubrication contamination feature vector, converting each component into a standardized feature with zero mean and unit variance.

[0036] A 4th-order Butterworth high-pass filter was applied to 60 seconds of dielectric constant data collected between 10:15:00 and 10:16:00, with a cutoff frequency of 100Hz. The standard deviation of high-frequency disturbances under healthy lubricating oil conditions was obtained from historical benchmark data and was 0.0008; therefore, the pulse amplitude threshold was... The filtered high-frequency perturbation sequence identified pulse events within a 60-second window. Table 3 shows the identification results of some pulse events.

[0037] Table 3. Results of High-Frequency Pulse Event Identification for Dielectric Constant (Partial): Within a 60-second window, 18 pulse events were identified. The frequency of these pulse events... 18 pulses per second. The amplitude sequence of 18 pulse events is as follows: Average pulse amplitude .

[0038] Calculate the mean of this amplitude sequence. Standard deviation Third-order central moments Fourth-order central moments Obtain the skewness parameter and kurtosis parameter .

[0039] The original lubrication contamination feature vector is .

[0040] Z-score normalization is performed on the feature vector using the mean vector of the healthy state sample set. and standard deviation vector To obtain the standardized feature vector .

[0041] Figure 3 The process of extracting the dielectric constant by high-frequency pulse perturbation is demonstrated.

[0042] Step 300: Decompose the tip position deviation and extract thermally relevant features to generate position-temperature correlation parameters; The three-dimensional position deviation of the tool tip obtained in step 100 is decomposed into axial and radial components. The average temperature sequence of the spindle is calculated. Correlation analysis is used to calculate the linear correlation coefficient between axial deviation and temperature. Based on the linear regression relationship between axial deviation and temperature, the thermal elongation coefficient is calculated, and the temperature-independent axial runout residual is extracted. A Fourier transform is performed on the radial deviation to extract the radial runout amplitude synchronized with the frequency.

[0043] It should be noted that the decomposition method of the above three-dimensional position deviation is as follows: Let the component of the three-dimensional position deviation of the blade tip measured by the laser rangefinder in the Cartesian coordinate system be... , , The Z-axis is parallel to the main spindle axis. The axial component is directly taken as the Z-axis component. The radial component is calculated as the deviation amplitude in the XY plane. .

[0044] It should be noted that the method for calculating the linear correlation coefficient between the axial deviation and temperature is as follows: Calculate the axial deviation sequence With spindle average temperature sequence Pearson correlation coefficient The closer the absolute value of the linear correlation coefficient is to 1, the stronger the linear correlation between axial deviation and temperature. The linear correlation coefficient is used to verify the existence of thermal elongation effect.

[0045] It should be noted that the above calculation method for the spindle average temperature sequence is as follows: assuming that there are [missing information - likely related to equipment or devices] deployed on the spindle. A temperature sensor, at time The measured values ​​of each sensor are The average spindle temperature is calculated as follows: ,in For the first temperature sensor at time The measured value, For the second temperature sensor at time The measured value, For the Mth temperature sensor at time... The measured value, For the number of temperature sensors, For summation index.

[0046] It should be noted that the input of the above Fourier transform is the radial deviation time series, and the output is the spectrum of the radial deviation. The amplitude at the frequency corresponding to the spindle rotation frequency is extracted from the spectrum as the radial jump amplitude of the rotation frequency synchronization.

[0047] Furthermore, the aforementioned spindle rotation frequency is read in real time by the spindle encoder of the CNC system. The spindle rotation frequency corresponds to the spindle speed divided by 60. When the spindle speed is 3000 rpm, the corresponding frequency is 50 Hz.

[0048] Furthermore, the method for extracting the amplitude at a specific frequency in the aforementioned spectrum is as follows: perform a discrete Fourier transform on the radial deviation time series to obtain a complex frequency domain sequence, and calculate the frequency index corresponding to the spindle rotation frequency. ,in For frequency index, Main spindle frequency The number of points in the Fourier transform. The sampling frequency is used to extract the complex amplitude at the frequency index corresponding to the spindle rotation frequency as the radial runout amplitude for frequency synchronization.

[0049] It should be noted that the above method for calculating the thermal elongation coefficient is as follows: Let the axial deviation sequence be... The principal shaft average temperature sequence is Establish a linear regression model ,in The coefficient of thermal elongation, For constant terms, This is the residual term. The parameters are estimated using the least squares method. and Thermal elongation coefficient The estimated value is ,in For covariance, The variance is given. The temperature-independent axial runout residual sequence is... .

[0050] Furthermore, the methods for calculating the covariance and variance mentioned above are as follows: Let the sample of the axial deviation sequence be... The temperature sequence sample is The covariance is calculated as follows: The variance is calculated as follows: ,in For the summation index, The sample size used to calculate covariance. For a moment The axial deviation measurement value, For a moment Spindle average temperature and These are the mean values ​​of axial deviation and temperature, respectively. The estimated value of the constant term is... .

[0051] The input to the aforementioned linear regression model is the axis-mean temperature series. The output is the predicted value of axial deviation. The training process employs a supervised learning model, using historical data of axial deviation and temperature sequences as training samples. The least squares method is used to solve for the mean square error loss function. Minimize parameters and ,in For the summation index, This represents the number of training samples for linear regression.

[0052] Furthermore, the time index in the above loss function The sampling time corresponding to the training samples, the number of samples The training time window length divided by the sampling interval; when the training duration is 2 hours and the sampling interval is 10 seconds, the number of samples is... .

[0053] Furthermore, the training samples collected during the complete thermal equilibrium process of the spindle under normal operating conditions were taken, with a data collection time of at least 2 hours, to cover the temperature change range of the spindle from cold start to hot stability, ensuring that the regression model can accurately characterize the linear relationship between temperature and axial deviation.

[0054] Using the synchronized three-dimensional position data in Table 2, the tool tip position deviation is decomposed into axial and radial components.

[0055] The average temperature of the four temperature sensors during the period from 10:15:00 to 10:16:00 was calculated as follows: Table 4 shows the location-temperature correlation data for this period.

[0056] Table 4. Decomposition of blade tip position deviation and temperature correlation data (partial): The Pearson correlation coefficient was calculated based on temperature and axial deviation data within a 60-second window. Mean axial deviation. μm, variance Average temperature ℃, variance Covariance Correlation coefficient This indicates that there is a strong linear correlation between axial deviation and temperature.

[0057] Estimating the coefficient of thermal elongation using the least squares method μm / ℃, constant term μm. The linear regression model is... The axial runout residual sequence was calculated, and the root mean square of the residuals was 0.28 μm, indicating a good model fit. A discrete Fourier transform was performed on the radial deviation. A spindle speed of 3000 rpm corresponds to a rotational frequency of 50 Hz, and the amplitude of the spectrum at 50 Hz is 1.57 μm, which is the radial runout amplitude at frequency synchronization.

[0058] Figure 4 The linear correlation analysis between the tool tip axial deviation and the average spindle temperature is shown.

[0059] Step 400: Construct a multi-physical quantity causal topology graph using the Granger causality test to identify the dominant causal path; A four-variable causal topological directed graph is constructed, with the four variables being: the lubrication contamination feature vector generated in step 200, the principal axis average temperature calculated in step 300, the axial deviation obtained in step 300, and the radial deviation obtained in step 300. Granger causality tests are used to calculate the causal strength of the three paths: lubrication contamination feature vector → temperature, temperature → axial deviation, and lubrication contamination feature vector → axial deviation. Dominant causal paths are identified based on the causal topological graph: if the lubrication contamination → temperature → axial deviation path is significant, it is determined to be the lubrication contamination dominant mode; if only the temperature → axial deviation path is significant, it is determined to be the external heat source dominant mode. The causal strength of each path is calculated as the lubrication contamination weight and the heat source weight.

[0060] It should be noted that the four variables mentioned above have different dimensions. The lubrication contamination feature vector has been standardized in step 200. The spindle average temperature, axial deviation, and radial deviation need to be Z-score standardized before performing the Granger causality test to eliminate the influence of the difference in dimensions on the addition of different terms in the autoregressive model.

[0061] It should be noted that the inputs to the Granger causality test described above are two standardized time series variables. and The output is the strength of the causal relationship. The specific execution method is as follows: Establish an autoregressive model ,in and For the summation index, and The lag order is... and For regression coefficients, For the residual. If the coefficient If it is significantly non-zero (as determined by the F-test), then it is considered... right Granger causality exists. The strength of causality is quantitatively represented by the F-statistic; a larger F-statistic indicates a stronger causal relationship.

[0062] Furthermore, the method for calculating the F-statistic mentioned above is as follows: Establish without including constraint model and contain The complete model is obtained, and the sum of squared residuals for both models is calculated separately. and The F-statistic is calculated as follows: ,in For the summation index, The lag order is... Let be the sample size for the Granger causality test. The number of constraints.

[0063] Furthermore, the calculation of the residual sum of squares mentioned above involves the time dimension, and the residual sum of squares of the constrained model is... The sum of squared residuals of the complete model is: in For time indexing, the time index starts from... Begin by ensuring that all lags have valid data. represents the sample size for the Granger causality test.

[0064] Furthermore, the aforementioned lag order and The Akaike Information Criterion is used to determine the lag order. The specific method is as follows: traverse different combinations of lag orders within the candidate range of 1 to 10, and calculate the Akaike Information Criterion value corresponding to each combination. ,in The number of model parameters, The order of the likelihood function is selected as the optimal lag order, which minimizes the value of the Akaike information criterion.

[0065] Furthermore, the above likelihood function values The calculation method under the Gaussian assumption is as follows: ,in Represents an exponential function. The estimated value of the residual variance , The sum of squared residuals of the model. Let be the sample size for the Granger causality test. This represents the number of model parameters.

[0066] It should be noted that the identification rule for the above-mentioned dominant causal path is as follows: Let the lubrication contamination feature vector be... The temperature is axial deviation is Path 1 → The causal strength of " is denoted as Path Two → The causal strength of " is denoted as Path Three → The causal strength of " is denoted as .like and (in If the threshold for causal significance is reached, then the mode is determined to be dominated by lubrication contamination, and the lubrication contamination weight is calculated as follows: .like and If the external heat source is dominant, the heat source weight is calculated as follows: .

[0067] Furthermore, the aforementioned causal significance threshold The critical value for the F-distribution at a significance level of 0.05 is set. The critical value corresponding to the causal significance threshold is determined by looking up the F-distribution table based on the degrees of freedom of the autoregressive model. In the weight calculation formula, when the denominator... When the value is less than 0.01, the denominator is set to 0.01 to avoid instability in numerical calculation.

[0068] A four-variable causal topology graph was constructed, with variables including the standardized lubrication contamination feature vector (taking the first principal component as a scalar), principal shaft mean temperature, axial deviation, and radial deviation. An augmented Dickey-Fuller test was performed on 5 minutes of historical data from 10:10:00 to 10:15:00. The ADF statistic p-values ​​for all variables were less than 0.05, indicating stationarity and no differencing was required. The optimal lag order was determined using the Akaike information criterion. , Perform Granger causality tests to calculate the causal strength of the three critical paths. Table 5 shows the results of the Granger causality tests.

[0069] Table 5. Granger causality test results and weight calculations: Causal significance threshold The F-distribution value, determined by looking up the degrees of freedom, is 4.26 (significance level 0.05). The causal strength of path one, "lubrication contamination → temperature". The causal strength of path two, "temperature → axial deviation". The causal strength of path three, "lubrication contamination → axial deviation". .because and This was determined to be the dominant lubrication contamination mode. The lubrication contamination weight was then calculated. Heat source weight .

[0070] Step 500: Calculate the lubrication contamination risk score and generate a comprehensive early warning index; An anomaly detection boundary is established on healthy lubricating oil data using a single-class support vector machine. The distance from the current lubrication contamination feature vector generated in step 200 to the boundary is calculated to generate a lubrication contamination risk score.

[0071] It should be noted that the input to the aforementioned single-class support vector machine is the lubrication contamination feature vector. The output is the decision function value. The data transfer process of a single-class support vector machine includes a kernel mapping unit and a linear decision unit. The kernel mapping unit passes the input feature vector through a kernel function. Mapped to a high-dimensional feature space, the linear decision unit calculates the distance from the input sample to the decision boundary in the high-dimensional space.

[0072] It should be noted that the above kernel function uses a radial basis function kernel, and the specific form of the kernel function is as follows: ,in Represents an exponential function. The squared Euclidean distance between two eigenvectors is given. This is the bandwidth parameter of the kernel function. Kernel mapping. The input feature vector is implicitly mapped from the original four-dimensional space to a high-dimensional feature space using the radial basis function kernel. Equivalent to the inner product of mapped vectors in a high-dimensional space ,in This indicates transpose.

[0073] Furthermore, the method for calculating the squared Euclidean distance is as follows: For four-dimensional feature vectors and The squared Euclidean distance is calculated as follows: ,in , , , respectively, feature vectors The 1st, 2nd, 3rd, and 4th components, , , , respectively, feature vectors The 1st, 2nd, 3rd, and 4th components, For summation index.

[0074] Furthermore, the aforementioned bandwidth parameters Set to the reciprocal of the feature vector dimension, i.e. ,in This represents the dimension of the lubrication contamination feature vector.

[0075] The training process of the aforementioned single-class support vector machine is as follows: collect a sample set of lubrication contamination feature vectors of the main axis under healthy lubrication conditions. The decision boundary for health status is learned using a single-class support vector machine algorithm, where... For the first The lubrication contamination feature vector of each training sample. Let be the number of training samples. The optimization objective of a single-class support vector machine is... The constraints are , ,in For the summation index, Indicates transpose. For kernel function mapping, Parameters used to control the degree of boundary relaxation.

[0076] The decision function is obtained after training. The training adopts an unsupervised learning mode, using only healthy state samples, and the optimization strategy uses the sequential minimum optimization algorithm to solve the constrained optimization problem.

[0077] Furthermore, the above constraints apply to all training sample indices. All of these hold true, ensuring that each training sample point is either inside the decision boundary or determined by slack variables. It is allowed to lie outside the boundary, and the sum of the slack variables is penalized by the regularization term in the objective function.

[0078] Furthermore, the number of samples in the aforementioned training sample set The training samples are at least 100. The data is collected continuously for 8 hours under the condition that the spindle is using new lubricating oil and running stably. A lubrication contamination feature vector is generated every 60 seconds as a training sample.

[0079] Furthermore, the above parameters The range of values ​​is ,parameter Simultaneously, the upper bound of the proportion of support vectors in the training samples and the lower bound of the proportion of outlier samples outside the boundary are controlled in this embodiment. Set to 0.1, parameter Setting it to 0.1 allows up to 10% of healthy samples to lie outside the decision boundary to enhance the model's generalization ability.

[0080] It should be noted that the above calculation method for lubrication contamination risk score is as follows: The current lubrication contamination feature vector is... Input the decision function and calculate its value. The lubrication contamination risk score is defined as follows: .when The time indicates that the current feature vector is outside the health boundary, and the risk score is positive, with the larger the value, the further away from the healthy state.

[0081] The single-class support vector machine was trained on 120 feature vector samples collected from historical healthy lubricating oil conditions, with radial basis function kernel parameters set. Boundary relaxation parameters After training, a decision boundary containing 32 support vectors is obtained. The current lubrication contamination feature vector at time 10:15:00 is then used. (After standardization) Input decision function.

[0082] Calculate the radial basis function kernel values ​​between the current feature vector and each support vector, for example, with the first support vector. The square of the Euclidean distance is: Kernel function value The contributions of all support vectors are aggregated to obtain the decision function value. .because The current feature vector is outside the health boundary, lubrication contamination risk score .

[0083] Step 600: Determine the early warning threshold and generate targeted maintenance instructions; Determine whether the combined index of the lubrication contamination weight calculated in step 400 and the lubrication contamination risk score calculated in step 500 exceeds the lubrication maintenance threshold, or whether the heat source weight calculated in step 400 exceeds the heat compensation threshold. Generate a targeted early warning signal: output a lubricating oil dehydration treatment command or a replacement command, or output heat compensation parameters to the CNC system control module.

[0084] It should be noted that the calculation method for the above-mentioned joint index is as follows: the lubrication and maintenance joint index is defined as follows: ,in The lubrication contamination weight calculated in step 400, The lubrication contamination risk score calculated in step 500. If (in If the lubrication maintenance threshold is reached, a lubricating oil dehydration or replacement instruction will be generated. (in If the threshold value is set to the thermal compensation threshold, then a thermal compensation parameter adjustment command will be generated.

[0085] Furthermore, the method for determining the above-mentioned lubricating oil dehydration treatment and replacement is as follows: if the lubricating contamination risk score... Less than the preset replacement threshold If so, a lubricating oil dehydration command is output, and the water in the lubricating oil is removed by a vacuum dehydration device; if If so, a lubricant replacement command will be output, where the replacement threshold is... Set to the risk score value corresponding to the upper limit of dehydration processing capacity.

[0086] Furthermore, the aforementioned replacement threshold The method for determining the upper limit of the processing capacity of the vacuum dehydration device is as follows: The risk score change before and after dehydration is tested on a lubricating oil sample with a known water content. The risk score before dehydration corresponding to the critical water content at which the risk score cannot be reduced to a safe range after dehydration is used as the replacement threshold. .

[0087] It should be noted that the above thermal compensation parameters are calculated based on the thermal elongation coefficient obtained in step 300. and current spindle temperature Calculate thermal elongation ,in The reference temperature is used. The thermal compensation parameter output is the axial position correction amount. The thermal compensation parameters are transmitted to the CNC system control module to correct the Z-axis coordinate of the tool path.

[0088] Furthermore, the aforementioned reference temperature The initial temperature of the spindle during cold start is set. The reference temperature is obtained by averaging the spindle temperature over the first 5 minutes after the CNC machining center starts, using the following formula: ,in For the summation index, This represents the number of temperature sampling points within a 5-minute time window. The sampling time.

[0089] In this embodiment, to improve the accuracy of causal relationship identification, a stationarity test is performed on the time series of each variable before performing the Granger causality test in step 400. If a variable's time series is non-stationary, it is differencing until stationary before performing the Granger causality test. Specifically, the augmented Dickey-Fuller test (ADF test) is used to determine the stationarity of the time series. The input to the augmented Dickey-Fuller test is the time series, and the output is the ADF statistic and its p-value. If the p-value of the ADF statistic is greater than 0.05, the series is considered non-stationary, and a first-order difference is performed on the series. Then, perform the augmented Dickey-Fuller test again until the sequence is stationary.

[0090] In this embodiment of the application, in order to adapt to the dynamic changes of the threshold under different operating conditions, the lubrication maintenance threshold in step 600 is... and thermal compensation threshold An adaptive adjustment mechanism is adopted. Specifically, this involves using combined indicators based on historical data within the most recent time window. and heat source weight The quantiles are set as follows: the lubrication maintenance threshold is set to the 95th quantile of historical data, and the thermal compensation threshold is set to the 90th quantile of historical data. When changes in operating conditions cause data distribution shifts, the thresholds are adaptively adjusted accordingly, thereby avoiding false alarms or missed alarms under different operating conditions with fixed thresholds.

[0091] Furthermore, the above quantile calculation method is as follows: Arrange the historical data sequence within the time window in ascending order; the 95th quantile corresponds to the position index in the sorted sequence. The value, the index of the 90th percentile is The value of, among which The number of historical data samples. This is for rounding up.

[0092] Furthermore, the length of the historical data time window is set to 24 hours. The historical data time window can cover multiple complete processing cycles to obtain stable statistical characteristics, and can also maintain an appropriate response speed to changes in operating conditions.

[0093] Lubrication contamination weight calculated in step 400 And the lubrication contamination risk score calculated in step 500 Calculate the combined indicators of lubrication and maintenance: An adaptive threshold mechanism was used to statistically analyze the joint indicator sequence of historical data from the past 24 hours. The 95th percentile was 0.52, indicating a lubrication maintenance threshold. .because This triggers a lubrication maintenance warning. Further assessment of the lubrication contamination risk score is then conducted. With replacement threshold Due to the relationship The system outputs a lubricating oil dehydration treatment command. Table 6 summarizes the key parameters and output commands for early warning judgment.

[0094] Table 6: Early Warning Judgment Results and Maintenance Command Output: At 10:15:05, the monitoring system output a "lubricating oil dehydration treatment" command to the maintenance terminal, prompting the operator to start the vacuum dehydration device to dehydrate the spindle lubrication system, with an estimated treatment time of 1.5 hours. The system recorded the trigger time, key parameter values, and causal analysis results of this warning event for subsequent threshold adaptive adjustment and anomaly pattern learning. Because the causal analysis clearly identified the root cause of the anomaly as lubricant contamination rather than an external heat source, unnecessary lubricant replacement or incorrect adjustment of thermal compensation parameters due to misjudgment was avoided, reducing maintenance costs and ensuring machining accuracy.

[0095] Figure 5 The time-series monitoring of lubrication and maintenance joint indicators, weighting parameters, and thresholds is presented.

[0096] The multi-sensor fusion-based monitoring method for CNC machining centers provided in this embodiment achieves coordinated sensing of three physical domains—lubrication, thermal state, and mechanical state—by fusing data from three types of sensors: dielectric constant, temperature field, and position deviation. This overcomes the limitation that a single sensor cannot fully reflect the equipment status.

[0097] To address the issue that trace amounts of moisture have a negligible impact on the average dielectric constant, this embodiment employs high-pass filtering in step 200 to extract high-frequency pulse disturbance characteristics of the dielectric constant. A four-dimensional feature vector is constructed using pulse event frequency, average pulse amplitude, skewness parameter, and kurtosis parameter, overcoming the problem of average value calculation masking the pulse disturbances caused by trace moisture. When trace amounts of moisture are mixed in the lubricating oil, water droplets cause high-frequency pulse fluctuations in the dielectric constant during oil flow. Although the impact on the average value is less than the measurement accuracy, high-pass filtering can extract these high-frequency disturbances, and pulse statistical features can quantify the frequency and intensity of the disturbances, thus enabling detection at the trace moisture (water content <0.1%) stage.

[0098] To address the issue of the coupling effect of thermal expansion and mechanical runout on the three-dimensional position deviation of the tool tip, this embodiment calculates the thermal expansion coefficient in step 300 through axial-radial position decomposition, temperature correlation analysis, and linear regression, thus separating thermal expansion from mechanical runout and overcoming the inability to distinguish the physical mechanisms in the three-dimensional deviation coupling. Thermal expansion is mainly manifested in the axial component and is linearly correlated with temperature. By establishing a linear regression model of axial deviation and temperature and extracting the residuals, the temperature-related thermal expansion effect can be decoupled from the temperature-independent mechanical runout effect, providing a clear physical quantity for subsequent causal analysis.

[0099] To address the core issue that increased temperature and positional deviation may originate from lubrication contamination or external heat sources, but with different maintenance requirements, this implementation method constructs a four-variable causal topology diagram in step 400, comprising a lubrication contamination feature vector, temperature, axial deviation, and radial deviation. It utilizes the Granger causality test to quantify the strength of two causal paths: lubrication contamination → temperature → axial deviation and external heat source → temperature → axial deviation. This overcomes the limitation of single-threshold warnings in distinguishing the root causes of anomalies. The Granger causality test analyzes the lag correlation of time series to determine the causal direction and strength between variables, identifying whether lubrication contamination causes temperature increases leading to axial deviation, or whether the temperature increase is caused by an external heat source while lubrication contamination characteristics show no abnormalities. By calculating the causal strength of each path and determining the dominant path, the root cause of the anomaly can be located.

[0100] In step 500, a single-class support vector machine is used to establish an anomaly detection boundary for the healthy lubricating oil state. A risk score is generated by combining the distance from the lubricating contamination feature vector to the boundary, which quantifies the deviation of the current lubrication state from the healthy baseline. In step 600, the lubricating contamination weights and the risk score are combined to form a joint index. A lubrication maintenance command is triggered only when the causal path indicates that lubricating contamination is dominant and the risk score exceeds a threshold, avoiding misjudgments based solely on a single index. When an external heat source is identified as the dominant factor, thermal compensation parameters are output instead of a lubrication replacement command, achieving separate decision-making for lubrication maintenance and thermal compensation.

[0101] Therefore, this implementation overcomes the technical defects of average value calculation masking trace moisture disturbances, three-dimensional deviation coupling failing to distinguish mechanisms, and single threshold failing to identify root causes by using multi-sensor fusion, high-frequency pulse disturbance extraction, position-temperature decoupling, multi-physical quantity causal topology analysis, and separation early warning decision-making. This allows for the triggering of correct maintenance strategies before trace moisture causes serious lubrication failure, avoiding the increased costs and unresolved problems caused by misjudging lubrication contamination as a heat source or vice versa.

[0102] The embodiments of the present invention have been described above. However, the embodiments are not limited to the specific implementation methods described above. The specific implementation methods described above are merely illustrative and not restrictive. Those skilled in the art can make more equivalent embodiments under the guidance of the present embodiments, and all of them are within the protection scope of the present embodiments.

Claims

1. A monitoring method for CNC machining centers based on multi-sensor fusion, characterized in that, Includes the following steps: Acquire high-sampling-rate raw data from the dielectric constant sensor of the lubrication circuit of the CNC machining center, data from the distributed temperature sensor array of the spindle, and three-dimensional position deviation data of the tool tip. Synchronize the multi-source sensor data to a unified time reference based on the timestamp alignment algorithm. High-pass filtering is performed on the dielectric constant time series to generate a high-frequency pulse perturbation sequence of dielectric constant. The pulse event frequency and average pulse amplitude within the sliding window are calculated. The skewness and kurtosis parameters of the pulse amplitude distribution are calculated using a probability density estimation algorithm to generate a lubrication contamination feature vector. The three-dimensional position deviation of the tool tip is decomposed into axial and radial components. The average temperature sequence of the spindle is calculated. Correlation analysis is used to calculate the linear correlation coefficient between axial deviation and temperature. The thermal elongation coefficient is calculated based on the linear regression relationship between axial deviation and temperature. A four-variable causal topological directed graph is constructed, consisting of lubrication contamination feature vector, temperature, axial deviation, and radial deviation. The Granger causality test is used to calculate the causal relationship strength of the three paths from lubrication contamination feature vector to temperature, temperature to axial deviation, and lubrication contamination feature vector to axial deviation. Based on the causal relationship strength, the dominant causal path is identified and the lubrication contamination weight and heat source weight are calculated. An anomaly detection boundary is established on healthy lubricating oil data using a single-class support vector machine, and the distance from the current lubricating contamination feature vector to the boundary is calculated to generate a lubricating contamination risk score. Determine whether the product of the lubrication contamination weight and the lubrication contamination risk score exceeds the lubrication maintenance threshold, or whether the heat source weight exceeds the heat compensation threshold, and generate a targeted early warning signal. The Granger causality test is performed as follows: For two standardized time series variables X and Y, an autoregressive model including X is established. The historical values ​​of Y and X are used as independent variables to predict the current value of Y. The sum of squared residuals of the complete model including X and the constrained model not including X is calculated. The F-statistic is calculated based on the difference between the sums of squared residuals of the two models. The value of the F-statistic is used as the strength of the causal relationship between X and Y.

2. The monitoring method for CNC machining centers based on multi-sensor fusion according to claim 1, characterized in that, The high-pass filter is a Butterworth high-pass filter with a cutoff frequency of 100Hz and a filter order of 4. The method for identifying the pulse event is as follows: set an amplitude threshold, and identify a pulse event when the absolute value of the high-frequency disturbance sequence exceeds the amplitude threshold. Count the number of pulse events within the time window. The pulse event frequency is equal to the number of pulse events divided by the time window length, and the average pulse amplitude is equal to the arithmetic mean of the amplitudes of all pulse events within the window. The amplitude threshold is set to 3 times the standard deviation of the high-frequency perturbation sequence under healthy lubricating oil conditions, and the time window length is set to 60 seconds. The standard deviation under healthy lubricating oil conditions is calculated by continuously collecting high-frequency disturbance sequence data for 1 hour under stable operating conditions with new lubricating oil on the spindle, and using the standard deviation of the high-frequency disturbance sequence data as the baseline standard deviation of the healthy state.

3. The monitoring method for CNC machining centers based on multi-sensor fusion according to claim 1, characterized in that, The probability density estimation adopts the kernel density estimation method, with the input being a pulse amplitude sample sequence, and the probability density function is estimated based on the Gaussian kernel function. Based on the estimated probability density function, the ratio of the third central moment to the cube of the standard deviation is calculated as the skewness parameter, and the ratio of the fourth central moment to the fourth power of the standard deviation is calculated and then 3 is subtracted as the kurtosis parameter. The lubrication contamination feature vector is a four-dimensional vector, including pulse event frequency, average pulse amplitude, skewness parameter, and kurtosis parameter; Z-score standardization is performed on the lubrication contamination feature vector to convert each component into a standardized feature with zero mean and unit variance.

4. The monitoring method for CNC machining centers based on multi-sensor fusion according to claim 1, characterized in that, The decomposition method of the three-dimensional position deviation is as follows: Let the components of the three-dimensional position deviation of the tool tip measured by the laser rangefinder in the Cartesian coordinate system be the X-axis component, the Y-axis component, and the Z-axis component. The Z-axis is parallel to the main axis. The axial component is directly taken as the Z-axis component. The radial component is calculated as the square root of the sum of the squares of the X-axis component and the squares of the Y-axis component. The calculation method for the spindle average temperature sequence is as follows: the arithmetic average of the measurement values ​​of multiple temperature sensors deployed on the spindle at the same time is performed. The linear correlation coefficient between the axial deviation and temperature is the Pearson correlation coefficient, calculated by dividing the covariance of the axial deviation sequence and the average temperature sequence of the main shaft by the square root of the product of the variance of the axial deviation and the variance of the temperature.

5. The monitoring method for CNC machining centers based on multi-sensor fusion according to claim 1, characterized in that, The method for calculating the coefficient of thermal elongation is as follows: A linear regression model is established, in which the axial deviation is expressed as the thermal elongation coefficient multiplied by the average temperature of the main shaft, plus a constant term and a residual term. The thermal elongation coefficient and the constant term are estimated using the least squares method. The estimated value of the thermal elongation coefficient is equal to the covariance of the axial deviation and the temperature divided by the variance of the temperature. The temperature-independent axial runout residual sequence is equal to the measured axial deviation minus the estimated value of the thermal elongation coefficient multiplied by the average temperature of the spindle and then minus the estimated value of the constant term; The training samples were collected during the complete thermal equilibrium process of the spindle under normal operating conditions, with a data collection time of at least 2 hours to cover the temperature change range of the spindle from cold start to hot stability.

6. The monitoring method for CNC machining centers based on multi-sensor fusion according to claim 1, characterized in that, The Granger causality test is performed by standardizing the time series of each variable before execution to eliminate dimensional differences. The method for calculating the F-statistic is as follows: establish a constrained model without independent variable X and a complete model with X, calculate the sum of squared residuals of the two models respectively, and the F-statistic is equal to the difference between the sum of squared residuals of the constrained model and the complete model divided by the number of constraints, then divided by the sum of squared residuals of the complete model divided by the degrees of freedom. The lag order of the autoregressive model is determined by the Akaike Information Criterion. Different combinations of lag orders are traversed within the candidate range, and the combination of orders that minimizes the value of the Akaike Information Criterion is selected as the optimal lag order.

7. The monitoring method for CNC machining centers based on multi-sensor fusion according to claim 1, characterized in that, The rule for identifying the dominant causal path is as follows: Let the causal strength from the lubrication contamination feature vector to temperature be the first causal strength, the causal strength from temperature to axial deviation be the second causal strength, and the causal strength from the lubrication contamination feature vector to axial deviation be the third causal strength. If the first causal strength is greater than the causal significance threshold and the third causal strength is greater than the causal significance threshold, then it is determined to be the lubrication contamination dominant mode. The lubrication contamination weight is equal to the product of the first causal strength and the third causal strength divided by the product of the first causal strength and the third causal strength plus the second causal strength. If the first causal strength is less than the causal significance threshold and the second causal strength is greater than the causal significance threshold, it is determined to be an external heat source dominant mode. The heat source weight is equal to the second causal strength divided by the product of the first causal strength and the third causal strength plus the second causal strength. The causal significance threshold is set as the critical value of the F distribution at a significance level of 0.

05.

8. The monitoring method for CNC machining centers based on multi-sensor fusion according to claim 1, characterized in that, The training process of the single-class support vector machine is as follows: Collect a sample set of lubrication contamination feature vectors of the spindle under healthy lubrication oil conditions, and use a single-class support vector machine algorithm to learn the decision boundary of the healthy state. The optimization objective is to minimize the square of the norm of the weight vector divided by 2 minus the regularization term of the bias parameter plus the relaxation variable. The constraint is that the inner product of the feature vectors after the weight vector is mapped to the kernel is greater than or equal to the bias parameter minus the relaxation variable and the relaxation variable is non-negative. The single-class support vector machine uses a radial basis function kernel. The kernel function is equal to an exponential function multiplied by the square of the Euclidean distance between two feature vectors with a negative bandwidth parameter. The bandwidth parameter is set to the reciprocal of the feature vector dimension. The method for calculating the lubrication contamination risk score is as follows: input the current lubrication contamination feature vector into the decision function to calculate the decision function value. The lubrication contamination risk score is equal to the larger of 0 and the negative value of the decision function. When the decision function value is less than 0, it indicates that the current feature vector is outside the health boundary.

9. The monitoring method for CNC machining centers based on multi-sensor fusion according to claim 1, characterized in that, The method for generating the targeted early warning signal is as follows: The lubrication maintenance joint index equals the lubrication contamination weight multiplied by the lubrication contamination risk score; If the combined lubrication maintenance index is greater than the lubrication maintenance threshold, a lubricating oil dehydration treatment instruction or a replacement instruction will be generated. If the lubrication contamination risk score is less than the replacement threshold, a dehydration treatment instruction will be output. If the lubrication contamination risk score is greater than or equal to the replacement threshold, a lubricating oil replacement instruction will be output. If the heat source weight is greater than the heat compensation threshold, the heat expansion is calculated based on the heat expansion coefficient and the current spindle temperature. The heat expansion is equal to the heat expansion coefficient multiplied by the difference between the current spindle temperature and the reference temperature. The heat compensation parameter output is the negative value of the heat expansion as the axial position correction amount. The lubrication maintenance threshold and heat compensation threshold adopt an adaptive adjustment mechanism. Based on the quantile of the joint index of historical data and the weight of heat source within the most recent time window, the lubrication maintenance threshold is set to the 95th quantile of historical data, and the heat compensation threshold is set to the 90th quantile of historical data.

10. A fault prediction system for CNC machining centers based on multi-sensor fusion, used to execute the monitoring method for CNC machining centers based on multi-sensor fusion as described in any one of claims 1-9, characterized in that, include: The data acquisition module is used to acquire data from the dielectric constant sensor of the lubricating oil circuit, data from the spindle temperature sensor array, and data from the three-dimensional position deviation of the tool tip; The data synchronization module is used to synchronize multi-source sensor data to a unified time reference; The feature extraction module is used to extract high-frequency pulse perturbation features of dielectric constant and position-temperature correlation parameters; The causal analysis module is used to construct a causal topology graph of multiple physical quantities and identify the dominant causal path; The risk assessment module is used to calculate the lubrication contamination risk score; The early warning and decision-making module is used to generate targeted maintenance instructions.

Images