A rotary bearing multi-working condition fault diagnosis method and system for a filling machine

Abstract

The application discloses a kind of filling machine slewing bearing multi-working condition fault diagnosis method and system, belong to filling machine slewing bearing fault diagnosis technical field.The method includes: multi-source heterogeneous data acquisition and signal pre-processing;Based on the orthogonal decomposition of feature space, working condition division and benchmark construction;Based on the health model construction and threshold determination of working condition correlation characteristics;Real-time online monitoring and hierarchical fault diagnosis.The application is by collecting vibration signal and process parameter, after feature extraction and orthogonal decomposition dimension reduction, adopts density peak value search method to identify working condition cluster center and establishes each working condition benchmark characteristic vector;By screening working condition correlation characteristics and establishing non-fixed health boundary and working condition switching transition zone, the quantitative evaluation of equipment health state is realized by calculating comprehensive health index;Through real-time working condition identification, trend analysis and multi-level alarm mechanism, accurate fault diagnosis and predictive maintenance of filling machine slewing bearing under multiple working conditions are realized.

Description

Technical Field

[0001] This invention belongs to the field of mechanical equipment condition monitoring and fault diagnosis technology, and relates to a multi-condition fault diagnosis method and system for the slewing bearing of a filling machine. Background Technology

[0002] Filling machines are core equipment on production lines in industries such as food and beverage, pharmaceuticals, and chemicals. Their slewing bearings, as key components that support the filling turntable and transmit rotational power, directly affect the stability and production efficiency of the entire production line. In existing technologies, fault diagnosis of filling machine slewing bearings mainly employs a single vibration threshold method. This involves installing vibration sensors at the bearing location, setting a fixed vibration amplitude threshold, and determining a fault when the detected vibration value exceeds the threshold.

[0003] Existing technologies have explored the monitoring and diagnosis of slewing bearings. For example, CN102183951A discloses a LabVIEW-based slewing bearing monitoring device that uses an accelerometer to collect vibration data and achieves online diagnosis through signal denoising and time-frequency feature extraction. However, this solution uses only a single vibration sensor, lacks synchronous acquisition of process parameters, and has fixed diagnostic rules, making it difficult to adapt to complex production environments with multiple operating conditions. CN111476339B proposes a system that integrates multiple feature extraction methods and uses PCA dimensionality reduction followed by SVM for fault identification, emphasizing the completeness and sparsity of features. However, this technology still relies on a single feature set, does not fully utilize the correlation between process parameters and multi-source heterogeneous data, and lacks specialized processing for smooth transitions between different operating conditions. CN115688018B further focuses on multi-condition bearing monitoring, using a torque sensor to collect vibration data and achieving feature fusion and SVM classification through multiple steps such as empirical mode decomposition and kernel principal component analysis. Despite the improved multi-condition identification capability, it still relies on offline feature fusion and fixed threshold judgment, making it difficult to achieve real-time online hierarchical fault diagnosis and dynamic threshold transition between operating conditions.

[0004] The common shortcomings of the above technologies include: relying on only a single or a few sensors, failing to achieve unified acquisition and synchronization of multi-source heterogeneous data; feature extraction and dimensionality reduction often adopt fixed procedures, lacking adaptive threshold construction for different working conditions; and the lack of a transition zone for the switching process between adjacent working conditions, resulting in false alarms or missed alarms during the switching period. Health assessment remains at the level of a single feature threshold, lacking a quantitative representation of a comprehensive health index. Summary of the Invention

[0005] The technical problems to be solved by this invention include at least one of the following: how to achieve synchronous acquisition, unified preprocessing and time alignment of multi-source heterogeneous data of filling machine slewing bearings; how to use orthogonal decomposition of feature space for adaptive dimensionality reduction while ensuring the integrity of feature information, so as to adapt to the feature distribution of different working conditions; how to automatically divide working conditions according to the density distribution of features after dimensionality reduction and build a health model with non-fixed thresholds for each working condition; how to build a transition zone during the switching between adjacent working conditions to achieve smooth fault identification and avoid false alarms or missed alarms during the switching period; how to weight and fuse the deviation degree of each related feature into a single comprehensive health index to achieve hierarchical alarm and support subsequent remaining life prediction.

[0006] To address the aforementioned technical problems, the present invention provides the following technical solutions.

[0007] A multi-condition fault diagnosis method for the slewing bearing of a filling machine includes the following steps:

[0008] Step S1: Multi-source heterogeneous data acquisition and signal preprocessing: Vibration sensors are installed at the outer ring raceway of the slewing bearing of the filling machine and at the bearing housing of the drive motor to collect equipment vibration signals. At the same time, a data interface is configured in the filling machine control system to extract process parameters such as filling speed, cap torque, turntable angle position, filling pressure, and filling volume in real time. The vibration signals are processed by anti-aliasing filtering, resampling, timestamp alignment, and sliding window segmentation. The process parameters are validated and outliers are removed to obtain a preprocessed multi-source heterogeneous dataset.

[0009] Step S2: Working condition division and benchmark construction based on orthogonal decomposition of feature space: Time-domain features, frequency-domain features, and time-frequency-domain features are extracted from the preprocessed vibration signal to construct a high-dimensional feature vector. The orthogonal decomposition of feature space is used to reduce the dimensionality of the high-dimensional feature vector to obtain a low-dimensional feature vector. In the low-dimensional feature vector, the density peak search method is used to identify the center of the dense region of the historical healthy sample distribution as the working condition cluster center, and the sample space is divided into several non-overlapping working condition clusters. For each working condition cluster, its geometric center is calculated as the benchmark feature vector of that working condition, and the average distance from the samples in the cluster to the benchmark feature vector is calculated as the dispersion index.

[0010] Step S3: Health Model Construction and Threshold Determination Based on Working Condition Correlation Features: For each identified working condition cluster, calculate the Pearson correlation coefficient between each vibration feature and the process parameter to screen correlation features, and retain vibration features with a correlation degree higher than a set threshold to form a correlation feature vector; stack the correlation feature vectors of multiple historical health samples row by row to form a correlation feature matrix; based on the correlation feature matrix, establish an independent non-fixed health boundary for each working condition cluster, the non-fixed health boundary is determined by adding or subtracting a certain number of standard deviations from the mean of the correlation features, the multiple being determined according to the dispersion index of the working condition; establish a transition zone between the non-fixed health boundaries of adjacent working conditions, when a working condition switching event is detected, perform weighted fusion based on the health boundaries of the source working condition and the target working condition to form a health judgment standard during the transition period; fuse the health status of each correlation feature into a single comprehensive health index, the comprehensive health index being calculated based on the feature deviation of each correlation feature relative to the non-fixed health boundary;

[0011] Step S4: Real-time online monitoring and graded fault diagnosis: Perform the same preprocessing procedure as in step S1 on the real-time collected data, extract the real-time high-dimensional feature vector and project it onto the low-dimensional feature vector constructed in step S2, calculate the weighted Euclidean distance from the real-time sample to the benchmark feature vector of each working condition to identify the current working condition category; input the real-time associated feature values ​​into the health model of the current working condition, calculate the comprehensive health index and the feature deviation of each associated feature; maintain the trend sequence of the comprehensive health index within the sliding time window, determine whether the equipment health status is continuously deteriorating or fluctuating instantaneously through trend analysis, and predict the remaining life of the equipment based on the deterioration rate; establish a multi-level alarm mechanism based on the range of the comprehensive health index.

[0012] In one embodiment of the present invention, in step S1, two axes of the vibration sensor are arranged in mutually perpendicular directions in the horizontal plane, and the other axis is arranged in the vertical direction, so as to fully sense the radial and axial vibration responses of the slewing bearing.

[0013] In one embodiment of the present invention, in step S1, the window length of the sliding window segment is set as an integer multiple of the rotation period of the slewing bearing, and the window movement step size is set as the sliding window overlap rate.

[0014] In one embodiment of the present invention, in step S2, the time-domain features include root mean square value, peak value index, waveform index, impulse index, margin index, skewness index, and kurtosis index; the frequency-domain features include rotational frequency amplitude, rotational frequency harmonic amplitude, characteristic frequency amplitude of slewing bearing inner ring fault, characteristic frequency amplitude of slewing bearing outer ring fault, characteristic frequency amplitude of rolling element fault, characteristic frequency amplitude of cage fault, and energy proportion of preset frequency bands; the time-frequency domain features include the energy of each frequency band and its energy entropy value obtained by wavelet packet decomposition.

[0015] In one embodiment of the present invention, in step S2, the density peak search method includes: calculating the local density of each sample point, wherein the local density is defined as the number of other sample points within a certain distance range of the sample point; calculating the distance from each sample point to the nearest sample point with higher local density; constructing a decision graph with local density as the horizontal axis and the distance as the vertical axis, and identifying the density peak point located in the upper right corner region as the center of the working condition cluster.

[0016] In one embodiment of the present invention, in step S2, for a sample located in the boundary region of two working condition clusters, its membership degree to the center of the adjacent working condition cluster is calculated. The membership degree is calculated by weighting the inverse of the distance from the sample to each center and stored in the form of a membership degree vector for subsequent transition zone calculation.

