CN109253882A - A kind of rotor crack fault diagnostic method based on variation mode decomposition and gray level co-occurrence matrixes - Google Patents

A kind of rotor crack fault diagnostic method based on variation mode decomposition and gray level co-occurrence matrixes Download PDF

Info

Publication number
CN109253882A
CN109253882A CN201811182815.0A CN201811182815A CN109253882A CN 109253882 A CN109253882 A CN 109253882A CN 201811182815 A CN201811182815 A CN 201811182815A CN 109253882 A CN109253882 A CN 109253882A
Authority
CN
China
Prior art keywords
rotor
image
gray level
matrix
formula
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.)
Granted
Application number
CN201811182815.0A
Other languages
Chinese (zh)
Other versions
CN109253882B (en
Inventor
钟志贤
焦博隆
王家园
刘翊馨
段戬
段一戬
祁雁英
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Guilin University of Technology
Original Assignee
Guilin University of Technology
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Guilin University of Technology filed Critical Guilin University of Technology
Priority to CN201811182815.0A priority Critical patent/CN109253882B/en
Publication of CN109253882A publication Critical patent/CN109253882A/en
Application granted granted Critical
Publication of CN109253882B publication Critical patent/CN109253882B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01M—TESTING STATIC OR DYNAMIC BALANCE OF MACHINES OR STRUCTURES; TESTING OF STRUCTURES OR APPARATUS, NOT OTHERWISE PROVIDED FOR
    • G01M13/00—Testing of machine parts
    • G01M13/04—Bearings
    • G01M13/045—Acoustic or vibration analysis

Landscapes

  • Physics & Mathematics (AREA)
  • Acoustics & Sound (AREA)
  • General Physics & Mathematics (AREA)
  • Investigating Or Analyzing Materials By The Use Of Ultrasonic Waves (AREA)
  • Image Processing (AREA)

Abstract

本发明提供一种基于变分模态分解和灰度共生矩阵的转子裂纹故障诊断方法，属于裂纹故障诊断领域。包括的步骤有：采集振动信号；将采集到的信号进行变分模态分解；对每个IMF分量生成对称极坐标图像；将对称极坐标图像转化为灰度图像；灰度图生成灰度共生矩阵，提取图像纹理特征作为特征参数；选取灰度共生矩阵的特征统计量熵；求取某一工况的第i个IMF分量；分别求对应状态下的样本的在4个方向的平均熵；采集待诊断样本若干，提取出特征向量；计算马氏距离；比较待测样本的d1,d2的大小，综合距离较小者的就是待诊断样本对应的状态。本发明克服了单纯采用灰度共生矩阵在转子裂纹故障信号处理的不足，能很好完成转子裂纹故障诊断。

The invention provides a rotor crack fault diagnosis method based on variational modal decomposition and gray level co-occurrence matrix, belonging to the field of crack fault diagnosis. The steps include: collecting vibration signals; performing variational modal decomposition on the collected signals; generating a symmetrical polar coordinate image for each IMF component; converting the symmetrical polar coordinate image into a grayscale image; generating grayscale co-occurrence from the grayscale image Matrix, extract image texture features as feature parameters; select the feature statistic entropy of the gray level co-occurrence matrix; find the i-th IMF component of a certain working condition; find the average entropy of the samples in the corresponding state in four directions; Collect a number of samples to be diagnosed, extract feature vectors; calculate Mahalanobis distance; compare the sizes of d 1 and d 2 of the samples to be tested, and the one with the smaller comprehensive distance is the state corresponding to the sample to be diagnosed. The invention overcomes the deficiency of simply using the gray-scale co-occurrence matrix in rotor crack fault signal processing, and can well complete the rotor crack fault diagnosis.

Description

It is a kind of to be diagnosed based on the rotor crack fault of variation mode decomposition and gray level co-occurrence matrixes Method
[technical field]
It is specifically a kind of to be based on variation mode decomposition and gray scale symbiosis square the present invention relates to rotary machinery fault diagnosis field The method of the rotor crack fault diagnosis of battle array.
[background technique]
With industrialization paces accelerate, rotating machinery using more and more extensive.Such mechanical equipment, as space flight is started Machine, wind-driven generator, rocket engine engine etc., play important in the industrial circles such as aviation, the energy, traffic, petrochemical industry Role.
It, can be before rotor be installed life's work, by super for rotor machining process or material self-defect Technology of acoustic wave reliably detects out that whether there are cracks.The core component rotor of rotating machinery is often because in the operating conditions such as load is changeable Under, it is influenced by vertiginous thermal stress and mechanical stress and generates stress and concentrate, development quickly is tired axis crackle. If handled without discovery in time, crackle will be because of being extended intensification, mutation by stress alternation so that cause fearful disconnected Split the accident with fatal crass.
Since structure is complicated and equipment lotus root connection property is strong for above-mentioned rotating machinery, early failure of rotor vibration signal shows as non- Linearly, it non-stationary property and is easily submerged in strong noise environment, significantly impacts the extraction of rotor fault feature.In recent years, Many scholars at home and abroad are directed to the detection and diagnosis of rotor fault signal, have done a large amount of research work.
The method of currently used mode decomposition has: wavelet decomposition, empirical mode decomposition, and local mean value is decomposed.Although this A little mode decomposition methods achieve certain progress in fault diagnosis field, but its there is also shortcomings.Wavelet decomposition It comforms and selects suitable wavelet basis function and Decomposition order in multi-wavelet bases.Empirical mode decomposition does not have mathematical theory basis, Algorithm is gradually decomposed using recurrence screening, can not backward error correction, and there are modal overlap phenomenons.Local mean value Decomposition iteration It is computationally intensive, it is also easy to generate end effect.
Variation mode decomposition (Variational Mode Decomposition, VMD), is a kind of predetermined scale Nonstationary random response method, this method can be divided into unstable sophisticated signal K preset amplitude-modulation frequency-modulation signal.Variation mould State, which is decomposed, uses frequency domain non-recursive fashion Solve problems, and Decomposition Accuracy is high and can preferably avoid modal overlap phenomenon.Gray scale The textural characteristics of co-occurrence matrix energy accurate description image, the texture image identification technology based on gray level co-occurrence matrixes is in medicine figure Picture, weather nephogram are widely applied in the identification of the grain of wood.
[summary of the invention]
Unstable, non-linear for cracked rotor vibration displacement signal, fault signature difficulty is extracted, and is lacked at failure initial stage Collection and arrangement to sample, the present invention provide a kind of rotor crack fault based on variation mode decomposition and gray level co-occurrence matrixes Diagnostic method.
In order to achieve the above object, the present invention is adopted the following technical scheme that
The critical speed of rotor-support-foundation system is measured first, then makes rotor-support-foundation system in 1.5 times of critical speeds of actual speed Under operating condition, following processing are carried out:
Step 1: the radial direction on current vortex sensor and infrared electro speed probe acquisition rotor at any is respectively adopted Vibration signal and rotor speed signal, acquisition have N number of sample under cracked rotor and flawless rotor two states；
Step 2: collected signal { x (n) } being subjected to variation mode decomposition, i.e., K are sought to the vibration signal of rotor Mode function uk(t), K IMF component is obtained；
Step 3: symmetrical polar image is generated to each IMF component；
Step 4: converting gray level image for symmetrical polar image；
Step 5: grayscale image generates gray level co-occurrence matrixes, extracts image texture characteristic as characteristic parameter；
Step 6: choosing the characteristic statistic entropy of gray level co-occurrence matrixes, take e1,e2,a3,a4Respectively represent 0 °, 45 °, 90 °, The entropy in 135 ° of directions；
Step 7: setting two operating conditions of rotor-support-foundation system as j, wherein j=1 represents flawless axis, and j=2 represents crackle axis, right In the feature vector of i-th of IMF component of a certain operating condition be just Hji=[ei1,ei2,ei3,ei4], it is 4 dimensional vectors；
Step 8: and then 0 ° of N number of sample under corresponding states is asked respectively, 45 °, 90 °, the average value of the entropy in 135 ° of directions, Corresponding states feature reference vectors are obtained, are as follows:
Step 9: acquisition to diagnostic sample several, also pass through the above process, extract feature vector Hji=[ei1, ei2,ei3,ei4]
Step 10: passing through formulaCalculate separately the feature of i-th of IMF of diagnostic data Mahalanobis distance between vector and the state reference feature vector of i-th of IMF of j state；
Step 11: calculating state comprehensive distance, the comprehensive distance of sample to be tested characteristic sequence and each state reference signal, And give corresponding weighting coefficient b1, bk, calculating formula are as follows: dj=b1·dj1+b2·dj2+…+bk·djk；
Step 12: comparing the d of sample to be tested1,d2Size, determine comprehensive distance smaller is exactly to diagnostic sample pair The state answered.
In the present invention, further, in the step (2), the step of variation mode decomposition, includes:
Step 2.1: initializationN, in formula, ukTo decompose K obtained modal components, ωkFor each mould The corresponding center frequency of state component, ^ are the frequency domain representation using Parseval/Plancherel Fourier's equilong transformation, and λ is Lagrange multiplier, n are the number of iterations；
Step 2.2: updatingWith
Step 2.3: update operator
Step 2.4: step 2.2~2.4 are repeated, until meet the condition of convergence,It gives Determine discrimination precision ∈ > 0, end loop exports K IMF component.
In the present invention, further, symmetrical polar image is converted in the step 4 specific method of gray level image Are as follows: symmetrical polar image array W is set, size is m × n, and the element value of a row b column is A(a,b), symmetrical polar figure Maximum value A as inmaxIt is corresponding with tonal gradation 255, minimum value AminIt is corresponding with tonal gradation 0, gray scale is converted by following formula Value is
Unit8 is 8 rounding operation symbols of no symbol, and operation result is the integer between 0~255.
In the present invention, further, in the step 5 method particularly includes: the picture for being B from the grey level on grayscale image First position (x, y) is set out, and statistics is s at a distance from it, and the pixel position (x+Dx, y+Dy) that grey level is C occurs simultaneously Frequency P (B, C, s, θ).
P (B, C, s, θ)=[(x, y), (x+Dx, y+Dy) | f (x, y)=B, f (x+Dx, y+Dy)=C] } (19)
In formula (19), B, C=0,1,2 ..., 255 be grey level, Dx, Dy be respectively both horizontally and vertically on position Setting offset, s is the step-length of grayscale image matrix, and θ is the direction that grayscale image matrix generates, take 0 °, 45 °, 90 °, 135 ° of 4 sides To,
Pass through normalization process process by formula (20)
Wherein, g (B, C) indicates that P (B, C) matrix passes through the matrix that Regularization process obtains.
In the present invention, further, entropy is calculated as follows in the step 6:
Take s=1, e1,e2,a3,a4Respectively represent 0 °, 45 °, 90 °, the entropy in 135 ° of directions.
Due to using above-mentioned technological means, compared with the prior art, the present invention has the following beneficial effects:
1. the method for the present invention combines the advantage of variation mode decomposition and gray level co-occurrence matrixes respectively, used variation mould State decomposition method decomposes vibration signal, effectively reduces reactive component and modal overlap, and arrangement entropy is combined to calculate letter Singly, the features such as noise resisting ability is strong, from each modal components middle (center) bearing fault signature of multiple dimensioned angle extraction, overcomes simple use Gray level co-occurrence matrixes can also complete rotor crack in noise jamming in the deficiency of rotor crack fault signal processing very well Fault diagnosis.
2. the present invention using variation mode decomposition to crackle rotor fault signal carry out mode decomposition, avoid well through Test the modal overlap phenomenon occurred when mode decomposition and local mean value decomposition, end effect phenomenon.
[Detailed description of the invention]
Fig. 1 rotor crack fault diagnostic method flow chart of the present invention.
Fig. 2 is variation mode decomposition algorithm flow chart.
Fig. 3 is the schematic diagram of symmetrical polar method.
Fig. 4 is symmetrical polar figure and grayscale image corresponding relationship.
Fig. 5 is 4 generation directions of gray level co-occurrence matrixes.
Fig. 6 is flawless rotor time-domain signal.
Fig. 7 is cracked rotor time-domain signal.
[specific embodiment]
The present invention is further illustrated below by way of specific embodiment and in conjunction with attached drawing.
Based on the method for the rotor crack fault of variation mode decomposition and gray level co-occurrence matrixes diagnosis, flow chart is shown in Fig. 1 It is shown, specific steps are as follows:
Step 1: the radial direction on current vortex sensor and infrared electro speed probe acquisition rotor at any is respectively adopted Vibration signal and rotor speed signal, sample frequency 2KHz acquire 512 points, and acquisition has cracked rotor and flawless rotor two N number of sample under kind state；
Step 2: collected signal { x (n) } being subjected to variation mode decomposition, collected signal { x (n) } is passed through into change Divide mode decomposition to obtain K IMF component, is set as ui, i=1,2, K.Each IMF includes that corresponding mode is special Reference breath, then extracting feature from K IMF, so that it may show original signal xt(n) feature.
Variation mode decomposition algorithm flow chart is shown in Fig. 2, in order to which mode decomposition method of the invention is better described, illustrates below The basic principle of variation mode decomposition:
Define an amplitude-modulation frequency-modulation signal is defined as:
uk(t)=Ak(t)cos[φk(t)] (1)
In formula (1), φkIt (t) is a non-decreasing phase function, φk(t) derivative is more than or equal to zero；Signal envelope value Ak(t) it is more than or equal to zero, φk(t) speed of the speed changed much larger than the change of instantaneous frequency and envelope value.In other words, exist In one sufficiently long region in [t- δ, t+ δ], 2 π of δ ≈/φ 'k(t), uk(t) it can be used as and possess certain temporal frequency domain and amplitude Pure harmonic signal.
Its detailed process is divided into two big steps
(1) under defined frame constraint, Variation Model is constructed.
1. by each mode function uk(t) Hilbert-Huang Transform is carried out, its respective unilateral frequency spectrum analytic signal is obtained:
2. for each mode uk(t), if its centre frequency isThe two is multiplied, is adjusted in estimating accordingly Frequency of heart, to filter out " base band ":
3. finally calculating square of demodulated signal time gradient L2- norm, and utilize the analysis of Gaussian smoothing demodulated signal Method obtains the finite bandwidth of each mode function.
Then, constraint variation condition are as follows:
In formula (4): f (t) is original input signal, { uk}={ u1…uk, { ωk}={ ω1…ωkRespectively indicate by The K IMF (Intrinsic Mode Function) of VMD and corresponding centre frequency.
(2) the solution of constraint variation model problem.
For convenience of solution, the constrained function of target is changed into free objective function optimization problem, passes through introducing Secondary penalty factor α and Lagrange operational form λ (t).Wherein secondary penalty factor α improves convergence under Gauusian noise jammer and protects Signal reconstruction precision is demonstrate,proved, Lagrange operational form λ (t) allows constraint condition with more stringency.
Then augmented Lagrange multiplier formula function are as follows:
It is alternately updated by two direction of multiplication operatorλn+1Seek formula augmented Lagrange multiplier formula function As a result, being the saddle point of formula (5).This method is multiplication operator alternating direction method (alternate direction method of Multipliers, ADMM).
Then update modeIt is exactly a minimum Solve problems, expression formula are as follows:
Under L2- norm, above formula is converted to by frequency domain using Parseval/Plancherel Fourier equilong transformation:
In formula (7) first item, according to Hermitain symmetry characteristic, with ω-ωkω is substituted, formula is transformed into non-negative The form of frequency separation integral:
Obtain the update of k-th mode, that is, the solution of optimization problem are as follows:
It is also minimum optimization problem, description that similarly center frequency domain, which updates, are as follows:
Formula (10) is for conversion into the form of non-negative frequency separation integral:
Then centre frequency is in frequency domain more new-standard cement are as follows:
Operator λn+1More new-standard cement are as follows:
In formula (13)For the Lagrange multiplier operator before update；τ is step-length.
Then iteration updates uk, ωk, formula until meeting the following condition of convergence:
∈ is given judgement precision (∈ > 0), is also preset convergence error.
IMF in frequency domain is done inverse fourier transform processing, removes the data of end effect, obtains IMF points in time domain Amount.
Step 3: symmetrical polar image being generated to each IMF component, Fig. 3 is the schematic diagram of symmetrical polar method.
For each IMF component uiIn, the amplitude of t moment is Ai, the amplitude at t+L moment is At+L, substitute into formula and be converted into Polar coordinate spacePoint
R (t) is polar radius in formula (15)~(17), and θ (t) is polar coordinates counterclockwise along the angle of initial line rotation Degree,It is polar coordinates clockwise along the angle of initial line rotation, AmaxFor component uiMaximum value, AminFor component uiMinimum Value, L are time interval, and n is plane of mirror symmetry number, and θ is initial rotation angle, and g is angle magnification factor.
Step 4, gray level image is converted by symmetrical polar image, extracts image texture characteristic as characteristic parameter.
Gray level image is a two-dimensional data matrix, and the subscript of matrix element is answered with its ranks coordinate pair in the picture, The value of element represents the brightness value of corresponding position, and taking tonal gradation is 256 grades, the maximum value A in symmetrical polar imagemaxWith Tonal gradation 255 is corresponding, minimum value AminIt is corresponding with tonal gradation 0.
If claiming polar coordinate image matrix W, size is m × n, and the element value of a row b column is A(a,b),
Being converted into gray value by formula (18) is
Unit8 is 8 rounding operation symbols of no symbol, and operation result is the integer between 0~255.
It is polar diagram and grayscale image corresponding relationship such as Fig. 4.
Step 5, gray level co-occurrence matrixes are generated according to grayscale image, extracts characteristic parameter.
From the pixel position that the grey level on grayscale image is B, (x, y), statistics are s, grey level at a distance from it The frequency P (B, C, s, θ) occurred simultaneously for the pixel position (x+Dx, y+Dy) of C.
P (B, C, s, θ)=[(x, y), (x+Dx, y+Dy) | f (x, y)=B, f (x+Dx, y+Dy)=C] } (19)
In formula (19), B, C=0,1,2 ..., 255 be grey level, Dx, Dy be respectively both horizontally and vertically on position Setting offset, s is the step-length of grayscale image matrix, and θ is the direction that grayscale image matrix generates, take 0 °, 45 °, 90 °, 135 ° of 4 sides To such as 4 generation directions that Fig. 5 is gray level co-occurrence matrixes.
Pass through normalization process process by formula (20).
Step 6, the characteristic statistic entropy of gray level co-occurrence matrixes is chosen.
Entropy is calculated by formula (21)
Take s=1, e1,e2,a3,a4Respectively represent 0 °, 45 °, 90 °, the entropy in 135 ° of directions.
Step 7, right if two operating conditions of rotor-support-foundation system are j (j=1 represents flawless rotor, and j=2 represents cracked rotor) In the feature vector of i-th of IMF component of a certain operating condition be just Hji=[ei1,ei2,ei3,ei4], it is 4 dimensional vectors.
Step 8,0 ° for then seeking N number of sample under corresponding states respectively, 45 °, 90 °, the average value of the entropy in 135 ° of directions, Corresponding states feature reference vectors are obtained, are as follows:
Step 9, acquisition to diagnostic sample several, also pass through the above process, extract feature vector Hji=[ei1, ei2,ei3,ei4]
Step 10, the state of the feature vector of i-th of IMF of diagnostic data and i-th of IMF of j state is calculated separately Mahalanobis distance (Mahalanobis distance) between reference feature vector, can effectively calculate unknown sample collection similarity. Var is the variance for calculating all corresponding elements of sampling feature vectors.
Step 11, the calculating of state comprehensive distance, the synthesis of sample to be tested characteristic sequence and each state reference signal are carried out Distance, and give corresponding weighting coefficient b1, bk, such as following formula:
dj=b1·dj1+b2·dj2+…+bk·djk (24)
Step 12, finally compare the d of sample to be tested1,d2Size, determine comprehensive distance smaller be exactly sample to be diagnosed This corresponding state.
Beneficial effects of the present invention are further illustrated below with reference to specific example.
Acquisition has crackle and each 9 groups of flawless state, and first 6 groups are used as sample, and latter 3 groups are used as to diagnostic signal, and acquisition has Vibration signal under crackle state and flawless state, Fig. 6 and Fig. 7 are respectively that the time-domain signal figure of flawless rotor and crackle turn The time-domain signal figure of son, takes K=4, L=3, g=π/6, θ=0,Weighting coefficient takes b1=0.03, b2=0.03, b3=0.03, b4=0.01.
Table 1 is system mode feature reference vectors；Table 2 is the variance from each element of sample signal extracted vector；Table 3 is To diagnostic sample condition diagnosing result.
It can be obtained from table 3, the feature reference vectors of 3 flawless axis samples and two states calculate comprehensive distance, warp Compare, obtaining wherein relatively short distance is respectively 1.492,1.735,1.394, corresponding with flawless Spindle Status.Equally, 3 split The feature reference vectors of line axis sample and two states calculate comprehensive distance, and obtaining relatively short distance is respectively 1.576,1.157, 0.995, it is also corresponding with crackle Spindle Status.The experimental results showed that 6 measured signal full diagnostics successes.
In order to enhance comparison, by flawless signal and Signal of Cracks, without variation mode decomposition process, ash is directly established It spends co-occurrence matrix and extracts feature vector, carry out mahalanobis distance judgement, distance results are calculated multiplied by coefficient 0.1.As a result such as table 4.It is connect it can be seen that No. 3 sample diagnostic errors of flawless rotor, No. 2 sample diagnostic errors of cracked rotor, and each group distance all compare Closely, not sensitive enough to fault signature.So the present invention overcomes merely using gray level co-occurrence matrixes in rotor crack fault signal The deficiency of processing can complete rotor crack fault diagnosis very well.
Table 1
Table 2
Table 3
Table 4
Above description is the detailed description for the present invention preferably possible embodiments, but embodiment is not limited to this hair Bright patent claim, it is all the present invention suggested by technical spirit under completed same changes or modifications change, should all belong to In the covered the scope of the patents of the present invention.

Claims (5)

1.一种基于变分模态分解和灰度共生矩阵的转子裂纹故障诊断方法，其特征在于：首先测得转子系统的临界转速，然后使转子系统在实际转速为1.5倍临界转速的工况下，进行下述处理：1. a rotor crack fault diagnosis method based on variational modal decomposition and gray-scale co-occurrence matrix, it is characterized in that: first measure the critical speed of the rotor system, then make the rotor system at the actual speed of 1.5 times the working condition of the critical speed Next, perform the following processing: 步骤1：分别采用电涡流传感器和红外光电转速传感器采集转子上一点处的径向振动信号和转子转速信号，采集有裂纹转子和无裂纹转子两种状态下的N个样本；Step 1: Use an eddy current sensor and an infrared photoelectric speed sensor to collect the radial vibration signal and the rotor speed signal at a point on the rotor, respectively, and collect N samples in two states of the rotor with cracks and the rotor without cracks; 步骤2：将采集到的信号{x(n)}进行变分模态分解，即对转子的振动信号寻求K个模态函数uk(t)，得到K个IMF分量；Step 2: Variational modal decomposition is performed on the collected signal {x(n)}, that is, K modal functions u k (t) are sought for the vibration signal of the rotor, and K IMF components are obtained; 步骤3：对每个IMF分量生成对称极坐标图像；Step 3: Generate a symmetrical polar coordinate image for each IMF component; 步骤4：将对称极坐标图像转化为灰度图像；Step 4: Convert the symmetric polar coordinate image to a grayscale image; 步骤5：灰度图生成灰度共生矩阵，提取图像纹理特征作为特征参数；Step 5: generate a grayscale co-occurrence matrix from the grayscale image, and extract image texture features as feature parameters; 步骤6：选取灰度共生矩阵的特征统计量熵，取e1，e2，a3，a4分别代表0°，45°，90°，135°方向的熵；Step 6: Select the feature statistic entropy of the gray level co-occurrence matrix, and take e 1 , e 2 , a 3 , and a 4 to represent the entropy in the directions of 0°, 45°, 90°, and 135°, respectively; 步骤7：设转子系统的两个工况为j，其中，j＝1代表无裂纹轴，j＝2代表裂纹轴，对于某一工况的第i个IMF分量的特征向量就为Hji＝[ei1，ei2，ei3，ei4]，为4维向量；Step 7: Let the two working conditions of the rotor system be j, where j=1 represents the crack-free axis, j=2 represents the cracked axis, and the eigenvector of the i-th IMF component for a certain working condition is H ji = [e i1 , e i2 , e i3 , e i4 ], is a 4-dimensional vector; 步骤8：然后分别求对应状态下的N个样本的0°，45°，90°，135°方向的熵的平均值，得到对应状态特征参考向量，为： Step 8: Then calculate the average value of the entropy in the 0°, 45°, 90°, and 135° directions of the N samples in the corresponding state, and obtain the corresponding state feature reference vector, which is: 步骤9：采集待诊断样本若干个，同样经过上述过程，提取出特征向量Hji＝[ei1，ei2，ei3，ei4]Step 9: Collect several samples to be diagnosed, and extract the feature vector H ji =[e i1 , e i2 , e i3 , e i4 ] through the same process as above 步骤10：通过式子分别计算待诊断数据的第i个IMF的特征向量与j状态的第i个IMF的状态参考特征向量之间的马氏距离；Step 10: Pass the formula Calculate the Mahalanobis distance between the eigenvector of the ith IMF of the data to be diagnosed and the state reference eigenvector of the ith IMF of the j state; 步骤11：计算状态综合距离，待测样本特征序列与各个状态参考信号的综合距离，并给予相应的加权系数b1，…，bk，计算式为：dj＝b1·dj1+b2·dj2+…+bk·djk；Step 11: Calculate the state comprehensive distance, the comprehensive distance between the feature sequence of the sample to be tested and each state reference signal, and give the corresponding weighting coefficients b 1 , . . . , b k , the calculation formula is: d j =b 1 ·d j1 +b 2 d j2 +…+b k d jk ; 步骤12：比较待测样本的d1，d2的大小，确定综合距离较小者的就是待诊断样本对应的状态。Step 12: Compare the sizes of d 1 and d 2 of the sample to be tested, and determine the state corresponding to the sample to be diagnosed with the smaller comprehensive distance. 2.根据权利要求1所述的故障诊断方法，其特征在于：所述步骤(2)中，变分模态分解的步骤包括：2. The fault diagnosis method according to claim 1, wherein in the step (2), the step of variational mode decomposition comprises: 步骤2.1：初始化式中，uk为分解得到的K个模态分量，ωk为各模态分量对应的频率中心，^为利用Parseval/Plancherel傅里叶等距变换的频域表示，λ为Lagrange乘子，n为迭代次数；Step 2.1: Initialization In the formula, u k is the K modal components obtained by decomposition, ω k is the frequency center corresponding to each modal component, ^ is the frequency domain representation using the Parseval/Plancherel Fourier equidistant transform, λ is the Lagrange multiplier, n is the number of iterations; 步骤2.2：更新和 Step 2.2: Update and 步骤2.3：更新算子 Step 2.3: Update the operator 步骤2.4：重复步骤2.2～2.4，直至满足收敛条件，给定判别精度∈＞0，结束循环，输出K个IMF分量。Step 2.4: Repeat steps 2.2 to 2.4 until the convergence conditions are met, Given the discrimination accuracy ∈>0, end the loop and output K IMF components. 3.根据权利要求2所述的故障诊断方法，其特征在于：所述步骤4中将对称极坐标图像转化为灰度图像的具体方法为：设对称极坐标图像矩阵W，尺寸大小为m×n，第a行第b列的元素值为A(a，b)，对称极坐标图像中的最大值Amax与灰度等级255对应，最小值Amin与灰度等级0对应，经过下式转化为灰度值为3. The fault diagnosis method according to claim 2, wherein the specific method for converting the symmetrical polar coordinate image into a grayscale image in the step 4 is: set the symmetrical polar coordinate image matrix W, and the size is m× n, the element value of row a and column b is A (a, b) , the maximum value A max in the symmetric polar coordinate image corresponds to gray level 255, and the minimum value A min corresponds to gray level 0. After the following formula Converted to grayscale value unit8为无符号8位取整运算符，运算结果为0～255之间的整数。unit8 is an unsigned 8-bit integer operator, and the result of the operation is an integer between 0 and 255. 4.根据权利要求1所述的故障诊断方法，其特征在于：所述步骤5中的具体方法为：从灰度图上的灰度级别为B的像元位置(x，y)出发，统计与其的距离为s，灰度级别为C的像元位置(x+Dx，y+Dy)同时出现的频度P(B，C，s，θ)。4 . The fault diagnosis method according to claim 1 , wherein the specific method in the step 5 is: starting from the pixel position (x, y) with the gray level B on the grayscale map, statistical The frequency P(B, C, s, θ) of the simultaneous occurrence of the pixel position (x+Dx, y+Dy) with a distance of s and a gray level of C. P(B，C，s，θ)＝{[(x，y)，(x+Dx，y+Dy)|f(x，y)＝B，f(x+Dx，y+Dy)＝C]} (19)P(B,C,s,θ)={[(x,y),(x+Dx,y+Dy)|f(x,y)=B,f(x+Dx,y+Dy)=C ]} (19) 式(19)中，B，C＝0，1，2，...，255为灰度级别，Dx，Dy分别为水平和垂直方向上的位置偏移量，s为灰度图矩阵的步长，θ为灰度图矩阵生成的方向，取0°，45°，90°，135°4个方向，In formula (19), B, C=0, 1, 2, ..., 255 are the gray levels, Dx, Dy are the position offsets in the horizontal and vertical directions, respectively, and s is the step of the gray image matrix. Length, θ is the direction generated by the grayscale image matrix, taking 4 directions of 0°, 45°, 90°, 135°, 按式(20)经过正规化处理过程According to the formula (20), the normalization process is carried out 其中，g(B，C)表示P(B，C)矩阵经过正则化处理过程得到的矩阵。Among them, g(B, C) represents the matrix obtained by the regularization process of the P(B, C) matrix. 5.根据权利要求1所述的故障诊断方法，其特征在于：所述步骤6按下式计算熵：5. fault diagnosis method according to claim 1, is characterized in that: described step 6 calculates entropy as follows: 取s＝1，e1，e2，a3，a4分别代表0°，45°，90°，135°方向的熵。Take s=1, e 1 , e 2 , a 3 , a 4 represent the entropy in the directions of 0°, 45°, 90°, and 135°, respectively.
CN201811182815.0A 2018-10-08 2018-10-08 A Rotor Crack Fault Diagnosis Method Based on Variational Mode Decomposition and Gray Co-occurrence Matrix Active CN109253882B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201811182815.0A CN109253882B (en) 2018-10-08 2018-10-08 A Rotor Crack Fault Diagnosis Method Based on Variational Mode Decomposition and Gray Co-occurrence Matrix

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201811182815.0A CN109253882B (en) 2018-10-08 2018-10-08 A Rotor Crack Fault Diagnosis Method Based on Variational Mode Decomposition and Gray Co-occurrence Matrix

Publications (2)

Publication Number Publication Date
CN109253882A true CN109253882A (en) 2019-01-22
CN109253882B CN109253882B (en) 2020-08-07

Family

ID=65045243

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201811182815.0A Active CN109253882B (en) 2018-10-08 2018-10-08 A Rotor Crack Fault Diagnosis Method Based on Variational Mode Decomposition and Gray Co-occurrence Matrix

Country Status (1)

Country Link
CN (1) CN109253882B (en)

Cited By (16)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109682561A (en) * 2019-02-19 2019-04-26 大连理工大学 A kind of automatic detection high-speed railway bridge free vibration responds the method to identify mode
CN109947047A (en) * 2019-03-28 2019-06-28 西安科技大学 A method for diagnosing unbalance faults of electric spindle
CN110134110A (en) * 2019-05-15 2019-08-16 哈尔滨工业大学 Rotor Crack Fault Detection Method Based on Interval Control Strategy
CN110595780A (en) * 2019-09-20 2019-12-20 西安科技大学 Bearing Fault Recognition Method Based on Vibration Gray Image and Convolutional Neural Network
CN111325722A (en) * 2020-02-17 2020-06-23 江苏诚印科技有限公司 Stamp image accurate identification method, stamp image identification processing method and stamp image identification system
CN111487046A (en) * 2020-02-27 2020-08-04 广西电网有限责任公司电力科学研究院 Fault diagnosis method for circuit breaker voiceprint and vibration entropy feature fusion
NL2024358B1 (en) * 2019-06-03 2020-12-08 Univ Anhui Sci & Technology Method for quantitatively evaluating dynamic quality of rolling bearing based on permutation entropy
CN112414715A (en) * 2020-11-05 2021-02-26 西安工程大学 Bearing fault diagnosis method based on mixed feature and improved gray level co-occurrence algorithm
CN112763904A (en) * 2020-12-29 2021-05-07 广州航天海特系统工程有限公司 Circuit breaker detection method, device, equipment and storage medium
CN112907524A (en) * 2021-02-04 2021-06-04 哈尔滨市科佳通用机电股份有限公司 Method for detecting fault of fire-proof plate of rail wagon based on image processing
CN113345399A (en) * 2021-04-30 2021-09-03 桂林理工大学 Method for monitoring sound of machine equipment in strong noise environment
CN113514246A (en) * 2021-04-23 2021-10-19 河北科技大学 Rotating mechanical system damage detection method, device and terminal
CN113984192A (en) * 2021-10-27 2022-01-28 广东电网有限责任公司佛山供电局 Transformer working state monitoring method and system based on sound signals
CN114548186A (en) * 2022-03-02 2022-05-27 浙江浙能技术研究院有限公司 A method for diagnosing faults of rotors of high-speed gearboxes of air compressors
CN114693650A (en) * 2022-03-31 2022-07-01 南通俊朗智能科技有限公司 Intelligent control method of mixing machine based on machine vision
CN117274117A (en) * 2023-11-23 2023-12-22 合肥工业大学 Frequency domain pseudo-color enhanced cardiomagnetic signal characteristic image generation method and storage medium

Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101593274A (en) * 2009-07-02 2009-12-02 浙江省电力公司 Texture-based Feature Extraction Method for Transmission Line Equipment
CN204594692U (en) * 2015-03-24 2015-08-26 桂林理工大学 The synchronous acoustical signal harvester of rolling bearing fault diagnosis experiment
CN104964821A (en) * 2015-05-22 2015-10-07 南京航空航天大学 Fault detection method and fault detection apparatus used for shafting device
CN105987809A (en) * 2015-02-10 2016-10-05 沈阳透平机械股份有限公司 Centrifugal-compressor semi-open-type impeller crack detection method based on random resonance
KR20170109872A (en) * 2016-03-22 2017-10-10 두산중공업 주식회사 Inspection method for a rotating component
US20180128711A1 (en) * 2012-10-26 2018-05-10 Acellent Technologies, Inc. System and method for monitoring the structural health of coupled bearings
CN108225762A (en) * 2016-12-15 2018-06-29 唐智科技湖南发展有限公司 A kind of single gear tooth crackle broken teeth fault identification diagnostic method
CN108414226A (en) * 2017-12-25 2018-08-17 哈尔滨理工大学 Fault Diagnosis of Roller Bearings under the variable working condition of feature based transfer learning

Patent Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101593274A (en) * 2009-07-02 2009-12-02 浙江省电力公司 Texture-based Feature Extraction Method for Transmission Line Equipment
US20180128711A1 (en) * 2012-10-26 2018-05-10 Acellent Technologies, Inc. System and method for monitoring the structural health of coupled bearings
CN105987809A (en) * 2015-02-10 2016-10-05 沈阳透平机械股份有限公司 Centrifugal-compressor semi-open-type impeller crack detection method based on random resonance
CN204594692U (en) * 2015-03-24 2015-08-26 桂林理工大学 The synchronous acoustical signal harvester of rolling bearing fault diagnosis experiment
CN104964821A (en) * 2015-05-22 2015-10-07 南京航空航天大学 Fault detection method and fault detection apparatus used for shafting device
KR20170109872A (en) * 2016-03-22 2017-10-10 두산중공업 주식회사 Inspection method for a rotating component
CN108225762A (en) * 2016-12-15 2018-06-29 唐智科技湖南发展有限公司 A kind of single gear tooth crackle broken teeth fault identification diagnostic method
CN108414226A (en) * 2017-12-25 2018-08-17 哈尔滨理工大学 Fault Diagnosis of Roller Bearings under the variable working condition of feature based transfer learning

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
刘军: ""发动机非线性转子裂纹检测与延缓扩展的研究"", 《计算机仿真》 *
钟志贤: ""主动电磁轴承-裂纹柔性转子系统动力学特性研究"", 《中国博士学位论文全文数据库》 *

Cited By (23)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109682561A (en) * 2019-02-19 2019-04-26 大连理工大学 A kind of automatic detection high-speed railway bridge free vibration responds the method to identify mode
CN109947047A (en) * 2019-03-28 2019-06-28 西安科技大学 A method for diagnosing unbalance faults of electric spindle
CN109947047B (en) * 2019-03-28 2021-03-23 西安科技大学 A method for diagnosing unbalance faults of electric spindle
CN110134110A (en) * 2019-05-15 2019-08-16 哈尔滨工业大学 Rotor Crack Fault Detection Method Based on Interval Control Strategy
NL2024358B1 (en) * 2019-06-03 2020-12-08 Univ Anhui Sci & Technology Method for quantitatively evaluating dynamic quality of rolling bearing based on permutation entropy
CN110595780A (en) * 2019-09-20 2019-12-20 西安科技大学 Bearing Fault Recognition Method Based on Vibration Gray Image and Convolutional Neural Network
CN111325722A (en) * 2020-02-17 2020-06-23 江苏诚印科技有限公司 Stamp image accurate identification method, stamp image identification processing method and stamp image identification system
CN111325722B (en) * 2020-02-17 2024-02-20 江苏诚印科技有限公司 Seal image accurate identification method and system and seal image identification processing method
CN111487046A (en) * 2020-02-27 2020-08-04 广西电网有限责任公司电力科学研究院 Fault diagnosis method for circuit breaker voiceprint and vibration entropy feature fusion
CN111487046B (en) * 2020-02-27 2021-10-26 广西电网有限责任公司电力科学研究院 Fault diagnosis method for circuit breaker voiceprint and vibration entropy feature fusion
CN112414715A (en) * 2020-11-05 2021-02-26 西安工程大学 Bearing fault diagnosis method based on mixed feature and improved gray level co-occurrence algorithm
CN112414715B (en) * 2020-11-05 2022-09-27 西安工程大学 Bearing fault diagnosis method based on mixed feature and improved gray level symbiosis algorithm
CN112763904A (en) * 2020-12-29 2021-05-07 广州航天海特系统工程有限公司 Circuit breaker detection method, device, equipment and storage medium
CN112907524A (en) * 2021-02-04 2021-06-04 哈尔滨市科佳通用机电股份有限公司 Method for detecting fault of fire-proof plate of rail wagon based on image processing
CN113514246A (en) * 2021-04-23 2021-10-19 河北科技大学 Rotating mechanical system damage detection method, device and terminal
CN113345399A (en) * 2021-04-30 2021-09-03 桂林理工大学 Method for monitoring sound of machine equipment in strong noise environment
CN113984192A (en) * 2021-10-27 2022-01-28 广东电网有限责任公司佛山供电局 Transformer working state monitoring method and system based on sound signals
CN113984192B (en) * 2021-10-27 2023-08-01 广东电网有限责任公司佛山供电局 Transformer working state monitoring method and system based on sound signals
CN114548186A (en) * 2022-03-02 2022-05-27 浙江浙能技术研究院有限公司 A method for diagnosing faults of rotors of high-speed gearboxes of air compressors
CN114693650A (en) * 2022-03-31 2022-07-01 南通俊朗智能科技有限公司 Intelligent control method of mixing machine based on machine vision
CN114693650B (en) * 2022-03-31 2024-05-24 无锡市慧航智能科技有限公司 Intelligent control method of mixer based on machine vision
CN117274117A (en) * 2023-11-23 2023-12-22 合肥工业大学 Frequency domain pseudo-color enhanced cardiomagnetic signal characteristic image generation method and storage medium
CN117274117B (en) * 2023-11-23 2024-02-02 合肥工业大学 Frequency domain pseudo-color enhanced magnetocardiogram signal characteristic image generation method and storage medium

Also Published As

Publication number Publication date
CN109253882B (en) 2020-08-07

Similar Documents

Publication Publication Date Title
CN109253882A (en) A kind of rotor crack fault diagnostic method based on variation mode decomposition and gray level co-occurrence matrixes
Chen et al. A meta-learning method for electric machine bearing fault diagnosis under varying working conditions with limited data
Xu et al. Fan fault diagnosis based on symmetrized dot pattern analysis and image matching
CN112304613A (en) Wind turbine generator bearing early warning method based on feature fusion
CN111238815B (en) Bearing fault identification method based on data enhancement under sample imbalance
CN107247231A (en) A kind of aerogenerator fault signature extracting method based on OBLGWO DBN models
CN111256993A (en) Method and system for diagnosing fault type of main bearing of wind turbine generator
CN110131109A (en) A wind turbine blade imbalance detection method based on convolutional neural network
CN112285554A (en) Information fusion-based demagnetization fault diagnosis method and device for permanent magnet synchronous motor
Huang et al. A Fault Diagnosis Approach for Rolling Bearing Based on Wavelet Packet Decomposition and GMM-HMM.
CN108444696A (en) A kind of gearbox fault analysis method
CN111412114B (en) Wind turbine generator impeller imbalance detection method based on stator current envelope spectrum
CN114549589B (en) Rotator vibration displacement measurement method and system based on lightweight neural network
CN108919116A (en) Ocean current generator imbalance stator current method for diagnosing faults based on MCCKAF-FFT-Softmax
Chen et al. An adversarial learning framework for zero-shot fault recognition of mechanical systems
CN117590989A (en) Motor rotating speed online estimation device and method based on neural network
CN109633270A (en) The identification of fault energy zone boundary and feature extracting method based on instantaneous spectrum entropy and noise energy difference
CN106053090A (en) Neighbor abnormal detection system of gas turbine
Cao et al. A novel method for detection of wind turbine blade imbalance based on multi-variable spectrum imaging and convolutional neural network
Yu et al. Threshing cylinder unbalance detection using a signal extraction method based on parameter-adaptive variational mode decomposition
Kharyton et al. Advanced processing of a blade vibratory response obtained with tip timing method using hyperparameter-free sparse estimation method
CN117951541A (en) A rolling bearing fault diagnosis method and system based on grey correlation
CN114781427B (en) Wind generating set antifriction bearing fault diagnosis system based on information fusion technology
Li et al. A signal based “W” structural elements for multi-scale mathematical morphology analysis and application to fault diagnosis of rolling bearings of wind turbines
CN115597868A (en) Bearing fault classification method based on wavelet signal processing and deep learning

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
GR01 Patent grant
GR01 Patent grant