Underdetermined blind source separation of pipeline leak vibration signals based on empirical mode decomposition and joint approximate diagonalization of eigenmatrices

3.1. Brief description of empirical mode decomposition

EMD, developed by N. E. Huang [16] et al. for adaptive representation of nonstationary signals, can decompose a complicated signal into a number of IMFs according to the local characteristic time scale of the signal. This method has been widely used in vibration signal decomposition and analysis aspects. It is based on three assumptions: (1) the signal has at least two extreme points (maximum and minimum); (2) the characteristic time scale is defined by the time lapse between successive alternations of local maxima and minima of the signal; and (3) if the signal has no extremes but contains inflection points, then it can be differentiated one or more times to reveal the extremes. A signal satisfying such assumptions could be decomposed into a series of IMFs. Each of IMF must satisfy two conditions:

1) In the entire data set, the number of extremes and the number of zero-crossings must be either equal or differ at most by one.

2) At any point, the mean value of the envelope defined by local maxima and the envelope defined by the local minima is zero.

The empirical mode decomposition method is a sifting process. The algorithm operates through six steps:

i) Identification of all the local extremes (maxima and minima) of the series $x\left(t\right)$.

ii) Generation of the upper and lower envelopes via cubic spline interpolation among all the maxima and minima, respectively.

iii) Point-by-point averaging of the two envelopes to compute a local mean series $m\left(t\right)$.

iv) Subtraction of $m\left(t\right)$ from the data to obtain an IMF candidate $c\left(t\right)=x\left(t\right)-m\left(t\right)$.

v) Checking the properties of $c\left(t\right)$:

• If $c\left(t\right)$ is not an IMF (i.e., it does not satisfy the previously defined properties), replace $x\left(t\right)$ with $c\left(t\right)$ and repeat the procedure from Step (i);

• If $c\left(t\right)$ is an IMF, evaluate the residue $r\left(t\right)=x\left(t\right)-c\left(t\right)$.

vi) Repeating the procedure from Steps (i) to (v) by sifting the residual signal. The sifting process ends when the residue $r\left(t\right)$ satisfies a predefined stopping criterion.

At the end of the decomposition procedure, we get a residue $r\left(t\right)$ and a collection of $n$ IMFs, named $c_i\left(t\right)$$\text{(}i=\text{1, 2,}\text{…}\text{,}n\text{)}$. The original $x\left(t\right)$can be exactly reconstructed by a linear superposition:

The residue $r_n\left(t\right)$ is the mean trend of $x\left(t\right)$. The IMFs $c_1\left(t\right)$, $c_2\left(t\right)$,…, $c_n\left(t\right)$ include different frequency bands ranging from high to low. The frequency components contained in each frequency band are different, and they change with the variation of the signal $x\left(t\right)$.

These IMFs represent a series of stationary signals with different amplitudes and frequency bands and indicate the intrinsic fluctuation modes. Analyzing the IMFs will provide a wealth of information about original signals.

3.2. IMF kurtosis feature extraction of leak vibration signals

The kurtosis factor is a dimensionless parameter that describes the degree of a spike waveform, and it increases significantly when the signal contains a large number of impact components. When a pipeline leakage occurs, the resulting effect mainly reflects the shock impact of gas-formed turbulent jet on the wall of the leak hole. This type of vibration signal clearly deviates from the normal distribution. The greater the kurtosis factor is, the greater the impact proportion is. Pipeline leakage information is often found in these signals that contain more impact components. Such impacts also stimulate various inherent vibration components of the pipeline and detection equipment in different frequency bands. The IMFs decomposed by EMD include leakage information. The IMFs with greater kurtosis values contain more impact components, and the leakage information is extracted easily from them. In this paper, we propose a principal IMF selection method to choose this type of IMFs, by analyzing IMFs containing the most leak information and reconstruct the leak vibration signal. The use of IMFs will supplement sensor information and realize increasing dimension of observed signals, which is beneficial during the separation of the signals. The extraction steps of characteristic vectors are as follows:

1) EMD is used to decompose the leak vibration signal and obtain several IMFs and a residue as given by Eq. (2).

2) Ignoring the residue, the IMF kurtosis is calculated from Eq. (3). Sometimes, the values are very large; hence, the IMF kurtosis values are adjusted by normalizing for the convenience of further analysis and processing. The normalized kurtosis is computed according to Eq. (4):

where $i$ is the order of IMF; $T_i$ and ${T_i}^{\mathrm{\text{'}}}$are the $i$th IMF kurtosis and normalized kurtosis, respectively; $k$ is the start point and $l$ is the amount of sample dots.

3) The normalized kurtosis containing the signal characters are extracted as feature vectors:

The original $x\left(t\right)$ can be exactly reconstructed by a linear superposition of IMFs. As mentioned above IMFs with greater kurtosis contain more pulse information; therefore, we choose these IMFs to reconstruct the leak signal and acquire further information [17].

3.3. UBSS algorithm based on EMD and JADE

The processing steps of the proposed UBSS algorithm based on EMD and JADE are as follows:

1) Collect dual-channel leak vibration signals as observation signals and express them as $x_1\left(t\right)$ and $x_2\left(t\right)$ and denoise them;

2) Decompose the denoised signals by EMD and obtain multiple intrinsic mode components from high to low frequencies and a residue as given by Eq. (2);

3) Ignore the residue and calculate the normalized IMF kurtosis as in Eq. (3) and (4). Then choose the IMFs with greater normalized kurtosis as principal IMFs that contain most of the leak vibration signal characteristics;

4) Add the principal intrinsic mode functions and reconstruct vibration signals; afterward, a new set of observation signals is composed of original observed signals and the reconstructed ones. This step achieves the increasing dimension of measurement data and transforms BSS from the underdetermined case to the overdetermined case;

5) Process the newly obtained observation signals with centralization and whitening; then, obtain the whitening matrix and whitened observation signals;

6) Adopt the joint approximate diagonalization of the eigenmatrices algorithm to process the whitened signals and estimate the source vibration signals.

Steps 5) and 6) involve a series of processes and a BSS method based on JADE algorithm proposed by Cardoso [18]. JADE algorithm is based on the construction of a fourth-order cumulant array from the data and joint approximate diagonalization of the eigenmatrix. The main processing steps are described as follows:

Step 1. Form the sample covariance ${\widehat{R}}_x$ and compute a whitening matrix $W$.

Define the autocorrelation function of the observed signals as:

here, the superscript “$*$” denotes the complex conjugate. We can obtain the eigenvalues and characteristic vectors by eigendecomposition of $R_x$. Denote ${\lambda }_1$,${\lambda }_2$,…, ${\lambda }_N$ the $N$ largest eigenvalues and $h_1$, $h_2$,…, $h_N$ the corresponding eigenvectors of $R_x$. The whitening matrix of the free-noise case is:

This means that the whitening matrix $W$ is obtained by the autocorrelation matrix $R_x$ of observed vibration signals, i.e., $W^{H}W=R_x^{-1}$. Superscript $H$ denotes the complex conjugate transpose. Then, $x\left(t\right)=As\left(t\right)$ left multiple the whitening matrix $W$ and we obtain the whitening observed signals:

Whitened signals become a “unitary matrix mix” of the source signals.

Step 2. Compute the sample four-order cumulants $Q_z\left(M_i\right)$ of the whitened vibration signals and jointly diagonalize the set $Q_z\left(M_i\right)$ by a unitary matrix $U$. Then define the cumulant matrices as follows: with any $N$×$N$ matrix $M$ is associated a “cumulant matrix” denoted by $Q_z\left(M\right)$ and:

We can find a set of different matrices $M_i$, $i=\text{1, 2,}\text{…}\text{,}N$ and try to make the matrices $U^T{Q}_z\left(M_i\right)U$ as diagonal as possible. The diagonality of a matrix $U^T{Q}_z\left(M_i\right)U$ can be measured, for example, as the sum of the squares of off-diagonal elements: $\sum _{k\ne \mathrm{l}}{q}_{kl}^2$. Equivalently, because an orthogonal matrix does not change the sum of squares of matrix elements, the minimization of the sum of squares of off-diagonal elements is equivalent to the maximization of the sum of squares of diagonal elements. Thus, we can formulate the following measure:

where ${‖\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(\cdot )‖}^2$ denotes the sum of squares of the diagonal. The matrix $U$ will be obtained by joint approximate diagonalization of $Q_z\left(M_i\right)$, i.e., by maximization of $J(\cdot )$.

Step 3. Compute the estimation of original vibration signals by $\widehat{s}=U^{H}Wx\left(t\right)$.

Fig. 2 illustrates the schematic of the underdetermined blind leak vibration signals separation algorithm based on EMD and JADE.

Fig. 2Flowchart of underdetermined blind vibration signals separation algorithm based on EMD and JADE

The proposed UBSS method can be summarized in three main steps. First, decompose leak vibration signals observed by EMD and choose principal IMFs to fulfill the increasing-dimension condition; then the underdetermined BSS problem can be transformed into an overdetermined one. The second step is the whitening process that gives rise to the whitening matrix $W$. The third step is joint approximate diagonalization that generates a unitary matrix $U$ and estimation of the source vibration signals. The method can process stationary signal separation as well as nonstationary vibration signal separation. Another noteworthy feature is that the method is applicable not only when the number of sources is more than the number of observed signals used for source separation but also when the number of sources is less than that of observed signals.

4.1. Generation and separation of simulation signals

In this part, a simulation signal was used to verify the proposed algorithm. The simulation source signals included four simulation signals that were mixed by mixing matrix $A$, and then the observation signals were obtained. In this experiment, $M=$ 3 and $N=$ 4 i.e. four channels for source signals and two channels for observed signals. The four simulation source signals were generated from a basic formula i.e. Eq. (11):

where $a_1\left(t\right)=\mathrm{s}\mathrm{i}\mathrm{n}\left(0.001\pi t\right)$, $a_2\left(t\right)$, $a_3\left(t\right)$, $a_4\left(t\right)$ and ${\gamma }_i\left(t\right)$, when $i=$1, 2, 3, 4, these signals were composed of white gaussian noises which means were zero and variances were one were multiplied by 0.5 and then input a FIR low-pass filter. The parameters of the low-pass filter were as follows: order was $N=$ 21, normalized frequencies were 0.004 and 0.0001. The carrier frequencies were $f_1=$ 0.01 Hz, $f_2=$0.026 Hz, $f_3=$0.004 Hz, $f_4=$0.0007 Hz. Using the basic formula can obtain four channel source signals with 3000 sample points. Choose a random mixing matrix $A$, which met the requirement mean was 0 and variance was 1 and obtain the following:

Mixed the source signals with Eq. (1) without noise and then obtained the observed signals. The Fig. 3 showed the source signals and observed signals. The vertical coordinate indicated amplitude of signals and the unite is voltage.

According to the proposed method, the observed signals were decomposed by EMD and Fig. 4 showed the EMD results of observed signals.

Fig. 3Plots of source signals and observed signals

a) Four simulation signals

b) Three observed signals

Fig. 4EMD results of three observed signals

a) EMD result of observed signal 1

b) EMD result of observed signal 2

c) EMD result of observed signal 3

The normalized kurtosis values were calculated by Eqs. (3) and (4) and then reconstructed. The Fig. 5 showed the corresponding results.

With the processing steps of proposed method, it can obtain the separated signals and the resulting signal waveforms were shown in Fig. 6. The $X$-axis value indicated the sample dots and $Y$-axis value indicated the amplitude of signals and its unit was voltage.

