Distributed fiber-optic sensing (DFOS), as a sophisticated microdeformation-detection technique [1-4], is extensively utilized in the realms of oil pipeline leak testing, speed monitoring, and perimeter guarding [5-7]. Over the decades, various DFOS configurations, involving the Sagnac, Michelson, and Mach-Zehnder interferometers [8-11] and optical time-domain reflectometer (OTDR) [12], have been designed for different application conditions and sensitivity levels. In terms of backscattering light, the OTDR technique has also developed into many novel configurations, such as the phase-sensitive OTDR (phi-OTDR) [13-15], polarization OTDR [16], and Brillouin OTDR [17], which have the advantages of spatial resolution and response sensitivity, compared to other DFOSs based on interferometry [18-22]. Especially for phi-OTDR, the waveform of distributed vibrations can be computed by discerning differential phases in the sensing system [23-25]. Due to this prominent technique, multiple physical parameters can be measured simultaneously and exactly, which has led to the realization of vibration detection in an ultralong fiber [26]. Nevertheless, because of its high sensitivity, the backscattering light is susceptible to the surrounding environment, in practice. In addition, random noises such as electrical noise and phase noise of the laser can deteriorate the phi-OTDR signal and make vibration detection more difficult. Thus, for locating disturbances and improving the signal-to-noise ratio (SNR) of phi-OTDR, many studies have been carried out. Some methods have added amplitude-based or frequency-modulated optical systems to the phi-OTDR, such as distributed Raman amplification [26], chirped pulse amplification [27], first-order Raman amplification [28], frequency-sweep pulse modulation [29], and so on. These improvements depend on the high performance of the introduced optical hardware, and are more adapted to conditions where SNR decrease is caused by long distance. On the other hand, the algorithmic improvement of signal processing, as a practical means, is also effective for raising sensitivity and spatial resolution. Moving average [30], continuous wavelet transform [31], temporal signal separation and determination algorithm [32], differential phase method [33], and adaptive two-dimensional bilateral filtering algorithm [34] for enhancing the SNR of phi-OTDR have been presented over the past decade. Due to the nonlinearity and nonstationarity of the phi-OTDR signal, though, these classical methods based on various smoothing windows cannot make great advancements. Recently, the methodology associated with empirical mode decomposition (EMD) [35-38], as an effective means of processing nonlinear and nonstationary signals, was introduced to denoise the phi-OTDR signal. It can decompose complex noisy signals into a series of inherent mode functions (IMFs) that are suitable for extraction of nonlinear and nonstationary features, which allows the signal processing of phi-OTDR systems to enter a new epoch.

In our view, an effective denoising method for a phi-OTDR signal not only should be suitable for nonlinear and nonstationary input, but also should demonstrate the difference between local signal segments. The essential purpose of a denoising method is to enhance the useful signal associated with the vibration, and suppress the noise. Therefore, we propose a novel signal processing method based upon component multiplication, which is in terms of the nonlinearity and nonstationarity of IMF components [38], and is not regarded as a traditional denoising method that separates the noisy component from the original signal. For suppression of the mode mixing of IMF components, we design a variant EMD-based combination algorithm to break down the raw signal. To generate a new waveform that manifests the features of phi-OTDR signals more clearly, we design the multiplication function according to the EMD theory [38, 39]. Real experimental data are used to verify the proposed methodology.

Initially we need to decompose the original signal into a set of components that contain the characteristics of the original signal, in both frequency and time domains. To improve the effect of signal decomposition, we design a novel algorithm, termed extreme-mean complementary empirical mode decomposition (ECEMD). This method combines the ensemble empirical mode decomposition (EEMD) [40] with an extreme-mean-based envelope processing method. It can effectively settle the problems of redundant noise and low computational efficiency in the traditional EMD-based method, theoretically. As a part of the ECEMD algorithm, the extreme-mean based EMD method is depicted briefly as follows.

1) Identify the maxima and minima of the original signal.

2) Calculate the midpoints between two neighboring maxima and minima.

3) Interpolate between the computed midpoints, ending up with the envelope of extreme means.

4) Subtract the calculated envelope of extreme means from the original signal, to acquire the detail signal.

5) Regard the detail signal as the original signal, and iterate on the residual.

If the detail signal meets the following conditions, it can be defined as an IMF: a) The envelope fitted by the mean of local maxima and minima must be zero, or below a threshold. b) The numbers of extrema and zero-crossings must differ by no more than one. Following the steps above, the original signal can be broken down into a sum of n IMFs and a residual.

Introducing the extreme-mean-based EMD method, the arrangement for ECEMD can be shown in the following:

1) Adding two complementary (positive and negative) Gaussian white noises n(x) and −n(x) to the original signal, to yield new artificial signals presented as:

2) Decomposing the artificial signals x₋₁ and x₁ by the extreme-mean-based EMD method, we can acquire a set of IMFs including c_−1i and c_1i respectively, where i expresses the ordinal number of decomposed IMFs.

