Fault feature extraction method based on EWT-SMF and MF-DFA for valve fault of reciprocating compressor

2.1.1. EWT

EWT is mainly divided into three steps:

Step 1. Perform a Fourier transform on the signal to get a spectrum diagram, and divide the spectrum by a specific algorithm to obtain a set of boundaries of the spectrum.

Reasonable spectral separation is the key to EWT. Gilles first proposed a spectral separation method based on the local maxima of the spectrum [5], but this method uses too few spectral information, and it is prone to modal aliasing when the wide modal and narrow modal neighbors or two large maxima on the spectrum belong to the same modality, etc. Then Gilles proposed a spectral space-based spectral separation method. The Gaussian kernel is used to scale the signal to obtain different scale spaces as shown in Eq. (1). The minimum value curve is obtained for the signal in the scale space, and then the position where the minimum value curve is larger than the threshold value is found, thereby determining the separation boundary. This method uses more information about the signal, so the resulting separation is better [19, 20]:

where:

The scale-space method finds local minimum curves in the scale space by searching different local minima described at multiple scales of the spectrum to achieve adaptive spectral separation:

From the spectrum diagram shown in Fig. 1(a), finding the boundaries in the spectrogram and dividing out the modality are similar to the fault diagnosis of rotating machinery that divides the spectrum by frequency doubling, and is to find the recesses in the spectrum, that is, the location of the local minimum. When we perform the scale-space transformation, the position and number of local minima will change. Let $x$ be the horizontal axis representing the position of the point; $s$ represents the number of scale transformations as the vertical axis, and the scaled space plane is constructed.

Fig. 1Example of modes in an spectrogram and its corresponding scale space representation

a) Simulation spectrogram

b) Spectrum support boundary

c) Scale-space curve

d) Screened scale-space curves based on threshold

After the transformation, the number of local minima described differently is a decreasing function with respect to the scale parameter $t$, and no new local minima will be generated as $t$ increases. For each scale transformation, the position of the new local minimum value obtained after the scale transformation is superimposed on the scale space plane. Therefore, the number of initial local minima can be written as $N_0$, and the position of each initial local minima will generate a scale-space curve $C_i\mathrm{}(i\in \left[1,N_0\right])$. The length of this curve depends on the number of occurrences of the $i$th local minimum in the entire scale space. Fig. 1(c) is an image of the scale space plane after scale transformation. Each initial local minimum value generates a curve in the scale space. However, not all local minima curves will remain. By using a classification algorithm and finding a suitable threshold, the minimum value curve is divided into two categories, so that a suitable scale-space curve is selected, and the position of the scale-space curve is finally defined as the separation boundaries of the spectrum as shown in Fig. 1(b) and Fig. 1(d).

Step 2. Based on the obtained boundary, use the empirical scale function and empirical wavelet function constructed by Meyer wavelet to obtain the orthogonal wavelet filter, and then EWT is performed on the spectrum to obtain a set of single-component or nearly-one-component AM-FM components.

The empirical scale function and empirical wavelet function of Meyer wavelet are as follows:

where $\beta \left(x\right)=x^4\left(35-84x+70x^2-20x^3\right)$, ${\tau }_n=\gamma {\omega }_n$, $\gamma <\mathrm{m}\mathrm{i}{\mathrm{n}}_n\left(\frac{{\omega }_{n+1}-{\omega }_n}{{\omega }_{n+1}+{\omega }_n}\right)$.

Then, as the wavelet analysis method, the inner product method is used to determine the detail coefficient and the approximate coefficient, respectively:

The empirical mode function can be expressed as:

Step 3. Analyze the AM-FM components using Hilbert transform to obtain the Hilbert spectrum.

After Gilles proposed EWT, Cao mentioned in [21] how empirical wavelet decomposition effectively separates the spectrum. In order to avoid erroneous segmentation caused by serious noise-induced internal modalities, Jinglong proposed that the signal can be denoised by wavelet-neighborhood coefficients to improve the SNR before EWT [22]. Mourad Kedadouche proposed an empirical wavelet transform method based on modal analysis (OMA) in the literature [23]. This paper uses the autoregressive moving average model to calculate the natural frequency and then spectrally separates them according to the natural frequency.

2.1.2. EWT modal optimization

