1 INTRODUCTION

One of the keys to monitoring and early warning of underground coal mine dynamic disaster through microseismic monitoring system is precision location. The seismic signals from the fracture of the underground rock are received by the micro-seismic sensors. The seismic source position inversion can be realized by combining the data of sensor coordinates and wave velocity of the rock. The arrangement of the sensor is restricted by the narrow space of the tunnel in the coal mine, therefore the error of the location inversion will be greatly increased.

As for micro-seismic monitoring the working face of underground coal mine, the sensors can only be installed around the crossheading (Ye et al., 2008). To avoid the sensors being arranged in the same horizontal plane, it needs to drill deep hole for installing and increasing the vertical distance between the points. However, issues including the collapse of the deep holes and the influence of the horizontal stress in the rock mass will affect the installation. Poor underground construction condition and deep hole installation means the increased construction difficulty and cost, meanwhile the drilling depth is limited. The surface safety monitoring of the mined-out area and the adjacent micro-seismic monitoring are also facing the same problems. Therefore, the problems of arrangement optimization of micro-seismic points and high precision location in ill-posed condition of micro-seismic measuring points need further research.

There are more than ten methods for source location and the location precision is between 5 m to 50 m (Lin et al., 2010; Dong et al., 2011). The linear location is one of the important methods. The source coordinates, the measuring points coordinates and the first arrival time can constitute the quadratic nonlinear equations. The nonlinear equations can be transformed into linear equations by reducing the dimension (Li and Zhu, 2009). Gaussian elimination method can be used to solve the equations (Li, 2008). Nevertheless, the solution error will be large when the condition number of linear equations coefficient matrix is large. When the order of magnitude of coefficient matrix elements has great difference, row or column equilibrium can be used to improve the matrix condition number (Liu, 1999). Linear location contains Geiger method, combined location method, relative location method, double residual method and etc. (Tian and Chen, 2002; Jin et al., 2014), nevertheless, there is no research on the ill-posed problems caused by measuring points arrangement. The error of location will be greatly increased when the seismic signals are only received by limited sensors.

