CN107391935A - The instantaneous Frequency Estimation method examined based on non-delayed cost function and Grubbs - Google Patents

Abstract

The invention discloses an instantaneous frequency estimation method based on a non-delayed cost function and Grubbs test, using short-time Fourier transform to convert the original signal into a time spectrum graph, using Canny detection algorithm to obtain multiple ridge bands, and using Grubbs test to exclude each The outliers of the ridge bands are superimposed to construct a synthetic ridge band with complete and clear edges, and the Grubbs test is used to exclude the outliers of the synthetic ridge bands, and the average value curve of the synthetic ridge band is calculated, and the mean value curve is smoothed to calculate the smoothed mean value curve In the confidence interval at the 95% confidence level, the smooth mean curve and its confidence interval are mapped to the target ridge, and the reference line and local search interval of the target ridge are obtained, and the target ridge is extracted using a non-delayed cost function. The invention is suitable for estimating the instantaneous frequency of complex multi-component frequency conversion signals, overcomes the defects of the traditional method in estimating the instantaneous frequency of mechanical vibration signals, has high accuracy and precision of the estimation result, and is convenient for engineering application.

Description

Instantaneous frequency estimation method based on non-delayed cost function and Grubbs test

technical field

The invention relates to the field of state monitoring and fault diagnosis of rotating machinery, in particular to an instantaneous frequency estimation method based on a non-delayed cost function and a Grubbs test.

Background technique

Due to the complexity of the working environment, rotating machinery often works under variable speed conditions. Instantaneous frequency estimation is an important prerequisite for evaluating the operating state and fault diagnosis of rotating machinery. The most commonly used instantaneous frequency estimation method is the one-step cost function method. The one-step cost function method can search for ridge points in the local frequency range, but the center point of the local frequency range depends on the position of the previous ridge point, which leads to a delay in the one-step cost function. In addition, the width of the local frequency search range is set arbitrarily according to experience, and the width at any moment is fixed and cannot change with time, which leads to the lack of sufficient adaptability of the one-step cost function. The above defects lead to the low accuracy and precision of the one-step cost function method in estimating the instantaneous frequency.

Contents of the invention

The problem to be solved by the present invention is to propose an instantaneous frequency estimation method based on a non-delayed cost function and a Grubbs test for the above deficiencies. Compared with the existing method, the present invention uses the mapped smooth mean curve as the reference line of the target ridge, and uses the mapped confidence interval as the local search interval of the target ridge, so the center point of the local frequency search range does not depend on the upper The location of a ridge point without any delay, the local frequency search range can be automatically set, the search bandwidth can be automatically changed with time, and the accuracy and precision of the instantaneous frequency estimation results are high.

For solving above technical problem, the present invention adopts following technical scheme: provide the instantaneous frequency estimation method based on non-delay cost function and Grubbs test, it is characterized in that, comprises the following steps:

Step 1: Use the short-time Fourier transform algorithm to convert the signal x(k) (k=1, 2, ..., N) into a time-spectrogram, where N represents the length of the signal;

Step 2: Select a local area with a high signal-to-noise ratio from the time-spectrogram, and use the Canny detection algorithm to convert the local area into a binary image. The binary image contains multiple ridge bands; a local area refers to at least two Ridge bands, the area where the signal-to-noise ratio is greater than 80dB;

Step 3: Use the Grubbs test algorithm to exclude outliers at the upper and lower edges of each ridge band;

Step 4: Superimpose the above-mentioned multiple ridges on one of the ridges with the most complete outline according to the kinematic proportional relationship between them to construct a synthetic ridge with complete and clear edges; the kinematic proportional relationship refers to the the transmission ratio between the corresponding machine components;

Step 5: using the Grubbs test algorithm to exclude the outliers at the upper and lower edges of the synthetic ridge band;

Step 6: Calculate the average value curve of the above-mentioned synthetic ridge band, adopt the five-point cubic smoothing algorithm to smooth the average value curve, obtain the smooth average value curve, and calculate the confidence interval of the smooth average value curve on the 95% confidence level;

