CN104122086B - Fault Diagnosis of Gearbox for Wind Turbine method based on Kriging model - Google Patents

Abstract

The invention discloses a kind of Fault Diagnosis of Gearbox for Wind Turbine method based on Kriging model, its step is as follows: gather multiple vibration signal sequences under each operating mode pattern of wind energy conversion system gear-box and calculate multiple fault-signal characteristic values; Set up corresponding relation and sample data table between operating mode pattern, diagnosis desired value; Set up variogram theoretical model according to sample data table; Build Kriging model based on variogram theoretical model; Detect the fault-signal characteristic value of vibration signal to be diagnosed and input Kriging model and obtain Kriging estimator, determine the operating mode pattern of vibration signal to be diagnosed according to the corresponding relation between Kriging estimator inquiry operating mode pattern, diagnosis desired value. The present invention have advantages of diagnostic result quick and precisely, nonlinear fitting is effective, use is flexible, amount of calculation is little, can lay the foundation for realizing the inline diagnosis of wind energy conversion system gearbox fault.

Description

Fault Diagnosis of Gearbox for Wind Turbine method based on Kriging model

Technical field

The present invention relates to the fault diagnosis technology of wind energy conversion system gear-box, be specifically related to a kind of wind energy conversion system gear based on Kriging modelBox fault diagnosis method.

Background technology

At present wind-powered electricity generation unit is main mainly with trunnion axis, comprise Stall Type, become oar type, variable-pitch variable-speed type, direct-driving type, permanent-magnet type,Half direct-driving types etc. are several, and adopting the megawatt-level wind unit of the double-fed asynchronous variable-speed generator of gear-box speedup is current wind-power marketMainstream model, wherein gear-box is made up of gear-box body, gear, power transmission shaft, bearing and bearing block. Along with China's wind-powered electricity generation dressThe put into operation accumulation of time of the sharp increase of machine capacity and wind energy conversion system, the unit outage event being caused by gearbox fault or damage happens occasionally.According to statistics, China's wind field gear-box spoilage is up to 40～50%, and individual brand set gear box turnover rate almost reaches100%. The direct and consequential damage bringing is thus also increasing, and the workload that attendant drops into related work also constantly rises.Taking 1.5MW large scale wind power machine gearbox fault as example, the dismounting of fault gear-box, transport and maintenance cost up to 1,000,000 yuan withUpper, be equivalent to 10% of wind energy conversion system total production cost, and this cost not yet comprises the production loss of compressor emergency shutdown several months. CauseThis, gear-box damages the safe operation of wind-powered electricity generation unit equipment and wind-powered electricity generation enterprise is constituted a serious threat.

Status monitoring and fault diagnosis are to reduce effective ways and the hand of wind energy conversion system gear-box and the operation of wind-powered electricity generation unit and maintenance costSection. The research of Fault Diagnosis of Gearbox for Wind Turbine method and application, can find incipient fault risk as early as possible, prevents degree of injuryDeterioration or the expansion of fault coverage, avoid contingency and casualties; Reduce breakdown maintenance number of times and downtime, increaseEquipment normal operating time and degree of reliability, improve utilization rate of equipment and installations; Reduce unnecessary part replacement, to reduce maintenance costWith downtime, save the expense of a large amount of operation and maintenances. Data shows, utilizes condition monitoring and fault diagnosis technology, canReduce accident rate 75%, reduce maintenance cost 25%～50%, make a profit and can reach 17:1 with output investment ratio.

Method for diagnosing faults is generally divided into traditional diagnosis method, Mathematical Diagnosis method and intelligent failure diagnosis method. At present, wind-powered electricity generationThe most frequently used method of set gear box fault diagnosis is that the vibration signal to collecting carries out spectrum analysis, is undertaken by artificial or machinePattern-recognition, thus judge whether to break down and failure cause, belong to the combination of traditional diagnosis method and Mathematical Diagnosis method.Along with the development of artificial intelligence, many intelligent diagnosing methods start to be applied to as fuzzy logic, expert system and neutral net etc.The fault diagnosis of wind power equipment system. But, because the relation between gearbox of wind turbine failure symptom and reason exists not trueQualitative, and existing various intelligent diagnosing methods are perfect not enough at present, cause current existing intelligent diagnosing method not only to be countedCalculation amount is larger, and real-time is poor, and fault False Rate is higher, is difficult to meet actual demands of engineering.

Kriging (Ke Lijin) method comes from a kind of interpolation method of South Africa geologist Krige Prospecting for Gold at first, through methodThe one optimum that the famous statistician Matheron of state theorizes and systematization grows up afterwards is gradually without bias estimation. The partyMethod is, in finite region, the value of regional change amount is carried out to nothing optimal estimation partially, limited by sample point to be estimated (or region)Some sample points certificates of having measured in neighborhood, consider the shape, size in sample spot (or region) and with point to be estimated (orRegion) the space structure information of the geometric properties such as spatial distribution position each other and variable, give one to each sampleFixed weight coefficient, estimates thereby adopt weighted mean method to treat the unknown quantity of estimating point (or region). Kriging method is initialBe applied to mining industry field, its range of application is well beyond original mining industry field, at meteorologic analysis, environment prison in recent yearsThe effect becoming more and more important is all being brought into play in the fields such as survey, image processing, mechanical engineering, structure optimization and Aero-Space.

Summary of the invention

The technical problem to be solved in the present invention is: for the above-mentioned technical problem of prior art, provide a kind of diagnostic result accurate fastReally, nonlinear fitting is effective, use is flexible, amount of calculation is little, can establish for realizing the inline diagnosis of wind energy conversion system gearbox faultThe Fault Diagnosis of Gearbox for Wind Turbine method based on Kriging model on basis.

In order to solve the problems of the technologies described above, the technical solution used in the present invention is:

A Fault Diagnosis of Gearbox for Wind Turbine method based on Kriging model, implementation step is as follows:

1) gather respectively the multiple vibration signal sequences of wind energy conversion system gear-box under each operating mode pattern, described each operating mode pattern bagDraw together normal condition, state of wear and the broken conditions of wind energy conversion system gear-box;

