Determination of the minimum distance between vibration source and fibre under existing optical vibration signals: a study

A signal emitting from a remote source and monitored in the presence of noise at two spatially separated receivers can be mathematically modelled as (18) {x1(n)=s(n)+n1(n)x2(n)=As(n−D)+n2(n) \left\{ {\matrix{ {{x_1}\left( n \right) = s\left( n \right) + {n_1}\left( n \right)} \hfill \cr {{x_2}\left( n \right) = As\left( {n - D} \right) + {n_2}\left( n \right)} \hfill \cr } } \right. where s(n) is the vibration source signal, x₁ (n) and x₂ (n) are vibration signals received by two receivers; A represents the attenuation coefficient; D is the TDOA between x₁ (n) and x₂ (n); n₁ (n) and n₂ (n) refer to zero-mean additive Gaussian noises, assuming that there is uncorrelation between noises.

GCC is the improvement on the basic cross-correlation method, and it can be used under the low SNR. GCC's principle is pre-filtering to eliminate noises and disturbance and then carrying out correlation operations.

According to the different weighting methods, GCC can be divided into four types, including Roth, SCOT, PHAT and ML. The ML algorithm can give the big weight to the frequency bands with the high SNR and the small weight to the bands with the small SNR, therefore it is the optimal filter statistically.

The generalised correlation between x₁ (n) and x₂ (n) is (19) Rx1x2(τ)=∫−∞∞ψg(f)Gx1x2(f)ej2πfτdf {R_{{x_1}{x_2}}}\left( \tau \right) = \int_{ - \infty }^\infty {{\psi _g}\left( f \right){G_{{x_1}{x_2}}}\left( f \right){e^{j2\pi f\tau }}df} For ML algorithm, ψ_g (f) is given (20) ψg(f)=1|Gx1x2(f)||γ12(f)|2[1−|γ12(f)|2] {\psi _g}\left( f \right) = {1 \over {\left| {{G_{{x_1}{x_2}}}\left( f \right)} \right|}}{{{{\left| {{\gamma _{12}}\left( f \right)} \right|}^2}} \over {\left[ {1 - {{\left| {{\gamma _{12}}\left( f \right)} \right|}^2}} \right]}} where |γ₁₂ (f)|² is the coherence function between x₁ (n) and x₂ (n).

The process of GCCSC is that the signals are autocorrelated and cross-correlated firstly, and then the results of autocorrelation and cross-correlation are used to carry out cross-correlation operation again, so as to improve resolution power and anti-noise performance.

Autocorrelation function of x₁ (n) is (21) Rx1x1(τ)=Rss(τ)+Rn1s(τ)+Rsn1(τ)+Rn1n2(τ). {R_{{x_1}{x_1}}}\left( \tau \right) = {R_{ss}}\left( \tau \right) + {R_{{n_1}s}}\left( \tau \right) + {R_{s{n_1}}}\left( \tau \right) + {R_{{n_1}{n_2}}}\left( \tau \right). The CC function between R_x1x2 (τ) and R_x1x1 (τ) is (22) RRR(τ)=E[Rx1x1(n)Rx1x2(n−τ)] {R_{RR}}\left( \tau \right) = E\left[ {{R_{{x_1}{x_1}}}\left( n \right){R_{{x_1}{x_2}}}\left( {n - \tau } \right)} \right] where R_x1x2 (n − τ) is calculated by Eq. (19). Ignoring cross-correlation function between noises and signals, Eq. (22) is simplified as (23) RRR(τ)=RRN(τ)+RRS(τ−D). {R_{RR}}\left( \tau \right) = {R_{RN}}\left( \tau \right) + {R_{RS}}\left( {\tau - D} \right). R_RN (τ) is secondary relation for noises, and R_RS (τ − D) is secondary relation for source signals. Considering R_RN (τ) = 0, there is maximum as τ in R_RR (τ) is D. D is the estimation value for TDOA.

To improve the estimation accuracy, two aspects of the secondary relation are improved: one is limiting the range of delay peak; the other is using three-point interpolation to improve the delay estimate.

Judging whether a data point is an outlier, it is necessary to select a suitable indicator. In this paper, the Deviation Square (DS) is used to judge whether Δt_i is an outlier. The definition is (24) DSi=(i2d2−Δti2v¯2−2Δtiv¯l0¯)2, D{S_i} = {\left( {{i^2}{d^2} - \Delta t_i^2{{\bar v}^2} - 2\Delta {t_i}\bar v\overline {{l_0}} } \right)^2}, where v¯ {\rm{\bar v}} and l0¯ \overline {{{\rm{l}}_0}} refer to the estimation values of wave velocity and distance calculated by Eq. (15). Figure 2 shows the geometrical intuitive significance of DS.

Fig. 2

According to Figure 2(a), when i is −1, 1, −6 and 6, the four DSi is relatively large. All DSi in Figure 2(a) are calculated and shown in Figure 2(b). In Figure 2(b), the four DSi marked with solid points are apparent outliers, and the rest 11 hollow points are more likely to be normal.

In Figure 2(b), the DS value of outliers is higher than normal points. The normal points gather in a specific region with smaller DS values, and the number is more. However, the outliers cluster into one or more classes in the region with larger DS values, and the number is less. Thus DS can be used to distinguish the normal points from outliers effectively.

This paper adopts the idea of eliminating outliers by clustering. The primary thought is as follows: First, all points are clustered into multiple classes according to the DS value. Second, the most likely outliers are determined by judging the number of points in the class and the distance between the classes. Third, and ared after removal of such points. Fourth, the new DS value is calculated. The above steps are repeated until the remaining data can be considered as one class. That is, no apparent outliers.

There are multiple clustering methods, including K-Means, hierarchical clustering, Clustering with Gaussian Mixture Models (GMM) and DBSCAN. However, this paper requires that only one class remain in the end, no matter how many classes there are at the beginning. K-Means and GMM cannot meet the requirement of changeable clustering number, and hierarchical clustering is not easy to determine the dividing layers. As a result, DBSCAN is used in this paper.

DBSCAN has two advantages: (1) High efficiency. It can remove multiple outliers at one time. Compared with other clustering methods, iterations are fewer. (2) The clustering number is reduced with the increase of iterations, and finally into one class.

Before using DBSCAN, it is necessary to set two parameters, one is the cluster radius ε, the other is the minimum number of points (MinPts), sometimes named as cluster density. ε indicates the distance between data in one class. MinPts specifies the minimum number of data contained in one class. How to set the two parameters is described in the next experiment.

Ricker wavelet is the typical waveform in vibration simulation. The time-domain expression is shown as follows (25) s(t)=(1−2π2fM2t2)e−π2fM2t2 s\left( t \right) = \left( {1 - 2{\pi ^2}f_M^2{t^2}} \right){e^{ - {\pi ^2}f_M^2{t^2}}} where t is time, and f_M represents the peak frequency. Fourier transform is applied to Equation (25) to obtain the expression of Ricker wavelet in the frequency domain (26) S(f)=2f2πfM3e−f2fM3 S\left( f \right) = {{2{f^2}} \over {\sqrt \pi f_M^3}}{e^{ - {{{f^2}} \over {f_M^3}}}} The oscillogram of Ricker wavelet with the center frequency of 50 Hz is shown as follows

Fig. 3

Figure 3 shows that the Ricker wavelet is not only similar to real vibration signals in the time domain but also in the frequency domain. Therefore, the Ricker wavelet is the desirable waveform of simulation vibration signals.

The vibration source data was composed of five Ricker wavelets with a frequency peak of 50 Hz. To improve the data quality, Ricker wavelet amplitude was expanded five times. After attenuation, the vibration source data was mixed with Gaussian noise with zero mean and 0.05 variance. The sampling frequency was set as 10000 Hz, the duration was 3 s, the wave velocity was 10000 m/s, and the minimum distance between the optical fibre and vibration source was 1 m. There were 19 fibre segments receiving vibration signals, and their serial numbers based on MMDFS were (−9, −8. . . −1, 0, 1. . . 8, 9).

At the same time, all data was stored as a matrix, with 30,000 rows and 19 columns. The stored data is shown in Figure 4.

Fig. 4

In Figure 4(b), each column refers to data received by one fibre segment. Among these columns, the 10th column is the datum fibre segment and has the best signal quality with an SNR of −1.05 db. But the first and 19th columns have the worst data quality with an SNR of −2.12 db. Five Ricker wavelets in signals of No. 10 column can be observed obviously, but the Ricker wavelets are extremely close to noises and challenging to distinguish in signals of No. 1 column. The simulation fibre data used in this paper is named as SimuFiberData.

The experiment used the method of GCCSC based on ML to estimate TDOA. To improve accuracy, three-point interpolation was conducted. Table 1 shows the simulation fiber data for making a comparison on the real TDOA with five methods, including CC, SCOT, ML, GCCSC based on SCOT and GCCSC based on ML. It is worth noting that the unit of TDOA is the sampling number, instead of the second. Mean Square Error and success rate of five methods in Table 1 are shown in Table 2.

TDOA estimation in Table 1 has two unreasonable situations: (1) TDOA value is negative. For example, the estimation value of No. −9 in SCOT is −2; (2) TDOA value dramatically exceeds the theoretical value. For example, the estimation value of No. 1 in SCOT is 13, but the actual value is 1.236. These belong to unsuccessful estimation values. Thus the success rate is also an important indicator to estimate TDOA.

The success rate of the five methods in Table 2 shows that GCCSC based on ML has the highest success rate than other methods. The mean square error estimated by GCCSC based on ML is the minimum and it has the highest success rate. For this reason, GCCSC based on ML was selected to do TDOA estimation.

To obtain satisfactory distance estimation results, not only the accurate TDOA method should be selected, but also the TDOA outliers should be removed reasonably. The entire detailed computational process of Simu-FiberData is shown in Table 3.

The 1st column in Table 3 is the serial numbering for TDOA. Successively, the 2nd column is the true TDOA; the 3rd column is the TDOA estimation value according to GCCSC based on ML. After the 3rd column, data in two columns are divided into a group. Each group represents one distance estimation. No. 1 column of each group is the serial numbering based on MMDFS. If the serial numbering is 0, it means data at this point have been eliminated from the distance estimation. No. 2 column in each group is the corresponding DS of data involved in TDOA estimation. The antepenultimate row is the mean square error of all DS values after each distance estimation, implying the effect of DS clustering. The penultimate row is the estimation value of wave velocity. The tailender row is the estimation value of the distance between the vibration source and optical fibre.

In the first estimation, wave velocity is 8508.390, the distance is 1.419, both of which belong to the reasonable value range. However, DS has the more significant mean square error and it is 105.451, indicating that the clustering result is not desirable. The TDOA value involved in estimation probably has outliers, which should be eliminated.

The primary step of eliminating outliers was to set up cluster radius ε and cluster density MinPts, both of which required pre-estimation. Moreover, the scatter diagram of DS can provide intuitional help for pre-estimation. As shown in Figure 5, the hollow points are normal points, but the solid points are suspected outliers. The distribution of the hollow points is relatively concentrated. The distance between two hollow points is about 10. The number of hollow points is about 10. The distribution of solid points is relatively scattered, and the distance between two solid points is >50. There were six solid points, among which five are within the range of 150–300, and one point is far away from others, therefor the five points were more likely to gather into a cluster. Cluster radius ε could be set slightly larger than the distance between two normal points. Cluster density MinPts was determined by the number of outliers in a cluster. Cluster radius ε was set as 10 and cluster density MinPts was set as 5 in the first time removing outlier.

Fig. 5

After removing outliers in the first time, the second estimation was conducted to obtain wave velocity as 9370.197 and distance as 0.958, close to the set value 1; DS mean square error was very small, that is 2.2577, indicating that the clustering result was desirable and outliers had been eliminated. Figure 5(b) shows the scatter diagram of all DS values in the No. 7 column in Table 3—the scatter diagram after removing outliers, showing that the DS values of residual points were maintained at 10 and there were no outliers. The whole calculation was over.

Fibre data were collected from the space in a college in Beijing, China. The space was an area of 800 m2 in a rectangle. The south-north length was about 20 m and the south-east length was about 40 m. The length of fibre was 300 m, paved around the site with two circles. Its buried depth was about 10 cm under the ground.

The fibre was wireline cable. NBX-S3000 from the Japanese Neubrex Company was used as the fibre collecting device. The sampling frequency was 10000 Hz. Fibre spatial resolution was 1.03 m, thus the entire fibre could be divided into 290 segments.

The excavation was carried out at 108.83 m of the fibre. A shovel was used to excavate on the ground. The excavation position was distant from the fibre about 1 m. The strongest excavation signal was at 108.83 m. Excavation signals could be distinguished at 99.59 m and 117.04 m. But it was difficult to distinguish farther excavation signals with the naked eyes. The vibration signals received from 99.59 m to 117.04 m were chosen as the analysis objects. According to MMDFS, the fibre segment at 108.83 m was the datum fibre segment with the serial numbering as 0. The serial numbering at 109.86 m, 117.04 m, 107.8 m and 99.59 m was 1, 8, −1 and −9, respectively. There were a total of 18 fibre segments, which had the sequence numbered as (−9, −8, −7, −6, −5, −4, −3, −2, −1, 0, 1, 2, 3, 4, 5, 6, 7, 8).

Practical fibre data were greatly affected by the background noises, and the SNR was low. The best SNR obtained from the strongest excavation signal at 108.83 m was 17.72 dB, the worst SNR was 9.49 dB. At 117.04 m, the best SNR was 8.45 dB; the worst SNR was −2.27 dB.

Relative to the simulation fibre data, practical fibre data were affected by noises and soil medium, resulting in greater deviation. As estimating TDOA, more unreasonable values than simulation fibre data were generated, sometimes accounting for >30% of the total data. Therefore, the main task of pre-treatment was to remove unreasonable TDOA.

Fig. 6

As shown in Figure 6, the fibre data contains four groups of apparent excavation signals circled with dotted lines. In distance estimation, four segments of data are partitioned. The length of each segment ranges from 130 to 200 sampling points. Excavation data in each group regards the No. 0 column as the benchmark to estimate TDOA with other columns and eliminate unreasonable TDOA, remaining 25 effective values.

Because the method adopted for estimating TDOA has strong noise resistance, this paper does not carry out de-noising processing on optical fibre data.

The complete details of distance estimation for excavation data are recorded in Table 4.

Differing from Table 3, because true TDOA is unknown, there is no real value column in Table 4. The rest is the same as Table 4. The original number ranges from 1 to 72. After pre-treatment, only 25 valid TDOA values remained. No. 2 column is the TDOA estimation value. The entire calculation conducted distance estimation five times. Outliers were eliminated four times. The setting of cluster radius ε and cluster density MinPts each time were listed in the last two rows of Table 4. Finally, the distance was 0.951 m, close to 1 m. The experiment proved that the above-mentioned distance estimation method was effective.

However, this method has the possibility of failure. Especially when there are many outliers with the DS values close to each other, they are easy to aggregate into a class, resulting in the deletion of the normal points. In this case, the result has a significant deviation.