J. Cent. South Univ. (2012) 19: 2297-2306

DOI: 10.1007/s11771-012-1275-z

Passive seismic velocity tomography on longwall mining panel based on simultaneous iterative reconstructive technique (SIRT)

N. Hosseini¹, K. Oraee², K. Shahriar³, K. Goshtasbi⁴

1. Department of Mining Engineering, South Tehran Branch, Islamic Azad University, Tehran 17776-13651, Iran;

2. Department of Management, University of Stirling, Stirling FK9 4LA, UK;

3. Department of Mining and Metallurgical Engineering, Amirkabir University of Technology, Tehran 15916-34311, Iran;

4. Department of Mining Engineering, Tarbiat Modares University, Tehran 14115-111, Iran

? Central South University Press and Springer-Verlag Berlin Heidelberg 2012

Abstract: Mining operation, especially underground coal mining, always has the remarkable risks of ground control. Passive seismic velocity tomography based on simultaneous iterative reconstructive technique (SIRT) inversion is used to deduce the stress redistribution around the longwall mining panel. The mining-induced microseismic events were recorded by mounting an array of receivers on the surface, above the active panel. After processing and filtering the seismic data, the three-dimensional tomography images of the p-wave velocity variations by SIRT passive seismic velocity tomography were provided. To display the velocity changes on coal seam level and subsequently to infer the stress redistribution, these three-dimensional tomograms into the coal seam level were sliced. In addition, the boundary element method (BEM) was used to simulate the stress redistribution. The results show that the inferred stresses from the passive seismic tomograms are conformed to numerical models and theoretical concept of the stress redistribution around the longwall panel. In velocity tomograms, the main zones of the stress redistribution around the panel, including front and side abutment pressures, and gob stress are obvious and also the movement of stress zones along the face advancement is evident. Moreover, the effect of the advance rate of the face on the stress redistribution is demonstrated in tomography images. The research result proves that the SIRT passive seismic velocity tomography has an ultimate potential for monitoring the changes of stress redistribution around the longwall mining panel continuously and subsequently to improve safety of mining operations.

Key words: longwall mining; passive seismic velocity tomography; simultaneous iterative reconstructive technique (SIRT); boundary element method; stress redistribution; ground control

1 Introduction

Although with development of the mechanized longwall mining method, underground coal mining has been improved both from production and productivity point of views, there are still certain risks in mining that can result in unacceptable level of safety. Generally, one of the most dangerous issues is related to roof fall and ground control. Accurate determination of the stress in surrounding rock mass of underground mining area is a complicated task and mostly depends on both the strength properties of the intact rock and the structural conditions of the rock mass [1]. However, determination of the nature of stress redistribution in the various phases of mining operation is necessary to understand the failure mechanism of the roof strata and the surrounding rock mass and also of the roof caving process. Accurate estimation of the front and the side abutment pressures is also important in designing the support system for longwall panel and side entries [2]. Achieving the optimum advance rate in a longwall face also, to a large extent, depends on the state of abutment pressures [3].

Several methods for estimation of stress around the longwall mining panel have been introduced. For this purpose, field measurement methods are time consuming, difficult and cause disruption in production. Theoretical methods, on the other hand, may be questionable in the accuracy of the obtained results.

In this work, the passive seismic velocity tomography based on simultaneous iterative reconstructive technique (SIRT) inversion, was applied to study the state of stress around the longwall mining panel. In passive tomography, the mining-induced microseismic events will be used as a seismic source. Therefore, the long-term and continuous monitoring of the stress changes will be possible.

Usually, during the pre-failure regime, the p-wave velocity increases linearly with stress at lower stress levels, and then plateaues at higher stress levels. Because the cracks and pore spaces are closures with increasing the stress, the p-wave velocity is related to stress [4].

2 Theoretical concept of stress redistribution around longwall mining panel

To analyze the stress redistribution in the longwall structure and ground strata, a simplified model based on composition of stress dynamic and displacement mechanics was used. The loads that are supported by coal seam prior to extraction, have to be transferred to solid adjacent areas, and based on the properties of mining area, the abutment pressure will be created [3]. In Fig. 1, the stress redistribution around the longwall panel, especially front and side abutment pressures, is shown. In this figure, the states of stress redistribution in sections a-a, b-b, and c-c are also demonstrated. Moreover, the overlaps of side and front abutment pressures are displayed on both side of the face. Generally, the situation of abutment pressure is different; particularly the location of maximum pressure in various longwall coal mines that cause the state of stress in rock mass surrounding panel and also geological features of each coalfield [3].