Due to the complex frequency spectrum components of the reciprocating compressor vibration signal, including a wide frequency range and frequency susceptibility to mutation, the AM-FM components obtained through the EWT cannot be regarded as a single component modality. So, there is not enough theoretical support for the Hilbert transform.

The EWT separates the spectrum based on the mathematical properties of the signal spectrum and does not incorporate physical meaning. Therefore, the selection of physical meaning energy is the first standard of the optimum mode. The modality obtained by EWT is sorted by the energy from the largest to the smallest, and the first two modalities are selected for the next step. For reciprocating compressor, the impact signal contains the most abundant fault information. Therefore, it is expected that the mode in which the impact signal is less stable can be selected from the above two optimum modalities. According to the physical characteristics of reciprocating compressor, reciprocating compressor vibration signals can be determined the location of the main impact, coupled with the impact signal amplitude, it can be compared whether the two modes contain more fault information. Taking the reciprocating compressor valve as an example, the position of the generated impact signal is determined by the pressure test signal and the key phase signal when the suction valve is opened. The amplitude of the impact signal position should be at least 50 % of the original signal impact position amplitude.

Because of the impact of reciprocating compressor during operation, the spectrum components in both the normal operation state and the fault state are very complicated, and it is difficult to use the same standard to process normal and fault signals. Therefore, in this paper, a comparative method is used to process the vibration signal of reciprocating compressor in the normal state by EWT to obtain the optimal modal state. And then the spectrum separating boundary of this mode is used for fault diagnosis of other state signals.

2.1.3. State-adaptive morphological filtering

For meaningful modalities obtained by EWT, the resulting modalities need to be filtered in order to reduce the effects of the random noise and other vibrations on the vibration signal of the diagnosis site and highlight the fault features.

Morphological filtering is mainly divided into two parts in one-dimensional signal processing. The first part is the selection of structural elements. The structural elements in the morphological operation are similar to the filtering window in general signal processing, and only the signal elements that match the size and shape of the structural elements be effectively preserved. The three elements of a structural element are shape, length, and height. Commonly used structural elements are flat, triangular and semi-circular. It is generally believed that the flat structure is conducive to maintaining the shape characteristics of the signal being processed, the semi-circular structure is suitable for filtering out the interference of random noise, and the triangular structure is suitable for filtering the interference of the impulse noise. Most commonly used morphological filters use structural elements of a single structural characteristic. This method has a good effect on smooth signals [24], but the reciprocating compressor vibration signals are the nonstationary and nonlinear signals, and the impact and non-impact waveform of the signal are very different. If only a single structural element is used for filtering, both the impact and non-impact waveform cannot be accommodated at the same time, resulting in filtering out some useful information and thus affecting the effect of fault diagnosis. So, this paper uses morphological filtering to process the signal and proposes a new state-adaptive morphological filtering method. In accordance with the characteristics of the reciprocating compressor vibration signal, morphological filtering of signals is performed using different structural elements for different states.

The following Fig. 2 describes the valve vibration signal of reciprocating compressor as an example.

One cycle signal of reciprocating compressor can be divided into four stages of expansion, suction, compression and exhaust in Fig. 2. During the expansion process, the suction valve is closed and the signal amplitude is small. In the suction phase, the suction valve opens at the point A and hits the lift limiter. At this point, the first large peak of impact waveform is formed. During the suction process, the suction valve plate is theoretically at rest, and the impact waveform that appears in the middle of the suction stage caused by the opening of the other exhaust valve (the double acting cylinder). When the suction process is completed, the suction valve begins to close at point B and hits the valve seat, forming a second impact waveform that begins with the closing of the suction valve. The impact waveform near point C is caused by the opening of the exhaust valve on the side of the cylinder. In summary, the information of the valve operation state mainly includes the impact waveform caused by the valve during the opening and closing of the valve. When the valve plate is in a static state, theoretically no vibration signal generated by the valve plate. The vibration signal at this stage is mainly composed of the cylinder pressure pulsation, the fluctuation of the air flow in the cylinder caused by the opening and closing of the valve, the vibration caused by the inertial force, and the vibration of other parts [25], so the vibration signal at this stage cannot effectively reflect the operation state of the valve. From the above analysis, it can be concluded that the ideal filtering effect is to maintain the impact vibration waveform caused by the opening and closing of the gas valve, reduce the impact of other vibrations such as cylinder pressure pulsation, air flow fluctuations in the stable stage of the valve, and suppress impact noise and random noise. The most intuitive method is to distinguish the amplitude of the vibration signal. Because it takes a certain duration before the valve plate is opened and closed until it is in close contact with the lift limiter or valve seat, it is necessary to segment the vibration signal, and then determine the state of each segment based on the vibration amplitude; after confirming the status , select the type of structural element and height of structural element corresponding to each state in accordance with the established criteria.

