Stochastic Modelling as a Tool for Seismic Signals Segmentation

Kucharczyk, Daniel; Wyłomańska, Agnieszka; Obuchowski, Jakub; Zimroz, Radosław; Madziarz, Maciej

doi:https://doi.org/10.1155/2016/8453426

Shock and Vibration

On this page

Abstract Introduction Conclusions References Copyright Related Articles

Research Article | Open Access

Volume 2016 | Article ID 8453426 | https://doi.org/10.1155/2016/8453426

Stochastic Modelling as a Tool for Seismic Signals Segmentation

Daniel Kucharczyk,¹Agnieszka Wyłomańska,¹Jakub Obuchowski,²Radosław Zimroz,²and Maciej Madziarz²

Academic Editor: Longjun Dong

Received15 Dec 2015

Revised28 Feb 2016

Accepted23 Mar 2016

Published08 May 2016

Abstract

In order to model nonstationary real-world processes one can find appropriate theoretical model with properties following the analyzed data. However in this case many trajectories of the analyzed process are required. Alternatively, one can extract parts of the signal that have homogenous structure via segmentation. The proper segmentation can lead to extraction of important features of analyzed phenomena that cannot be described without the segmentation. There is no one universal method that can be applied for all of the phenomena; thus novel methods should be invented for specific cases. They might address specific character of the signal in different domains (time, frequency, time-frequency, etc.). In this paper we propose two novel segmentation methods that take under consideration the stochastic properties of the analyzed signals in time domain. Our research is motivated by the analysis of vibration signals acquired in an underground mine. In such signals we observe seismic events which appear after the mining activity, like blasting, provoked relaxation of rock, and some unexpected events, like natural rock burst. The proposed segmentation procedures allow for extraction of such parts of the analyzed signals which are related to mentioned events.

1. Introduction

Real world signals are modeled in order to extract information about the analyzed phenomena. However, very often the observations exhibit properties adequate to nonstationary processes. The nonstationarity can be manifested in different ways. One of the possibilities is the coexistence of two or more processes and in fact the analyzed time series is a mixture of different stochastic systems. There are many examples of such phenomena. For example, we can observe that for a given time period the analyzed time series has properties related to one process whereas for other time intervals it completely changes the character, so it exhibits properties related to different process. In this case it switches between two (sometimes more) processes which can not be described by the same model or in less radical case the parameters of the model change. Therefore, the first step of experimental data analysis is the segmentation, that is, extraction of such parts of the signal which have homogeneous properties.

There are many segmentation methods in the literature. Most of them were invented and dedicated for specific applications. Signal segmentation might be done in different domains. Most of the methods are based on the properties of the signal in time domain. We should mention here the statistical techniques applied to analyzed time series [1, 2] and methods which take into consideration appropriate models describing the signal [3–5]. There are also methods which are based on the properties of the signal in other domains [6–11]. In the literature one can also find segmentation techniques based on the time series decomposition [12–14].

The segmentation problem can be found in several applications. We only mention here a few of them like condition monitoring [9, 14], machine performance analysis [15], experimental physics [1, 16, 17], biomedical signals (like ECG signals) [7, 18–23], speech [24–26], and econometrics [27, 28].

The problem of extraction of parts of the time series with homogeneous properties is also important in analysis of seismic signals [8, 29–34]. One of the most important issues is to identify, quantify, and understand seismic activity by seismic signal analysis. Mining operations in underground mine produce many different seismic events related to mining room and pillar technology, blasting technology, complicated geomechanical conditions in deep (>1000 m) underground mine, heavy duty machinery systems, and so forth [11].

During regular mining activities in the underground mine, seismic events might be related to events induced in somehow predictable way by mining activity (blasting procedures, provoked relaxation of rock) and unexpected, as natural rock burst. Apart from these two main sources of events, other activities in the mine might also initiate seismic signal recording procedure (e.g., mobile machine moving nearby the sensor). Obviously, significance of each type of recorded signal is very different; its shape in time domain, energy, and frequency structure (spectrum of the signal) will be different. Moreover, depending on the physical background of the event, seismic signal might contain very complex structure that in fact is a sequence of impulses related to sequence of events [11].

In this paper we propose two novel methods of seismic signal segmentation. Both of them are based on the statistical properties of analyzed time series and take into consideration models appropriate for different parts of the analyzed signals. In the first method we assume that the noisy part of the signal (not related to the seismic event) is Gaussian distributed whereas the nonnoisy part (that can be related to seismic event) is described by other distributions. By using this assumption we can extract those parts of the signal that can be described by Gaussian distribution and finally the procedure leads to extraction of those segments which may be considered as seismic events. The second proposed approach also takes into consideration the statistical properties of the underlying time series. However, here we assume that the noisy as well as the nonnoisy parts are described by autoregressive model of order 1. The only difference between those segments is in the parameters of the autoregressive model. In this method we propose to use the Markov regime-switching model in order to obtain the probabilities of being in one of the mentioned regimes. On the basis of the estimated probabilities we can easily segment the signal in order to detect the parts with homogeneous properties. It is worth mentioning that the Markov regime-switching models (called also hidden Markov models) were initially used as segmentation methods in speech recognition [35], but they have been also exploited in other contexts, including as diverse areas as software reliability [36], electromagnetic field measurements [37], or electricity markets [27]. Other interesting applications of Markov regime-switching model can be found in [38–42]. This model was also used in diagnostics [43–47]. We compare the proposed segmentation methods to the classical technique of seismic signal analysis called STA/LTA ratio method that take into consideration also the properties of the time series in time domain [48]. As it was mentioned, the segmentation of seismic signal is a first step of experimental data analysis. It is especially important in the problem of locating a single event [49–51].