2) calculate respectively multiple fault-signal characteristic values of described vibration signal sequence;

3) according to fault-signal characteristic value corresponding to multiple vibration signal sequences under each operating mode pattern, set up by membership functionCorresponding relation and sample data table between operating mode pattern, diagnosis desired value, each list item of described sample data table comprises workThe data of condition pattern, multiple fault-signal characteristic value, diagnosis desired value three types;

4), according to the some or all of data in described sample data table, adopt least square method to carry out variogram theoretical modelMatching, obtain parameter value and the variogram theoretical model of variogram theoretical model;

5), based on described variogram theoretical model, set up according to Kriging method the input that each fault-signal characteristic value formsApproximate function between vector, diagnosis desired value obtains the Kriging mould being made up of the calculation expression of Kriging estimatorType, the input quantity of described Kriging model be the input that forms of each fault-signal characteristic value that vibration signal to be diagnosed is corresponding toAmount, the output quantity of described Kriging model is the Kriging estimator corresponding with diagnosis desired value;

6) detect the wind energy conversion system gear-box vibration signal to be diagnosed under in working order, each of vibration signal to be diagnosed described in calculatingFault-signal characteristic value, described in the input vector that each fault-signal characteristic value of vibration signal described to be diagnosed is formed is inputtedKriging model, obtains the Kriging estimator that Kriging model is exported, and inquires about described work according to described Kriging estimatorCorresponding relation between condition pattern, diagnosis desired value, thus determine operating mode pattern corresponding to vibration signal to be diagnosed.

Preferably, described multiple fault-signal characteristic value comprises Power Spectral Entropy, wavelet packet energy entropy, kurtosis, the degree of bias, correlation dimensionAnd box counting dimension.

Preferably, described step 2) in specifically refer to the Power Spectral Entropy of calculating vibration signal sequence according to formula (1);

$G_F=-\underset{i=1}{\overset{N}{\Σ}}{q}_i{\mathrm{logq}}_i---\left(1\right)$

In formula (1), G_FRepresent the Power Spectral Entropy of vibration signal sequence, q_iRepresent that the power spectrum of i vibration signal is in vibrationShared proportion in summed power spectrum in burst, N represents the quantity of vibration signal in vibration signal sequence.

Preferably, described step 2) in specifically refer to the wavelet packet energy entropy that calculates vibration signal sequence according to formula (2);

$H_F=-\underset{i=1}{\overset{N}{\Σ}}{E}_i{\mathrm{logE}}_i---\left(2\right)$

In formula (2), H_FRepresent the wavelet packet energy entropy of vibration signal sequence, E_iThe frequency band energy that represents i vibration signal existsShared proportion in frequency band energy summation in vibration signal sequence, N represents the quantity of vibration signal in vibration signal sequence.

Preferably, described step 2) in specifically refer to the kurtosis of calculating vibration signal sequence according to formula (3);

$K_q=\frac{1}{N}\underset{i=1}{\overset{N}{\Σ}}\frac{x_i^4}{{\left(\left({\mathrm{\Σ}}_{i=1}^N{x}_i^2\right)/N\right)}^2}-3---\left(3\right)$

In formula (3), K_qRepresent the kurtosis of vibration signal sequence, { x_i(1≤i≤N) expression vibration signal sequence, N represents vibrationThe quantity of vibration signal in burst.

Preferably, described step 2) in specifically refer to the degree of bias of calculating vibration signal sequence according to formula (4);

$K_p=\frac{1}{N}\frac{\left({\mathrm{\Σ}}_{i=1}^N\right|x_i^3\left|\right)/N}{{\left({\mathrm{\Σ}}_{i=1}^N{x}_i^2\right)}^{1.5}}---\left(4\right)$

In formula (4), K_pRepresent the degree of bias of vibration signal sequence, { x_i(1≤i≤N) expression vibration signal sequence, N represents vibrationThe quantity of vibration signal in burst.

Preferably, described step 2) in to calculate the detailed step of correlation dimension of vibration signal sequence as follows:

2.1) by vibration signal sequence { x_i(i=1,2 ..., n) carry out phase space reconfiguration and obtain nm vector y_iThe vector sequence of composition{y_i}(i＝1,2,…,nm)；

2.2) calculate described vector sequence { y_i(i=1,2 ..., nm) in arbitrary reference point to the distance of all the other points;

2.3) respectively by described vector sequence { y_i(i=1,2 ..., nm) in any two reference points as current vector pair, calculate work asThe front right distance of vector, compares distance right current vector and measurement vector to related given parameters r, asThis distance of fruit is less than given parameters r and judges that current vector is to being incidence vector pair, otherwise judges that current vector is to being non-passConnection vector pair;

2.4) calculate the right correlation integral of each incidence vector according to formula (5);

$C_m\left(r\right)=\frac{1}{{\mathrm{nm}}^2}\underset{i,j=0}{\overset{nm}{\Σ}}H(r-r_{ij})---\left(5\right)$

In formula (5), C_m(r) represent the described right correlation integral of each incidence vector, nm represents vector sequence { y_i}(i＝1,2,…,nm)Middle vector y_iQuantity, r represents to weigh vector to related given parameters, r_ijRepresent vector sequence { y_i(i=1,2 ..., nm) inAs right i the reference point of current vector, to the distance between j reference point, H (x) is jump function;

2.5) according to the correlation dimension of formula (6) calculating vibration signal sequence;

$D_2\left(m\right)=\underset{\mathrm{\ϵ}\→0}{\lim }\frac{{\mathrm{lnC}}_m\left(r\right)}{\ln \left(r\right)}---\left(6\right)$

In formula (6), D₂(m) correlation dimension of expression vibration signal sequence, ε represents the little square box covering for fractalThe length of side, C_m(r) represent the correlation integral that each incidence vector is right, r represents to weigh vector to related given parameters.

Preferably, described step 2) in specifically refer to the box counting dimension that calculates vibration signal sequence according to formula (7);

