Features based on instantaneous frequency for seismic signals clustering

2.1. Segmentation procedure

Analyzed seismic recordings may consist of a few separate events. Thus, to analyze single excitation, segmentation procedure is recommended.

The segmentation procedure used in this paper is based on time-frequency analysis and is focused on detection of significant changes of the amplitude spectrum [19]. Moreover, this procedure trims the signals, thus signal portions with near-zero amplitude are cut-off. The main assumption of the segmentation was that each segment consists of one seismic event. However, some inconveniences may occur, when two separate events overlap significantly. It might be difficult to recognize event consisting of e.g. blasting and roof collapse, or blasting and rockburst.

2.2. Empirical mode decomposition

An intrinsic mode function (IMF) is defined as a function which satisfies following conditions: number of zero crossings and local extrema must either be equal or differ at most by one, and at any point mean value of the envelope defined by the local maxima and that defined by the local minima is zero [9].

It is known, that signal’s instantaneous frequency is meaningful if and only if the signal is symmetric with respect to the local zero mean and has the same number of zero crossings and extrema [9, 20]. It has to be monocomponent signal [15].

It was shown that IF obtained from a function satisfying IMF conditions almost always is correct and in worst condition it is still consistent with the physics of the system studied [9].

Empirical mode decomposition (EMD) is an algorithm that decomposes basic signal into several functions that satisfies IMF conditions, and a Residuum which is lower (in terms of amplitude) than the predetermined value, or it is a monotonic function. EMD is an iterative, data driven procedure. We make use of stopping criterion proposed in [21].

2.3. Calculation of an instantaneous frequency

Instantaneous phase for a real-valued function $f\left(t\right)$ is the argument of its analytic representation. The analytic representation $f_a$ may be received by $f_a\left(t\right)=f\left(t\right)+j\widehat{f}\left(t\right)$, where $j$ is imaginary unit, and $\widehat{f}$ is function’s Hilbert transform, which is described as a convolution of the signal with a $h\left(t\right)=1/\pi t$ function.

Instantaneous frequency is a time derivative of unwrapped instantaneous phase.

2.4. Signals’ analyzed features

It has been noticed that in case of seismic signals, when dealing with such low noise ratio, 1st IMF well resembles the original signal. Thus it was decided to represent segments’ by the IF of their first IMFs and to base features that well distinguish signals on it.

2.5. Clustering methods

Features’ results need to be subjected to clustering methods in order to show their potential of seismic signal discrimination. In this paper we compare results of two clustering methods: $k$-means, and agglomerative hierarchical algorithm.

2.5.2. Agglomerative hierarchical clustering

The analyzed features have been also tested using a deterministic method, namely agglomerative hierarchical clustering. AHCA begins with locating each observation into single-element cluster, and merges 2 nearest (in terms of specific measure) groups as it proceeds. In this paper the distance measure proposed by Ward is used [25]:

2

$d_w\left(r,s\right)=\frac{n_r{n}_s{\left(d\left({\stackrel{-}{x}}_r,{\stackrel{-}{x}}_s\right)\right)}^2}{n_r+n_s},$

where ${\stackrel{-}{x}}_r=\left({\sum }_{i=1}^{n_r}{x}_{ri}\right)/n_r$, and $x_{ri}$, $i=$ 1,…, $n_r$ are all elements which belong to cluster $r$, and $d\left(\bullet ,\bullet \right)$ is Euclidean metric.

Fig. 1Diagram of signal analysis’ procedure

2.6. Silhouettes method

Clustering procedures require an a priori assumed number of clusters. So there is a need to estimate the optimal number of clusters in analyzed data. Visual analysis of clustering procedure results can be supported/dismissed by e.g. silhouettes [26]. Their construction is based on the concept of average dissimilarities between an object and other elements in clusters. One can construct a silhouette statistic, which for a fixed number of clusters $k$ assigns to $i$th observation following value:

where $a\left(x_i\right)$ is the average distance to cotenants in cluster where $i$th case is allocated and $b\left(x_i\right)$ is a distance to nearest neighboring cluster’s center. The final output of the silhouettes method is $SC=\left({\sum }_{i=1}^{n}s\left(x_i\right)\right)/n$ (silhouette coefficient) and this value is therefore interpreted according to Table 1. The value of the criterion varies from 0 to 1.

Average value of silhouettes over all data (Silhouette Coefficient, $SC$) is a measure of how properly the clustering is performed. Table 1 represents interpretation of $SC$ proposed by Kaufman and Rousseeuw [18].