Micro-seismic source location is an inversion problem based on the data of the first arrival time, wave velocity, etc. In 1960s, Tikhonov and Arsenin (1977) and Phillips independently proposed the regularization method to solve the mathematic inverse problems and obtain stable approximate solution. Regularization methods include the truncated singular value decomposition (SVD), improved truncated SVD, generalized truncated SVD, the attenuation of the SVD, the minimum entropy principle, the maximum entropy principle, quadratic constrained least square method and truncated total least square (Hansen, 2007). At present, Tikhonov method is widely used. Regularization parameter is one of the most important parameters for the regularization theory and the choice of the regularization parameter is the basis of regularization method. A good regularization parameter λ can generate maximum equilibrium between right-hand side disturbance error and the regularization solution error (Wang, 2007). The λ controls the weight which minimizes the solution norm and the residual norm, balancing the instability and smooth. The selection of regularization parameter is divided into two categories. One is the error estimation which is based on the principle of deviation (Cheng, 2010). The main problem of this method will be an over regularized solution with huge error when an oversize parameter is chosen. The other is extracting necessary information from the right-hand side, including the L-curve criterion, the generalized cross validation (GCV) and approximate optimal criterion. GCV method may mistake correlated noise for effective signal. L-curve can distinguish the correlation interference excellently well (Hansen and O'Leary., 1993; Hansen, 2001).

Pang (2004) found that Tikhonov regularized solution converges much faster than least square and regularized solution when using Tikhonov regularization method to iterative calculation for the ill-posed problems in the process of the micro-seismic linear equations location. Zeng et al. (2015) used the relationship between regularization parameter λ and cut-off wave number ωc to obtain the improved regular low pass filter function and apply to the conversion of the vertical derivative of the potential field successfully.

The most convenient graphical analysis tool for discrete ill-posed problems is the L-curve which includes all valid regularization parameters. L-curve is a compromise method to minimize the solution norm and residual error norm, and shows the dependency of these quantities on the regularization parameter (Hansen 1997). This method is firstly applied by Miller (1970), Lawson and Hanson (1974) to deal with the least square problems. L-curve plays an essential role in the connection of Tikhonov regularization and ill-posed problems. Only if the limited sensors receive the vibration signals and the monitoring points are approximately on the same plane, the matrix structured by monitoring points coordinates, arrival time and etc. may be ill-posed. Tikhonov regularization method based on L-curve criterion can improve the solution precision of the ill-posed matrix.

This paper focuses on the existing problems of the sensors arranged around underground roadway and approximately on the same plane in micro-seismic monitoring of confined space, analyses the factors which influence the source location precision, and puts forward pretreatment before solving regularized solution for ill-posed problems on the basis of measuring points arrangement optimization and the characteristics of field test data. In order to reduce the gap of magnitude order between coefficient matrix A and b, centralization algorithm for measuring points original coordinate is put forward. The row balance pretreatment operator matrix is constructed to reduce the condition number of coefficient matrix A. Firstly, the monitoring data is carried out to SVD after centralization and row balance pretreatment. Then condition number is calculated by L-curve criterion. Finally the applicability of the Tikhonov regularized solution is validated.

2 THE THEORETICAL BASIS

The unknown source coordinates are assumed as (x, y, z), the coordinates of the sensor i are x_i, y_i, z_i. The distance between source and the sensor i is calculated as follows:

(1)

where K_i=x_i² + y_i² + z_i².

The precise coordinates of seismic source can be calculated by turning quadratic nonlinear equation into linear equation. r_{i, 1} is assumed to be the difference between r_i and r₁, and r₁ is the distance between source and the S₁ sensor that first received signal (Li and Zhu, 2009).

(2)

where x_{i, 1}=x_i -x₁, y_{i, 1}=y_i -y₁, z_{i, 1}=z_i -z₁, r₁=(t₁ -t₀)v, r_{i, 1}=(t_i -t₀)v -(t₁ -t₀)v=(t_i -t₁)v. As shown in Fig. 1, t₁ is the time when seismic wave arrives to S₁, t_i is the time when seismic wave arrives to S_i, t₀ is the occurrence time of micro-seismic event and is an unknown quantity, r₁ is the function of t₀, v is the wave velocity.

(3)

Formula (3) can be simplified as formula (4).

(4)

For linear system of equations, A is non-singular matrix and its condition number is calculated as (Li, 2008):

(5)

The linear system of equations is ill-posed when the condition number of A is much greater than 1, cond (A) ≫ 1. The greater the A's condition number is the more ill-posed the linear equation group will be. The solution deviation will be greater by the conventional method (Li, 2008; Xiao, 2003). This study applies formula cond (A)_∞ to calculate condition number of linear equation group coefficient matrix.

3 THE SOLUTION FOR ILL-POSED PROBLEMS

Measuring points arrangement is the key factor that whether the coefficient matrix A is ill-posed or not. Firstly, the measuring points arrangement should be optimized so as to realize high precision monitoring of rock mass activity of the target area. And then the methods of pretreatment and regularization are adopted to solve the still ill-posed problems after optimizing.

3.1 Linear Independent Optimization

To achieve three-dimension location of the seismic source, more than four measuring points shouldn't be installed on the same horizontal plane, otherwise the z coordinates of these points will be equal. A will be a singular matrix detA=0 and A^-1 won't exist, which would lead to non-solution or infinite solution of the equation system. Therefore, detA > 0 should be ensured.

Approximate correlation of the row or column vector of the coefficient matrix A, which means the angle between the vectors is very small, is the direct cause of the ill-posed matrix A. The measuring points should be set at unequal interval in three directions of space and arranged randomly so as to build non-singular coefficient matrix A. Ultimately linear system of equations can be solved by direct method. The formulas a_i ≠ qa_j and c_i ≠ qc_j should be satisfied, where a_i and a_j are row vector, c_i and c_j are column vector, and q is arbitrary real number. The condition number tends to be positive infinity when two arbitrary rows or columns of the coefficient matrix A are proportional and tend to be overlapped completely.

The condition number of matrix A can be reduced effectively by increasing the installation depth difference among sensors for monitoring a wide enough area under the condition that the sensors plane coordinates (x_i, y_i) are ascertained.

3.2 Initial Coordinates Centralization

The equations are ill-posed if a small change of the constant term b causes great change of the solution of the equation AX=b. Supposing b is perturbed by value δ_b, the perturbation of corresponding solution can be calculated as (Liu, 1999):

(6)

(7)

For a certain micro-seismic event, r_{i, 1}, x_{i, 1}, y_{i, 1} and z_{i, 1} are fixed values which are determined by the relative position relation of the measuring points, and aren't affected by the points original coordinates value. The right-hand side of formula (4) can be calculated as follows:

(8)

It can be seen from the above formula that the initial coordinate value of x₁, y₁, z₁ has a huge impact on b. If 8 bit integer geodetic coordinates are adopted to indicate x₁ and y₁, it would cause the huge gap of magnitude order between matrix A and b. Therefore, it needs to transform b to constant term, like A, which wouldn't be affected by initial coordinates, and ensure the solution error will be decreased during regularization solving.

In this study, the following centralization method for the given initial coordinates is put forward.

(9)

(10)

where i=1, …, n. The linear equations are related to the five point coordinates, and therefore only the five x_i coordinates participating the calculation are taken into arithmetic mean. The mean value is zero and the variance matrix is invariant after transformation. The same centralized method is used to treat y_i and z_i, and calculate the right-hand side b with the transformed x_i', y_i', z_i'.

3.3 Row Balance Pretreatment

The difference among x_{i, 1}, y_{i, 1} and z_{i, 1} in magnitude order is huge because of the wide monitoring area and simple rapid installation to save cost. Row balance pretreatment to the formula (4) by the following method is proposed so as to reduce the condition number of the coefficient matrix. Non-singular matrix P ∈ R^n×n is selected and the formula (4) is changed into the equivalent equations:

(11)

(12)

(13)

Formula (12) is put forward by Liu (1999), denoted as m_max for short. In this study, the diagonal matrix construction method based on 2-norm is put forward, denoted as m₂ as follows:

(14)

where ai is the row vector of matrix A, and A=[a₁, a₂, …, a_n]^T, ‖PA‖₂=1. The best construction method would be chosen at the following data processing part by comparing the practical application effect of the two methods.

3.4 Regularization Inversion

A in formula (4) is a bounded linear and compact operator in Euclidean space. The equation has the unique solution when detA ≠ 0. When the condition number is too big, the problem will be unsuitable or ill-posed. To solve this kind of problem, it needs to seek a stable approximate solution in the case of little change in b. The regularized solution x_λ (Tikhonov and Arsenin, 1977) is defined as

(15)

where λ is the regularization parameter which controls the weight of relative boundary constraint and the residual norm minimization. L is the boundary limitation matrix. x^* is the estimated initial value. Tikhonov method is a direct method, the regularization model is to solve the following least square problem (ELDEN 1977).

(16)

The SVD of A is calculated as follows:

(17)

where u_i is row vector in U, v_i is row vector in V and i is row index.

The standard form of the Tikhonov regularization model is L=I_n, x^*=0 and the regularized solution is calculated by the following formula (Hansen 2001):

(18)

(19)

3.5 To Calculate λ by L-curve Criterion

‖x_λ‖₂ and ‖Ax_λ -b‖₂ play an important role in the judgment of λ. The residual error ‖Ax_λ -b‖₂ will be too large induced by over-regularization of the solution. However, it will be much smaller when there is less regularization of the solution. Nevertheless, the regularized solution will be controlled by data error, which will lead ‖x_λ‖₂ to be great. L-curve is a compromise curve when two quantities need to be controlled.

Take a set of constants as λ, and substitute them into formula (18) and (19) to solve x_λ and b -Ax_λ. L-curve constitutes x coordinate ‖Axλ -b‖₂ and y coordinate ‖x_λ‖₂ as shown in Fig. 2.

L-curve criterion means the regularization parameter corresponding to the maximum curvature at the L corner is the optimal regularization parameter. The exact optimal regularization parameter can't be obtained directly from the graph. The optimal regularization parameter in Fig. 2 is 0.039833 and its exact numerical solving formula is defined as (Hansen and O'Leary, 2001):

