CN104122086A - Kriging model based wind turbine gearbox fault diagnosis method - Google Patents

Abstract

The invention discloses a Kriging model based wind turbine gearbox fault diagnosis method which includes the steps of acquiring multiple vibration signal sequences of a wind turbine gearbox in different modes of working conditions and calculating multiple fault signal feature values; building a corresponding relation between the modes of working conditions and diagnosis target values and building a sample data list; building a variogram theoretical model according to the sample data list; building a Kriging model on the basis of the variogram theoretical model; detecting the fault signal feature values of to-be-diagnosed vibration signals and inputting the Kriging model to acquire a Kriging estimator, and inquiring the corresponding relation between the modes of working conditions and the diagnosis target values according to the estimator so as to determine the modes of working conditions of the to-be-diagnosed vibration signals. The Kriging model based wind turbine gearbox fault diagnosis method has the advantages of rapid diagnosis, accurate diagnosis result, good nonlinear fitting effect, use flexibility and small calculation amount, and can lay a foundation for online fault diagnosis of wind turbine gearboxes.

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 case, be specifically related to a kind of Fault Diagnosis of Gearbox for Wind Turbine method based on Kriging model.

Background technology

Wind-powered electricity generation unit is main mainly with transverse axis at present, comprise that Stall Type, change oar type, variable-pitch variable-speed type, direct-driving type, permanent-magnet type, half direct-driving type etc. are several, adopting the megawatt-level wind unit of the double-fed asynchronous variable-speed generator of gear case speedup is the mainstream model of current wind-power market, and wherein gear case is made up of gear case body, gear, transmission shaft, bearing and bearing seat.Along with the put into operation accumulation of time of the sharp increase of China's installed capacity of wind-driven power 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 case spoilage is up to 40～50%, and individual brand set gear box turnover rate almost reaches 100%.The direct and consequential damage bringing is thus also increasing, and the workload that maintainer drops into related work also constantly rises.Taking 1.5MW large scale wind power machine gearbox fault as example, the dismounting of fault gear case, transport and maintenance cost, up to more than 1,000,000 yuan, are 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.Therefore, gear case 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 means of wind energy conversion system gear case and the operation of wind-powered electricity generation unit and maintenance cost.The research of Fault Diagnosis of Gearbox for Wind Turbine method and application, can find incipient fault risk as early as possible, prevents the deterioration of degree of injury or the expansion of fault coverage, avoids mishap and casualties; Reduce breakdown maintenance number of times and stop time, increase equipment normal operating time and degree of reliability, improve plant factor; Reduce unnecessary part replacement, to reduce maintenance cost and stop time, save the expense of a large amount of operation and maintenances.Data shows, utilizes condition monitoring and fault diagnosis technology, can reduce accident rate 75%, reduces maintenance cost 25%～50%, makes 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, the most frequently used method of gearbox of wind turbine fault diagnosis is that the vibration signal to collecting carries out spectrum analysis, carry out pattern-recognition by artificial or machine, thereby 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 the fault diagnosis of wind power equipment system as fuzzy logic, expert system and neural network etc.But, because the relation between gearbox of wind turbine failure symptom and reason exists uncertain, and existing various intelligent diagnosing methods are perfect not enough at present, not only calculated amount is larger to cause current existing intelligent diagnosing method, real-time is poor, and fault False Rate is higher, be difficult to meet actual demands of engineering.

Kriging (Ke Lijin) method comes from a kind of interpolation method of South Africa sand smeller Krige Prospecting for Gold at first, and the one optimum that the French famous statistician Matheron of process theorizes and systematization grows up afterwards is gradually without bias estimation.The method is, in limited area, the value of regional change amount is carried out to nothing optimal estimation partially, by some sample points certificates of having measured in the limited neighborhood of sample point to be estimated (or region), consider the shape, size in sample spot (or region) and the space structure information with the geometric properties such as point to be estimated (or region) space distribution position each other and variable thereof, give certain weight coefficient to each sample, estimate thereby adopt method of weighted mean to treat the unknown quantity of estimating point (or region).Kriging method is applied to mining industry field at first, its range of application, well beyond original mining industry field, is all bringing into play in fields such as meteorologic analysis, environmental monitoring, image processing, mechanical engineering, structure optimization and Aero-Space the effect becoming more and more important in recent years.