[0017] In one embodiment of the present invention, in step S3, the threshold value of the Pearson correlation coefficient is in the range of 0.7 to 0.9; a multicollinearity test is performed on the features in the associated feature vector, and when the absolute value of the Pearson correlation coefficient between two features is greater than the multicollinearity judgment threshold, the feature with a higher degree of correlation with the process parameters is retained.

[0018] In one embodiment of the present invention, in step S3, the width of the transition band is determined based on the distance between the reference feature vectors of the two operating conditions and the switching rate. Within the transition band, the health boundary is a weighted fusion result of the source operating condition and the target operating condition boundary, and the weight is determined based on the progress ratio of the current moment in the transition process.

[0019] In one embodiment of the present invention, in step S4, when the real-time sample falls into the fuzzy boundary region between two working conditions, the health models of the two adjacent working conditions are activated simultaneously to enter the dual-working-condition monitoring mode, and the comprehensive health index under the two working conditions is calculated respectively. The comprehensive health index with the lower value is taken as the current judgment result.

[0020] This invention also provides a multi-condition fault diagnosis system for the slewing bearing of a filling machine, comprising:

[0021] The data acquisition unit is used to collect equipment vibration signals at the outer ring raceway of the slewing bearing of the filling machine and the bearing housing of the drive motor, and to extract process parameters such as filling speed, cap torque, turntable angle position, and filling pressure from the filling machine control system in real time.

[0022] The signal preprocessing unit receives the vibration signal and process parameters from the data acquisition unit, performs anti-aliasing filtering, resampling, timestamp alignment and sliding window segmentation on the vibration signal, performs validity verification and outlier removal on the process parameters, and outputs the preprocessed multi-source heterogeneous dataset.

[0023] The feature extraction and dimensionality reduction unit receives the preprocessed multi-source heterogeneous dataset from the signal preprocessing unit, extracts time-domain features, frequency-domain features, and time-frequency-domain features from the vibration signal to construct a high-dimensional feature vector, and uses the feature space orthogonal decomposition method to perform dimensionality reduction processing on the high-dimensional feature vector to reduce the dimensionality and output a low-dimensional feature vector.

[0024] The working condition division unit receives the low-dimensional feature vector output by the feature extraction and dimensionality reduction unit, and uses the density peak search method to identify the center of the dense region of the historical health sample distribution as the working condition cluster center, and divides the sample space into several non-overlapping working condition clusters.

[0025] The benchmark construction unit receives the working condition cluster division result from the working condition division unit, calculates the geometric center of each working condition cluster as the benchmark feature vector of that working condition, and calculates the average distance from the samples within the cluster to the benchmark feature vector as the dispersion index.

[0026] The health model construction unit receives the baseline feature vector and dispersion index from the baseline construction unit, calculates the Pearson correlation coefficient between each vibration feature and the process parameters to screen related features, retains vibration features with a correlation degree higher than a set threshold to form a related feature vector; stacks the related feature vectors of multiple historical health samples row by row to form a related feature matrix; based on the related feature matrix, establishes an independent non-fixed health boundary for each working condition cluster, establishes a transition zone between the non-fixed health boundaries of adjacent working conditions, and outputs health model parameters.

[0027] The comprehensive health index calculation unit receives the health model parameters from the health model construction unit and integrates the health status of each related feature into a single comprehensive health index.

[0028] The real-time monitoring unit is connected to the feature extraction and dimensionality reduction unit, the benchmark construction unit, the health model construction unit, and the comprehensive health index calculation unit, respectively. It is used to project the real-time collected data onto the low-dimensional feature vector after preprocessing and feature extraction, calculate the weighted Euclidean distance from the real-time sample to the benchmark feature vector of each working condition to identify the current working condition category, and simultaneously activate the health models of two adjacent working conditions in the fuzzy area of ​​the working condition boundary, and take the comprehensive health index with the lower value as the current judgment result.

[0029] The trend analysis and prediction unit is used to receive the comprehensive health index from the comprehensive health index calculation unit, maintain the trend sequence of the comprehensive health index within the sliding time window, determine whether the health status of the equipment is continuously deteriorating or fluctuating instantaneously through trend analysis, and predict the remaining lifespan of the equipment based on the deterioration rate.

[0030] The alarm unit is connected to the comprehensive health index calculation unit and the trend analysis and prediction unit, respectively, and is used to establish a multi-level alarm mechanism based on the range of the comprehensive health index.

[0031] Compared with the prior art, the present invention has the following beneficial effects:

[0032] This invention achieves high-precision time synchronization of vibration signals and process parameters through unified multi-source data acquisition and anti-aliasing preprocessing; based on orthogonal decomposition of feature space, it retains only a few orthogonal features with a cumulative contribution rate of not less than 85%, effectively compressing the feature dimension while maintaining key operating condition information; it uses a density peak search method to naturally divide operating conditions and construct corresponding benchmark feature vectors, enabling the diagnostic model to adapt to different operating states; it uses weighted fusion to form a transition zone for thresholds of adjacent operating conditions, significantly reducing the false alarm rate during operating condition switching; and it obtains a unified comprehensive health index and implements hierarchical alarms by weighted summarization of the deviation of all related features, improving the sensitivity of early fault identification and the reliability of diagnosis. Attached Figure Description

[0033] Figure 1 This is a flowchart illustrating the overall process of the fault diagnosis method of the present invention.

[0034] Figure 2 The figure shows the experimental results of the sliding window parameter optimization of the present invention; wherein: Figure 2 A represents the experimental results of the sliding window overlap rate; Figure 2 B represents the window length test result;

[0035] Figure 3 This is a peak density decision distribution diagram for the present invention.

[0036] Figure 4 This is a block diagram of the modular architecture of the fault diagnosis system of the present invention. Detailed Implementation

[0037] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the invention.

[0038] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments.

[0039] like Figure 1 As shown, the present invention provides a multi-condition fault diagnosis method for the slewing bearing of a filling machine, comprising the following steps:

[0040] Step S1: Multi-source heterogeneous data acquisition and signal preprocessing;

[0041] Step S2: Working condition division and benchmark construction based on orthogonal decomposition of feature space;

[0042] Step S3: Construction of a health model based on working condition correlation features and determination of thresholds;

[0043] Step S4: Real-time online monitoring and graded fault diagnosis.

[0044] Step S1 further includes the following steps:

[0045] Step S1-1: Arrange triaxial MEMS accelerometers at the outer raceway of the slewing bearing of the filling machine and at the bearing housing of the drive motor. The triaxial MEMS accelerometers are also used as vibration sensors to collect vibration signals from the equipment.

[0046] Preferably, the sampling frequency of the triaxial MEMS accelerometer is set based on the highest passing frequency of the slewing bearing to ensure that the bearing fault characteristic frequency and its harmonic components can be captured. Simultaneously, a data interface is configured in the filling machine control system to extract process parameters such as filling speed, cap torque, turntable angle position, and filling pressure in real time.

[0047] Preferably, the triaxial MEMS accelerometer is fixed to the metal surface by magnetic attraction or adhesive bonding and connected to the edge computing gateway via a shielded cable.

[0048] Furthermore, two axes of the triaxial MEMS accelerometer are arranged in mutually perpendicular directions in the horizontal plane, and the other axis is arranged in the vertical direction, so as to comprehensively sense the radial and axial vibration response of the slewing bearing.

[0049] Steps S1-2: Anti-aliasing filtering is performed on the raw vibration signal output from the triaxial MEMS accelerometer to remove frequency components higher than half the sampling frequency, preventing spurious frequency folding during spectrum analysis. Resampling is then performed to unify the data rate, ensuring data from different sensors and time periods have the same temporal resolution. Time stamp alignment of process parameters is then performed to establish a strict correspondence between vibration data and process data in the time dimension, with an alignment accuracy no less than that of a single sampling period.

[0050] Furthermore, a sliding window is used to segment the vibration signal. The window length is set as an integer multiple of the slewing bearing rotation period to ensure that each data segment contains complete information about the rotational machinery cycle and eliminate spectral leakage caused by truncation. The window movement step size is set based on the sliding window overlap rate, preferably 50%, with the corresponding window movement step size being 50% of the window length.

[0051] Schematic, the sliding window overlap rate is calculated as O=(LS) / L×100%, where: O is the sliding window overlap rate, ranging from 30% to 70%, with a typical value of 50%; L is the sliding window length, in units of sampling points, and is 3 to 5 times the number of sampling points corresponding to the rotation cycle of the slewing bearing; S is the window movement step size, in units of sampling points.

[0052] like Figure 2 As shown, through comparative experiments with multiple sets of different overlap rates, when the overlap rate is less than 30%, the continuity of features between adjacent windows is broken, and the early fault missed rate increases significantly; when the overlap rate is greater than 70%, the computational load increases significantly but the improvement in diagnostic accuracy is limited.