The rest of the paper is organized as follows. In Section 2 we present in detail the proposed segmentation methods and we remind the readers with the short-term-average and long-term-average ratio technique and the classical method of seismic signal analysis. Next, in Section 3 we check their effectiveness for simulated time series imitating seismic signals. In order to check versatility of invented techniques we analyze two different signals: with single seismic event and with multiple events. We compare the results obtained for proposed segmentation techniques to those obtained for STA/LTA ratio method. In Section 4 we analyze the real seismic signals acquired in a copper-ore underground mine. Last section contains conclusions.

2. Methodology

In this section we introduce two novel methods of seismic signal segmentation. Both of them are based on the stochastic properties of analyzed signal in time domain. In both proposed methods we assume that the noisy part of the signal can be modeled by using Gaussian-based system. However, in the first approach we assume the noisy part is a sequence of independent identically distributed random variables whereas, in the second one, the process is supposed to be driven by autoregressive model of order 1. Moreover, in the first proposed method we do not assume the distribution of the segment related to nonnoisy part (which can be related to seismic event) whereas in the second approach we take into consideration the nonnoisy part driven by another independent autoregressive process of order 1. The rationale behind the assumption of distribution of the noisy part of the signal is driven by the fact that the recorded real signal is characterized by a discrete time series (the discreteness of the time series is related to the sensitivity of the sensors) for which the proposed Gaussian distribution is of the central limit theorem which says that any distribution with finite second moment can be approximated by Gaussian one. For many real data analyses this approach is taken as well. According to the assumption of AR(1) model: in the autocorrelation function (ACF) and partial autocorrelation function (PACF) of the noise for real signal we observe behavior adequate to autoregressive models; that is, the ACF has sinusoidal behavior and PACF is larger than zero only for few lags. The simplest autoregressive model is AR model of order 1. If the order is larger then the calculations are more complicated and actually if we take the larger model the results will be very similar to those we obtained. In the last subsection we recall the classical method of seismic signal segmentation based on the STA/LTA ratio technique.

2.1. Cumulative Sum of Gaussian Probability Density Functions

In this part, we introduce a new seismic signal segmentation method based on the use of a suitable test statistic , defined as a cumulative sum of zero-mean Gaussian probability density functions; that is, where is the value of Gaussian probability density function given at point , with a zero mean () and relatively small variance () given as follows: where denotes the value of the observed signal at time . In other words, we assume that the noisy part of the seismic record is considered to be Gaussian distributed according to the probability density function given in (2), whereas the nonnoisy part is assumed to be driven by some unknown distribution. Note though that the assumption about the nonnoise is not relevant in the presented approach, as there is a one-to-one correspondence between the noise segments of the signal and the segments affected by the seismic event; that is, identifying the noisy segments of the record at hand gives us a direct indication about the location of the seismic events belonging to the nonnoisy segments.

From the above definition of the test statistics given in (1), it is clear that, for each point belonging to the noise part of the signal, the value of the Gaussian probability density function given at that point is strictly positive assuming that the true noisy part is generated according to (2). In case when the observation comes from the segment governed by the seismic event the value of (2) is most likely to be neither equal nor close to zero. By taking these into account, it becomes apparent that the behavior of in signal (nonnoisy) segments is highly different from its behavior in noise intervals. To put this into a more analytical perspective, the slope of the curve whose values are calculated in accordance with (1) is neither positive nor zero depending on whether the time interval is considered to be noise or the proper seismic signal segment. The aforementioned behavior of the used statistics in a simulated signal record is presented in Figure 1. To this end, we have not described the way of choosing the proper value of the variance parameter defined in (2), as it depends on the sensitivity of the seismometer used for recording seismic signals. The parameter can be characterized as a seismometer driven parameter that can be easily derived by taking the sample standard deviation from the recorded noise segments.

(a)

(b)

(c)

After plotting the values of test statistics for each point , it is clearly visible that the period of times characterized by the constant value of is affected by seismic events; see Figure 1. In order to extract the periods of times (segments) where the seismic events have taken place, we propose to use the vector of differences between consecutive observations taken from the original input seismic record and mark the periods of time characterized by relatively small increments as the period of time driven by a seismic event.

Having obtained the vector of differences, we check whether the values are equal to 0 or not in order to identify the period of time where the seismic events occur. To remove sudden drops and peaks in the output vector, we follow the procedure given below.(1)When the derivative is negative then we check the signal until the derivative is positive. If the space between them is less than steps then we fill the output vector with ones at this location.(2)When the derivative is positive then we check the signal until the derivative changes sign. If the size of this gap is less than steps then we fill the output vector with zeros at this location. The last step is to remove segments shorter than seconds (under the assumption that the sampling interval of the record is equal to ms, equivalently all segments containing less than 200 observations will be removed), since they are too short to analyze and can be some artefacts rather than true events. The final result of the procedure can be seen in Figure 1, whereas the intermediate results of applying the smoothing procedure to the initial vector of conditional probabilities are shown in Figure 1(c).

2.2. Markov Regime-Switching Models

In this section, we consider Markov regime-switching (MRS) model with application to the seismic signal segmentation problem. The general idea behind MRS model is to represent the observed stochastic behavior of a specific time series described by at least two separate states (regimes) each driven by different underlying stochastic process. For the purpose of discriminating the seismic signal between noise and proper signal segments, it is enough to assume that the underlying time series (the recorded signal at hand) can be fully described by two separated regimes; that is, the first regime is responsible for characterizing the period of time where there is no seismic event and the second one represents the behavior of the seismic signal governed by the seismic event.

