Intelligent monitoring method for mechanical vibration in bearing operation process

Abstract

The application belongs to the technical field of intelligent monitoring, and particularly relates to a mechanical vibration intelligent monitoring method in a bearing operation process, comprising the following steps: carrying out multi-layer wavelet packet decomposition on a bearing original vibration signal, extracting a target frequency band sub-signal and an instantaneous energy sequence thereof according to a kurtosis maximization principle; determining a non-extensive entropy parameter from a kurtosis coefficient, and determining a candidate box width search interval in combination with a quartile range and a weighting factor; calculating a penalty Tsallis entropy value containing a model complexity penalty term, and optimizing a maximum penalty Tsallis entropy and an optimal box width through a coarse search and a fine search combination strategy; constructing a three-dimensional state feature vector and calculating a Mahalanobis distance square, and determining a bearing operation state transition through comparison with a monitoring threshold value. The application realizes self-adaptive adjustment of a monitoring scale according to signal pulse characteristics, and improves the accuracy and stability of bearing early state transition perception in a strong background noise environment.

Description

Technical Field

[0001] This invention relates to the field of intelligent monitoring technology. More specifically, this invention relates to an intelligent monitoring method for mechanical vibration during bearing operation. Background Technology

[0002] As a core supporting component in modern rotating machinery, the operating status of rolling bearings is crucial to the safety, stability, and efficiency of the entire mechanical system. Under long-term and variable operating conditions, bearings are susceptible to damage from fatigue, wear, and impact. Failure can lead to equipment downtime, production interruptions, and potentially safety accidents and economic losses. Real-time and accurate condition monitoring and fault diagnosis of bearings during operation have significant engineering application value. Among various monitoring methods, mechanical vibration signal analysis has become one of the most important methods for bearing condition monitoring due to its rich information on equipment operating status and ease of online acquisition. However, because bearings are affected by internal mechanical structure interactions and external environmental noise during actual operation, the collected mechanical vibration signals often exhibit highly nonlinear, non-stationary, and time-varying characteristics. Especially when bearings are in the early stages of degradation or when their operating state is just beginning to change, the weak characteristic information indicating state changes is easily submerged in background noise and modulation signals, posing challenges to vibration signal analysis and feature extraction.

[0003] To extract sensitive features representing bearing operating conditions from vibration signals, nonlinear dynamic analysis methods based on information entropy have been applied to the field of mechanical condition monitoring. Information entropy can assess the complexity and uncertainty of signal time series, demonstrating unique advantages in revealing the implicit dynamic evolution laws of the system. Existing condition monitoring methods, when estimating signal probability density or partitioning state space, often rely on fixed empirical formulas or manual presets for key analysis parameters, lacking a mechanism to optimize based on the transient energy distribution characteristics of the signal itself. This results in the inability to guarantee the precision and rationality of feature extraction under different operating conditions. Conventional feature extraction indices often fail to simultaneously address the identification of system nonlinear characteristics and control model complexity, easily reducing sensitivity and anti-interference capabilities to subtle changes in bearing state due to information redundancy or single-dimensionality. There is an urgent need for a monitoring method that can combine the inherent statistical characteristics of the signal, optimize the information feature extraction scale, and integrate multi-dimensional sensitive feature indices to improve the accuracy and stability of bearing operating state transformation perception. Summary of the Invention

[0004] To address the technical problems of bearing vibration signals being easily overwhelmed by noise due to their high nonlinearity, and the lack of adaptive optimization mechanisms due to the reliance on manually preset monitoring parameters, which makes it difficult to balance feature precision and model complexity, resulting in insufficient accuracy and stability in sensing bearing state changes, this invention provides an intelligent monitoring method for mechanical vibration during bearing operation.

[0005] This invention provides an intelligent monitoring method for mechanical vibration during bearing operation, comprising: performing multi-level wavelet packet decomposition on the acquired original bearing vibration signal to obtain sub-signals of each frequency band; extracting the target frequency band sub-signals according to the kurtosis maximization principle; calculating the instantaneous energy sequence of the target frequency band sub-signals; determining the non-extensive entropy parameter from the kurtosis coefficient of the instantaneous energy sequence through a monotonic function relationship; using the reference box width calculated based on the interquartile range of the instantaneous energy sequence as a reference, expanding the reference box width using a weighting factor proportional to the absolute deviation of the non-extensive entropy parameter and 1, and determining the candidate box width search interval; for any candidate box width within the candidate box width search interval, performing histogram statistical normalization on the instantaneous energy sequence to obtain... The discrete probability distribution is used to calculate the Tsallis entropy value under the current candidate bin width using non-extensive entropy parameters. The penalized Tsallis entropy value is obtained by subtracting the model complexity penalty term proportional to the number of data bins. Within the candidate bin width search interval, a search strategy combining coarse search and fine search is used to optimize the maximum penalized Tsallis entropy and the corresponding optimal bin width. A three-dimensional state feature vector is constructed, consisting of the maximum penalized Tsallis entropy, the optimal bin width, the ratio of the instantaneous energy sequence standard deviation, and the non-extensive entropy parameters. The squared Mahalanobis distance between the three-dimensional state feature vector and the center of the preset health state feature cluster is calculated. When the squared Mahalanobis distance is greater than the monitoring threshold, the bearing operating state is determined to have changed.

[0006] By adopting the above technical solution, this invention introduces an adaptive determination mechanism for non-extensive entropy parameters and combines it with a Tsallis entropy optimization strategy based on model complexity penalty to construct a three-dimensional state feature vector integrating energy distribution, scale resolution, and statistical evolution characteristics. The kurtosis coefficient of the instantaneous energy sequence is used to dynamically determine the value of the non-extensive entropy parameters, and this is used to simultaneously expand the candidate bin search range. This allows the monitoring system to automatically adjust the observation scale according to the inherent impulse characteristics of the signal, resolving the contradiction between the depth of information mining and the model's generalization ability under fixed parameters in traditional methods. Furthermore, by using the squared Mahalanobis distance to correlate and fuse three-dimensional heterogeneous features, including non-extensive entropy parameters, dimensional differences and feature redundancy are eliminated, improving the accuracy and stability of bearing early state transition perception under strong background noise.

[0007] Preferably, the step of performing multi-level wavelet packet decomposition on the acquired original bearing vibration signal to obtain sub-band signals includes: using a vibration acceleration sensor to collect vibration acceleration data of the bearing as the original vibration signal; and using the db4 wavelet basis function to perform three-level wavelet packet decomposition on the original vibration signal to reconstruct the sub-band signals of eight independent frequency bands.

[0008] Preferably, the step of extracting the target frequency band sub-signal according to the kurtosis maximization principle and calculating the instantaneous energy sequence of the target frequency band sub-signal includes: calculating the fourth-order normalized central moments of the frequency band sub-signals of the eight independent frequency bands to obtain their respective kurtosis values, and selecting the frequency band sub-signal corresponding to the maximum kurtosis value as the target frequency band sub-signal; calculating the Teager-Kaiser energy value point by point by subtracting the product of the values ​​of the previous and next adjacent sampling points from the square of the current sampling point value of the target frequency band sub-signal, and constructing the instantaneous energy sequence from all the calculated energy values.

