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 1Two-sensor model of TDOA estimation.

Download

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

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

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

Eq. (2) can be rewritten as

The auto-correlation function has the property

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

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

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

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 2Gaussian signal, α-stable distribution signal and their dynamic sample variance.

Download

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

Download

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.

Download

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 5The absolute value for the difference between the positive and negative counts in each microseismic record.

Download

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

and

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(τ).

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

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

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

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

because

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

Then any order moments of Y are limited and

This can be simplified to

where C_e=0.57721566… is a Euler constant. Then

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.

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 6Comparison of the TDOA estimation results of GCC, FLOC and PHAT-GCC.

Download

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.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 7The influence of different α to the TDOA estimation result.

Download

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.