Rotary machinery rub-impact fault identification method

Abstract

The invention provides a rubbing fault identification method for a rotating machine. The rubbing fault identification method comprises the following steps: acquiring a vibration signal of the rotating machine; a Hankel matrix of the vibration signals is constructed; calculating a covariance matrix of the Hankel matrix; calculating a mean value of each column in the covariance matrix and constructing a reconstruction signal R by using the mean value; performing binary quantization on the reconstructed signal R according to a 1D-LBP method by taking a local mean value of the in-window signal as a criterion to obtain a local texture signal; and carrying out spectrum analysis on the local texture signal and extracting a prominent frequency component in a spectrum to carry out rotating machinery rub-impact fault feature extraction and fault identification. According to the rotating machine rub-impact fault identification method, noise can be effectively suppressed, and effective extraction of rub-impact fault feature information and accurate fault identification can be realized.

Description

Technical Field

[0001] The invention belongs to the technical field of fault diagnosis, and in particular relates to a method for identifying rubbing faults of rotating machinery. Background Art

[0002] Rotating machinery is widely used in fields such as aviation and machinery. Effective monitoring of its operating status and accurate diagnosis of faults are crucial for minimizing economic losses and preventing major safety incidents. To improve the performance of large equipment such as aircraft engines, it is often necessary to reduce the rotor-stationary clearance, which significantly increases the likelihood of rotor-stationary rubbing. Rubbing can cause excessive vibration in the equipment and even lead to catastrophic accidents. Therefore, research on rubbing fault diagnosis methods is of great practical significance.

[0003] To accurately extract rub-impact fault signature information, researchers have used multiple approaches to identify and assess rub-impact faults, including signal decomposition (such as variational mode decomposition and empirical mode decomposition), noise reduction algorithms (such as wavelet threshold denoising and singular value denoising), and information fusion (such as principal component analysis). These signal analysis methods offer their own advantages in processing rub-impact fault signals, but they also have limitations. For example, the number of decomposition layers and the determination of the center frequency in variational mode decomposition are crucial for accurately extracting fault signatures. The recent rise of local binary patterns (LBP) and its extended one-dimensional local binary pattern (1D-LBP) method have provided new insights into fault signature extraction. The 1D-LBP algorithm can exploit texture feature extraction to mine information about equipment operating status. However, current research based on the 1D-LBP algorithm primarily focuses on directly classifying faults using the acquired vibration signals. However, the 1D-LBP algorithm itself is sensitive to noise; even small changes in neighboring elements can lead to significant differences in local structure. However, in actual engineering, there is inevitably a lot of noise, which greatly increases the difficulty of texture feature extraction based on the 1D-LBP algorithm.

[0004] Therefore, how to improve the noise resistance of 1D-LBP, achieve accurate signal quantization, and accurately extract fault information from multiple angles has become an urgent problem to be solved. Summary of the Invention

[0005] In view of this, the present invention provides a method for identifying rubbing faults in rotating machinery to solve the problems existing in the prior art.

[0006] The present invention provides a method for identifying a rubbing fault in a rotating machine, comprising:

[0007] S1: Get the vibration signal U of the rotating machinery, U=(u1,u2,u3,…u N ), where N represents the length of the vibration signal;

[0008] S2: Use the vibration signal U to construct the Hankel matrix H of the vibration signal U with a delay step of 1 between each row vector. m×n :

[0009]

[0010]

[0011] In the formula, u (i,j) = u i+j-1 , where i = 1, 2, 3..., m; j = 1, 2, 3..., n; N = m + n - 1; n < N, and m and n are the number of rows and columns of the Hankel matrix H m×n respectively;

[0012] S3: Calculate the covariance matrix C m×n of the Hankel matrix H n×n ;

[0013] S4: Calculate the mean of each column in the covariance matrix C n×n and construct the reconstructed signal R using the mean, where R = [R(1), R(2),..., R(J),..., R(n)], J = 1, 2, 3..., n, and R(J) is the mean of the J-th column in C n×n ;