[0009] By adopting the above technical solution, this invention selects the target frequency band based on the principle of kurtosis maximization, ensuring that the subsequent analysis object is always the signal component with the richest impact characteristics; in conjunction with the Teager-Kaiser energy operator to extract instantaneous energy from the target frequency band, and by utilizing its nonlinear and nonlocal characteristics, it amplifies the transient energy mutations in the signal caused by weak pitting, and enhances the nonlinear dynamic characteristics implicit in the original vibration sequence.

[0010] Preferably, the step of determining the non-extensive entropy parameter from the kurtosis coefficient of the instantaneous energy sequence through a monotonic function relationship includes: extracting the sample mean and sample standard deviation of the instantaneous energy sequence; calculating the ratio of the centered fourth moment to the fourth power of the sample standard deviation to obtain the kurtosis coefficient of the instantaneous energy sequence; substituting the kurtosis coefficient into the monotonic function relationship to calculate the corresponding non-extensive entropy parameter, wherein the monotonic function relationship is such that the non-extensive entropy parameter is equal to 0.5 times the kurtosis coefficient power of the natural constant e.

[0011] By adopting the above technical solution, this invention establishes an exponential mapping relationship between the kurtosis coefficient of the instantaneous energy sequence and the non-extensive entropy parameter, enabling the entropy measure to adaptively evolve according to the degree of deviation of the signal from the Gaussian distribution. When the bearing generates an impact pulse that causes an increase in kurtosis, the automatically increasing value of the non-extensive entropy parameter enhances the sensitivity of the entropy system to extreme fluctuation events, thereby achieving a precise assessment of the non-extensive characteristics of the system in a physical sense.

[0012] Preferably, the step of using the reference bin width calculated based on the interquartile range of the instantaneous energy sequence as a reference includes: sorting the data of the instantaneous energy sequence from smallest to largest, calculating the difference between the 75th percentile value and the 25th percentile value to obtain the interquartile range; calculating the reference bin width using the Freedman-Diaconis rule, wherein the reference bin width is equal to a constant 2 multiplied by the interquartile range and then multiplied by the length of the instantaneous energy sequence. Power of 1.

[0013] Preferably, the step of expanding the baseline bin width using a weighting factor proportional to the absolute deviation of the non-extensive entropy parameter and 1 to determine the candidate bin width search interval includes: calculating the absolute value of the difference between the non-extensive entropy parameter and the constant 1; multiplying the absolute value by a preset scaling coefficient and adding an offset constant to obtain the weighting factor; determining the quotient of the baseline bin width divided by the weighting factor as the lower bound of the candidate bin width search interval; and determining the product of the baseline bin width multiplied by the weighting factor as the upper bound of the candidate bin width search interval, thus forming a closed interval as the candidate bin width search interval.

[0014] Preferably, the step of performing histogram statistical normalization on the instantaneous energy sequence to obtain a discrete probability distribution for any candidate bin width within the candidate bin width search interval, calculating the Tsallis entropy value under the current candidate bin width using a non-extensive entropy parameter, and subtracting a model complexity penalty term proportional to the number of data bins to obtain a penalized Tsallis entropy value includes: dividing the data coverage of the instantaneous energy sequence into several non-overlapping data bins with their heads and tails connected according to the candidate bin width; counting the number of sample points falling in each data bin and dividing it by the total sequence length to obtain the discrete probability distribution, and obtaining the total number of data bins; calculating the Tsallis entropy value based on the non-extensive entropy parameter and the discrete probability distribution, and subtracting the model complexity penalty term from the Tsallis entropy value to obtain the penalized Tsallis entropy value, wherein the model complexity penalty term is equal to a preset constant regularization factor multiplied by the total number of data bins.

[0015] By adopting the above technical solution, this invention establishes an evaluation mechanism similar to an information criterion by subtracting a model complexity penalty term proportional to the number of data bins from the Tsallis entropy calculation. This regularization method effectively constrains the fragmentation tendency of histogram partitioning, forcing the system to consider the simplicity of the statistical model while acquiring information gain, and preventing overfitting risks caused by excessively small bin widths at the algorithm level.

[0016] Preferably, the step of using a search strategy combining coarse and fine search within the candidate bin width search interval to optimize the maximum penalty Tsallis entropy and the corresponding optimal bin width includes: generating multiple bin width sample points proportionally in logarithmic coordinates between the lower and upper bounds of the candidate bin width search interval and calculating the corresponding penalty Tsallis entropy value for each; recording the bin width sample point with the largest penalty Tsallis entropy value among the bin width sample points as the extreme value center point; if the extreme value center point is not located on the boundary, extracting the extreme value center point and the adjacent left and right bin width sample points to construct a univariate quadratic polynomial function, setting the first derivative equal to 0 to solve for the root of the independent variable as the optimal bin width, and substituting the optimal bin width into the univariate quadratic polynomial function for calculation to obtain the maximum penalty Tsallis entropy; if the extreme value center point is located on the boundary, the extreme value center point and the corresponding penalty Tsallis entropy value are taken as the optimal bin width and the maximum penalty Tsallis entropy.

[0017] By adopting the above technical solution, this invention employs a strategy combining coarse search with logarithmic coordinates and fine search with quadratic polynomial fitting to optimize the optimization process of the non-convex penalty Tsallis entropy function. The logarithmic distribution ensures dense sampling within small-scale intervals, while polynomial interpolation provides the analytical roots of the optimal solution. This significantly reduces computational overhead while maintaining feature extraction accuracy, thus meeting the real-time requirements of online monitoring in industrial settings.

[0018] Preferably, the calculation of the squared Mahalanobis distance between the three-dimensional state feature vector and the preset health state feature cluster center includes: extracting a sample set composed of multiple historical three-dimensional state feature vectors during the normal operation phase of the bearing; calculating the mean vector of the sample set as the preset health state feature cluster center; calculating the feature covariance matrix of the sample set and inverting it to obtain the inverse covariance matrix; subtracting the preset health state feature cluster center from the measured current three-dimensional state feature vector to obtain the difference vector; and calculating the continuous matrix product of the transpose of the difference vector, the inverse covariance matrix, and the difference vector to obtain the squared Mahalanobis distance.

[0019] Preferably, determining the change in bearing operating state when the squared Mahalanobis distance is greater than the monitoring threshold includes: determining the monitoring threshold based on the chi-square distribution critical value with 3 degrees of freedom at a preset confidence level; comparing the squared Mahalanobis distance with the monitoring threshold; and triggering an early warning signal when the squared Mahalanobis distance is greater than the monitoring threshold to determine that the bearing operating state has changed from a normal state to a degraded state.

[0020] The technical solution of the present invention has the following beneficial technical effects:

[0021] This invention extracts target frequency band sub-signals and their instantaneous energy sequences through multi-layer wavelet packet decomposition and the principle of kurtosis maximization, thereby enhancing the core features of bearing operating state evolution while suppressing background noise interference. By establishing a monotonic function mapping relationship between the kurtosis coefficient and non-extensive entropy parameters, and utilizing a nonlinear expansion mechanism based on interquartile range and weighting factors to rationally define the candidate box width search interval, combined with an optimization strategy integrating coarse and fine search, and a penalized Tsallis entropy calculation method incorporating a model complexity penalty term, the computational efficiency and accuracy of optimal box width optimization are improved, and the risk of overfitting in probability density statistical modeling is effectively reduced. Furthermore, this invention constructs a three-dimensional state feature vector integrating the maximum penalized Tsallis entropy, the ratio of optimal box width to standard deviation, and non-extensive entropy parameters. Mahalanobis distance is used to measure the deviation of the current feature vector from the center of a preset healthy state feature cluster to determine state transitions, eliminating the influence of dimensional differences and correlations among multi-dimensional features, deeply characterizing the inherent dynamic distribution law of vibration signals, thereby improving the accuracy and stability of the entire bearing operation monitoring process. Attached Figure Description

[0022] Figure 1 This is a flowchart of the intelligent monitoring method for mechanical vibration during bearing operation in this invention; Figure 2 This is a schematic diagram of the Teager-Kaiser instantaneous energy sequence of the target frequency band sub-signal; Figure 3 This is a schematic diagram illustrating the evolution of the candidate bin width search interval with the non-extensive entropy parameter; Figure 4 This is a schematic diagram of the Mahalanobis distance squared over the entire lifetime under different monitoring schemes. Detailed Implementation

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

[0024] This invention discloses an intelligent monitoring method for mechanical vibration during bearing operation, referring to... Figure 1 This includes steps S1-S3: S1, extract the target frequency band sub-signal and determine the candidate box width search range.

[0025] The acquired original bearing vibration signal is decomposed into multi-level wavelet packet decomposition to obtain sub-signals of each frequency band. The target frequency band sub-signals are extracted according to the principle of kurtosis maximization, and the instantaneous energy sequence of the target frequency band sub-signals is calculated. The non-extensive entropy parameter is determined by the kurtosis coefficient of the instantaneous energy sequence through a monotonic function relationship. The reference box width calculated based on the interquartile range of the instantaneous energy sequence is used as the reference, and the reference box width is expanded using a weighting factor that is proportional to the absolute deviation of the non-extensive entropy parameter and 1 to determine the candidate box width search interval.

[0026] The original vibration signal of the bearing is acquired using a piezoelectric accelerometer. To highlight signal characteristics, the original vibration signal undergoes multi-level wavelet packet decomposition, for example, a three-level wavelet packet decomposition, selecting db4 as the wavelet basis function to obtain eight frequency band sub-signals. Subsequently, the kurtosis value of each of the eight frequency band sub-signals is calculated using a kurtosis calculation function, and the frequency band sub-signal with the largest kurtosis value is selected as the target frequency band sub-signal. Next, the instantaneous energy sequence is calculated using a feature enhancement algorithm. For example, the target frequency band sub-signal can be subjected to a Hilbert transform to obtain an analytic signal and the amplitude squared can be extracted, or the Teager-Kaiser energy operator can be used to generate the instantaneous energy sequence. Then, the kurtosis coefficient of the instantaneous energy sequence is calculated using a statistical module, and the non-extensive entropy parameter q is calculated based on the monotonic function relationship, for example, the formula is q equal to 0.5 times the kurtosis coefficient power of the natural constant e. Further, the 75% and 25% position values ​​of the instantaneous energy sequence are calculated and the difference is obtained to obtain the interquartile range, and the width of the reference box is calculated using the Freedman-Diaconis rule. Finally, the absolute value of the difference between the non-extensive entropy parameter q and the constant 1 is calculated, and the absolute value of the natural constant e raised to the power of this value is used as a weighting factor. Alternatively, the absolute value can be mapped to a weighting factor through a linear transformation, such as multiplying the absolute value by a preset scaling factor and adding a bias constant as a weighting factor. The preset scaling factor is preferably in the range of 0.05 to 0.5, and in this embodiment, it is most preferably 0.2. The bias constant is preferably in the range of 0.5 to 2.0, and in this embodiment, it is most preferably a constant 1.

[0027] It should be noted that the preset scale adjustment coefficient and bias constant are set within a certain range to ensure that the weighting factor can appropriately respond to the deviation of the non-extensive entropy parameter. If the preset scale adjustment coefficient is too large, it will cause the candidate bin width search interval to expand excessively, increasing meaningless computational overhead and potentially introducing interfering bin widths; if it is too small, the candidate bin width search interval will not expand sufficiently, easily missing the optimal resolution bin width when the device experiences heavy-tailed impact conditions. Maintaining the bias constant around 1 ensures that, under ideal Gaussian distribution conditions (i.e., when the absolute deviation approaches 0), the candidate bin width search interval can smoothly retreat and converge to near the reference bin width.

[0028] The candidate box width search interval is generated by dividing the baseline box width by the weighting factor and using the baseline box width multiplied by the weighting factor as the upper bound of the candidate box width search interval.

[0029] In an optional embodiment, the acquired original bearing vibration signal is decomposed into multi-level wavelet packet decomposition to obtain sub-signals of each frequency band, including: using a vibration acceleration sensor to collect the vibration acceleration data of the bearing as the original vibration signal; using the db4 wavelet basis function to perform 3-level wavelet packet decomposition on the original vibration signal to reconstruct the sub-signals of 8 independent frequency bands.

[0030] In practical implementation, during the stage of acquiring the original vibration signal of the bearing, a high-frequency piezoelectric vibration accelerometer with a sensitivity of 100mV / g to 500mV / g is preferably used. This sensor is threaded or magnetically attached to the vertical bearing area of ​​the bearing housing. To ensure high-frequency coverage and resolution of the signal, the sampling frequency of the data acquisition card is set. The sampling frequency ranges from 12000Hz to 25600Hz; for example, in this embodiment, the preferred sampling frequency is 12800Hz. The length L of a single data acquisition is set to 2048 to 8192 sampling points; for example, in this embodiment, the data length is 4096 sampling points. Subsequently, the db4 wavelet basis function from the Daubechies wavelet family is called to perform a 3-level discrete wavelet packet decomposition on the acquired 4096 original sampled data sequences. This decomposition process uses an orthogonal mirror filter bank to equally divide the signal's frequency band within the low-frequency to high-frequency range (0 to 6400Hz) into 8 frequency bands, each with a bandwidth of 800Hz. Then, the corresponding reconstruction coefficients are used to perform single-branch reconstruction on the signals at each node, resulting in 8 frequency band sub-signal sequences, each with a length of 4096. , where i=1,2,...,8, n=1,2,...,4096.

[0031] In an optional embodiment, the target frequency band sub-signals are extracted according to the kurtosis maximization principle, and the instantaneous energy sequence of the target frequency band sub-signals is calculated, including: calculating the fourth-order normalized central moments of the frequency band sub-signals of the eight independent frequency bands to obtain their respective kurtosis values, and selecting the frequency band sub-signal corresponding to the maximum kurtosis value as the target frequency band sub-signal; calculating the Teager-Kaiser energy value point by point by subtracting the product of the values ​​of the previous and next adjacent sampling points from the square of the current sampling point value of the target frequency band sub-signal, and constructing an instantaneous energy sequence from all the calculated energy values.

