# Stochastic Modelling as a Tool for Seismic Signals Segmentation

Authors: Daniel Kucharczyk; Agnieszka Wyłomańska; Jakub Obuchowski; Radosław Zimroz; Maciej Madziarz
Journal: Shock and Vibration
 (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.


---

## Body

## 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 Section2 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 statisticMLt, defined as a cumulative sum of zero-mean Gaussian probability density functions; that is,(1)MLt=∑k=1tϕxt,t=1,2,…,T, where ϕ(xk) is the value of Gaussian probability density function given at point xk, with a zero mean (μ=0) and relatively small variance (σ<1) given as follows:(2)ϕxk=12πσexp⁡-xk22σ2, where xt denotes the value of the observed signal at time t. 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 statisticsMLt given in (1), it is clear that, for each point xk 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 MLt 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.Figure 1
The application of the proposed method to the simulated data record using the signal described in Section3. (a) Detection of the seismic signal in the recorded dataset. (b) Sequence of the test statistics MLt. (c) The application of the smoothing procedure used for segmentation purposes based on MLt statistics. The grey curves depict the nonsmoothed differences between consecutive values of the test statistics, whereas the red curve is used for marking the final segments.
(a)
(b)
(c)After plotting the values of test statisticsMLt for each point xt, it is clearly visible that the period of times characterized by the constant value of MLt 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 than5 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 than5 steps then we fill the output vector with zeros at this location. The last step is to remove segments shorter than 0.4 seconds (under the assumption that the sampling interval of the record is equal to 2ms, 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 chainRt. In general, the Markov chain Rt is described by the transition matrix P containing the probabilities(3)pij=PRt+1=j∣Rt=iof switching from regime i at time t to regime j at time t+1. Assuming that there are only two regimes, i,j=1,2, we have(4)P=pij=p11p12p21p22. 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:(5)Xt=αRt+βRtXt-1+σRtϵt, where ϵt are assumed to be Gaussian distributed N(0,1), that is, with parameters μ and σ equal zero and unity, respectively. In the second case with independent regimes, the underlying time series Xt is defined as(6)Xt=Xt,1,ifRt=1,Xt,2,ifRt=2, where at least one regime is defined as AR(1) model given in the form(7)Xt,i=αi+βiXt-1,i+σiϵt,i.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, θ(0)=ϕi(0),P(0),ρi(0) for i=1,2, where ϕi=(αi,βi,σi), denotes parameters of ith regime distribution and ρi(0)≡P(R1=i) 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 calledE-part estimates conditional probabilities of being in regime i at time t, since switching mechanism Rt is not directly observable and we can only derive expected value of the process EIRt=i∣x1,…,xT;θ which is mentioned conditional probability PRt=i∣x1,…,xT;θ, where x1,…,xT denotes the realizations of the process Xt. The results from this step are used in theM-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 θ(n) from previous iteration ((n) denote the nth iteration of the procedure) and a vector xt=(x1,…,xt) of realization of the process Xt up to given point time t. This part is divided into two steps.(i) Filtering. At this step we derive the probabilities with given time t in the iterative equation. Let t=1,2,…,T; then(8)PRt=i∣xt;θn=PRt=i∣xt-1;θnfxt∣Rt=i;xt-1;θn∑j=12PRt=j∣xt-1;θnfxt∣Rt=j;xt-1;θn,PRt+1=i∣xt;θn=∑j=12pjinPRt=j∣xt;θn, where fxt∣Rt=i;xt-1;θn is the conditional density of the regime i process at time t. The starting point of this iteration is given as PR1=i∣x0;θn=ρi(n).(ii) Smoothing. In this part we are calculating conditional probabilities with given xT; this means that we are conditioning with whole history of the process xt. This step goes backward, t=T-1, T-2,…,1 and iterates on(9)PRt=i∣xT;θn=∑j=12PRt=i∣xt;θnPRt+1=j∣xT;θnpijnPRt+1=j∣xt;θn.To calculate probabilities in (9) we need a derivation of f(xt∣Rt=i;xt-1;θ(n)). For the first model which is AR(1) process with switching parameters (dependent regimes) and the same Gaussian random noise we use conditional distribution formula given as(10)fxt∣Rt=i;xt-1;θn=12πσinexp⁡-xt-βinxt-1-αin22σin2. For the model described in (6), with two independent random variables we use simple formula(11)fxt∣Rt=i;xt-1;θn=gixt, where gi(x) is a density function of a chosen random variable in regime i, since this random variable is independent and identically distributed.The calculation becomes more complicated when at least one of the regimes is given byAR(1) process and is independent of the other ones. As the regime process Rt is latent and unobserved, the way of conditioning Xt with given Xt-1 makes it dependent on the whole history of the process and in this case all possible paths of this process up to time t should be considered. This gives us 2t 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 xt-1 in formula (10) with their expectation x~t-1=EXt-1∣xt-1;θn, which contains all information about process xt until time t-1. This value can be computed with the following formula:(12)x~t-1,i=EXt,i∣xt;θn=PRt=i∣xt;θnxt+PRt≠i∣xt;θnαin+1-βinEXt-1,i∣xt-1;θn. Using formula (12) one can rapidly reduce computational complexity of this algorithm without losing quality of the estimation.M-Step. The M-step estimates both regime distribution parameters θ(n+1) using the modified maximum likelihood method, the same as in standard algorithm of maximum likelihood method of estimation where we use log-likelihood function ∑t=1Tln⁡f(xt,θ(n)), 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 ith regime parameters leads to the following formulas:(13)α^i=∑t=2TPRt=i∣xT;θnxt-β^ixt-1∑t=2TPRt=i∣xT;θn,β^i=B1-B2B3,B1=∑t=2TPRt=i∣xT;θnxt-1xt-xt-1,B2=∑t=2TPRt=i∣xT;θnxt-xt-1∑t=2TPRt=i∣xT;θn·∑t=2TPRt=i∣xT;θnxt-1,B3=∑t=2TPRt=i∣xT;θnxt-12∑t=2TPRt=i∣xT;θn-∑t=2TPRt=i∣xT;θnxt-12,σ^i2=∑t=2TPRt=i∣xT;θnxt-α^i-β^ixt-12∑t=2TPRt=i∣xT;θn.For the model defined in (6) where at least one of the regimes is independent of AR(1) process we can use the estimators given above, but we need to replace observation xt-1 with its expectation x~t-1,i 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(14)μ^=∑t=2TxtPRt=i∣xT;θn∑t=2TPRt=i∣xT;θn,σ^2=∑t=2Txt-μ^2PRt=i∣xT;θn∑t=2TPRt=i∣xT;θn.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 than5 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 than5 steps then we fill the output vector with zeros at this location. The last step is to remove segments shorter than 0.4 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.Figure 2
The application of the proposed method to the simulated signal from Section3. (a) Detection of the seismic signal in the recorded dataset. (b) Sequence of the smoothed conditional probabilities of observing the second regime PRt=2∣xt;θn (red line) including the nonsmoothed vector of PRt=2∣xt;θn marked as grey line.
(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 τA (activation threshold), the consecutive recorded values are marked to be event-driven until the ratio falls below another predefined value τD, called the deactivation threshold. In general the STA/LTA ratio algorithm can be described as follows:(15)SLRt=1/ns∑k=t-nstCFxk1/nl∑k=t-nltCFxk,t=1,2,…,T, where ns and nl 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 τA and deactivation τD thresholds correspond to 1.4 and 1.2, respectively.Figure 3
The application of the proposed method to a simulated signal described in Section3. (a) Detection of the signal segment in the recorded dataset. (b) Sequence of the test statistics SLRt, where the short and long time windows are equal to 0.3 s and 1.5 s, respectively; the activation τA and deactivation τD thresholds correspond to 1.4 and 1.2, respectively.
(a)
(b)


## 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 statisticMLt, defined as a cumulative sum of zero-mean Gaussian probability density functions; that is,(1)MLt=∑k=1tϕxt,t=1,2,…,T, where ϕ(xk) is the value of Gaussian probability density function given at point xk, with a zero mean (μ=0) and relatively small variance (σ<1) given as follows:(2)ϕxk=12πσexp⁡-xk22σ2, where xt denotes the value of the observed signal at time t. 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 statisticsMLt given in (1), it is clear that, for each point xk 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 MLt 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.Figure 1
The application of the proposed method to the simulated data record using the signal described in Section3. (a) Detection of the seismic signal in the recorded dataset. (b) Sequence of the test statistics MLt. (c) The application of the smoothing procedure used for segmentation purposes based on MLt statistics. The grey curves depict the nonsmoothed differences between consecutive values of the test statistics, whereas the red curve is used for marking the final segments.
(a)
(b)
(c)After plotting the values of test statisticsMLt for each point xt, it is clearly visible that the period of times characterized by the constant value of MLt 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 than5 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 than5 steps then we fill the output vector with zeros at this location. The last step is to remove segments shorter than 0.4 seconds (under the assumption that the sampling interval of the record is equal to 2ms, 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 chainRt. In general, the Markov chain Rt is described by the transition matrix P containing the probabilities(3)pij=PRt+1=j∣Rt=iof switching from regime i at time t to regime j at time t+1. Assuming that there are only two regimes, i,j=1,2, we have(4)P=pij=p11p12p21p22. 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:(5)Xt=αRt+βRtXt-1+σRtϵt, where ϵt are assumed to be Gaussian distributed N(0,1), that is, with parameters μ and σ equal zero and unity, respectively. In the second case with independent regimes, the underlying time series Xt is defined as(6)Xt=Xt,1,ifRt=1,Xt,2,ifRt=2, where at least one regime is defined as AR(1) model given in the form(7)Xt,i=αi+βiXt-1,i+σiϵt,i.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, θ(0)=ϕi(0),P(0),ρi(0) for i=1,2, where ϕi=(αi,βi,σi), denotes parameters of ith regime distribution and ρi(0)≡P(R1=i) 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 calledE-part estimates conditional probabilities of being in regime i at time t, since switching mechanism Rt is not directly observable and we can only derive expected value of the process EIRt=i∣x1,…,xT;θ which is mentioned conditional probability PRt=i∣x1,…,xT;θ, where x1,…,xT denotes the realizations of the process Xt. The results from this step are used in theM-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 θ(n) from previous iteration ((n) denote the nth iteration of the procedure) and a vector xt=(x1,…,xt) of realization of the process Xt up to given point time t. This part is divided into two steps.(i) Filtering. At this step we derive the probabilities with given time t in the iterative equation. Let t=1,2,…,T; then(8)PRt=i∣xt;θn=PRt=i∣xt-1;θnfxt∣Rt=i;xt-1;θn∑j=12PRt=j∣xt-1;θnfxt∣Rt=j;xt-1;θn,PRt+1=i∣xt;θn=∑j=12pjinPRt=j∣xt;θn, where fxt∣Rt=i;xt-1;θn is the conditional density of the regime i process at time t. The starting point of this iteration is given as PR1=i∣x0;θn=ρi(n).(ii) Smoothing. In this part we are calculating conditional probabilities with given xT; this means that we are conditioning with whole history of the process xt. This step goes backward, t=T-1, T-2,…,1 and iterates on(9)PRt=i∣xT;θn=∑j=12PRt=i∣xt;θnPRt+1=j∣xT;θnpijnPRt+1=j∣xt;θn.To calculate probabilities in (9) we need a derivation of f(xt∣Rt=i;xt-1;θ(n)). For the first model which is AR(1) process with switching parameters (dependent regimes) and the same Gaussian random noise we use conditional distribution formula given as(10)fxt∣Rt=i;xt-1;θn=12πσinexp⁡-xt-βinxt-1-αin22σin2. For the model described in (6), with two independent random variables we use simple formula(11)fxt∣Rt=i;xt-1;θn=gixt, where gi(x) is a density function of a chosen random variable in regime i, since this random variable is independent and identically distributed.The calculation becomes more complicated when at least one of the regimes is given byAR(1) process and is independent of the other ones. As the regime process Rt is latent and unobserved, the way of conditioning Xt with given Xt-1 makes it dependent on the whole history of the process and in this case all possible paths of this process up to time t should be considered. This gives us 2t 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 xt-1 in formula (10) with their expectation x~t-1=EXt-1∣xt-1;θn, which contains all information about process xt until time t-1. This value can be computed with the following formula:(12)x~t-1,i=EXt,i∣xt;θn=PRt=i∣xt;θnxt+PRt≠i∣xt;θnαin+1-βinEXt-1,i∣xt-1;θn. Using formula (12) one can rapidly reduce computational complexity of this algorithm without losing quality of the estimation.M-Step. The M-step estimates both regime distribution parameters θ(n+1) using the modified maximum likelihood method, the same as in standard algorithm of maximum likelihood method of estimation where we use log-likelihood function ∑t=1Tln⁡f(xt,θ(n)), 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 ith regime parameters leads to the following formulas:(13)α^i=∑t=2TPRt=i∣xT;θnxt-β^ixt-1∑t=2TPRt=i∣xT;θn,β^i=B1-B2B3,B1=∑t=2TPRt=i∣xT;θnxt-1xt-xt-1,B2=∑t=2TPRt=i∣xT;θnxt-xt-1∑t=2TPRt=i∣xT;θn·∑t=2TPRt=i∣xT;θnxt-1,B3=∑t=2TPRt=i∣xT;θnxt-12∑t=2TPRt=i∣xT;θn-∑t=2TPRt=i∣xT;θnxt-12,σ^i2=∑t=2TPRt=i∣xT;θnxt-α^i-β^ixt-12∑t=2TPRt=i∣xT;θn.For the model defined in (6) where at least one of the regimes is independent of AR(1) process we can use the estimators given above, but we need to replace observation xt-1 with its expectation x~t-1,i 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(14)μ^=∑t=2TxtPRt=i∣xT;θn∑t=2TPRt=i∣xT;θn,σ^2=∑t=2Txt-μ^2PRt=i∣xT;θn∑t=2TPRt=i∣xT;θn.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 than5 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 than5 steps then we fill the output vector with zeros at this location. The last step is to remove segments shorter than 0.4 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.Figure 2
The application of the proposed method to the simulated signal from Section3. (a) Detection of the seismic signal in the recorded dataset. (b) Sequence of the smoothed conditional probabilities of observing the second regime PRt=2∣xt;θn (red line) including the nonsmoothed vector of PRt=2∣xt;θn marked as grey line.
(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 τA (activation threshold), the consecutive recorded values are marked to be event-driven until the ratio falls below another predefined value τD, called the deactivation threshold. In general the STA/LTA ratio algorithm can be described as follows:(15)SLRt=1/ns∑k=t-nstCFxk1/nl∑k=t-nltCFxk,t=1,2,…,T, where ns and nl 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 τA and deactivation τD thresholds correspond to 1.4 and 1.2, respectively.Figure 3
The application of the proposed method to a simulated signal described in Section3. (a) Detection of the signal segment in the recorded dataset. (b) Sequence of the test statistics SLRt, where the short and long time windows are equal to 0.3 s and 1.5 s, respectively; the activation τA and deactivation τD thresholds correspond to 1.4 and 1.2, respectively.
(a)
(b)


## 3. Simulations

In this section we describe the results of performed simulation study of the methods presented in Section2. 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μ=0 and σ≪1;(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:(16)Xt=Ae-λtsin⁡ωt+θ, where Xt denotes the instantaneous amplitude of the modeled microseismic event at time t; A 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.Figure 4
Microseismic event simulated using the formula of the dumped sine wave given in (16) with parameters A=20,  λ=4,  ω=80π, and θ=0, where the dominant frequency is 40Hz and the sample rate is 2ms.To show the performance of the methods presented in the preceding section, that is, Section2, 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.Figure 5
(a) Simulated time series containing single seismic event generated according to (16) with parameters A=20,  λ=4, ω=80π, and θ=0, with signal-to-noise ratio equal to 25. (b) Simulated time series consisting of several seismic events where each seismic event is driven by different set of parameters used in (16).
(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 A=20,  λ=4,  ω=80π, and θ=0; 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 MLt, 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.Figure 6
The comparison of segments produced by different methods applied to the simulated multievents dataset.Looking at the results presented in Figure6, 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 MLt and MRS statistics produced almost exactly the same segments, whereas the reference SLRt 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, ns and nl, both denoting the number of observations used for calculating appropriate moving averages; τD and τA 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 SLRt and it is driven by the choice of the characteristic function CF. In the above considered case, the characteristic function has been defined as(17)CFxk=xk2,k=1,2,…,n. On the contrary the other methods require only one parameter. In case of MLt, 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 0.1353. Similarly to MLt 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 Figure5(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 MLt 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 MLt 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 MLt 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 0.4 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 SLRt are used directly to estimate the final segments. The values of parameters used for calculating the appropriate statistics are as follows: σ=0.0570 for MLt; the short ns and long nl time windows are equal to 0.3 s and 1.5 s, respectively; the activation τA and deactivation τD thresholds correspond to 1.05 and 0.8, respectively, in case of SLRt approach; and the probability threshold for MRS model is equal to 0.4. The performance of the discussed methods is shown in Figure 8.Figure 7
The application of the discussed methods to the multievent time series, presented in Figure5(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 MLt 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 SLRt, where the short and long time windows are equal to 0.3 s and 1.5 s, respectively; the activation τA and deactivation τD thresholds correspond to 1.05 and 0.8, respectively.
(a)
(b)
(c)Figure 8
The comparison of segments produced by different methods applied to the simulated multievents dataset.


## 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 Figure9) and measures vibration acceleration in the range of 0.001–100 m/s2 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.Figure 9
The accelerometer located on the mining corridor roof.Figure 10
The analyzed real seismic signal with single event.Figure 11
The analyzed real seismic signal with multiple events.The application of the discussed methods to real dataset consisting of single seismic event with the final segments is presented in Figure12. In case of MLt statistics the selected value of the parameter σ was set to 0.0570 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 0.4. 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 τA and deactivation τD thresholds corresponding to 1.6 and 0.3, respectively.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 onMLt 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 SLRt, where the short and long time windows are equal to 0.3 s and 1.2 s, respectively; the activation τA and deactivation τD thresholds corresponding to 1.6 and 0.3, respectively.
(a)
(b)
(c)As a result of applying all methods to the real data, presented in Figure10, we can conclude that all methods produce quite consistent results. As in the case of the simulated datasets, the obtained segment calculated based on SLRt 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 MLt-MLt-1=0 to retrieve the information about the segments for which MLt 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 MLt.Figure 13
The comparison of the discussed segmentation methods applied to the real data consisting of a single seismic event.In the last part of this section, we briefly discuss the results obtained for the real data, presented in Figure11, where several seismic events of different magnitude have been recorded. The application of the smoothing procedures for MLt and MRS methods along with the values of SLRt 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 MLt 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 MLt and MRS those segments are removed after applying the smoothing procedure. Apart from that we can observe that MLt 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.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 onMLt 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 SLRt, where the short and long time windows are equal to 0.3 s and 1.2 s, respectively; the activation τA and deactivation τD thresholds correspond to 1.1 and 0.8, respectively.
(a)
(b)
(c)Figure 15
The comparison of the discussed segmentation methods applied to the real data containing several seismic events.


## 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.



---

*Source: 8453426-2016-05-08.xml*