In mining operation process, the abutment pressures shift in the direction of the face advance. The representative model of traveling the abutment pressure is shown in Fig. 2. The magnitude of the peak stress and its location within the rib are given in Table 1.

The abutment pressure increase in the longwall face is dependent on the geometry of the gob area. Peak abutment pressures increase as the shape of the mined area changes from the slim rectangular opening at the start of mining to a square shape with the advance equal to the longwall face length [2].

It should be mentioned that the magnitude of the peak front abutment pressure is related to the geological properties of immediate roof and main roof, and usually will occur at corners or center of the face [5]. Of course, the mine design and mining operation, especially the advance rate of the face have a significant effect on the location of maximum front abutment pressure [1-2]. The side abutment pressure is due to the portion of roof over the void that is not supported by the caved material or gob [3].

3 Tomography

RADON was the first to establish the tomography framework [6]. He proved that an infinite number of rays passing through a two-dimensional object at an infinite number of angles could generate a perfectly tomography image of the object. This theory was also applied for three-dimensional objects [6]. Tomography includes body divided into the grid cells called pixel in two- dimensional mode, or cubes called voxel in three- dimensional mode to estimate the characteristics of the body in all pixels or voxels [6-7].

Fig. 1 General form of stress redistribution around longwall mining panel [3]

Fig. 2 Concept of abutment pressures movement [3]

Table 1 Abutment pressure distribution

Tomography based on the type of source used, is classified as “active” and “passive” [8]. In active tomography, the seismic wave is artificially created. For example, hammer strikes against roof and rib bolts, and controlled explosions may be used as active source in underground mining tomography. However, in passive tomography, the mining-induced seismic events, such as vibration of drilling operation or coal cutter machine, are used as a source. Usually, using the passive source in mining tomography is preferred. Because passive source really is a part of mining operations, there is no disruption in the mine production, and it also will not impose additional costs. However, in use of the passive source, uncontrolled features of wave including the location, occurrence time, and magnitude are the most challenges.

The velocity tomography relies on a simple relationship; the velocity along a seismic wave is equal to the raypath distance divided by the time to travel between the source and the receiver. Therefore, the time is calculated by integral of the inverse velocity (slowness), from the source to the receiver [9], as shown in Eqs. (1)-(3):

                                                                               (1)

                                                                    (2)

                                                        (3)

where v is the velocity; l is the distance between the source and the receiver; t is the travelling time; p is the slowness; N is the total number of rays; M is the number of voxels; t_i is the travel time of the i-th ray; p_j is the slowness of the j-th voxel; l_ij is the distance of the i-th ray in the j-th voxel.

In passive seismic velocity tomography, the source location and also the raypath must be calculated. For this purpose, the initial velocity model that is developed based on field measured data, is used. Although velocity, distance, and travel time in each individual voxel (or pixel) are unknown, by the initial velocity model, distance and time along the raypath can be calculated. Of course, calculation of the distance in each voxel (or pixel) is not difficult, while the calculation of the velocity and the time are complicated [9]. If the time, distance, and slowness for each voxel are arranged into the matrix form, the velocity based on inverse theory is calculated by

                                                                 (4)

where T is the travel time per ray matrix (1×N); L is the distance per ray per voxel matrix (N×M); P is the slowness per grid cell matrix (1×M).

Usually, the inverse problems are either underdetermined (more voxel than rays), or overdetermined (more rays than voxels). The most effective way to solve such problems is iterative process [10]. SIRT is one of the famous methods of iterative process and includes the following steps [11]:

Step 1: Ray tracing.

Step 2: Calculation of the ray distance in voxels that ray passes through.

Step 3: Calculation of the residual time for ray (observed time minus calculated time) based on slowness distribution and save the results.

Step 4: Repeat steps 1 to 3 for all rays.

Step 5: Modification of the slowness in each cell with regards to all the passing rays.

Repeat steps 1-5 until the residual time is less than an acceptable amount; sometimes the solution based on the number of iterative will be limited. Slowness in each cell is modified by

                    (5)

Equation (6) shows the matrix form of Eq. (5):

                        (6)

where p^ák+1? and p^ák? are the slowness vectors (M×1) after the (k+1)-th and the k-th iteration, respectively; ψ_ij=

 is the diagonal matrix (M×M);  is the transpose matrix (N×M);  is the