[0032] In practice, to screen target frequency band sub-signals containing fault impact characteristics, kurtosis values, i.e., fourth-order normalized central moments, are calculated for each of the eight frequency band sub-signals. The specific calculation process involves extracting the sample mean and standard deviation of the frequency band sub-signal sequence, summing the fourth power of the difference between each point's value and the mean, dividing by the total number of points, and then dividing by the fourth power of the standard deviation to obtain the kurtosis value. Assuming that the kurtosis value of a certain high-frequency node reaches the maximum of 5.6, which is higher than the normal Gaussian distribution of 3, the reconstructed sub-signal corresponding to that frequency band is selected as the target frequency band sub-signal. In the instantaneous energy sequence calculation stage, a nonlinear, nonlocal Teager-Kaiser energy operator is introduced to enhance the weak impact characteristics. Specifically, it iterates through the index n of the target frequency band sub-signal sequence, where n = 2, 3, ..., 4095, and subtracts the value of the previous adjacent sampling point from the square of the current sampling point value. Value of the next adjacent sampling point The product of these factors is used to calculate the Teager-Kaiser energy values ​​at each discrete time point. This results in an instantaneous energy sequence of length 4094, which amplifies the amplitude abrupt changes caused by weak transient impacts, laying the data foundation for subsequent extraction of non-extensive entropy measures. Combined with... Figure 2 As shown in the figure, the instantaneous energy sequence curve of the Teager-Kaiser is displayed. The horizontal axis is the sampling number n, the vertical axis is the instantaneous energy value, the curve baseline is the noise floor, and the discrete spikes correspond to the energy mutation pulses generated by the bearing failure impact.

[0033] In an optional embodiment, the non-extensive entropy parameter is determined from the kurtosis coefficient of the instantaneous energy sequence through a monotonic function relation, including: extracting the sample mean and sample standard deviation of the instantaneous energy sequence; calculating the ratio of the centered fourth moment to the fourth power of the sample standard deviation to obtain the kurtosis coefficient of the instantaneous energy sequence; substituting the kurtosis coefficient into the monotonic function relation to calculate the corresponding non-extensive entropy parameter q, where the monotonic function relation is q equal to 0.5 times the kurtosis coefficient power of the natural constant e.

[0034] In practical implementation, to determine the non-extensive entropy parameter q that can identify the non-Gaussian and fractal characteristics of the system, it is necessary to extract high-order statistical features from the previously obtained instantaneous energy sequence of length 4094. First, the accumulator is initialized to calculate the sample mean of the instantaneous energy sequence; for example, in this embodiment, the measured sample mean is... Simultaneously extract the sample standard deviation, for example, measured as... Then, the sequence is centered by subtracting the sample mean from each element to obtain a mean-free sequence. The elements of this mean-free sequence are then raised to the fourth power and summed, divided by the total sequence length of 4094 to obtain the centered fourth moment. This centered fourth moment is then divided by the fourth power of the sample standard deviation to calculate the kurtosis coefficient K of the instantaneous energy sequence. Assuming that under slightly worn operating conditions, the instantaneous energy sequence deviates more significantly from the Gaussian distribution, and internal abrupt pulses cause the calculated kurtosis coefficient K to reach 3.85.

[0035] After obtaining the kurtosis coefficient K, a specific exponential transformation model is used to map this coefficient to a non-extensive entropy parameter q, representing the degree to which the system deviates from extensive statistical mechanics. In this mapping process, as a preferred implementation, the selected monotonic function relationship is: where e is the base of the natural logarithm. This is the kurtosis coefficient. It is obtained by the microprocessor performing floating-point calculations. The result is 6.855. This mapping mechanism ensures that when the bearing damage intensifies, leading to increased signal impact and a heavy-tailed instantaneous energy distribution causing an increase in the kurtosis coefficient, the calculated non-extensive entropy parameter q increases exponentially. Since the entropy system emphasizes the contribution of high-probability events and extreme fluctuations when the non-extensive entropy parameter q is greater than 1, this nonlinear transformation method can adaptively adjust the sensitivity of the entropy measure to pulses generated by local faults. Furthermore, the entire calculation process requires only basic algebraic instructions, which is highly beneficial for the real-time performance of parameter updates within edge computing devices and for engineering implementation.

[0036] In an optional embodiment, using a reference bin width calculated based on the interquartile range of the instantaneous energy sequence as a reference, the method includes: sorting the data of the instantaneous energy sequence from smallest to largest, calculating the difference between the 75th percentile value and the 25th percentile value to obtain the interquartile range; and calculating the reference bin width using the Freedman-Diaconis rule, where the reference bin width is equal to a constant 2 multiplied by the interquartile range and then multiplied by the length of the instantaneous energy sequence. Power of 1.

[0037] In an optional embodiment, the baseline bin width is expanded using a weighting factor proportional to the absolute deviation of the non-extensive entropy parameter q from 1, and the candidate bin width search interval is determined. This includes: calculating the absolute value of the difference between the non-extensive entropy parameter q and the constant 1; multiplying the absolute value by a preset scaling factor and adding an offset constant to obtain the weighting factor; determining the quotient of the baseline bin width divided by the weighting factor as the lower bound of the candidate bin width search interval; and determining the product of the baseline bin width multiplied by the weighting factor as the upper bound of the candidate bin width search interval, thus forming a closed interval as the candidate bin width search interval.

[0038] In practical implementation, when determining the bin width for the histogram probability distribution, the first step is to extract the statistic representing the degree of data dispersion, namely the interquartile range (IQR). Specifically, the 4094-bit instantaneous energy sequence data obtained earlier is sorted using a quicksort algorithm, arranged in ascending order of numerical values. Then, the third quartile at the 75th percentile and the first quartile at the 25th percentile are located. Assuming the third quartile is measured to be 0.28 and the first quartile to be 0.12, the difference between the two is the IQR: IQR = 0.28 - 0.12 = 0.16. Next, the baseline bin width is calculated using the Freedman-Diaconis criterion. Given a total length of 4094, its The power equals 0.0625. Substituting the interquartile range into the formula... That is, the constant 2 multiplied by 0.16 and then multiplied by 0.0625, the baseline bin width is calculated to be 0.02. This baseline value can be used as the optimal bin width estimate for achieving probability density smoothing under ideal conditions.

[0039] To enable the bin width search space to adapt to unsteady and impact-heavy tail characteristics under varying operating conditions, this scheme introduces a weighted expansion mechanism linked to the non-extensive entropy parameter q. First, the absolute value of the difference between parameter q and scalar 1 in the current state is calculated. For example, if the non-extensive entropy parameter q obtained above is 6.855, then this absolute value is 5.855. Then, the absolute value of this value is raised to the power of the natural constant e as a weighting factor, or the absolute value is multiplied by a preset scale adjustment coefficient, preferably 0.2 in this embodiment, and then a bias constant is added, preferably constant 1 in this embodiment, to calculate a dimensionless weighting factor of 2.171. Finally, using the calculated baseline bin width as the anchor point, a multiplication-division nonlinear boundary expansion is performed to construct a closed interval for the candidate bin width search. The lower bound of the candidate bin width search interval is obtained by dividing the baseline bin width 0.02 by the weighting factor 2.171, resulting in a value of 0.0092. The upper bound of the candidate bin width search interval is obtained by multiplying the baseline bin width of 0.02 by a weighting factor of 2.171, resulting in a value of 0.0434. The microprocessor defines the candidate bin width search interval output to subsequent optimization algorithms as a closed interval [0.0092, 0.0434]. This strategy of using the deviation of the non-extensive entropy parameter q for nonlinear boundary expansion enables the monitoring system to provide a coarse-grained data domain with appropriate resolution in both the early and late stages of a fault. Figure 3 As shown in the figure, the curve of the candidate box width search interval varies with the non-extensive entropy parameter q. When the non-extensive entropy parameter q is greater than 1, the upper bound of the candidate box width search interval increases monotonically with q, while the lower bound decreases monotonically with q. This successfully realizes the nonlinear expansion of the candidate box width search interval under the fault heavy-tail condition.

