Time difference of arrival estimation of microseismic signals based on alpha-stable distribution

Microseismic signals are generally considered to follow the Gauss distribution. A comparison of the dynamic characteristics of sample variance and the symmetry of microseismic signals with the signals which follow α-stable distribution reveals that the microseismic signals have obvious pulse characteristics and that the probability density curve of the microseismic signal is approximately symmetric. Thus, the hypothesis that microseismic signals follow the symmetric α-stable distribution is proposed. On the premise of this hypothesis, the characteristic exponent α of the microseismic signals is obtained by utilizing the fractional low-order statistics, and then a new method of time difference of arrival (TDOA) estimation of microseismic signals based on fractional low-order covariance (FLOC) is proposed. Upon applying this method to the TDOA estimation of Ricker wavelet simulation signals and real microseismic signals, experimental results show that the FLOC method, which is based on the assumption of the symmetric α-stable distribution, leads to enhanced spatial resolution of the TDOA estimation relative to the generalized cross correlation (GCC) method, which is based on the assumption of the Gaussian distribution.

1 Introduction

Microseismic monitoring technology has been widely applied to mine rock burst monitoring, oil and gas field fracturing monitoring, reservoir seismic monitoring, slope stability evaluation and so on. Seismic source location is one of the key technologies used (Zhao et al., 2017). The conventional seismic source localization method usually first needs to pick up the P arrival time of multi-channel seismic signals and then calculate the time difference of arrival (TDOA) of the signals to solve the equation to obtain the source location (Schwarz et al., 2016). As a result, the accuracy of the calculated TDOA directly affects the accuracy of the seismic source location. However, in the process of actual operation, the first arrival time of the microseismic signals is not obvious, and there is much external noise (Jia et al., 2015). Therefore, it becomes very difficult to determine the time difference between waves from the same seismic source.

The basic problem the TDOA solves is to measure and estimate the TDOA between waves from the same seismic source accurately and rapidly. Since the classic article on TDOA written by Knapp and Carter was published in 1976, this problem has perpetually been a research focus in the field of international signal processing (Knapp et al., 1976). The common method of TDOA includes the generalized cross correlation method (Knapp et al., 1976; Souden et al., 2010; Jin et al., 2013), the phase spectrum estimation method (Youn et al., 1982; Qiu et al., 2012), the generalized bi-spectral estimation method (Hinich et al., 1992; Hou et al., 2013), the adaptive estimation method (Gedalyahu et al., 2010; Salvati and Canazza, 2013), the energy method based on Hilbert–Huang transform (Sun et al., 2016) and so on. These methods have been widely used in many fields. However, the vast majority of these methods assume that the signals and noises follow the Gaussian distribution. In the case of non-Gaussian distributions, their algorithm shows serious degradation in spatial resolution and does not function anymore (Ma and Nikias, 1996; Cornelis et al., 2010; Park et al., 2011).

Noisy microseismic signals have conspicuous non-stationary characteristics, such as impulsiveness and randomness; therefore, they belong to the category of non-Gaussian signals. If the microseismic signal is simulated by Gaussian signal, it is inevitable that the TDOA algorithm will have serious performance degradation. To solve the problem in theory, we intend to introduce the α-stable distribution to describe the microseismic signal and noise in the distribution model. The α-stable distribution model has achieved excellence in the field of non-Gaussian signal processing, such as seismic inversion, speech de-noising and enhancement, sound source localization and mechanical fault diagnosis (Li and Yu, 2010; Yue et al., 2012; Zhang et al., 2014). However, its use is rare in research projects and publications on TDOA estimation of microseismic signals from the same seismic source. This paper intends to describe the characteristics of microseismic signals and noise with α-stable distributions, studies the impact of non-Gaussian noise on the spatial resolution of TDOA and proposes an improved TDOA algorithm based on fractional low-order covariance (FLOC). Compared with the traditional TDOA algorithm, this improved algorithm could inhibit both the Gaussian noise and the α-stable distribution noise.

2 The TDOA model

2.1 The basic model of TDOA

In order to facilitate data processing, the three component time travel curve of microseismic is first transformed into a set of energy gradient time travel curve (He et al., 2016). The basic model of TDOA is shown in Fig. 1. The original microseismic signal is represented by s(n). It spreads to the two seismic geophones S₁ and S₂ through the rock stratum. Due to different propagation paths, the time at which the signals arrive at the geophones are different.