[0053] The results of the two sets of control experiments show that... Figure 2 In the sliding window overlap rate test, when the overlap rate is less than 30%, the early fault missed rate increases significantly and the diagnostic accuracy decreases significantly. This is because the continuity of fault features between adjacent analysis windows is broken. When the overlap rate exceeds 70%, the diagnostic accuracy only slightly improves, but the relative computational load of the algorithm increases sharply, and the consumption of computational resources increases significantly. The overlap range of 30% to 70% is a reasonable range for parameter selection. Among them, a 50% overlap rate can keep the computational overhead at a reasonable level while ensuring the continuity of fault features and reducing the early fault missed rate. Figure 2The window length test (B) uses the slewing bearing rotation cycle as a benchmark. Fault characteristics can be stably displayed within 3 to 5 rotation cycles. When the window length is less than 3 times the rotation cycle, the fault characteristic capture rate is low, and the fault impact information cannot be completely recorded. When the window length exceeds 3 times the rotation cycle, the fault location error continues to increase, the location accuracy gradually deteriorates, and the overall performance index gradually decreases after the 3-cycle position. Therefore, selecting 3 times the rotation cycle as the window length can not only completely capture the fault impact signal, but also avoid the problem of decreased fault location accuracy caused by excessively long windows. Based on the two sets of test data, the present invention preferably uses a parameter combination of 50% window overlap rate and 3 times the rotation cycle window length.

[0054] Steps S1-3: Outlier removal is performed on the aligned multi-source data. Specifically, the statistics of vibration signals within each data segment are calculated. If the peak value or energy value of a certain data segment exceeds the historical normal range by several times, the segment is identified as an outlier segment affected by a sudden impact and is removed. For process parameters, validity is verified by setting a physically reasonable range, and erroneous values ​​caused by sensor failure or communication interruption are removed. After outlier removal, a clean multi-source heterogeneous dataset is obtained. Further, the continuity of the data segments after outlier removal is checked. If the number of consecutively missing data segments exceeds a set threshold, the time period is marked as a data missing interval and is not included in subsequent modeling.

[0055] Steps S1-4: The preprocessed vibration data and process data are locally cached and initially compressed at the edge computing gateway, and then uploaded to the cloud platform via wired or wireless network.

[0056] The edge computing gateway performs preliminary feature extraction on the raw vibration signal, uploading only the feature values ​​and some key waveform data to reduce network bandwidth consumption. Correspondingly, the cloud platform receives and stores data from multiple filling machines, forming a distributed device database. The cloud platform performs secondary verification on the uploaded data, checking data integrity, timestamp continuity, and feature value rationality. Data that passes verification is added to the historical sample database, while data that fails verification is returned to the edge computing gateway for re-collection.

[0057] Step S2 further includes the following steps:

[0058] Step S2-1: Extract multi-dimensional feature values ​​from the preprocessed vibration signal to construct a high-dimensional feature vector.

[0059] In the time domain, the root mean square value, peak value, waveform value, impulse value, margin value, skewness value, and kurtosis value of each data segment are calculated. These statistics reflect the overall energy level, impact characteristics, and probability distribution of the vibration signal.

[0060] In the frequency domain, the time-domain signal is converted to the frequency domain by fast Fourier transform, and the amplitude at the frequency shift and its harmonics, the characteristic frequency amplitude of the inner ring fault of the slewing bearing, the characteristic frequency amplitude of the outer ring fault of the slewing bearing, the characteristic frequency amplitude of the rolling element fault, and the energy proportion of the preset frequency band are extracted.

[0061] In the time-frequency domain, wavelet packet decomposition is used to decompose the signal into multiple frequency bands, and the energy and energy entropy of each band are calculated to reflect the frequency structure changes of the signal at different resolutions. The aforementioned time-domain, frequency-domain, and time-frequency-domain feature values ​​are arranged in order to form a high-dimensional feature vector describing the current equipment state. Simultaneously, process parameters such as filling speed, cap torque, turntable angle position, filling pressure, and filling volume are normalized and used as normalized process parameters for subsequent correlation analysis.

[0062] The feature extraction process is completed at two levels: the edge computing gateway and the cloud platform. The edge computing gateway is responsible for extracting simple time-domain features in real time, while the cloud platform is responsible for extracting complex frequency-domain and time-frequency-domain features.

[0063] Step S2-2: Due to information redundancy, dimensional differences, and inter-feature correlations in high-dimensional feature vectors, directly using them for work condition classification would result in high computational complexity and the classification results would be affected by redundant information. This invention employs an orthogonal decomposition method of feature space to reduce the dimensionality of high-dimensional features.

[0064] Specifically, a large number of samples from historical healthy operation phases are collected, with each sample corresponding to a high-dimensional feature vector. If the number of samples is n and the feature dimension of each sample is m, then an n-row m-column data matrix is ​​formed.

[0065] First, the data matrix is ​​centered by subtracting the sample mean of each feature dimension from the value of that dimension, so that the mean of each dimension is 0.

[0066] Subsequently, the covariance structure of the centered data matrix is ​​calculated, constructing an m x m covariance matrix. The characteristic equation of this covariance matrix is ​​solved to obtain m eigenvalues ​​and their corresponding m orthogonal eigenvectors. These eigenvalues ​​are then sorted in descending order of magnitude, and the first k orthogonal eigenvectors are selected to form a low-dimensional projection space. The selection of k is based on the cumulative eigenvalue contribution rate.

[0067] In a preferred embodiment, the cumulative eigenvalue contribution rate In the formula: η k The cumulative eigenvalue contribution rate of the first k eigenvalues, expressed as a percentage, with a value not less than 85%; λ i Let λ be the i-th eigenvalue, arranged in descending order, i.e., λ1≥λ2≥…≥λ m ≥0. k is the number of orthogonal feature vectors selected, ranging from 3 to 8; m is the dimension of the original high-dimensional feature vector, with a value of 24.

[0068] Through comparative analysis of multiple sets of different thresholds, it was found that when the cumulative feature value contribution rate is below 85%, crucial operating condition variation information is lost, leading to a significant decrease in the accuracy of operating condition segmentation. When the cumulative feature value contribution rate is above 95%, the dimension of the low-dimensional feature vector usually needs to be increased to more than 10, increasing the computational load but providing only limited improvement in segmentation accuracy. Therefore, in this invention, the 85% threshold can reduce computational complexity while retaining most of the effective operating condition information.

[0069] The 24-dimensional features extracted by this invention include 7 time-domain features, 10 frequency-domain features, and 7 time-frequency-domain features, covering the typical characteristic patterns of slewing bearing faults.