The switching mechanism between the states is assumed to be an unknown unobserved Markov chain . In general, the Markov chain is described by the transition matrix containing the probabilitiesof switching from regime at time to regime at time . Assuming that there are only two regimes, , we have In theory, the MRS models can be divided into two groups depending on whether the regimes under consideration are independent of each other or not. In the first case only the model parameters change depending on the state process values, while in the second case the individual regimes are driven by independent process. Following the notation used by [52], the first case is commonly described by a parameter switching time series of the form given as follows: where are assumed to be Gaussian distributed , that is, with parameters and equal zero and unity, respectively. In the second case with independent regimes, the underlying time series is defined as where at least one regime is defined as AR(1) model given in the form

For the purpose of presenting the usage of Markov regime-switching models in the area of seismic data segmentation, we limit ourselves to the case where the model with independent regimes is taken into consideration. Moreover, we would assume that both regimes are driven by autoregressive models of order 1, in brief AR(1).

As one might expect, the calibration procedure of MRS models is the most difficult part as the regimes are only latent and hence not directly observable. In this paper we follow the iterative calibration procedure introduced by [52], where the vector of initial parameters, for , where , denotes parameters of th regime distribution and is stationary probability in Markov chain. This procedure is divided into two parts consisting of estimation of latent observation of Markov process and estimation of parameters of given processes.

The first step of this iterative procedure called E-part estimates conditional probabilities of being in regime at time , since switching mechanism is not directly observable and we can only derive expected value of the process which is mentioned conditional probability , where denotes the realizations of the process . The results from this step are used in the M-step, where the vector of parameters is calculated using the modified maximum likelihood method. After both steps likelihood function is calculated and when it reaches the local maximum the procedure ends.

The detailed description of the procedure is as follows.

E-Step. To calculate probabilities in this step we need a vector from previous iteration ( denote the th iteration of the procedure) and a vector of realization of the process up to given point time . This part is divided into two steps.

(i) Filtering. At this step we derive the probabilities with given time in the iterative equation. Let ; then where is the conditional density of the regime process at time . The starting point of this iteration is given as .

(ii) Smoothing. In this part we are calculating conditional probabilities with given ; this means that we are conditioning with whole history of the process . This step goes backward, , and iterates on

To calculate probabilities in (9) we need a derivation of . For the first model which is process with switching parameters (dependent regimes) and the same Gaussian random noise we use conditional distribution formula given as For the model described in (6), with two independent random variables we use simple formula where is a density function of a chosen random variable in regime , since this random variable is independent and identically distributed.

The calculation becomes more complicated when at least one of the regimes is given by process and is independent of the other ones. As the regime process is latent and unobserved, the way of conditioning with given makes it dependent on the whole history of the process and in this case all possible paths of this process up to time should be considered. This gives us possible ways and it increases exponentially with the size of the sample. When probability of 10 consecutive observations of one regime is negligible then we can use just 10 last observations. This method was proposed by [53], but this still makes computation very intensive and we cannot use this assumption for our data. In order to avoid the computational burden of the algorithm, we would modify it following the logic described in [52], where authors suggest to replace the value of in formula (10) with their expectation , which contains all information about process until time . This value can be computed with the following formula: Using formula (12) one can rapidly reduce computational complexity of this algorithm without losing quality of the estimation.

M-Step. The -step estimates both regime distribution parameters using the modified maximum likelihood method, the same as in standard algorithm of maximum likelihood method of estimation where we use log-likelihood function , but each component of this sum is weighted with corresponding probability (9). In particular, for the model defined by (6) explicit formulas for the estimates can be obtained by setting the partial derivatives of the log-likelihood function to zero. For MRS process given by (5), estimation of th regime parameters leads to the following formulas:

For the model defined in (6) where at least one of the regimes is independent of process we can use the estimators given above, but we need to replace observation with its expectation given as in (12).

For independent random variables it is enough to use the maximum likelihood estimators of their parameters with the same modification as in the case of calculating the estimators for model (5). In case of the autoregressive model of order 1 with Gaussian noise, the estimators are given by

For the purpose of discriminating the seismic signal record between the noisy segments and the seismic events, we use the vector of probabilities given as in (9). In addition to that we introduce the so-called smoothing procedure responsible for creating proper segments; that is, depending on the characteristics of the seismic record one might expect that the observed vector of probabilities will not behave in a smooth way, and there might be some periods of time where the probabilities of being in the appropriate regime vary over the time so frequently causing gaps between the seismic events. The procedure itself is based on the vector consisting of elements calculated as a difference between consecutive observations. Having obtained such a vector, we check the nonzero values to identify the period of time where the seismic events occur. To remove sudden drops and peaks in the output, we follow the procedure given below:(1)When the derivative is negative then we check the signal until the derivative is positive. If the space between them is less than steps then we fill the output vector with ones at this location.(2)When the derivative is positive then we check the signal until the derivative changes sign. If the size of this gap is less than steps then we fill the output vector with zeros at this location. The last step is to remove segments shorter than seconds, since they are too short to analyze and can be some artefacts rather than true events. The final result of the procedure applied to simulated signal can be seen in Figure 2.

(a)

(b)

2.3. Short-Term-Average and Long-Term-Average (STA/LTA) Ratio

