A Strong Robust Equipment Operating State Monitoring Method Based on Multivariate Shrinkage Control Chart

Abstract

The present invention discloses a strong robust device operation state monitoring method based on a multivariate shrinkage control chart, including: First, group the collected vibration signals at equal time intervals in chronological order; use the grouped signals to construct a Hankel matrix; screen the singular value features by incremental singular entropy; calculate the L1 median instead of the mean, and calculate the H statistic of the multivariate uniformly weighted moving control chart; adopt a covariance shrinkage operator based on M-estimation to shrink the covariance of the H statistic; construct a T 2 statistic and calculate the control line by using the kernel density estimation method; construct a control chart, and judge the monitoring state by comparing the statistic with the control line. If the statistic exceeds the control line, it is considered that the device has a fault, so as to achieve the purpose of device state monitoring. The bearing is a typical component in a rotating machine. In this embodiment of the present invention, the bearing is taken as the object, and a set of publicly available vibration data of the bearing full-life test is used to verify the usability and generality of the present invention.

Description

Technical Field

[0001] The present invention belongs to the field of condition monitoring of rotating machinery equipment, and particularly relates to a strong-robust signal construction and equipment operation state monitoring method. Background Art

[0002] In the field of industrial equipment, rotating machinery, as the main body or a key component of industrial machinery, its stable and reliable operation is the guarantee for the safe operation of the entire equipment. However, rotating machinery is extremely vulnerable to damage and defects. Once a failure occurs during the working process, it is very likely to have a significant impact on the operation of the overall production equipment, resulting in huge economic losses or even major accidents. Therefore, it is of great significance to monitor the operation state of rotating machinery equipment and carry out early fault warning.

[0003] For the condition monitoring of rotating machinery, common monitoring methods include vibration analysis method, acoustic emission method, temperature analysis method, etc. Among them, since the vibration signals have great differences in different parts and different degrees of fault manifestations, and the vibration signals are relatively intuitive and have clear physical meanings, the vibration analysis method is currently a relatively common and commonly used monitoring method.

[0004] The feature extraction of signals is a key step in equipment condition monitoring, and the selection of features plays a crucial role in monitoring the degradation process of equipment. Traditional feature indicators include time-domain features, frequency-domain features, time-frequency domain features, energy features, etc. Considering that in rotating machinery, typical components such as gears and bearings are difficult to identify early fault signals in a large environmental noise, traditional feature indicators cannot well display the degradation state of equipment performance. This patent proposes a feature extraction method based on the Hankel matrix, and uses incremental singular entropy to select appropriate singular values as monitoring features.

[0005] Multivariate statistical process control is a commonly used method for quantitatively analyzing the characteristics of target parameters using control charts, and is one of the important methods of modern quality management. This method mainly includes two steps: one is to obtain the control line according to the healthy process data; the other is to monitor the subsequent process based on the obtained control line. Considering that there is a long-term small change in the occurrence of faults in rotating machinery equipment, this patent uses the MHWMA (Multivariate Homogeneous Weighted Moving) control chart to monitor statistical indicators. The MHWMA control chart takes into account historical data and can effectively monitor small mean drifts in the process, but lacks effective planning for healthy data and has more false alarms. This patent uses a covariance shrinkage operator based on M-estimation to shrink the data covariance, improve the robustness of monitoring features, and uses KDE (Kernel Density Estimation) to estimate the control line, reduce the false alarm rate, and improve the monitoring accuracy. Summary of the Invention

[0006] In view of the above, the present invention provides a strongly robust method for monitoring the operating state of equipment based on a multivariate shrinkage control chart. The object of the present invention is to extract degradation characteristics and monitor the state of rotating mechanical equipment.

[0007] To achieve the above object, the present invention provides the following technical means:

[0008] A strongly robust method for monitoring the operating state of equipment based on a multivariate shrinkage control chart, comprising:

[0009] S1. First, group the collected vibration signals at equal time intervals in chronological order;

[0010] S2. Construct a Hankel matrix sequence S from the obtained signal x(t) according to the sample sequence with a step size of L = Fs / r i , where Fs is the sampling frequency and r is the rotational speed of the rotating shaft.

[0011]

[0012] where N = L and i is the number of Hankel matrix sequences.