Step 7: Map the above smooth mean curve and its confidence interval to the target ridge line according to the kinematic proportional relationship between the smooth mean curve and the target ridge line to be estimated;

Step 8: Use the mapped smooth mean curve as the reference line of the target ridge, and use the mapped confidence interval as the local search interval of the target ridge;

Step 9: Use the non-delayed cost function to search for ridge points in the local search interval corresponding to each moment, determine the instantaneous frequency corresponding to each moment, and finally obtain the instantaneous frequency on the entire time interval.

Further, the short-time Fourier transform algorithm in the step 1 includes the following steps:

1) Perform a short-time Fourier transform on the signal x(k):

,

TF(t, f) represents the short-time Fourier transform result of the signal x(k), t represents the time factor, f represents the scale factor, and the function w(z) represents the window function whose independent variable is z;

2) Calculate the time spectrum of the signal x(k):

,

Spectrogram(t, f) represents the time spectrum of x(k).

Further, the Canny detection algorithm in the step 2 includes the following steps:

1) Use a Gaussian filter to smooth the image f(x, y) to eliminate noise and irrelevant details in the image:

,

,

g(x, y) represents the smoothed image, G(x, y) represents the two-dimensional Gaussian filter, x represents the time point of the image, y represents the frequency point of the image, the symbol * represents the convolution calculation, and σ represents the Gaussian standard Poor, σ=1 in the present invention;

2) Calculate the magnitude and direction of the g(x, y) intensity gradient

,

,

,

M(x, y) represents the magnitude of the intensity gradient, θ(x, y) represents the direction of the intensity gradient, g _x (x, y) represents the partial derivative of g(x, y) to x, g _y (x, y) represents the partial derivative of g(x, y) with respect to y;

3) Use non-maximum suppression to eliminate spurious edges:

In the gradient direction θ(x, y), if the non-zero gradient value is greater than the two adjacent gradient values, the non-zero gradient value remains unchanged, otherwise, the non-zero gradient value is set to zero;

4) Use two thresholds of different sizes to filter the above-mentioned image after eliminating false edges. The two thresholds are denoted as T ₁ and T ₂ , T ₁ < T ₂ , the image obtained by T ₁ is denoted as I ₁ , and the image obtained by T 1 is denoted as I 1 . ₂ The obtained image is recorded as I ₂ ; in the present invention, T ₁ =0.0063, T ₂ =0.0156;

5) Eliminate weak edges that are not connected with strong edges from I ₂ , and then connect the edges in I ₁ and I ₂ to form continuous edges.

Further, in the step 3, the Grubbs test algorithm comprises the following steps:

1) For the signal x _n (n=1, 2, ..., N), establish the Grubbs test statistic , Represents the sample mean, σ represents the sample standard deviation, and N represents the sample length;

2) Set the significance level to α, according to the probability formula Determine g ₀ (N, α), g ₀ (N, α) represents that the data length is N, and the significance level is the Grubbs critical value corresponding to α, and α=0.05 among the present invention;

3) if , (b=1, 2, …,N), then remove x _b .

Further, the non-delay cost function in said step 9 includes the following steps:

1) The local search interval FB _k corresponding to the kth moment is defined as

,

f _k (pmc) represents the value of the mapped smooth mean curve at the kth moment, Represents half of the confidence interval of the mapped smooth mean curve at the kth moment, and m represents the length of the target ridge;

2) The non-delayed cost function CF _k corresponding to the kth moment is defined as:

,

,

f _k (i) represents the frequency value taken within the range of FB _k , TF(t _k , f _k ) represents the value of TF(t, f) at the kth moment, t _k represents the value of t at the kth moment value, f _k represents the value of f at the kth moment, and e _k represents the weight factor.

Further, the relative error is ≤0.561%, and the average relative error is ≤0.046%.

The present invention adopts the above technical scheme, compared with the prior art, the present invention has the following advantages:

1) The present invention has real-time performance: the present invention takes the mapped synthetic ridge band smooth mean curve as a reference line, can instantly determine the center point of the local frequency search range at the current moment, avoids dependence on the previous ridge point, and eliminates time Delayed, real-time.

2) The present invention is self-adaptive: the present invention utilizes the confidence interval of the mapped synthetic ridge band smooth mean curve as the local frequency search interval, and can adaptively determine the local frequency search range corresponding to each moment, and the search bandwidth can follow the Time changes automatically, and there is no need to rely on experience to set the search bandwidth, thereby eliminating errors caused by human factors.

3) The experimental results show that: the maximum relative error between the estimated instantaneous frequency obtained by the present invention and the measured value is 0.561%, and the average relative error is 0.046%. Compared with the result of the one-step cost function method, the maximum relative error reduces by 96.58% , the average relative error is reduced by 97.85%.

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

Description of drawings

Accompanying drawing 1 is the flowchart of the instantaneous frequency estimation method based on non-delay cost function and Grubbs test in the embodiment of the present invention;

Accompanying drawing 2 is the planetary gearbox vibration signal in the embodiment of the present invention;

Accompanying drawing 3 is the time spectrum of planetary gearbox vibration signal in the embodiment of the present invention;

Accompanying drawing 4 is the local region with high signal-to-noise ratio selected from the time-spectrum diagram in the embodiment of the present invention;

Accompanying drawing 5 is the edge of the partial image area detected by Canny algorithm in the embodiment of the present invention;

Accompanying drawing 6 is the result after adopting Grubbs test algorithm to eliminate the abnormal point of each ridge band in the embodiment of the present invention;

Accompanying drawing 7 is the synthetic ridge band (the lowermost ridge band is the synthetic ridge band) formed by superimposing the kinematic proportional relationship between the ridge bands in the embodiment of the present invention;

Accompanying drawing 8 is the result after adopting Grubbs test algorithm in the embodiment of the present invention to eliminate the abnormal point of synthetic ridge band;

Accompanying drawing 9 is the mean smooth curve and 95% confidence interval thereof of synthetic ridge band in the embodiment of the present invention;

Accompanying drawing 10 is the mean smooth curve and its confidence interval that are mapped in the embodiment of the present invention;

Accompanying drawing 11 is the instantaneous frequency estimation value in the embodiment of the present invention.

detailed description

Embodiment, as shown in Figure 1, based on the instantaneous frequency estimation method of non-delay cost function and Grubbs test, comprises the following steps:

Step 1: Use the short-time Fourier transform algorithm to convert the signal x(k) (k=1, 2, ..., N) into a time-spectrogram, where N represents the length of the signal;

Step 2: Select a local area with a high signal-to-noise ratio from the time-spectrogram, and use the Canny detection algorithm to convert the local area into a binary image. The binary image contains multiple ridge bands; a local area refers to at least two Ridge bands, the area where the signal-to-noise ratio is greater than 80dB;

Step 3: Use the Grubbs test algorithm to exclude outliers at the upper and lower edges of each ridge band;

Step 4: Superimpose the above-mentioned multiple ridges on one of the ridges with the most complete outline according to the kinematic proportional relationship between them to construct a synthetic ridge with complete and clear edges; the kinematic proportional relationship refers to the the transmission ratio between the corresponding machine components;

Step 5: using the Grubbs test algorithm to exclude the outliers at the upper and lower edges of the synthetic ridge band;

Step 6: Calculate the average value curve of the above-mentioned synthetic ridge band, adopt the five-point cubic smoothing algorithm to smooth the average value curve, obtain the smooth average value curve, and calculate the confidence interval of the smooth average value curve on the 95% confidence level;

Step 7: Map the above smooth mean curve and its confidence interval to the target ridge line according to the kinematic proportional relationship between the smooth mean curve and the target ridge line to be estimated;

Step 8: Use the mapped smooth mean curve as the reference line of the target ridge, and use the mapped confidence interval as the local search interval of the target ridge;