$D_0=\underset{\mathrm{\ϵ}\→0}{\lim }\frac{\ln N\left(\mathrm{\ϵ}\right)}{\ln N(1/\mathrm{\ϵ})}---\left(7\right)$

In formula (7), D₀Represent the box counting dimension of vibration signal sequence, ε represents the limit of the little square box covering for fractalLong, 1/ ε represents the inverse of ε, and N () represents that utilizing the little square box of the corresponding length of side to carry out signal covers required box number.

Preferably, described step 4) detailed step as follows:

4.1) select the some or all of list item in described sample data table, by the multiple fault-signal features in each list itemValue is common forms an input vector as a sample point X, observation station in the variogram survey region shown in use formula (8)Between spatial coherence, corresponding variogram value γ (d) when obtaining Euclidean distance d value and changing;

$\mathrm{\γ}\left(d\right)=\frac{1}{2N\left(d\right)}\underset{r=1}{\overset{N\left(d\right)}{\Σ}}{\[y\left(X_i^r\right)-y\left(X_j^r\right)\]}^2---\left(8\right)$

In formula (8), γ (d) represents variogram value, and X represents sample point, and y (X) is the response observation of sample point X, and d is sampleThis between Euclidean distance,WithFor any two sample points in Euclidean distance d, N (d) is sample point pairNumber;

4.2) when Euclidean distance d value variation, corresponding variogram value γ (d), as input, adopts formula according to least square method(9) Gauss model shown in carries out matching, the variogram theory that obtains the parameter value of Gauss model and be made up of this Gauss modelModel;

$\mathrm{\γ}\left(d\right)=\left\{\begin{array}{cc}0& d=0\\ c_0+c\left(1-e^{-\frac{d^2}{a^2}}\right)& d>0\end{array}\right.---\left(9\right)$

In formula (9), d represents Euclidean distance, and γ (d) represents variogram value corresponding to Euclidean distance d, c₀, c, a be respectivelyThe parameter value of Gauss model, wherein c₀The piece gold that represents Gauss model, c represents the sagitta of Gauss model,Represent Gauss modelRange.

Preferably, described step 5) detailed step as follows:

5.1), for the sample data of given research object, utilize formula (10) and formula (11) to carry out the power of calculative determination observationHeavy coefficient w;

$\left\{\begin{array}{c}\mathrm{\Γ}={\left\{\mathrm{\γ}(x^{\left(i\right)}-x^{\left(j\right)})\right\}}_{ij},\\ \mathrm{\γ}=\{\mathrm{\γ}(x^{\left(1\right)}-x),...,\mathrm{\γ}(x^{\left(n\right)}-x)\}\·\end{array}\right.---\left(10\right)$

$w={\left(\mathrm{\γ}+I\frac{1-I^T{\mathrm{\Γ}}^{-1}\mathrm{\γ}}{I^T{\mathrm{\Γ}}^{-1}I}\right)}^T{\mathrm{\Γ}}^{-1}---\left(11\right)$

In formula (10) and formula (11), γ represents the variogram value of utilizing formula (8) to calculate, and Γ represents by variogramThe matrix that calculated value γ forms as element, x⁽ⁱ⁾And x^(j)Represent any two sample points of sample data, w represents observationWeight coefficient, n is sample point number used in calculating, I represents n dimension unit column vector, the transposition computing of T representing matrix;

5.2), according to the weight coefficient w of described observation, set up according to Kriging method defeated that each fault-signal characteristic value formsIncoming vector, the approximate function of diagnosis between desired value, obtain the calculation expression structure of Kriging estimator shown in formula (12)The Kriging model becoming, the input quantity of described Kriging model is each fault-signal characteristic value that vibration signal to be diagnosed is correspondingThe input vector forming, the output quantity of described Kriging model is the Kriging estimator corresponding with diagnosis desired value;

$\hat{z}\left(x\right)=\underset{i=1}{\overset{n}{\Σ}}{w}_{i}z\left(x^{\left(i\right)}\right),---\left(12\right)$

In formula (12),Represent the output quantity of Kriging model, w_iRepresent i the component of weight coefficient w, z (x⁽ⁱ⁾)Represent the input quantity of Kriging model, sample point number used during n represents to calculate.

The Fault Diagnosis of Gearbox for Wind Turbine method that the present invention is based on Kriging model has following advantage: the present invention is by gathering windThe multiple vibration signal sequences of power machine gear-box under each operating mode pattern are also calculated multiple fault-signal characteristic values; Set up operating mode mouldCorresponding relation and sample data table between formula, diagnosis desired value; Set up variogram theoretical model according to sample data table; BaseBuild Kriging model in variogram theoretical model; Detect the fault-signal characteristic value input Kriging of vibration signal to be diagnosedModel obtains Kriging estimator, true according to the corresponding relation between Kriging estimator inquiry operating mode pattern, diagnosis desired valueThe operating mode pattern of vibration signal surely to be diagnosed, thus judge whether wind energy conversion system gear-box breaks down, the type of fault and fault tightHeavy degree. The present invention can be applied to the fault diagnosis of gearbox of wind turbine, so that according to the actual motion of wind energy conversion system gear-boxSituation, obtains diagnostic result comparatively accurately fast, can under the operation of wind energy conversion system transmission system, carry out fault diagnosis to its gear-box,Can reflect timely and accurately gearbox fault situation, not only can, according to running state of gear box arrangement maintenance, avoid currentThe planned shutdown of wind energy turbine set maintenance, and can find as early as possible incipient fault risk, prevents the deterioration of degree of injury or fault coverageExpand, avoid contingency and casualties, reduce breakdown maintenance number of times and downtime, increase equipment normal operating time andDegree of reliability, improves utilization rate of equipment and installations, and the expense of saving a large amount of operation and maintenances, has improved wind-powered electricity generation enterprise operation economy; ThisInvention can overcome that at present the existing existing real-time of Fault Diagnosis of Gearbox for Wind Turbine technology is poor, fault False Rate is higherProblem, meets the actual needs of gearbox of wind turbine fault diagnosis preferably, have diagnostic result quick and precisely, Nonlinear QuasiClose advantage effective, that use is flexible, amount of calculation is little, can lay the foundation for realizing the inline diagnosis of wind energy conversion system gearbox fault.