Fig. 2Acceleration vibration signal of reciprocating compressor valve

The second part of morphological filtering is the morphological operations including four basic operators: morphological erosion, morphological expansion, morphology opening, and morphology closure. Let the original signal $f\left(n\right)$ be a discrete function defined on $F$$(\mathrm{1,2},\dots ,N-1)$, the structural element $g\left(m\right)$ is a discrete function defined on $G(\mathrm{1,2},\dots ,M-1)$, and $N\ge M$, then [26].

The $g\left(m\right)$ erosion of $f\left(n\right)$ is defined as:

The $g\left(m\right)$ expansion of $f\left(n\right)$ is defined as:

The $g\left(m\right)$ open operation of $f\left(n\right)$ is defined as:

The $g\left(m\right)$ closed operation of $f\left(n\right)$ is defined as:

where $\mathrm{\Theta }$, $\oplus$, ○, ● are operators for erosion, expansion, opening, and closing operations, respectively.

Four basic operations in morphology can be used to extract different signal feature information. But if only a single morphological operator is used, only information of a particular shape in the signal can be extracted. Therefore, a combination of basic operators is needed to more comprehensively extract the information from the signal. The commonly combinations of morphological operators filters. include the mean filter of corrosion and expansion in Eq. (12), the mean filter of open and closed operation in Eq. (13), the mean filter of open and closed-closed operation in Eq. (14), and the differential filter of erosion and expansion in Eq. (15):

Through repeated data comparison, this paper selects the mean filter of open and closed operation to filter and obtain the final filtering result.

2.2.1. Method description

The main calculation steps of MF-DFA for non-stationary time series $x_k$ is as follows:

1) Determine of signal outlines and segments:

Divide $Y\left(i\right)$ from positive and negative directions into non-overlapping $2N_s=int\left(N/s\right)$ segments of equal length $s$.

2) The local trend of each interval is calculated by least square fitting:

where $y_{\upsilon }$ is the trend of the polynomial data of the $\upsilon$ segment.

3) Define the $q$ order fluctuation function of the sequence:

4) The bilogarithmic function of $F_q\left(s\right)$ and $s$ is analyzed to determine the scaling property of fluctuation function, as $s$ increases, $F_q\left(s\right)$ increases in power law relationships:

where the range of the generalized Hurst exponent is $h\left(q\right)\in \left(\mathrm{0,1}\right)$ and the long-range correlation of the sequence can be judged by $h\left(q\right)$.

5) Plug $F^2(\upsilon ,s)=\left[Y\right(\upsilon s)-Y((\upsilon -1)s){]}^2$ into Eq. (20):

6) Suppose that $N_{\mathrm{s}}=N/s$:

7) The scaling exponent $\tau \left(q\right)$ is defined by the partition function $Z_q\left(s\right)$:

8) The relationship between $h\left(q\right)$ and the scaling exponent $\tau \left(q\right)$, the singular exponent $\alpha$ and multi-fractal spectrum $f\left(\alpha \right)$:as follows:

2.2.2. MF-DFA parameters

By observing the multi-fractal singular spectrum obtained by the MF-DFA method, it can be seen that the following important parameters can quantitatively express the probability distribution ratio and unevenness degree of fractal structure shown in Fig. 3: the spectral width $∆\alpha ={\alpha }_{max}-{\alpha }_{min}$describes the unevenness of probability measure distribution in the whole fractal structure. The peak value $f\left({\alpha }_o\right)$ of the multi-fractal spectrum describes the rate at which the number of identical probability units changes with the observation scale. The singular value ${\alpha }_o$ corresponds to the peak value $f\left({\alpha }_o\right)$ of multifractal spectrum. The difference $∆f=f\left({\alpha }_{max}\right)-f\left({\alpha }_{min}\right)$ in the fractal dimension of maximum and minimum probability subset, describes the difference between the number of elements in the maximum and minimum probability subsets [17].

Fig. 3Typical multi-fractal singular spectrum