[0070] Schematic, the 24-dimensional feature X extracted by this invention is X=[T1,T2,T3,T4,T5,T6,T7,F1,F2,F3,F4,F5,F6,F7,F8,F9,F 10 [W1, W2, W3, W4, W5, W6, W7]. The seven time-domain features are: T1 root mean square value, T2 peak value, T3 waveform value, T4 pulse value, T5 margin value, T6 skewness value, and T7 kurtosis value; the ten frequency-domain features are: F1 slewing bearing rotational frequency amplitude, F2 second harmonic amplitude, F3 third harmonic amplitude, F4 slewing bearing inner ring fault characteristic frequency amplitude, F5 slewing bearing outer ring fault characteristic frequency amplitude, F6 rolling element fault characteristic frequency amplitude, F7 cage fault characteristic frequency amplitude, F8 low-frequency energy proportion, F9 mid-frequency energy proportion, F... 10 The proportion of high-frequency band energy; the seven time-frequency domain features are: W1 first frequency band energy, W2 second frequency band energy, W3 third frequency band energy, W4 fourth frequency band energy, W5 fifth frequency band energy, W6 sixth frequency band energy, and W7 wavelet packet energy entropy.

[0071] The original high-dimensional feature vectors are projected one by one onto the low-dimensional feature vectors spanned by the selected k orthogonal feature vectors. The projection operation is achieved by calculating the inner product of the original vector and each orthogonal feature vector, resulting in the k-dimensional feature representation after dimensionality reduction.

[0072] The following example, using high-speed, large-capacity filling operations, illustrates the specific implementation process of orthogonal decomposition of the feature space.

[0073] According to the requirements of step S2-2, the first three orthogonal eigenvectors are selected to form a low-dimensional projection space. In practical applications, the cumulative eigenvalue contribution rate should not be less than 85%, and in typical cases, the cumulative contribution rate of the first three orthogonal eigenvectors can reach 85% to 92%.

[0074] In an illustrative application scenario, operational data is collected from the slewing bearing of a 36,000 bottles / hour filling machine on a food and beverage production line. The equipment is in a stable operating state. The characteristic values ​​for each dimension are illustrated below (actual values ​​may vary depending on the equipment model and operating conditions):

[0075] X=[0.82,3.15,1.24,4.72,5.89,0.12,2.98,0.35,0.18,0.09,0.07,0.05 ,0.04,0.03,0.21,0.52,0.27,0.42,0.31,0.15,0.08,0.03,0.01,1.26];

[0076] The time-domain characteristic T1 (root mean square value) is 0.82g, which is consistent with the normal vibration level of the slewing bearing under high-speed filling conditions.

[0077] The time-domain characteristic T7 (kurtosis index) is 2.98, which is close to the kurtosis value of 3 of a normal distribution, indicating that the equipment has no obvious impact failure.

[0078] The frequency domain characteristic F1 (rotational frequency amplitude) is 0.35, corresponding to the rotational frequency characteristic of the slewing bearing at a rotational speed of 15 r / min;

[0079] Frequency domain characteristics F8~F 10 The values ​​are 0.21, 0.52, and 0.27 respectively, and the sum of the three is 1.0, indicating that the vibration energy is mainly concentrated in the mid-frequency range, which is consistent with the frequency distribution characteristics of a healthy bearing.

[0080] The time-frequency domain feature W7 (wavelet packet energy entropy) is 1.26, indicating that the frequency structure of the vibration signal is relatively uniform and there are no obvious abnormal frequency components.

[0081] By performing orthogonal decomposition of the feature space on historical health samples, three orthogonal feature vectors, V1, V2, and V3, are obtained. Each vector contains 24 elements, corresponding to the weight coefficients of the original 24-dimensional features. The three vectors are pairwise orthogonal, meaning that the theoretical value of the inner product of any two vectors is 0.

[0082] Schematic, the orthogonal eigenvector V1 (first principal direction, with a typical contribution rate of about 50% to 55%) mainly reflects the overall changes in the filling machine speed and load, and is most sensitive to process parameters such as filling speed and filling volume.

[0083] V1=[0.32,0.28,0.15,0.22,0.19,0.08,0.12,0.41,0.35,0.27,0.18,0.1 2,0.09,0.06,0.38,0.25,0.17,0.21,0.16,0.11,0.07,0.04,0.02,0.14].

[0084] The orthogonal eigenvector V2 (second principal direction, with a typical contribution rate of about 20% to 25%) mainly reflects the impact characteristics of the vibration signal and is most sensitive to early pitting and spalling faults in bearings.

[0085] V2=[0.12,0.35,0.09,0.42,0.38,0.15,0.45,0.08,0.12,0.18,0.32,0.2 8,0.24,0.19,0.11,0.15,0.22,0.17,0.21,0.26,0.31,0.27,0.18,0.33].

[0086] The orthogonal eigenvector V3 (the third principal direction, with a typical contribution rate of about 10% to 15%) mainly reflects the frequency structure changes of the vibration signal and is most sensitive to faults such as poor lubrication and improper installation.

[0087] V3=[0.18,0.12,0.35,0.08,0.11,0.42,0.09,0.15,0.21,0.28,0.12,0.1 7,0.22,0.27,0.25,0.38,0.32,0.31,0.26,0.21,0.16,0.11,0.07,0.24].

[0088] Orthogonality verification:

[0089] V1·V2=0.32×0.12+0.28×0.35+...+0.14×0.33=0.0002≈0.

[0090] V1·V3=0.32×0.18+0.28×0.12+...+0.14×0.24=0.0001≈0.

[0091] V2·V3=0.12×0.18+0.35×0.12+...+0.33×0.24=0.0003≈0.

[0092] Theoretically, V1·V2=0, V1·V3=0, and V2·V3=0. However, in actual calculations, due to the precision limitations of floating-point operations and the randomness of sample data, the inner product result is usually a very small value close to 0 (e.g., on the order of 10). -3 (The following is an example of an error that is acceptable within the scope of engineering applications.)

[0093] The projection operation is achieved by calculating the inner product of the original high-dimensional eigenvector X and each orthogonal eigenvector. The inner product operation is the sum of the element-wise multiplication of two vectors, and the formula is as follows: , where Y i Let X be the i-th eigenvalue after dimensionality reduction. jV is the j-th element of the original feature vector. i,j It is the j-th element of the i-th orthogonal eigenvector.

[0094] For the above illustrative data, the calculations are as follows: Y1=X·V1≈4.70; Y2=X·V2≈7.80; Y3=X·V3≈3.35.

[0095] The 3D feature vector after dimensionality reduction is Y = [Y1, Y2, Y3] ≈ [4.70, 7.80, 3.35].

[0096] The original 24-dimensional features were compressed into 3-dimensional features, reducing the data volume by approximately 87.5%. If the cumulative feature value contribution rate is 89.2%, then approximately 89% of the effective information is retained. Y1≈4.70 corresponds to the typical load level under high-speed, high-capacity filling conditions; Y2≈7.80 indicates that the equipment has no obvious impact failures; Y3≈3.35 indicates that the equipment's frequency structure is normal and the lubrication condition is good.

[0097] By calculating the weighted Euclidean distance between the 3D feature vector and the baseline feature vector of each working condition, the current working condition category of the equipment can be identified.

[0098] Preferably, directions with eigenvalues ​​close to 0 or negative indicate that the data in that direction has almost no variation and are therefore discarded; only directions with positive and significant eigenvalues ​​are retained. The selection of orthogonal eigenvectors also needs to consider physical interpretability. If an orthogonal direction corresponds to a clear physical meaning, such as primarily reflecting changes in rotational speed or load fluctuations, then that direction is preferentially retained.

[0099] Steps S2-3: In the reduced-dimensional feature vector, historical healthy samples naturally cluster into several dense regions based on their distribution density. This invention uses a density peak search method to identify the centers of these dense regions. Specifically, the local density of each sample point is calculated, and the local density is defined as the number of other sample points within a certain distance of that sample point.

[0100] Preferably, the local density calculation method for sample points is as follows: In the formula: ρ s Let d be the local density of the s-th sample point, taking the value of a non-negative integer; st d is the Euclidean distance between the s-th sample point and the t-th sample point; c The cutoff distance is the 1% to 2% quantile of the distance between all sample points; χ(x) is the indicator function, which takes the value 1 when x < 0 and 0 when x ≥ 0; n' is the total number of historical healthy samples.

[0101] When the cutoff distance is set to the 1% quantile, the local density becomes overly sensitive to changes in sample distribution, easily generating too many false cluster centers. When the cutoff distance is set to the 3% quantile, the local density becomes too smooth, making it impossible to distinguish adjacent dense regions. The 1% to 2% quantile range can effectively avoid false clustering and missed identification while ensuring the accuracy of cluster center identification.

[0102] Simultaneously, the distance from each sample point to the nearest sample point with higher local density is calculated. If the sample point itself has the highest local density, its distance is taken as the distance to the farthest sample point.

[0103] A decision map is constructed with local density as the horizontal axis and distance as the vertical axis. Sample points located in the upper right corner of the decision map are identified. These points have a local density that is significantly higher than that of the surrounding samples and are far away from other high-density points. These points are density peak points and serve as the centers of the initial working condition cluster.

[0104] Centered on each density peak point, the sample space is divided into several non-overlapping operating condition clusters based on the distance between samples and local density changes. During the division process, each sample is assigned to the operating condition cluster represented by the nearest density peak point. However, if the distance of a sample to the center of the cluster exceeds several times the average dispersion index of the cluster, the sample is considered an outlier and is temporarily excluded from any operating condition cluster. Outliers are marked separately for manual verification.

[0105] Each operating condition cluster corresponds to a standard operating state of the filling machine, such as high-speed filling, low-speed adjustment, start-stop transition, cleaning and maintenance, and standby. Furthermore, for samples located at the boundary between two operating condition clusters, their membership degree to the center of the adjacent cluster is calculated. The membership degree is calculated by weighting the inverse of the distance from the sample to each center; the closer the distance, the higher the membership degree. Soft partitioning is performed based on the membership degree magnitude to ensure the continuity and smoothness of operating condition discrimination and avoid hard misclassification of samples at the boundary. The soft partitioning results are stored in the form of a membership degree vector for subsequent transition zone calculations.

[0106] like Figure 3 As shown, in one application scenario, the decision graph is constructed based on 1200 historical healthy operation samples of the slewing bearing of a filling machine. The horizontal axis represents the local density ρ of the sample points. s The vertical axis represents the relative distance δ between sample points. s All sample points are calculated based on their ρ. s and δ s The values, when mapped to a two-dimensional plane, exhibit the following distribution characteristics:

[0107] The upper right corner area contains 5 isolated high-density points, which are density peak points, corresponding to the 5 standard operating conditions of the filling machine;

[0108] Central region: contains a large number of ordinary sample points, forming dense clusters around their respective density peak points;

[0109] The lower left corner area contains a small number of discrete points, i.e., outliers, which correspond to instantaneous anomalies or data acquisition errors during equipment operation.

[0110] Schematic, five density peak points can be identified in the decision map. The parameters and corresponding operating conditions of each point are as follows:

[0111] Peak point A (ρ) s ≈110, δ s ≈5.2): Corresponding to high-speed, large-capacity filling conditions, it has the largest number of samples, the highest local density, and is far from other conditions;

[0112] Peak point B (ρ) s ≈85, δ s ≈4.7): This corresponds to medium-speed, medium-capacity filling conditions, with the second largest sample size.