In this subsection we briefly recall the idea of the method proposed in [48], called short-term-average and long-term-average (STA/LTA) trigger method. The underlying idea of the STA/LTA method is to evaluate in a continuous fashion the value of characteristic function (CF) of a seismic window in two consecutive moving-time windows in order to detect the seismic event. The characteristic function used for calculation purpose can be defined as energy, absolute amplitude, or envelope function of the microseismic trace. Irrespective of the definition of the characteristic function (CF), the short time window (STA) is supposed to measure the instantaneous amplitude of the seismic signal, whereas the long time window (LTA) provides information about the temporal amplitude of seismic noise. When the ratio of both, the short-term-average and the long-term-average, exceeds a predefined value (activation threshold), the consecutive recorded values are marked to be event-driven until the ratio falls below another predefined value , called the deactivation threshold. In general the STA/LTA ratio algorithm can be described as follows: where and denote the number of observations to be available in the short and long time windows, respectively. In Figure 3 we present the results of applying the STA/LTA ratio method to simulated seismic signal with the characteristic function (CF) defined in terms of signal energy. In addition to that we have assumed that the short and long time windows are equal to 0.3 s and 1.5 s, respectively; the activation and deactivation thresholds correspond to and , respectively.

(a)

(b)

3. Simulations

In this section we describe the results of performed simulation study of the methods presented in Section 2. As a reminder, we present the results for three different segmentation methods, namely,(i)the novel approach where the cumulative sum of Gaussian probability density function is used to identify the noisy segments of the seismic record: this method assumes that the observations from the noisy segments of the data are independently and identically distributed according to the Gaussian distribution with parameters and ;(ii)the method where the Markov regime-switching model is applied to the recorded seismic signal: in this case we consider the model with two different independent regime models, where both regimes are defined as autoregressive model of order 1, in brief AR(1), with different set of model parameters;(iii)the most well known, widely used method called STA/LTA ratio trigger algorithm, where the short-term-average (STA) window and the long-term-average (LTA) of the characteristic value of observations are calculated to detect the onset time of the seismic event wave. For the purpose of determining the effectiveness of the methods, we have decided to generate several simulated datasets, each consisting of at least one microseismic event generated using a formula of the dumped sine wave, defined as follows: where denotes the instantaneous amplitude of the modeled microseismic event at time ; is the initial amplitude of the signal envelope; denotes the decay constant; is the phase angle; and denotes the angular frequency. The example of the simulated microseismic signal generated based on the formula given above is shown in Figure 4.

To show the performance of the methods presented in the preceding section, that is, Section 2, we consider the following datasets:(i)the dataset sampled at 2 ms consisting of 2000 observations with the specified signal-to-noise ratio of 5, where the single seismic event is located in the middle of the recorded sample,(ii)the dataset sampled at 2 ms of the length of 10 seconds with several randomly located seismic events with uniformly generated signal-to-noise ratios ranging from 1 to 10. Both described datasets are presented in Figure 5.

(a)

(b)

As mentioned before we start our performance analysis by considering the simulated seismic signal record sampled with 2 ms frequency with the signal-to-noise ratio of 5, where the record sample consists of a single seismic event located in the middle of the record sample. The corresponding microseismic event has been generated according to the formula given in (16) with parameters , , , and ; see Figure 5(a). The final segments obtained by applying different segmentation methods discussed throughout this paper are shown in Figure 6. The intermediate results for the presented methods, that is, the cumulative sum of Gaussian , Markov regime-switching models with both independent autoregressive models of order 1, and the short-term-average and long-term-average (STA/LTA), are presented along Section 2 in Figures 1, 2 and 3, respectively.

Looking at the results presented in Figure 6, one may easily observe that by using considered methods we can properly identify the beginning of the segments. On the other hand the period of time where the end of the true seismic event has taken place is not uniquely identified by the discussed methods. We see that and MRS statistics produced almost exactly the same segments, whereas the reference statistics slightly deviates from the rest. This might be explained by the fact that at least one of the parameters used in the STA/LTA model (see (15)) has been wrongly selected. As a reminder, the possible set of parameters available while using STA/LTA methods consists of four different parameters, namely, and , both denoting the number of observations used for calculating appropriate moving averages; and are used as thresholds for defining the starting and ending points for each segments. Apart from that, there is also another degree of freedom while calculating and it is driven by the choice of the characteristic function . In the above considered case, the characteristic function has been defined as On the contrary the other methods require only one parameter. In case of , the expert is asked to provide a single number which is supposed to describe the variation of the noise segments. For the purpose of this simulation study, the first two hundred observations have been taken into account while calculating the value of parameter. The final value of was set to . Similarly to method, the discussed Markov regime-switching model only requires selecting appropriate threshold associated with observing the time series being in the second regime (the regime that has been driven by the seismic event). In all of our considered cases the value of the probability threshold has been set up to 0.4.

In the next part of this section, we will describe the results obtained for the simulated multievents time series; see Figure 5(b). The time series consists of five different seismic events, each of them characterized by different set of initial parameters used for generating the events according to formula given in (16). Similarly to the previous case, the intermediate results for each segmentation method are shown in Figure 7. Figure 7(a) depicts the performance of method before and after applying the smoothing procedure. The nonsmoothed version of the statistics is shown as grey line, whereas the final results obtained by applying smoothing procedure are presented in red. Comparing the results produced by method with the MRS approach, we can observe that the underlying test statistics used for determining the final segments are different. In case of the MRS method (Figure 7(b)) the nonsmoothed version the underlying test statistics is almost the same as the results obtained by applying smoothing procedure. It does not hold in the case of method where the characteristics of the nonsmoothed segments are much more volatile and there are several spikes/gaps in the estimated segments. However, after applying the smoothing procedure we obtain the same results as in the case of the MRS method. In addition to that the segments appearing one by one within the small time interval may cause some problems with the proper segmentation of the data; see Figures 7(a), 7(b), and 8. In the presented case, the fourth and fifth events are falling under these circumstances. Both presented methods, instead of producing two separated segments, just give us a single segment that has been created by merging two consecutive seismic events into the single one. The reason for that is the application of the smoothing procedure that simply assumes that the segments shorter than seconds are not considered as a true event, since they are too short to analyze and can be some artifacts rather than true events. In the case of the STA/LTA algorithm both events have been identified properly leading to the satisfactory precision of the obtained segments; see Figure 7(c). It is worth noting that the STA/LTA algorithm does not require any smoothing procedure as the obtained underlying statistics are used directly to estimate the final segments. The values of parameters used for calculating the appropriate statistics are as follows: for ; the short and long time windows are equal to 0.3 s and s, respectively; the activation and deactivation thresholds correspond to 1.05 and 0.8, respectively, in case of approach; and the probability threshold for MRS model is equal to . The performance of the discussed methods is shown in Figure 8.