Summary of the invention

The technical problem to be solved in the present invention is: for the above-mentioned technical matters of prior art, provide a kind of diagnostic result quick and precisely, nonlinear fitting is effective, use is flexible, calculated amount is little, can be for realizing the Fault Diagnosis of Gearbox for Wind Turbine method based on Kriging model that the inline diagnosis of wind energy conversion system gearbox fault lays the foundation.

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 case under each operating mode pattern, described each operating mode pattern comprises normal condition, state of wear and the broken conditions of wind energy conversion system gear case;

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

3) according to fault-signal eigenwert corresponding to multiple vibration signal sequences under each operating mode pattern, set up corresponding relation and the sample data table between operating mode pattern, diagnosis desired value by membership function, each list item of described sample data table comprises the data of operating mode pattern, multiple fault-signal eigenwert, 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 the matching of variogram theoretical model, obtain parameter value and the variogram theoretical model of variogram theoretical model;

5) based on described variogram theoretical model, the approximate function of setting up between input vector, the diagnosis desired value that each fault-signal eigenwert forms according to Kriging method obtains the Kriging model being made up of the calculation expression of Kriging estimator, the input quantity of described Kriging model is the input vector of each fault-signal eigenwert formation that vibration signal to be diagnosed is corresponding, and the output quantity of described Kriging model is and the corresponding Kriging estimator of diagnosis desired value;

6) detect the wind energy conversion system gear case vibration signal to be diagnosed under in working order, each fault-signal eigenwert of vibration signal to be diagnosed described in calculating, the input vector that each fault-signal eigenwert of vibration signal described to be diagnosed is formed is inputted described Kriging model, obtain the Kriging estimator of Kriging model output, inquire about the corresponding relation between described operating mode pattern, diagnosis desired value according to described Kriging estimator, thereby determine operating mode pattern corresponding to vibration signal to be diagnosed.

Preferably, described multiple fault-signal eigenwert comprises Power Spectral Entropy, wavelet packet energy entropy, kurtosis, the degree of bias, correlation dimension and 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}{\mathrm{\Σ}}}{q}_i\log q_i---\left(1\right)$

In formula (1), G _frepresent the Power Spectral Entropy of vibration signal sequence, q _ishared proportion in the summed power spectrum of the power spectrum that represents i vibration signal in vibration signal sequence, 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}{\mathrm{\Σ}}}{E}_i\log E_i---\left(2\right)$