[0113] Peak point C (ρ) s ≈60, δ s ≈4.1): This corresponds to low-speed, small-capacity filling conditions, with a relatively small sample size;

[0114] Peak point D (ρ) s ≈40, δ s ≈3.8): This corresponds to a relatively small sample size for cleaning and maintenance conditions.

[0115] Peak point E (ρ) s ≈30, δ s ≈5.5): Corresponding to the standby condition, the number of samples is the smallest, the local density is the smallest, but the relative distance is the largest.

[0116] The five peak points mentioned above all meet the characteristics of having a significantly higher local density than the surrounding samples and being far away from other high-density points, and therefore were selected as the initial working condition cluster centers.

[0117] Centered on each density peak point, each sample point is assigned to the operating condition cluster represented by the nearest density peak point, forming 5 non-overlapping operating condition clusters:

[0118] Cluster A: Corresponds to high-speed, large-capacity filling conditions, with a sample size accounting for approximately 9% to 10% of the total healthy samples;

[0119] Cluster B: Corresponds to medium-speed, medium-volume filling conditions, accounting for approximately 25% to 35% of the sample size;

[0120] Cluster C: Corresponds to low-speed, small-capacity filling conditions, accounting for approximately 15% to 25% of the sample size;

[0121] Cluster D: Corresponds to cleaning and maintenance conditions, accounting for approximately 5% to 10% of the sample size;

[0122] Cluster E: Corresponds to standby mode, with a sample size of approximately 3% to 8%.

[0123] Samples whose distance from the cluster center exceeds 2.5 times the average dispersion index of the cluster are marked as outliers. Outliers typically account for 0.5% to 2% of the total number of healthy samples, are not included in any condition cluster, and are marked separately for manual verification.

[0124] For a sample located at the boundary between two operating condition clusters, its membership degree to the center of the adjacent operating condition cluster is calculated. The membership degree is calculated by weighting the inverse of the distance from the sample to each center. For example, if a sample is 0.8 from the center of cluster A and 1.2 from the center of cluster B, its membership degree to cluster A is (1 / 0.8) / [(1 / 0.8)+1 / 1.2]=0.6, and its membership degree to cluster B is 0.4. The soft partitioning result of this sample is stored in the form of a membership degree vector [0.6,0.4,0,0,0] for subsequent transition zone calculations, ensuring the continuity and smoothness of operating condition discrimination.

[0125] Step S2-4: For each established cluster of operating conditions, calculate its geometric center as the reference feature vector for that operating condition. This reference feature vector represents the typical characteristic pattern of the equipment health status under that operating condition.

[0126] When calculating the geometric center, the arithmetic mean of the low-dimensional feature vectors of all samples within the cluster is taken to obtain the mean vector. Each component of this mean vector represents the average projection value of the operating condition on the corresponding orthogonal feature direction. Simultaneously, the average distance from all samples within the cluster to the geometric center is calculated as a dispersion index, reflecting the range of normal equipment fluctuations under that operating condition. Further, the variance contribution of each operating condition cluster in each feature direction in the original high-dimensional feature space is calculated to form a feature weight vector. The variance contribution is obtained by calculating the proportion of the variance of each feature within that operating condition cluster to the total variance of all features. This baseline feature vector, the dispersion index, and the feature weight vector together constitute the basic parameters of the operating condition discrimination rule base. When a new real-time sample enters the system, it is projected onto the low-dimensional feature vector, and its weighted Euclidean distance to the baseline feature vectors of each operating condition is calculated. The weighted Euclidean distance is determined based on the weights of the low-dimensional feature values ​​(i.e., the normalized value of the contribution rate of each principal direction), ensuring that feature directions with high information content have a greater weight in the distance calculation.

[0127] Continue Figure 3The collected data will be used to further illustrate this. After identifying five initial operating condition clusters using the density peak search method and completing cluster assignment for all samples, the geometric center of each independent operating condition cluster is first calculated in the 3D low-dimensional orthogonal space constructed in step S2-2. This geometric center is the reference feature vector for the corresponding operating condition. The calculation process involves taking the arithmetic mean of the low-dimensional feature vectors of all samples within the cluster dimension by dimension. The resulting mean vector is the reference feature vector for that operating condition, with each component corresponding to the average projection value of the equipment on each orthogonal feature direction under that operating condition.

[0128] Taking the high-speed, high-capacity filling operation cluster as an example, the cluster contains 112 healthy operating samples (accounting for approximately 9.3% of the total number of healthy samples). Each sample is a 3-dimensional low-dimensional feature vector. The first component of the reference vector is obtained by summing the projection values ​​of the first dimension of all samples and dividing by 112, which corresponds to the typical load level under the operation condition. Similarly, the second and third components are calculated, which correspond to the typical impact characteristics and frequency structure characteristics, respectively.

[0129] The reason for using the arithmetic mean to calculate the geometric center is that each dimension is independent in the orthogonal feature space. The arithmetic mean is the optimal cluster center estimate in the least squares sense, which can represent the typical characteristic pattern of the equipment health status under the working condition to the greatest extent and eliminate the influence of random fluctuations of individual samples.

[0130] After obtaining the baseline feature vectors for each working condition, the dispersion index of each working condition cluster is further calculated.

[0131] The magnitude of the dispersion index is directly related to the stability of the operating conditions. Generally, high-speed, high-capacity filling operations have relatively large load fluctuations during production, resulting in a larger dispersion index (e.g., 0.20 to 0.30); while in standby conditions, the equipment experiences almost no load fluctuations, leading to a smaller dispersion index (e.g., 0.10 to 0.20). Specific values ​​vary depending on the equipment and production environment and must be calculated based on actual historical data.

[0132] The dispersion index serves two purposes: first, it acts as a boundary threshold for real-time operating condition identification. When the weighted Euclidean distance from a real-time sample to the benchmark center of a certain operating condition is less than 1.5 times the dispersion of that operating condition, the sample is determined to clearly belong to that operating condition; second, it provides a basis for the dynamic adjustment of the health boundary in the subsequent step S3. Operating conditions with larger dispersion will adopt a wider health boundary.

[0133] Calculate the high-dimensional feature weight vector for each operating condition cluster. This weight vector is calculated in the original 24-dimensional high-dimensional feature space by calculating the proportion of the variance of each original feature within that operating condition cluster to the total variance of all features.

[0134] A larger variance contribution indicates a higher degree of variation for that feature under that operating condition, encompassing richer information about equipment status, and thus a larger weight. For example, under high-speed filling conditions, features related to rotational speed, such as frequency amplitude, typically have a larger variance contribution; while under standby conditions, features reflecting foundation vibration levels, such as root mean square (RMS) values, typically have a larger variance contribution. The specific weight values ​​for each feature need to be calculated and determined based on the variance analysis results within the actual operating condition cluster.

[0135] The aforementioned baseline feature vector, dispersion index, and low-dimensional feature value weights together constitute the basic parameters of the working condition discrimination rule base. The low-dimensional feature value weights are calculated by normalizing the contribution rates of the first three feature values ​​in step S2-2. Illustratively, if the contribution rates of the three main directions are 52.7%, 23.5%, and 13.0%, the normalized weights are approximately 0.59, 0.26, and 0.15, respectively, with a sum of 1. The actual weight values ​​depend on the specific historical sample decomposition results.

[0136] When new real-time samples enter the system, they first undergo a preprocessing step S1 to extract a 24-dimensional high-dimensional feature vector. Then, the orthogonal feature vectors obtained in step S2-2 are projected onto a 3-dimensional low-dimensional feature vector to obtain a real-time low-dimensional feature vector. Subsequently, the system calculates the weighted Euclidean distance from this real-time feature vector to all five benchmark feature vectors for different operating conditions. The aforementioned low-dimensional feature value weights are introduced during the distance calculation process, allowing the feature directions with high information content to occupy a larger proportion in the distance calculation, further amplifying the feature differences between different operating conditions.

[0137] Finally, based on the principle of minimum weighted distance, real-time samples are assigned to the nearest working condition cluster to complete real-time working condition identification. If the minimum weighted Euclidean distance from a real-time sample to the center of all working condition references is greater than twice the dispersion of the corresponding working condition, the sample is determined to belong to an unknown working condition or an abnormal state.

[0138] Step S3 further includes the following steps:

[0139] Step S3-1: For each identified operating condition cluster, analyze the correlation between the vibration characteristics of the slewing bearing and the process parameters, and screen out the correlation characteristics that are sensitive to changes in the operating condition and can reflect the health status of the equipment.

[0140] Specifically, the correlation between each vibration characteristic and process parameters such as filling speed, cap torque, turntable angle position, and filling pressure is calculated, and the correlation is measured by the Pearson correlation coefficient.

[0141] The Pearson correlation coefficient is a commonly known algorithm in statistics and a routine technique for analyzing the correlation between fault characteristics and process parameters. This invention does not improve the Pearson correlation coefficient algorithm itself, but applies it to the feature association screening stage of multi-condition fault diagnosis of slewing bearings in filling machines. The specific calculation method will not be elaborated here.