[0040] S2, finding the maximum penalty Tsallis entropy and the optimal bin width.

[0041] For any candidate bin width within the candidate bin width search interval, the instantaneous energy sequence is histogram-normalized to obtain a discrete probability distribution. The Tsallis entropy value under the current candidate bin width is calculated using the non-extensive entropy parameter. The penalized Tsallis entropy value is obtained by subtracting the model complexity penalty term proportional to the number of data bins. Within the candidate bin width search interval, a search strategy combining coarse search and fine search is used to find the maximum penalized Tsallis entropy and the corresponding optimal bin width.

[0042] The algorithm iterates through each candidate bin width within the search interval. The number of bins is obtained by subtracting the minimum value from the maximum value of the instantaneous energy sequence and dividing by the current candidate bin width. A histogram statistical function is used to count the number of sample points in each bin, and this count is divided by the total sequence length to obtain the discrete probability distribution sequence within each bin. The Tsallis entropy value for the current candidate bin width is calculated according to the definition of Tsallis entropy. This is done by summing the terms of the discrete probability distribution sequence to the power of q, subtracting this sum from the constant 1, and then dividing the result by the difference between q and 1. To control model complexity, a preset constant regularization factor is set. In practice, the preferred value range for the preset constant regularization factor is 0.0001 to 0.005, and in this embodiment, the optimal value is 0.0005. It should be noted that this regularization factor is used to balance the model's fit and complexity. If the value is too large, it will lead to over-penalization, causing the system to tend to select very few data bins and lose local features; if the value is too small, the penalty will be insufficient, and the system may divide the data bins into fragmented bins, thus falling into overfitting. The model complexity penalty term is obtained by multiplying the preset constant regularization factor by the number of data bins. The penalized Tsallis entropy value is obtained by subtracting the model complexity penalty term from the current Tsallis entropy value. Subsequently, a search strategy combining coarse and fine search is adopted. Within the candidate bin width search interval, nodes are extracted with a step size of 5% of the total interval length to calculate the corresponding penalized Tsallis entropy value and find the node containing the initial maximum value. The interval length is then expanded to the left and right by 5% from this node as the fine search range. By finding the extreme point that maximizes the penalized Tsallis entropy value, this extreme point is extracted as the optimal bin width, and the maximum penalized Tsallis entropy corresponding to this extreme point is recorded.

[0043] In an optional embodiment, for any candidate bin width within the candidate bin width search interval, the instantaneous energy sequence is histogram-normalized to obtain a discrete probability distribution. The Tsallis entropy value under the current candidate bin width is calculated using a non-extensive entropy parameter, and a model complexity penalty term proportional to the number of data bins is subtracted to obtain a penalized Tsallis entropy value. This includes: dividing the data coverage of the instantaneous energy sequence into several non-overlapping data bins with their beginning and end connected according to the candidate bin width; counting the number of sample points falling in each data bin and dividing it by the total sequence length to obtain a discrete probability distribution, and obtaining the total number of data bins; calculating the Tsallis entropy value based on the non-extensive entropy parameter and the discrete probability distribution, and subtracting the model complexity penalty term from the Tsallis entropy value to obtain the penalized Tsallis entropy value. The model complexity penalty term is equal to a preset constant regularization factor multiplied by the total number of data bins.

[0044] In practical implementation, when evaluating the fitness of a candidate bin width within the candidate bin width search interval, for example, when selecting a candidate bin width of 0.025, it is necessary to construct the discrete probability mass distribution function of the instantaneous energy sequence at that resolution. First, find the maximum value in the instantaneous energy sequence, for example... And the minimum value, for example The total span of the data is determined to be 1.2. Using the current candidate bin width of 0.025 as the step size, the total span of 1.2 is divided from bottom to top into data bins that are connected end-to-end and do not overlap. The total number of data bins generated at this point is equal to the span divided by the bin width and rounded up. Next, the entire instantaneous energy sequence of length 4094 is traversed, and each sampling point is accumulated into its corresponding bin using a hash mapping, counting the frequency of samples falling into each of the i data bins. The frequency of each bin is divided by the total sequence length 4094 to obtain the discrete probability value corresponding to each data bin. This generates a normalized set of discrete probability distributions. The set satisfies the constraint that the sum of probabilities is 1.

[0045] Based on this set of discrete probability distributions, the non-zero discrete probability values ​​and the previously calculated non-extensive entropy parameter q=6.855 are substituted into the Tsallis entropy core formula for nonlinear combination operations. The computational unit performs nonlinear combination operations on all non-zero probability terms. The algorithm performs a power of 6.855 operation and sums the results. Assuming the sum is 0.012, subtracting this sum from the constant 1 yields a numerator of 0.988. Simultaneously, the denominator, q minus 1, equals 5.855. Dividing the numerator by the denominator yields the base Tsallis entropy value for this specific bin width, which is 0.1687. To reduce the algorithm's tendency to generate excessively small bin widths when seeking the maximum entropy, leading to fragmented overfitting in probability statistics, this scheme introduces a penalty mechanism. Specifically, a preset constant regularization factor, such as 0.0005, is multiplied by the total number of generated data bins (48), resulting in a model complexity penalty term of 0.024. Finally, the base Tsallis entropy value of 0.1687 is subtracted from the penalty term 0.024, resulting in the output penalized Tsallis entropy value of 0.1447. This achieves a trade-off between local feature detection and the overall generalization ability of the histogram.

[0046] In an optional embodiment, a search strategy combining coarse and fine search is employed within the candidate bin width search interval to optimize the maximum penalized Tsallis entropy and the corresponding optimal bin width. This includes: generating multiple bin width sample points proportionally in logarithmic coordinates between the lower and upper bounds of the candidate bin width search interval and calculating the corresponding penalized Tsallis entropy values ​​for each; recording the bin width sample point with the largest penalized Tsallis entropy value as the extreme value center point; if the extreme value center point is not located on the boundary, extracting the extreme value center point and the adjacent left and right bin width sample points to construct a univariate quadratic polynomial function, setting the first derivative equal to 0 to solve for the root of the independent variable as the optimal bin width, and substituting the optimal bin width into the univariate quadratic polynomial function for calculation to obtain the maximum penalized Tsallis entropy; if the extreme value center point is located on the boundary, the extreme value center point and the corresponding penalized Tsallis entropy value are taken as the optimal bin width and the maximum penalized Tsallis entropy.

