2.1.

Time-Frequency Analysis

Time-frequency diagrams obtained by time-frequency analysis show the combined information in time domain and frequency domain, which directly reflect change of the frequency components of the signals with time. Common time-frequency transform methods include short-time Fourier transform (STFT), continuous wavelet transform, and S-transform. Because spectrograms obtained by STFT have succeeded in practical application in the fields of music, sonar, radar,¹² and speech processing¹³ and has a rich image color information that is helpful to image recognition, we choose STFT as the time-frequency analysis method in this paper.

2.1.2.

Spectrogram

After the signal in time domain in each segment is transformed to the frequency domain by fast Fourier transform (FFT), the amplitude frequency is rotated by 90 deg, and the spectral amplitude mapped to a gray level (0 to 255) value. The higher the amplitude is, the darker the corresponding region. Time-frequency representation of a signal is referred to as spectrogram.

Due to the randomness of the vibration signal and the complexity of its frequency spectrum performance, we use the Hamming window function to segment the signal into the small frames. The Hamming window function can reduce the discontinuities of the signal at the edge of frame and eliminate high-frequency interference and energy leaks. It is formulated with the following equation:

Eq. (2)

$$w(n)=0.54-0.46\cos [2n\pi /(N-1)],0\le n\le N-1,$$

where $N$ is the length of the window.

By using high-resolution FFT to process the vibration signal, we found that the lowest frequency of its spectrum peak was around 25 Hz. Therefore, the frequency resolution should be lower than 25 Hz. With a system sampling rate of 5000 Hz, the maximum value of the number of FFT could be calculated as 200. To obtain a high-frequency resolution and facilitate the computing, the number of FFT is set as 128 and the time length of each frame is 64 ms.

Take a digging signal as an example. The total time length of the signal is 2.4 s. After using STFT to transform the time-domain signal to the frequency domain, the time-frequency representation of the digging signal is shown in Fig. 1.

Fig. 1

An example of a knocking signal: (a) the original time-domain waveform of the knocking signal and (b) the spectrogram of the knocking signal.

2.1.3.

Spectral subtraction

Spectral subtraction¹⁴ is the most common method for dealing with wide-band noise. Its principle is obtaining the spectrum of the pure speech by subtracting the estimated noise amplitude spectrum from the amplitude spectrum of the noisy signal. The effect is equivalent to carrying out a type of equalization to the noisy signal in the transform domain. Compared with other noise reduction methods, spectral subtraction has the least constraint, the most direct physical meaning, the least amount of calculation, and the desired effect.

Assuming that the correlation between noise in the noisy signal and the vibration signal itself is zero, that is, the two are independent of each other and additive in the frequency domain, the additive model of the signal can be expressed as

Eq. (3)

$$X_w(\omega )=S_w(\omega )+N_w(\omega ),$$

where $X_w(\omega )$ is the short-time Fourier spectra of the small frame of noisy signal, $S_w(\omega )$ and $N_w(\omega )$ are the short-time Fourier spectra of noise component and effective component of the signal frame, respectively.

The magnitude squared of the short-time Fourier spectrum gives the short-time power spectrum

Eq. (4)

$${|Y_w(\omega )|}^2={|S_w(\omega )|}^2+{|N_w(\omega )|}^2+S_w(\omega )N_w^{*}(\omega )+S_w^{*}(\omega )N_w(\omega ).$$

Since $s(m)$ and $n(m)$ are independent of each other, the cross-correlation statistical mean is 0. The short-time power spectrum of the pure signal is estimated as

Eq. (5)

$${|{\widehat{S}}_{\omega }(\omega )|}^2={|Y_{\omega }(\omega )|}^2-E[{|N_{\omega }(\omega )|}^2],$$

where $S_{\omega }(\omega )$ is the short-time power spectrum estimation of the pure signal.

The phase is restored and then the inverse Fourier transform is used to obtain the denoised time-domain signal

Eq. (6)

$$\widehat{s}(m)=\mathrm{IFFT}[|{\widehat{s}}_w(\omega )|e^{i\phi (\omega )}].$$

Since the noise is locally stationary, to estimate the noise power spectrum, it is assumed that the noise at no vibration is close to the noise at the time of vibration. Thus, the noise power spectrum at the time of vibration can be estimated from the short-time power spectrum of the silent frame without vibration according to Welch’s method.

After using spectral subtraction to denoise the original knocking signal in Fig. 1, the time-frequency presentation of the denoised signal is shown in Fig. 2.

Fig. 2

The effect of spectral subtraction on the vibration signal: (a) the time-domain waveform of the knocking signal after noise reduction and (b) the spectrogram of the knocking signal after noise reduction.

It can be seen from Fig. 2(a) that the signal quality improves after noise reduction, and meanwhile the signal intensity and time-domain characteristics of the effective signals are not weakened. Figure 2(b) shows the broadband noise in noisy signals is significantly suppressed and other noise with multiple harmonics generated by the acquisition process is eliminated. By this method, adverse effects from noise in vibration signals on the spectrogram classification are reduced.