(a)

(b)

(c)

Figure 7

The application of the discussed methods to the multievent time series, presented in Figure 5(b). The final segments have been marked by using bold red-dashed curves. (a) The application of the smoothing procedure used for segmentation purposes based on statistics. The grey curve depicts the nonsmoothed segments. (b) The application of the smoothing procedure to the vector of probabilities of being in the second regime with the boundary threshold set to 0.4. The grey curve denotes the nonsmoothed vector of probabilities. (c) Sequence of the test statistics , where the short and long time windows are equal to 0.3 s and 1.5 s, respectively; the activation and deactivation thresholds correspond to 1.05 and 0.8, respectively.

4. Real Data Analysis

In this section we apply the discussed segmentation techniques to the real dataset provided by KGHM Cuprum. The seismic vibration signals analyzed in this paper were acquired in a copper-ore underground mine using a prototype data acquisition system independent of the commercial seismic system in the mine. The triaxial accelerometer is installed on the mining corridor roof (see Figure 9) and measures vibration acceleration in the range of 0.001–100 m/s² with frequency range of 0.5–400 Hz. Sampling frequency is set to 1250 Hz. A rigid connection between sensor and rock is ensured, since the sensor is mounted using gypsum binder. Primarily, the sensor was located 20 m from the mining face and the distance was increasing up to 140 m, due to advance of mining works. In this paper we analyze 2 exemplary signals with sources located approximately 110–130 m from the sensor [54, 55]. The analyzed signals are presented in Figures 10 and 11.

The application of the discussed methods to real dataset consisting of single seismic event with the final segments is presented in Figure 12. In case of statistics the selected value of the parameter was set to and it is equivalent to the sample standard deviation calculated based on the 400 observations. The value of the probability threshold has not been changed and it stays at the level of . The corresponding set of parameters for STA/LTA is specified by the short and long time windows equal to 0.3 s and 1.2 s, respectively, with the activation and deactivation thresholds corresponding to 1.6 and 0.3, respectively.

(a)

(b)

(c)

Figure 12

The application of the discussed methods to real dataset consisting of single seismic event. The final segments have been marked by using bold red-dashed curves. (a) The application of the smoothing procedure used for segmentation purposes based on statistics. The grey curve depicts the nonsmoothed segments. (b) The application of the smoothing procedure to the vector of probabilities of being in the second regime with the boundary threshold set to 0.4. The grey curve denotes the nonsmoothed vector of probabilities. (c) Sequence of the test statistics , where the short and long time windows are equal to 0.3 s and 1.2 s, respectively; the activation and deactivation thresholds corresponding to 1.6 and 0.3, respectively.

As a result of applying all methods to the real data, presented in Figure 10, we can conclude that all methods produce quite consistent results. As in the case of the simulated datasets, the obtained segment calculated based on is shorter than the segments produced by the competitive models. More detailed analysis of the intermediate as well as the final results is presented in Figures 12 and 13, respectively. For Figures 12(a) and 12(b) we present the application of the smoothing procedure to the different kind of underlying test statistics that at the end leads to the conclusions. In Figure 12(a) the underlying test statistics is defined as to retrieve the information about the segments for which does not vary over the time, that is, constant. As it can be observed there are several periods of times that are too short in order to be classified as a “true” segment. The idea behind smoothing procedure is to remove such segments. The final segments obtained by applying the smoothing procedure are shown in Figures 12(a) and 12(b). With regard to Figure 12(b) the characteristics of the underlying test statistics (i.e., the vector of probabilities of being in the second regime with the boundary threshold set to 0.4) are more stable comparing to the underlying test statistics calculated based on .

In the last part of this section, we briefly discuss the results obtained for the real data, presented in Figure 11, where several seismic events of different magnitude have been recorded. The application of the smoothing procedures for and methods along with the values of is presented in Figures 14(a), 14(b), and 14(c), respectively. The final segments are depicted in Figure 15. In general we can observe that the STA/LTA approach has a tendency to overestimate the number of segments, as some of them are wrongly classified. In particular, the reader is advised to take a look at Figure 14 and verify the depicted segments around 20 seconds of the record. In general it can be explained by the fact that STA/LTA method requires few more parameters that need to be carefully selected before applying it to the real data. In contrast, the novel methods like and MRS are less demanding in terms of the number of parameters required for segments determining. Also it is worth noting that the wrongly classified segments located around 20 seconds of the records are visible in all presented methods. However, for methods like and MRS those segments are removed after applying the smoothing procedure. Apart from that we can observe that and MRS have a tendency to merge several seismic events appearing one after another in a relatively short period of time into the single blocks, whereas STA/LTA approach correctly recognizes most of the events. Such behavior is observed between 6 and 8 seconds of the recorded signal; see Figures 14 and 15.

(a)

(b)

(c)

Figure 14

