Arch. Min. Sci., Vol. 57 (2012), No 1, p. 139155 Electronic version (in color) of this paper is available: http://mining.archives.pl DOI 10.2478/v10267-012-0010-9 NAVID HOSSEINI*, KAZEM ORAEE**, KOUROSH SHAHRIAR***, KAMRAN GOSHTASBI**** PASSIVE SEISMIC VELOCITY TOMOGRAPHY AND GEOSTATISTICAL SIMULATION ON LONGWALL MINING PANEL TOMOGRAFIA PASYWNA POLA PRDKOCI I SYMULACJE GEOSTATYSTYCZNE W OBRBIE POLA CIANOWEGO Generally, the accurate determination of the stress in surrounding rock mass of underground mining area has an important role in stability and ground control. In this paper stress redistribution around the longwall face has been studied using passive seismic velocity tomography based on Simultaneous Iterative Reconstructive Technique (SIRT) and Sequential Gaussian Simulation (SGS). The mining-induced microseismic events are used as a passive source. Since such sources are used, the ray coverage is insufficient and in order to resolve this deficiency, the wave velocity is estimated in a denser network and by the SGS method. Consequently the three-dimensional images of wave velocity are created and sliced into the coal seam. To analyze the variations of stress around the panel during the study period, these images are interpreted. Results show that the state of stress redistribution around the longwall panel can be deduced from these velocity images. In addition, movements of the stressed zones, including front and side abutments and the goaf area, along the longwall face are evident. The applied approach illustrated in this paper can be used as a useful method to monitoring the stress changes around the longwall face continuously. This can have significant safety implications and contribute to improvements in operational productivity. Keywords: Longwall Mining, Passive Seismic Tomography, SIRT, SGS, Stress Redistribution Dokladne okrelenie napre górotworu w warstwach otaczajcych wyrobiska podziemne ma podstawowe znaczenie dla stabilnoci i zabezpieczenia powierzchni. W artykule tym zbadano rozklady napre wokól przodka cianowego przy wykorzystaniu tomografii pola prdkoci, w oparciu o techniki rekonstrukcji przy równoczesnych iteracjach i sekwencyjnej symulacji Gaussa (SSG). Drobne wydarzenia mikrosejsmiczne wykorzystane zostaly jako ródla bierne. Z uwagi na wykorzystanie takich ródel, zasig promieni jest niewystarczajcy, dlatego te prdko fali okrelana jest przy uyciu gstszej siatki * DEPARTMENT OF MINING ENGINEERING, SCIENCE AND RESEARCH BRANCH, ISLAMIC AZAD UNIVERSITY, TEHRAN, IRAN; firstname.lastname@example.org ** STIRLING UNIVERSITY, STIRLING, SCOTLAND, UK; email@example.com *** DEPARTMENT OF MINING AND METALLURGICAL ENGINEERING, AMIRKABIR UNIVERSITY OF TECHNOLOGY, TEHRAN, IRAN; firstname.lastname@example.org **** DEPARTMENT OF MINING ENGINEERING, TARBIAT MODARES UNIVERSITY, TEHRAN, IRAN; email@example.com i w oparciu o metody sekwencyjnej symulacji Gaussa. W rezultacie otrzymujemy trójwymiarowe obrazy prdkoci fali, które nastpnie ,,narzucane" s warstwami na poklad wgla. Okrelenie naprenia wokól pola w trakcie badania wymaga interpretacji tych obrazów. Wyniki wskazuj, e rozklad stanu naprenia wokól ciany mona okreli na podstawie obrazów prdkoci. Ponadto, uwidaczniaj si ruchy stref podlegajcych napreniom, w tym warstw ssiadujcych przednich i bocznych oraz obszaru samego wyrobiska wzdlu przodka ciany. Zastosowane podejcie moe zosta wykorzystane jako skuteczna metoda biecego monitorowania zmian napre w rejonie ciany. Ma to powane znaczenie z punktu widzenia bezpieczestwa pracy, a co za tym idzie przyczynia si do podniesienia wydobycia. Slowa kluczowe: wybieranie cianowe, tomografia pola prdkoci, metody rekonstrukcji przy równoczesnych iteracjach, symulacja sekwencyjna Gaussa, rozklad napre Introduction Mining operation, especially underground coal mining, always has the remarkable risks of ground control. Evaluation the state of stress and it's comparison with bearing capacity of rock mass can provide a criterion of stability. Moreover in mining process, the stress redistribution has a significant effect on the preparation, development, and extraction (Peng, 2006). In a longwall coal mine, the fracturing-property of coal seam increases with increasing the stress which in result, improves the productivity of operation, but excess in stress not only does not help to improve the production but also the mining operation will be dramatically impaired (Jeremic, 1985). Therefore, the knowledge of stress redistribution around the longwall panel has an important role in reforming and improving of mine condition. The conventional methods of stress estimation are based on static field measurements and/or analytical calculation but each method has some limitations (Peng, 2006, 2008; Jeremic, 1985). Actually, the mechanical status evaluation of rock mass can be relied on the dynamic properties. Although the theoretical concept of wave propagation through the rock mass is associated with some complications, but its implementation as an engineering tool, have unique advantages than other methods. In this paper, the passive seismic velocity tomography is 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 variations will be possible. In several studies the ground failure process by analysis of microseismic events is predicted (Eneva et al., 1996, 1998; Iannacchione et al., 2005b). Field and laboratory studies show that the seismic tomography is a useful tool to analysis the stress changes in rocks. Such as; the stress distribution consideration in the laboratory samples (Meglis et al., 2004; Chapman, 2004; Nur & Simmons, 1969; Eberhart-Phillips et al., 1989; Scott et al. 1994; Mitra & Westman, 2009; Chow et al., 1995; Nakagawa et al., 2000), the pillar stability (Young & Maxwell, 1992; Friedel et al., 1996; Maxwell & Young, 1993; Scott et al., 1999, 2004; Manthei, 1997; Obert & Duvall, 1967; Obert, 1941), the state of stress and stability around the entries, tunnels, and other underground spaces (Watanabe & Sassa, 1996; Maxwell & Young, 1995, 1996; Wright et al., 2000; Friedel et al., 1997; Glazer & Lurka, 2007; Zuo et al., 2009; Gale et al., 2001; Cook, 1963; Prugger & Gendzwill, 1993; Scott et al., 2004), the rock burst (Scott et al., 1997; Brady, 1977; Li et al., 2007; Alcott et al., 1998; Lurka 2008), and analysis the stress in coal mines (Kormendi et al., 1986; Westman et al., 1996; Hebblewhite & Simpson, 2001; Mason, 1981; Hanna & Haramy, 1998; Hanson et al., 2002; Harris & Akintunde, 2005; Iannacchione et al., 2005a; Wong et al., 1989; Westman et al., 2008; Luxbacher et al., 2008a, 2008b; Luo et al., 2008, 2009; Lurka & Swanson, 2009). Usually, during the pre-failure regime, the p-wave velocity increases linearly with stress at lower stress levels, and then plateaus at higher stress levels. Because the cracks and pore spaces are closure with increases the stress and therefore the p-wave velocity is related to stress (Young & Maxwell, 1992). Although the use of passive seismic tomography has many advantages, it is also associated with some limitations. Lack of control on the location, magnitude and other features of seismic events are the most important limitations. Because the high scattering and non-uniform of mining-induced microseismic events causes inadequate ray coverage in some area and thus in many cases, the quality of tomographic images to interpretations and stress inferring would not be appropriate (Luxbacher et al., 2008a; Luxbacher et al., 2008b). In this paper, firstly by using tomography with Simultaneous Iterative Reconstructive Technique (SIRT) inversion method, the location and velocity of seismic events are calculated. Then, to reach a better ray coverage, the p-wave velocity in a denser network by using the Sequential Gaussian Simulation (SGS) method is estimated, and finally the variations of p-wave velocity in rock mass are imaged. Stress Redistribution around the Longwall Mining Panel To analysis the stress redistribution in the longwall structure and ground strata, a simplified model based on composition of stress dynamic and displacement mechanic is used. The loads that 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 (Jeremic, 1985). In Figure 1 the stress redistribution around the longwall panel, especially front and side abutment pressures are shown. In addition, the overlap of side and front abutment pressure are displayed on both side Fig. 1. General form of stress redistribution around longwall mining panel (Peng, 2006) of the face. Generally, the situation of abutment pressure is different; particularly the location of maximum pressure in various longwall coal mines that cause to the state of stress in rock mass surrounding panel and also geological features of each coalfield (Peng, 2006, 2008; Jeremic, 1985). In Figure 2, zones of stress around longwall mining panel and also, in Figure 3 the distribution of abutment pressure in section AA and BB of Figure 2, are shown. As seen in Figure 3, along section B-B, the peak of vertical stress on the longwall panel is occurred in a short distance of front the face, and immediately behind the face reaches to zero. Therefore, the vertical stress concentration problem at behind the face usually does not exist; meanwhile the shields have a high vertical support capacity. Since roof strata after the fractured are still attached to each other, so usually the vertical stress in goaf is less than average of overburden pressure (Jeremic, 1985). Section A-A shows the distribution of side abutment pressure. The maximum of vertical pressure is occurred behind the face, where a large length of the roof is not supported by shields. In this area, abutment pressure is relatively uniform (Peng, 2006; Jeremic, 1985). Fig. 2. The zones of abutment pressures around longwall mining panel Fig. 3. Abutment pressures distribution for sections A-A and B-B of Figure 2; X shows the average of in-situ stress Distribution of abutment pressures around the longwall panel is depend on various parameters such as the properties of coal seam and surrounding strata, mine depth, panel geometry, and type of the support system (Peng, 2006). However, the overburden pressure has the most effect on the magnitude of abutment pressure (Jeremic, 1985). Thus, with increasing the depth of mine, the greater abutment pressure will be expected. Tomography Radon was the first to establish the tomography framework. He proved that with passing an infinite number of rays through a two-dimensional object at an infinite number of angles could generate a perfectly tomography image of the object. This theory also was applied for three- dimensional objects (Cox, 1999). Tomography includes body dividing 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 (Lee & Pereyra, 1993; Cox, 1999). Tomography based on the type of used source, is classified as "active" and "passive" (Swanson et al., 1992). 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 and hence no disruption in the mine production, and also will not impose additional costs (Hardy, 2003; Luxbacher et al., 2008a). But in use of the passive source, uncontrolled features of wave including the location, occurrence time, and magnitude are the most challenges (Glazer & Lurka, 2007). 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 (Luxbacher et al., 2008b) as shown in Equations 1 to 3: v= t= l ® vt = l t dl = R S (1) R1 S pdl (2) ti = å pj lij j =1 (i = 1,..., N ) (3) where, v is the velocity, l is the distance between the source and the receiver, t is the travel time, p is the slowness, N is the total number of rays, M is the number of voxels, ti is the travel time of the i th ray, pj is the slowness of the j th voxel, and lij 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, but 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 (Luxbacher et al., 2008b). If the time, distance, and slowness for each voxel arranging into the matrix form, the velocity based on inverse theory is calculated as Equation 4: -1 T = LP ® P = L T (4) where, T is the travel time per ray matrix (1 × N), L is the distance per ray per voxel matrix (N × M), and 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 (Manthei, 1997; Menke, 1987). SIRT is one of the famous methods of iterative process and include the following steps (Tarantola, 1987): 1. Ray tracing. 2. Calculation the ray distance in voxels that ray passes through them. 3. Calculation the residual time for ray (observed time minus calculated time) based on slowness distribution and save the results. 4. Repeat steps 1 to 3 for all rays. 5. Modification the slowness in each cell with regard to all the passing rays. Repeat steps 1 to 5 continue 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 Equation 5: pj k +1 = pj k + å ( ei ) M j lij å (lij) å lij lij N (5) Equation 6 shows the matrix form of Equation 5: -1 k +1 =p T +y . Q . A . e (6) where, p k +1 and p k are the slowness vector (M × 1) after the k + 1 th and the k th iterative N respectively, y jj = é å lijù ê i ú ë û is the diagonal matrix (M × M ), ij = [lij]T is the transpose matrix 2 M (N × M ), Aii = éå (lij) ù is the diagonal matrix (N × N), e k is the residual vector (N × 1) ê j ú ë û after k th iterative, N is the total number of rays (number of equations), M is the number of voxels (number of parameters), pjk is the slowness of the j th voxel after the k th iterative, and (ei) k is the residual of the i th ray after the k th iterative (Santamarina & Fratta, 2005). -1 Case Study The used data in this study were collected by National Institute for Occupational Safety and Health (NIOSH) at an underground coal mine in western United Stats, over an eighteen days period (MSHA, 2010). In this mine, the retreating longwall mining method is employed. The coal seam thickness is ranging from 2.6 to 3.0 m, and the mine depth is approximately 350 m. The coal seam is overlain and underlain by sandstone with 2 m and 35 m thickness respectively. The length and width of interested panel are 5490 and 250 m respectively. The adjacent panel at tailgate side already is mined, but the mining operation in adjacent panel at headgate side is not started yet. During microseismic events recording, the face approximately 431 m, averaging daily 24 m is retreated. To record the mining-induced microseismic events, sixteen receivers on the surface above the active longwall panel are mounted. The raw data includes 172632 p-wave arrival times, and 11696 microseismic events over the course of this study. The velocity and the frequency of the seismic events are averagely 3600 m/s and 30 Hz respectively. The number of event recorded per day is ranging from 100 to 800 events, with average daily 650 events (MSHA, 2010; Luxbacher et al., 2008a; Luxbacher et al., 2008b). In Figure 4 the location of mounted receivers on the surface, raypaths, and also longwall panel geometry are depicted. Fig. 4. The location of receivers (red circulars), ray paths, and panel geometry To avoid creating artificial anomaly in seismic velocity tomograms, the events that is 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 tends to both converge and diverge is slowly, so the solution is relatively stable (Nolet, 2008). Voxel size 15 m × 15 m × 15 m is considered. This size was determined to be sufficiently small to ascertain the general stress trend, but sufficiently large that low and high velocity artifacts would not disrupt interpretation of the tomograms. The ideal size of voxel is equal to the wavelength of seismic wave and for smallest size the half wavelength has been suggested (Hsieh, 2003). Since the average p-wave velocity for this data set is 3600 m/s with a typical frequency content of 30 Hz, hence 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 is not detectable. However each voxels with 15 m size is traversed by more than 1000 rays. SIRT is an iterative technique; the algorithm must have an initial velocity value to perturb the first iteration (Nolet, 2008). The initial velocity model allows the inversion to be calculated more efficiently and accurately. The initial velocity model was provided with the raw data from NIOSH (MSHA, 2010; Luxbacher et al., 2008a), and a one-dimensional interpretation of the model is displayed in Figure 5. Anisotropy is defined as variations in the properties of material with direction of measurement (Obert & Duvall, 1967). In tomography, anisotropy refers to variation of p-wave velocity when parallel or perpendicular to the bedding layers are measured. The anisotropy vector is assumed normal to the dipping layers of the initial velocity model, and according NIOSH provided data as [-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, ranges from 0.8 to 1.2, with the goal of minimizing the travel time residuals from the inversion (Cox, 1999). Thus, the selected magnitude of anisotropy is 1.1 to minimize the travel time residuals, or in other words to improves the model that have a better fits on the data. This value indicates the seismic wave velocity that is measured orthogonally to the anisotropy vector is 1.1 times faster than parallel to the vector. In SIRT tomographic inversion the raypath can assume both straight and curved. The straight ray based on straight distance between the source and the receiver is simply calculated, while the calculation of curved ray is need to bending calculation according to Snell's law (Tarantola, 1987; Watanabe et al., 1999). Although, the root-meansquare of travel time residuals for the straight ray assumption is actually smaller than for the curved ray assumption, but the sum of the residuals for the curved ray assumption is significantly smaller. Moreover the higher sum of the residuals for the straight ray assumption demonstrates which the straight ray algorithm consistently underestimates the length of raypath. Therefore, the curved ray assumption for these data was appropriated, and each tomogram with 10 curved ray iterations was generated. The smoothing constant of 0.02 was applied in all directions. Smoothing replaces the velocity at a node by a weighted average of velocity at that node and the surrounding nodes. Fig. 5. The one-dimensional interpretation of the initial velocity model; coal seam is shown by bold black line at approximately 1700 m Geostatistical Simulation Although various geostatistical estimation methods can estimate with minimum error, but the estimated results have a less variations than initial data (measured data); in other word the estimated results are smoothed. However the geostatistical simulation method generates the consistent data with regionalized variations, so that the histogram and the variogram of simulated data will be quite similar to initial data (Deustch, 2002). SGS is an applicable and flexible geostatistical simulation technique that is used in this study to simulate the p-wave velocity with aim of better ray coverage. To use SGS, the drift from data should be removed and also the data should be normal distribution. After drift removing and converting the data of each day to the standard normal distribution, N(0,1), the variogram by using normalized data is drawn and also for simulation, a network is selected. Then, data on this network is simulated and histogram based on mean and variance of simulation results is drawn. Next, a number is extracted from the histogram randomly, and then the data conditioning is done. This process was carry out for the whole network and finally by using the reverse conversion the data are returned to the initial situation (Cassiani et al., 1998; Deustch, 2002; Lee & Xu, 2000). The wave velocity changes in rock mass have the spatial structure and can be defined as regionalized variations (Deustch, 2002). To simulate the velocity in this study, firstly by statistical tests on the data of each day, the situation of data is analyzed and then the drift from data is removed. Next, the data of each day using Hermite's polynomial functions is converted to standard normal distribution. Variogram is a tool for structural analysis and continual evaluation of regionalized variations (Cassiani et al., 1998). On the data of each day, the non-directional variography is performed and the bi-structural spherical model is fitted. For example, in Figure 6 the variogram of sixth day is plotted. Fig. 6. The variogram of 6th day with fitted bi-structural spherical model In drawn variograms, the nugget effect is very low, because the tomography creates a similarity in the neighboring nods. Based on variography parameters, particularly range effect, the radius of estimation ellipse is determined. Since in the output of tomography algorithm, the velocity is obtained as point samples and each point represents a block with 15 m × 15 m × 15 m dimension, this size has been selected as a block dimensions. Then the velocity variables that have been converted to standard normal distribution, based on variography of structural analysis model, within block model space, by using SGS are simulated10 times. Therefore, for each block, 10 groups of simulated data are available. The simulated data for each day have also standard normal distribution, because the used data in simulation had standard normal distribution. In the next step, using reverse conversion and initial situation, all simulated data is returned to original dimensions. To control the simulation results, the variogram and the histogram of the simulated data and initial data should be similar. For example, the variograms and histograms of 10 simulations with variogram and histogram of initial data for sixth day are shown in Figure 7 and Figure 8 respectively. As can be seen, the comparison of variograms and histograms of initial data and simulated data shows that the simulated results have high similarity with initial data. Fig. 7. Variograms of initial data (bold black line) and 10 simulated data for 6th day Next, the E-type block model for each day is created. In this model, the final average value for each block by averaging the simulated results is calculated (Cassiani et al., 1998; Deustch, 2002; Lee & Xu, 2000). Thus, an average model of interested area is obtained. Fig. 8. Histograms of initial data (Real) and 10 simulated data for 6th day Results and Discussion Using velocity data obtained from tomography and geostatistical simulation, the three-dimensional images of velocity variations for each days of study period are created. Then, to study the velocity variation around the longwall panel, these three-dimensional images into the coal seam level (depth of 350 m) are sliced. In Figure 9, the p-wave velocity images in coal seam level, face location, and longwall panel geometry for eighteen days of study course are shown. Overall, the variation of p-wave velocity in rock mass is related to the state of stress. Although this relationship is complicated but study the velocity variation, especially during a time period, can make the proper understanding of stress redistribution. Increasing the stress in rock mass caused the closure of fractures and pore spaces, and thus the p-wave velocity increases. Therefore high velocity zones in tomogram images are represented the high stress areas. Of course such inferring the state of stress is also associated with some complications. Because increasing the stress cause the development of new fractures in rock mass and thus the p-wave velocity will be decreased (Marcak, 2008). In other word, the variations of velocity in very low Fig. 9. The velocity imaging in coal seam level with longwall panel geometry and face location for eighteen days study course and very high stress zones is somewhat similar, and the interpretation of these areas only with help of engineering judgments and general understanding of the state of stress rely on theoretical models is possible. As seen in Figure 9, the state of stress redistribution around the longwall mining panel is depended on the advance rate of the face. With decreases the advance rate of the face, the pattern of stress around the longwall panel is changed. The effect of this stress increases is the fractures expansion which is demonstrated in the tomograms of following days. Therefore maintaining a uniform trend in the advance rate of the face can be effective in controlling the stress changes. In most days of the study period, a high-velocity zone in front of the face is evident which moving along with the advance of the face. This zone represents the front abutment pressure. Moreover, another high-velocity zone is on the tailgate that represents the side abutment pressure in tailgate. In addition, on the headgate, a low-velocity zone is seen that shows side abutment pressure in headgate. Since an un-mined panel is located in other side of headgate, so the rock mass structure in this area is strong and pillar system is well undergo the applied loads. Thus the high stress zone in headgate is not created. While on the other side of tailgate, the already mined panel is located that is filled with caved materials. The very high-velocity zone in the corners of face, especially at tailgate side, is seen that is matched with theoretical concept of stress distribution around the longwall panel. Additionally in goaf area, a relatively high-velocity zone is seen. Based on the structure of longwall mining, the goaf is consists of fractured rocks (caved materials) and the wave attenuation rate in goaf area is high. However with adequately ray coverage due to using SGS, the study of stress redistribution in goaf is also possible. Rock burst is one of the dangerous phenomena in underground mining that caused by the rock over stressed. With continuous study of velocity variations around the longwall mining panel and inference the stress redistribution on tomographic images, the state of stress concentration will understand; so that the stress concentration zones are recognized and the potential of rock burst to be determined. Conclusions Results of this study show the passive seismic velocity tomography is an appropriate technology for continuous monitoring the stress redistribution around the longwall mining panel. The insufficient ray coverage, especially in a far area of the face, is the most important limitation of this technique, which in this study by using of SGS has improved. Moreover, applying the SIRT to inversion, location of events and tomography is an appropriate approach for determination the stress distribution around the panel. In tomographic images that are created by SIRT and SGS, the front and side abutment pressures and also the stress redistribution on goaf area are clearly evident. The estimated state of stress from the velocity images is matched acceptably with theoretical models. The process of stress redistribution around the longwall panel with study the tomographic images of successive days are understandable. So that the effect of advance rate of the face on stress redistribution and also movement of the stress zone along with the face is obvious. However, the most important limitation of this approach is related to use of the passive source and lack of control on the wave parameters; which with more suitable configuration of receivers array on the surface, particularly the use of more receivers, can be improved. The research conclusion proves that despite of some limitations, this approach is a significant tool in improving of the safety and performance of longwall mining method.
Archives of Mining Sciences – de Gruyter
Published: Oct 29, 2012