Stochastic Modelling as a Tool for Seismic Signals Segmentation

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.


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][4][5].There are also methods which are based on the properties of the signal in other domains [6][7][8][9][10][11]. In the literature one can also find segmentation techniques based on the time series decomposition [12][13][14].
The problem of extraction of parts of the time series with homogeneous properties is also important in analysis of seismic signals [8,[29][30][31][32][33][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 regimeswitching 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][39][40][41][42].This model was also used in diagnostics [43][44][45][46][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][50][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.

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.

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 ML  , 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 ( = 0) and relatively small variance ( < 1) 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 ML  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 ML  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.
After plotting the values of test statistics ML  for each point   , it is clearly visible that the period of times characterized by the constant value of ML  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 5 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 5 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 2 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).
Shock and Vibration

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 P containing the probabilities of switching from regime  at time  to regime  at time  + 1.
Assuming that there are only two regimes, ,  = {1, 2}, 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 N(0, 1), 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,  (0) = ( (0)  , P (0) ,  (0)  ) for  = 1, 2, where   = (  ,   ,   ), denotes parameters of th regime distribution and  (0)  ≡ ( 1 = ) 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 (I   = |  1 , . . .,   ; ) which is mentioned conditional probability (  =  |  1 , . . .,   ; ), where  1 , . . .,   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 x  = ( 1 , . . .,   ) 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  = 1, 2, . . ., ; then where (  |   = ; x −1 ;  () ) 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 x  ; this means that we are conditioning with whole history of the process   .This step goes backward,  =  − 1,  − 2, . . ., 1 and iterates on To calculate probabilities in (9) we need a derivation of (  |   = ; x −1 ;  () ).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 ) . ( 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 AR(1) process and is independent of the other ones.As the regime process   is latent and unobserved, the way of conditioning   with given  −1 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 2  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  −1 in formula (10) with their expectation x−1 = ( −1 | x −1 ;  () ), which contains all information about process   until time −1.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  (+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 ∑  =1 ln (  ,  () ), 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 AR(1) process we can use the estimators given above, but we need to replace observation  −1 with its expectation x−1, 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 Shock and Vibration 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 5 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 5 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.

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 1.4 and 1.2, respectively.

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  = 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 shortterm-average (STA) window and the long-termaverage (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.
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  = 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 ML  , 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 ML  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 CF.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 ML  , 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 ML  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 ML  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 ML  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 ML  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   are used directly to estimate the final segments.The values of parameters used for calculating the appropriate statistics are as follows:  = 0.0570 for ML  ; 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, in case of   approach; and the probability threshold for MRS model is equal to 0.4.The performance of the discussed methods is shown in Figure 8.

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 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 ML  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   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  defined as ML  − ML −1 = 0 to retrieve the information about the segments for which ML  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 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 ML  and MRS 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 ML  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 ML  and MRS those segments are removed after applying the smoothing procedure.Apart from that we can observe that ML  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.

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

Figure 1 :
Figure 1: The application of the proposed method to the simulated data record using the signal described in Section 3. (a) Detection of the seismic signal in the recorded dataset.(b) Sequence of the test statistics ML  .(c) The application of the smoothing procedure used for segmentation purposes based on ML  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.

Figure 2 :
Figure 2: The application of the proposed method to the simulated signal from Section 3. (a) Detection of the seismic signal in the recorded dataset.(b) Sequence of the smoothed conditional probabilities of observing the second regime (  = 2 | x  ;  () ) (red line) including the nonsmoothed vector of (  = 2 | x  ;  () ) marked as grey line.

Figure 3 :Figure 4 :
Figure 3: The application of the proposed method to a simulated signal described in Section 3. (a) Detection of the signal segment in the recorded dataset.(b) 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.4 and 1.2, respectively.

Figure 5 :Figure 6 :
Figure 5: (a) Simulated time series containing single seismic event generated according to(16) with parameters  = 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).

Figure 7 :
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 ML  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.

Figure 8 :
Figure 8: The comparison of segments produced by different methods applied to the simulated multievents dataset.

Figure 13 :
Figure 13: The comparison of the discussed segmentation methods applied to the real data consisting of a single seismic event.

Figure 14 :Figure 15 :
Figure14: 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 ML  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.
) and12(b).With regard to Figure12(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 ML  .