(20)

(21)

where k is the curvature of L-curve. k is regarded as the vertical coordinate and λ is regarded as the horizontal coordinate in Fig. 2. The λ corresponding to curvature maximum value is the optimal parameter.

3.6 Processing Procedure of Ill-posed Problems

To avoid ill-posed problems from the beginning, the linear independent optimization should be done in the stage of measuring point arrangement. The available points will be decreased as the working face advancing. Pretreatment and regularization are necessary when the coefficient matrix composed of the remaining points or the optimized points inevitably become ill-posed problems. The specific data processing steps are as follows:

(1) Select the minimum first arrival time point as the base reference, and select another 4 points to form Eq.(4).

(2) To calculate the coefficient matrix A, and solve cond (A)_∞. If cond (A)_∞ < 10, then skip to step 4. If cond (A)_∞≥10, skip to step 3.

(3) Implement centralization and row balance pretreatment.

(4) Implement regularization treatment, and obtain the regularized approximate solution.

4 FIELD TESTING 4.1 Measuring Points Arrangement

A certain mine has extremely thick coal seam and its mining method is fully mechanized top caving. For the sake of timely monitoring rock burst to predict mining dynamic hazard, the micro-seismic monitoring system is installed, with eleven triaxial sensors installed in roof and the floor of the track transportation tunnel, as shown in Fig. 3. The black filled dots represent the sensors installed in roof and the squares represent sensors installed in floor. The calibration blasting hole is arranged in the floor.