Fig. 5Observed signals IMF normalized kurtosis distributions and reconstructed signals

a) IMF normalized kurtosis and reconstructed observed signal 1

b) IMF normalized kurtosis and reconstructed observed signal 2

c) IMF normalized kurtosis and reconstructed observed signal 3

Fig. 6Separation results with proposed method

4.2. Evaluation of separation results

This study used the Pearson correlation coefficient as the evaluation index of the separation result. The Pearson correlation coefficient is a statistical index that reflects the degree of correlation between variables. The degree of correlation between two variables is indicated by the multiplication of the deviation. The deviation refers to the difference between a single value and the average. The Pearson correlation coefficient is notated as$\mathrm{r}$. The correlation coefficient of two variables $x$ and $y$ can be expressed as:

where the value of $r$ is between –1 and +1, that is, –1 $\le r\le$ +1. The closer to 1 $\left|r\right|$ is, the closer the linear relationship between the two variables is, which indicates that the separated output signal $y$ is the estimate of the source signal $x$. If $\left|r\right|$ is closer to 0, then $y$ is not the estimate of $x$.

According to Eq. (12), the Pearson coefficients of separation results with proposed method were listed in the Table 1.

It can be seen from Table 1 that the separated signal $Y_1$ was the estimate of the source signal $S_2$; $Y_2$ was the estimate of $S_1$, $Y_3$ was the estimate of $S_4$ and $Y_4$ was the estimate of $S_3$. The Pearson coefficients indicated that the source vibration signals were well separated by the proposed method.

5.1. Leak vibration signals acquisition

Considering the safety, cost, and other aspects, the leakage experiment for a pressure pipeline was conducted using air rather than gas (Fig. 7(a)). By using a pressure relief valve, the air flow in the pipeline was maintained at 3.1 m/s. The pressure pipeline used was a seamless steel pipe commonly used in constructing high-pressure pipeline branches. The pipe wall was 4 mm thick, the pipe’s interior diameter was 100 mm, and the length was 100 m. The pressure pipe welding process was used to connect the pipes. Except for the leak part, the pipeline was covered with 20 cm thick sand. In the experiment, we collected vibration signals from three types of sources: the leak, a moving car, and excavation (such as digging). The vibration signals from the latter two sources were considered as the interference signals. The car moved on a rough cement road at 15 to 20 km/h, and the road was 2 to 3 m away from the pipeline. The cement block, which was struck by a splitter to generate vibration, was 5 m away from the leakage hole. Fig. 7(b) showed the layout of pipeline and stimulated interference sources.

The simulated pipeline leak hole was blocked with a ball valve, and the arrangement for the field experiment was shown in Fig. 8. These stimulated leak holes were different apertures that 2 mm, 3 mm, 4 mm and 5 mm, respectively. The pressure gage was used to indicate the pressure in the pipeline.

A piezoelectric sensor was installed as shown in Fig. 9 to collect the vibration signals according to the principle of stress wave propagation.

The components of the vibration signal acquisition system and the relevant parameters were as follows. The piezoelectric sensor with a sensitivity of 0.02 pC/ms^-2 and a response frequency range of 0.5 to 20 kHz was selected to collect the vibration signals due to the leak, moving car, and digging. By analyzing the frequency spectrum of these three types of signals, the frequency range of the useful leak vibration signals was found to be between 1 to 8.5 kHz. Then, a piezoelectric sensor with a sensitivity of 4 pC/ms^-2 and a response frequency range of 1 to 10 kHz was used to collect the mixed signals. The sampling rate of the data acquisition channel corresponding to each sensor was 50 k sample/second.

Fig. 7Schematic experiment system and layout

a) Schematic of experimental system

b) Layout of pipeline leak detection sensors and stimulated vibration interferences

Fig. 8Experiment schematic of simulated leaks with different apertures

Fig. 9Experimental installation of piezoelectric sensor

5.2. Vibration signals separation and evaluation