The entire procedure of signal preparation before clustering is presented in (Fig. 1).

3.1. Data preparation

In this paper we analyze 98 signals (segments) which are obtained from 51 seismic recordings.

In Fig. 4 one can find an exemplary result of the segmentation procedure described in [19]. Trimming property is easy to notice. What also could be noted is that near $t=$3 s there is a segment that consists of couple single events. When single excitations are overlapping, segmentation procedure may have difficulties in division.

Fig. 3Typical recording of a seismic signal (natural event)

After the segmentation single-event signals are subjected to EMD algorithm and, as it was mentioned, for each segment we compute instantaneous frequency of a first IMF.

Fig. 4Typical recording of a seismic signal (mining face blasting) with applied segmentation

In Fig. 5. we can see an exemplary segment, first four IMFs and corresponding IFs. It can be noticed that IF of raw signal is meaningless, since negative frequencies can be noticed [15]. Also some negative values might be seen on IF of 4th IMF. It is due to possible imperfection of IF from IMF. There exist methods that can deal with that problem [27, 28]. However, we assume that problem of negative IMFs’ IF is in this case marginal and we haven’t considered it further yet.

Afterwards, 2 features are extracted from IFs, that were proposed by the Authors in [5].

First of the discrimination features is a vector of probability density function values of the 1st IMF’s IF, derived from Rozenblatt-Parzen kernel density estimator with Gaussian kernel [22, 23], in fixed number of points (i.e. 512), located uniformly between 0 and 625 Hz (as the sampling frequency of gathered data is 1250 Hz). In Fig. 6, we present the segments’ first feature representations.

The second feature, which is specific condensation of the previous one, represents each seismic signal segment by its dominant frequency of 1st IMF’s IF. With such reduction we lose significant part of information, however it is simpler to decide whether two clusters shall be cotenants or not. In Fig. 7 we include values of 2nd feature presented on a histogram.

Fig. 5Typical seismic signal with 5 first IMFs and corresponding instantaneous frequencies

a) Signal

b) IF of signal

c) 1 IMF

d) IF of 1 IMF

e) 2 IMF

f) IF of 2 IMF

g) 3 IMF

h) IF of 3 IMF

i) 4 IMF

j) IF of 4 IMF

Fig. 6Probability density function estimators for all segments

3.2. Silhouette results

Each feature was analyzed by silhouettes in order to find optimal number of classes. The results are presented in a table below.

As it can be seen from Table 2, best results of $SC$ for distribution of IF, are acquired with two-cluster division (for both clustering methods). $SC$ for $k=\text{2}$ dominate over other options. All results may be interpreted as “reasonable structures” (Table 1). The differences between clustering methods for $k\in \left\{\text{2,}\text{3,}\text{4,}\text{5,}\text{6}\right\}$ are not significant, i.e. each value of $SC$ is between 0.5 and 0.7.

Fig. 7Histogram of 2nd feature values

Fig. 8Results of IF density estimator – 98 segments grouped into 2 clusters by both algorithms

a)$K$-means cluster 1/2

c) AHCA cluster 1/2

b)$K$-means cluster 2/2

d) AHCA cluster 2/2

Illustration of clustering into 2 groups obtained using k-means and AHCA is presented in Fig. 8. Results of AHCA and k-means are mostly coinciding. First cluster (for both methods) possess dominant shape of densities: medium-dispersal with mode located around 300 Hz. Additionally, couple of visible outliers are located in this cluster, however their number is negligible ($k$-means: 6, AHCA: 4). Cluster 1 (especially after AHCA) seems to be well build. However, cotenants in the 2nd cluster seems to be too much diversified in terms of dispersion and dominant frequency. It is possible that clustering with another number of groups would be much better in terms of visual inspection.

Table 3 shows that the second feature has $k=\text{2}$ as an optimal number of clusters. Results for all sets from $k=\text{2}$ to $k=\text{7}$ are evaluated as “strong structures” (Table 1). Number of clusters that was evaluated as the least is $k=\text{4}$ ($k$-means) and $k=\text{5}$ (AHCA). On the other hand, in Fig. 9 one can notice several modes in distribution of the feature. Similar to method presented in [29], we divide the frequency axis in order to separate each significant mode.

In Fig. 9 results of k-means and AHCA clustering are contained. Segments which are on the left side of the line have been allocated in the first cluster, and these on the right side – in the second one. The line does not indicate any specific frequency. Both of the methods separated the segments into 2 groups in the same way. Although, 4 clear modes can be noticed in the histogram. Moreover, 6 segments with mode located between 230 and 280 Hz could even build another group.