The application of the discussed methods to real dataset containing several seismic events. The final segments have been marked by using bold red-dashed curves. (a) The application of the smoothing procedure used for segmentation purposes based on statistics. The grey curve depicts the nonsmoothed segments. (b) The application of the smoothing procedure to the vector of probabilities of being in the second regime with the boundary threshold set to 0.4. The grey curve denotes the nonsmoothed vector of probabilities. (c) Sequence of the test statistics , where the short and long time windows are equal to 0.3 s and 1.2 s, respectively; the activation and deactivation thresholds correspond to 1.1 and 0.8, respectively.

5. Conclusions

In this paper novel seismic signal segmentation procedures were proposed. The purpose of new methods is to isolate the seismic signal of a single event from collection of events that appear after the mining activity, like blasting procedures, provoked relaxation of rock, and some events that are unexpected, like natural rock burst. The proposed techniques are based on stochastic properties of the examined signals in time domain. Both invented methods assume the noisy part of the seismic signal can be described by the Gaussian distribution. However, in the first approach we assume the dataset constitutes sample of independent identically distributed random variables while in the second approach we assume short time dependency of the time series. Moreover, in the second case we also specify the model of the nonnoisy part of the signal which is related to the seismic event. We applied the proposed methods to the simulated as well as real seismic signals acquired in a copper-ore underground mine. We compared the results obtained by application of two new segmentation techniques to the classical method of seismic signal segmentation based on the STA/LTA ratio. We should mention that the proposed segmentation methods give similar results and especially for the multievents signals they seem to be more robust as the classical technique. After analysis of real seismic signals we can conclude that the STA/LTA approach has a tendency to overestimate the number of segments and moreover some of them are wrongly classified. At the end we should mention that although the presented procedures were dedicated to seismic signals segmentation they can be also applied to other applications.

Competing Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

References