[0142] An illustrative Pearson correlation coefficient threshold is set as follows: vibration features with an absolute Pearson correlation coefficient greater than 0.7 are retained, while features with an absolute value less than 0.7 are considered redundant and are discarded.

[0143] Vibration features with a correlation level higher than a set threshold are retained, while features with a correlation level lower than the threshold are considered redundant and discarded, forming the correlation feature vector for this operating condition. This correlation feature vector retains only features closely related to the current operating condition, allowing the health model to focus more on fault information relevant to the operating condition and reducing the interference of irrelevant features on fault judgment.

[0144] Preferably, after the correlation feature screening is completed, p effective correlation features are obtained under this working condition. For a single observation sample, its values ​​on the p effective correlation features are arranged in a fixed order to form a 1×p-dimensional correlation feature vector, which represents the coordinate position of the current sample in the health feature space of this working condition.

[0145] Preferably, the threshold value is set based on the statistical distribution of historical data, and preferably a value between 0.7 and 0.9.

[0146] Furthermore, a multicollinearity test is performed on the features in the associated feature vector. A threshold of 0.9 for multicollinearity is preferred. If the Pearson correlation coefficient between two features exceeds this threshold, the feature with the higher correlation to the process parameters is retained to avoid feature duplication that could lead to model instability. The dimensionality of the associated feature matrix varies depending on the complexity of the operating conditions; 3 to 5 features are retained for simple conditions, and 8 to 12 features are retained for complex conditions.

[0147] For a historical set of healthy samples, the associated feature vectors corresponding to each sample are stacked row-wise to form an n×p-dimensional associated feature matrix. Each row of this matrix corresponds to the associated feature vector of a sample, and each column corresponds to the value distribution of an associated feature across all historical samples. This associated feature matrix serves as the statistical basis for subsequent calculations of the mean, standard deviation, and non-fixed health boundary of each associated feature.

[0148] Step S3-2: Based on the correlation feature matrix of historical health samples, establish an independent health boundary for each operating condition cluster.

[0149] Unlike traditional fixed threshold methods, this invention sets a non-fixed health boundary for each associated feature under each working condition. Specifically, for each associated feature under each working condition, the mean and standard deviation of historical healthy samples are calculated, and the mean plus or minus a certain multiple of the standard deviation is used as the initial health boundary. The multiple is determined based on the dispersion index of historical samples under that working condition.

[0150] Schematic illustration: The non-fixed health boundary calculation method is as follows In the formula: μ f σ is the sample mean of the associated feature f under this working condition; f Let f be the sample standard deviation of the correlation feature under this operating condition; k f This is the boundary multiplier coefficient, with a value range of 2 to 3; the left boundary is the lower limit of health, and the right boundary is the upper limit of health.

[0151] Preferably, when the operating condition dispersion index is higher than the average dispersion of the entire sample, k f Take 3; when the operating condition dispersion index is lower than the average dispersion of the entire sample, k f Take 2.

[0152] Larger factors are used for conditions with high dispersion, and smaller factors are used for conditions with low dispersion. For example, a factor of 3 is used for conditions with dispersion index higher than the average dispersion of the entire sample, and a factor of 2 is used for conditions with dispersion index lower than the average.

[0153] Step S3-3: Considering the gradual change in speed and load during the switching of operating conditions in the filling machine, directly using the strict health boundary of the target operating condition for judgment can easily lead to transient false alarms during the transition period. This invention establishes a transition zone between the health boundaries of adjacent operating conditions. When the system determines that the equipment is in the process of switching operating conditions, it performs transition calculations based on the health boundaries of the source and target operating conditions to form a health judgment standard for the transition period.

[0154] Specifically, the system identifies operating condition switching events by monitoring the rate of change of process parameters. When the change in filling speed or cap torque per unit time exceeds a set threshold, the system determines that an operating condition switch has begun. The identifiers of the source and target operating conditions are recorded, and the health boundaries of the two conditions are extracted. The width of the transition zone is determined based on the distance between the baseline feature vectors of the two operating conditions and the switching rate; the greater the distance and the slower the switching, the wider the transition zone. Within the transition zone, the health boundary is a weighted fusion of the source and target operating condition boundaries, with the weights determined based on the current progress ratio during the transition process.

[0155] Schematic representation: The weighted fusion calculation method for the health boundary of the transition zone is as follows In the formula: B t B is the temporary health boundary at transition time t; s B represents the health boundary of the source operating condition.g α represents the health boundary of the target operating condition; α is the transition progress coefficient, with a value ranging from 0 to 1.

[0156] The method for calculating the transition progress coefficient is as follows: , where v current v represents the current filling speed. s v is the reference speed for the source operating condition. g The target operating condition reference speed is used. When α=0, the source operating condition boundary is fully adopted; when α=1, the target operating condition boundary is fully adopted.

[0157] The duration of the transition zone is dynamically determined based on the switching rate, with a typical value of 10 to 30 seconds. The transition zone ends when the process parameters stabilize within the allowable deviation range of the target operating condition baseline value, and this stable state continues for a preset confirmation time (preferably 5 seconds). In this case, the transition zone ends, and the system switches to the target operating condition health model.

[0158] The progress ratio is calculated based on the relative position of the current process parameters with the source and target operating conditions and their benchmark process parameters.

[0159] For example, if the current filling speed has completed 60% of the speed change from the source condition to the target condition, then the target condition boundary weight is 60% and the source condition boundary weight is 40%.

[0160] Steps S3-4: To facilitate on-site maintenance personnel's intuitive understanding of equipment health status, this invention integrates the health status of multiple related features into a single comprehensive health index. Specifically, for each related feature, the feature deviation of its current value relative to the health boundary of that operating condition is calculated. The feature deviation is calculated with the center of the health boundary as the reference, and the ratio of the distance from the current value to the boundary to the half-width of the boundary is used as the feature deviation index.

[0161] The characteristic deviation is calculated as follows In the formula, D f x represents the feature deviation of the associated feature f, with a value range of [0,2]. f μ is the current value of the associated feature f; f The health boundary center of the associated feature f under this working condition; k f σ f Let f be the half-width of the health boundary associated with the feature f.

[0162] When D f When =0, it indicates that the eigenvalue is at the center of the health boundary, and the equipment is in optimal condition; D f When D = 1, it indicates that the characteristic value has reached the health boundary and the equipment status is abnormal; when D f When D > 1, it indicates that the characteristic value exceeds the health boundary, and the equipment has a potential for failure. fWhen the value is greater than 2, it is truncated to 2 to avoid excessive influence of individual extreme characteristic values ​​on the comprehensive health index.

[0163] A feature deviation of 0 indicates the area is at the center of the health boundary, a feature deviation of 1 indicates the area is touching the health boundary, and a feature deviation greater than 1 indicates the area is outside the health boundary. The feature deviations of each associated feature are weighted and fused according to the weight of that feature in the associated feature matrix to obtain the comprehensive health index.

[0164] Comprehensive Health Index In the formula: HI is the comprehensive health index, ranging from 0 to 100; p is the number of associated features under this working condition; w f Let f be the weight of the associated feature f, with values ​​ranging from 0 to 1, and satisfying the following conditions: ;D f denoted as the feature deviation of the associated feature f.

[0165] weight w f The weight of a feature is directly proportional to the significance of its change in historical failure samples. The significance of the change is obtained by calculating the ratio of the feature value change before and after the failure. The more significant the change, the higher the weight of the feature, making the comprehensive health index more sensitive to early failures.

[0166] An illustrative comprehensive health index classification standard is as follows: HI≥85 indicates a healthy state, 70≤HI<85 indicates mild deterioration, 50≤HI<70 indicates moderate deterioration, and HI<50 indicates severe deterioration.

[0167] The comprehensive health index is normalized to a range of 0 to 100, where 100 indicates the equipment is in perfect health and 0 indicates the equipment is in a critical fault state. Preferably, the feature weights are determined based on the significance of the feature's changes in historical fault samples. The significance of the change is obtained by calculating the ratio or slope of the change in the feature value before and after the fault. The more significant the change, the higher the weight of the feature, making the comprehensive health index more sensitive to early faults. The update frequency of the comprehensive health index is synchronized with the data acquisition window, achieving near real-time quantitative assessment of the equipment's health status.

[0168] Step S4 further includes the following steps:

[0169] Step S4-1: The system continuously collects vibration and process data within a set time window. The length of the time window is the same as the length of the sliding window in step S1.

[0170] The real-time acquired data undergoes the same preprocessing procedure as in step S1, extracting real-time high-dimensional feature vectors and projecting them onto the low-dimensional orthogonal space constructed in step S2. The weighted Euclidean distance from the real-time sample to the baseline feature vectors of each operating condition is calculated, using the same method as in steps S2-4. The current operating condition category is identified based on the minimum weighted distance principle. If a real-time sample falls into the ambiguity region at the boundary between two operating conditions (i.e., the distance difference to the baseline feature vectors of two adjacent operating conditions is less than the set tolerance), the health models of the two adjacent operating conditions are activated simultaneously, entering a dual-operating-condition monitoring mode.

