CN109934206A - A kind of rotary machinery fault diagnosis method under non-stationary operating condition - Google Patents
A kind of rotary machinery fault diagnosis method under non-stationary operating condition Download PDFInfo
- Publication number
- CN109934206A CN109934206A CN201910276504.9A CN201910276504A CN109934206A CN 109934206 A CN109934206 A CN 109934206A CN 201910276504 A CN201910276504 A CN 201910276504A CN 109934206 A CN109934206 A CN 109934206A
- Authority
- CN
- China
- Prior art keywords
- frequency
- index
- domain
- value
- time
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
- 238000003745 diagnosis Methods 0.000 title claims abstract description 45
- 238000000034 method Methods 0.000 title claims abstract description 35
- 230000002123 temporal effect Effects 0.000 claims abstract description 30
- 238000001228 spectrum Methods 0.000 claims abstract description 17
- 230000001360 synchronised effect Effects 0.000 claims abstract description 16
- 230000006835 compression Effects 0.000 claims abstract description 13
- 238000007906 compression Methods 0.000 claims abstract description 13
- 238000005457 optimization Methods 0.000 claims abstract description 11
- 239000000284 extract Substances 0.000 claims abstract description 9
- 239000002245 particle Substances 0.000 claims abstract description 7
- 238000012706 support-vector machine Methods 0.000 claims description 16
- 238000009826 distribution Methods 0.000 claims description 14
- 238000012545 processing Methods 0.000 claims description 11
- 230000009466 transformation Effects 0.000 claims description 9
- 238000000605 extraction Methods 0.000 claims description 7
- 230000004927 fusion Effects 0.000 claims description 7
- 238000004458 analytical method Methods 0.000 claims description 6
- 238000005452 bending Methods 0.000 claims description 6
- 238000004364 calculation method Methods 0.000 claims description 6
- 230000003595 spectral effect Effects 0.000 claims description 6
- 238000012360 testing method Methods 0.000 claims description 6
- 238000012549 training Methods 0.000 claims description 6
- 230000001133 acceleration Effects 0.000 claims description 3
- 238000010606 normalization Methods 0.000 claims description 3
- 238000006073 displacement reaction Methods 0.000 claims description 2
- 239000011159 matrix material Substances 0.000 claims 1
- 238000005516 engineering process Methods 0.000 description 5
- 230000006870 function Effects 0.000 description 4
- 230000008901 benefit Effects 0.000 description 2
- 230000008859 change Effects 0.000 description 2
- 238000004040 coloring Methods 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 235000013399 edible fruits Nutrition 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 241000208340 Araliaceae Species 0.000 description 1
- 235000005035 Panax pseudoginseng ssp. pseudoginseng Nutrition 0.000 description 1
- 235000003140 Panax quinquefolius Nutrition 0.000 description 1
- 230000002159 abnormal effect Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000002405 diagnostic procedure Methods 0.000 description 1
- 230000002349 favourable effect Effects 0.000 description 1
- 238000007306 functionalization reaction Methods 0.000 description 1
- 235000008434 ginseng Nutrition 0.000 description 1
- 230000036541 health Effects 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 238000012544 monitoring process Methods 0.000 description 1
- 238000007500 overflow downdraw method Methods 0.000 description 1
- 230000008447 perception Effects 0.000 description 1
- 230000035945 sensitivity Effects 0.000 description 1
- 230000019491 signal transduction Effects 0.000 description 1
Abstract
The invention discloses the rotary machinery fault diagnosis methods under a kind of non-stationary operating condition, N rank synchronous compression is carried out to original vibration signal first to convert to obtain high-precision time-frequency spectrum, then time-frequency energy accumulating region is found out by region of interest (ROI) extracted region algorithm, extract temporal signatures, temporal frequency domain feature and time-frequency image feature comprehensively disclose the operating status of equipment, then optimization of the optimized parameter realization to supporting vector machine model is filtered out using modified particle swarm optiziation, the independent tentative diagnosis result in each domain is obtained by SVM, the failure modes result in three domains is finally carried out information using D-S evidence theory to merge to obtain final judging result.Rotary machinery fault diagnosis method provided by the invention, largely improves the accuracy rate of fault diagnosis, while improving the robustness of judging result.
Description
Technical field
The present invention relates to the rotary machinery fault diagnosis sides under Diagnosis Technique field, especially non-stationary operating condition
Method.
Background technique
With the rapid development of science and technology, mechanical equipment in modern industry is more and more diversified, functionalization, complication.Rotation
Favourable turn tool (such as gas compressor, wind turbine and aero-engine) is one of most important equipment in modern industry application,
Industry especially is being manufactured, carrying out production using the even abnormal equipment of any failure will inevitably impact on product
Final mass leads to weight huge economic loss, or even can cause work safety accident, seriously threatens the life security of personnel, because
This needs to carry out timely health monitoring and fault diagnosis to rotating machinery, prevents trouble before it happens.
Previous research is largely the fault diagnosis to rotating machinery under steady operating condition, but in modern project scene,
Rotating machinery working environment variation multiplicity, operating condition is (such as Wind turbines) also complicated and changeable: when varying load, when variable speed and wink
State phenomenon etc., while being influenced by interacting between noise, signaling pathways and inner body, vibration signal is modulated
The interference layer by layer of frequency modulation width, noise and other vibration sources etc. shows plyability, non-linear, non-stationary.
The time-frequency for the non-stationary signal that nonlinear change is presented in existing Time-Frequency Analysis Method processing at any time at present is differentiated
Rate is poor, be easy to appear mode obscure, spurious energy the problems such as, while using the fault signature in single piece of information source solve it is practical
It is easy to expose diagnostic message incomplete drawback when engineering problem, and then shadow is caused to the accuracy and robustness of diagnostic result
It rings, it is therefore desirable to carry out the research of the signal processing and method for diagnosing faults of the Non-stationary vibration signal of rotating machinery, to improve
Precision, sensitivity, the precision of prediction and the reliability of early warning of rotating machinery fault perception of time-frequency representation.
Summary of the invention
In view of the above-mentioned problems, the object of the present invention is to provide rotary machinery fault diagnosis method under a kind of non-stationary operating condition,
This method is based on the transformation of N rank synchronous compression and obtains high-precision time-frequency spectrum, and by the fault diagnosis side of Multi-source Information Fusion
Method is to improve the precision and reliability of fault diagnosis.
In order to achieve the above objectives, the invention adopts the following technical scheme:
A kind of rotary machinery fault diagnosis method under non-stationary operating condition, the method includes the steps of:
1) a certain state H of equipment is obtained by the vibration acceleration sensor being mounted on slewingiUnder (1≤i≤k)
Vibration signal X, the state HiIt can be the state under normal condition, various typical fault types, randomly select several groups work
For sample data, remaining each group is as sample data to be detected；
2) transformation of N rank synchronous compression is carried out to original signal X, obtains high-precision time-frequency spectrum；
3) temporal signatures E is extracted1, temporal frequency domain feature E2, time-frequency image feature E3；
4) support vector machines (SVM) prediction model in each single domain feature of improved population (PSO) algorithm optimization is utilized；
5) optimized parameter C and g are inputted in SVM classifier, selects suitable kernel function, training dataset makes it have most
Excellent feature is then applied to test data set and obtains failure modes result and classification accuracy, uses the SVM model after optimization point
It Li Yong not temporal signatures E1, temporal frequency domain feature E2With time-frequency image feature E3Training sample be trained make it have respectively
Most there is feature, is then respectively applied to the test sample in domain, carries out independent tentative diagnosis and obtain failure modes result and rotation
Make a connection classification accuracy P of the tool under different characteristic domainSVMi, i.e., each domain independence tentative diagnosis as a result, this classification accuracy PSVMi
It is Basic probability assignment function, other WSVMefExpression virtual condition is HeIt is mistaken for state HfProbability (1≤e, f≤k)；
6) based on D-S evidence theory to time domain E1, temporal frequency domain E2With time-frequency image domain E3Local diagnosis knot in three domains
Fruit carries out information fusion；
7) according to three domain information convergence analysis probability distribution PDSiProvide the final judging result of fault diagnosis, i.e. PDSiValue is most
Corresponding bearing operating status is last diagnostic state when big, it is hereby achieved that rotating machinery has fault-free and failure mode
Judgement.
2. the rotary machinery fault diagnosis method under a kind of non-stationary operating condition according to claim 1, feature exist
In: the transformation of N rank synchronous compression is carried out to obtain high-precision time-frequency spectrum to original signal X in the step 2), is specifically converted
It is as follows:
If the FM amplitude modulation signal X under non-stationary operating condition is defined as: f (τ)=A (τ) ei2πφ(τ), in which: A (τ) indicates wink
When amplitudeφ (τ) indicates instantaneous phaseA
(τ) and φ (τ) are higher order functionality, φ(k)(t) k local derviation of the φ about t is indicated, then the vibration under the non-stationary operating condition
Signal is converted by the synchronous transformed time-frequency Energy distribution of extraction of high-order by formula (1):
Each parameter is respectively as follows: in the formula (1)
Indicate the Short Time Fourier Transform of signal:
Indicate N rank instantaneous Frequency Estimation operator:
ωf(t, η) indicates synchronous compression instantaneous frequency:
Indicate N contrast operator:
3. the rotary machinery fault diagnosis method under a kind of non-stationary operating condition according to claim 1, feature exist
In: the step 3) is subdivided into following two step:
3-1): using region of interest (region of interest, ROI) processing technique, i.e., extracted and calculated using ROI region
The characteristic area of method quick obtaining most worthy, i.e., in the crestal line region of time-frequency energy accumulating, to promote follow-up data processing
Speed；
3-2): extracting temporal signatures E1Index has 13: absolute mean, peak value, root-mean-square value, root amplitude, variance, peak
- peak value, degree of skewness index, kurtosis index, peak index, waveform index, pulse index, margin index, the coefficient of variation；It extracts
11 temporal frequency domain feature E2Index: average frequency, square frequency, root mean square frequency, frequency variance, frequency standard be poor, spectral peak
Index of stability, and according to the frequency feature of equipment different conditions, frequency domain is divided equally into 5 frequency bands, each band bandwidth is Bf, point
The relative power spectrum energy of each frequency band is not calculated；Extract time-frequency image feature E3Index has: color second moment, Hu not displacement,
Gray level co-occurrence matrixes, brightness index, Tamura texture roughness index；The calculation expression of each index is as follows:
Absolute mean:Peak value: xp=max | xi|；
Root-mean-square value:Root amplitude:
Variance:Peak-to-peak value: xp-p=max (xi)-min(xi)；
Degree of skewness index:Kurtosis index:
Peak index:Waveform index:
Pulse index:Margin index:
The coefficient of variation:Average frequency:
Square frequency:Root mean square frequency:
Frequency variance:
Frequency standard is poor:
Spectral peak index of stability:
First band relative energy:
Second band relative energy:
Third frequency band relative energy:
4th frequency band relative energy:
5th frequency band relative energy:
F indicates the frequency of signal in each characteristic index calculation formula, and p (f) indicates the power spectrum of signal, BfIndicate 1/
5 frequency range values, FsIndicate highest frequency value；
Color second momentWherein, pi,jIndicate i-th of the face of j-th of pixel of color image
Colouring component, N indicate the number of pixels in image；
Bending moment does not share 7 to Hu, is typically chosen first 2 not bending moments: M1=η20+η02,Its
Middle ηpqIndicate normalized center away from；
Gray level co-occurrence matrixesWherein P (i, j, d, θ) indicates the gray scale from gray level image
Rank be i pixel position (x, y) set out, statistics with its distance be d, pixel position (x+ Δ x, the y+ Δ y) that grey level is j
The frequency P (i, j, d, θ) occurred simultaneously=and [(x, y), (x+ Δ x, y+ Δ y)] | f (x, y)=i, f (x+ Δ x, y+ Δ y)=
j}；
Brightness I=(r+g+b)/3, wherein r, g and b are respectively red, green and blue 3 single Color Channels in input picture
Pixel value；
Tamura texture roughness FcrsReflect the severe degree of variation of image grayscale, the texture primitive size the big, feels
Feel more coarse, calculates in image centered on coordinate (x, y), the active window (such as 1 × 1,2 × 2 ..., 32 that size is 2k × 2k
× 32), in window each pixel gray averageWherein (i, j) indicates neighborhood
The position of middle pixel；G (i, j) indicates the gray value for being located at (i, j) pixel；K determines the pixel coverage for participating in mean value computation；For
Each pixel, the average gray difference calculated separately between the active window that it is not overlapped on both horizontally and vertically are other
Are as follows: Ek,v=| Ak(x,y+2k-1)-Ak(x,y-2k-1) | and Ek,h=| Ak(x+2k-1,y)-Ak(x-2k-1, y) |, for each picture
Member finds optimal size Skbest=2k, its average gray value E value is made to reach maximumThat
Roughness is defined as each pixel optimal size S in entire imagekbestAverage value:M
Picture traverse and height are respectively indicated with N.
4. the rotary machinery fault diagnosis method under a kind of non-stationary operating condition according to claim 1, feature exist
In: it is predicted in the step 4) using the support vector machines (SVM) under each single domain feature of improved population (PSO) algorithm optimization
Model, the optimization screen suitable penalty factor and nuclear parameter g, described two ginsengs using improved population PSO algorithm
Number C and g is to influence two important parameters of svm classifier performance, and the step 4) is subdivided into following four step:
4-1) parameter initialization: setting Population Size, itself Studying factors c1 and sociology department factor c2, the number of iterations K,
The maximum value and minimum value of inertia weight ω, the value range of parameter C and g；
The position and speed of primary 4-2) is randomly generated: calculating the fitness size of each particle, sequence is found a
Body extreme point and global extreme point calculates the average fitness of every generation population；
4-3) update particle rapidity and location information；
It 4-4) examines termination condition: judging whether to meet maximum number of iterations (or meeting required precision), if meeting
Penalty factor and nuclear parameter g are exported, otherwise return step b.
5. the rotary machinery fault diagnosis method under a kind of non-stationary operating condition according to claim 1, feature exist
Be subdivided into: the step 6) the following three steps:
6-1) calculate normalization factor K value:
6-2) calculate a certain state H when time domain and temporal frequency domain Fusion FeaturesiProbability distribution:With
6-3) the result for after the same method merging the time domain and temporal frequency domainAnd time-frequency image
Property field merges again, and then calculates and diagnose operating status H by three domain evidence fusionsiProbability distribution obtain
Compared with the existing technology, the present invention has following remarkable advantage:
1. the present invention proposes to carry out letter to the characteristic of rotating machines vibration signal under non-stationary operating condition using the transformation of n rank synchronous compression
Number analysis, can preferably reflect that non-linear fast-changing multi -components are believed by carrying out high-order approximation simultaneously to amplitude and phase
Number, the time-frequency Energy distribution result that this method obtains more connects with video Energy distribution ideally compared with the existing technology
Closely, really reflect the trend that instantaneous frequency and amplitude change over time, the resolution ratio of time-frequency representation is higher, non-linear more handling
Component FM amplitude modulation signal has distinctive advantage, provides high-precision time-frequency spectrum for fault diagnosis；
2. the method for diagnosing faults of multi-domain characteristics fusion proposed by the present invention can make full use of time domain, temporal frequency domain and time-frequency
Feature in three domains of image area, the opposite characteristic index being used alone in each domain are diagnosed, and can be obtained more valuable
Judge information, improve the reliability of fault diagnosis result；Before image characteristics extraction, propose to utilize region of interest
(ROI) crestal line regional scope is recognized accurately in extracted region algorithm, can effectively accelerate the speed of subsequent image processing；In addition should
Method also proposes the mode by two layer diagnosis, solves the problems, such as that single domain feature diagnostic method accuracy rate is low, stability is poor, protects
The accuracy of fault diagnosis is demonstrate,proved.
Characteristic of rotating machines vibration signal analysis and fault diagnosis of this method suitable for fault diagnosis field.
Detailed description of the invention
Fig. 1 is the rotary machinery fault diagnosis method flow chart under a kind of non-stationary operating condition of the present invention；
The time-frequency spectrum of Fig. 2 difference time frequency processing technology；
Fig. 3 extracts the region of time-frequency energy accumulating using region of interest (ROI) extracted region algorithm.
Specific embodiment
Present invention will be further explained below with reference to the attached drawings and specific embodiments.
As shown in Figure 1, the rotary machinery fault diagnosis method under a kind of non-stationary operating condition provided by the invention, comprising following
Step:
1) a certain state H of equipment is obtained by the vibration acceleration sensor being mounted on slewingiUnder (1≤i≤k)
Vibration signal X, the state HiIt can be the state under normal condition, various typical fault types, randomly select several groups work
For sample data, remaining each group is as sample data to be detected；
2) transformation of N rank synchronous compression is carried out to original signal X, obtains high-precision time-frequency spectrum, specific transformation is as follows:
If the FM amplitude modulation signal X under non-stationary operating condition is defined as: f (τ)=A (τ) ei2πφ(τ), in which: A (τ) indicates wink
When amplitudeφ (τ) indicates instantaneous phaseA
(τ) and φ (τ) are higher order functionality, φ(k)(t) k local derviation of the φ about t is indicated, then the vibration under the non-stationary operating condition
Signal is converted by the synchronous transformed time-frequency Energy distribution of extraction of high-order by formula (1):
Each parameter is respectively as follows: in the formula (1)
Indicate the Short Time Fourier Transform of signal:
Indicate N rank instantaneous Frequency Estimation operator:
ωf(t, η) indicates synchronous compression instantaneous frequency:
Indicate N contrast operator:
As shown in Fig. 2, the figure is only a schematic diagram, to show the different resulting time-frequency spectrums of time frequency processing technology
Difference, Fig. 2 (a) are Short Time Fourier Transform (STFT) treated time-frequency figure, it can be seen that time-frequency figure is very fuzzy, and there are moulds
State aliasing, time-frequency energy accumulating are poor, and Fig. 2 (b) is the N rank synchronous compression time-frequency that converts that (wherein when order n=3) treated
Figure, it can be seen that the figure time-frequency energy accumulating, crestal line is clear, and precision is high, can preferably provide for Fault Diagnosis of Rotating Equipment Based
Reliable Time-Frequency Information.The frequency structure that Fig. 2 shows is relatively easy, the displaying of only one of situation, actual vibration letter
Number multiple bar chart time-frequency energy accumulating band may may be showed comprising complicated frequency structure；
3) on the basis of the time-frequency spectrum that the step 2) obtains, temporal signatures E is extracted1, temporal frequency domain feature E2, when
Frequency characteristics of image E3, it can specifically be subdivided into two steps:
3-1): using region of interest (region of interest, ROI) processing technique, i.e., extracted and calculated using ROI region
The characteristic area of method quick obtaining most worthy, i.e., in the crestal line region of time-frequency energy accumulating, to promote follow-up data processing
Speed；If Fig. 3 is also schematic diagram, the region in figure between two yellow curves is after Fig. 2 (b) N rank synchronous compression converts
The region of interest that extracts of time-frequency figure, that is, the region of time-frequency energy accumulating, the few time-frequency energy except two yellow curves
It is low to ignore in subsequent feature extraction, feature extraction is carried out just for the part in region, therefore can be significantly
Promote the precision of calculating speed and fault diagnosis.ROI region extraction algorithm is prior art, and details are not described herein.
3-2): extracting temporal signatures E1Index has 13: absolute mean, peak value, root-mean-square value, root amplitude, variance, peak
- peak value, degree of skewness index, kurtosis index, peak index, waveform index, pulse index, margin index, the coefficient of variation；It extracts
11 temporal frequency domain feature E2Index: average frequency, square frequency, root mean square frequency, frequency variance, frequency standard be poor, spectral peak
Index of stability, and according to the frequency feature of equipment different conditions, frequency domain is divided equally into 5 frequency bands, each band bandwidth is Bf, point
The relative power spectrum energy of each frequency band is not calculated；The calculation expression of each index is as follows:
Absolute mean:Peak value: xp=max | xi|；
Root-mean-square value:Root amplitude:
Variance:Peak-to-peak value: xp-p=max (xi)-min(xi)；
Degree of skewness index:Kurtosis index:
Peak index:Waveform index:
Pulse index:Margin index:
The coefficient of variation:Average frequency:
Square frequency:Root mean square frequency:
Frequency variance:
Frequency standard is poor:
Spectral peak index of stability:
First band relative energy:
Second band relative energy:
Third frequency band relative energy:
4th frequency band relative energy:
5th frequency band relative energy:
F indicates the frequency of signal in each characteristic index calculation formula, and p (f) indicates the power spectrum of signal, BfIndicate 1/
5 frequency range values, FsIndicate highest frequency value；
Color second momentWherein, pi,jIndicate i-th of the face of j-th of pixel of color image
Colouring component, N indicate the number of pixels in image；
Bending moment does not share 7 to Hu, is typically chosen first 2 not bending moments: M1=η20+η02,Its
Middle ηpqIndicate normalized center away from；
Gray level co-occurrence matrixesWherein P (i, j, d, θ) indicates the gray scale from gray level image
Rank be i pixel position (x, y) set out, statistics with its distance be d, pixel position (x+ Δ x, the y+ Δ y) that grey level is j
The frequency P (i, j, d, θ) occurred simultaneously=and [(x, y), (x+ Δ x, y+ Δ y)] | f (x, y)=i, f (x+ Δ x, y+ Δ y)=
j}；
Brightness I=(r+g+b)/3, wherein r, g and b are respectively red, green and blue 3 single Color Channels in input picture
Pixel value；
Tamura texture roughness FcrsReflect the severe degree of variation of image grayscale, the texture primitive size the big, feels
Feel more coarse, calculates in image centered on coordinate (x, y), the active window (such as 1 × 1,2 × 2 ..., 32 that size is 2k × 2k
× 32), in window each pixel gray averageWherein (i, j) indicates neighborhood
The position of middle pixel；G (i, j) indicates the gray value for being located at (i, j) pixel；K determines the pixel coverage for participating in mean value computation；For
Each pixel, the average gray difference calculated separately between the active window that it is not overlapped on both horizontally and vertically are other
Are as follows: Ek,v=| Ak(x,y+2k-1)-Ak(x,y-2k-1) | and Ek,h=| Ak(x+2k-1,y)-Ak(x-2k-1, y) |, for each picture
Member finds optimal size Skbest=2k, its average gray value E value is made to reach maximumSo
Roughness is defined as each pixel optimal size S in entire imagekbestAverage value:M and
N respectively indicates picture traverse and height.
4) support vector machines (SVM) prediction model under each single domain feature of improved population (PSO) algorithm optimization is utilized；
Suitable penalty factor is screened using improved population PSO algorithm and nuclear parameter g, described two parameters are shadows
Ring two important parameters of svm classifier performance；
4-1) parameter initialization.Population Size, itself Studying factors c1 and sociology department factor c2, the number of iterations K are set,
The maximum value and minimum value of inertia weight ω, the value range of parameter C and g；
The position and speed of primary 4-2) is randomly generated.The fitness size of each particle is calculated, sequence is found a
Body extreme point and global extreme point calculates the average fitness of every generation population；
4-3) update particle rapidity and location information；
4-4) examine termination condition.Judge whether to meet maximum number of iterations (or meeting required precision), if meeting
Penalty factor and nuclear parameter g are exported, otherwise return step b.
5) optimized parameter C and g are inputted in SVM classifier, selects suitable kernel function, training dataset makes it have most
Excellent feature is then applied to test data set and obtains failure modes result and classification accuracy, uses the SVM model after optimization point
It Li Yong not temporal signatures E1, temporal frequency domain feature E2With time-frequency image feature E3Training sample be trained make it have respectively
Most there is feature, is then respectively applied to the test sample in domain, carries out independent tentative diagnosis and obtain failure modes result and rotation
Make a connection classification accuracy P of the tool under different characteristic domainSVMi, i.e., each domain independence tentative diagnosis as a result, this classification accuracy PSVMi
It is Basic probability assignment function, other WSVMefExpression virtual condition is HeIt is mistaken for state HfProbability (1≤e, f≤k),
It is as shown in table 1:
6) based on D-S evidence theory to time domain E1, temporal frequency domain E2With time-frequency image domain E3Local diagnosis knot in three domains
Fruit carry out information fusion, be subdivided into the following three steps:
6-1) calculate normalization factor K value:
6-2) calculate a certain state H when time domain and temporal frequency domain Fusion FeaturesiProbability distribution:With
6-3) the result for after the same method merging the time domain and temporal frequency domainAnd time-frequency image
Property field merges again, and then calculates and diagnose operating status H by three domain evidence fusionsiProbability distribution obtain
7) according to three domain information convergence analysis probability distribution PDSiProvide the final judging result of fault diagnosis, i.e. PDSiValue is most
Corresponding bearing operating status is last diagnostic state when big, it is hereby achieved that rotating machinery has fault-free and failure mode
Judgement.
Claims (5)
1. the rotary machinery fault diagnosis method under a kind of non-stationary operating condition, which is characterized in that comprise the steps of:
1) a certain state H of equipment is obtained by the vibration acceleration sensor being mounted on slewingiVibration under (1≤i≤k)
Dynamic signal X, the state HiIt can be the state under normal condition, various typical fault types, randomly select several groups as sample
Notebook data, remaining each group is as sample data to be detected；
2) transformation of N rank synchronous compression is carried out to original signal X, obtains high-precision time-frequency spectrum；
3) temporal signatures E is extracted1, temporal frequency domain feature E2, time-frequency image feature E3；
4) support vector machines (SVM) prediction model in each single domain feature of improved population (PSO) algorithm optimization is utilized；
5) optimized parameter C and g are inputted in SVM classifier, selects suitable kernel function, training dataset makes it have optimal spy
Sign, is then applied to test data set and obtains failure modes result and classification accuracy, sharp respectively using the SVM model after optimization
With temporal signatures E1, temporal frequency domain feature E2With time-frequency image feature E3Training sample be trained to make it have respectively and most have
Then feature is respectively applied to the test sample in domain, carries out independent tentative diagnosis and obtain failure modes result and whirler
Classification accuracy P of the tool under different characteristic domainSVMi, i.e., each domain independence tentative diagnosis as a result, this classification accuracy PSVMiIt is
Basic probability assignment function, in addition WSVMefExpression virtual condition is HeIt is mistaken for state HfProbability (1≤e, f≤k)；
6) based on D-S evidence theory to time domain E1, temporal frequency domain E2With time-frequency image domain E3Local diagnosis result in three domains into
Row information fusion；
7) according to three domain information convergence analysis probability distribution PDSiProvide the final judging result of fault diagnosis, i.e. PDSiWhen value is maximum
Corresponding bearing operating status is last diagnostic state, it is hereby achieved that rotating machinery has sentencing for fault-free and failure mode
It is disconnected.
2. the rotary machinery fault diagnosis method under a kind of non-stationary operating condition according to claim 1, it is characterised in that: institute
It states in step 2) and the transformation of N rank synchronous compression is carried out to obtain high-precision time-frequency spectrum to original signal X, specific transformation is as follows:
If the FM amplitude modulation signal X under non-stationary operating condition is defined as: f (τ)=A (τ) ei2πφ(τ), in which: A (τ) indicates instantaneous vibration
Widthφ (τ) indicates instantaneous phaseA(τ)
It is higher order functionality, φ with φ (τ)(k)(t) k local derviation of the φ about t is indicated, then the vibration under the non-stationary operating condition is believed
It number extracts transformed time-frequency Energy distribution by high-order is synchronous and is converted by formula (1):
Each parameter is respectively as follows: in the formula (1)
Indicate the Short Time Fourier Transform of signal:
Indicate N rank instantaneous Frequency Estimation operator:
ωf(t, η) indicates synchronous compression instantaneous frequency:
Indicate N contrast operator:
3. the rotary machinery fault diagnosis method under a kind of non-stationary operating condition according to claim 1, it is characterised in that: institute
It states step 3) and is subdivided into following two step:
3-1): use region of interest (region of interest, ROI) processing technique, i.e., it is fast using ROI region extraction algorithm
Speed obtains the characteristic area of most worthy, i.e., in the crestal line region of time-frequency energy accumulating, to promote the speed of follow-up data processing；
3-2): extracting temporal signatures E1Index has 13: absolute mean, peak value, root-mean-square value, root amplitude, variance, peak-peak
Value, degree of skewness index, kurtosis index, peak index, waveform index, pulse index, margin index, the coefficient of variation；Extract 11
Temporal frequency domain feature E2Index: average frequency, square frequency, root mean square frequency, frequency variance, frequency standard is poor, spectral peak is stable
Index, and according to the frequency feature of equipment different conditions, frequency domain is divided equally into 5 frequency bands, each band bandwidth is Bf, count respectively
Calculate the relative power spectrum energy of each frequency band；Extract time-frequency image feature E3Index has: color second moment, Hu not displacement, gray scale
Co-occurrence matrix, brightness index, Tamura texture roughness index；The calculation expression of each index is as follows:
Absolute mean:Peak value: xp=max | xi|；
Root-mean-square value:Root amplitude:
Variance:Peak-to-peak value: xp-p=max (xi)-min(xi)；
Degree of skewness index:Kurtosis index:
Peak index:Waveform index:
Pulse index:Margin index:
The coefficient of variation:Average frequency:
Square frequency:Root mean square frequency:
Frequency variance:
Frequency standard is poor:
Spectral peak index of stability:
First band relative energy:
Second band relative energy:
Third frequency band relative energy:
4th frequency band relative energy:
5th frequency band relative energy:
F indicates the frequency of signal in each characteristic index calculation formula, and p (f) indicates the power spectrum of signal, BfIndicate 1/5 frequency range
Value, FsIndicate highest frequency value；
Color second momentWherein, pi,jIndicate i-th of color point of j-th of pixel of color image
Amount, N indicate the number of pixels in image；
Bending moment does not share 7 to Hu, is typically chosen first 2 not bending moments: M1=η20+η02,Wherein ηpq
Indicate normalized center away from；
Gray level co-occurrence matrixesWherein P (i, j, d, θ) indicates the grey level from gray level image
It sets out for the pixel position (x, y) of i, (x+ Δ x, y+ Δ y) is simultaneously with its pixel position that distance is d, grey level is j for statistics
The frequency P (i, j, d, θ) of appearance=and [(x, y), (x+ Δ x, y+ Δ y)] | f (x, y)=i, f (x+ Δ x, y+ Δ y)=j }；
Brightness I=(r+g+b)/3, wherein r, g and b are respectively red, green and blue 3 single Color Channel pixels in input picture
Value；
Tamura texture roughness FcrsReflect variation of image grayscale severe degree, the texture primitive size the big, feel more
It is coarse, calculate image in centered on coordinate (x, y), size be 2k × 2k active window (such as 1 × 1,2 × 2 ..., 32 ×
32), in window each pixel gray averageWherein (i, j) is indicated in neighborhood
The position of pixel；G (i, j) indicates the gray value for being located at (i, j) pixel；K determines the pixel coverage for participating in mean value computation；For every
A pixel, the average gray difference calculated separately between the active window that it is not overlapped on both horizontally and vertically are respectively as follows:
Ek,v=| Ak(x,y+2k-1)-Ak(x,y-2k-1) | and Ek,h=| Ak(x+2k-1,y)-Ak(x-2k-1, y) |, for each pixel,
Find optimal size Skbest=2k, its average gray value E value is made to reach maximumIt is so thick
Rugosity is defined as each pixel optimal size S in entire imagekbestAverage value:M and N
Respectively indicate picture traverse and height.
4. the rotary machinery fault diagnosis method under a kind of non-stationary operating condition according to claim 1, it is characterised in that: institute
Support vector machines (SVM) prediction model utilized under each single domain feature of improved population (PSO) algorithm optimization in step 4) is stated,
The optimization i.e. using improved population PSO algorithm screen suitable penalty factor and nuclear parameter g, described two parameter C and
G is to influence two important parameters of svm classifier performance, and the step 4) is subdivided into following four step:
4-1) parameter initialization: setting Population Size, itself Studying factors c1 and sociology department factor c2, the number of iterations K, inertia
The maximum value and minimum value of weights omega, the value range of parameter C and g；
The position and speed of primary 4-2) is randomly generated: calculating the fitness size of each particle, individual pole is found in sequence
Value point and global extreme point, calculate the average fitness of every generation population；
4-3) update particle rapidity and location information；
It 4-4) examines termination condition: judging whether to meet maximum number of iterations (or meeting required precision), exported if meeting
Penalty factor and nuclear parameter g, otherwise return step b.
5. the rotary machinery fault diagnosis method under a kind of non-stationary operating condition according to claim 1, it is characterised in that: institute
State step 6) be subdivided into the following three steps:
6-1) calculate normalization factor K value:
6-2) calculate a certain state H when time domain and temporal frequency domain Fusion FeaturesiProbability distribution:With
6-3) the result for after the same method merging the time domain and temporal frequency domainWith time-frequency image feature
Domain is merged again, and then is calculated and diagnosed operating status H by three domain evidence fusionsiProbability distribution obtain
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910276504.9A CN109934206A (en) | 2019-04-08 | 2019-04-08 | A kind of rotary machinery fault diagnosis method under non-stationary operating condition |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910276504.9A CN109934206A (en) | 2019-04-08 | 2019-04-08 | A kind of rotary machinery fault diagnosis method under non-stationary operating condition |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN109934206A true CN109934206A (en) | 2019-06-25 |
Family
ID=66989528
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201910276504.9A Pending CN109934206A (en) | 2019-04-08 | 2019-04-08 | A kind of rotary machinery fault diagnosis method under non-stationary operating condition |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN109934206A (en) |
Cited By (21)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110441060A (en) * | 2019-08-07 | 2019-11-12 | 贵州电网有限责任公司 | One kind being based on mechanical failure diagnostic method combined of multi-sensor information |
| CN110595811A (en) * | 2019-09-11 | 2019-12-20 | 浙江工业大学之江学院 | Method for constructing health state characteristic diagram of mechanical equipment |
| CN110988990A (en) * | 2019-12-10 | 2020-04-10 | 北京化工大学 | High-precision seismic attribute inversion method |
| CN111012316A (en) * | 2020-01-18 | 2020-04-17 | 四川知周光声医疗科技有限公司 | Image reconstruction system of photoacoustic mammary gland |
| CN111012317A (en) * | 2020-01-18 | 2020-04-17 | 四川知周光声医疗科技有限公司 | Photoacoustic mammary gland image reconstruction method and system |
| CN111044902A (en) * | 2019-12-31 | 2020-04-21 | 朗斯顿科技(北京)有限公司 | Motor fault diagnosis method based on current and voltage signals |
| CN111120348A (en) * | 2019-12-25 | 2020-05-08 | 中国石化销售股份有限公司华南分公司 | Centrifugal pump fault early warning method based on support vector machine probability density estimation |
| CN111289251A (en) * | 2020-02-27 | 2020-06-16 | 湖北工业大学 | Rolling bearing fine-grained fault identification method |
| CN111337263A (en) * | 2020-02-12 | 2020-06-26 | 中国民航大学 | Fault diagnosis method for engine turbine disk |
| CN111651728A (en) * | 2020-05-29 | 2020-09-11 | 中铁二院工程集团有限责任公司 | Method for identifying non-stationarity and non-stationarity of actually measured wind speed |
| CN112504434A (en) * | 2020-11-06 | 2021-03-16 | 常州大学 | System and method for measuring relative movement speed of object and air sound wave attenuation coefficient |
| CN112505641A (en) * | 2020-11-12 | 2021-03-16 | 南京理工大学 | Radar interference signal identification method based on characteristic parameter extraction |
| CN112710446A (en) * | 2020-12-21 | 2021-04-27 | 北京和中普方新能源科技有限公司 | Judgment method and system for vibration test of electric vehicle battery system and storage medium |
| CN112860183A (en) * | 2021-01-07 | 2021-05-28 | 西安交通大学 | Multisource distillation-migration mechanical fault intelligent diagnosis method based on high-order moment matching |
| CN112966632A (en) * | 2021-03-19 | 2021-06-15 | 浙江中自庆安新能源技术有限公司 | Fault identification method and system based on vibration signal imaging |
| CN113704922A (en) * | 2021-09-01 | 2021-11-26 | 江苏师范大学 | Method for predicting surface roughness based on sound vibration and texture features |
| CN113804425A (en) * | 2021-08-29 | 2021-12-17 | 西北工业大学 | Early friction instability fault identification method for sleeve gear connection structure |
| CN114165474A (en) * | 2022-02-11 | 2022-03-11 | 蘑菇物联技术（深圳）有限公司 | Method, apparatus and computer storage medium for detecting a fault condition of an air compressor |
| CN115562133A (en) * | 2022-11-10 | 2023-01-03 | 浙江大学 | Intelligent gateway of axial plunger pump |
| CN116402825A (en) * | 2023-06-09 | 2023-07-07 | 华东交通大学 | Bearing fault infrared diagnosis method, system, electronic equipment and storage medium |
| CN117076869A (en) * | 2023-10-16 | 2023-11-17 | 广东石油化工学院 | Time-frequency domain fusion fault diagnosis method and system for rotary machine |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US6253175B1 (en) * | 1998-11-30 | 2001-06-26 | International Business Machines Corporation | Wavelet-based energy binning cepstal features for automatic speech recognition |
| CN105718961A (en) * | 2016-02-15 | 2016-06-29 | 哈尔滨理工大学 | Rotating machinery intelligent fault diagnosis method based on CEEMD and image texture features |
| CN108318249A (en) * | 2018-01-24 | 2018-07-24 | 广东石油化工学院 | A kind of method for diagnosing faults of bearing in rotating machinery |
-
2019
- 2019-04-08 CN CN201910276504.9A patent/CN109934206A/en active Pending
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US6253175B1 (en) * | 1998-11-30 | 2001-06-26 | International Business Machines Corporation | Wavelet-based energy binning cepstal features for automatic speech recognition |
| CN105718961A (en) * | 2016-02-15 | 2016-06-29 | 哈尔滨理工大学 | Rotating machinery intelligent fault diagnosis method based on CEEMD and image texture features |
| CN108318249A (en) * | 2018-01-24 | 2018-07-24 | 广东石油化工学院 | A kind of method for diagnosing faults of bearing in rotating machinery |
Non-Patent Citations (1)
| Title |
|---|
| DUONG-HUNG PHAM等: ""High-Order Synchrosqueezing Transform for Multicomponent Signals Analysis—With an Application to Gravitational-Wave Signal"", 《IEEE TRANSACTIONS ON SIGNAL PROCESSING》 * |
Cited By (33)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110441060A (en) * | 2019-08-07 | 2019-11-12 | 贵州电网有限责任公司 | One kind being based on mechanical failure diagnostic method combined of multi-sensor information |
| CN110595811A (en) * | 2019-09-11 | 2019-12-20 | 浙江工业大学之江学院 | Method for constructing health state characteristic diagram of mechanical equipment |
| CN110988990A (en) * | 2019-12-10 | 2020-04-10 | 北京化工大学 | High-precision seismic attribute inversion method |
| CN111120348A (en) * | 2019-12-25 | 2020-05-08 | 中国石化销售股份有限公司华南分公司 | Centrifugal pump fault early warning method based on support vector machine probability density estimation |
| CN111044902B (en) * | 2019-12-31 | 2022-04-26 | 朗斯顿科技(北京)有限公司 | Motor fault diagnosis method based on current and voltage signals |
| CN111044902A (en) * | 2019-12-31 | 2020-04-21 | 朗斯顿科技(北京)有限公司 | Motor fault diagnosis method based on current and voltage signals |
| CN111012316A (en) * | 2020-01-18 | 2020-04-17 | 四川知周光声医疗科技有限公司 | Image reconstruction system of photoacoustic mammary gland |
| CN111012317B (en) * | 2020-01-18 | 2022-10-25 | 中川新迈科技有限公司 | Photoacoustic mammary gland image reconstruction method and system |
| CN111012316B (en) * | 2020-01-18 | 2022-10-28 | 中川新迈科技有限公司 | Image reconstruction system of photoacoustic mammary gland |
| CN111012317A (en) * | 2020-01-18 | 2020-04-17 | 四川知周光声医疗科技有限公司 | Photoacoustic mammary gland image reconstruction method and system |
| CN111337263A (en) * | 2020-02-12 | 2020-06-26 | 中国民航大学 | Fault diagnosis method for engine turbine disk |
| CN111289251A (en) * | 2020-02-27 | 2020-06-16 | 湖北工业大学 | Rolling bearing fine-grained fault identification method |
| CN111651728A (en) * | 2020-05-29 | 2020-09-11 | 中铁二院工程集团有限责任公司 | Method for identifying non-stationarity and non-stationarity of actually measured wind speed |
| CN111651728B (en) * | 2020-05-29 | 2023-03-28 | 中铁二院工程集团有限责任公司 | Method for identifying non-stationarity and non-stationarity of actually measured wind speed |
| CN112504434A (en) * | 2020-11-06 | 2021-03-16 | 常州大学 | System and method for measuring relative movement speed of object and air sound wave attenuation coefficient |
| CN112504434B (en) * | 2020-11-06 | 2022-08-26 | 常州大学 | System and method for measuring relative movement speed of object and air sound wave attenuation coefficient |
| CN112505641B (en) * | 2020-11-12 | 2024-04-05 | 南京理工大学 | Radar interference signal identification method based on characteristic parameter extraction |
| CN112505641A (en) * | 2020-11-12 | 2021-03-16 | 南京理工大学 | Radar interference signal identification method based on characteristic parameter extraction |
| CN112710446A (en) * | 2020-12-21 | 2021-04-27 | 北京和中普方新能源科技有限公司 | Judgment method and system for vibration test of electric vehicle battery system and storage medium |
| CN112860183A (en) * | 2021-01-07 | 2021-05-28 | 西安交通大学 | Multisource distillation-migration mechanical fault intelligent diagnosis method based on high-order moment matching |
| CN112966632A (en) * | 2021-03-19 | 2021-06-15 | 浙江中自庆安新能源技术有限公司 | Fault identification method and system based on vibration signal imaging |
| CN112966632B (en) * | 2021-03-19 | 2023-12-12 | 浙江中自庆安新能源技术有限公司 | Vibration signal imaging-based fault identification method and system |
| CN113804425A (en) * | 2021-08-29 | 2021-12-17 | 西北工业大学 | Early friction instability fault identification method for sleeve gear connection structure |
| CN113804425B (en) * | 2021-08-29 | 2022-06-14 | 西北工业大学 | Early friction instability fault identification method for sleeve gear connection structure |
| CN113704922B (en) * | 2021-09-01 | 2024-02-02 | 江苏师范大学 | Method for predicting surface roughness based on sound vibration and texture characteristics |
| CN113704922A (en) * | 2021-09-01 | 2021-11-26 | 江苏师范大学 | Method for predicting surface roughness based on sound vibration and texture features |
| CN114165474B (en) * | 2022-02-11 | 2022-05-06 | 蘑菇物联技术（深圳）有限公司 | Method, apparatus and computer storage medium for detecting a fault condition of an air compressor |
| CN114165474A (en) * | 2022-02-11 | 2022-03-11 | 蘑菇物联技术（深圳）有限公司 | Method, apparatus and computer storage medium for detecting a fault condition of an air compressor |
| CN115562133A (en) * | 2022-11-10 | 2023-01-03 | 浙江大学 | Intelligent gateway of axial plunger pump |
| CN116402825A (en) * | 2023-06-09 | 2023-07-07 | 华东交通大学 | Bearing fault infrared diagnosis method, system, electronic equipment and storage medium |
| CN116402825B (en) * | 2023-06-09 | 2023-08-11 | 华东交通大学 | Bearing fault infrared diagnosis method, system, electronic equipment and storage medium |
| CN117076869A (en) * | 2023-10-16 | 2023-11-17 | 广东石油化工学院 | Time-frequency domain fusion fault diagnosis method and system for rotary machine |
| CN117076869B (en) * | 2023-10-16 | 2024-01-26 | 广东石油化工学院 | Time-frequency domain fusion fault diagnosis method and system for rotary machine |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN109934206A (en) | A kind of rotary machinery fault diagnosis method under non-stationary operating condition | |
| CN106950047B (en) | The visual analysis method of vibration acceleration signal frequency spectrum | |
| CN107003663B (en) | The monitoring of device with movable part | |
| Zhao et al. | Intelligent fault diagnosis of multichannel motor–rotor system based on multimanifold deep extreme learning machine | |
| CN107483014B (en) | A kind of photovoltaic panel failure automatic detection method | |
| CN104680524B (en) | A kind of leafy vegetable disease screening method | |
| CN107368854A (en) | A kind of circuit breaker failure diagnostic method based on improvement evidence theory | |
| CN104408459A (en) | Image identification method applied to power equipment monitoring | |
| CN104596766B (en) | Early fault determining method and device for bearing | |
| KR20120027733A (en) | Rotating machinery fault diagnostic method and system using support vector machines | |
| CN110334764A (en) | Rotating machinery intelligent failure diagnosis method based on integrated depth self-encoding encoder | |
| CN112304614B (en) | Intelligent fault diagnosis method for end-to-end rolling bearing by adopting multi-attention mechanism | |
| EP3788328B1 (en) | System and process for pattern matching bearing vibration diagnostics | |
| CN110163075A (en) | A kind of multi-information fusion method for diagnosing faults based on Weight Training | |
| CN107588947B (en) | Wind turbine generator transmission chain state monitoring method considering operation condition and information simplification | |
| CN103674234A (en) | State early warning method and system for abnormal vibration of wind generating set | |
| CN109255333A (en) | A kind of large-scale wind electricity unit rolling bearing fault Hybrid approaches of diagnosis | |
| Cai et al. | A novel improved local binary pattern and its application to the fault diagnosis of diesel engine | |
| CN104778696B (en) | A kind of image border hierarchical detection method based on visual pathway azimuth sensitivity | |
| CN103984952A (en) | Method for diagnosing surface crack fault of blade of wind turbine generator of electric power system based on machine vision image | |
| CN108335019A (en) | The fault diagnosis recognition methods of fan yaw system | |
| CN107481260A (en) | A kind of region crowd is detained detection method, device and storage medium | |
| CN111353400A (en) | Whole scene vibration intensity atlas analysis method based on visual vibration measurement | |
| CN102760293B (en) | Image quality evaluation method based on distance matrix | |
| CN109918791A (en) | A kind of nuclear plant digital master control room operator human reliability analysis method |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| WD01 | Invention patent application deemed withdrawn after publication |
Application publication date: 20190625 |
|
| WD01 | Invention patent application deemed withdrawn after publication |