Step 9: Use the non-delayed cost function to search for ridge points in the local search interval corresponding to each moment, determine the instantaneous frequency corresponding to each moment, and finally obtain the instantaneous frequency on the entire time interval.

The short-time Fourier transform algorithm in step 1 includes the following steps:

1) Perform a short-time Fourier transform on the signal x(k):

,

TF(t, f) represents the short-time Fourier transform result of the signal x(k), t represents the time factor, f represents the scale factor, and the function w(z) represents the window function whose independent variable is z;

2) Calculate the time spectrum of the signal x(k):

,

Spectrogram(t, f) represents the time spectrum of x(k).

Further, the Canny detection algorithm in the step 2 includes the following steps:

1) Use a Gaussian filter to smooth the image f(x, y) to eliminate noise and irrelevant details in the image:

,

,

g(x, y) represents the smoothed image, G(x, y) represents the two-dimensional Gaussian filter, x represents the time point of the image, y represents the frequency point of the image, the symbol * represents the convolution calculation, and σ represents the Gaussian standard Poor, σ=1 in the present invention;

2) Calculate the magnitude and direction of the g(x, y) intensity gradient

,

,

,

M(x, y) represents the magnitude of the intensity gradient, θ(x, y) represents the direction of the intensity gradient, g _x (x, y) represents the partial derivative of g(x, y) to x, g _y (x, y) represents the partial derivative of g(x, y) with respect to y;

3) Use non-maximum suppression to eliminate spurious edges:

In the gradient direction θ(x, y), if the non-zero gradient value is greater than the two adjacent gradient values, the non-zero gradient value remains unchanged, otherwise, the non-zero gradient value is set to zero;

4) Use two thresholds of different sizes to filter the above-mentioned image after eliminating false edges. The two thresholds are denoted as T ₁ and T ₂ , T ₁ < T ₂ , the image obtained by T ₁ is denoted as I ₁ , and the image obtained by T 1 is denoted as I 1 . ₂ The obtained image is recorded as I ₂ ; in the present invention, T ₁ =0.0063, T ₂ =0.0156;

5) Eliminate weak edges that are not connected with strong edges from I ₂ , and then connect the edges in I ₁ and I ₂ to form continuous edges.

The Grubbs test algorithm in step 3 includes the following steps:

1) For the signal x _n (n=1, 2, ..., N), establish the Grubbs test statistic , Represents the sample mean, σ represents the sample standard deviation, and N represents the sample length;

2) Set the significance level to α, according to the probability formula Determine g ₀ (N, α), g ₀ (N, α) represents that the data length is N, and the significance level is the Grubbs critical value corresponding to α, and α=0.05 among the present invention;

3) if , (b=1, 2, …,N), then remove x _b .

The non-delayed cost function in step 9 includes the following steps:

1) The local search interval FB _k corresponding to the kth moment is defined as

,

f _k (pmc) represents the value of the mapped smooth mean curve at the kth moment, Represents half of the confidence interval of the mapped smooth mean curve at the kth moment, and m represents the length of the target ridge;

2) The non-delayed cost function CF _k corresponding to the kth moment is defined as:

,

,

f _k (i) represents the frequency value taken within the range of FB _k , TF(t _k , f _k ) represents the value of TF(t, f) at the kth moment, t _k represents the value of t at the kth moment value, f _k represents the value of f at the kth moment, and e _k represents the weight factor.

The performance of the algorithm described in the present invention is verified by using the vibration data of the fan turbine planetary gearbox.

The vibration data is collected from the gearbox housing close to the planetary gear train, the data length N=2736825, and the sampling frequency fs=5000 Hz.

The collected vibration data of the planetary gearbox are shown in Fig. 2.

The short-time Fourier transform algorithm is used to convert the vibration data of the planetary gearbox shown in Figure 2 into a time-spectrum diagram, and the resulting time-spectrum diagram is shown in Figure 3.