[0171] The boundary ambiguity region refers to the region where the difference in the weighted Euclidean distance from the real-time sample to the two adjacent working condition reference feature vectors is less than a set tolerance, and the distance from the real-time sample to both working condition reference feature vectors is less than a preset multiple of the corresponding working condition dispersion index.

[0172] In dual-condition monitoring mode, the comprehensive health index is calculated for both conditions, and the lower comprehensive health index is taken as the current judgment result to ensure the conservatism and reliability of the monitoring. Dual-condition monitoring mode continues until the real-time sample clearly falls within a certain condition cluster, after which it exits.

[0173] Step S4-2: Input the real-time correlated feature values ​​into the health model of the current operating condition, and calculate the comprehensive health index according to the method described in Step S3-4. At the same time, calculate the feature deviation of each correlated feature relative to the health boundary of the operating condition.

[0174] When the overall health index falls below the set alarm threshold or any single feature exceeds the health boundary, the system records the timestamp of the abnormal event, the current operating condition category, the overall health index value, and the name of the feature that exceeded the limit.

[0175] Furthermore, the system maintains a comprehensive health index trend sequence within a sliding time window. The length of the comprehensive health index trend sequence is set according to the equipment degradation rate, preferably between 1 hour and 24 hours.

[0176] The overall health index trend sequence is used to determine whether the equipment's health status is continuously deteriorating or experiencing momentary fluctuations. Alarms are only issued for continuously deteriorating conditions to avoid misjudgments caused by sudden shocks. Continuous deterioration is determined by a monotonically decreasing overall health index across multiple consecutive time windows or by the overall health index falling below an alarm threshold across multiple consecutive windows. Monotonically decreasing is determined by calculating the difference between the overall health indexes of adjacent windows; if all consecutive differences are negative, it is considered a monotonically decreasing condition. Momentary fluctuations are determined by calculating the variance of the overall health index trend sequence; if the variance exceeds a set threshold and the overall health index quickly recovers to the normal range, it is considered a momentary fluctuation.

[0177] The maintenance of the comprehensive health index trend series adopts the first-in, first-out principle, where old data is removed when new data enters, thus maintaining a constant series length.

[0178] The method for calculating the equipment degradation rate is as follows: In the formula: v d The rate of equipment degradation is expressed as the comprehensive health index per hour, and the value is negative; t w This is the timestamp for the w-th time window, in hours; HI w The comprehensive health index for the w-th time window; The average time of the time window; is the average value of the comprehensive health index; N is the length of the comprehensive health index trend series, with a value ranging from 10 to 20.

[0179] Through trend analysis and comparison of multiple sets of sequences with different lengths, when the sequence length is less than 10, the calculated degradation rate fluctuates significantly, and the false positive rate increases markedly; when the sequence length is greater than 20, the computational load increases and the sensitivity to early degradation decreases. Considering both computational stability and early degradation detection capability, a sequence length of 15 time windows is optimal.

[0180] Remaining lifetime prediction method: Based on the degradation rate, the time when the equipment reaches the severe degradation threshold can be predicted. The prediction method is as follows: , among which, HI current This is the current comprehensive health index, in which T r To predict remaining lifespan (in hours), v d This refers to the rate of equipment degradation.

[0181] Step S4-3: This invention establishes a multi-level alarm mechanism.

[0182] To illustrate, when the comprehensive health index is between 70 and 85, it is considered mild degradation, and the system issues a warning, suggesting that maintenance personnel increase monitoring frequency and pay attention to changes in relevant process parameters. The equipment can continue to operate at this time, but close monitoring is necessary. When the comprehensive health index is between 50 and 70, it is considered moderate degradation, and the system issues a maintenance alarm, suggesting planned maintenance and preparation of spare parts. The equipment should be shut down for inspection within the nearest maintenance window. When the comprehensive health index is below 50 or a sudden change occurs in the vibration characteristic amplitude, it is considered severe degradation, and the system issues an emergency shutdown alarm, suggesting immediate shutdown for inspection to prevent secondary damage.

[0183] like Figure 4 As shown, in another embodiment of the present invention, a multi-condition fault diagnosis system for the slewing bearing of a filling machine is illustrated, comprising:

[0184] The data acquisition unit is used to collect equipment vibration signals at the outer ring raceway of the slewing bearing of the filling machine and the bearing housing of the drive motor, and to extract process parameters such as filling speed, cap torque, turntable angle position, and filling pressure from the filling machine control system in real time.

[0185] The signal preprocessing unit receives the vibration signal and process parameters from the data acquisition unit, performs anti-aliasing filtering, resampling, timestamp alignment and sliding window segmentation on the vibration signal, performs validity verification and outlier removal on the process parameters, and outputs the preprocessed multi-source heterogeneous dataset.

[0186] The feature extraction and dimensionality reduction unit receives the preprocessed multi-source heterogeneous dataset from the signal preprocessing unit, extracts time-domain features, frequency-domain features, and time-frequency-domain features from the vibration signal to construct a high-dimensional feature vector, and uses the feature space orthogonal decomposition method to perform dimensionality reduction processing on the high-dimensional feature vector to reduce the dimensionality and output a low-dimensional feature vector.

[0187] The working condition division unit receives the low-dimensional feature vector output by the feature extraction and dimensionality reduction unit, and uses the density peak search method to identify the center of the dense region of the historical health sample distribution as the working condition cluster center, and divides the sample space into several non-overlapping working condition clusters.

[0188] The benchmark construction unit receives the working condition cluster division result from the working condition division unit, calculates the geometric center of each working condition cluster as the benchmark feature vector of that working condition, and calculates the average distance from the samples within the cluster to the benchmark feature vector as the dispersion index.

[0189] The health model construction unit receives the baseline feature vector and dispersion index from the baseline construction unit, calculates the Pearson correlation coefficient between each vibration feature and the process parameters to screen related features, retains vibration features with a correlation degree higher than a set threshold to form a related feature vector; stacks the related feature vectors of multiple historical health samples row by row to form a related feature matrix; based on the related feature matrix, establishes an independent non-fixed health boundary for each working condition cluster, establishes a transition zone between the non-fixed health boundaries of adjacent working conditions, and outputs health model parameters.

[0190] The comprehensive health index calculation unit receives the health model parameters from the health model construction unit and integrates the health status of each related feature into a single comprehensive health index.

[0191] The real-time monitoring unit is connected to the feature extraction and dimensionality reduction unit, the benchmark construction unit, the health model construction unit, and the comprehensive health index calculation unit, respectively. It is used to project the real-time collected data onto the low-dimensional feature vector after preprocessing and feature extraction, calculate the weighted Euclidean distance from the real-time sample to the benchmark feature vector of each working condition to identify the current working condition category, and simultaneously activate the health models of two adjacent working conditions in the fuzzy area of ​​the working condition boundary, and take the comprehensive health index with the lower value as the current judgment result.

[0192] The trend analysis and prediction unit is used to receive the comprehensive health index from the comprehensive health index calculation unit, maintain the trend sequence of the comprehensive health index within the sliding time window, determine whether the health status of the equipment is continuously deteriorating or fluctuating instantaneously through trend analysis, and predict the remaining lifespan of the equipment based on the deterioration rate.

[0193] The alarm unit is connected to the comprehensive health index calculation unit and the trend analysis and prediction unit, respectively, and is used to establish a multi-level alarm mechanism based on the range of the comprehensive health index.

[0194] It will be apparent to those skilled in the art that the present invention is not limited to the details of the exemplary embodiments described above, and that the invention can be implemented in other specific forms without departing from its spirit or essential characteristics. Therefore, the embodiments should be considered in all respects as exemplary and non-limiting, and the scope of the invention is defined by the appended claims rather than the foregoing description. Thus, all variations falling within the meaning and scope of equivalents of the claims are intended to be included within the present invention. No reference numerals in the claims should be construed as limiting the scope of the claims.

[0195] Furthermore, it should be understood that although this specification describes embodiments, not every embodiment contains only one independent technical solution. This narrative style is merely for clarity. Those skilled in the art should consider the specification as a whole, and the technical solutions in each embodiment can also be appropriately combined to form other embodiments that can be understood by those skilled in the art.

Claims