This leads us to believe that visual-best results may not be consistent with that obtained from the automatic algorithm.

Fig. 9Location of the mode – 98 segments partitioned into 2 clusters by k-means and AHCA

3.3. Best visual results

Recall that the silhouette method indicated that 2 is the optimal number of clusters for both instantaneous frequency features. However, visual inspection of these results (Fig. 8, Fig. 9) pointed out that the obtained clusters consist of segments that differ significantly in terms of dominant mode and density of instantaneous frequency. Thus, an attempt to find the number of clusters which would better fulfill the main assumptions of good clustering (similar cotenants and different neighboring clusters) has been made.

By empirical tests we decided, that 5 is the best number of classes for density of IF. Lower quantity makes cotenants too different, and the larger make neighboring clusters too similar.

Clustering results of AHCA (Fig. 11) seems to be consistent with those derived by the k-means (Fig. 10). Such phenomenon indicates high robustness of a discrimination feature. As it can be seen, the algorithm distributed the segments with respect to theirs IF distribution and dominant frequency value. Instantaneous frequencies concentrate in one of 4 different modes 140 Hz, 170 Hz, 200 Hz (2 groups) and 300 Hz. They manifest different dispersions within the clusters. First set (140 Hz, 5/5 in Fig. 10, 1/5 in Fig. 11) has the highest dispersion, also in couple of its segments multimodal feature can be seen. Set which has the 2nd lowest dispersion (200 Hz, 3/5 in Fig. 10, 2/5 in Fig. 11), shares modal frequency with the set whose dispersion is largest (200 Hz, 4/5 in Fig. 10, 4/5 in Fig. 11). 3rd (170 Hz, 2/5 in Fig. 10, 3/5 in Fig. 11) and 5th (300 Hz, 1/5 in Fig. 10, 5/5 in Fig. 11) clusters have similar dispersions, but different modes.

Fig. 10Results of clustering using IF density estimator – 98 segments partitioned into 5 clusters, k-means

a) Probability of frequency (cluster 1/5)

b) Probability of frequency (cluster 2/5)

c) Probability of frequency (cluster 3/5)

d) Probability of frequency (cluster 4/5)

e) Probability of frequency (cluster 5/5)

Fig. 11Results of clustering using IF density estimator – 98 segments partitioned into 5 clusters, AHCA

a) Probability of frequency (cluster 1/5)

b) Probability of frequency (cluster 2/5)

c) Probability of frequency (cluster 3/5)

d) Probability of frequency (cluster 4/5)

e) Probability of frequency (cluster 5/5)

One can notice couple outliers (Fig. 10 cluster 1/5 and 3/5, Fig. 11 cluster 2/5 and 5/5). They correspond to 6 segments with mode located between 230 and 280 Hz in Fig. 9. Thus considering them as outliers in terms of this feature may guide to repeat consideration for the second discrimination feature.

The AHCA four-cluster results seems to be reliable (Fig. 12). The clusters are similar to those obtained from first feature. The clusters are concentrated in 140 Hz, 170 Hz, 200 Hz, and 300 Hz (outliers are allocated in this set). Two groups of segments with mode located around 200 Hz from previous discrimination feature are now merged into one set.

Fig. 12Results of clustering using mode location – 98 segments partitioned into 4 clusters, AHCA

However, when $k$-means is applied couple of times with 4 as a number of clusters it can be noticed, that in many cases outliers form a separate set (see Fig. 13). Thus, there must be such a cluster that contains segments with significantly different modes (here, clusters with modes at 140 Hz and 170 Hz), which seems to be an inappropriate result (due to analysis of Fig. 9). Thus, we consider 5 clusters to improve clustering results.

The results of AHCA clustering into 5 groups divide the cluster with modes centered around 300 Hz into 2 clusters: segments with dominant frequency around 300 Hz and 6 outliers mentioned previously.

Fig. 13Results of clustering using mode location – 98 segments partitioned into 4 clusters, k-means

In Fig. 14 we can see $k$-means clustering. The segments with modes at 170 Hz are slightly mixing with the 200 Hz segments. The algorithm tends to return a bit random assignment especially for outliers and segments represented by modes between 200 and 170 Hz. However, the general results are satisfying, since $k$-means does not merge all of the segments with modes between 170 and 200 Hz into one cluster.

The fact that 2 different clustering methods return us similar (or even the same) results proves reliability and robustness of the considered features.

Fig. 14Second feature under the k-means, 5 clusters