Select a local area with a high signal-to-noise ratio from the time-spectrum diagram shown in Figure 3, and the obtained local area is shown in Figure 4.

The Canny detection algorithm is used to detect the edge of the local area shown in Figure 4, and the obtained image edge is shown in Figure 5.

The Grubbs test algorithm is used to eliminate the abnormal points of each ridge band in Figure 5, and the results are shown in Figure 6.

According to the kinematic proportional relationship between the ridges, each ridge is superimposed on one of the ridges with the most complete outline, and the synthetic ridge constructed is shown in Figure 7 (the lowermost ridge is the synthetic ridge) .

The Grubbs test algorithm is used to eliminate the abnormal points of the synthetic ridge belt, and the result is shown in Fig. 8.

The smoothed mean curves of the synthetic ridge bands and their 95% confidence intervals were calculated, and the results are shown in Fig. 9.

According to the kinematic proportional relationship between the smooth mean curve and the target ridge, the smooth mean curve and its confidence interval are mapped to the target ridge, and the results are shown in Figure 10.

Using the non-delayed cost function to search for the ridge point of the target ridge, the obtained instantaneous frequency curve is shown in Figure 11.

Show through many experiments, the maximum relative error between the estimated value of instantaneous frequency obtained by the present invention and the measured value is 0.561%, and the average relative error is 0.046%, and the difference between the estimated value of instantaneous frequency obtained by the one-step cost function method and the measured value is 0.561%. The maximum relative error between them is 16.39%, and the average relative error is 2.14%. The maximum relative error of the present invention is reduced by 96.58%, and the average relative error is reduced by 97.85%.

According to the experimental results, it is believed after analysis that:

1) The traditional one-step cost function needs to rely on the position of the previous ridge point when determining the center point of the current search interval, and there is a time delay phenomenon. The present invention uses the mapped smooth mean curve as a reference line to instantly determine the current search interval The central position of , does not depend on the last ridge point at all, so it has real-time performance.

2) The traditional one-step cost function method lacks adaptability, and needs to manually set the search interval, and the search width is fixed, which inevitably brings errors. The present invention uses the mapped smooth mean curve confidence interval to automatically determine the local search Interval, the search bandwidth can change automatically as time changes without manual participation, so it is adaptive.

3) Compared with the traditional one-step cost function method, the present invention has high precision and accuracy.

Those skilled in the art should realize that the above-mentioned specific embodiments are only exemplary, and are intended to enable those skilled in the art to better understand the content of the present invention, and should not be construed as limiting the protection scope of the present invention. The improvements made in the technical solution of the invention all fall into the protection scope of the present invention.

Claims (6)

1. the instantaneous Frequency Estimation method examined based on non-delayed cost function and Grubbs, it is characterised in that including following step Suddenly：

Step 1：Using Short Time Fourier Transform algorithm by signal x（k）（k=1, 2, …,N）Time-frequency spectrum is converted to, N is represented The length of signal；

Step 2：One piece of regional area having compared with high s/n ratio is chosen from time-frequency spectrum, should using Canny detection algorithms Regional area is converted into bianry image, and bianry image includes more vallate bands；

Step 3：Excluded using Grubbs check algorithms per exceptional value of the vallate with lower edges；

Step 4：Above-mentioned more vallate bands are most complete according to the mutual kinematics proportionate relationship wherein profile that is added to Ridge band on, structure one have complete sharp edge synthesis ridge band；

Step 5：Above-mentioned exceptional value of the synthesis ridge with lower edges is excluded using Grubbs check algorithms；

Step 6：Calculate the Mean curve of above-mentioned synthesis ridge band, using 5 points three times smoothing algorithm Mean curve is smoothly located Reason, obtains smooth Mean curve, calculates the confidential interval of the smooth Mean curve in 95% confidence level；

Step 7：By above-mentioned smooth Mean curve and its confidential interval according between smooth Mean curve and target crestal line to be estimated Kinematics proportionate relationship be mapped on target crestal line；