1. A rotary bearing multi-working condition fault diagnosis method for a filling machine, characterized in that, Includes the following steps: Step S1: Multi-source heterogeneous data acquisition and signal preprocessing: Vibration sensors are installed at the outer ring raceway of the slewing bearing of the filling machine and at the drive motor bearing housing to collect equipment vibration signals. At the same time, a data interface is configured in the filling machine control system to extract the following process parameters in real time: filling speed, cap torque, turntable angle position, filling pressure, and filling volume. The vibration signals are subjected to anti-aliasing filtering, resampling, timestamp alignment, and sliding window segmentation. The process parameters are validated and outliers are removed to obtain the preprocessed multi-source heterogeneous dataset. Step S2: Working condition division and benchmark construction based on orthogonal decomposition of feature space: Time-domain features, frequency-domain features, and time-frequency-domain features are extracted from the preprocessed vibration signal to construct a high-dimensional feature vector. The orthogonal decomposition of feature space is used to reduce the dimensionality of the high-dimensional feature vector to obtain a low-dimensional feature vector. In the low-dimensional feature vector, the density peak search method is used to identify the center of the dense region of the historical healthy sample distribution as the working condition cluster center, and the sample space is divided into several non-overlapping working condition clusters. For each working condition cluster, its geometric center is calculated as the benchmark feature vector of that working condition, and the average distance from the samples in the cluster to the benchmark feature vector is calculated as the dispersion index. Step S3: Health model construction and threshold determination based on working condition correlation features: For each identified working condition cluster, calculate the Pearson correlation coefficient between each vibration feature and the process parameter to screen the correlation features, and retain the vibration features with a correlation degree higher than the set threshold to form the correlation feature vector. A correlation feature matrix is ​​constructed by stacking the correlation feature vectors of multiple historical health samples row by row. Based on the correlation feature matrix, an independent non-fixed health boundary is established for each operating condition cluster. The non-fixed health boundary is determined by adding or subtracting a certain number of standard deviations from the mean of the correlation features, with the multiple determined according to the dispersion index of the operating condition. A transition zone is established between the non-fixed health boundaries of adjacent operating conditions. When an operating condition switching event is detected, a health judgment standard for the transition period is formed by weighted fusion based on the health boundaries of the source and target operating conditions. The health status of each correlation feature is fused into a single comprehensive health index, which is calculated by weighting the feature deviation of each correlation feature relative to the non-fixed health boundary. Step S4: Real-time online monitoring and graded fault diagnosis: Perform the same preprocessing procedure as in step S1 on the real-time collected data, extract the real-time high-dimensional feature vector and project it onto the low-dimensional feature vector constructed in step S2, calculate the weighted Euclidean distance from the real-time sample to the benchmark feature vector of each working condition to identify the current working condition category; input the real-time associated feature values ​​into the health model of the current working condition, and calculate the comprehensive health index and the feature deviation of each associated feature; Maintain the trend sequence of the comprehensive health index within the sliding time window, determine whether the health status of the equipment is continuously deteriorating or fluctuating instantaneously through trend analysis, and predict the remaining lifespan of the equipment based on the deterioration rate; establish a multi-level alarm mechanism based on the range of the comprehensive health index.

2. The rotary bearing multi-condition fault diagnosis method for a filling machine according to claim 1, characterized in that, In step S1, two axes of the vibration sensor are arranged in a horizontal plane perpendicular to each other, and the other axis is arranged in a vertical direction to fully sense the radial and axial vibration responses of the slewing bearing.

3. The rotary bearing multi-condition fault diagnosis method for a filling machine according to claim 2, characterized in that, In step S1, the window length of the sliding window segment is set as an integer multiple of the rotation period of the slewing bearing, and the window movement step size is set as the overlap rate of the sliding window.

4. The multi-condition fault diagnosis method for the slewing bearing of a filling machine according to claim 3, characterized in that, In step S2, the time-domain features include root mean square value, peak value, waveform value, impulse value, margin value, skewness value, and kurtosis value; the frequency-domain features include rotational frequency amplitude, rotational frequency harmonic amplitude, characteristic frequency amplitude of slewing bearing inner ring fault, characteristic frequency amplitude of slewing bearing outer ring fault, characteristic frequency amplitude of rolling element fault, characteristic frequency amplitude of cage fault, and energy proportion of preset frequency bands; the time-frequency domain features include the energy of each frequency band and its energy entropy value obtained by wavelet packet decomposition.

5. The multi-condition fault diagnosis method for the slewing bearing of a filling machine according to claim 4, characterized in that, In step S2, the density peak search method includes: calculating the local density of each sample point, wherein the local density is defined as the number of other sample points within a certain distance of the sample point; calculating the distance from each sample point to the nearest sample point with higher local density; constructing a decision graph with local density as the horizontal axis and the distance as the vertical axis, and identifying the density peak point located in the upper right corner region as the center of the working condition cluster.

6. The multi-condition fault diagnosis method for the slewing bearing of a filling machine according to claim 5, characterized in that, In step S2, for a sample located at the boundary of two working condition clusters, its membership degree to the center of the adjacent working condition cluster is calculated. The membership degree is calculated by weighting the inverse of the distance from the sample to each center and stored in the form of a membership degree vector for subsequent transition zone calculation.

7. The multi-condition fault diagnosis method for the slewing bearing of a filling machine according to claim 6, characterized in that, In step S3, the threshold value of the Pearson correlation coefficient is between 0.7 and 0.9; a multicollinearity test is performed on the features in the associated feature vector, and when the absolute value of the Pearson correlation coefficient between two features is greater than the multicollinearity judgment threshold, the feature with a higher degree of correlation with the process parameters is retained.

8. The multi-condition fault diagnosis method for the slewing bearing of a filling machine according to claim 7, characterized in that, In step S3, the width of the transition zone is determined based on the distance between the reference feature vectors of the two operating conditions and the switching rate. Within the transition zone, the health boundary adopts the weighted fusion result of the source operating condition and the target operating condition boundary, and the weight is determined based on the progress ratio of the current moment in the transition process.

9. The multi-condition fault diagnosis method for the slewing bearing of a filling machine according to claim 8, characterized in that, In step S4, when the real-time sample falls into the fuzzy boundary region between the two working conditions, the health models of the two adjacent working conditions are activated simultaneously to enter the dual-working-condition monitoring mode, and the comprehensive health index under the two working conditions is calculated respectively. The comprehensive health index with the lower value is taken as the current judgment result.

10. A multi-condition fault diagnosis system for the slewing bearing of a filling machine, characterized in that, include: The data acquisition unit is used to collect equipment vibration signals at the outer ring raceway of the slewing bearing and the drive motor bearing housing of the filling machine, and to extract the following process parameters from the filling machine control system in real time: filling speed, cap torque, turntable angle position, filling pressure, and filling volume. The signal preprocessing unit receives the vibration signal and process parameters from the data acquisition unit, performs anti-aliasing filtering, resampling, timestamp alignment and sliding window segmentation on the vibration signal, performs validity verification and outlier removal on the process parameters, and outputs the preprocessed multi-source heterogeneous dataset. The feature extraction and dimensionality reduction unit receives the preprocessed multi-source heterogeneous dataset from the signal preprocessing unit, extracts time-domain features, frequency-domain features, and time-frequency-domain features from the vibration signal to construct a high-dimensional feature vector, and uses the feature space orthogonal decomposition method to perform dimensionality reduction processing on the high-dimensional feature vector to reduce the dimensionality and output a low-dimensional feature vector. The working condition division unit receives the low-dimensional feature vector output by the feature extraction and dimensionality reduction unit, and uses the density peak search method to identify the center of the dense region of the historical health sample distribution as the working condition cluster center, and divides the sample space into several non-overlapping working condition clusters. The benchmark construction unit receives the working condition cluster division result from the working condition division unit, calculates the geometric center of each working condition cluster as the benchmark feature vector of that working condition, and calculates the average distance from the samples within the cluster to the benchmark feature vector as the dispersion index. The health model construction unit receives the baseline feature vector and dispersion index from the baseline construction unit, calculates the Pearson correlation coefficient between each vibration feature and the process parameters to screen related features, and retains vibration features with a correlation degree higher than a set threshold to form a related feature vector. The correlation feature vectors of multiple historical health samples are stacked row by row to form a correlation feature matrix; based on the correlation feature matrix, an independent non-fixed health boundary is established for each working condition cluster, and a transition zone is established between the non-fixed health boundaries of adjacent working conditions to output the health model parameters. The comprehensive health index calculation unit receives the health model parameters from the health model construction unit and integrates the health status of each related feature into a single comprehensive health index. The real-time monitoring unit is connected to the feature extraction and dimensionality reduction unit, the benchmark construction unit, the health model construction unit, and the comprehensive health index calculation unit, respectively. It is used to project the real-time collected data onto the low-dimensional feature vector after preprocessing and feature extraction, calculate the weighted Euclidean distance from the real-time sample to the benchmark feature vector of each working condition to identify the current working condition category, and simultaneously activate the health models of two adjacent working conditions in the fuzzy area of ​​the working condition boundary, and take the comprehensive health index with the lower value as the current judgment result. The trend analysis and prediction unit is used to receive the comprehensive health index from the comprehensive health index calculation unit, maintain the trend sequence of the comprehensive health index within the sliding time window, determine whether the health status of the equipment is continuously deteriorating or fluctuating instantaneously through trend analysis, and predict the remaining lifespan of the equipment based on the deterioration rate. The alarm unit is connected to the comprehensive health index calculation unit and the trend analysis and prediction unit, respectively, and is used to establish a multi-level alarm mechanism based on the range of the comprehensive health index.

Images