In formula (2), H _frepresent the wavelet packet energy entropy of vibration signal sequence, E _ishared proportion in the frequency band energy summation of the frequency band energy that represents i vibration signal 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}{\mathrm{\Σ}}}\frac{x_i^4}{{(\left({\mathrm{\Σ}}_{i=1}^N{x}_i^2\right)/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 the quantity of vibration signal in vibration signal sequence.

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 the quantity of vibration signal in vibration signal sequence.

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 _ivector sequence { the y of composition _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 point are as current vector pair, calculate the current right distance of vector, by distance right current vector with weigh vector related given parameters r is compared, if this distance is less than given parameters r, judges that current vector is to being incidence vector pair, otherwise judge that current vector is to being dereferenced vector pair;

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

$C_m\left(r\right)=\frac{1}{{\mathrm{nm}}^2}\underset{i,j=0}{\overset{\mathrm{nm}}{\mathrm{\Σ}}}H(r-r_{\mathrm{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) and middle vector y _iquantity, r represents to weigh vector to related given parameters, r _ijrepresent vector sequence { y _i(i=1,2 ..., nm) in as right i the reference point of current vector to the distance between j reference point, H (x) is step 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{\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 length of side of the little square box covering for fractal, 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 length of side of the little square box covering for fractal, and 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, form an input vector as a sample point X using common the multiple fault-signal eigenwerts in each list item, spatial coherence in variogram survey region shown in use formula (8) between observation station, corresponding variogram value γ (d) while obtaining the variation of Euclidean distance d value;

$\mathrm{\γ}\left(d\right)=\frac{1}{2N\left(d\right)}\underset{r=1}{\overset{N\left(d\right)}{\mathrm{\Σ}}}{[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 observed reading of sample point X, d be sample point between Euclidean distance, with for any two sample points in Euclidean distance d, N (d) is sample point pair number;

4.2) when Euclidean distance d value variation, corresponding variogram value γ (d) is as input, adopt the Gauss model shown in formula (9) to carry out matching, the variogram theoretical model that obtains the parameter value of Gauss model and formed by this Gauss model according to least square method;

$\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 respectively the parameter value of Gauss model, wherein c ₀the piece gold that represents Gauss model, c represents the sagitta of Gauss model, represent the range of Gauss model.

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 weight coefficient w of calculative determination observed reading;

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

$\left\{\begin{array}{c}w={(\mathrm{\γ}+I\frac{1-I^T{\mathrm{\Γ}}^{-1}\mathrm{\γ}}{I^T{\mathrm{\Γ}}^{-1}I})}^T{\mathrm{\Γ}}^{-1},\\ \mathrm{\λ}=\frac{1-I^T{\mathrm{\Γ}}^{-1}\mathrm{\γ}}{I^T{\mathrm{\Γ}}^{-1}I}.\end{array}\right.---\left(\text{11}\right)$

In formula (10) and formula (11), γ represents to utilize the variogram value that shown in formula (8), variogram theoretical model calculates, and Γ represents the matrix being made up of as element variogram calculated value γ, x ⁽ⁱ⁾and x ^(j)any two sample points that represent sample data, w represents the weight coefficient of observed reading, and n is sample point number used in calculating, and I represents n dimension unit column vector, the transposition computing of T representing matrix.

5.2) according to the weight coefficient w of described observed reading, set up the approximate function between input vector, the diagnosis desired value that each fault-signal eigenwert forms according to Kriging method, obtain the Kriging model of the calculation expression formation of Kriging estimator shown in formula (12), the input quantity of described Kriging model is the input vector of each fault-signal eigenwert formation that vibration signal to be diagnosed is corresponding, and the output quantity of described Kriging model is and the corresponding Kriging estimator of diagnosis desired value;

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

In formula (12), the output quantity that represents Kriging model, wi represents 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 multiple vibration signal sequences under each operating mode pattern of wind energy conversion system gear case and calculating multiple fault-signal eigenwerts; 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; The fault-signal eigenwert input Kriging model that detects vibration signal to be diagnosed obtains 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, thereby judge whether wind energy conversion system gear case breaks down, the type of fault and the fault order of severity.The present invention can be applied to the fault diagnosis of gearbox of wind turbine, so that according to the practical operation situation of wind energy conversion system gear case, obtain fast diagnostic result comparatively accurately, can under the operation of wind energy conversion system kinematic train, carry out fault diagnosis to its gear case, can reflect timely and accurately gearbox fault situation, not only can overhaul according to running state of gear box arrangement, avoid the planned shutdown maintenance of current wind energy turbine set, and can find as early as possible incipient fault risk, prevent the deterioration of degree of injury or the expansion of fault coverage, avoid mishap and casualties, reduce breakdown maintenance number of times and stop time, increase equipment normal operating time and degree of reliability, improve plant factor, save the expense of a large amount of operation and maintenances, improve wind-powered electricity generation enterprise operation economy, the present invention can overcome the existing problem that the existing real-time of Fault Diagnosis of Gearbox for Wind Turbine technology is poor, fault False Rate is higher at present, meet preferably the actual needs of gearbox of wind turbine fault diagnosis, have advantages of diagnostic result quick and precisely, nonlinear fitting is effective, use is flexible, calculated amount 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 framed structure schematic diagram of the Fault Diagnosis of Gearbox for Wind Turbine device of application embodiment of the present invention method.

Embodiment

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 case under each operating mode pattern, each operating mode pattern comprises normal condition, state of wear and the broken conditions of wind energy conversion system gear case; In the present embodiment, specifically gather the vibration signal of gear case by being laid in the vibration signal sensor of wind energy conversion system gear case appropriate location check point, the fault-signal of gear case gathers by vibration acceleration sensor, and its sample frequency is 8192Hz, efficiently sampling time >=2s.

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

3) according to fault-signal eigenwert corresponding to multiple vibration signal sequences under each operating mode pattern, set up corresponding relation and the sample data table between operating mode pattern, diagnosis desired value by membership function, each list item of sample data table comprises the data of operating mode pattern, multiple fault-signal eigenwert, 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 matching of variogram theoretical model, obtain parameter value and the variogram theoretical model of variogram theoretical model;

5) based on variogram theoretical model, the approximate function of setting up between input vector, the diagnosis desired value that each fault-signal eigenwert forms according to Kriging method 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 eigenwert formation that vibration signal to be diagnosed is corresponding, and the output quantity of Kriging model is and the corresponding Kriging estimator of diagnosis desired value;

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

To sum up step 1)～6), in the present embodiment by gathering multiple vibration signal sequences under each operating mode pattern of wind energy conversion system gear case and calculating multiple fault-signal eigenwerts, 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, the fault-signal eigenwert input Kriging model that detects vibration signal to be diagnosed obtains Kriging estimator, according to Kriging estimator inquiry operating mode pattern, corresponding relation between diagnosis desired value is determined the operating mode pattern of vibration signal to be diagnosed, by operating mode pattern, can determine the state (normal condition that wind energy conversion system gear case is current, or state of wear, or broken conditions), current existing method for diagnosing faults and the existing real-time of technology are solved poor, the problem that fault False Rate is higher, lay a good foundation for realizing the inline diagnosis of wind energy conversion system gearbox fault, there is diagnostic result quick and precisely, nonlinear fitting is effective, use flexibly, the advantage that calculated amount is little.

In view of the uncertainty of wind turbine gearbox vibration signal, non-stationary and complicacy, for Efficient Characterization fault-signal feature, the multiple fault-signal eigenwerts in the present embodiment comprise Power Spectral Entropy, wavelet packet energy entropy, kurtosis, the degree of bias, correlation dimension and box counting dimension.Wherein Power Spectral Entropy and wavelet-packet energy entropy can embody the complicacy that signal energy distributes, i.e. signal energy distribution is simpler, and 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 and coefficient of skewness index mainly reflect the size of impact energy, and therefore, in the time that gear breaks down, vibrational energy can have a greater change, and now kurtosis coefficient and the coefficient of skewness can characterize gear distress feature preferably.Box counting dimension and correlation dimension are simple, conventional fractal dimensions, effectively faults signal characteristic, box counting dimension and correlation dimension minimum under gear normal condition conventionally, and under gear tooth breakage state box counting dimension and correlation dimension maximum.

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}{\mathrm{\Σ}}}{q}_i\log q_i---\left(1\right)$

In formula (1), G _frepresent the Power Spectral Entropy of vibration signal sequence, shared proportion in the summed power spectrum of the power spectrum that qi represents i vibration signal in vibration signal sequence, N represents the quantity of vibration signal in vibration signal sequence.If vibration signal sequence { x _ithe discrete Fourier transformation 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}{\mathrm{\Σ}}}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 quantity of vibration signal in vibration signal sequence.

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}{\mathrm{\Σ}}}{E}_i\log E_i---\left(2\right)$