In Table 1, S₁, S₂, …, S₁₁ are the sensors or measuring points code number, x_i, y_i, z_i are the ith measuring point three-dimensional coordinates, and t_i is the first arrival time of P wave. The three-dimensional coordinates of calibration blasting are 3899.89, 5097.48 and 1009.30. S₁ is the nearest base point to the source. The two sensors S_i and S_j are selected randomly, and the distances between the two sensors to calibration blasting are d_i and d_j. The P wave first arrival time is t_i and t_j respectively. Therefore, the wave velocity can be calculated by

(22)

Four groups of floor measuring points combinations and four of top measuring points combinations are selected. The largest interval, the measuring points of crossing the track transportation tunnel, and the smallest interval are gained from the above groups, which can calculate the corresponding velocity v_ij. The arithmetic mean velocity of all v_ij is 4.35 m·ms^-1, as shown in Table 2.

4.2 The Solution of Heavily Ill-posed Problems

Through designing different inclination angle and different drill hole depth, the three-dimensional direction and differential optimal arrangement is achieved. Unfortunately, it is still limited by the narrow space of the roadway. The maximum distance among z coordinates in bottom is only 1.64 m, and the maximum distance among z coordinates in the roof points is 10.56 m, whereas the maximum distance interval value of y direction measuring points is 174 m. Analysis should be made on the ill-posed problems caused by space layout. The monitoring system layout hasn't gone through linear independent optimization, therefore S₁, S₃, S₇, S₁₀ are almost in the same horizontal plane. When the floor points participate in the seismic source inversion location, there appear severe ill-posed problems in the data of the following 4 groups. Which include combination one (S₁, S₃, S₄, S₇, S₁₁), combination two (S₁, S₃, S₉, S₁₀, S₁₁), combination three (S₁, S₃, S₆, S₇, S₁₀) and combination four (S₁, S₄, S₇, S₉, S₁₁). S₁ is the nearest base point to source. The equation solutions are shown in Table 3. Each combination adopts separately the diagonal matrix of m₂, m_max structure to balance, and the condition numbers are calculated, as shown in Fig. 3. The calculating methods include Gaussian elimination (Gaussian for short), regularization, regularization after centralization (Centralization for short) and regularization after centralization and row balance (Pretreatment for short). The solution error is shown in Fig. 5.

The three kinds of condition numbers of all combinations are basically decreasing as shown in Fig. 4. The condition numbers are largest before row balance pretreatment, and significantly reduced while being calculated with the diagonal matrix constructed by m₂. The solution errors calculated by four algorithms are also decreasing, as shown in Fig. 5. The error of Gaussian elimination solution is the biggest, and the source inversion effect by pretreatment and regularization is the best.