Brief description of the drawings

Fig. 1 is the implementing procedure schematic diagram of embodiment of the present invention method.

Fig. 2 is the frame structure schematic diagram of the Fault Diagnosis of Gearbox for Wind Turbine device of application embodiment of the present invention method.

Detailed description of the invention

As shown in Figure 1, the implementation step of the Fault Diagnosis of Gearbox for Wind Turbine method of the present embodiment based on Kriging model is as follows:

1) gather respectively the multiple vibration signal sequences of wind energy conversion system gear-box under each operating mode pattern, each operating mode pattern comprises windNormal condition, state of wear and the broken conditions of power machine gear-box; In the present embodiment, specifically by being laid in wind energy conversion system gear-boxThe vibration signal sensor of appropriate location test point gathers the vibration signal of gear-box, and the fault-signal of gear-box adds by vibrationVelocity sensor gathers, and its sample frequency is 8192Hz, efficiently sampling time >=2s.

2) calculate respectively multiple fault-signal characteristic values of vibration signal sequence;

3) according to fault-signal characteristic value corresponding to multiple vibration signal sequences under each operating mode pattern, set up by membership functionCorresponding relation and sample data table between operating mode pattern, diagnosis desired value, each list item of sample data table comprises operating mode mouldThe data of formula, multiple fault-signal characteristic value, diagnosis desired value three types;

4), according to the some or all of data in sample data table, adopt least square method to carry out the plan of variogram theoretical modelClose, obtain parameter value and the variogram theoretical model of variogram theoretical model;

5) based on variogram theoretical model, according to Kriging method set up input vector that each fault-signal characteristic value forms,Approximate function between diagnosis desired value obtains the Kriging model being made up of the calculation expression of Kriging estimator,The input quantity of Kriging model is the input vector of each fault-signal characteristic value formation that vibration signal to be diagnosed is corresponding, KrigingThe output quantity of model is the Kriging estimator corresponding with diagnosis desired value;

6) detect the wind energy conversion system gear-box vibration signal to be diagnosed under in working order, calculate each fault of vibration signal to be diagnosedSignal characteristic value, the input vector input Kriging model that each fault-signal characteristic value of vibration signal to be diagnosed is formed,To the Kriging estimator of Kriging model output, between Kriging estimator inquiry operating mode pattern, diagnosis desired valueCorresponding relation, thereby determine operating mode pattern corresponding to vibration signal to be diagnosed.

To sum up step 1)～6), in the present embodiment, believe by gathering multiple vibrations of wind energy conversion system gear-box under each operating mode patternNumber sequence is also calculated multiple fault-signal characteristic values; Set up corresponding relation and sample data between operating mode pattern, diagnosis desired valueTable; Set up variogram theoretical model according to sample data table; Build Kriging model based on variogram theoretical model; DetectThe fault-signal characteristic value input Kriging model of vibration signal to be diagnosed obtains Kriging estimator, estimates according to KrigingCorresponding relation between amount inquiry operating mode pattern, diagnosis desired value is determined the operating mode pattern of vibration signal to be diagnosed, by operating mode mouldFormula, can determine the current state of wind energy conversion system gear-box (normal condition, or state of wear, or broken conditions), has solved currentExisting method for diagnosing faults and the problem that the existing real-time of technology is poor, fault False Rate is higher, for realizing wind energy conversion system gearCase on-line fault diagnosis is laid a good foundation, have diagnostic result quick and precisely, nonlinear fitting is effective, use is flexible, calculateMeasure little advantage.

In view of the uncertainty of wind turbine gearbox vibration signal, non-stationary and complexity, for Efficient Characterization fault-signal feature,Multiple fault-signal characteristic values in the present embodiment comprise Power Spectral Entropy, wavelet packet energy entropy, kurtosis, the degree of bias, correlation dimension and boxDimension. Wherein Power Spectral Entropy and wavelet-packet energy entropy can embody the complexity that signal energy distributes, and signal energy distributes simplerSingle, energy is more concentrated, and entropy is less; Otherwise signal is more complicated, energy is overstepping the bounds of propriety loose, and entropy is larger. Kurtosis coefficient andCoefficient of skewness index mainly reflects the size of impact energy, and therefore, in the time that gear breaks down, vibrational energy can have a greater change,Now kurtosis coefficient and the coefficient of skewness can characterize gear distress feature preferably. Box counting dimension and correlation dimension are simple, conventionalFractal dimension, effectively faults signal characteristic, box counting dimension and correlation dimension minimum under gear normal condition conventionally, andBox counting dimension and correlation dimension maximum under gear tooth breakage state.

In the present embodiment, step 2) in specifically refer to the Power Spectral Entropy of calculating vibration signal sequence according to formula (1);

$G_F=-\underset{i=1}{\overset{N}{\Σ}}{q}_i{\mathrm{logq}}_i---\left(1\right)$

In formula (1), G_FRepresent the Power Spectral Entropy of vibration signal sequence, q_iRepresent that the power spectrum of i vibration signal is in vibrationShared proportion in summed power spectrum in burst, N represents the quantity of vibration signal in vibration signal sequence. If vibration letterNumber sequence { x_iThe discrete Fourier transform of (1≤i≤N) is X (i), its power spectrum be S (i)=| X (i) |²/ (2 π N), i spectrum value q_iExpression formula suc as formula shown in (1-1).

$q_i=S\left(i\right)/\left(\underset{i=1}{\overset{N}{\Σ}}S\left(i\right)\right)---(1-1)$

In formula (1-1), S (i) represents the power spectrum of i vibration signal, and N represents the number of vibration signal in vibration signal sequenceAmount.

In the present embodiment, step 2) in specifically refer to the wavelet packet energy entropy that calculates vibration signal sequence according to formula (2);

$H_F=-\underset{i=1}{\overset{N}{\Σ}}{E}_i{\mathrm{logE}}_i---\left(2\right)$