2.2.

Convolutional Neural Network

CNN is a hierarchical neural network composed of a sequence of layers. A typical model usually consists of several convolutional layers where image contents are represented as a set of feature maps obtained after convolving the input with a variety of filters, which are learned during the training phase. Pooling layers are introduced after convolutional layers to reduce the image size and accumulate maximum activation features from convolutional feature maps. Furthermore, CNNs may also contain fully connected (FC) layers where each neuron of the input layer is connected with every neuron in the layer. A sequence of convolutional, pooling, and FC layers form a feature extraction pipeline that models the input data in abstract form. Finally, a soft-max layer performs the final classification task based on this representation.

The proposed CNN model, shown in Fig. 3, has a structure similar to that used in AlexNet. The network contains five convolutional layers, two FC layers, and a soft-max layer. Each convolutional layer is followed by an ReLU and a max-pooling layer. The ReLU is a nonlinear function expressed as $f(x)=\max (0,x)$. It is applied as an activation function after the convolutions. Max-pooling following the ReLU passes on the maximum value in each $2\times 2$ block. The first convolutional layer takes the $224\times 224\times 3$ image and applies $1285\times 5$ filters. Through an ReLU and a $2\times 2$ max-pooling layer, the volume of the original image becomes $110\times 110\times 128$. The next two convolutional layers both have $1283\times 3$ filters and are followed by an ReLU and $2\times 2$ max-pooling, resulting in a $26\times 26\times 128$ image volume. The resulting volume is fed into the following two layers, each of which consists of $643\times 3$ filters, an ReLU, and a $2\times 2$ max-pooling layer. After all convolution and pooling, the $4\times 4\times 64$ image volume is flattened into a $1024\times 1$ vector. The first and second FC layers reduce the size of the feature vector to 256 and 4 in turn, and then the soft-max classifier takes the $4\times 1$ vector and outputs the final result.

Fig. 3

Visual depiction of the Image Network architecture (the soft-max classifier is not pictured).

2.3.

Support Vector Machine

SVM¹⁵ has advantages in the small sample, nonlinear, and high-dimensional expression and is widely used in pattern recognition. Its basic idea is that there is always an optimal hyperplane to separate a binary dataset through the support vector. Assuming that the given $N$-classes training sample set is inseparable in the low-dimensional space, it is expressed as $(x_i,y_i)$, where $i=\mathrm{1,2},\dots ,n$, $n$ is the training sample number, $x_i$ represents the input training samples, and $y_i$ represents the class markers of data samples. The input samples are projected using nonlinear mapping onto a high-dimensional space $R^n$ to construct an optimal hyperplane as

In the binary classification problem, classification results are determined by the symbol value after data samples to be measured are input to the classification function. To deal with $N$-classes classification problems (where $N>2$), the combinatorial classifier is often used. This type of classifier is a combination of all possible binary subclassifiers, which are $N(N-1)/2$ in total. Subsequently, the voting method is used to determine the class with the most votes, which is used as the class of the input data sample.

4.1.

Algorithm Execution Process

The flowchart of vibration signal identification in distributed optical fiber sensing system is shown in Fig. 5.

Fig. 5

The classification algorithm flowchart.

First, a differential operation and a normalization operation are performed on each data sample as a preprocessing step to obtain the normalized vibration waveform without DC offset. Subsequently, we apply spectral subtraction to reduce the noise in the vibration signals. The denoised signals are transformed into spectrograms by STFT. And then all spectrograms are randomly divided into a training set and a testing set, the proportion of which are 75% and 25%, respectively. Next, the proposed CNN is trained on the training set. In the training process, spectrograms are flattened into feature vectors after all convolutional, pooling, and FC layers. There are three following methods used for classification:

4.2.

Vibration Event Simulation Experiments

A segment of a 300-m-long optical fiber, led out in the vibration position of 20 km in the sensing fiber, was buried in the shallow soil of 10-cm depth. Striking the ground above the buried sensing fiber with shovels simulated digging events. Walking on the same position simulated walking events. Driving a car through this position simulated vehicle-passing events. Shaking the sensing fiber exposed on the ground simulated fiber-damage events.

Two hundred trials were carried out in each of the four vibration modes in the vibration position and around the vibration position, and data samples collected in the nearest five consecutive position nodes were stored. Through this method, 4000 data samples were obtained after all trials. Considering the fact that the system should respond to vibration events within 3 s, the duration and the total sampling point number of a single data sample were set as 2.4 s and 4800, respectively.

4.3.

Time-Frequency Analysis

We used spectral subtraction to reduce the noise in the vibration signals with the estimated noise power spectrums that were estimated by the average power spectrums of the optical signals within the last 0.6 s before the occurrence of vibration events in the corresponding position nodes. Then the denoised vibration signal waveforms were transformed into spectrograms by STFT. During the STFT, both the length of window function and the number of FFT point were 128, and the overlapping rate was 50%. The resolution of spectrograms was $224\times 224$, and the transformation process was encoded in Python and takes 0.3 s for each spectrogram. The typical spectrograms of all types of vibration signals after noise reduction are shown in Fig. 6.