Step 8：Reference line using the smooth Mean curve after mapping as target crestal line, using the confidential interval after mapping as mesh Mark the Local Search section of crestal line；

Step 9：Ridge point is searched in the Local Search section corresponding to each moment using non-delayed cost function, it is determined that each Instantaneous frequency corresponding to moment, finally obtain the instantaneous frequency on whole time interval.

2. the instantaneous Frequency Estimation method according to claim 1 examined based on non-delayed cost function and Grubbs, its It is characterised by, Short Time Fourier Transform algorithm comprises the following steps in the step 1：

1）To signal x（k）Carry out Short Time Fourier Transform：

,

TF (t, f) representation signal x（k）Short Time Fourier Transform result, t represents time factor, and f represents scale factor, function W (z) represents window function of the independent variable as z；

2）Calculate signal x（k）Time-frequency spectrum：

,

Spectrogram (t, f) represents x（k）Time-frequency spectrum.

3. the instantaneous Frequency Estimation method according to claim 1 examined based on non-delayed cost function and Grubbs, its It is characterised by, Canny detection algorithms comprise the following steps in the step 2：

1) image f (x, y) is smoothed using Gaussian filter, eliminates noise and Extraneous details in image：

,

,

Image after g (x, y) representatives are smooth, G (x, y) represent 2-d gaussian filterses device, the time point of x representative images, and y is represented The Frequency point of image, symbol * represent convolutional calculation, and it is poor that σ represents Gauss standard；

2) amplitude and direction of g (x, y) intensity gradient are calculated

,

,

,

M (x, y) represents the amplitude of intensity gradient, and θ (x, y) represents the direction of intensity gradient, g_x(x, y) represents g (x, y) to x Partial derivative, g_y(x, y) represents partial derivatives of the g (x, y) to y；

3) false edge is eliminated using non-maximum coercion Acts：

On gradient direction θ (x, y), if non-zero gradient value is more than two adjacent Grad, the non-zero gradient value is kept It is constant, otherwise, the non-zero gradient value zero setting；

4) image after above-mentioned elimination false edge is filtered using two threshold values of different sizes, two threshold values are designated as T₁ And T₂, T₁<T₂, by T₁Obtained image is designated as I₁, by T₂Obtained image is designated as I₂；

5) from I₂It is middle to reject the weak edge not being connected with strong edge, then connect I₁And I₂In edge formed continuous boundary.

4. the instantaneous Frequency Estimation method according to claim 1 examined based on non-delayed cost function and Grubbs, its It is characterised by：Grubbs check algorithms comprise the following steps in the step 3：

1）To signal x_n（n=1, 2, …,N）, establish Grubbs test statistics,Representative sample average, σ Representative sample standard deviation, N representative sample length；

2）Significance is set as α, according to new probability formulaDetermine g₀(N, α), g₀(N, α) generation Table data length is N, and significance is the Grubbs critical values corresponding to α；

If 3), （b=1, 2, …,N）, then x is rejected_b。

5. the instantaneous Frequency Estimation method according to claim 1 examined based on non-delayed cost function and Grubbs, its It is characterised by：Non-delayed cost function comprises the following steps in the step 9：

1）Local Search section FB corresponding to k-th of moment_kIt is defined as

,

f_k(pmc) value of the smooth Mean curve after mapping k-th of moment is represented,It is bent to represent the smooth average after mapping For line confidential interval in the half of k-th of moment width, m represents the length of target crestal line；

2）Non-delayed cost function CF corresponding to k-th of moment_kIt is defined as：

,

,

f_k(i) represent in FB_kIn the range of the frequency values that are taken, TF (t_k, f_k) represent values of the TF (t, f) k-th of moment, t_k Represent values of the t k-th of moment, f_kRepresent values of the f k-th of moment, e_kRepresent weight factor.

6. the instantaneous Frequency Estimation method according to claim 1 examined based on non-delayed cost function and Grubbs, its It is characterised by：Relative error≤0.561%, average relative error≤0.046%.