In formula (2), H _frepresent the wavelet packet energy entropy of vibration signal sequence, E _ishared proportion in the frequency band energy summation of the frequency band energy that represents i vibration signal in vibration signal sequence, N represents the quantity of vibration signal in vibration signal sequence.If A _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}{\mathrm{\Σ}}}{A}_{\mathrm{ik}}^2---(2-1)$

In formula (2-1), E _ishared proportion in the frequency band energy summation of the frequency band energy that represents i vibration signal in vibration signal sequence, 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}{\mathrm{\Σ}}}\frac{x_i^4}{{(\left({\mathrm{\Σ}}_{i=1}^N{x}_i^2\right)/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 the quantity of vibration signal in vibration signal sequence.

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 the quantity of vibration signal in vibration signal sequence.

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 _ivector sequence { the y of composition _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 point are as current vector pair, calculate the current right distance of vector, by distance right current vector with weigh vector related given parameters r is compared, if this distance is less than given parameters r, judges that current vector is to being incidence vector pair, otherwise judge that current vector 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{\mathrm{nm}}{\mathrm{\Σ}}}H(r-r_{\mathrm{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) and middle vector y _iquantity, r represents to weigh vector to related given parameters, r _ijrepresent vector sequence { y _i(i=1,2 ..., nm) in as right i the reference point of current vector to the distance between j reference point, H (x) is step function, its expression formula is suc as formula shown in (5-1);

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

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 length of side of the little square box covering for fractal, 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 length of side of the little square box covering for fractal, and 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, and its computing method are not limited to formula (6) and the formula (7) of above-mentioned record.

The present embodiment, according to fault-signal eigenwert corresponding to multiple vibration signal sequences under each operating mode pattern, is set up corresponding relation and the sample data table between operating mode pattern, diagnosis desired value by membership function.The present embodiment utilizes the vibration signal characteristics value of existing fault mode and corresponding failure, set up fault mode and the corresponding relation (membership function of fault mode feature) of diagnosing desired value span, wherein diagnosing desired value is in order to diagnose and to judge fault mode, adopt simple absolute figure to carry out quantization signifying to possessed feature degree, its sign function is the membership function of this fault mode feature, and its span is the interval of this membership function.For example, fault mode adopts simple piecewise function to represent with the corresponding relation of diagnosis target span, a numerical value of respectively corresponding its degree of membership interval of the diagnosis desired value of normal condition, state of wear, broken conditions, here the degree of membership interval that might as well establish above-mentioned three kinds of fault modes is respectively (0,1), (1,2), (2,3), the expectation value of diagnosis desired value is the intermediate value of degree of membership interval, 0.5,1.2 and 2.5, as shown in table 1.The data of the sample data table obtaining in the present embodiment 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, multiple fault-signal eigenwerts in each list item (Power Spectral Entropy, wavelet packet energy entropy, kurtosis, the degree of bias, correlation dimension, box counting dimension) are formed to an input vector jointly as a sample point X, spatial coherence in variogram survey region shown in use formula (8) between observation station, corresponding variogram value γ (d) while obtaining the variation of Euclidean distance d value;

$\mathrm{\γ}\left(d\right)=\frac{1}{2N\left(d\right)}\underset{r=1}{\overset{N\left(d\right)}{\mathrm{\Σ}}}{[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 observed reading of sample point X, d be sample point between Euclidean distance, with for any two sample points in Euclidean distance d, N (d) is sample point pair number;

4.2) when Euclidean distance d value variation, corresponding variogram value γ (d) is as input, adopt the Gauss model shown in formula (9) to carry out matching, the variogram theoretical model that obtains the parameter value of Gauss model and formed by this Gauss model according to least square method;

$\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 respectively the parameter value of Gauss model, wherein c ₀the piece gold that represents Gauss model, c represents the sagitta of Gauss model, represent the range of Gauss model.

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 weight coefficient w of calculative determination observed reading;

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

$\left\{\begin{array}{c}w={(\mathrm{\γ}+I\frac{1-I^T{\mathrm{\Γ}}^{-1}\mathrm{\γ}}{I^T{\mathrm{\Γ}}^{-1}I})}^T{\mathrm{\Γ}}^{-1},\\ \mathrm{\λ}=\frac{1-I^T{\mathrm{\Γ}}^{-1}\mathrm{\γ}}{I^T{\mathrm{\Γ}}^{-1}I}.\end{array}\right.---\left(\text{11}\right)$

In formula (9) and formula (10), γ represents to utilize the variogram value that shown in formula (8), variogram theoretical model calculates, and Γ represents the matrix being made up of as element variogram calculated value γ, x ⁽ⁱ⁾and x ^(j)any two sample points that represent sample data, w represents the weight coefficient of observed reading, and n is sample point number used in calculating, and I represents n dimension unit column vector, the transposition computing of T representing matrix.

5.2) according to the weight coefficient w of observed reading, set up the approximate function between input vector, the diagnosis desired value that each fault-signal eigenwert forms according to Kriging method, obtain the Kriging model of the calculation expression formation of Kriging estimator shown in formula (12), the input quantity of Kriging model is the input vector of each fault-signal eigenwert formation that vibration signal to be diagnosed is corresponding, and the output quantity of Kriging model is and the corresponding Kriging estimator of diagnosis desired value;

$\hat{z}=\left(x\right)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{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), in given survey region Ω, n observation station and corresponding observed reading are { (x ⁽ⁱ⁾, z (x ⁽ⁱ⁾)) | x ⁽ⁱ⁾∈ Ω, i=1,2 ..., n}, waits arbitrarily to estimate the Kriging estimator of an x (x ∈ Ω) can represent by the linear combination of this n sample observations.

As shown in Figure 2, device involving vibrations acceleration sensing unit 1, signal transmitting device 2, data acquisition unit 3 and the analyzing and diagnosing unit 4 of the Fault Diagnosis of Gearbox for Wind Turbine method of application based on Kriging model.Vibration acceleration sensing unit 1 picks up in real time the vibration acceleration signal on wind energy conversion system gear case representative point and exports to signal transmitting device 2, signal transmitting device 2 converts the vibration acceleration signal of input the voltage signal (or current signal of standard) of standard to and exports data acquisition unit 3 to, data acquisition unit 3 carries out the voltage/current signals of standard to export to analyzing and diagnosing unit 4 after detection and filtering, analyzing and diagnosing unit 4 adopts special application software to carry out analyzing and processing to gathered vibration signal, extract the variation characteristic of vibration signal, carrying out sample data processes and variogram matching, and utilize the variogram parameter value obtaining to set up Kriging model, malfunction by this Model Diagnosis gear case is also exported diagnostic result.In the present embodiment, it is the general vibration acceleration sensor of blue Deco skill BZ1163 that vibration acceleration sensing unit 1 adopts model; It is vibration acceleration signal conditioning and the collector of blue Deco skill BZ2704 that signal transmitting device 2 and data acquisition unit 3 adopt model; Analyzing and diagnosing unit 4 adopts computer realization, and computing machine is loaded with by LabVIEW and programs to realize operational modal analysis.The configuration of computing machine is as follows: CPU adopts the CPU of Intel, CPU frequency is 2.93GHz, inside saves as 2GB, and hard disk is 500GB, computing machine in the present embodiment can control 2 passage～256 port numbers and adopt the sampling of running simultaneously, and met the measurement demand of hyperchannel, high precision, 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, and all technical schemes belonging under thinking of the present invention all belong to protection scope of the present invention.It should be pointed out that for those skilled in the art, some improvements and modifications without departing from the principles of the present invention, these improvements and modifications also should be considered as protection scope of the present 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 case under each operating mode pattern, described each operating mode pattern comprises normal condition, state of wear and the broken conditions of wind energy conversion system gear case;

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

3) according to fault-signal eigenwert corresponding to multiple vibration signal sequences under each operating mode pattern, set up corresponding relation and the sample data table between operating mode pattern, diagnosis desired value by membership function, each list item of described sample data table comprises the data of operating mode pattern, multiple fault-signal eigenwert, 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 the matching of variogram theoretical model, obtain parameter value and the variogram theoretical model of variogram theoretical model;

5) based on described variogram theoretical model, the approximate function of setting up between input vector, the diagnosis desired value that each fault-signal eigenwert forms according to Kriging method obtains the Kriging model being made up of the calculation expression of Kriging estimator, the input quantity of described Kriging model is the input vector of each fault-signal eigenwert formation that vibration signal to be diagnosed is corresponding, and the output quantity of described Kriging model is and the corresponding Kriging estimator of diagnosis desired value;

6) detect the wind energy conversion system gear case vibration signal to be diagnosed under in working order, each fault-signal eigenwert of vibration signal to be diagnosed described in calculating, the input vector that each fault-signal eigenwert of vibration signal described to be diagnosed is formed is inputted described Kriging model, obtain the Kriging estimator of Kriging model output, inquire about the corresponding relation between described operating mode pattern, diagnosis desired value according to described Kriging estimator, thereby 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 eigenwerts comprise 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}{\mathrm{\Σ}}}{q}_i\log q_i---\left(1\right)$

In formula (1), G _frepresent the Power Spectral Entropy of vibration signal sequence, q _ishared proportion in the summed power spectrum of the power spectrum that represents i vibration signal in vibration signal sequence, 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}{\mathrm{\Σ}}}{E}_i\log E_i---\left(2\right)$

