CN110458073A - A Feature Extraction Method of Optical Fiber Vibration Signal Based on MEEMD-Hilbert and Multilayer Wavelet Decomposition - Google Patents

Abstract

The present invention relates to the feature extraction method of optical fiber vibration signal based on MEEMD-Hilbert and multilayer wavelet decomposition, is a kind of method that distributed optical fiber vibration signal is carried out feature extraction, belongs to the field of signal processing and machine learning, is characterized in that adopting following steps: (1) Determine the vibration signal after introducing the white noise signal; (2) Determine the first IMF component sequence set; (3) Perform time-delay space reconstruction; (4) Determine the permutation entropy; (5) Determine the remaining components; ( 6) Calculate the Hilbert transform of the sequence; (7) Determine the analytical signal of the sequence and perform autocorrelation processing; (8) Discrete wavelet transform; (9) Calculate the average energy on different frequency bands; (10) Calculate each frequency band The average energy ratio above. The invention has higher time-frequency resolution, and provides a method with obvious effect for feature extraction of optical fiber vibration signals.

Description

A fiber vibration signal characteristic based on MEEMD-Hilbert and multilayer wavelet decomposition extraction method

technical field

The invention relates to the fields of signal processing and machine learning, and mainly relates to a method for feature extraction of distributed optical fiber vibration signals.

Background technique

At present, for the feature extraction of distributed optical fiber vibration signals, it is mainly based on time domain and frequency domain feature extraction or wavelet analysis. Signal time-domain features are divided into short-term features and long-term sequence features. The short-term feature is to extract the envelope information of the vibration signal, which mainly includes the rising edge time, falling edge time and amplitude peak value. The long-term series feature is the time series change feature of the signal over a long period of time, which mainly includes the number of zero crossings, the number of vibration segments, and the average amplitude peak value, etc. However, only extracting the time-domain and frequency-domain features of the signal can only represent part of the information of the vibration signal, which makes it impossible to improve the efficiency and effect of recognition when the signal is subsequently recognized. Wavelet decomposition is to decompose the signal into multiple layers and use the energy distribution of the signal in each sub-band as the feature of the signal, but the feature of the wavelet component will be mixed with low-frequency noise features, and with the increase of the number of decomposition layers, the feature dimension will become smaller. The higher the value, the higher the computational burden.

Vibration signal feature extraction is mainly used in the field of optical fiber early warning. There are already some relatively mature feature extraction methods, such as parameter analysis, gene cycle, Fourier analysis, wavelet decomposition, etc. However, there are many types of intrusion signals and the degree of danger of each type is different. Faced with a variety of complex intrusion types, accurate extraction of stable and representative target vibration source signal characteristics is one of the core links of the optical fiber early warning system. Sensitive to distinguishing effective features of optical fiber vibration signals in noisy environments, high recognition rate, convenient and accurate feature extraction algorithm, providing strong support for the field of optical fiber early warning.

Contents of the invention

In view of the problems existing in the above-mentioned prior art, the technical problem to be solved by the present invention is to provide a method for feature extraction of optical fiber vibration signals based on MEEMD-Hilbert and multi-layer wavelet decomposition, the specific process of which is shown in FIG. 1 .

The implementation steps of the technical solution are as follows:

(1) Determine the vibration signal after introducing the white noise signal and

Introducing white noise signals n _p (t) and -n _p (t) into the vibration signal x(t), we get

In the formula, and Indicates the vibration signal after introducing the white noise signal, x(t) indicates the original optical fiber vibration signal, n _p (t) and -n _p (t) indicate the white noise signal, a _p indicates the amplitude of the noise signal introduced for the pth time, p=1, 2, . . . , N _noise , where N _noise represents the total number of noises introduced.

(2) Determine the first IMF component sequence set I ₁ (t):

right and Perform EMD decomposition separately to obtain the first set of component sequence sets and The components with the same subscript number in the sequence set are summed, accumulated, and averaged to obtain I ₁ (t):

In the formula, I ₁ (t) represents the first IMF component sequence set, and Represents the first group of IMF component sequence sets, p=1, 2, . . . , N _noise , where N _noise represents the total number of noises introduced.

(3) Delay space reconstruction:

Delay space reconstruction is performed on the digital sequence I ₁ (q) corresponding to I ₁ (t), and the following sequence is obtained:

In the formula, In ( _n ) represents the sequence after delay space reconstruction, q=1, 2, ..., n, ..., N, N represents the total length of the sequence, Δn represents the sequence delay length, m represents Spatial reconstruction dimensionality.