In formula (2), H_FRepresent the wavelet packet energy entropy of vibration signal sequence, E_iThe frequency band energy that represents i vibration signal existsShared proportion in frequency band energy summation in vibration signal sequence, N represents the quantity of vibration signal in vibration signal sequence. IfA_ikFor the reconstruction coefficients of wavelet packet, the frequency band energy E of i vibration signal_iExpression formula suc as formula shown in (2-1).

$E_i=\underset{k=1}{\overset{n}{\Σ}}{A}_{ik}^2---(2-1)$

In formula (2-1), E_iInstitute in the frequency band energy summation of the frequency band energy that represents i vibration signal in vibration signal sequenceThe proportion accounting for, the number that n is discrete point, A_ikFor the reconstruction coefficients of wavelet packet.

In the present embodiment, step 2) in specifically refer to the kurtosis of calculating vibration signal sequence according to formula (3);

$K_q=\frac{1}{N}\underset{i=1}{\overset{N}{\Σ}}\frac{x_i^4}{{\left(\right({\mathrm{\Σ}}_{i=1}^N{x}_i^2)/N)}^2}-3---\left(3\right)$

In formula (3), K_qRepresent the kurtosis of vibration signal sequence, { x_i(1≤i≤N) expression vibration signal sequence, N represents vibrationThe quantity of vibration signal in burst.

In the present embodiment, step 2) in specifically refer to the degree of bias of calculating vibration signal sequence according to formula (4);

$K_p=\frac{1}{N}\frac{\left({\mathrm{\Σ}}_{i=1}^N\right|x_i^3\left|\right)/N}{{\left({\mathrm{\Σ}}_{i=1}^N{x}_i^2\right)}^{1.5}}---\left(4\right)$

In formula (4), K_pRepresent the degree of bias of vibration signal sequence, { x_i(1≤i≤N) expression vibration signal sequence, N represents vibrationThe quantity of vibration signal in burst.

In the present embodiment, step 2) in to calculate the detailed step of correlation dimension of vibration signal sequence as follows:

2.1) by vibration signal sequence { x_i(i=1,2 ..., n) carry out phase space reconfiguration and obtain nm vector y_iThe vector sequence of composition{y_i}(i＝1,2,…,nm)；

2.2) compute vectors sequence { y_i(i=1,2 ..., nm) in arbitrary reference point to the distance of all the other points;

2.3) respectively by vector sequence { y_i(i=1,2 ..., nm) in any two reference points as current vector pair, calculate currentThe distance that vector is right, compares distance right current vector and measurement vector, if should to related given parameters rDistance is less than given parameters r and judges that current vector is to being incidence vector pair, otherwise judges that current vector is to vowing as dereferencedIt is right to measure;

2.4) calculate the right correlation integral of each incidence vector according to formula (5);

$C_m\left(r\right)=\frac{1}{{\mathrm{nm}}^2}\underset{i,j=0}{\overset{nm}{\Σ}}H(r-r_{ij})---\left(5\right)$

In formula (5), C_m(r) represent the correlation integral that each incidence vector is right, nm represents vector sequence { y_i}(i＝1,2,…,nm)Middle vector y_iQuantity, r represents to weigh vector to related given parameters, r_ijRepresent vector sequence { y_i(i=1,2 ..., nm) inAs right i the reference point of current vector, to the distance between j reference point, H (x) is jump function, its expression formula asShown in formula (5-1);

$H\left(x\right)=\left\{\begin{array}{cc}1,& x>0,\\ 0,& x\≤0.\end{array}\right.---\left(5-1\right)$

2.5) according to the correlation dimension of formula (6) calculating vibration signal sequence;

$D_2\left(m\right)=\underset{\mathrm{\ϵ}\→0}{\lim }\frac{{\mathrm{lnC}}_m\left(r\right)}{\ln \left(r\right)}---\left(6\right)$

In formula (6), D₂(m) correlation dimension of expression vibration signal sequence, ε represents the little square box covering for fractalThe length of side, C_m(r) represent the correlation integral that each incidence vector is right, r represents to weigh vector to related given parameters.

In the present embodiment, step 2) in specifically refer to the box counting dimension that calculates vibration signal sequence according to formula (7);

$D_0=\underset{\mathrm{\ϵ}\→0}{\lim }\frac{\ln N\left(\mathrm{\ϵ}\right)}{\ln N(1/\mathrm{\ϵ})}---\left(7\right)$

In formula (7), D₀Represent the box counting dimension of vibration signal sequence, ε represents the limit of the little square box covering for fractalLong, 1/ ε represents the inverse of ε, and N () represents that utilizing the little square box of the corresponding length of side to carry out signal covers required box number.

It should be noted that, the correlation dimension and the box counting dimension that calculate vibration signal sequence are conventional vibration signal characteristics at present,Its computational methods are not limited to formula (6) and the formula (7) of above-mentioned record.

The present embodiment is according to fault-signal characteristic value corresponding to multiple vibration signal sequences under each operating mode pattern, by degree of membership letterNumber is set up corresponding relation and the sample data table between operating mode pattern, diagnosis desired value. The present embodiment utilizes existing fault modeWith the vibration signal characteristics value of corresponding failure, set up fault mode and the corresponding relation (fault mode of diagnosing desired value spanThe membership function of feature), wherein diagnosing desired value is possessed feature degree to be adopted fault mode in order to diagnose and to judgeSimple absolute figure carries out quantization signifying, and its sign function is the membership function of this fault mode feature, its spanBe the interval of this membership function. For example, fault mode adopts simple with the corresponding relation of diagnosis target spanPiecewise function represents, respectively corresponding its degree of membership interval of the diagnosis desired value of normal condition, state of wear, broken conditionsA numerical value, the degree of membership interval that might as well establish above-mentioned three kinds of fault modes is here respectively (0,1), (1,2), (2,3), diagnosisThe desired value of desired value is the median of degree of membership interval, 0.5,1.2 and 2.5, and as shown in table 1. In the present embodimentThe data of the sample data table arriving are specifically as shown in table 1.