In formula (2), H _frepresent the wavelet packet energy entropy of vibration signal sequence, E _ishared proportion in the frequency band energy summation of the frequency band energy that represents i vibration signal 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}{\mathrm{\Σ}}}\frac{x_i^4}{{(\left({\mathrm{\Σ}}_{i=1}^N{x}_i^2\right)/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 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 _i(1≤i≤N) expression vibration signal sequence, N represents 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 _i(i=1,2 ..., n) carry out phase space reconfiguration and obtain nm vector y _ivector sequence { the y of composition _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 point are as current vector pair, calculate the current right distance of vector, by distance right current vector with weigh vector related given parameters r is compared, if this distance is less than given parameters r, judges that current vector is to being incidence vector pair, otherwise judge that current vector is to being dereferenced vector pair;

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

$C_m\left(r\right)=\frac{1}{{\mathrm{nm}}^2}\underset{i,j=0}{\overset{\mathrm{nm}}{\mathrm{\Σ}}}H(r-r_{\mathrm{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) and middle vector y _iquantity, r represents to weigh vector to related given parameters, r _ijrepresent vector sequence { y _i(i=1,2 ..., nm) in as right i the reference point of current vector to the distance between j reference point, H (x) is step 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{\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 length of side of the little square box covering for fractal, 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 length of side of the little square box covering for fractal, and 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 method based on Kriging model described in any one in claim 1～8, it 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, form an input vector as a sample point X using common the multiple fault-signal eigenwerts in each list item, spatial coherence in variogram survey region shown in use formula (8) between observation station, corresponding variogram value γ (d) while obtaining the variation of Euclidean distance d value;

$\mathrm{\γ}\left(d\right)=\frac{1}{2N\left(d\right)}\underset{r=1}{\overset{N\left(d\right)}{\mathrm{\Σ}}}{[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 observed reading of sample point X, d be sample point between Euclidean distance, with for any two sample points in Euclidean distance d, N (d) is sample point pair number;

4.2) when Euclidean distance d value variation, corresponding variogram value γ (d) is as input, adopt the Gauss model shown in formula (9) to carry out matching, the variogram theoretical model that obtains the parameter value of Gauss model and formed by this Gauss model according to least square method;

$\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, and c0, c, a are respectively the parameter value of Gauss model, wherein c ₀the piece gold that represents Gauss model, c represents the sagitta of Gauss model, represent the range of Gauss model.

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 weight coefficient w of calculative determination observed reading;

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

$\left\{\begin{array}{c}w={(\mathrm{\γ}+I\frac{1-I^T{\mathrm{\Γ}}^{-1}\mathrm{\γ}}{I^T{\mathrm{\Γ}}^{-1}I})}^T{\mathrm{\Γ}}^{-1},\\ \mathrm{\λ}=\frac{1-I^T{\mathrm{\Γ}}^{-1}\mathrm{\γ}}{I^T{\mathrm{\Γ}}^{-1}I}.\end{array}\right.---\left(\text{11}\right)$

In formula (10) and formula (11), γ represents to utilize the variogram value that shown in formula (8), variogram theoretical model calculates, and Γ represents the matrix being made up of as element variogram calculated value γ, x ⁽ⁱ⁾and x ^(j)any two sample points that represent sample data, w represents the weight coefficient of observed reading, and n is sample point number used in calculating, and I represents n dimension unit column vector, the transposition computing of T representing matrix.

5.2) according to the weight coefficient w of described observed reading, set up the approximate function between input vector, the diagnosis desired value that each fault-signal eigenwert forms according to Kriging method, obtain the Kriging model of the calculation expression formation of Kriging estimator shown in formula (12), the input quantity of described Kriging model is the input vector of each fault-signal eigenwert formation that vibration signal to be diagnosed is corresponding, and the output quantity of described Kriging model is and the corresponding Kriging estimator of diagnosis desired value;

$\hat{z}=\left(x\right)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{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.