Rotating machinery fault diagnosis method based on improved continuous variational modal decomposition

Abstract

The method for diagnosing a rotating machine fault based on improved continuous variational modal decomposition includes: S1: collecting vibration signals; S2: setting the variation range and increasing step length of the balance parameter α; S3: performing continuous variational modal decomposition SVMD with each step value of α increasing; S4: calculating the energy position index value EP of the modal component obtained in each step; S5: finding the minimum EP index value of each step and the serial number of the corresponding target modal component; S6: drawing a relationship curve between α and the minimum EP index value of each step and the serial number of the target modal component; S7: finding the optimal value of α; S8: finding the serial number of the target modal component corresponding to the optimal value of α; S9: performing SVMD with the optimal value of α, selecting the optimal target modal component; S10: extracting the fault characteristic frequency, and realizing the diagnosis of the rotating machine fault. The present invention overcomes the problem that the balance parameters are difficult to determine and the target modal components are difficult to select in the existing SVMD, and can conveniently and effectively realize the diagnosis of the rotating machine fault.

Description

Technical Field

[0001] The invention belongs to the technical field of mechanical fault diagnosis, and in particular relates to a rotating machinery fault diagnosis method based on improved continuous variational modal decomposition. Background Art

[0002] Rotating machinery and equipment usually contain rotating parts such as bearings and gears, which have a vital impact on the normal operation of rotating machinery and equipment. However, in practice, these parts often have to operate under dynamic loads, alternating loads, and even overload conditions, and are prone to various types of fault damage, such as pitting, cracks, spalling, and wear, which seriously affect the performance and economic benefits of rotating machinery and equipment. Therefore, it is necessary to take effective technical measures to monitor the condition and diagnose faults of rotating machinery and equipment, so as to detect damaged parts as early as possible, and repair or replace them in time to prevent them from happening and reduce the maintenance cost of rotating machinery and equipment.

[0003] A common method for state monitoring and fault diagnosis of rotating machinery is to collect vibration signals of the equipment, process and analyze the vibration signals, identify the fault characteristics, and then diagnose the fault status of the equipment. However, since there are usually many rotating parts in rotating machinery, the vibration signals excited by different rotating parts will be coupled and superimposed with each other, and the vibration signals excited by damaged parts will usually be greatly attenuated on the transmission path from the excitation source to the vibration signal collection point. Coupled with the interference of background noise, it is not easy to identify weak fault characteristics from the vibration signals collected by the sensor. Research in this area is also an important part of the research on fault diagnosis of rotating machinery.

[0004] At present, a common method for extracting fault features from vibration signals of rotating machinery is variational mode decomposition (VMD). The VMD algorithm was proposed by scholars Dragomiretskiy and Zosso in 2014. It decomposes the vibration signal into a series of narrowband intrinsic mode components with different center frequencies. It has a complete mathematical theoretical support and good noise resistance. It has received extensive attention and research in the field of rotating machinery fault diagnosis. However, VMD also has some shortcomings, such as the difficulty in optimizing the number of modal decomposition layers and balance parameters, and the difficulty in selecting the target modal components in the decomposition results. For this reason, scholars Nazari and Sakhaei proposed the continuous variational mode decomposition (SVMD) algorithm in 2020 based on VMD. SVMD does not need to know the number of modal decomposition layers when executing, and can adaptively decompose the signal into a series of modal components. Compared with VMD, SVMD lacks the parameter of the number of modal decomposition layers, so it is easier to execute and has higher computational efficiency. But even so, SVMD still has the problems of difficult setting of balance parameters and difficult selection of target modal components in the decomposition results. These two problems have a crucial impact on the performance of SVMD, especially in the processing and analysis of rotating machinery vibration signals. Solving these two problems is the key to extracting fault features and realizing fault diagnosis. Summary of the invention

[0005] The present invention aims to overcome the shortcomings of the prior art continuous variational modal decomposition in the extraction of fault features of rotating machinery vibration signals, and proposes a rotating machinery fault diagnosis method based on improved continuous variational modal decomposition.

[0006] The rotating machinery fault diagnosis method based on improved continuous variational modal decomposition of the present invention comprises the following steps:

[0007] S1: On rotating machinery, s is the sampling frequency, and a vibration signal x(t) with a length of L is collected;

[0008] S2: Set the range of the balance parameter α to [α min ,α max ], let α be from α min Start with a step size of s α Increase, and the value when it increases to the i-th step is:

[0009] α i =α min +(i-1)·s α (I)

[0010] Where i = 1, 2, 3, ..., N s , N s is the total number of steps that α increases, and N s =(αmax -α min ) / s α +1;

[0011] S3: With α i As a balance parameter, continuous variational mode decomposition SVMD is performed to obtain multiple modal components u ik (t)(k=1,2,…,K i ), where k is the serial number of the modal component, K i is the number of modal components obtained by the continuous variational mode decomposition (SVMD) in the i-th step;

[0012] S4: Calculate the modal component u ik Energy position index value EP of (t) i (k);

[0013] S5: From K i Energy position index EP i Find the minimum value minEP in (k) i and the serial number of its corresponding target modal component i ;

[0014] S6: According to the final N s minEP i Value, draw α i With minEP i and draw the relationship curve between i With ON i The relationship curve between

[0015] S7: Using α i With minEP i The relationship curve between minEP and α is used to find out the relationship between minEP and α when the equilibrium parameter α increases. i The α value corresponding to the minimum value is used as the optimal value of the balance parameter α opt ;

[0016] S8: Using α i With ON i The relationship curve between i The value is α opt Corresponding ON i The value of ON T ;

[0017] S9: With α opt As the value of the balance parameter α, perform continuous variational mode decomposition SVMD to obtain multiple modal components, from which the ONth modal component is selected. T The modal component is the optimal target modal component u opt (t);

[0018] S10: Calculate the optimal target modal component u opt The square envelope spectrum of (t) is used to extract the fault characteristic frequency and perform fault diagnosis on rotating machinery.

[0019] Preferably, in step S4 and step S5, the energy position index value calculation of the modal component includes the following steps:

[0020] S2-1: For a certain modal component u(t), calculate its square envelope spectrum SES(f), where f represents the frequency;

[0021] S2-2: Normalize the square envelope spectrum SES(f) to obtain the normalized square envelope spectrum NS(f). The specific method is:

[0022]

[0023] Wherein, max[SES(f)] represents the maximum amplitude in the square envelope spectrum SES(f);

[0024] S2-3: Calculate the energy concentration index EC of the fault characteristic frequency component using the normalized square envelope spectrum NS(f);

[0025] S2-4: Calculate the position accuracy index PA of the fault characteristic frequency using the normalized square envelope spectrum NS(f);

[0026] S2-5: Calculate the energy position index EP corresponding to the modal component u(t), the formula is:

[0027]

[0028] Where p is the adjustment coefficient, and 0<p≤1, β is the balance coefficient, and the calculation formula of β is:

[0029]

[0030] Further preferably, in step S10 and step S2-1, the calculation of the square envelope spectrum of the modal component includes the following steps:

[0031] S3-1: For a certain modal component u(t), calculate its square envelope signal u SE (t), the formula is:

[0032]

[0033] Where j is the imaginary unit, represents Hilbert transform, |·| represents the modulus of a complex number;

[0034] S3-2: Calculate uSE The square envelope spectrum SES(f) of (t) is given by:

[0035]

[0036] in, represents Fourier transform.

[0037] Further preferably, in step S2-3 and step S2-5, the calculation of the energy concentration index EC of the fault characteristic frequency component includes the following steps:

[0038] S4-1: rearrange the amplitudes of the normalized square envelope spectrum NS(f) from large to small to obtain a rearranged amplitude sequence SNS(n) (n=1,2,…,N), where N represents the number of frequency points;

[0039] S4-2: Calculate the average value of the difference between the first amplitude and the subsequent M amplitudes in SNS(n) to obtain the energy concentration index EC, the formula is:

[0040]

[0041] Further preferably, in step S2-4 and step S2-5, the calculation of the position accuracy index PA of the fault characteristic frequency includes the following steps:

[0042] S5-1: Find the frequency value f corresponding to the maximum amplitude in the normalized square envelope spectrum NS(f) ma ;

[0043] S5-2: Calculate the position accuracy index PA, the formula is:

[0044]

[0045] Where C is the number of rotating parts in the rotating machinery, f j (j=1,2,…,C) is the theoretical fault characteristic frequency of each rotating component.

[0046] The positive effects achieved by the present invention are as follows: the present invention utilizes the continuous variational modal decomposition method to avoid the problem that the number of modal decomposition layers in the existing variational modal decomposition method is difficult to determine; the proposed energy position index is utilized to overcome the problem that the target mode is difficult to select in the existing continuous variational modal decomposition method; the improved continuous variational modal decomposition method can be used to obtain the optimal target modal component containing sufficient and complete fault feature information, and combined with the square envelope spectrum analysis, the fault feature frequency can be conveniently and effectively extracted to realize the fault diagnosis of rotating machinery. BRIEF DESCRIPTION OF THE DRAWINGS

[0047] Figure 1It is a flow chart of the implementation of the present invention;

[0048] Figure 2 A time domain waveform diagram of a rolling bearing vibration signal in an embodiment of the present invention;

[0049] Figure 3 Graph 1 is a square envelope spectrum of a rolling bearing vibration signal in an embodiment of the present invention.

[0050] Figure 4 In the embodiment of the present invention, in the continuous variational mode decomposition, the balance parameter α i and the minimum energy position index minEP i The relationship curve between , where i = 1, 2, 3, ..., 200;

[0051] Figure 5 In the embodiment of the present invention, in the continuous variational mode decomposition, the balance parameter α i The modal component number corresponding to the minimum energy position index is ON i The relationship curve between , where i = 1, 2, 3, ..., 200;

[0052] Figure 6 is a time domain waveform diagram of the optimal target modal component in an embodiment of the present invention;

[0053] Figure 7 is the square envelope spectrum of the optimal target modal component in the embodiment of the present invention. DETAILED DESCRIPTION

[0054] The present invention will be further described below in conjunction with the accompanying drawings and embodiments.

[0055] Reference Figure 1 , a rotating machinery fault diagnosis method based on improved continuous variational modal decomposition includes the following steps:

[0056] S1: On rotating machinery, s is the sampling frequency, and a vibration signal x(t) with a length of L is collected;

[0057] S2: Set the range of the balance parameter α to [α min ,α max ], let α be from α min Start with a step size of s α Increase, and the value when it increases to the i-th step is:

[0058] α i =α min +(i-1)·s α (I)

[0059] Where i = 1, 2, 3, ..., N s , Ns is the total number of steps that α increases, and N s =(α max -α min ) / s α +1;

[0060] S3: With α i As a balance parameter, continuous variational mode decomposition SVMD is performed to obtain multiple modal components u ik (t)(k=1,2,…,K i ), where k is the serial number of the modal component, K i is the number of modal components obtained by the continuous variational mode decomposition (SVMD) in the i-th step;

[0061] S4: Calculate the modal component u ik Energy position index value EP of (t) i (k);

[0062] S5: From K i Energy position index EP i Find the minimum value minEP in (k) i and the serial number of its corresponding target modal component i ;

[0063] S6: According to the final N s minEP i Value, draw α i With minEP i and draw the relationship curve between i With ON i The relationship curve between

[0064] S7: Using α i With minEP i The relationship curve between minEP and α is used to find out the relationship between minEP and α when the equilibrium parameter α increases. i The α value corresponding to the minimum value is used as the optimal value of the balance parameter α opt ;

[0065] S8: Using α i With ON i The relationship curve between i The value is α opt Corresponding ON i The value of ON T ;

[0066] S9: With α opt As the value of the balance parameter α, perform continuous variational mode decomposition SVMD to obtain multiple modal components, from which the ONth modal component is selected. TThe modal component is the optimal target modal component u opt (t);

[0067] S10: Calculate the optimal target modal component u opt The square envelope spectrum of (t) is used to extract the fault characteristic frequency and perform fault diagnosis on rotating machinery.

[0068] In step S4 and step S5, the calculation of the energy position index value of the modal component includes the following steps:

[0069] S2-1: For a certain modal component u(t), calculate its square envelope spectrum SES(f), where f represents the frequency;

[0070] S2-2: Normalize the square envelope spectrum SES(f) to obtain the normalized square envelope spectrum NS(f). The specific method is:

[0071]

[0072] Wherein, max[SES(f)] represents the maximum amplitude in the square envelope spectrum SES(f);

[0073] S2-3: Calculate the energy concentration index EC of the fault characteristic frequency component using the normalized square envelope spectrum NS(f);

[0074] S2-4: Calculate the position accuracy index PA of the fault characteristic frequency using the normalized square envelope spectrum NS(f);

[0075] S2-5: Calculate the energy position index EP corresponding to the modal component u(t), the formula is:

[0076]

[0077] Where p is the adjustment coefficient, and 0<p≤1, β is the balance coefficient, and the calculation formula of β is:

[0078]

[0079] In step S10 and step S2-1, the calculation of the square envelope spectrum of the modal component includes the following steps:

[0080] S3-1: For a certain modal component u(t), calculate its square envelope signal u SE (t), the formula is:

[0081]

[0082] Where j is the imaginary unit, represents Hilbert transform, |·| represents the modulus of a complex number;

[0083] S3-2: Calculate u SE The square envelope spectrum SES(f) of (t) is given by:

[0084]

[0085] in, represents Fourier transform.

[0086] In step S2-3 and step S2-5, the calculation of the energy concentration index EC of the fault characteristic frequency component includes the following steps:

[0087] S4-1: rearrange the amplitudes of the normalized square envelope spectrum NS(f) from large to small to obtain a rearranged amplitude sequence SNS(n) (n=1,2,…,N), where N represents the number of frequency points;

[0088] S4-2: Calculate the average value of the difference between the first amplitude and the subsequent M amplitudes in SNS(n) to obtain the energy concentration index EC, the formula is:

[0089]

[0090] In step S2-4 and step S2-5, the calculation of the position accuracy index PA of the fault characteristic frequency includes the following steps:

[0091] S5-1: Find the frequency value f corresponding to the maximum amplitude in the normalized square envelope spectrum NS(f) ma ;

[0092] S5-2: Calculate the position accuracy index PA, the formula is:

[0093]

[0094] Where C is the number of rotating parts in the rotating machinery, f j (j=1,2,…,C) is the theoretical fault characteristic frequency of each rotating component.

[0095] The above invention is applied to the vibration signal processing of rolling bearings. The bearing data comes from the Bearing Data Center of Western Reserve University. The data used is the vibration data of the outer ring of the rolling bearing. The specific fault is a single point damage with a diameter of 0.021 inches and a depth of 0.011 inches at the 6 o'clock position of the outer ring. The load applied by the loading motor is 3HP. The output speed of the drive motor measured by the encoder is 1721r / min, and the corresponding rotation frequency is f r =1721 / 60=28.6833Hz, so the theoretical fault characteristic frequency of the inner ring of the rolling bearing is calculated to be f ir=5.4152·f r =155.3260Hz, the theoretical fault characteristic frequency of the outer ring is f or =3.5848·f r =102.8240Hz, the theoretical fault characteristic frequency of the rolling element is f ba =4.7135·f r =135.1989 Hz. The sampling frequency of the rolling bearing vibration data is 12000 Hz. The steps of selecting the target mode using the present invention are as follows.

[0096] S1: Select a vibration signal x(t) with a length of L=4000 from the data set corresponding to the rolling bearing. Its time domain waveform is as follows: Figure 2 As shown, its square envelope spectrum is as follows Figure 3 As shown, it can be seen that there are many large-amplitude interference components, and the characteristic frequency component of the rolling bearing fault indicated by f=102Hz is not prominent, making it difficult to determine that the rolling bearing has suffered a fault damage;

[0097] S2: Set the range of the balance parameter α to [50,10000], starting from 50 and increasing in step size s α =50, the value when it increases to the i-th step is:

[0098] α i =50+(i-1)×50

[0099] Where i = 1, 2, 3, ..., N s , N s is the total number of steps that α increases, and N s =200;