Table 1: sample data table.

In the present embodiment, step 4) detailed step as follows:

4.1) select the some or all of list item in sample data table, by the multiple fault-signal characteristic value (merits in each list itemRate spectrum entropy, wavelet packet energy entropy, kurtosis, the degree of bias, correlation dimension, box counting dimension) jointly form an input vector as a sampleThis some X, the spatial coherence in the variogram survey region shown in use formula (8) between observation station, obtains Euclidean distance dCorresponding variogram value γ (d) when value changes;

$\mathrm{\γ}\left(d\right)=\frac{1}{2N\left(d\right)}\underset{r=1}{\overset{N\left(d\right)}{\Σ}}{\[y\left(X_i^r\right)-y\left(X_j^r\right)\]}^2---\left(8\right)$

In formula (8), γ (d) represents variogram value, and X represents sample point, and y (X) is the response observation of sample point X, and d is sampleThis between Euclidean distance,WithFor any two sample points in Euclidean distance d, N (d) is sample point pairNumber;

4.2) when Euclidean distance d value variation, corresponding variogram value γ (d), as input, adopts formula according to least square method(9) Gauss model shown in carries out matching, the variogram theory that obtains the parameter value of Gauss model and be made up of this Gauss modelModel;

$\mathrm{\γ}\left(d\right)=\left\{\begin{array}{cc}0& d=0\\ c_0+c\left(1-e^{-\frac{d^2}{a^2}}\right)& d>0\end{array}\right.---\left(9\right)$

In formula (9), d represents Euclidean distance, and γ (d) represents variogram value corresponding to Euclidean distance d, c₀, c, a be respectivelyThe parameter value of Gauss model, wherein c₀The piece gold that represents Gauss model, c represents the sagitta of Gauss model,Represent Gauss modelRange.

In the present embodiment, step 5) detailed step as follows:

5.1), for the sample data of given research object, utilize formula (10) and formula (11) to carry out the power of calculative determination observationHeavy coefficient w;

$\left\{\begin{array}{c}\mathrm{\Γ}={\left\{\mathrm{\γ}(x^{\left(i\right)}-x^{\left(j\right)})\right\}}_{ij},\\ \mathrm{\γ}=\{\mathrm{\γ}(x^{\left(1\right)}-x),...,\mathrm{\γ}(x^{\left(n\right)}-x)\}\·\end{array}\right.---\left(10\right)$

$w={(\mathrm{\γ}+I\frac{1-I^T{\mathrm{\Γ}}^{-1}\mathrm{\γ}}{I^T{\mathrm{\Γ}}^{-1}I})}^T{\mathrm{\Γ}}^{-1}---\left(11\right)$

In formula (9) and formula (10), γ represents the variogram value of utilizing formula (8) to calculate, and Γ represents by variogram meterThe matrix that calculation value γ forms as element, x⁽ⁱ⁾And x^(j)Any two sample points that represent sample data, w represents the power of observationHeavy coefficient, n is sample point number used in calculating, I represents n dimension unit column vector, the transposition computing of T representing matrix;

5.2) according to the weight coefficient w of observation, according to Kriging method set up input that each fault-signal characteristic value forms toAmount, the approximate function of diagnosis between desired value, obtain that the calculation expression of Kriging estimator shown in formula (12) formsKriging model, the input quantity of Kriging model be each fault-signal characteristic value that vibration signal to be diagnosed is corresponding form defeatedIncoming vector, the output quantity of Kriging model is the Kriging estimator corresponding with diagnosis desired value;

$\hat{z}\left(x\right)=\underset{i=1}{\overset{n}{\Σ}}{w}_{i}z\left(x^{\left(i\right)}\right),---\left(12\right)$

In formula (12),Represent the output quantity of Kriging model, w_iRepresent i the component of weight coefficient w, z (x⁽ⁱ⁾)Represent the input quantity of Kriging model, sample point number used during n represents to calculate. Known according to formula (11), given study areaIn the Ω of territory, n observation station and corresponding observation are { (x⁽ⁱ⁾,z(x⁽ⁱ⁾))|x⁽ⁱ⁾∈ Ω, i=1,2 ..., n}, waits arbitrarily to estimate an x (x ∈ Ω)Kriging estimatorCan represent by the linear combination of this n sample observations.

As shown in Figure 2, the device involving vibrations acceleration of the Fault Diagnosis of Gearbox for Wind Turbine method of application based on Kriging modelSensing unit 1, signal transmitting device 2, data acquisition unit 3 and analyzing and diagnosing unit 4. Vibration acceleration sensing unit 1 picks up in real timeGet the vibration acceleration signal on wind energy conversion system gear-box representative point and export to signal transmitting device 2, signal transmitting device 2 by inputVibration acceleration signal converts the voltage signal (or current signal of standard) of standard to and exports data acquisition unit 3, data toCollector 3 carries out the voltage/current signals of standard after detection and filtering, to export to analyzing and diagnosing unit 4, analyzing and diagnosing unit 4Adopt special application software to carry out analyzing and processing to gathered vibration signal, extract the variation characteristic of vibration signal, carry out sampleNotebook data is processed and variogram matching, and utilizes the variogram parameter value obtaining to set up Kriging model, examines by this modelThe malfunction of broken teeth roller box is also exported diagnostic result. In the present embodiment, it is blue moral that vibration acceleration sensing unit 1 adopts modelThe general vibration acceleration sensor of science and technology BZ1163; It is blue Deco skill that signal transmitting device 2 and data acquisition unit 3 adopt modelVibration acceleration signal conditioning and the collector of BZ2704; Analyzing and diagnosing unit 4 adopts computer realization, computer be loaded with byLabVIEW programs to realize operational modal analysis. The configuration of computer is as follows: CPU adopts the CPU of Intel, CPUDominant frequency is 2.93GHz, inside saves as 2GB, and hard disk is 500GB, and the computer in the present embodiment can be controlled 2 passage～256Port number is adopted the sampling of running simultaneously, and has met the measurement demand of multichannel, high accuracy, high speed Dynamic Signal.