J. Gajda, G. Sikora, and A. Wylomanska, “Regime variance testing—a quantile approach,” Acta Physica Polonica B, vol. 44, no. 5, p. 1015, 2013.
View at: Publisher Site | Google Scholar
R. S. Tsay, “Outliers, level shifts, and variance changes in time series,” Journal of Forecasting, vol. 7, no. 1, pp. 1–20, 1988.
View at: Publisher Site | Google Scholar
M. Lopatka, C. Laplanche, O. Adam, J.-F. Motsch, and J. Zarzycki, “Non-stationary time-series segmentation based on the Schur prediction error analysis,” in Proceedings of the IEEE/SP 13th Workshop on Statistical Signal Processing, pp. 251–256, Novosibirsk, Russia, July 2005.
View at: Publisher Site | Google Scholar
R. Makowski and R. Zimroz, “A procedure for weighted summation of the derivatives of reflection coefficients in adaptive Schur filter with application to fault detection in rolling element bearings,” Mechanical Systems and Signal Processing, vol. 38, no. 1, pp. 65–77, 2013.
View at: Publisher Site | Google Scholar
R. Makowski and R. Zimroz, “New techniques of local damage detection in machinery based on stochastic modelling using adaptive Schur filter,” Applied Acoustics, vol. 77, pp. 130–137, 2014.
View at: Publisher Site | Google Scholar
C. Li, M. Liang, and T. Wang, “Criterion fusion for spectral segmentation and its application to optimal demodulation of bearing vibration signals,” Mechanical Systems and Signal Processing, vol. 64-65, pp. 132–148, 2015.
View at: Publisher Site | Google Scholar
C. Bhagavatula, A. Jaech, M. Savvides et al., “Automatic segmentation of cardiosynchronous waveforms using cepstral analysis and continuous wavelet transforms,” in Proceedings of the 19th IEEE International Conference on Image Processing (ICIP '12), pp. 2045–2048, Orlando, Fla, USA, October 2012.
View at: Publisher Site | Google Scholar
T. D. Popescu and D. Aiordachioaie, “Signal segmentation in time-frequency plane using Renyi entropy—application in seismic signal processing,” in Proceedings of the Conference on Control and Fault-Tolerant Systems, pp. 312–317, Nice, France, October 2013.
View at: Publisher Site | Google Scholar
J. Obuchowski, A. Wyłomańska, and R. Zimroz, “The local maxima method for enhancement of time–frequency map and its application to local damage detection in rotating machines,” Mechanical Systems and Signal Processing, vol. 46, no. 2, pp. 389–405, 2014.
View at: Publisher Site | Google Scholar
J. Urbanek, T. Barszcz, R. Zimroz, and J. Antoni, “Application of averaged instantaneous power spectrum for diagnostics of machinery operating under non-stationary operational conditions,” Measurement, vol. 45, no. 7, pp. 1782–1791, 2012.
View at: Publisher Site | Google Scholar
R. Zimroz, M. Madziarz, G. Zak, A. Wylomanska, and J. Obuchowski, “Seismic signal segmentation procedure using time-frequency decomposition and statistical modelling,” Journal of Vibroengineering, vol. 17, no. 6, pp. 3111–3121, 2015.
View at: Google Scholar
M. Musaruddin and R. Zivanovic, “Signal segmentation of fault records based on Empirical Mode Decomposition,” in Proceedings of the IEEE Region 10 Conference (TENCON '11), pp. 138–143, Bali, Indonesia, November 2011.
View at: Publisher Site | Google Scholar
H. Azami and S. Sanei, “Automatic signal segmentation based on singular spectrum analysis and imperialist competitive algorithm,” in Proceedings of the 2nd International e-Conference on Computer and Knowledge Engineering (ICCKE '12), pp. 50–55, Mashhad, Iran, October 2012.
View at: Publisher Site | Google Scholar
J. A. Crossman, H. Guo, Y. L. Murphey, and J. Cardillo, “Automotive signal fault diagnostics—part I: signal fault analysis, signal segmentation, feature extraction and quasi-optimal feature selection,” IEEE Transactions on Vehicular Technology, vol. 52, no. 4, pp. 1063–1075, 2003.
View at: Publisher Site | Google Scholar
A. Wylomanska and R. Zimroz, “Signal segmentation for operational regimes detection of heavy duty mining mobile machines-a statistical approach,” Diagnostyka, vol. 15, no. 2, pp. 33–42, 2014.
View at: Google Scholar
K. Burnecki, A. Wyłomańska, A. Beletskii, V. Gonchar, and A. Chechkin, “Recognition of stable distribution with Lévy index α close to 2,” Physical Review E, vol. 85, no. 5, Article ID 056711, 2012.
View at: Publisher Site | Google Scholar
B. Tóth, F. Lillo, and J. D. Farmer, “Segmentation algorithm for non-stationary compound Poisson processes,” The European Physical Journal B, vol. 78, no. 2, pp. 235–243, 2010.
View at: Publisher Site | Google Scholar
H. Azami, K. Mohammadi, and B. Bozorgtabar, “An improved signal segmentation using moving average and Savitzky-Golay filter,” Journal of Signal and Information Processing, vol. 3, pp. 39–44, 2012.
View at: Publisher Site | Google Scholar
J. Terrien, G. Germain, C. Marque, and B. Karlsson, “Bivariate piecewise stationary segmentation; improved pre-treatment for synchronization measures used on non-stationary biological signals,” Medical Engineering & Physics, vol. 35, no. 8, pp. 1188–1196, 2013.
View at: Publisher Site | Google Scholar
R. V. Andreão, B. Dorizzi, and J. Boudy, “ECG signal analysis through hidden Markov models,” IEEE Transactions on Biomedical Engineering, vol. 53, no. 8, pp. 1541–1549, 2006.
View at: Publisher Site | Google Scholar
H. J. L. M. Vullings, M. H. G. Verhaegen, and H. B. Verbruggen, “Automated ECG segmentation with dynamic time warping,” in Proceedings of the 20th Annual International Conference of the IEEE Engineering in Medicine and Biology Society, vol. 1, pp. 163–166, Hong Kong, November 1998.
View at: Publisher Site | Google Scholar
P. Micó, M. Mora, D. Cuesta-Frau, and M. Aboy, “Automatic segmentation of long-term ECG signals corrupted with broadband noise based on sample entropy,” Computer Methods and Programs in Biomedicine, vol. 98, no. 2, pp. 118–129, 2010.
View at: Publisher Site | Google Scholar
S. Choi and Z. Jiang, “Comparison of envelope extraction algorithms for cardiac sound signal segmentation,” Expert Systems with Applications, vol. 34, no. 2, pp. 1056–1069, 2008.
View at: Publisher Site | Google Scholar
B. Lovell and B. Boashash, “Segmentation of non-stationary signals with applications,” in Proceedings of the International Conference on Acoustics, Speech, and Signal Processing (CASSP '88), Cat. No.88CH2561-9, New York, NY, USA, April 1988.
View at: Google Scholar
V. Khanagha, K. Daoudi, O. Pont, and H. Yahia, “Phonetic segmentation of speech signal using local singularity analysis,” Digital Signal Processing, vol. 35, pp. 86–94, 2014.
View at: Publisher Site | Google Scholar
R. Makowski and R. Hossa, “Automatic speech signal segmentation based on the innovation adaptive filter,” International Journal of Applied Mathematics and Computer Science, vol. 24, no. 2, pp. 259–270, 2014.
View at: Publisher Site | Google Scholar | Zentralblatt MATH
J. Janczura, “Pricing electricity derivatives within a Markov regime-switching model: a risk premium approach,” Mathematical Methods of Operations Research, vol. 79, no. 1, pp. 1–30, 2014.
View at: Publisher Site | Google Scholar
J. Janczura and R. Weron, “Goodness-of-fit testing for the marginal distribution of regime-switching models with an application to electricity spot prices,” AStA Advances in Statistical Analysis, vol. 97, no. 3, pp. 239–270, 2013.
View at: Publisher Site | Google Scholar
X. Xu, L. Dou, C. Lu, and Y. Zhang, “Frequency spectrum analysis on micro-seismic signal of rock bursts induced by dynamic disturbance,” Mining Science and Technology, vol. 20, no. 5, pp. 682–685, 2010.
View at: Publisher Site | Google Scholar
C. H. Chen, “On a segmentation algorithm for seismic signal analysis,” Geoexploration, vol. 23, no. 1, pp. 35–40, 1984.
View at: Publisher Site | Google Scholar
J. E. Gaby and K. R. Anderson, “Hierarchical segmentation of seismic waveforms using affinity,” Geoexploration, vol. 23, no. 1, pp. 1–16, 1984.
View at: Publisher Site | Google Scholar
T. D. Popescu, “Signal segmentation using changing regression models with application in seismic engineering,” Digital Signal Processing, vol. 24, pp. 14–26, 2014.
View at: Publisher Site | Google Scholar
A. Leśniak and Z. Isakow, “Space-time clustering of seismic events and hazard assessment in the Zabrze-Bielszowice coal mine, Poland,” International Journal of Rock Mechanics and Mining Sciences, vol. 46, no. 5, pp. 918–928, 2009.
View at: Publisher Site | Google Scholar
Z. Haiyang, H. Zhonghua, and C. Guoxing, “Numerical modeling on the seismic responses of a large underground structure in soft ground,” Journal of Vibroengineering, vol. 17, no. 2, pp. 802–815, 2015.
View at: Google Scholar
L. R. Rabiner, “A tutorial on hidden Markov models and selected applications in speech recognition,” Proceedings of the IEEE, vol. 77, no. 2, pp. 257–286, 1989.
View at: Publisher Site | Google Scholar
N. Ravishanker, Z. Liu, and B. K. Ray, “NHPP models with Markov switching for software reliability,” Computational Statistics & Data Analysis, vol. 52, no. 8, pp. 3988–3999, 2008.
View at: Publisher Site | Google Scholar | MathSciNet
P. Bieńkowski, K. Burnecki, J. Janczura, R. Weron, and B. Zubrzak, “A new method for automated noise cancellation in electromagnetic field measurement,” Journal of Electromagnetic Waves and Applications, vol. 26, no. 8-9, pp. 1226–1236, 2012.
View at: Publisher Site | Google Scholar
L. Antwarg, L. Rokach, and B. Shapira, “Attribute-driven hidden markov model trees for intention prediction,” IEEE Transactions on Systems, Man and Cybernetics Part C, vol. 42, no. 6, pp. 1103–1119, 2012.
View at: Publisher Site | Google Scholar
X. Cui, J. Xue, X. Chen et al., “Hidden markov acoustic modeling with bootstrap and restructuring for low-resourced languages,” IEEE Transactions on Audio, Speech and Language Processing, vol. 20, no. 8, pp. 2252–2264, 2012.
View at: Publisher Site | Google Scholar
C.-L. Liu, W.-H. Hsaio, C.-H. Lee, and H.-C. Chi, “An HMM-based algorithm for content ranking and coherence-feature extraction,” IEEE Transactions on Systems, Man, and Cybernetics: Systems, vol. 43, no. 2, pp. 440–450, 2013.
View at: Publisher Site | Google Scholar
X. Ma, D. Schonfeld, and A. A. Khokhar, “Video event classification and image segmentation based on noncausal multidimensional hidden Markov models,” IEEE Transactions on Image Processing, vol. 18, no. 6, pp. 1304–1313, 2009.
View at: Publisher Site | Google Scholar | MathSciNet
A. Sloin and D. Burshtein, “Support vector machine training for improved hidden Markov modeling,” IEEE Transactions on Signal Processing, vol. 56, no. 1, pp. 172–188, 2008.
View at: Publisher Site | Google Scholar | MathSciNet
T. Boutros and M. Liang, “Detection and diagnosis of bearing and cutting tool faults using hidden Markov models,” Mechanical Systems and Signal Processing, vol. 25, no. 6, pp. 2102–2124, 2011.
View at: Publisher Site | Google Scholar
M. Dong and D. He, “A segmental hidden semi-Markov model (HSMM)-based diagnostics and prognostics framework and methodology,” Mechanical Systems and Signal Processing, vol. 21, no. 5, pp. 2248–2266, 2007.
View at: Publisher Site | Google Scholar
M. Ge, R. Du, and Y. Xu, “Hidden Markov model based fault diagnosis for stamping processes,” Mechanical Systems and Signal Processing, vol. 18, no. 2, pp. 391–408, 2004.
View at: Publisher Site | Google Scholar
Z. Li, Z. Wu, Y. He, and C. Fulei, “Hidden Markov model-based fault diagnostics method in speed-up and speed-down process for rotating machinery,” Mechanical Systems and Signal Processing, vol. 19, no. 2, pp. 329–339, 2005.
View at: Publisher Site | Google Scholar
Q. Miao and V. Makis, “Condition monitoring and classification of rotating machinery using wavelets and hidden Markov models,” Mechanical Systems and Signal Processing, vol. 21, no. 2, pp. 840–855, 2007.
View at: Publisher Site | Google Scholar
R. V. Allen, “Automatic earthquake recognition and timing from single traces,” Bulletin of the Seismological Society of America, vol. 68, no. 5, pp. 1521–1532, 1978.
View at: Google Scholar
L. Dong and X. Li, “A microseismic/acoustic emission source location method using arrival times of PS waves for unknown velocity system,” International Journal of Distributed Sensor Networks, vol. 2013, Article ID 307489, 8 pages, 2013.
View at: Publisher Site | Google Scholar
X. Li and L. Dong, “An efficient closed-form solution for acoustic emission source location in three-dimensional structures,” AIP Advances, vol. 4, no. 2, Article ID 027110, 2014.
View at: Publisher Site | Google Scholar
L. Dong, X. Li, and G. Xie, “An analytical solution for acoustic emission source location for known P wave velocity system,” Mathematical Problems in Engineering, vol. 2014, Article ID 290686, 6 pages, 2014.
View at: Publisher Site | Google Scholar
J. Janczura and R. Weron, “Efficient estimation of Markov regime-switching models: an application to electricity spot prices,” AStA Advances in Statistical Analysis, vol. 96, no. 3, pp. 385–407, 2012.
View at: Google Scholar
R. Huisman and C. de Jong, Option Formulas for Mean-Reverting Power Prices Withspikes, ERIM Report Series Reference no. ERS 2002-96, F & A, 2002.
M. Polak, J. Obuchowski, M. Madziarz, A. Wyłomańska, and R. Zimroz, “Time-varying group delay as a basis for clustering and segmentation of seismic signals,” Journal of Vibroengineering, vol. 18, no. 1, pp. 267–275, 2016.
View at: Google Scholar
J. Sokolowski, J. Obuchowski, M. Madziarz, A. Wylomanska, and R. Zimroz, “Seismic signals discrimination based on instantaneous frequency,” Vibroengineering Procedia, vol. 6, pp. 200–205, 2015.
View at: Google Scholar

Copyright

Copyright © 2016 Daniel Kucharczyk et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

PDF Download Citation

Download other formats

Order printed copies

Views

1435

Downloads

972

Citations