[0013] S3. Perform SVD (singular value decomposition) on the obtained Hankel matrix sequence S i to obtain the corresponding eigen-sequence S n (i, j). Where i is the number of eigenvectors obtained from the Hankel matrix decomposition and j is the Hankel matrix sequence number.

[0014] S4. Calculate the incremental singular entropy of each sequence of S n (i), and take the singular value sequence with an incremental singular entropy greater than 0 and increasing and an increment greater than the singular value sequence of the next lower order of the maximum increment as the monitoring feature S s (i, j). Where i is the number of features and j is the number of Hankel sequences.

[0015] S5. Select the monitoring features under normal conditions to calculate the L1 median u of S s (i), and calculate the multivariate uniform weighted moving control chart statistic H and the covariance sequence σ of all samples MM . H .

[0016]

[0017] H i = wS s (:, i)+(1 - w)μ MM #(3)

[0018]

[0019] where w is the smoothing coefficient and σ0 is the monitoring feature S s the covariance of (i, j); μ MM is S m the final value of, where S m is the quantity set for obtaining μ MM ;

[0020] S6. Solve the covariance shrinkage weight β based on M-estimation and calculate the shrunk covariance ∑ sr . Where p is the number of feature variables of S s (i), I is the unit diagonal matrix, where tr(σ H ) is the trace of the feature covariance sequence and p is the feature dimension;

[0021] S7. Calculate the T 2 statistic according to the above obtained results.

[0022]

[0023] where u MM is the L1 median and H is the multivariate homogeneous weighted moving average (MHWMA) control chart statistic.

[0024] S8. Estimate the distribution of the normal data T 2 using kernel density estimation to obtain the control line.

[0025] S9. Plot the T 2 graph and analyze it according to the judgment criterion using the obtained control line to judge the operating state of the device.

[0026] As a further improvement of the present invention, the step S4 specifically includes:

[0027] Take the average of the healthy part of the feature sequence S n (i, j) obtained in step S3 to obtain the feature sequence S n2 (i), i = 1, 2,..., L / 2. Calculate its incremental singular entropy according to the following formula:

[0028]

[0029] Curvature represents the sensitivity to state changes and can be used to reflect the state changes of the incremental singular entropy sequence. Select appropriate singular value feature numbers. Since it is a discrete sequence here, use the difference to approximate the derivative:

[0030]

[0031] y″ = ΔE i+1 -2ΔE i +ΔE i-1 #(8)

[0032] y′ = min(ΔE i+1 -ΔE i , ΔE i -ΔE i-1 ) #(9)

[0033] Calculate the curvature increment ΔCE j = C j -C j+1 , take ΔCE j > 0 and take the first term of the curvature increment as the cut-off with a magnitude one order smaller to obtain the monitoring feature sequence S s (i, j).

[0034] As a further improvement of the present invention, the step S6 specifically includes:

[0035] The calculation formula of the contraction operator β is as follows:

[0036]

[0037]

[0038]

[0039]

[0040] Therefore, the required contraction covariance is:

[0041]

[0042] where σ H is the covariance of the feature sequence S s (i, j), p is the feature dimension, and I is the unit diagonal matrix.

[0043] The beneficial effect of the present invention is that a strong robust device operation state monitoring method based on a multivariate contraction control chart can effectively extract the performance degradation index of rotating mechanical equipment and complete the function of state monitoring, and the monitoring effect is good, and it can detect early faults of the equipment. It is a method that can be applied to industry BRIEF DESCRIPTION OF THE DRAWINGS

[0044] In order to more clearly illustrate the existing solutions and the technical solutions in the embodiments of the present invention, the following will give the relevant drawings for the description of the existing solutions and the technical description of the embodiments and briefly introduce them. Obviously, the following described drawings are only some embodiments recorded in the present invention. For those of ordinary skill in the art, other drawings can be obtained according to the following drawings without creative labor.

[0045] Figure 1Schematic diagram of the process of a strong robust device operation status monitoring method based on a multivariate shrinkage control chart according to the present invention.

[0046] Figure 2 Complete trend chart of bearing degradation monitoring characteristic indicators in the embodiment.

[0047] Figure 3 Robust multivariate control chart and detailed display chart in the embodiment. Specific embodiments

[0048] To more clearly describe the technical solutions and advantages of the present invention, the present invention will be further described in detail below in conjunction with the accompanying drawings and embodiments. Obviously, the specific embodiments described herein are only used to explain the present invention and are not limited to this example. The present invention is implemented through the following technical solutions.