(4) Determine the permutation entropy S(q) of I ₁ (q):

In the formula, g=1, 2, ..., k represents the type of sequence number, m is the spatial reconstruction dimension,

(5) Determine the remaining component r(t):

Set the threshold of S(q), when S(q) is lower than the threshold, judge I ₁ (t) as a non-abnormal signal, and remove it as the first IMF component from the original signal x(t), that is r(t)=x(t)-I ₁ (t), the residual component r(t) is obtained.

Repeat steps (1)-(5) for r(t), and then get the IMF components I(t)=I ₁ (t), I ₂ (t), ..., I _s (t), where s is the IMF The number of components.

(6) Calculate the Hilbert transform of I(t):

In the formula, Represents the Hilbert transform of I(t), I(t) represents the decomposed IMF component, t represents time, and τ represents the time interval.

(7) Determine the analytical signal g(t) of I(t) and perform autocorrelation processing:

In the formula, I(t) represents the decomposed IMF component, Represents the Hilbert transform of I(t), j represents the imaginary number unit, g(t) represents the analytical signal, y(t) represents the signal after autocorrelation processing, z(t) represents the autocorrelation signal template, and t represents time , τ represents the time interval.

(8) Discrete wavelet transform:

In the formula, α represents the scale factor, ψ(t) represents the mother wavelet function, τ represents the displacement information provided by the mother wavelet function, t represents time, y(t) represents the signal after autocorrelation processing, ξ _i (n) represents the discrete wavelet Transformed wavelet coefficients, i represents the current i-th frequency band, n represents the current n-th frequency band. Therefore, under different scales and displacements, multi-level wavelet decomposition is realized, and the schematic diagram of the process is shown in Figure 3.

(9) Calculate the average energy on different frequency bands:

In the formula, ξ _i (n) represents the wavelet coefficient on each frequency band, i represents the current i-th frequency band, n represents the current n-th frequency band, and N represents the length of wavelet coefficients in each frequency band.

(10) Calculate the average energy ratio Ω _i in each frequency band:

In the formula, Ω _i represents the average energy ratio of each frequency band, E _sum represents the sum of energy of all frequency bands, and E _i represents the energy of the current frequency band.

The average energy ratio obtained from the highest frequency coefficient is not significantly different from other coefficients. In order to reduce the number of redundant features and eliminate the correlation between features, the energy ratio of the highest frequency coefficient is removed. Therefore, the finally obtained eigenvector is e={Ω ₁ , Ω ₂ , Ω ₃ , Ω ₄ , Ω ₅ } ^T .

The present invention has the advantage over prior art:

(1) The present invention applies the multi-layer wavelet decomposition method to the feature extraction of optical fiber vibration signals, extracts the essential features of the vibration signals in different frequency bands, fully considers the detailed features of the vibration signals in each frequency band, and can reflect the time The feature information such as domain mutation position and corresponding frequency has high time-frequency resolution.

(2) the present invention combines wavelet decomposition technology and MEEMD-Hilbert technology, it is applied in the feature extraction of optical fiber vibration signal, has obtained obvious feature extraction effect compared with prior art, illustrates and utilizes the present invention to carry out optical fiber vibration signal Feature extraction can achieve better results.

Description of drawings

In order to better understand the present invention, further description will be made below in conjunction with the accompanying drawings.

Fig. 1 is the flow chart of the steps of establishing the fiber vibration signal feature extraction method based on MEEMD-Hilbert and multi-layer wavelet decomposition;

Fig. 2 is to establish the flow chart of the feature extraction method of optical fiber vibration signal based on MEEMD-Hilbert and multi-layer wavelet decomposition;

Fig. 3 is a schematic diagram of multilayer wavelet decomposition process;

Fig. 4 is a schematic diagram of energy ratio of wavelet coefficients of four kinds of optical fiber vibration signals;

specific implementation plan

The present invention will be described in further detail below through examples of implementation.

In this implementation case, four typical vibration signals of electric drill, pedestrian trampling, excavator excavation, and vehicle passing were selected for experiments. Among them, the electric drill signal is generated by mechanical operations, and its inherent rotation frequency can generate continuous and uninterrupted mechanical signals, and the amplitude and frequency are relatively regular and stable. The number of acquisitions of each type of vibration signal is 30 times, and the sampling frequency is 2KHz, corresponding to four vibration signals, and there are 120 sets of experimental data in total.

The overall flow of the optical fiber vibration signal feature extraction algorithm provided by the present invention is shown in Figure 1, and the specific steps are as follows:

(1) Determine the vibration signal after introducing the white noise signal and