Fig. 6

The typical spectrograms of all types of denoised vibration signals: (a) the digging signal, (b) the walking signal, (c) the vehicle-passing signal, and (d) the damaging signal.

4.4.

Network Training

Two separate datasets comprised of 4000 spectrograms from original signals and 4000 spectrograms from denoised signals, respectively. Each dataset was assigned into a training set and a testing set. For each vibration mode, 750 spectrograms generated in 150 trials were selected as part of the training set and other 250 sets were taken as testing data.

We used training sets from two datasets to train two CNNs. The training processes were performed in Caffe, running on NVIDIA GTX TITAN X Pascal GPU with 12-GB onboard memory. Each training process was an iterative process. With every iteration, the loss and accuracy were obtained, and the hyperparameters were changed to optimize for the next training. We set the batch size, the dropout value, and the initial learning rate as 100, 0.5, and 0.001, respectively, and run the training process for 40 epochs. Each network was tested on the testing phase every 20 iterations. The training took around 40 min on each dataset.

Figure 7 shows the train loss and test accuracy curves of two CNNs trained by different datasets. The CNN1 was trained and tested on the dataset in which spectrogram samples were transformed from original vibration signals, while the CNN2 was trained and tested on spectrograms of denoised vibration signals. In CNN1 training process, the training loss, which represented the network error, decreased slowly until hovering around 0.65 after 25 epochs, and the best testing accuracy was only 55%. On the contrary, the test accuracy of the CCN2 grew up to 88%, and a training loss of 0.30 was achieved after 30 epochs. The figure indicates that the noise reduction process on vibration signals based on the spectral subtraction improved remarkably the classification performance of the CNN2.

Fig. 7

The train loss and test accuracy curves for CNN1 and CNN2, which were trained by original spectrograms and enhanced spectrograms after noise reduction, respectively.

5.1.

Classification Based on Convolutional Neural Network

The CNN2 trained on spectrograms generated from denoised vibration signals was proved to have a good performance on vibration event recognition for distributed optical fiber sensing systems. The results of the CNN2 on the test set are shown in Table 1.

Table 1

Recognition accuracy of vibration signals using the CNN2.

Each row in Table 1 represents the recognition outcomes of each type of vibration signals on the testing set from the CNN2. It can be seen that the CNN2 classified four types of vibration signals with 97.2%, 96.0%, 74.0%, 84.8% accuracy on the testing set, respectively. These results indicate that this type of CNNs can be trained to classify very different types of vibration events, especially distinguish between instantaneous effects (digging and walking) and long-time effects (vehicle passing and damaging), but have a difficulty in accurately recognizing the vehicle signals. Though using a CNN is a more complicated way of carrying out this classification, these results demonstrate that this method successfully yields better results than other methods^8–9 and avoids manual feature selection that takes much more time in new application scenarios.

5.2.

Classification Based on Convolutional Neural Network and Support Vector Machine

In the RCNN proposed by Girshick et al.,¹⁶ an SVM classifier is applied to replace the soft-max classifier in the CNN. Girshick et al. extracted a 4096-dimensional feature vector extracted from each mean-subtracted $227\times 227$ RGB image using a CNN containing five convolutional layers and two connected layers, and then performed classification with an SVM classifier comprised of a set of category-specific linear SVMs. This method improved mean average precision by more than 30% compared to the performance of traditional CNN approach on object recognition. It have been proved that an SVM has a better performance than a soft-max classifier, especially in solving problems of small sample size.

We first extracted the output vectors of the last FC layer in the CNN2 as feature vectors in the SVM classifier. Considering that each feature vector only had a size of $4\times 1$, we used the radial basis function as the kernel function in the SVM classifier. The classification results are shown in Table 2. From Tables 1 and 2, we can see that classification results of the soft-max classifier and the nonlinear SVM classifier were approximately the same, except that the recognition rate of vehicle signals rose up to 79.2%. This was because the outputs of the last FC layer in a CNN were highly abstracted and reflected the classification results to a great degree. Thus, it was inefficient to replace the soft-max classifier with a complex classifier.

Table 2

Recognition accuracy of vibration signals using a nonlinear SVM classifier.

Then, we extracted the input vectors of the first FC layer in the CNN2 as feature vectors in the SVM classifier. The number of feature parameters was 1600, more than the sample number of each type of vibration events. Thus, we chose the linear kernel function as the SVM’s kernel function. Table 3 shows the classification performance of the linear SVM classifier. It indicates that this classification method achieved a classification accuracy of 93.3%, much better than other methods. Meanwhile, the successful application of a linear SVM proved that the input vectors of the first FC layer had a good expression on spectrogram features. And we only needed to optimize the parameter $C$ when using a linear SVM.

Table 3

Recognition accuracy of vibration signals using a linear SVM classifier.