[0047] In practical implementation, when searching for the optimal bin width within the set candidate bin width search interval [0.0092, 0.0434], to reduce the excessive time consumption of global exhaustive traversal, which could prevent meeting the real-time industrial response requirements, a coarse search based on a logarithmic grid is performed. The calculation module takes the logarithm to base 10 for the lower and upper bounds, with values ​​of approximately -2.036 and -1.363 respectively. Multiple interpolation points are generated at equal intervals within this logarithmic scale space; for example, in this embodiment, 10 interpolation points are preferably generated. The inverse exponential operation maps these interpolation points back to the coordinate system, obtaining 10 box width sample points that are logarithmically increasing. Since the perturbation of information resolution is more pronounced in small bin width intervals, a logarithmic distribution can achieve both small-scale dense sampling and large-scale sparse sampling. The penalized Tsallis entropy values ​​corresponding to these 10 sample points are calculated in parallel, and the maximum value is found through comparison. Assuming the 6th sample point... For example, the penalty Tsallis entropy value of 0.155 generated at a bin width of 0.022 is the highest, then... It is used as a local approximate extreme point and recorded as the extreme center point.

[0048] Then, Lagrange parabolic interpolation is initiated for a fine-grained search. The controller first checks the extreme center point. Does it fall within the boundary of the candidate bin width search interval? If the extreme value center does not fall within the boundary of the candidate bin width search interval, i.e., it is not equal to... or extreme value center point For example, a bin width of 0.022 and a corresponding entropy value of 0.155, and its left adjacent bin width sample point. For example, a bin width of 0.018 and a corresponding entropy value of 0.142, and the sample point of the right adjacent bin width. For example, the bin width is 0.027 and the corresponding entropy value is 0.151. By simultaneously solving the system of coordinates (0.018, 0.142), (0.022, 0.155), and (0.027, 0.151) for the (bin width, entropy value) of these three samples, the quadratic polynomial equation in one variable can be analytically solved. The coefficients a, b, and c in the equation are used. To obtain the pole coordinates, the stationary point formula is employed. Calculate the analytical root where the first derivative of the parabola reaches zero. For example, if the root of the independent variable is calculated to be 0.0235, this value is the optimal bin width obtained through the search. Substitute this optimal bin width of 0.0235 back into the quadratic polynomial equation to calculate the dependent variable value of 0.158 at the peak, which is the maximum penalized Tsallis entropy. It should be noted that if the extreme center point locked by the coarse search falls on the left or right endpoint of the interval, it indicates that the optimal point is at the boundary or forms a monotonic surface. In this case, the extreme center point of the boundary and the corresponding penalized Tsallis entropy value are directly used as the optimal bin width and the maximum penalized Tsallis entropy. This hierarchical joint optimization architecture only requires a small number of probability distribution operations to identify the optimal parameters.

[0049] S3, construct a three-dimensional state feature vector and determine the bearing operating state transition.

[0050] A three-dimensional state feature vector is constructed, consisting of the maximum penalty Tsallis entropy, the ratio of the optimal box width to the standard deviation of the instantaneous energy sequence, and the non-extensive entropy parameter. The squared Mahalanobis distance between the three-dimensional state feature vector and the center of the preset health state feature cluster is calculated. When the squared Mahalanobis distance is greater than the monitoring threshold, the bearing operating state is determined to change.

[0051] In practice, the standard deviation of the instantaneous energy sequence is first calculated. The optimal bin width obtained in the previous optimization step is then divided by this standard deviation to obtain a dimensionless ratio. Subsequently, the maximum penalty Tsallis entropy obtained from the optimization, this ratio, and the non-extensive entropy parameter q are concatenated sequentially to generate a three-dimensional state feature vector in a one-dimensional array format. To assess the deviation between the current state and the healthy state, a sample set of historical three-dimensional state feature vectors from pre-collected pure normal states is extracted. The mean vector of this sample set is calculated along the column direction and used as the center of the preset healthy state feature cluster. Next, the covariance matrix of the sample set is calculated, and its inverse is obtained. The center of the preset healthy state feature cluster is subtracted from the currently generated three-dimensional state feature vector to obtain the difference vector. The squared Mahalanobis distance, representing the state deviation, is obtained by calculating the continuous matrix product of the transpose, inverse, and difference vector of the difference vector. Finally, the pPF function of the chi-square distribution in the stats module of the scipy library is called to determine the monitoring threshold, and the greater than logical operator is used to determine whether the squared Mahalanobis distance is greater than the monitoring threshold. If the returned boolean value is true, an out-of-bounds warning electrical signal is output by calling the system interface, thereby determining that the bearing's operating state has changed.

[0052] In an optional embodiment, calculating the squared Mahalanobis distance between the three-dimensional state feature vector and the center of the preset health state feature cluster includes: extracting a sample set composed of multiple historical three-dimensional state feature vectors during the normal operation phase of the bearing; calculating the mean vector of the sample set as the center of the preset health state feature cluster; calculating the feature covariance matrix of the sample set and performing an inverse operation to obtain the inverse covariance matrix; subtracting the center of the preset health state feature cluster from the measured current three-dimensional state feature vector to obtain the difference vector; and calculating the continuous matrix product of the transpose of the difference vector, the inverse covariance matrix, and the difference vector to obtain the squared Mahalanobis distance.

[0053] In practical implementation, to perform unsupervised early warning of bearing failure status, a parametric model needs to be established based on baseline data of the equipment under normal operating conditions. During the initial health cycle of equipment operation, historical three-dimensional state feature vectors containing the maximum penalized Tsallis entropy, the ratio of optimal box width to instantaneous energy sequence standard deviation, and the non-extensive entropy parameter q are periodically extracted and saved. After accumulating and collecting M sets of normal feature vectors with sufficient statistics, for example, in this embodiment, 200 sets of historical three-dimensional state feature vectors are preferably accumulated and collected. The arithmetic mean of each dimension is calculated offline and combined into a column vector with a dimension of 3×1 as the center of the preset health status feature cluster. ,For example Subsequently, the variance and cross-correlation among the three features are calculated using centered matrix multiplication, generating a 3×3 symmetric feature covariance matrix. The Gaussian-Jordan elimination method is then used to invert this matrix, yielding the inverse covariance matrix resident in memory.

[0054] Once the real-time closed-loop monitoring process begins, and the system receives the raw vibration segment from the latest acquisition cycle, it immediately performs online calculations to extract and generate the three-dimensional state feature vector for the current testing cycle. Assuming the measured value at a certain moment is... To reduce the heterogeneity of units of measurement across different dimensions and eliminate redundant correlations, the state offset at a given moment is calculated using a matrix norm index. The computational unit obtains the difference vector using subtraction, which is the current 3D state feature vector minus the preset health state feature cluster center. The transpose of this difference vector is then multiplied sequentially by a pre-cached inverse covariance matrix, and finally multiplied by the difference vector itself to obtain a dimensionless scalar, which is the squared Mahalanobis distance. Let's assume the result of this matrix multiplication is 12.45.

[0055] In an optional embodiment, determining a change in bearing operating state when the squared Mahalanobis distance is greater than a monitoring threshold includes: determining a monitoring threshold based on a chi-square distribution critical value with 3 degrees of freedom at a preset confidence level; comparing the squared Mahalanobis distance with the monitoring threshold; and triggering an early warning signal when the squared Mahalanobis distance is greater than the monitoring threshold to determine that the bearing operating state has evolved from a normal state to a degraded state.