[0014] S5: Based on the local mean of the in-window signal as the criterion, perform binary quantization on the reconstructed signal R according to the 1D-LBP method to obtain the local texture signal LTS, LTS = [LTS(​​​​​​​​​​​​​​​​​​​​​​

[0018]

[0019] The vector The formula is as follows:

[0020]

[0021] Where u (i,j) is the Hankel matrix H m×n The data in row i and column j of the , i = 1, 2, 3…, m; j = 1, 2, 3…, n;

[0022] S32: For the Hankel matrix H m×n Each element of is decentralized, and the formula is as follows:

[0023]

[0024] Where u′ (i,j) is the Hankel matrix H m×n The element u in row i and column j of (i,j) Decentralized data;

[0025] S33: Construct the covariance matrix C n×n ,in,

[0026]

[0027] The covariance matrix C n×n The data c in row I and column J (I,J) The calculation formula is as follows:

[0028]

[0029] Where u′ (i,I) is the Hankel matrix H m×n The decentralized data in row i and column i, u′ (i,J) is the Hankel matrix H m×n The decentralized data in the i-th row and J-th column of , where I = 1, 2, …, n and J = 1, 2, …, n.

[0030] Further preferably, in S5, the calculation method of the local texture signal LTS is as follows:

[0031] S51: Calculate the reconstructed signal R in the moving window signal r (p,q) The mean r p And the reconstructed signal R is the signal r in the moving window (p,q) and the signal r in the moving window (p,q) The mean r pThe difference b (p,q) , where the signal r in the moving window (p,q) The mean r p The calculation formula is as follows:

[0032]

[0033] Difference b (p,q) The calculation formula is as follows:

[0034] b (p,q) =r (p,q) -r p ;

[0035] Where r (p,q) is the qth vibration signal in the pth moving window, k is the window length of the moving window, k=2, 3, ... 256, q=1, 2, ... k, p=1, 2, ... n-k+1, n is the length of the reconstructed signal R;

[0036] S52: The difference b is calculated according to the 1D-LBP method. (p,q) Perform binary quantization to obtain the binary signal y(b (p,q) ), where b (p,q) >0,y(b (p,q) ) is equal to 1; when b (p,q) ≤0,y(b (p,q) ) is equal to 0;

[0037] S53: The binarized signal y(b (p,q) ) is converted back into a decimal sequence to obtain a local texture signal LTS, LTS = [LTS(1), LTS(2), ...LTS(p), ...LTS(n-k+1)], Where p = 1, 2,…n-k+1.

[0038] The rotating machinery rub-impact fault identification method provided by the present invention constructs a Hankel matrix of a signal before extracting texture features of the rub-impact fault according to the 1D-LBP method; enhances rub-impact fault information and reduces noise according to the covariance matrix of the Hankel matrix; and simultaneously, when screening effective component signals of the covariance matrix, no longer relies on signal evaluation indicators in conventional methods, but instead constructs a new signal sequence based on the mean of each column of signals in the covariance matrix; and by combining the newly obtained signal sequence with the 1D-LBP algorithm, feature extraction and fault identification of the rub-impact fault are achieved from the perspective of texture feature extraction.

[0039] The rotating machinery rub fault identification method provided by the present invention integrates a one-dimensional local binary pattern and a Hankel matrix covariance matrix. By establishing the covariance matrix of the Hankel matrix, the one-dimensional time series is expanded to a high-dimensional space, which can reduce noise while highlighting the periodic characteristics; by constructing a new signal sequence using the average value of the elements in each column of the covariance matrix, the uncertainty problem of screening effective component signals according to signal evaluation indicators can be solved; by combining 1D-LBP with the covariance matrix of the Hankel matrix, when quantizing the signal, the traditional quantization strategy of "overall vibration signal" is converted into a method of quantizing the signal Hankel matrix covariance matrix, which can solve the problem that 1D-LBP is easily affected by noise and further highlight weak fault information; by introducing 1D-LBP into rub fault identification, rub fault information can be extracted from the perspective of local texture feature extraction, and rub faults can be accurately identified based on the 1D-LBP method. BRIEF DESCRIPTION OF THE DRAWINGS

[0040] In order to more clearly illustrate the technical solutions in the embodiments of the present invention, the following briefly introduces the drawings required for use in the description of the embodiments. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without creative work.

[0041] Figure 1 A flowchart of a method for identifying rubbing faults in rotating machinery provided by the present invention;

[0042] Figure 2 This is a corresponding signal diagram in the process of identifying rubbing faults using the method provided by the present invention. In the figure, (a) and (b) are the original vibration signal U and its spectrum, respectively; (c) and (d) are the time domain and spectrum of the reconstructed signal R obtained based on the mean of the covariance matrix, respectively; (e) and (f) are the time domain and spectrum of the signal LTS obtained after binary quantization of the reconstructed signal R according to the 1D-LBP method, respectively; (g) is a local enlarged view of (f);

[0043] Figure 3 The uncertainty results of the classic method of determining the effective components of the covariance matrix by parameter index are shown in Figure 1. (a) and (b) are the uncertainty results of the method of determining the effective components of the covariance matrix by parameter index. Figure 2 The original vibration signal in (a) is obtained by screening the effective components of the covariance matrix based on information entropy and performing 1D-LBP based on the reconstructed signal, and the signal time domain and its spectrum are obtained. (c) and (d) are respectively obtained by screening the effective components of the covariance matrix based on the waveform factor and performing 1D-LBP based on the reconstruction, and (e) is a local enlarged view of (d). DETAILED DESCRIPTION

[0044] To make the objectives, technical solutions, and advantages of the present invention clearer, the present invention will be described in detail below with reference to the accompanying drawings and specific embodiments. Here, it should be noted that in order to avoid obscuring the present invention due to unnecessary details, only the processing steps closely related to the solution of the present invention are shown in the drawings, while other details less relevant to the present invention are omitted.

[0045] The present invention provides a method for identifying the rubbing fault of a rotating machine, comprising the following steps:

[0046] S1: Obtain the vibration signal U of the rotating machine, U = (u 1, u 2, u 3, …u N ), where N represents the length of the vibration signal;

[0047] S2: Use the vibration signal U to construct the Hankel matrix H of the vibration signal U with a delay step of 1 between each row vector m×n :

[0048]

[0049] In the formula, u (i,j) = u i+j-1 , where i = 1, 2, 3…, m; j = 1, 2, 3…, n; N = m + n - 1; n < N, and m and n are the number of rows and columns of the Hankel matrix H m×n respectively;

[0050] S3: Calculate the covariance matrix C m×n of the Hankel matrix H n×n ;

[0051] S4: Calculate the mean value of each column in the covariance matrix C n×n and construct the reconstructed signal R, R = [R(1), R(2), …R(J), …R(n)], J = 1, 2, 3…, n, and R(J) is equal to the mean value of the Jth column in C n×n . The calculation formula of R(J) is as follows:

[0052]

[0053] In the formula, c (In,J) represents the data in the Ith row and Jth column of the covariance matrix C n×n ;

[0054] S5: Taking the local mean of the signal within the window as a criterion, the reconstructed signal R is binary quantized according to the 1D-LBP method to obtain a local texture signal LTS, LTS = [LTS(1), LTS(2), ...LTS(p), ...LTS(n-k+1)], where k is the window length of the moving window, k = 2, 3, ...256, p = 1, 2, ...n-k+1, and n is the length of the reconstructed signal R;

[0055] S6: Perform spectrum analysis on the local texture signal LTS and extract prominent frequency components in the spectrum to perform rotating machinery rubbing fault feature extraction and fault identification.

[0056] Among them, the Hankel matrix H is calculated in S3 m×n The covariance matrix C n×n The specific steps are as follows:

[0057] S31: Calculate the Hankel matrix H m×n The mean of each column and constructing a vector using the mean of each column Among them, the Hankel matrix H m×n The mean of the jth column The calculation formula is as follows:

[0058]

[0059] The vector The formula is as follows:

[0060]

[0061] Where u (i,j) is the Hankel matrix H m×n The data in row i and column j, i = 1, 2, 3…, m; j = 1, 2, 3…, n);

[0062] S32: For the Hankel matrix H m×n Each element of is decentralized, and the formula is as follows:

[0063]

[0064] Where u′ (i,j) is the Hankel matrix H m×n The element u in row i and column j of (i,j) Decentralized data;

[0065] S33: Construct the covariance matrix C n×n ,in,

[0066]

[0067] The covariance matrix C n×n The data c in row I and column J (I,J) The calculation formula is as follows:

[0068]

[0069] Where u′ (i,I) is the Hankel matrix H m×n The decentralized data in row i and column i, u′ (i,J) is the Hankel matrix H m×n The decentralized data in the i-th row and J-th column of , where I = 1, 2, …, n and J = 1, 2, …, n.

[0070] Among them, in S5, the calculation method of the local texture signal LTS is as follows:

[0071] S51: Calculate the reconstructed signal R in the moving window signal r (p,q) The mean r p And the reconstructed signal R is the signal r in the moving window (p,q) and the signal r in the moving window (p,q) The mean r p The difference b (p,q) , where the signal r in the moving window (p,q) The mean r p The calculation formula is as follows:

[0072]

[0073] Difference b (p,q) The calculation formula is as follows:

[0074] b (p,q) =r (p,q) -r p ;

[0075] Where r (p,q) is the qth vibration signal in the pth moving window, k is the window length of the moving window, k=2, 3, ... 256, q=1, 2, ... k, p=1, 2, ... n-k+1, n is the length of the reconstructed signal R;

[0076] S52: The difference b is calculated according to the 1D-LBP method. (p,q) Perform binary quantization to obtain the binary signal y(b (p,q) ), where b (p,q) >0,y(b (p,q) ) is equal to 1; when b (p,q) ≤0,y(b (p,q) ) is equal to 0;

[0077]

[0078] S53: The binarized signal y(b (p,q) ) is converted back into a decimal sequence to obtain a local texture signal LTS to characterize the local feature information of the vibration signal, LTS = [LTS(1), LTS(2), ...LTS(p), ...LTS(n-k+1)], Where p = 1, 2,…n-k+1.

[0079] Example 1:

[0080] In order to demonstrate the effectiveness of the proposed method for identifying rubbing faults in rotating machinery, the rubbing faults in Table 1 are randomly analyzed. The rubbing position of this group of data is vertically above the turbine casing, the acceleration sensor is installed vertically below the turbine casing, the collision degree is moderate, the rotation speed is 1515.96 r / min, and the corresponding rotation frequency is 25.27 Hz (f r =1515.96 / 60=25.27).

[0081] Table 1 Test data information table

[0082]

[0083] The vibration signal of the rotating machinery in Case 1 is obtained and processed using the method proposed in the present invention, wherein in S5 , the moving window length k=8.

[0084] Among them, Figure 2 As shown in the figure, (a) and (b) are the original vibration signal U and its spectrum, respectively; (c) and (d) are the time domain and spectrum of the reconstructed signal R obtained based on the mean of the covariance matrix, respectively; (e) and (f) are the time domain and spectrum of the signal LTS obtained after binary quantization of the reconstructed signal R according to the 1D-LBP method, respectively; (g) is a local enlarged view of (f), wherein the spectrum of the signal in this embodiment is normalized.

[0085] analyze Figure 2 You can find:

[0086] (1) Spectrum of the original vibration signal ( Figure 2 (b)) High-frequency signature information is very prominent, while low-frequency signature information is relatively weak. In engineering practice, when monitoring the operating status of equipment such as aircraft engines, identification and judgment are usually based on the low-frequency signature of the signal. Therefore, monitoring rubbing faults directly based on the signal spectrum is likely to result in missed diagnosis.

[0087] (2) The spectrum of the signal constructed based on the mean of the covariance matrix ( Figure 2 (d)), the noise interference is reduced; at the same time, the relative amplitude of the low-frequency component is greatly improved; but only 2f r (51.27Hz) This frequency component makes it difficult to accurately determine the fault type;

[0088] (3) The spectrum of the signal obtained after binary quantization of the reconstructed signal R according to the 1D-LBP method ( Figure 2 (g)), the low-frequency feature information is further enriched, and the frequency components related to the rubbing fault are extracted, including 2f r (51.27Hz),5.5f r (136.7Hz),10f r (249Hz),14f r (349.1Hz) and 19f r (481Hz), based on the integer multiples and fractional multiples of these rotation frequencies, it can be determined that a friction fault has occurred in the equipment.

[0089] From (1) to (3), it can be found that the method of constructing signals based on the mean of the covariance matrix can effectively suppress noise while highlighting the low-frequency characteristics of the signal; combining it with 1D-LBP can further enrich the fault feature information in the low-frequency band of the spectrum, realize the accurate identification of rubbing faults, and is beneficial to practical engineering applications.

[0090] Analysis shows that the main reason is that the signal sequence constructed based on the mean of the column vector of the covariance matrix reduces the influence of noise and highlights the periodic characteristics, while avoiding the randomness problem of selecting component signals based on signal evaluation indicators; at the same time, the 1D-LBP method can better characterize the local texture characteristics of the signal and realize the accurate extraction of equipment operation status information.

[0091] Example 2:

[0092] In order to further demonstrate the effectiveness and superiority of the method proposed in this invention, the friction fault in Table 1 is taken as an example for analysis. Figure 3 The uncertainty result diagram of the classic method of determining the effective components of the covariance matrix by parameter indicators is shown. Among them, the characteristic parameter indicators select the more representative information entropy and waveform factor. Only when the component signal is selected, it is performed according to the signal evaluation index (the mean of the column vector is taken in the present invention). The rest of the steps are completely consistent with the method proposed in the present invention. For comparative verification and analysis, the friction fault in Table 1 of Example 1 is still selected as an example for analysis. The results are as follows Figure 3 shown. Figure 3 In the figure, (a) and (b) are the Figure 2The original vibration signal in (a) is obtained by screening the effective components of the covariance matrix based on information entropy and performing 1D-LBP based on the reconstructed signal, and the signal time domain and its spectrum are obtained. (c) and (d) are respectively obtained by screening the effective components of the covariance matrix based on the waveform factor and performing 1D-LBP based on the reconstruction, and (e) is a local enlarged view of (d).

[0093] analyze Figure 3 You can find:

[0094] (1) When information entropy is used as a parameter indicator, the signal spectrum obtained is ( Figure 3 (b) It is difficult to extract the frequency multiples or fractional frequency multiples of the rotation frequency in the low frequency band, and it is also difficult to accurately determine the rubbing fault;

[0095] (2) When the waveform factor is used as a parameter indicator, the signal spectrum obtained is ( Figure 3 (d), (e), some feature information including 2f can be extracted in the low frequency band r (51.27Hz),3f r (75.68Hz),5f r (122.1Hz);

[0096] Combining (1) and (2), we can find that since the manifestation of fault characteristics in the signal is not certain but has a certain degree of randomness, there will be great uncertainty when selecting component signals based on signal evaluation indicators that reflect certain characteristic information of the fault signal.

[0097] contrast Figure 2 (g) and Figure 3 (b) and (e) show that the method of constructing signals based on the mean of the covariance matrix proposed in the present invention can reflect the characteristics of the average level of the signal. Compared with the method of selecting effective components for reconstruction based on signal evaluation indicators such as information entropy and waveform factor, it solves the uncertainty problem existing in selecting sensitive fault component signals based on signal evaluation indicators, and can extract richer low-frequency information.

[0098] The method proposed in the present invention can effectively suppress noise and highlight the low-frequency characteristics of the signal; reduce the uncertainty of component selection; extract rub-impact fault information from the perspective of texture feature extraction, and realize accurate identification of rub-impact faults.

[0099] It should be noted that the purpose of disclosing the embodiments is to facilitate a further understanding of the present invention. However, those skilled in the art will appreciate that various substitutions and modifications are possible without departing from the spirit and scope of the present invention and the appended claims. Therefore, the present invention should not be limited to the contents disclosed in the embodiments, and the scope of protection claimed by the present invention shall be determined by the scope defined in the claims.

Claims

1. A method for identifying rubbing faults in rotating machinery, characterized in that: include: S1: Obtain the vibration signal U of the rotating machinery, U=(u 1, u 2, u 3, …u N ), where N represents the length of the vibration signal; S2: Using the vibration signal U, with 1 as the delay step between each row vector, construct the Hankel matrix H of the vibration signal U m×n : where \(u\) (i,j) = \(u\) i+j-1 , where \(i = 1, 2, 3, \cdots, m\); \(j = 1, 2, 3, \cdots, n\); \(N = m + n - 1\); \(n \lt N\), and \(m\) and \(n\) are the number of rows and columns of the Hankel matrix \(H\) m×n respectively. S3: Calculate the Hankel matrix H m×n The covariance matrix C n×n ; S4: Calculate the covariance matrix C n×n The mean of each column in C is used to construct the reconstructed signal R, where R = [R(1), R(2), ... R(J), ... R(n)], J = 1, 2, 3 ..., n, R(J) is C n×n The mean of the Jth column in ; S5: Taking the local mean of the signal within the window as a criterion, the reconstructed signal R is binary quantized according to the 1D-LBP method to obtain a local texture signal LTS, LTS = [LTS(1), LTS(2), ...LTS(p), ...LTS(n-k+1)], where k is the window length of the moving window, k = 2, 3, ...256, p = 1, 2, ...n-k+1, and n is the length of the reconstructed signal R; S6: Perform spectrum analysis on the local texture signal LTS and extract prominent frequency components in the spectrum to perform rotating machinery rubbing fault feature extraction and fault identification.

2. The rotating machinery rubbing fault identification method according to claim 1, characterized in that: S3 calculates the Hankel matrix H m×n The covariance matrix C n×n The steps are as follows: S31: Calculate the Hankel matrix H m×n The mean of each column and constructing a vector using the mean of each column Among them, the Hankel matrix H m×n The mean of the jth column The calculation formula is as follows: The vector The formula is as follows: Where u (i,j) is the Hankel matrix H m×n The data in row i and column j of the , i = 1, 2, 3…, m; j = 1, 2, 3…, n; S32: For the Hankel matrix H m×n Each element of is decentralized, and the formula is as follows: Where u′ (i,j) is the Hankel matrix H m×n The element u in row i and column j of (i,j) Decentralized data; S33: Construct the covariance matrix C n×n ,in, The covariance matrix C n×n The data c in row I and column J (I,J) The calculation formula is as follows: Where u′ (i,I) is the Hankel matrix H m×n The decentralized data in row i and column i, u′ (i,J) is the Hankel matrix H m×n The decentralized data in the j-th row and J-th column of , where I = 1, 2, …, n and J = 1, 2, …, n.

3. The rotating machinery rubbing fault identification method according to claim 1, characterized in that: In S5, the calculation method of the local texture signal LTS is as follows: S51: Calculate the reconstructed signal R in the moving window signal r (p,q) The mean r p And the reconstructed signal R is the signal r in the moving window (p,q) and the signal r in the moving window (p,q) The mean r p The difference b (p,q) , where the signal r in the moving window (p,q) The mean r p The calculation formula is as follows: Difference b (p,q) The calculation formula is as follows: b (p,q) ＝r (p,q) -r p ； Where r (p,q) is the qth vibration signal in the pth moving window, k is the window length of the moving window, k=2, 3, ... 256, q=1, 2, ... k, p=1, 2, ... n-k+1, n is the length of the reconstructed signal R; S52: The difference b is calculated according to the 1D-LBP method. (p,q) Perform binary quantization to obtain the binary signal y(b (p,q) ), where b (p,q) >0,y(b (p,q) ) is equal to 1; when b (p,q) ≤0,y(b (p,q) ) is equal to 0; S53: The binarized signal y(b (p,q) ) is converted back into a decimal sequence to obtain the local texture signal LTS, Where p = 1, 2,…n-k+1.

Images