Figure 1 Two-sensor model of TDOA estimation.

If the microseismic acquisition system is discrete, the signals received by the geophones S₁ and S₂ can be expressed as

$$\begin{array}{}\text{(1)}& \left\{\begin{array}{l}{x}_{\mathrm{1}}\left(n\right)={\mathit{\lambda }}_{\mathrm{1}}s\left(n\right)+b_{\mathrm{1}}\left(n\right)\\ x_{\mathrm{2}}\left(n\right)={\mathit{\lambda }}_{\mathrm{2}}s\left(n-D\right)+b_{\mathrm{2}}\left(n\right)\end{array}\right.,\end{array}$$

In the above equation, s(n) represents the original signal. D is the time delay value. b₁(n) and b₂(n) are the external Gaussian noises. In addition, s(n), b₁(n) and b₂(n) are uncorrelated.

Correlation analysis is commonly used to calculate the TDOA estimation of two signals. In the case that the substance of the problem is not affected and the calculation is simplified, we take ${\mathit{\lambda }}_{\mathrm{1}}={\mathit{\lambda }}_{\mathrm{2}}=\mathrm{1}$. Then, the cross correlation function of the two microseismic signals x₁(n) and x₂(n) can be represented as

$$\begin{array}{ll}\text{(2)}& {\displaystyle }& {\displaystyle }{R}_{x_{\mathrm{1}}{x}_{\mathrm{2}}}\left(\mathit{\tau }\right)=E\left[x_{\mathrm{1}}\left(n\right)x_{\mathrm{2}}\left(n+\mathit{\tau }\right)\right]{\displaystyle }& {\displaystyle }=R_{\text{ss}}\left(\mathit{\tau }-D\right)+R_{s{b}_{\mathrm{1}}}\left(\mathit{\tau }-D\right)+R_{s{b}_{\mathrm{2}}}\left(\mathit{\tau }\right)+R_{b_{\mathrm{1}}{b}_{\mathrm{2}}}\left(\mathit{\tau }\right),\end{array}$$

where R_ss(⋅) represents the auto-correlation function of the original signal. R_pq(⋅) is the cross-correlation function of the two signals p and q. It is assumed that s(n), b₁(n) and b₂(n) are unrelated Gaussian noises. Then

$$\begin{array}{}\text{(3)}& R_{s{b}_{\mathrm{1}}}\left(\mathit{\tau }-D\right)=R_{s{b}_{\mathrm{2}}}\left(\mathit{\tau }\right)=R_{b_{\mathrm{1}}{b}_{\mathrm{2}}}\left(\mathit{\tau }\right)=\mathrm{0}.\end{array}$$

Eq. (2) can be rewritten as

$$\begin{array}{}\text{(4)}& R_{x_{\mathrm{1}}{x}_{\mathrm{2}}}\left(\mathit{\tau }\right)=R_{\text{ss}}\left(\mathit{\tau }-D\right).\end{array}$$

The auto-correlation function has the property

$$\begin{array}{}\text{(5)}& \left|R_{\text{ss}}\left(\mathit{\tau }-D\right)\right|\le R_{\text{ss}}\left(\mathrm{0}\right).\end{array}$$

Therefore, R_ss is maximized when $\mathit{\tau }-D=\mathrm{0}$. Thus, the TDOA estimation between the two seismic geophones can be expressed by the maximum of R_ss(τ−D).

$$\begin{array}{}\text{(6)}& \widehat{D}=\text{arg}\left\{\underset{\mathit{\tau }}{max}\left[R_{x_{\mathrm{1}}{x}_{\mathrm{2}}}\left(\mathit{\tau }-D\right)\right]\right\}\end{array}$$

When noisy signal follows the Gaussian distribution, the above method can estimate the time delay by detecting the peak position of the cross-correlation function of the signals x₁(n) and x₂(n). Based on this observation, the authors propose the generalized cross correlation method (Hertz and Azaria, 1985; Kang et al., 2012; Zhang et al., 2015), the phase spectrum estimation method (Harada, 2014; Choudhuri et al., 2016), the adaptive estimation method (Carrier and Got, 2014; Wang et al., 2017) and so on to implement the TDOA estimation, which significantly improves the anti-noise property, estimation accuracy and resolution of the algorithm. However, these algorithms are all based on the second-order statistics and the assumption that the noises follow the Gaussian distribution. In the process of microseismic monitoring, noisy signals are non-Gaussian and their pulse characteristics are obvious (Xu et al., 2015; Jia et al., 2016). Therefore, Eq. (3) is not suitable for them and it is necessary to introduce other models to describe the distribution characteristics of the microseismic signals and noises, and establish a new method to estimate the time delay.

2.2 The model of α-stable distribution

In the process of microseismic monitoring, external noises are composed of man-made noises, mechanical vibration, etc. The common characteristics of these noises are that their time domain waveforms have conspicuous pulse characteristics, the energy diminishes from low to high frequencies and their corresponding probability density functions have a thicker tail than that of Gaussian signals. In the field of signal processing, this type of non-Gaussian noise is usually described by the α-stable distribution model.

The α-stable distribution is a random signal model that can be applied to an extensive range of problems. Except for a few specific situations, there is no uniform probability density function expression; therefore, Eq. (7) is used to express it (Shao and Nikias, 1993).

$$\begin{array}{}\text{(7)}& {\displaystyle }& {\displaystyle }\mathit{\phi }\left(m\right)=\exp \left\{\mathrm{j}\mathit{\delta }m-\mathit{\gamma }{\left|m\right|}^{\mathit{\alpha }}\left[\mathrm{1}+\mathrm{j}\mathit{\beta }\text{sgn}\left(m\right)\mathit{\omega }\left(m,\mathit{\alpha }\right)\right]\right\},\text{(8)}& {\displaystyle }& {\displaystyle }\mathit{\omega }\left(m,\mathit{\alpha }\right)=\left\{\begin{array}{ll}\tan {\displaystyle \frac{\mathit{\alpha }\mathit{\pi }}{\mathrm{2}}},& \mathit{\alpha }\ne \mathrm{1}\\ {\displaystyle \frac{\mathit{\pi }}{\mathrm{2}}}\log \left|m\right|,& \mathit{\alpha }=\mathrm{1}\end{array}\right.,\text{(9)}& {\displaystyle }& {\displaystyle }sgn\left(m\right)=\left\{\begin{array}{ll}\mathrm{1},& m>\mathrm{0}\\ \mathrm{0},& m=\mathrm{0}\\ -\mathrm{1},& m<\mathrm{0}\end{array}\right.,\end{array}$$

where φ(m) is the characteristic function of the probability density, α represents the characteristic exponent. Smaller values of α result in thicker tails of the probability density function. β is the skew parameter, representing the deviation degree of signals. It is a symmetric α-stable distribution signal when β=0, which is also called the SαS distribution. γ is the scale parameter, representing the dispersion degree of signal around the location parameters, which is similar to the variance in the Gaussian distribution. δ is the location parameter, which is similar to the mean or mid-value in the Gaussian distribution.

We can infer from Eq. (7) that the corresponding eigenfunction is the same as when α=2. That is to say, the Gaussian distribution is a special case of the α-stable distribution. When $\mathrm{0}<\mathit{\alpha }<\mathrm{2}$, Eq. (7) represents the eigenfunction of the signals following the non-Gaussian distribution, which is also called the fractional lower-order α-stable distribution.

2.3 Non-Gaussian property of microseismic signals

The difference in determining a signal between the Gaussian distribution and the α-stable distribution is that the latter has stronger pulse characteristics. Due to the existence of the pulse, the secondary moment of the observation data that follow the α-stable distribution is not convergent, and there is no limited high-order moment above the second order. However, the observation data that follow the Gaussian distribution have both stable secondary moment and limited high-order moment (Sun and Qiu, 2008). Therefore, whether or not the signal follows the Gaussian distribution can be determined from whether or not the sample variance of the observed data is convergent.

If x_i, $i=\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{\dots }N$ represents the observed data sequence and N represents the sample number of observed data. The dynamic sample variance of the first $k\left(\mathrm{1}\le k\le N\right)$ observed data is defined as

$$\begin{array}{}\text{(10)}& {\displaystyle }& {\displaystyle }{S}_k^{\mathrm{2}}={\displaystyle \frac{\mathrm{1}}{k}}\sum _{i=\mathrm{1}}^k{\left(x_i-\stackrel{\mathrm{‾}}{x},\right)}^{\mathrm{2}}\text{(11)}& {\displaystyle }& {\displaystyle }\stackrel{\mathrm{‾}}{x}={\displaystyle \frac{\mathrm{1}}{k}}\sum _{i=\mathrm{1}}^k{x}_i.\end{array}$$

With the continuous increase of k, if $S_k^{\mathrm{2}}$ converges to a certain value, the observed data sequence follows the Gaussian distribution. Otherwise, it follows the α-stable distribution. To illustrate the changes of the dynamic sample variance of the Gaussian signals and the α-stable distribution signals, three sets of random data are produced for comparison. The sample length of the three sets of data are all 1000 (Fig. 2), and Fig. 2a is a Gaussian signal. This means that α=2.0, β=0, γ=1, δ=0. Figure 2b is a random signal that follows the α-stable distribution, and α=1.6, β=0, γ=1, δ=0. Figure 2c is another random signal following the α-stable distribution, and α=1.2, β=0, γ=1, δ=0. The figures (Fig. 2d–e) are the dynamic sample variances corresponding to signals (Fig. 2a–c), respectively.

Figure 2 Gaussian signal, α-stable distribution signal and their dynamic sample variance.

Figure 3 (a) Microseismic wave and (b) its dynamic sample variance.

A comparison of the waveform characteristic of the signals (Fig. 2a–c) shows that with the gradual decrease of the characteristic exponent α, the pulse characteristic of signals is enhanced. The signal (Fig. 2a) follows the Gaussian distribution. Its pulse characteristic is not obvious, and its dynamic sample variance converges to a stable value. The characteristic exponent α of the signal (Fig. 2b) is 1.6. It has a strong pulse characteristic. Its dynamic sample variance springs stepwise and does not converge to a stable value with an increase in sample points. The characteristic exponent α of the signal (Fig. 2c) is 1.2. Its pulse characteristic is more obvious. The step amplitude of the dynamic sample variance increases sharply; therefore, it is more difficult to converge to a stable value.

Figure 4 (a) The α-stable densities with different skew parameter; (b) the α-stable densities of microseismic signals.

We select a measured microseismic wave and calculate its dynamic sample variance according to Eq. (10) (Fig. 3). It shows that the microseismic signal's dynamic sample variance jumps stepwise and also does not converge to a stable value. Thus, one can conclude that the microseismic signal follows the fractional lower-order α-stable distribution. Through the analysis of a large number of seismic signals and the calculation of characteristic exponents, current literature (Yue et al., 2013) shows that the characteristic exponent α of a seismic signal is less than 2, usually between 1.8458 and 1.9301.

2.4 The determination of symmetry property of microseismic signal

Before the parameter estimation of the α-stable distribution, we should determine whether the distribution of the signal is symmetric. The methods for identifying symmetry are listed below:

1. Draw the probability density curve of the sample sequence and observe the symmetry

2. Count the number of positive and negative values in the sample sequence. If the number of positive and negative values are approximately same, the signal is symmetric.

Figure 4a shows that when the skew parameter β=0, the probability density curve is symmetric; when β=0.8, the probability density curve is right-skewed; when $\mathit{\beta }=-\mathrm{0.8}$, the probability density curve is left-skewed. Figure 4b shows five probability density curves of microseismic signals from the same seismic source. It is obvious that these curves are symmetric. As a result, the distribution of microseismic signal is considered symmetric.

For further validation of the symmetry of microseismic signal, we randomly select 30 signals from the microseismic records in different places, truncate the continuous 3000 sampling points of each signal and then count the number of positive and negative values. The absolute value for the difference between the numbers of positive and negative is shown in Fig. 5.

Figure 5 The absolute value for the difference between the positive and negative counts in each microseismic record.

According to the data in Fig. 5, we can use estimate maximum likelihood estimator for parameters μ (difference of data number) and δ (standard error), and μ=1.8667 and δ=26.8356 are obtained. Compared with the 3000 of sample data, the microseismic signal is approximately considered symmetric.

In conclusion, the microseismic signal follows the symmetric α-stable distribution, which is also called the SαS distribution. Because the α-stable distribution does not have limited secondary and high-order moment, the above TDOA method is based on the assumption that the secondary moment (or high-order moment) and the Gaussian noise shows serious performance degradation. It is necessary to do some research on the new TDOA algorithm based on the low-order statistics.

3 The improved TDOA estimation algorithm

3.1 The TDOA estimation based on FLOC

According to the study of Sect. 2.4, the microseismic signals and noises are more consistent with the α-stable distribution. Therefore, this paper is intended to describe the characteristics of microseismic signals and noises by α-stable distribution. ON this basis, an improved TDOA algorithm based on fractional low-order covariance (FLOC) is proposed. Compared with the traditional TDOA algorithm, the improved algorithm has a good suppression effect on the noise of α-stable distribution noise and Gauss noise.

In the case that the noise follows the α-stable distribution, an existing study (Ma and Nikias, 1995) puts forward a TDOA algorithm based on fractional low-order covariance (FLOC). The FLOC of two signals x_i(t) ($i=\mathrm{1},\mathrm{2}$) is defined as

$$\begin{array}{ll}\text{(12)}& {\displaystyle }& {\displaystyle }{R}_{\text{d}}\left(\mathit{\tau }\right)=E\left[x_{\mathrm{2}}{\left(t\right)}^{<A>}{x}_{\mathrm{1}}{\left(t+\mathit{\tau }\right)}^{<B>}\right],{\displaystyle }& {\displaystyle }\mathrm{0}\le A<{\displaystyle \frac{\mathit{\alpha }}{\mathrm{2}}},\mathrm{0}\le B<{\displaystyle \frac{\mathit{\alpha }}{\mathrm{2}}},\mathrm{0}<\mathit{\alpha }\le \mathrm{2},\end{array}$$

and

$$\begin{array}{}\text{(13)}& x^{<c>}={\left|x\right|}^c\text{sgn}\left(x\right),\text{sgn}\left(x\right)=\left\{\begin{array}{ll}\mathrm{1},& x>\mathrm{0}\\ \mathrm{0},& x=\mathrm{0}\\ -\mathrm{1},& x<\mathrm{0}\end{array}\right.,\end{array}$$

where A and B represent the fractional low-order exponents of the two input signals x_i(t) ($i=\mathrm{1},\mathrm{2}$), respectively. τ is the translation relative to the signal x₁(t) when calculating FLOC. The TDOA estimation can be obtained by detecting the peak of the function R_d(τ).

$$\begin{array}{}\text{(14)}& D=-\arg \left\{\underset{\mathit{\tau }}{max}\left[R_{\text{d}}\left(\mathit{\tau }\right)\right]\right\}\end{array}$$

The FLOC algorithm can be used for the TDOA estimation of microseismic signals. If the two microseismic signal samples are x_i(n) ($i=\mathrm{1},\mathrm{2}$; $n=\mathrm{1},\mathrm{2},\mathrm{\dots },N)$, Eq. (12) can be expressed by

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\widehat{R}}_{\text{d}}\left(\mathit{\tau }\right)={\displaystyle \frac{\mathrm{1}}{N}}\sum _{n=\mathrm{1}}^N{\left|x_{\mathrm{2}}\left(n\right)\right|}^A{\left|x_{\mathrm{1}}\left(n+\mathit{\tau }\right)\right|}^B\cdot \text{sgn}\left[x_{\mathrm{2}}\left(n\right)x_{\mathrm{1}}\left(n+\mathit{\tau }\right)\right],\\ \text{(15)}& {\displaystyle }& {\displaystyle }\mathrm{0}\le A<{\displaystyle \frac{\mathit{\alpha }}{\mathrm{2}}},\mathrm{0}\le B<{\displaystyle \frac{\mathit{\alpha }}{\mathrm{2}}},\mathrm{0}<\mathit{\alpha }\le \mathrm{2}.\end{array}$$

The TDOA estimation can be obtained by detecting the peak of the function${\widehat{R}}_{\text{d}}\left(\mathit{\tau }\right)$.

$$\begin{array}{}\text{(16)}& \widehat{D}=-\arg \left\{\underset{\mathit{\tau }}{max}\left[{\widehat{R}}_{\text{d}}\left(\mathit{\tau }\right)\right]\right\}\end{array}$$

The TDOA method based on FLOC requires very few calculations and its real-time implementation is simple. However, the α parameter needs to be estimated in advance; otherwise, the FLOC algorithm will have serious performance degradation and will lead to incorrect results when A and B are greater than α∕2.

3.2 The estimation of the characteristic exponent α

For the random variable X, which follows the α-stable distribution, the fractional lower-order moment is defined as E(|X|^p), $\mathrm{0}<p<\mathit{\alpha }\le \mathrm{2}$. p is the order of fractional lower-order moment. From the Zolotarev theorem (Zolotarev, 1966) we obtain

$$\begin{array}{}\text{(17)}& {\displaystyle }& {\displaystyle }E\left({\left|X\right|}^p\right)=C\left(p,\mathit{\alpha }\right){\mathit{\gamma }}^{p/\mathit{\alpha }},\text{(18)}& {\displaystyle }& {\displaystyle }C\left(p,\mathit{\alpha }\right)={\displaystyle \frac{{\mathrm{2}}^p\mathrm{\Gamma }\left(\frac{p+\mathrm{1}}{\mathrm{2}}\right)\mathrm{\Gamma }\left(\mathrm{1}-\frac{p}{\mathit{\alpha }}\right)}{\sqrt{\mathit{\pi }}\mathrm{\Gamma }\left(\mathrm{1}-\frac{p}{\mathrm{2}}\right)}},\end{array}$$

where α represents the characteristic exponent, γ represents the scale parameter and Γ(⋅) represents the gamma function.

If the random variable X follows the SαS distribution, a study has found that there is a negative-order moment in the SαS distribution (Ma and Nikias,1995). Equation (17) can then be changed to

$$\begin{array}{}\text{(19)}& E\left({\left|X\right|}^p\right)=C\left(p,\mathit{\alpha }\right){\mathit{\gamma }}^{p/\mathit{\alpha }},-\mathrm{1}<p<\mathit{\alpha }\le \mathrm{2},\end{array}$$

because

$$\begin{array}{}\text{(20)}& E\left({\left|X\right|}^p\right)=E\left(e^{p\log \left|X\right|}\right)=C\left(p,\mathit{\alpha }\right){\mathit{\gamma }}^{p/\mathit{\alpha }}.\end{array}$$

Eq. (20) is continuous at the point p=0 after the introduction of negative-order moment. If Y=log |X|, E(e^pY) is the moment-generating function of Y and

$$\begin{array}{}\text{(21)}& E\left(e^{pY}\right)=\sum _{k=\mathrm{0}}^{\mathrm{\infty }}E\left(Y^k\right){\displaystyle \frac{p^k}{k\mathrm{!}}}=C\left(p,\mathit{\alpha }\right){\mathit{\gamma }}^{p/\mathit{\alpha }}.\end{array}$$

Then any order moments of Y are limited and

$$\begin{array}{}\text{(22)}& E\left(Y^k\right)={\displaystyle \frac{d^k}{d{p}^k}}\left[C\left(p,\mathit{\alpha }\right){\mathit{\gamma }}^{p/\mathit{\alpha }}\right]\left|{}_{p=\mathrm{0}}\right..\end{array}$$

This can be simplified to

$$\begin{array}{}\text{(23)}& E\left(Y\right)=C_{\text{e}}\left({\displaystyle \frac{\mathrm{1}}{\mathit{\alpha }}}-\mathrm{1}\right)+{\displaystyle \frac{\mathrm{1}}{\mathit{\alpha }}}\log \mathit{\gamma },\end{array}$$

where C_e=0.57721566… is a Euler constant. Then

$$\begin{array}{}\text{(24)}& \text{Var}\left(Y\right)=E\left\{{\left[Y-E\left(Y\right)\right]}^{\mathrm{2}}\right\}={\displaystyle \frac{{\mathit{\pi }}^{\mathrm{2}}}{\mathrm{6}}}\left({\displaystyle \frac{\mathrm{1}}{{\mathit{\alpha }}^{\mathrm{2}}}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\right).\end{array}$$

For the microseismic signals Y_i ($i=\mathrm{1},\mathrm{2},\mathrm{\dots },N$, N is the sampling number), and the mean and variance can be obtained by Eqs. (25) and (26), respectively.

$$\begin{array}{}\text{(25)}& {\displaystyle }& {\displaystyle }\stackrel{\mathrm{‾}}{Y}={\displaystyle \frac{\mathrm{1}}{N}}\sum _{i=\mathrm{1}}^N{Y}_i\text{(26)}& {\displaystyle }& {\displaystyle }\text{Var}\left(Y\right)={\displaystyle \frac{\mathrm{1}}{N}}\sum _{i=\mathrm{1}}^N{\left(Y_i-\stackrel{\mathrm{‾}}{Y}\right)}^{\mathrm{2}}\end{array}$$

Plugging in the value gained from Eq. (26) into Eq. (24), we can obtain an estimated value of α. Then, the value of α is plugged into Eq. (23), and we obtain the value of γ.

Figure 6 Comparison of the TDOA estimation results of GCC, FLOC and PHAT-GCC.

3.3 Algorithm procedures

If x_i(n) ($i=\mathrm{1},\mathrm{2}$; $n=\mathrm{1},\mathrm{2},\mathrm{\dots },N$) represents the sample of two microseismic signals from the same seismic source, the TDOA algorithm is shown below:

Step 1:

for a given sequence of discrete signal x₁(n) and x₂(n), calculate their characteristic exponents α₁ and α₂ according to Eqs. (24), (25) and (26);

Step 2:

assigning $A=\frac{\mathrm{0.95}\times {\mathit{\alpha }}_{\mathrm{1}}}{\mathrm{2}}$ and $B=\frac{\mathrm{0.95}\times {\mathit{\alpha }}_{\mathrm{2}}}{\mathrm{2}}$ to find that $\mathrm{0}\le A<\frac{{\mathit{\alpha }}_{\mathrm{1}}}{\mathrm{2}}$, $\mathrm{0}\le B<\frac{{\mathit{\alpha }}_{\mathrm{2}}}{\mathrm{2}}$, 0.95 is an empirical value;

Step 3:

add the Hanning window to x₁(n) and x₂(n), and set the window lengths to max(size(x₁(n))) and max(size(x₂(n))). The cross-correlation function ${\widehat{R}}_{\text{d}}\left(\mathit{\tau }\right)$ of x₁(n) and x₂(n) is calculated according to Eq. (15);

Step 4:

detect the peak of the function ${\widehat{R}}_{\text{d}}\left(\mathit{\tau }\right)$. Then, the TDOA estimation $\widehat{D}$ can be obtained.

4 Simulation and analysis

The signals Ricker1 and Ricker2 used in the simulation are two Ricker wavelets. Their spectral peak frequency is 25 Hz. The sampling frequency is 1 kHz, and the number of sampling points is 1000. The time delay between the two Ricker wavelets is set to 70 ms (Fig. 6a). The generalized signal-to-noise ratio (GSNR) is defined in Eq. (27) and used to describe the power ratio of signal and noise (Ma and Nikias, 1996).

$$\begin{array}{}\text{(27)}& \text{GSNR}=\mathrm{10}\lg {\displaystyle \frac{{\mathit{\sigma }}_{\text{s}}}{\mathit{\gamma }}}.\end{array}$$

In the equation, σ_s represents the signal power, and γ represents the noise figure of the α-stable distribution.

4.1 Experiment 1

The spatial resolution on TDOA estimation of the generalized cross correlation (GCC), PHAT-GCC (phase transfer–generalized correlation) method based on the Gaussian distribution and the FLOC method based on the non-Gaussian distribution are compared and verified. α-stable distribution noises to the two Ricker wavelets are added. The parameter of the α-stable distribution (α, β, γ, δ) is set to (1.2, 0, 1, 0). Because of the randomness of α-stable distribution noises, the two noises are independent of each other. The two Ricker wavelets with added noises are shown in Fig. 6b and c. In the case of α-stable distribution noises, the TDOA estimation results of the GCC, PHAT-GCC method and the FLOC method when GSNR = 0 dB and GSNR = 15 dB are shown in Fig. 6d–j.

It is evident from Fig. 6d, f, h and j that the GCC and PHAT-GCC method shows serious performance degradation when GSNR = 0 dB and GSNR = 15 dB. There are several peak positions in the curve so that the correct result is difficult to get. However, the FLOC method has a strong anti-interference ability. The peak appears at 70 ms and can estimate the time delay correctly.

4.2 Experiment 2

The influence of different α to the TDOA estimation results are verified. The two noises are generated randomly when α takes different values between 0 and 2 and are added to Ricker1 and Ricker2, respectively. When GSNR = 0 dB, the TDOA estimation result of the two signals with noises obtained by the FLOC method is shown in Fig. 7. Figure 7a–d show the waveforms of Ricker1 and Ricker2 with different noise signals added to these waveforms. Figure 7a^*–d^* show the corresponding TDOA estimation results.

Figure 7 The influence of different α to the TDOA estimation result.

It is evident from Fig. 7 that a smaller α (which implies that there is a stronger pulse of noise) corresponds to better performance in the TDOA estimation of the FLOC algorithm. When α is close or equal to 2, the performance is degraded (Fig. 7d^*), but it still can obtain the correct TDOA estimation result: 70 ms. Therefore, the FLOC method performs well irrespective of the noise and follows the Gaussian distribution or the α-stable distribution.

5 Case study

To verify the effectiveness of the FLOC method for TDOA estimation of real microseismic signals, we select eight microseismic signals from the same seismic source to do the experiment. The eight signals come from the ISS microseismic monitoring system of a coal mine in central China. Seismic geophones are laid along the mining roadway every 50 m in the system. The frequency bandwidth of the seismic geophones is between 3 and 2000 Hz. The data acquisition frequency is 1 kHz. For convenience of comparing and analyzing the experiment results, the first 2000 sampling points of each waveform are picked as the data object. The P arrival time of each microseismic signal is recorded manually, and the time delay between any two of the microseismic signals as a reference of the experimental result is calculated.

As an example, the microseismic signals in roadway nos. 2 and 7 are selected to explain the result. The waveforms of the microseismic signals after interception are shown in Fig. 8a and b. The time delay between the two signals obtained by manual method is 19 ms. First, when the microseismic signal follows the Gaussian distribution, the PHAT-GCC method, which, of the GCC methods, performs best, is chosen for the TDOA estimation of microseismic signals from the same seismic source. The result is shown in Fig. 8c. Second, when the microseismic signal follows the α-stable distribution, FLOC is used for the TDOA estimation. The characteristic exponent α of the two picked signals are calculated according to Eqs. (23), (24) and (25). α₂=1.802, α₇=1.835 are obtained as results. According to Step 2 in Sect. 3.3, assign A=0.85595 and B=0.87163. The TDOA estimation result obtained by the FLOC method is shown in Fig. 8d.

Figure 8 The comparison of TDOA estimation results based on PHAT-GCC and FLOC.

Figure 8c and d show that the two methods both obtain the correct result, 19 ms, but the peak of the FLOC method is sharper than the GCC method. This implies that the FLOC method performs better.

Each of the eight microseismic signals is considered to be a set of data following the α-stable distribution. Their characteristic exponent α values are calculated according to Eqs. (23), (24) and (25) and are shown in the table (Table 1). The values are between 1.802 and 1.913 (Fig. 9). It can be seen that the characteristic exponent α values of all of the signals are less than 2. According to the data in Table 1, we can use estimate maximum likelihood estimator for parameters μ (difference of data number) and δ (standard error), and μ=1.8550, δ=0.0377 can be obtained.

Table 1 The characteristic exponent α of microseismic signal.

We can obtain 28 pairs of microseismic signals by the pair combination of the eight signals in Table 1. The comparison of TDOA estimations obtained by the PHAT-GCC, FLOC and manual method is shown in the table (Table 2).

Table 2 The comparison of TDOA estimation results of microseismic signal.

Figure 9 Comparison of the characteristic exponent α of eight microseismic signals.

An analysis of tables (Tables 1, 2) indicates that the pulse of actual microseismic signal is stronger than the one following the Gaussian distribution. Because the characteristic exponent of the actual microseismic signal is less than 2, it is considered to be a signal following the α-stable distribution. Based on this observation, we can say that the spatial resolution of the FLOC method is better than the PHAT-GCC method at TDOA estimation of microseismic signals.

6 Conclusions

Through the analysis of the convergence of dynamic sample variance, the microseismic signal with noises is shown to follow the α-stable distribution. The analysis of the symmetry of probability density curve of the sample sequence proves that the microseismic signal is approximately symmetric. Therefore, it is more reasonable to regard the microseismic signal with noises as the α-stable distribution signal.

Because of the absence of second-order statistics of α-stable distribution, one cannot obtain optimal or correct estimation values via the traditional TDOA method based on the Gaussian distribution.

Microseismic monitoring data obtained from a coal mine in central China are used for TDOA estimation based on the GCC method and the FLOC method to study cases when the microseismic signals follow the Gaussian distribution and the α-stable distribution. In the comparison of the results and the time delay obtained manually, we observe that the FLOC method performs better than the traditional GCC method irrespective of whether the noise follows the Gaussian distribution or the α-stable distribution. This method is suitable for the TDOA estimation of microseismic signals from the same seismic source.