[0056] In practical implementation, based on the law of large numbers, the squared Mahalanobis distance of multiple variables in a stable controlled system follows a chi-square law. The probability distribution has a feature dimension of 3, so the system's degrees of freedom are set to 3. To determine the acceptable false alarm tolerance level, a pre-set confidence level is selected. For example, in this embodiment, the preferred pre-set confidence level is 95%, i.e., a significance level of 0.05. The chi-square distribution's ppf function in the stats module of the scipy library is called to calculate the critical value with 3 degrees of freedom and a 95% confidence level. This is obtained by inverse integration of the density function, yielding the critical right-tailed quantile constant 7.815, which is then locked as the monitoring threshold. Subsequently, the logic gate comparator uses a greater-than logical operator to compare the real-time Mahalanobis distance squared result (12.45) with the monitoring threshold 7.815. Since the former is greater than the pre-set constant limit, a Boolean value of true is returned. The system automatically outputs an out-of-bounds warning signal through the system interface, sending a warning signal to the central control room, indicating that the bearing surface has undergone health degradation changes such as pitting fatigue.

[0057] To verify the effectiveness of this scheme in practical engineering, the experiment used a full-life-cycle vibration dataset of rolling bearings from a heavy machinery test bench. The accelerometer sampling frequency was set to 12800 Hz, and the length of a single data acquisition was 4096 sampling points. Three comparison conditions were established: Scheme 1 was the complete scheme for penalized non-extensive entropy multidimensional feature extraction and Mahalanobis distance monitoring proposed in this invention; Scheme 2 was a traditional Shannon entropy monitoring scheme using fixed box width and fixed parameters; and Scheme 3 was an ablation scheme based on the scheme of this invention, removing the Teager-Kaiser energy operator. All schemes used the same normal baseline data to build the model and uniformly used the squared Mahalanobis distance value of 7.815 as the early warning monitoring threshold for operational state transitions.

[0058] During the first 100 hours of normal break-in period of the bearing operation, the average squared Mahalanobis distance calculated by Scheme 1, Scheme 2, and Scheme 3 were 2.15, 2.68, and 2.41, respectively, all remaining stably below the preset monitoring threshold. When an early minor pitting failure occurred at the 115th hour of operation, the non-extensive entropy parameter calculated by Scheme 1 rose to 5.85, corresponding to a Mahalanobis distance squared of 12.45 for the three-dimensional state characteristics, successfully exceeding the threshold of 7.815 and triggering an early warning. At the same time, the Mahalanobis distance squared of Scheme 2 was only 4.12, only reaching 8.05 and triggering a warning when the damage worsened at the 132nd hour. Although Scheme 3 reached the alarm threshold of 7.92 at the 121st hour, compared to the complete scheme of this invention, its alarm time was delayed by 6 hours, and the characteristic values ​​fluctuated significantly during the previous transition phase.

[0059] Further statistical analysis of 50 independent full-lifecycle test samples revealed that Scheme 1 achieved an average early warning time of 25.4 hours for early faults, with a false alarm rate controlled at 1.2%. Scheme 2 had an average early warning time of only 8.5 hours, but a false alarm rate as high as 6.4%. Scheme 3 had an average early warning time of 16.2 hours and a false alarm rate of 3.5%. In terms of feature sensitivity measurement at the initial stage of minor faults, Scheme 1 achieved a feature relative variation rate of 310%, significantly better than Scheme 2's 45% and Scheme 3's 120%. This demonstrates that the complete scheme of this invention can accurately identify bearing degradation characteristics under strong background noise. Figure 4 As shown in the figure, the Mahalanobis distance squared lifetime evolution curves of three comparative schemes are displayed. The data evolution trend fully verifies the timeliness and stability advantages of the present invention in early state transition warning.

[0060] Figure 2This is a schematic diagram of the Teager-Kaiser instantaneous energy sequence of the target frequency band sub-signal. The dense fluctuations in the figure represent the random noise floor of the vibration signal; the vertically upward spikes represent the transient impact characteristics generated when the bearing is damaged. The image shows that the spikes are discretely distributed on the time axis and significantly higher than the noise floor, proving that the energy operator successfully amplifies the weak impact component. Observing the amplitude changes of the pulses, it can be found that the energy corresponding to the fault impact is much greater than that of the normal fluctuations, which corresponds to the physical phenomenon described in the specific implementation method of amplifying the early pitting corrosion characteristics of the bearing using nonlinear mapping.

[0061] Figure 3 This diagram illustrates the evolution of the candidate bin width search interval with respect to the non-extensive entropy parameter. The horizontal dashed lines represent the baseline bin width calculated based on the interquartile range. The upper solid line represents the upper bound of the candidate bin width search interval, and the lower solid line represents the lower bound. The envelope between these two lines represents the candidate search space for bin width optimization. The diagram shows that as the non-extensive entropy parameter increases, the solid lines representing the upper and lower bounds gradually widen, demonstrating that the algorithm successfully achieves adaptive expansion of the bin width optimization range without incorrectly limiting the search scale to a single static interval. Observing the trend of the boundary changes, it can be found that the further the parameter deviates from the constant 1, the larger the span of the candidate bin width search interval. This corresponds to the mathematical logic described in the specific implementation method of using a weighting factor to nonlinearly expand the candidate bin width search interval to match the heavy-tailed distribution characteristics.

[0062] Figure 4 This diagram illustrates the lifetime evolution of the squared Mahalanobis distance under different monitoring schemes. The horizontal dotted lines represent the monitoring thresholds for fault detection. The solid line represents the monitoring curve for Scheme 1; the long dashed line represents the monitoring curve for Scheme 2; and the dotted line represents the monitoring curve for Scheme 3. The intersection of the curve and the monitoring threshold line represents the moment the algorithm triggers an early warning. The image shows that the solid line crosses the monitoring threshold first around 115 hours, and its rate of increase is significantly higher than other comparative curves. This demonstrates the sensitivity advantage of the multi-dimensional feature vector in early fault detection. Observing the curve fluctuations during normal operation, it can be seen that all three monitoring curves run smoothly below the dotted line without generating false alarms. This corresponds to the effect described in the specific implementation method of improving the accuracy and stability of state transition perception through three-dimensional state feature fusion.

[0063] It should be noted that those skilled in the art can make various modifications and improvements without departing from the concept of this invention, and these modifications and improvements are all within the scope of protection of this invention.

Claims

1. A method for intelligent monitoring of mechanical vibration during bearing operation, characterized in that, include: S1: Perform multi-level wavelet packet decomposition on the acquired original bearing vibration signal to obtain sub-signals of each frequency band. Extract the target frequency band sub-signals according to the kurtosis maximization principle and calculate the instantaneous energy sequence of the target frequency band sub-signals. Determine the non-extensive entropy parameter from the kurtosis coefficient of the instantaneous energy sequence through a monotonic function relationship. Using the benchmark bin width calculated based on the interquartile range of the instantaneous energy sequence as a benchmark, the benchmark bin width is expanded using a weighting factor proportional to the absolute deviation of the non-extensive entropy parameter and 1 to determine the candidate bin width search interval; S2: For any candidate bin width within the candidate bin width search interval, the instantaneous energy sequence is statistically normalized using histogram to obtain a discrete probability distribution. The Tsallis entropy value under the current candidate bin width is calculated using the non-extensive entropy parameter, and the penalized Tsallis entropy value is obtained by subtracting the model complexity penalty term proportional to the number of data bins; within the candidate bin width search interval, a search strategy combining coarse search and fine search is used to optimize the maximum penalized Tsallis entropy and the corresponding optimal bin width; S3: A three-dimensional state feature vector is constructed, consisting of the maximum penalized Tsallis entropy, the ratio of the optimal bin width to the standard deviation of the instantaneous energy sequence, and the non-extensive entropy parameter. The squared Mahalanobis distance between the three-dimensional state feature vector and the center of the preset health state feature cluster is calculated. When the squared Mahalanobis distance is greater than the monitoring threshold, the bearing operating state is determined to change.