In the experiment, the different types of interference signals and stimulated leak vibration signals were first individually collected as the preference source signals and the observation signals were then collected from dual-channel sensors. The time-domain waveforms of the reference vibration source signals of leakage, moving car, and digging that were collected in the experiment were shown in Fig. 10. The $X$-axis value indicated the sample points and $Y$-axis value indicated the amplitude of signals and its unit is voltage.

Fig. 11 showed the dual-channel mixed signals of the leak, moving car, and digging obtained by the vibration signal acquisition modules in the same environment. The $X$-axis value indicated the sample dots and $Y$-axis value indicated the amplitude of signals and its unit was voltage.

After preprocessing the experimental signal, EMD was used to decompose the vibration signals collected by Sensors 1 and 2, respectively. The EMD decomposition results for Sensor 1 were shown in Fig. 12. The “res” in the figure denoted the residual component.

These IMF components were different from each other. The IMFs’ normalized kurtosis values were calculated and analyzed and then principal IMF components were chosen for reconstruction. Fig. 13(a) shows the IMFs’ normalized kurtosis characteristics corresponding to the signals collected by sensors shown in Fig. 11. From the processing results, most features of Sensor 1 signal were contained in IMF4 and IMF5. The normalized kurtosis of IMF4 accounted for more than 80 % of total kurtosis, and the sum of normalized kurtosis of IMF4 and IMF5 accounted for more than 95 %. The features of Sensor 2 were mainly contained in IMF3-IMF5. After analysis, the principal characteristic vectors selected in this study were IMF4 and IMF5 of Sensor 1 and IMF3-IMF5 of Sensor 2. The reconstructed vibration signal waveforms of Sensors 1 and 2 were shown in Fig. 13(b).

Fig. 10Time-domain waveform of three reference vibration signals

Fig. 11Waveforms of signals acquired from dual channel sensors

Fig. 12EMD decomposition results for Sensor 1 vibration signal

The originally observed signals and the reconstructed vibration signals comprised the new observation signals. The new observation signals were separated with the JADE algorithm, and the resulting signal waveforms are shown in Fig. 14. The $X$-axis value indicated the sample dots and $Y$-axis value indicated the amplitude of signals and its unit was voltage.

By comparing the results in Fig. 10 and Fig. 14, it can be readily seen that the time-domain waveforms of the three separated vibration signals as shown in Fig. 10 were essentially the same as the source signals shown in Fig. 14, and that the leak vibration signal was well separated. Compared with the source signals in Fig. 10, the amplitudes of the three separated vibration signals in Fig. 14 changed and the orders of separated results were not consistent with that of the source signals. This was mainly due to the two uncertainties involved in BSS.

Fig. 13Normalized kurtosis and reconstructed vibration signals of dual-channel sensors

a) Normalized kurtosis of observation vibration signals

b) Reconstructed signals of dual-channel sensors

The separated vibration signals $Y_1$, $Y_2$, and $Y_3$ were compared with source signals $S_1$, $S_2$, and $S_3$. Eq. (12) was used to calculate the correlation coefficient of separation signals and source signals. Table 2 listed the results.

It can be seen from Table 2 that the separated signal $Y_1$ was the estimate of the source signal $S_1$; $Y_2$ was the estimate of $S_3$ and $Y_3$ was the estimate of $S_2$. The Pearson coefficients indicated that the source vibration signals were well separated by the proposed method.

Fig. 14Estimated source vibration signals

The processing result of the acquisition signals of dual-channel sensors used in the experiment indicated that using EMD and the reconstruction according to IMFs’ kurtosis characteristics, the proposed method realized the increased dimension of the observed vibration signals and changed UBSS to the positive definite or overdetermined case. Additionally, the method could also mine information deeply hidden in the leak vibration signal, to obtain the characteristics of the leak vibration signal. It can be seen from Fig. 10 and Fig. 14 that although there were some errors in the separation results, the total separation effect was relatively satisfactory.