diagonal matrix (N×N); e^ák? is the residual vector (N×1) after k-th iterative; N is the total number of rays (number of equations), M is the number of voxels (number of parameters),  is the slowness of the j-th voxel after the k-th iterative, and  is the residual of the i-th ray after the k-th iterative [12].

Resolution is one of the most important parameters in evaluating the tomograms. Resolution is typically dependent on three factors: the distribution of the velocity, the source-receiver geometry, and the wavelength [13]. Higher velocity contrasts will cause more ray bending and more uneven ray coverage, which can adversely affect resolution. Also, source-receiver geometry will affect ray coverage. Some studies indicate that resolution is related to wavelength with theoretical resolution about one wavelength [14], although this relationship is not very well known [15].

The reliability of a tomogram can be described by resolution matrix. MENKE [10] gave a simple description of resolution matrix calculation. First, the inverse problem is defined as

Gu=d                                  (7)

where G is the data kernel (ray distance per cell); u is a matrix of model parameters (slowness), and d is a matrix of data (travel times).

The solution to this problem is

u_est=G_-gd                                      (8)

where u_est is the estimate of model parameters, and G_-g is the inverse of G.

It follows that:

u_est=G_-gd_obs=G_-g[Gu_true]=[G_-gG]u_true=Ru_true                     (9)

where d_obs is the observed data; u_true is the true, but unknown set of model parameters; R is the resolution matrix [10].

Equation (9) illustrates that resolution relates the estimated model parameters with the true model parameters. In tomography, the starting model is comparable to u_true, and the calculated model is comparable to u_est, so a comparison of the two models gives an idea of resolution.

4 Case study

The seismic data used in this work were collected by National Institute for Occupational Safety and Health (NIOSH) in a coal longwall mine in western United States [16]. The mining-induced microseismic events during eighteen days were recorded by sixteen mounted receivers on the surface, above the active panel. The locations of longwall face in each day, panel geometry and location of the mounted receivers on the surface are shown in Fig. 3. In this mine, the average annual mining production is 7.5×10⁶ t. The overburden is approximately 350 m, and the coal seam thickness ranges from 2.6 m to 3.0 m. The coal seam is overlain and underlain by sandstone with 2 m and 35 m in thickness, respectively. The length and width of interested panel are approximately 5 490 m and 250 m, respectively. The adjacent panel in tailgate already is mined, but the mining operation in the adjacent panel in headgate has not yet started. The mining operation in the panel under study is retreated longwall mining method. Over the course of data recording, the face is retreated for about 431 m, so about 24 m per day in average [9, 16-17].

Fig. 3 Locations of longwall face in each day (lines), panel geometry and locations of receivers (circle)

The collected raw data during eighteen days period consist of 172 632 p-wave arrival times, and 11 696 microseismic events. The recorded events typically have a frequency around 30 Hz. The number of recorded events per day depends on the amount of mining operation, ranging about 100 to more than 800 events per day, with a mean of 650 [16].

To avoid creating artificial anomaly in seismic velocity tomograms, the events recorded by less than teen of sixteen receivers are omitted. In tomography study and determination of the data location, the SIRT algorithm is employed. Generally, for such study, the SIRT is an appropriate method, because the solution tending to both converge and diverge is slowly, so the solution is relatively stable [18]. Voxel size of 15 m ×  15 m × 15 m is considered. This size is small enough to interpretate the velocity changes, and also is large enough to avoid creating artificial anomaly. The ideal size of voxel is equal to the wavelength of seismic wave and for smallest size the half wavelength has been suggested [19]. In this work, the average wavelength of recorded events is 120 m and therefore the ideal size of voxels is 120 m; but the expected features of velocity in this voxel size are not detectable. However, each voxel with 15 m in size is traversed by more than 1 000 rays. A generated three-dimensional model of p-wave velocity in rock mass based on geophysical field data that are collected by NIOSH [9, 16-17] is used as the initial velocity model. SIRT is an iterative method and the initial velocity value in first iteration must perturb [18]. The one-dimensional interpretation of the initial velocity model is displayed in Fig. 4.

Fig. 4 One-dimensional interpretation of initial velocity model (Coal seam is shown by bold black line at approximately     1 700 m elevation)