Introducing white noise signals n _p (t) and -n _p (t) into the vibration signal x(t), we get

In the formula, and Indicates the vibration signal after introducing the white noise signal, x(t) indicates the original optical fiber vibration signal, n _p (t) and -n _p (t) indicate the white noise signal, a _p indicates the amplitude of the noise signal introduced for the pth time, p=1, 2, . . . , N _noise , where N _noise represents the total number of noises introduced. In this example, the amplitudes of the noise signals introduced are 1.22, 1.37, 0.15, 0.87, 4.32, 1.27, 3.98, 2.61, 1.95, 0.21, and the total number of noises N _noise is set to 10.

(2) Determine the first IMF component sequence set I ₁ (t):

right and Perform EMD decomposition separately to obtain the first set of component sequence sets and The components with the same subscript number in the sequence set are summed, accumulated, and averaged to obtain I ₁ (t):

In the formula, I ₁ (t) represents the first IMF component sequence set, and Indicates the first group of IMF component sequence sets, and p indicates the noise introduced for the pth time.

(3) Delay space reconstruction:

Delay space reconstruction is performed on the digital sequence I ₁ (q) corresponding to I ₁ (t), and the following sequence is obtained:

In the formula, In ( _n ) represents the sequence after delay space reconstruction, q=1, 2, ..., n, ..., N, N represents the total length of the sequence, Δn represents the sequence delay length, m represents Spatial reconstruction dimensionality. In this example, the spatial reconstruction dimension m is 6.

(4) Determine the permutation entropy S(q) of I ₁ (q):

In the formula, g=1, 2, ..., k represents the type of serial number,

(5) Determine the remaining component r(t):

Set the threshold of S(q), when S(q) is lower than the threshold, judge I ₁ (t) as a non-abnormal signal, and remove it as the first IMF component from the original signal x(t), that is r(t)=x(t)-I ₁ (t), the residual component r(t) is obtained. In this case, set the threshold of S(q) to 0.6.

Repeat steps (1)-(5) for r(t) to obtain IMF components I(t)=I ₁ (t), I ₂ (t), . . . , I ₅ (t) in sequence.

(6) Calculate the Hilbert transform of I(t):

In the formula, Represents the Hilbert transform of I(t), I(t) represents the decomposed IMF component, t represents time, and τ represents the time interval.

(7) Determine the analytical signal g(t) of I(t) and perform autocorrelation processing:

In the formula, I(t) represents the decomposed IMF component, Represents the Hilbert transform of I(t), j represents the imaginary number unit, g(t) represents the analytical signal, y(t) represents the signal after autocorrelation processing, z(t) represents the autocorrelation signal template, and t represents time , τ represents the time interval.

(8) Discrete wavelet transform:

In the formula, α represents the scale factor, ψ(t) represents the mother wavelet function, τ represents the displacement information provided by the mother wavelet function, t represents time, y(t) represents the signal after autocorrelation processing, ξ _i (n) represents the discrete wavelet Transformed wavelet coefficients, i represents the current i-th frequency band, n represents the current n-th frequency band. Therefore, under different scales and displacements, multi-level wavelet decomposition is realized, and the schematic diagram of the process is shown in Figure 3. In this implementation case, the value of α is 0.03.

(9) Calculate the average energy on different frequency bands:

In the formula, ξ _i (n) represents the wavelet coefficient on each frequency band, i represents the current i-th frequency band, n represents the current n-th frequency band, and N represents the length of wavelet coefficients in each frequency band. In this case, the value of N is 10.

(10) Calculate the average energy ratio Ω _i in each frequency band:

In the formula, Ω _i represents the average energy ratio of each frequency band, E _sum represents the sum of energy of all frequency bands, and E _i represents the energy of the current frequency band.

The average energy ratio obtained from the highest frequency coefficient is not significantly different from other coefficients. In order to reduce the number of redundant features and eliminate the correlation between features, the energy ratio of the highest frequency coefficient is removed. Therefore, the finally obtained eigenvector is e={Ω ₁ , Ω ₂ , Ω ₃ , Ω ₄ , Ω ₅ } ^T .

In order to verify the effectiveness of the present invention for feature extraction of optical fiber vibration signals, an experiment for feature extraction of optical fiber vibration signals was carried out on the present invention. In order to visualize the differences between each feature, their energy ratio was calculated. The experimental results are shown in Figure 4. It can be seen from Figure 4 that the energy proportions of the wavelet coefficients of the four different optical fiber vibration signals are significantly different, and using the energy ratio as the eigenvector of the vibration signal has a clear degree of discrimination. This shows that the fiber optic vibration signal feature extraction algorithm established by the present invention is effective, provides a better vibration signal feature extraction method for establishing an accurate distributed fiber optic vibration signal classification model, and is more suitable for practical use.