[0049] The bearing is a typical component in rotating machinery, and the example of the present invention takes this as the target object.

[0050] The data comes from the full-life test of the bearing conducted by the Intelligent Equipment Maintenance Center of the University of Cincinnati in the United States. This test records all the data of the bearing from normal operation to failure in chronological order. This example selects the data recorded in Test 2. This data contains a total of 984 data files, each file contains 20,480 data points, and the vibration signals of the bearing from the healthy state to the obvious failure are recorded every 10 minutes. The recording time is from 10:32:39 on February 12, 2004 to 6:22:39 on February 19, 2004. As shown in the appendix Figure 1 As shown, in this example, a strong robust device operation status monitoring method based on a multivariate shrinkage control chart includes the following steps:

[0051] S1. Classify the target bearing files in 984 data files by data number, and 984 sample sets can be obtained.

[0052] S2. Construct the Hankel matrix sequence S of the obtained signal x(t) according to the sample sequence with a step size of L = Fs / r, where Fs is the sampling frequency and r is the rotational speed of the rotating shaft. In the example, Fs is 20 kHz and r is 2000 RPM, then the length of the Hankel matrix window function is 600. i

[0053] S3. Perform SVD (singular value decomposition) on the obtained Hankel matrix sequence S i to obtain the corresponding eigen-sequence S n(i, j). Here, i is the number of features obtained from the Hankel matrix decomposition, and j is the serial number of the Hankel matrix sequence. Considering the large amount of data in this example, each sample has 34 Hankel matrices, and there are many singular value decomposition results. In this example, the average value of the 34 feature matrices after singular value decomposition of each sample is taken as the feature S of each sample. n (i, j).

[0054] S4. Calculate S n For each increment singular entropy of the (i) sequence, select the singular value sequence with an increment singular entropy greater than 0 and increasing, and the increment greater than the singular value sequence with the largest increment by one order of magnitude as the monitoring feature S s (i, j). Here, i is the number of features, and j is the number of Hankel sequences. Select the monitoring features under normal conditions to calculate S s The L1 median u of S(i) MM , and calculate the multivariate uniformly weighted moving control chart statistic and covariance sequence σ of all samples H .

[0055]

[0056] H i = wS s (:, i) + (1 - w)μ MM

[0057]

[0058] where w is the smoothing coefficient, taken as 0.1 in this example, and σ0 is the covariance of the monitoring feature S(i, j). s (i, j).

[0059] S6. Solve the covariance shrinkage weight β based on M-estimation, and calculate the shrunk covariance ∑ sr . Among them p is the number of feature variables of S s (i), and I is the identity diagonal matrix.

[0060] S7. According to the results obtained above, calculate the T2 statistic.

[0061]

[0062] where u MM is the L1 median of the first 400 samples, and H is the statistic of the MHWMA control chart.

[0063] Draw the trend chart of the T 2 statistic, as Figure 2 shown. Obviously, the first half of the data is relatively stable, and the first 400 samples can be used as the healthy data under normal operation of the bearing.

[0064] S8. Estimate the distribution of normal data T using kernel density estimation to obtain the control line UCL. 2

[0065] S9. Plot the T 2 image as shown in Figure 3 and analyze it according to the judgment criteria using the obtained control line to judge the operating state of the device.

[0066] As Figure 3 shown in the detailed display of the robust multivariate control chart, it is not difficult to find that almost all the bearing performance degradation indicators calculated from the data with sample numbers 1 - 531 are below the control line; almost all the bearing performance degradation indicators calculated from sample 539 onwards are above the control line; the bearing performance degradation indicators calculated from the data of samples 532 - 538 fluctuate above and below the control line. Then, based on this, it can be judged that: the bearing had a fault omen at 4:13:39 on February 16, 2004, and started to fail at 5:23:39 on February 16, 2004.

[0067] In summary, the strong robust device operating state monitoring method based on multivariate shrinkage control chart of the present invention can effectively extract the performance degradation indicators of rotating mechanical equipment and complete the function of state monitoring with good monitoring effect, and can detect early faults of the device, which is a method that can be applied to industry.

[0068] The reference terms "one embodiment", "example", "specific example", etc. described in this specification mean that the specific features, materials, structures or characteristics of the embodiment are included in at least one embodiment or example of the present invention. In this specification, the schematic meanings of the above terms do not necessarily refer to the same embodiment or example. Moreover, the specific features, structures, materials or characteristics described can be combined in a suitable manner in any one or more embodiments or examples.