Anisotropy is defined as variations in the properties of material with direction of measurement [5]. In tomography, anisotropy refers to that the variation of p-wave velocity parallel or perpendicular to the bedding layers is measured. The anisotropy vector is assumed normal to the dipping layers of the initial velocity model, and according to Ref. [16], [-0.068, 0.057, 0.996] is defined. The magnitude of anisotropy refers to the ratio of the velocity measured orthogonally to the anisotropy vector to the velocity that is measured parallel to the anisotropy vector. The magnitude of anisotropy is determined experimentally by inverting the data, ranging from 0.8 to 1.2, with the goal of minimizing the travel time residuals from the inversion [6]. Thus, the selected magnitude of anisotropy is 1.1 to minimize the travel time residuals, in other words, to improve the model that has a better fits on the data. This value indicates that the seismic wave velocity that is measured orthogonally to the anisotropy vector is 1.1 times faster than that parallel to the vector.

In SIRT tomographic inversion, the raypath can be assumed to be straight or curved. The straight ray based on straight distance between the source and the receiver is simply calculated, while the calculation of curved ray is needed to conduct bend calculation according to SNELL’s law [11, 14]. Although, the root-mean-square of travelling time residuals for the straight ray assumption is actually smaller than that for the curved ray assumption, the sum of the residuals for the curved ray assumption is significantly small. Moreover, the higher sum of the residuals for the straight ray assumption demonstrates that the straight ray algorithm consistently underestimates the length of raypath. Therefore, the curved ray assumption for these data is appropriated, and each tomogram with 10 curved ray iterations is generated.

The smoothing constant of 0.02 is applied in all directions. Smoothing replaces the velocity at a node by a weighted average of velocity at that node and the surrounding nodes.

After data processing and evaluating the resolution and sharpness of tomograms based on resolution matrix, the three-dimensional tomography images for each day of study period are provided. To investigate the state of stress redistribution, these images into the coal seam level (depth of 350 m) are sliced. Because the wave propagation velocity through rock mass is related to the stress level, the changes of stress redistribution are inferable from seismic velocity tomograms. Since the use of passive seismic velocity tomography in detecting the state of stress is an innovative technique and it is in development stage, the numerical modeling is employed to verify the tomography results.

5 Numerical modeling

Although the analytical method is a useful tool in geotechnical problems, it also has some limitations. These methods should be used in a domain of assumptions that have been developed based on them. These assumptions usually include the elastic, homogeneous, and isotropic behaviors, while these features do not often exist in the actual problems. In this condition, the use of numerical methods can be an appropriate solution. The boundary element method (BEM) is one of the conventional numerical methods for the continuum domains. In BEM, only the boundaries of the continuum need to be discretized. Also, if the medium extends to infinity, which is common in problems in geomechanics, no artificial boundaries are required. The BEM automatically satisfies far-field conditions. In the BEM, the solution is approximated at the boundaries while equilibrium and compatibility are exactly satisfied in the interior of the medium [19].

To model the stress redistribution around the longwall mining panel, the LaModel [20] as a BEM code was used. LaModel was developed by NIOSH and it is a useful tool for modeling stresses and convergences in underground coal mines as well as optimizing pillar sizes. This software simulates the overburden as a stack of homogeneous isotropic layers with frictionless interfaces, and with each layer having the identical elastic modulus, Poisson ratio, and thickness. This “homogeneous stratification” formulation does not require specific material properties for each individual layer, and yet it still provides a realistic suppleness to the overburden that is not possible with the classic, homogeneous isotropic elastic overburden used in previous boundary element formulations. The LaModel modeling process utilizes three different programs in series: LamPre, LaModel, and LamPlt. LamPre is used to input the project parameters such as the seam thickness, overburden depth, material properties, and the seam geometry. This analysis utilized a mine with one seam and two mining stages. The default overburden and rock mass properties were accepted with the exception of using a 25 m lamination layer thickness. This thickness was optimized experimentally to match existing theory. The default overburden and rock mass properties included the values which are listed in Table 2. Next, the seam geometry and boundary conditions were entered into LamPre.

The required seam geometry parameters included the element width (the width of each grid cell in the seam geometry graphical interface screen), the number of elements in the x-direction, and the number of elements in the y-direction. The element width was chosen to be  3 m with 210 elements chosen for the x-direction and 200 elements for the y-direction. Through this input screen, the seam thickness was set to be 2.7 m and the overburden depth was set to be 350 m [9, 16-17]. Next, the seam boundary conditions were set to be rigid on the head and tail sides and rigid on the top and bottom sides of the model. The final stage of the LamPre parameter input process calculates material properties for the coal and for the gob based on the previous inputs. After these material properties are calculated, the user is prompted to input the layout for each mining stage. In this modeling, a longwall mine includes the panel, two side entries and also chain pillars are defined.

Table 2 LamPre input parameters [16]