Claims (1)

1. the present invention is characterized in that: (1) determine the vibration signal after introducing the white noise signal; (2) determine the first IMF component sequence set; (3) carry out time-delay space reconstruction; (4) determine the permutation entropy; ( 5) Determine the residual component; (6) Calculate the Hilbert transform of the sequence; (7) Determine the analytical signal of the sequence and perform autocorrelation processing; (8) Discrete wavelet transform; (9) Calculate the average energy on different frequency bands; (10) Calculate the average energy ratio on each frequency band, specifically including the following ten steps: Step 1: Determine the vibration signal after introducing the white noise signal and Introducing white noise signals n _p (t) and -n _p (t) into the vibration signal x(t), we get In the formula, and Indicates the vibration signal after introducing the white noise signal, x(t) indicates the original optical fiber vibration signal, n _p (t) and -n _p (t) indicate the white noise signal, a _p indicates the amplitude of the noise signal introduced for the pth time, p=1, 2,..., N _noise , N _noise represents the total number of noises introduced; Step 2: Determine the first IMF component sequence set I ₁ (t): right and Perform EMD decomposition separately to obtain the first set of component sequence sets and The components with the same subscript number in the sequence set are summed, accumulated, and averaged to obtain I ₁ (t): In the formula, I ₁ (t) represents the first IMF component sequence set, and Indicates the first set of IMF component sequence sets, p=1, 2, ..., N _noise , N _noise represents the total number of times noise is introduced; Step 3: Perform delay space reconstruction: Delay space reconstruction is performed on the digital sequence I ₁ (q) corresponding to I ₁ (t), and the following sequence is obtained: In the formula, In ( _n ) represents the sequence after delay space reconstruction, q=1, 2, ..., n, ..., N, N represents the total length of the sequence, Δn represents the sequence delay length, m represents Spatial reconstruction dimension; Step 4: Determine the permutation entropy S(q) of I ₁ (q): In the formula, g=1, 2, ..., k represents the type of sequence number, m is the spatial reconstruction dimension, Step 5: Determine the residual component r(t): Set the threshold of S(q), when S(q) is lower than the threshold, judge I ₁ (t) as a non-abnormal signal, and remove it as the first IMF component from the original signal x(t), that is r(t)=x(t)-I ₁ (t), get the residual component r(t); Repeat steps (1)-(5) for r(t), and then get the IMF components I(t)=I ₁ (t), I ₂ (t), ..., I _s (t), where s is the IMF the number of components; Step 6: Calculate the Hilbert transform of I(t): In the formula, Represents the Hilbert transform of I(t), I(t) represents the decomposed IMF component, t represents time, and τ represents the time interval; Step 7: Determine the analytical signal g(t) of I(t) and perform autocorrelation processing: In the formula, I(t) represents the decomposed IMF component, Represents the Hilbert transform of I(t), j represents the imaginary number unit, g(t) represents the analytical signal, y(t) represents the signal after autocorrelation processing, z(t) represents the autocorrelation signal template, and t represents time , τ represents the time interval; Step 8: Discrete wavelet transform: In the formula, α represents the scale factor, ψ(t) represents the mother wavelet function, τ represents the displacement information provided by the mother wavelet function, t represents time, y(t) represents the signal after autocorrelation processing, ξ _i (n) represents the discrete wavelet Transformed wavelet coefficients, i represents the current i-th frequency band, n represents the current n-th frequency band, therefore, under different scales and displacements, multi-level wavelet decomposition is realized; Step 9: Calculate the average energy on different frequency bands: In the formula, ξ _i (n) represents the wavelet coefficient on each frequency band, i represents the current i-th frequency band, n represents the current n-th frequency band, and N represents the length of the wavelet coefficient in each frequency band; Step 10: Calculate the average energy ratio Ω _i in each frequency band: In the formula, Ω _i represents the average energy ratio of each frequency band, E _sum represents the sum of energy of all frequency bands, and E _i represents the energy of the current frequency band; The average energy ratio obtained by the highest frequency coefficient is not significantly different from other coefficients. In order to reduce the number of redundant features and eliminate the correlation between features, the energy ratio of the highest frequency coefficient is removed. Therefore, the final eigenvector is e={Ω ₁ , Ω ₂ , Ω ₃ , Ω ₄ , Ω ₅ } ^T .