[0100] S3: With α i As a balance parameter, continuous variational mode decomposition SVMD is performed to obtain multiple modal components u ik (t)(k=1,2,…,K i ), where k is the serial number of the modal component, K i is the number of modal components obtained by the continuous variational mode decomposition (SVMD) in the i-th step;

[0101] S4: Calculate the modal component u ik Energy position index value EP of (t) i (k), comprising the following steps:

[0102] S4-1: For modal component u ik (t), calculate its square envelope spectrum SES ik (f), where f represents frequency;

[0103] S4-2: Squared Envelope Spectrum SESik (f) Perform normalization processing to obtain the normalized square envelope spectrum NS ik (f) The specific method is:

[0104]

[0105] Among them, max[SES ik (f)] represents the square envelope spectrum SES ik (f) the maximum amplitude;

[0106] S4-3: Using the normalized square envelope spectrum NS ik (f) Calculate the energy concentration index EC of the fault characteristic frequency component ik , which in turn comprises the following steps:

[0107] S4-3-1: Normalized square envelope spectrum NS ik The amplitudes of (f) are rearranged from large to small to obtain the rearranged amplitude sequence SNS ik (n)(n=1,2,…,2000);

[0108] S4-3-2: Calculate SNS ik The energy concentration index EC is obtained by taking the average value of the difference between the first amplitude and the following 10 amplitudes in (n), and the formula is:

[0109]

[0110] S4-4: Using the normalized square envelope spectrum NS ik (f) Calculate the position accuracy index PA of the fault characteristic frequency ik ; It also includes the following steps;

[0111] S4-4-1: Find the normalized square envelope spectrum NS ik The frequency value f corresponding to the maximum amplitude in (f) ik_ma ;

[0112] S4-4-2: Calculate the position accuracy index PA ik , considering the three rotating parts of the rolling bearing, the inner ring, the outer ring and the rolling elements, then PA ik The calculation formula is:

[0113] PA ik =(|f ik_ma -f ir |·|f ik_ma -f or |·|f ik_ma -f ba |) 1 / 3

[0114] Among them, f ir 、f or and f ba They are respectively the theoretical fault characteristic frequencies of the inner ring, outer ring and rolling element of the rolling bearing calculated above.

[0115] S4-5: Calculate the modal component u ik (t) The corresponding energy position index EP ik , the formula is:

[0116]

[0117] Where p is the adjustment coefficient, and p=1, β is the balance coefficient, and the calculation formula of β is:

[0118]

[0119] S5: From K i Energy position index EP i Find the minimum value minEP in (k) i and the serial number of its corresponding target modal component i ;

[0120] S6: According to the final N s minEP i Value, draw α i With minEP i The relationship curve between Figure 4 As shown, draw α i With ON i The relationship curve between Figure 5 As shown;

[0121] S7: According to Figure 4 α shown i With minEP i From the relationship curve between minEP and α, we can see that as the equilibrium parameter α increases, minEP i The α value corresponding to the minimum value is 6650, which is taken as the optimal value of the balance parameter α opt =6650;

[0122] S8: According to Figure 5 According to α i With ON i The relationship curve between α i The value is α opt =6650, the corresponding ON i The value is ON T =8;

[0123] S9: With α opt= 6650 as the value of the balance parameter α, perform continuous variational mode decomposition SVMD, obtain multiple modal components, and select the ONth modal component from them. T =8 modal components are the optimal target modal components u opt (t), such as Figure 6 As shown;

[0124] S10: Calculate the optimal target modal component u opt The square envelope spectrum of (t) is as follows: Figure 7 As shown in the figure, the extracted fault characteristic frequency is f = 102 Hz, which is consistent with the theoretical fault characteristic frequency f of the outer ring of the rolling bearing. or =102.824Hz, indicating that the outer ring of the rolling bearing has indeed been damaged.

[0125] The contents described in the embodiments of this specification are merely an enumeration of the implementation forms of the inventive concept. The protection scope of the present invention should not be regarded as limited to the specific forms described in the embodiments. The protection scope of the present invention also extends to equivalent technical means that can be conceived by those skilled in the art based on the inventive concept.