The output file of the LamPre was used as an input for the LaModel. The state of stress was calculated by the LaModel and finally by using the LamPlt, and the results were displayed. The numerical plots of the stress redistribution around the longwall panel for eighteen days of study with two-dimensional tomograms are shown in Figs. 5-7. The efficiency of passive seismic tomography in inferring the stress redistribution will be confirmed by comparison between the numerical plots and tomograms.

6 Results and discussion

Based on mining-induced seismic data and SIRT inversion algorithm, the three-dimensional seismic velocity tomograms for each day of study period are provided. To display the velocity changes in coal seam level and subsequently to deduce the stress  redistribution, these tomograms into the coal seam level (depth of 350 m) are sliced. Two-dimensional images of velocity variations in coal seam level with longwall panel geometry and face location for eighteen days of study course, with numerical plots of the state of stress are shown in Figs. 5-7. As can be seen, the face from northeast to southwest is retreated.

Generally, the p-wave velocity variations are linked to the changes of stress in rock mass. With the increase of stress, rock gradually is more compact, the structural discontinuities are closed and hence the velocity of wave propagation increases. However, it should be noted that if the stress exceeds a certain limit, causing new fractures expansion in the rock mass, then the velocity of wave propagation decreases. This is a major challenge to infer the stress by using seismic tomography, because in tomograms, the very high and very low stress almost have the same images and the diagnosis of these areas is difficult. In this situation, the experience and engineering judgments of interpreter can be helpful.

Fig. 5 Two-dimensional seismic velocity tomograms and numerical plots of stress (1-6 d)

Numerical plots show that the front abutment is the widest on the tailgate side and shows a maximum value at the center of the face. The pillars adjacent to the gob on the tailgate side have yielded around the edges and now show a high stress concentration at their centers. The magnitude of the tailgate-side abutment is greater than the magnitude of the headgate-side abutment and extends beyond the face due to the adjacent mined-out panel. Behind the face line, the stress level greatly reduces and in gob area, a fairly stressed zone is illustrated. In addition, abutment pressures move along the face with advance of the mining operation.

In tomograms, similar to numerical plots, a high velocity zone is seen on the tailgate which is based on relationship between stress and wave propagation velocity, representing a high stressed zone as the side abutment pressure. Also, a high velocity zone is depicted in the front face that represents the front abutment pressure. On most days, this zone is biased toward the tailgate. As the theoretical model (Fig. 1) is expressed and also numerical plots are shown, the intersection of the tailgate and face is an over stressed area. In addition, a low velocity zone, behind the face line (or around the face line) is displayed as a stress reduction behind the longwall face that is described in theoretical concept of stress redistribution around the longwall mining panel (Fig. 1). A fairly high velocity zone is also illustrated on the headgate that is known as the side abutment pressure on headgate. The state of stress in gob area is complex. Since the gob area contains broken material, the wave propagation velocity in this area is very low. However, with advancement of longwall operation, the gob areas are gradually compacted and hence, the wave velocity increases. This situation with the tomographic images of successive days is evident. However, in numerical plots, the state of stress in gob area is uniform. Moreover, in tomograms similar to theoretical concepts and numerical plots, the stress zones shift in the direction of the face advance.

Fig. 6 Two-dimensional seismic velocity tomograms and numerical plots of stress (7-12 d)

Fig. 7 Two-dimensional seismic velocity tomograms and numerical plots of stress (13-18 d)

In numerical modeling, many parameters are assumed constant, which is likely not to be matched with the actual conditions. Sometimes, this situation will affect the numerical modeling results. For example, based on theoretical concepts, the stress redistribution around the longwall panel is related to the advance rate of the face, but it is not presented in numerical plots. However, in tomography, the actual conditions of filed are considered and this situation is demonstrated. So, with decreasing the advance rate of the face (6th, 9th, and especially 13th days), the pattern of stress around the longwall panel is changed and generally increases. The effect of this stress increase is the fracture that is demonstrated in the tomograms of following days. Of course, such application of passive seismic velocity tomography to infer the stress redistribution is an innovational approach which is still in development stage and inevitably faces with some restrictions.

Generally, the changes of stress redistribution around the longwall mining panel by comparison of the tomograms of successive days are demonstrated. Moreover, the highly stressed zones that are usually the prone area of rockburst, roof fall, and uncontrolled failures are detectable in tomograms.

7 Conclusions

1) In passive seismic tomograms, the main stress redistribution zones around the panel, including front and side abutment pressures, and gob stress are detectable, and the stress zones shift in the direction of the face advance.