[0069] Although the embodiments of the present invention have been shown and described, those of ordinary skill in the art can understand that various changes, modifications, substitutions and variations can be made to these embodiments without departing from the principles and purposes of the present invention, and the scope of the present invention is defined by the claims and their equivalents.​

Claims

1. A strong robust method for monitoring the operating state of equipment based on a multivariate shrinkage control chart, characterized in that, The method includes the following steps: S1. First, group the collected vibration signals at equal time intervals in chronological order; S2. Construct a Hankel matrix sequence S from the obtained signal x(t) with a step size L = Fs / r, where Fs is the sampling frequency and r is the rotational speed of the rotating shaft; i , where Fs is the sampling frequency and r is the rotational speed of the rotating shaft; S3. Perform SVD singular value decomposition on the obtained Hankel matrix sequence S i to obtain the corresponding eigen-sequence S n (i, j), where i is the number of eigenvectors obtained from the Hankel matrix decomposition and j is the Hankel matrix sequence number; S4. Calculate S n (i) The incremental singular entropy of each sequence. Select the singular value sequence with an incremental singular entropy greater than 0 and increasing, and the increment greater than the singular value sequence of the next lower order of magnitude of the maximum increment as the monitoring feature S s (i,j), where i is the number of features and j is the number of Hankel sequences; S5. Select the monitoring features in the normal state to calculate S s The L1 median u of (i) MM , and calculate the multivariate uniformly weighted moving control chart statistic H for all samples i and the covariance sequence σ H ; S6. Solve for the covariance shrinkage weight β based on M-estimation and calculate the shrunk covariance Σ sr ; where I is the unit diagonal matrix, where tr(σ H ) is the trace of the characteristic covariance sequence, and p is the characteristic dimension; S7. Calculate T based on the result obtained in S6 2 Statistic where u MM is the median of L1, and H is the statistic of the multivariate homogeneous weighted moving average control chart; S8. Use kernel density estimation to estimate the distribution of health data T 2 to obtain a control line; S9. Draw T 2 image, and analyze it according to the judgment criterion by using the obtained control line to judge the operation state of the device.

2. The strong robust device operating state monitoring method based on a multivariate shrinkage control chart according to claim 1, wherein, Matrix sequence S i is as follows: where N = L.

3. A strong robust device operating state monitoring method based on a multivariate shrinkage control chart according to claim 1, characterized in that, Calculate the multivariate uniformly weighted moving control chart statistic H for all samples i and the covariance sequence σ H as follows: H i = wS s (:, i)+(1 - w)μ MM (2) Among them, S s (:, i) is the above-mentioned monitoring feature sequence, and μ MM is the median of the feature sequence L1 of S s (:, i), and S m is a quantity set for obtaining μ MM ; where w is the smoothing coefficient and σ0 is the monitoring feature S s (i,j) covariance 4. A strong robust device operating state monitoring method based on a multivariate shrinkage control chart according to claim 1, characterized in that, For the normal operation characteristic sequence S obtained in step S3 n The average value of (i, j) is obtained to get the characteristic sequence S n2 (i), i = 1, 2, …, L / 2, and its incremental singular entropy is calculated according to the following formula: Curvature represents the sensitivity to state changes and can be used to reflect the state changes of the incremental singular entropy sequence. Select appropriate singular value feature numbers. Since this is a discrete sequence here, use the difference to approximate the derivative: y″ = ΔE i+1 -2ΔE i +ΔE i-1 (8) y′=min(ΔE i+1 -ΔE i ,ΔE i -ΔE i-1 ) (9) Calculate the curvature increment ΔCE j = C j - C j+1 , take ΔCE j > 0 and take the first term of the curvature increment as the cut-off with a smaller order of magnitude to obtain the monitoring feature sequence S s (i, j), where C j is the value in the required curvature C sequence.

5. A strong robust device operating state monitoring method based on a multivariate shrinkage control chart according to claim 1, characterized in that, In step S6, the calculation formula of the contraction weight β is as follows: where p is the feature dimension obtained in step S4, n is the number of columns and the length of the features obtained in S4, and γ is the spherical parameter, where: Therefore, the required contraction covariance is: where σ H is the covariance of the feature sequence S s (i, j), and I is the identity diagonal matrix.

Images