The above is only the preferred embodiment of the present invention, and protection scope of the present invention is also not only confined to above-described embodiment, allThe technical scheme belonging under thinking of the present invention all belongs to protection scope of the present invention. It should be pointed out that common for the artTechnical staff, some improvements and modifications without departing from the principles of the present invention, these improvements and modifications also should be considered as thisThe protection domain of invention.

Claims (10)

1. the Fault Diagnosis of Gearbox for Wind Turbine method based on Kriging model, is characterized in that implementation step is as follows:

1) gather respectively the multiple vibration signal sequences of wind energy conversion system gear-box under each operating mode pattern, described each operating mode pattern bagDraw together normal condition, state of wear and the broken conditions of wind energy conversion system gear-box;

2) calculate respectively multiple fault-signal characteristic values of described vibration signal sequence;

3) according to fault-signal characteristic value corresponding to multiple vibration signal sequences under each operating mode pattern, set up by membership functionCorresponding relation and sample data table between operating mode pattern, diagnosis desired value, each list item of described sample data table comprises workThe data of condition pattern, multiple fault-signal characteristic value, diagnosis desired value three types;

4), according to the some or all of data in described sample data table, adopt least square method to carry out variogram theoretical modelMatching, obtain parameter value and the variogram theoretical model of variogram theoretical model;

5), based on described variogram theoretical model, set up according to Kriging method the input that each fault-signal characteristic value formsApproximate function between vector, diagnosis desired value obtains the Kriging mould being made up of the calculation expression of Kriging estimatorType, the input quantity of described Kriging model be the input that forms of each fault-signal characteristic value that vibration signal to be diagnosed is corresponding toAmount, the output quantity of described Kriging model is the Kriging estimator corresponding with diagnosis desired value;

6) detect the wind energy conversion system gear-box vibration signal to be diagnosed under in working order, each of vibration signal to be diagnosed described in calculatingFault-signal characteristic value, described in the input vector that each fault-signal characteristic value of vibration signal described to be diagnosed is formed is inputtedKriging model, obtains the Kriging estimator that Kriging model is exported, and inquires about described work according to described Kriging estimatorCorresponding relation between condition pattern, diagnosis desired value, thus determine operating mode pattern corresponding to vibration signal to be diagnosed.

2. the Fault Diagnosis of Gearbox for Wind Turbine method based on Kriging model according to claim 1, is characterized in that:Described multiple fault-signal characteristic value comprises Power Spectral Entropy, wavelet packet energy entropy, kurtosis, the degree of bias, correlation dimension and box counting dimension.

3. the Fault Diagnosis of Gearbox for Wind Turbine method based on Kriging model according to claim 2, is characterized in that:Described step 2) in specifically refer to the Power Spectral Entropy of calculating vibration signal sequence according to formula (1);

$G_F=-\underset{i=1}{\overset{N}{\Σ}}{q}_i{\mathrm{logq}}_i---\left(1\right)$

In formula (1), G_FRepresent the Power Spectral Entropy of vibration signal sequence, q_iRepresent that the power spectrum of i vibration signal is in vibrationShared proportion in summed power spectrum in burst, N represents the quantity of vibration signal in vibration signal sequence.

4. the Fault Diagnosis of Gearbox for Wind Turbine method based on Kriging model according to claim 2, is characterized in that:Described step 2) in specifically refer to the wavelet packet energy entropy that calculates vibration signal sequence according to formula (2);

$H_F=-\underset{i=1}{\overset{N}{\Σ}}{E}_i{\mathrm{logE}}_i---\left(2\right)$

In formula (2), H_FRepresent the wavelet packet energy entropy of vibration signal sequence, E_iThe frequency band energy that represents i vibration signal existsShared proportion in frequency band energy summation in vibration signal sequence, N represents the quantity of vibration signal in vibration signal sequence.

5. the Fault Diagnosis of Gearbox for Wind Turbine method based on Kriging model according to claim 2, is characterized in that:Described step 2) in specifically refer to the kurtosis of calculating vibration signal sequence according to formula (3);

$K_q=\frac{1}{N}\underset{i=1}{\overset{N}{\Σ}}\frac{x_i^4}{{\left(\right({\mathrm{\Σ}}_{i=1}^N{x}_i^2)/N)}^2}-3---\left(3\right)$

In formula (3), K_qRepresent the kurtosis of vibration signal sequence, { x_iRepresent vibration signal sequence, wherein 1≤i≤N, N tableShow the quantity of vibration signal in vibration signal sequence.

6. the Fault Diagnosis of Gearbox for Wind Turbine method based on Kriging model according to claim 2, is characterized in that:Described step 2) in specifically refer to the degree of bias of calculating vibration signal sequence according to formula (4);

$K_p=\frac{1}{N}\frac{\left({\mathrm{\Σ}}_{i=1}^N\right|x_i^3\left|\right)/N}{{\left({\mathrm{\Σ}}_{i=1}^N{x}_i^2\right)}^{1.5}}---\left(4\right)$

In formula (4), K_pRepresent the degree of bias of vibration signal sequence, { x_iRepresent vibration signal sequence, wherein 1≤i≤N, N tableShow the quantity of vibration signal in vibration signal sequence.

7. the Fault Diagnosis of Gearbox for Wind Turbine method based on Kriging model according to claim 2, is characterized in that:Described step 2) in to calculate the detailed step of correlation dimension of vibration signal sequence as follows:

2.1) by vibration signal sequence { x_iCarry out phase space reconfiguration and obtain nm vector y_iVector sequence { the y of composition_i, Qi ZhongzhenMoving burst { x_iMiddle i=1,2 ..., n, vector sequence { y_iMiddle i=1,2 ..., nm;

2.2) calculate described vector sequence { y_iIn arbitrary reference point to the distance of all the other points, wherein vector sequence { y_iIni＝1,2,…,nm；

2.3) respectively by described vector sequence { y_iIn any two reference points as current vector pair, wherein vector sequence { y_i}Middle i=1,2 ..., nm, calculates the right distance of current vector, and distance right current vector and measurement vector are given relatedDetermine parameter r and compare, if this distance is less than given parameters r, judges that current vector is to being incidence vector pair, otherwise sentenceVector before settled is to being dereferenced vector pair;

2.4) calculate the right correlation integral of each incidence vector according to formula (5);

$C_m\left(r\right)=\frac{1}{{\mathrm{nm}}^2}\underset{i,j=0}{\overset{nm}{\Σ}}H(r-r_{ij})---\left(5\right)$

In formula (5), C_m(r) represent the described right correlation integral of each incidence vector, nm represents vector sequence { y_iMiddle vector y_iQuantity, r represents to weigh vector to related given parameters, r_ijRepresent vector sequence { y_iIn as current vector rightI reference point is to the distance between j reference point, and H (x) is jump function; Wherein, vector sequence { y_iMiddle i=1,2 ..., nm;

2.5) according to the correlation dimension of formula (6) calculating vibration signal sequence;

$D_2\left(m\right)=\underset{\mathrm{\ϵ}\→0}{\lim }\frac{\ln C_m\left(r\right)}{ln\left(r\right)}---\left(6\right)$

In formula (6), D₂(m) correlation dimension of expression vibration signal sequence, ε represents the little square box covering for fractalThe length of side, C_m(r) represent the correlation integral that each incidence vector is right, r represents to weigh vector to related given parameters.

8. the Fault Diagnosis of Gearbox for Wind Turbine method based on Kriging model according to claim 2, is characterized in that:Described step 2) in specifically refer to the box counting dimension that calculates vibration signal sequence according to formula (7);

$D_0=\underset{\mathrm{\ϵ}\→0}{\lim }\frac{\ln N\left(\mathrm{\ϵ}\right)}{\ln N(1/\mathrm{\ϵ})}---\left(7\right)$

In formula (7), D₀Represent the box counting dimension of vibration signal sequence, ε represents the limit of the little square box covering for fractalLong, 1/ ε represents the inverse of ε, and N () represents that utilizing the little square box of the corresponding length of side to carry out signal covers required box number.

9. according to the Fault Diagnosis of Gearbox for Wind Turbine side based on Kriging model described in any one in claim 1～8Method, is characterized in that, described step 4) detailed step as follows:

4.1) select the some or all of list item in described sample data table, by the multiple fault-signal features in each list itemValue is common forms an input vector as a sample point X, observation station in the variogram survey region shown in use formula (8)Between spatial coherence, corresponding variogram value γ (d) when obtaining Euclidean distance d value and changing;

$\mathrm{\γ}\left(d\right)=\frac{1}{2N\left(d\right)}\underset{r=1}{\overset{N\left(d\right)}{\Σ}}{\[y\left(X_i^r\right)-y\left(X_j^r\right)\]}^2---\left(8\right)$

In formula (8), γ (d) represents variogram value, and X represents sample point, and y (X) is the response observation of sample point X, and d is sampleThis between Euclidean distance,WithFor any two sample points in Euclidean distance d, N (d) is sample point pairNumber;

4.2) when Euclidean distance d value variation, corresponding variogram value γ (d), as input, adopts formula according to least square method(9) Gauss model shown in carries out matching, the variogram theory that obtains the parameter value of Gauss model and be made up of this Gauss modelModel;

$\mathrm{\γ}\left(d\right)=\left\{\begin{array}{cc}0& d=0\\ c_0+c(1-e^{-\frac{d^2}{a^2}})& d>0\end{array}\right.---\left(9\right)$

In formula (9), d represents Euclidean distance, and γ (d) represents variogram value corresponding to Euclidean distance d, c₀, c, a be respectivelyThe parameter value of Gauss model, wherein c₀The piece gold that represents Gauss model, c represents the sagitta of Gauss model,Represent Gauss modelRange.

10. the Fault Diagnosis of Gearbox for Wind Turbine method based on Kriging model according to claim 9, is characterized in that,Described step 5) detailed step as follows:

5.1), for the sample data of given research object, utilize formula (10) and formula (11) to carry out the power of calculative determination observationHeavy coefficient w;

$\left\{\begin{array}{c}\mathrm{\Γ}={\left\{\mathrm{\γ}(x^{\left(i\right)}-x^{\left(j\right)})\right\}}_{ij},\\ \mathrm{\γ}=\{\mathrm{\γ}(x^{\left(1\right)}-x),...,\mathrm{\γ}(x^{\left(n\right)}-x)\}.\end{array}\right.---\left(10\right)$

$w={(\mathrm{\γ}+I\frac{1-I^T{\mathrm{\Γ}}^{-1}\mathrm{\γ}}{I^T{\mathrm{\Γ}}^{-1}I})}^T{\mathrm{\Γ}}^{-1}---\left(11\right)$

In formula (10) and formula (11), γ represents the variogram value of utilizing formula (8) to calculate, and Γ represents by variogramThe matrix that calculated value γ forms as element, x⁽ⁱ⁾And x^(j)Represent any two sample points of sample data, w represents observationWeight coefficient, n is sample point number used in calculating, I represents n dimension unit column vector, the transposition computing of T representing matrix;

5.2), according to the weight coefficient w of described observation, set up according to Kriging method defeated that each fault-signal characteristic value formsIncoming vector, the approximate function of diagnosis between desired value, obtain the calculation expression structure of Kriging estimator shown in formula (12)The Kriging model becoming, the input quantity of described Kriging model is each fault-signal characteristic value that vibration signal to be diagnosed is correspondingThe input vector forming, the output quantity of described Kriging model is the Kriging estimator corresponding with diagnosis desired value;

$\hat{z}\left(x\right)=\underset{i=1}{\overset{n}{\Σ}}{w}_{i}z\left(x^{\left(i\right)}\right),---\left(12\right)$

In formula (12),Represent the output quantity of Kriging model, w_iRepresent i the component of weight coefficient w, z (x⁽ⁱ⁾)Represent the input quantity of Kriging model, sample point number used during n represents to calculate.