2) The simplifying assumptions in analytical and numerical modeling, reduce the accuracy of the results. However, the passive seismic velocity tomograms are provided based on actual conditions and as a result, show the changes of stress redistribution in more detail.

3) Tomography images show that the pattern of stress changes around the longwall mining panel is dependent on the conditions of previous days. However in numerical plots, this dependence is not seen.

4) In tomograms, the gob area in the tailgate side, as expected, is more compacted. Because in most days, more stresses are applied on tailgate. However, in numerical stress plots, the gob area in the headgate side is more compacted based on the theoretical concepts.

References

[1] PENG S S. Coal mine ground control [M]. 3rd ed. Englewood, USA: Society for Mining, Metallurgy, and Exploration, Inc. (SME), 2008.

[2] PENG S S. Longwall mining [M]. 2nd ed. Englewood, USA: Society for Mining, Metallurgy, and Exploration, Inc. (SME), 2006.

[3] JEREMIC M L. Strata mechanics in coal mining [M]. Rotterdam, Netherlands: A. A. BALKEMA, 1985.

[4] YOUNG R P, MAXWELL S C. Seismic characterization of a highly stressed rock mass using tomographic imaging and induced seismicity [J]. Journal of Geophysical Research, 1992, 97(B9): 12361-12373.

[5] BRADY B H G, BROWN E T. Rock mechanics for underground mining [M]. London: George Allen & Unwin, 2004.

[6] COX M. Static corrections for seismic reflection surveys [M]. Tulsa: Society of Exploration Geophysicists, 1999.

[7] LEE W H K, PEREYRA V. Mathematical introduction to seismic tomography, inseismic tomography: Theory and practice [M]. London: Chapman and Hall, 1993: 9-22.

[8] SWANSON P, ESTEY L, BOLER F, BILLINGTON S. Accuracy and precision of microseismic event locations in rock burst research studies [R]. Report of investigations (U.S. Bureau of Mines) No. 9395, NIOSH, 1992.

[9] LUXBACHER K, WESTMAN E, SWANSON P, KARFAKIS M. Three-dimensional time-lapse velocity tomography of an underground longwall panel [J]. International Journal of Rock Mechanics and Mining Sciences, 2008, 45(4): 478-485.

[10] MENKE W. Geophysical data analysis: Discrete inverse theory [M]. 2nd ed. San Diego: International Geophysics Series, Harcourt Brace Jovanovich, 1989.

[11] TARANTOLA A. Inverse problem theory: Methods for data fitting and model parameter estimation [M]. Amsterdam: Elsevier, 1987.

[12] SANTAMARINA J C, FRATTA D. Discrete signals and inverse problem: An introduction for engineers and scientists [M]. New York: Wiley, 2005.

[13] SOLDATI G, BOSHI L. The resolution of whole earth seismic tomographic models [J]. Geophysical Journal International, 2005, 161(1): 143-153.

[14] WATANABE T, MATSUOKA T, ASHIDA Y. Seismic travel-time tomography using fresnel volume approach [C]// Proceedings of the 69th SEG Annual Meeting. Houston: SEG, 1999: 1402-1405.

[15] TSELENTIS G A, SERPETSIDAKI A, MARTAKIS N, SOKOS E, PARASKEVOPOULOS P, KAPOTAS S. Local high-resolution passive seismic tomography and kohonen neural networks- Application at the rio-antirio strait, central greece [J]. Geophysics Journal, 2007, 72(4): 93-106.

[16] MSHA. Data retrieval system [EB/OL]. http://www.msha.gov/drs/ drshome.htm. [2007-24-10]. 2010.

[17] LUXBACHER K, WESTMAN E, SWANSON P. Time-lapse tomography of a longwall panel: A comparison of location schemes [C]// Proceedings of the 26th International Conference on Ground Control in Mining. Morgantown, West Virginia University: ICGCM, 2008: 217-225.

[18] NOLET G. A breviary of seismic tomography [M]. UK: Cambridge University Press, 2008.

[19] BOBET A. Numerical methods in geomechanics [J]. The Arabian Journal for Science and Engineering, 2010: 35(1B): 27-48.

[20] NIOSH. LAMODEL software-stress and displacement calculations [EB/OL]. [2003-03-09]. http://origin.cdc.gov/niosh/mining/ products/product54.htm.

(Edited by DENG Lü-xiang)

Received date: 2011-10-12; Accepted date: 2012-02-01

Corresponding author: N. Hosseini, PhD; Tel: +98-9124711842; E-mail: navid.hosseini@gmail.com