2. The intelligent monitoring method for mechanical vibration during bearing operation according to claim 1, characterized in that, The step of performing multi-level wavelet packet decomposition on the acquired original bearing vibration signal to obtain sub-signals of each frequency band includes: using a vibration acceleration sensor to collect vibration acceleration data of the bearing as the original vibration signal; and using the db4 wavelet basis function to perform three-level wavelet packet decomposition on the original vibration signal to reconstruct the sub-signals of eight independent frequency bands.

3. The intelligent monitoring method for mechanical vibration during bearing operation according to claim 2, characterized in that, The step of extracting the target frequency band sub-signal according to the kurtosis maximization principle and calculating the instantaneous energy sequence of the target frequency band sub-signal includes: calculating the fourth-order normalized central moments of the frequency band sub-signals of the eight independent frequency bands to obtain their respective kurtosis values, and selecting the frequency band sub-signal corresponding to the maximum kurtosis value as the target frequency band sub-signal; calculating the Teager-Kaiser energy value point by point by subtracting the product of the values ​​of the previous and next adjacent sampling points of the current sampling point from the square of the current sampling point value of the target frequency band sub-signal, and constructing the instantaneous energy sequence from all the calculated energy values.

4. The intelligent monitoring method for mechanical vibration during bearing operation according to claim 1, characterized in that, The method of determining the non-extensive entropy parameter from the kurtosis coefficient of the instantaneous energy sequence through a monotonic function relationship includes: extracting the sample mean and sample standard deviation of the instantaneous energy sequence; calculating the ratio of the centered fourth moment to the fourth power of the sample standard deviation to obtain the kurtosis coefficient of the instantaneous energy sequence; substituting the kurtosis coefficient into the monotonic function relationship to calculate the corresponding non-extensive entropy parameter, wherein the monotonic function relationship is such that the non-extensive entropy parameter is equal to 0.5 times the kurtosis coefficient power of the natural constant e.

5. The intelligent monitoring method for mechanical vibration during bearing operation according to claim 1, characterized in that, The step of using the reference bin width calculated based on the interquartile range of the instantaneous energy sequence as a benchmark includes: sorting the data of the instantaneous energy sequence from smallest to largest; calculating the difference between the 75th percentile value and the 25th percentile value to obtain the interquartile range; and calculating the reference bin width using the Freedman-Diaconis rule, wherein the reference bin width is equal to a constant 2 multiplied by the interquartile range and then multiplied by the length of the instantaneous energy sequence. Power of 1.

6. The intelligent monitoring method for mechanical vibration during bearing operation according to claim 5, characterized in that, The step of expanding the baseline bin width using a weighting factor proportional to the absolute deviation of the non-extensive entropy parameter and 1 to determine the candidate bin width search interval includes: calculating the absolute value of the difference between the non-extensive entropy parameter and the constant 1; multiplying the absolute value by a preset scaling coefficient and adding an offset constant to obtain the weighting factor; determining the quotient of the baseline bin width divided by the weighting factor as the lower bound of the candidate bin width search interval; and determining the product of the baseline bin width multiplied by the weighting factor as the upper bound of the candidate bin width search interval, thus forming a closed interval as the candidate bin width search interval.

7. The intelligent monitoring method for mechanical vibration during bearing operation according to claim 1, characterized in that, The process of performing histogram-based statistical normalization on the instantaneous energy sequence for any candidate bin width within the candidate bin width search interval to obtain a discrete probability distribution, calculating the Tsallis entropy value for the current candidate bin width using a non-extensive entropy parameter, and subtracting a model complexity penalty term proportional to the number of data bins to obtain a penalized Tsallis entropy value includes: dividing the data coverage of the instantaneous energy sequence into several non-overlapping data bins with consecutive ends according to the candidate bin width; counting the number of sample points falling in each data bin and dividing by the total sequence length to obtain the discrete probability distribution and the total number of data bins; calculating the Tsallis entropy value based on the non-extensive entropy parameter and the discrete probability distribution, and subtracting the model complexity penalty term from the Tsallis entropy value to obtain the penalized Tsallis entropy value, wherein the model complexity penalty term is equal to a preset constant regularization factor multiplied by the total number of data bins.

8. The intelligent monitoring method for mechanical vibration during bearing operation according to claim 1, characterized in that, The method of using a search strategy combining coarse and fine search within the candidate bin width search interval to optimize the maximum penalty Tsallis entropy and the corresponding optimal bin width includes: generating multiple bin width sample points proportionally in logarithmic coordinates between the lower and upper bounds of the candidate bin width search interval and calculating the corresponding penalty Tsallis entropy value for each; recording the bin width sample point with the largest penalty Tsallis entropy value as the extreme value center point; if the extreme value center point is not located on the boundary, extracting the extreme value center point and the adjacent left and right bin width sample points to construct a univariate quadratic polynomial function, setting the first derivative equal to 0 to solve for the root of the independent variable as the optimal bin width, and substituting the optimal bin width into the univariate quadratic polynomial function for calculation to obtain the maximum penalty Tsallis entropy; if the extreme value center point is located on the boundary, the extreme value center point and the corresponding penalty Tsallis entropy value are taken as the optimal bin width and the maximum penalty Tsallis entropy.

9. The intelligent monitoring method for mechanical vibration during bearing operation according to claim 1, characterized in that, The calculation of the squared Mahalanobis distance between the three-dimensional state feature vector and the center of the preset health state feature cluster includes: extracting a sample set composed of multiple historical three-dimensional state feature vectors during the normal operation phase of the bearing; calculating the mean vector of the sample set as the center of the preset health state feature cluster; calculating the feature covariance matrix of the sample set and inverting it to obtain the inverse covariance matrix; subtracting the center of the preset health state feature cluster from the measured current three-dimensional state feature vector to obtain the difference vector; and calculating the continuous matrix product of the transpose of the difference vector, the inverse covariance matrix, and the difference vector to obtain the squared Mahalanobis distance.

10. The intelligent monitoring method for mechanical vibration during bearing operation according to claim 9, characterized in that, The method of determining a change in bearing operating state when the squared Mahalanobis distance is greater than the monitoring threshold includes: determining the monitoring threshold based on the chi-square distribution critical value with 3 degrees of freedom at a preset confidence level; comparing the squared Mahalanobis distance with the monitoring threshold; and triggering an early warning signal when the squared Mahalanobis distance is greater than the monitoring threshold to determine that the bearing operating state has changed from a normal state to a degraded state.

Images