In the four combinations, the condition numbers are large without pretreatment, of which the minimum value is 6062.69 and the maximum is 14834.2. Consequently, all of the equations are ill-posed. The minimum of source location error is 2411.34 m and the maximum is 10558.34 m when solved by Gaussian elimination method, and the results aren't significant. The solution error has been improved when Tikhonov regularization algorithm is adopted without pretreatment. The minimum of source location error is 2370.15 m and the maximum is 2514.29 m, nevertheless the error is still great and the data has no applicable value.

By analyzing the characteristics of the initial coefficient matrix A and b, it can be found that the absolute value of the element in A is 10² and the absolute value of the element in b is 10⁶. The initial analysis is that large difference of the order of magnitude causes the great solution error. The four groups data is processed by the centralization pretreatment, and the matrix A has no change therefore the condition numbers haven't decreased. The order of magnitude of righthand side b decreases significantly, and the absolute value of maximum order of magnitude of element is 10⁴, and most of them are 10³. If the data after centralization is calculated with Gaussian, the location error can't be improved. However, the location result is greatly improved when Tikhonov regularization method is adopted, and the error of combination one and four is reduced to less than 10 m, the minimum location error is only 7.28 m and the average value is 14.4 m. The results have reached the level of practical application.

After centralization pretreatment, row balance is also used to analyze the data of the four groups, it turns out that cond (A)_∞ is improved significantly, and the maximum value is down from 14834.2 to 10818. As is shown in Fig. 4, the condition numbers are calculated with m₂ and are smaller than that calculated with mmax in combination one, two and four. Therefore, diagonal matrix constructed with 2-norm, which is put forward in this study to execute row balance pretreatment, has a better effect. Nevertheless, the equations are still ill-posed and should be regularized. The regularized solution after pretreatment is improved significantly in combination one and two. Combination one has the smallest error of 3.09 m, and the error decreases by 0.1 m in combination three. Only in combination four, the error increases by 1.06 m. The location precision is further improved compared to regularization after centralization, and the average error is decreased to 11.6 m. The measuring points arrangement, signal quality, average velocity and etc will influence the solution error.

5 CONCLUSIONS

Linear independent optimization is needed in measuring point arrangement when linear method is adopted to invert for seismic source. Four or more than four measuring points shouldn't be set on the same horizontal plane, and the measuring point should be set at unequal interval in three dimensional direction and be arranged randomly. The installation depth difference among sensors should be increased.

The initial value of the measuring points is the most important factor that affects the error level of the source location. The centralization pretreatment is proposed based on the characteristics of the micro-seismic data. The order of magnitude of constant term b decreased greatly while coefficient matrix A isn't affected, and the location error is reduced considerably. The results don't have any significance before centralization, whereas the location results reach the level of practical application after centralization.

The 2-norm of the matrix A's row vector is proposed to construct the diagonal matrix P and impose row balance pretreatment. The condition number of matrix A is decreased effectively, nevertheless, the improvement is limited. It is necessary to implement regularization processing when the equations are heavily ill-posed and can't reach the optimization expectation.

The location error is still large, in case the micro-seismic monitoring data solutions are solved directly. The method of pretreating the data before regularization can greatly improve the solution precision. The minimum location error is 3.09 m. There are three groups of source location inversion error that is less than 10 m. The result reaches the requirement of actual engineering application.

In general, the basis to obtain high precision source solution is optimal arrangement of measuring points, high quality sampling signal, velocity model optimization, etc. In this paper, the available and useful information is extracted from the invalid data that is unable to be processed by conventional method. In the case that the valid signal is monitored by limited points, the accurate source location can also be obtained. This paper provides a new algorithm to solve the ill-posed problems of micro-seismic source location, which promotes the development of the micro-seismic monitoring technology.

ACKNOWLEDGMENTS

The research presented in this paper is supported by Natural Science Foundation of China (51574250), the national key research and development program (2016YFC0801405), and the Science and Technology Development Plan of Shandong Province (2014GSF120017).