3) Repeating steps 1 and 2 n times, n pairs of IMFs c_−1i, c_1i, c_−2i, c_2i, ··· c_−Ni, , c_Ni can be obtained. Computing the ensemble means of these IMFs, the results are demonstrated as:

4) The final outcome is calculated as the mean of c_i and c_−i.

To compare different EMD-based approaches, we summarize their features in terms of noise addition and mean envelope. Theoretically speaking, the mean-envelope calculation is the kernel part of the EMD-based method, which determines its adaptability to the nonlinearity and nonstationarity of the original signal directly. Extreme-mean EMD improves the traditional EMD on the type of mean envelope, which can be regarded as using the envelope of extremum means to replace the mean of extremum envelopes in the procedure of traditional EMD. This suppresses mode mixing of IMF components effectively. On the other hand, noise addition, as an effective method of restraining mode mixing, is extensively applied in improvements of the traditional EMD-based algorithm, such as the EEMD method. Nevertheless, due to the extra noise, the SNR of some IMF components is sacrificed. To alleviate this issue, we introduce the regime that positive and negative noisy IMF components complement each other, which can neutralize the extra noise and reduce the computing time effectively. In addition, the complementary noisy IMF components are more beneficial for suppressing the strong noise by multiplying each other, because mode mixing is restrained effectively. As demonstrated in Table 1, the ECEMD method makes improvements in the aspects of mean envelope and noise addition that will lead to comprehensive improvement of mode-mixing suppression, as well as SNR enhancement.

According to the EMD theory [35-41], decomposed IMF components can be considered as a series of independent signal components, which all involve the essential segments of the original signal, but which do not all involve the noise segments of the one in the time domain [39, 40]. The essential segments are caused by the strong vibration, due to knocking, climbing and so on, while the noise segments are only caused by the noise. At one time point, only if the values of all useful IMFs are not zero can the signal point be considered essential; otherwise, it is considered to be noise, as shown in Fig. 1 [39-41]. In other words, the strong vibration can act on all of the IMF components, which results in the values of the corresponding point of all the IMF components being unequal to zero, but the noise cannot.

Following this point, we can enhance the SNR by multiplying the useful IMF components in the time domain, which can increase the SNR by several orders of magnitude. Detailed steps are organized as follows:

1) Calculate the correlation coefficient ρ_i for the original signal with each IMF component via Eq. (4),

2) Compute the test function t_i of the correlation coefficient: . Check the critical value t_α/2 in terms of degrees of freedom N − 2 and given significance α. If |t_i| ≥ t_α/2, the IMF components corresponding to ρ_i will be reserved, and termed as , , ···, . Otherwise, they will be eliminated.

3) Multiply the reserved IMF components to obtain the final artificial signal, which can be considered as the denoised signal, when the original non-noisy waveform mainly consists of a constant signal x_a(t), such as the phi-OTDR signal.

Reviewing the above steps that process the original phi-OTDR signal, we can acquire the flow chart for the proposed method, depicted in Fig. 2.

We simulate a single pulse to test the proposed method. After adding Gaussian white noise with SNR of 10 dB, the original waveform is completely unrecognizable in the artificial signal, as shown in Fig. 3(b). We perform the ECEMD method on the noisy signal, which generates the IMFs numbered from 1 to 10. Selecting the useful IMF components according to the correlation coefficient, we reserve IMFs 1-5. Multiplying them, the outcome can be plotted as Fig. 3(c), which indicates the location of the original pulse clearly and removes the noise at the other positions. Although the processed waveform is not a strictly denoised signal, the location of original pulse decorated with background noise can be identified clearly. Additionally, from the reserved IMFs (Fig. 4), it can be seen that all of them have a value unequal to zero at the position of the original pulse exactly, which powerfully verifies the EMD theory that only when the values of all useful IMFs are not zero is the signal segment considered essential. Finally, to compare our approach to a traditional denoising method, we apply the db2 wavelet denoising algorithm to process the artificial signal. In Fig. 3(d), it can be seen that the peak at 1.8 km confuses the vibration location seriously, and might be considered as a new vibration. In contrast, by our method the vibration point is positioned only at 5 km exactly, and hardly any confusing peaks exist. Moreover, the processed waveform is highly similar to the original one, which makes the SNR of phi-OTDR system enhance by several orders of magnitude.

To investigate the positioning performance of the proposed method for different SNRs, we mix the original signal (as shown in Fig. 3(a)) with continuously changing noise, in steps of 0.5 db. Our approach, the EMD-based method, the wavelet method, and the EEMD method are used to denoise the artificial signal and locate the vibration. Figure 5(a) demonstrates the location results of the four methods. One can observe that the proposed method shows a higher accuracy of vibration detection compared to the other approaches. With the SNR from 4.5 dB to 20 dB, our approach can position the pulse at 5 km exactly, while the EMD-based methods and conventional wavelet method cannot give accurate outcomes in cases with SNR below 13 dB. It verifies the effectiveness of the proposed approach for single-vibration detection. Moreover, in Fig. 5(b) we demonstrate the SNR after processing, with SNR from 4.5 to 10 dB, for the four approaches. Our method is again much superior to the other methods. One can see that even if the initial SNR is only 4.5 dB, after denoising it can exceed 35 dB, which is enough to give an exact position of the pulse.