Claims

1. Rotating machinery fault diagnosis method based on improved continuous variational modal decomposition, The following steps are involved: S1: On rotating machinery, s is the sampling frequency, and a vibration signal x(t) with a length of L is collected; S2: Set the range of the balance parameter α to [α min ,α max ], let α be from α min Start with a step size of s α Increase, and the value when it increases to the i-th step is: a i =a min +(i-1)·s α (I) Where i = 1, 2, 3, ..., N s , N s is the total number of steps that α increases, and N s =(α max -α min ) / s α +1; S3: With α i As a balance parameter, continuous variational mode decomposition SVMD is performed to obtain multiple modal components u ik (t), where k is the serial number of the modal component, k = 1, 2, ..., K i , K i is the number of modal components obtained by the continuous variational mode decomposition (SVMD) in the i-th step; S4: Calculate the modal component u ik Energy position index value EP of (t) i (k); S5: From K i Energy position index EP i Find the minimum value minEP in (k) i and the serial number of its corresponding target modal component i ; S6: According to the final N s minEP i Value, draw α i With minEP i and draw the relationship curve between α i With ON i The relationship curve between S7: Using α i With minEP i The relationship curve between minEP and α is used to find out the relationship between minEP and α when the equilibrium parameter α increases. i The α value corresponding to the minimum value is used as the optimal value of the balance parameter α opt ; S8: Using α i With ON i The relationship curve between i The value is α opt Corresponding ON i The value of ON T ; S9: With α opt As the value of the balance parameter α, perform continuous variational mode decomposition SVMD to obtain multiple modal components, from which the ONth modal component is selected. T The modal component is the optimal target modal component u opt (t); S10: Calculate the optimal target modal component u opt The square envelope spectrum of (t) is used to extract the fault characteristic frequency and perform fault diagnosis on rotating machinery; In step S4 and step S5, the calculation of the energy position index value of the modal component includes the following steps: S2-1: For a certain modal component u(t), calculate its square envelope spectrum SES(f), where f represents the frequency; S2-2: Normalize the square envelope spectrum SES(f) to obtain the normalized square envelope spectrum NS(f). The specific method is: Wherein, max[SES(f)] represents the maximum amplitude in the square envelope spectrum SES(f); S2-3: Calculate the energy concentration index EC of the fault characteristic frequency component using the normalized square envelope spectrum NS(f); S2-4: Calculate the position accuracy index PA of the fault characteristic frequency using the normalized square envelope spectrum NS(f); S2-5: Calculate the energy position index EP corresponding to the modal component u(t), the formula is: Where p is the adjustment coefficient, and 0<p≤1, β is the balance coefficient, and the calculation formula of β is: In step S2-3 and step S2-5, the calculation of the energy concentration index EC of the fault characteristic frequency component includes the following steps: S4-1: rearrange the amplitudes of the normalized square envelope spectrum NS(f) from large to small to obtain a rearranged amplitude sequence SNS(n), n = 1, 2, ..., N, where N represents the number of frequency points; S4-2: Calculate the average value of the difference between the first amplitude and the subsequent M amplitudes in SNS(n) to obtain the energy concentration index EC, the formula is: In step S2-4 and step S2-5, the calculation of the position accuracy index PA of the fault characteristic frequency includes the following steps: S5-1: Find the frequency value f corresponding to the maximum amplitude in the normalized square envelope spectrum NS(f) ma ; S5-2: Calculate the position accuracy index PA, the formula is: Where C is the number of rotating parts in the rotating machinery, f j is the theoretical fault characteristic frequency of each rotating component, j=1,2,…,C.

2. The rotating machinery fault diagnosis method based on improved continuous variational modal decomposition as claimed in claim 1, It is characterized in that In step S10 and step S2-1, the calculation of the square envelope spectrum of the modal component includes the following steps: S3-1: For a certain modal component u(t), calculate its square envelope signal u SE (t), the formula is: Where j is the imaginary unit, represents Hilbert transform, |·| represents the modulus of a complex number; S3-2: Calculate u SE The square envelope spectrum SES(f) of (t) is given by: in, represents Fourier transform.

Images