Multiple-vibration detection is the advantage of the phi-OTDR system, which increases the demand for its robustness against noise. To test the ability of the proposed method to process multiple vibrations, we produce an original signal that consists of two pulses, as shown in Fig. 6(a). Gaussian white noise with SNR of 10 dB is again mixed with the original signal. Performing the proposed method on the artificial signal, in Fig. 6(c) we can acquire a denoised signal presenting only two spike pulses, which can indicate the positions and amplitudes of the two pulses exactly. As a comparison, the outcome of the wavelet algorithm is far inferior to that of our method. In Fig. 6(d), the noise from 3 km to 7 km is still preserved, which causes confusion as to the number and positions of the original pulses.

Being similar to the single-pulse detection, continuous noise with a SNR step of 0.5 dB is mixed with the original signal. The outcomes of the four methods are shown in Fig. 7. Because the original signal is more complex, our method cannot provide an accurate location in the case of SNR below 6 dB, but it still demonstrates more robustness against noise than do the other three methods. The two positions (3 km and 7 km) of the original pulses are detected accurately in most situations, which verifies the effectiveness of our method for multiple-vibration detection.

The experimental setup of a general phi-OTDR system is schematically illustrated in Fig. 8. A narrow-linewidth laser launches continuous-wave light with a linewidth of 100 Hz into the optical fiber. By means of the isolator, the reflected light can be effectively prevented from returning to the laser. Then, an optical coupler separates 1% of the light for a local reference, and 99% of the light for the probe. An acoustic-optic modulator (AOM) driven by a function generator modulates the light into a pulse chain with 1.8 kHz repetition rate and 100 ns pulse duration. An erbium-doped fiber amplifier (EDFA) amplifies the light pulse that is injected into an optical sensing fiber 3 km in length. The local light and backscattered light are combined by a coupler, and the coherent signal is detected by a balanced photodetector (BPD). We use a data acquisition card (DAQ) to sample the output signal at a rate of 100 MHz. To simulate the vibration, we use the cylindrical piezo transducer (PZT) as a vibration source, which can simulate vibrational frequencies from Hz to kHz.

Based upon the phi-OTDR system, we perform the experiment of vibration detection for verifying the effect of proposed method in a real environment. A vibration is implemented by the PTZ at 1750 m, which induces a raising pulse in the strong background noise at the corresponding position. Our method and wavelet algorithm are applied to process noisy signal. In Fig. 9, it can be observed that even in a real environment, the outcome by our approach also has an essential difference from the one by traditional wavelet method, the line pulse almost only exists at the position of the vibration, which is beneficial for identifying the vibration position. In contrast, the result using conventional wavelet method still reserves some noise from 0 m to 1500 m, which impacts the identification of the vibration pulse seriously. On the other hand, in the results of two method, compared with the spick pulse by the wavelet algorithm, the line one is more suitable for indicating the positions of the pulses in the horizontal coordinate. In addition, the signal pressing times of our method and wavelet algorithm are 4.2 s and 4.4 s respectively, which are very close.

To test the capability of detecting multiple vibrations in practical applications, we generate two vibrations at 1550 m and 2250 m in the phi-OTDR system. As shown in Fig. 10(a), the responses are submerged in the strong noise. The proposed method and the wavelet algorithm are employed to remove the noise. From Fig. 10(a), it can be seen that the noise seriously impacts the determination of the vibration positions, and at least five peaks can be suspected as pulses caused by vibrations, which directly induces the incorrect results of the wavelet algorithm. In Fig. 10(c) 970 m, 1800 m, and 2000 m are mistaken for positions of vibrations, and 2250 m is neglected; only 1550 m is a correct vibration position. However, with our method, although the noise impacts some IMF components seriously, the useful IMFs can also be employed to determine the vibration positions correctly. 1550 m and 2250 m are detected successfully, and the line pulses caused by strong noise are not enough to impact the ones induced by actual vibrations in Fig. 10(b). The signal processing times of our method and the wavelet algorithm are 4.8 and 4.7 s respectively. Undoubtedly, the proposed approach is very suitable for detecting multiple vibrations with the phi-OTDR system.

In this paper, we demonstrate a novel method to process a phi-OTDR signal mixed with strong background noise. Unlike conventional denoising methods for removing the noise from the original signal, our method aims to create an effective waveform signal that can indicate the positions of vibrational pulses and suppress strong noise. It is not a traditional denoising method based on subtraction. In this algorithm, we employ ECEMD to break down the noisy signal into IMF components, and generate the waveform signal by multiplying selected IMFs components. The simulations and experiments based on real data verify the validity of this method for detecting single or multiple vibrations. We believe that our approach can establish a new idea for denoising the raw signal from a phi-OTDR system.