RPPLlED GEoPHvsIcs ELSEVIER Journal of Applied Geophysics 34 (1995) l-21 3-D tomographic imaging of anomalous conditions in a deep silver mine M.J. Friedel a, M.J. Jackson a, D.F. Scott b, T.J. Williams b, M.S. Olson a a U.S. Bureau ofMines, Twin Cities Research Center, 5629 Minnehaha Ave. S., Minneapolis, MN 55417, USA b U.S. Bureau of Mines, Spokane Research Center, E. 31.5 Montgomery Avenue, Spokane, WA 99207, USA Received 16 July 1994; accepted 20 April 1995 Abstract Mining-induced stress-field changes pose both safety and economic hazards. In an effort aimed at developing technology for mitigating such hazards, the Bureau of Mines together with Hecla Mining Company conducted an active 3-D seismic tomographic investigation of anomalous rockmass conditions in a large underground, high-grade, remnant ore pillar, at the Lucky Friday mine near Mullan, ID. IRoughly 2400 P-wave traveltime measurements, were simultaneously inverted to obtain a velocity distribution. The resulting velocity structure appears extremely heterogeneous ( 1 S-6.0 km/s) and well correlated with mechan- ical models indicating the transfer of stress in direct response to mining. Regions of anomalous ground (intense fracturing or high-stress) were identified using threshold probabilities; the minimum velocity regions surrounding drifts indicate a zone of stress relief that extends up to three drift diameters into the rockmass, while regions of maximum velocity above and below mined-out portions of veins indicate the likelihood of concentrated stress. Active tomographic imaging provides the engineer with a flexible tool for routine underground 3-D monitoring of mechanical conditions in large mine structures. 1. Introduction 1.1. Problem signijicance The ability to “see”’ and quantify risks associated with existing or developing anomalous rockmass con- ditions in three dimensions gives au engineer the poten- tial to improve both the safety and economics of underground mining. One of the immediate benefits of such a technology is the identification of regions sub- jected to elevated stress. This is particularly important since catastrophic and violent failure of rock surround- ing underground openings, i.e. rockburst, is believed to occur in response to (either an over-stressed structure or the passing of a seismic wave (Blake and Cuvelier, 1990). The safety problems associated with rockbursts are serious, sometimes causing injuries and deaths. 0926-9851/95/$09.50 0 1995 Elsevier Science B.V. All rights reserved SSD10926-985 1(95)00007-O Also, the associated clean up and rock support costs, and lost production time due to rockbursts are immense. Although a variety of countermeasures are in use at the Lucky Friday and other burst-prone mines, no method has proven entirely reliable. Conventional rockburst countermeasures can be passive, e.g. obser- vation of temporal/spatial sequences of rockmass fail- ure using microseismic monitoring, or active, e.g. excavation sequencing, application of counter stresses through bolt installation, and/or destressing through a drill and blast sequence. Despite the use of these coun- termeasures, rockbursts continue, often with no micro- seismic precursors, (Blake and Cuvelier, 1990). As mining proceeds to greater depths, rockbursts will remain a persistent safety hazard and an impediment to production. Presently, there remains the need to develop a simple, rapid, and cost-effective tool for peri- 2 M.J. Friedel et al. /Journal of Applied Geophysics 34 (1995) 1-21 odic imaging of mine structures under the full range of stress conditions. One such promising “high-resolu- tion” geophysical technology is active 3-D seismic tomography. This paper describes the results of a study designed to evaluate the feasibility of imaging stress conditions in a deep mine through the use of seismic tomography. 1.2. Rationale The rationale for using seismic velocity tomograms to aid in mapping mine-induced, stress-related phe- nomena is based largely on laboratory observations showing that P-wave velocity increases under loading (Thill, 1989). The increase in velocity is related fun- damentally to the closure of void space, e.g. pores and cracks. In general, the rate of velocity increase is non- linear and greatest with an early incremental increase in stress. As stress increases further, the rate at which velocity increases is reduced in response to the for- mation of new cracks (yield point) parallel to loading. These observations suggest that regions of high veloc- ity are likely to indicate zones of high stress concentra- tion, whereas low-velocity regions indicate zones of stress relief. 1.3. Background To date, a number of workers have investigated the application of seismic tomography to assist in solving mine-induced stress-related problems. The earliest application was in 1981, when Mason guided com- pressional waves (P-waves) along a mine roof to locate the position of adjacent coal pillars (Mason, 198 1) . A few years later, seismic tomography was used to mon- itor large-scale changes in stress behind an advancing longwall (Kormendi et al., 1986). Similar experiments were conducted in Hungary and Poland to detect anom- alous bump prone regions, i.e. locate stress concentra- tion, in a coal panel (Ivansson, 1985; Dubinski and Dworak, 1989). United States and Canadian researchers have applied improved methods of traveltime tomography to eval- uate mechanical conditions in metal and nonmetal mines (Young and Hill, 1985; Friedel et al., 1992, 1993, 1994; McCreary et al., 1992; Young and Max- well, 1992; Jackson et al., 1993). Friedel et al. (1993, 1994) demonstrated that yield pillar behavior in coal could be monitored as the longwall advanced, throug the use of seismic refraction tomography. In the fir; 3-D tomographic mining application, passive-sourc images were computed from induced seismicity usin a simultaneous linear inversion procedure for bot source location and velocity structure (Maxwell an Young, 1992). This paper describes the application c 3-D active source seismic tomography for assessin anomalous conditions in a large, heterogeneous, ren nant mine pillar at the Lucky Friday mine. 2. Lucky Friday mine 2.1. Local geology The Lucky Friday (LF) mine, owned and operate by Hecla Mining Company, is located in the Coeu d’ Alene Mining District near Mullan, ID. The two pti mary rock formations at the mine are the argillitic S Regis, and the Revette, composed primarily of vitreou quartzite with thinly ( 1 cm to tens of cm) interbedde sericitic quartzite and argillite (Scott, 1993). Th steeply southward-dipping (about 80“)) E-W-trend ing, LF vein is a silver-lead-zinc deposit that continue to depths of more than 1.62 km. In the region of study the LF vein bifurcates into what are known locally a the hanging and footwall veins. Unconfined laborator sonic measurements conducted on intact samples a vitreous quartzite revealed that P-wave velocities aver age about 4.45 km/s. The study region, known as the “remnant ore pi1 lar,” represents a block of ground that is vertical1 defined by previous mining activity (Fig. 1) delimite by the 4900 and 5 100 levels (the level designatio indicates the approximate depth in feet from the SUI face). Note that at some locations mining activit (stoping) has also occurred both below the 4901 (underhand cut-and-fill) and above the 5100 leve (overhand cut-and-fill). In all cases, these features ar backfilled with a low-modulus sand mix (200 psi) consequently, they are likely to act as localized stres risers within this block of ground. The width of this remnant structure is defined by th lateral extent of crosscuts on the two levels, e.g. 4900 87, 93 and 95; and 5100: haulage, 88, 92, while it length is defined by the distance spanning the outermos crosscuts (Fig. 2). In this region, bedding strikes E. h4.J. Friedel et al. /Journal ofApplied Geophysics 34 (1995) 1-21 95 crosscut ?-_---IO Scala. m Sandfilled Stope VE??Z??Z Fig. 1. Vertical cross-section through the main Lucky Friday vein. Region between levels is the study region and is known as the “remnant ore pillar.” W, essentially parallel to the vein, and dips SSW at approximately 60”. The N40W trending South Control fault terminates the western edge of the LF vein. At the time of this survey, mining activity using underhand cut and fill was engaged at the 5500 level, a cover depth of 1700 m, and proceeding to greater depths. 2.2. In-situ stress To assess the general influence of mining activities on the local stress distribution (and thus tomographic outcome), a finite-element model was used to solve the plane strain problem gjiven local mining geometries and far field stress conditions for the 20,000 E and 20,050 E sections of the remnant structure. Noteworthy mine features include the 4900 and 5 100 lateral haulage drifts and ramp. The location of geologic structures, such as the veins and SC fault, were interpolated between levels for the stress modeling. Previous mining activity asso- ciated with stoping has left portions of the vein back- filled with a low-modulus sand composite at top and bottom locations of this structure. The plane strain condition assumes that features are infinite in extent normal to the plane of investigation. In reality, however, the continuity of mining, e.g. drifts and stopes, and geologic features, e.g. veins and fault, exists only over relatively short distances ( < 50 m). The feature that most grossly violates the plane strain condition is the 93 crosscut in the 20,000 E section. Also, displacements are assumed to be continuous sin- gle valued functions that are independent of the direc- tion perpendicular to the solution domain. Obviously, such a displacement field precludes the possibility of gaps, e.g. joints, faults, etc. To introduce fractured regions, such as around drifts and where rockbursting occurred, the adjacent elemental modulus was decreased to half that of the host rock. Other simpli- fying rockmass assumptions include a massive homo- geneous and anisotropic host (principal component oriented 60” from the horizontal, i.e. in the direction of bedding), homogeneous and isotropic veins and sand- fill (no attempt was made to account for the fault), all subject to far field stresses (Whyatt et al., 1993). The far field stresses were assigned based on prior in situ M.J. Friedel et al. /Journal of Applied Geophysics 34 (1995) 1-21 w W W W w W w II I 20,200 N 0 Spad lOCatiOnS -___-----_______-. Argillite fractures 15 meters 0 I I I I Geophone 75 ,$ _- - - $5 92 1 60 ??Scale locations (5100 level ( 0 Vein I -.. 4i Fig. 2. Lucky Friday mine geometry and geologic structure on (a) the 4900 and (b) the 5 100 levels. 6 M.J. Friedel et ul. /Journal ofApplied Geophysics 34 (1995) 1-21 tests, i.e. the vertical and horizontal stresses were set to lithostatic and two times this value, respectively (stress ratio of 2). Fig. 3 shows the finite-element grid and resulting maximum principal compressional stress distributions ( 10-70 MPa) for the two vertical sections; 20,000 E and 20,050 E. In general, stress appears elevated between the 4900 and 5100 levels with a region of maximum stress ( > 50 MPa) occurring immediately above the sandfilled stopes. Conversely, a region of minimum stress is present across the lower portion of these figures where the ramp and backfilled stopes are present. The interpretation is that the sandfill cannot support any stress, and it is shed to the adjacent more competent vein. In the 20,000 E section (Fig. 3a), the zone of stress concentration is more pervasive, extending from the 5 1 OO-level sandfill to the end of the 4900 level 93 crosscut. In this case, the high stress follows along the dip of the vein until it is influenced by the 93 crosscut, where it turns southward, eventually meeting up with the end of this crosscut. Other anomalies present in the 20,000 E section include an intermediate stress region that wraps from the lower side of the 4900 lateral to the far side of the 93 crosscut. Superimposed on this hem- ispherical shape are several disconnected and higher stress features. Note that these features bound a low- stress region directly below and adjacent to this same drift. Another smaller ( l-2 drift diameters) interme- diate zone of stress occurs to the south of and above the 5 100 crosscut. In the 20,050 E section (Fig. 3b), the zone of max- imum stress concentration ( > 50 MPa) also follows along the dip of the vein ( > 50 MPa); however, it is more localized. Intermediate stresses extend up and outward between the upper (4900 level) and lower (5 100 level) sandfilled stopes. Below the interburden, a low-stress region wraps around the bottom and through the 5 100 level stopes toward the north along the outer perimeter up to the 4900 lateral. Assuming that velocity is a sensitive indicator of stress, the general trends indicated by high and low compressional fea- tures should appear in the tomogram as high and low velocity anomalies, respectively. 2.3. Microseismic measurements In previous studies, a fair correlation was shown to exist between compressional velocity and location of seismicity in active mining regions (Maxwell and Young, 1992; McCreary et al., 1992). In a rockmass, the seismicity is attributed to a sudden release of accu- mulated strain energy. The release of strain energy is then accompanied by the propagation of an elastic wave. Depending on the magnitude of seismic energy and/or its incident angle, this may critically alter the stress state across structural heterogeneities located in the mine vicinity. With this in mind, the Lucky Friday microseismic database (recorded using an MP250 monitoring system) was queried for all events occur- ring in the proximity of the mine structure during the week of the tomographic survey, as well as one week before and after. Only those events with an estimated source location accuracy of less than 6 m were saved for comparative evaluation with the velocity tomo- grams. 3. Data acquisition 3.1. Ray path coverage Perhaps the most important tomographic survey con- siderations are to achieve a high density of seismic ray paths and a wide range of viewing angles through the rockmass. During this 3-D tomographic survey, more than 2400 seismic traces were recorded for ray paths between drifts on common levels, as well as between levels. The global survey represents the compilation of smaller 2-D cross pillar surveys, and level-to-level sur- veys. For example, on the 4900 level, seismic energy traversed pillars defined by the 87-haulage-93, and 93- haulage-95 combinations (Fig. 2a). Similarly, on the 5 100 level seismic energy traversed a pillar defined by the 92-haulage-ramp combinations (Fig. 2b). Follow- ing a cross pillar survey at the upper 4900 level, receiv- ers were left in place, and seismic energy was excited at locations on the lower 5100 level. In total, the trav- eltime response for 1650 interburden seismic ray paths (between the 4900 and 5100 levels) were recorded. The fact that source-receiver geometries surrounding a pillar approximated a horizontal plane, and recogniz- ing that the velocity gradients are probably greatest in this plane, made it possible to invert the traveltime data 2-dimensionally for a same-day evaluation of the gross velocity structure. M.J. Friedel et al. /Journal ofApplied Geophysics 34 (1995) 1-21 I 25 Scale, m I Fig. 3. Modeled maximum compressional stress distribution and mine related features for (a) the 20,000 E and (b) the 20,050 E vertical sections. Modeled assuming a plane strain condition. M.J. Friedel et al. /Journal of Applied Geophysics 34 (1995) l-21 Fig. 4. 3-D seismic ray path orientations assuming straight rays. The individual tomographic surveys were combined, into a global 3-D data set, representing seismic energy that traveled along intersecting ray paths throughout the remnant structure. Fig. 4 shows a stereographic pro- jection of the seismic ray path orientations, which are restricted to inclinations between zero and 10” (sub- horizontal) and between 50 and 90” (steeply inclined). While this appears to indicate a lack of moderately inclined viewing angles, it should be noted that this is compensated to some degree by curvature of rays in the presence of velocity heterogeneities (completeness of the curved ray path coverage is shown later). 3.2. Digital recording Seismic energy was introduced into the rockmass by impacting rock bolts along mine drifts using an 8-lb sledge-hammer as the source. The first arriving seismic energy traveled from the point of impact, across the mine structure, to 24 single-component 100 Hz geo- phones. Each of the geophones was screwed into a rock bolt at varying locations along the mine drift (see Fig. 2a,b). A digital instantaneous floating-point seis- mograph enabled rapid, multichannel data collection and storage, with a total dynamic range of 120 dB . The signal character was optimized by stacking an average of 20 records (each sampled at 0.1 ms intervals, result- ing in a Nyquist frequency of 5 kHz), and applying analog 5 12 Hz highpass and 2 kHz anti-aliasing filters. Fig. 5 shows a typical 24-channel common source gather (following a hammer blow to station R7, in crosscut 95 recorded while surveying the 93-4900-95 d lb i0 40 50 TIME, ms Fig. 5. Typical common source gather: hammer blow (source) at station R7, geophones (receivers) mounted to rock bolts in 93-49- 95 pillar walls. pillar. Reciprocity checks were conducted to test the repeatability of timing. The complete 3-D seismic data acquisition phase took about three g-hour shifts, using a 3-person team. An additional day was required to determine the spatial coordinates for each source and receiver location. The actual amount of time required for a tomographic study depends on a number of factors including the number of seismic records, manpower, number of recording channels, takeout spacing, ambient noise, and other mine related issues, e.g. drilling, blasting, mucking and hauling of the ore, and rockbursts. Experience suggests that optimum tomographic survey conditions are achieved while recording during “off” production time, e.g. during a maintenance shift. Alternatively, quality seismic data can also be recorded in regions that are a minimum of 200 m away from drilling activities. 4. Data reduction and analysis The tomographic information is organized in an ordered format: label, source and receiver coordinates, and traveltimes. The tomographic data file is typically constructed using a spreadsheet. The spreadsheet facil- itates basic data reduction and analysis from which improper coordinate and/or traveltime information, outliers, and statistical trends can be discerned. For any tomographic subset, e.g. 2-D pillar or level-to-level survey, the following five steps are implemented: ( 1) pick traveltimes, (2) combine traveltime and source-receiver coor- dinates, (3) compute apparent velocities, M.J. Friedel et al. /Journal of Applied Geophysics 34 (1995) 1-21 9 40 r- -I 0 25 !jO 75 100 125 150 APPARENT DISTANCE, m Fig. 6. Global time-distance plot; apparent velocities of 1.0 and 3.8 km/s are attributed to fracturing around drifts and intact rock of the interburden, respectively. (4) generate time-distance scatterplots and velocity histograms, and (5) compute statistics. Problems associated with assembling tomographic data sets of this magnitude typically include importing wrong coordinate and /or traveltime groups, i.e. from a common source gather, and traveltime mispicks. While automatic traveltime pickers are commercially availa- ble, manual. selection has provided maximum conti- nuity and fidelity. The tradeoff, however, is that it may take up to one full day for the identification of 2400 traveltimes. Following the removal of outliers identi- fied in the scattetplclts and velocity histograms, the complete 3-D data file is assembled by collating each of the smaller tomographic subsets. 4.1. Time-distance scatterplot After the tomographic data are compiled, a scatter- plot is generated. The scatterplot shows the traveltime (dependent variable) as a function of distance (inde- pendent variable). Since the actual propagation path is not known in advance, the distance between any given source-receiver pair is assumed again to be a straight path. Obvious outliers are detected by visual inspec- tion, and time picks a.nd coordinates are re-examined for possible errors. “Best” straight lines are fit to the data, in the least-squares sense, to identify apparent velocity patterns. For tthe Lucky Friday traveltime data set (Fig. 6), two linear segments may be discerned: a steeply sloping trend for path lengths less than about 8 m, and a more gently sloping trend for longer paths. In a manner analogous to 1-D surface refraction studies, the segments may be associated with structural ele- ments ( “refractors”) in the rock mass, although the interpretation differs somewhat due to the 3-D nature of the data. The apparent velocities (inverse slopes) for the two refractors are 1.0 and 3.8 km/s, respec- tively. Ray paths shorter than 8 m correspond to energy traveling along or near the surfaces of the drifts, where sources and receivers are closely spaced. Thus the slower refractor represents an average global response to travel through the stress relieved (fractured) zone surrounding underground drifts. Using standard refrac- tion equations and simplifying assumptions (Telford et al., 1984)) with a best-fit zero-offset delay of 5.2 ms, the average distance from a drift wall to competent rock is determined to be about 2.5 m. 4.2. Summary statistics Dividing each source-receiver distance by its cor- responding traveltime yielded an apparent velocity for each ray path. Basic statistics for the global data set indicated the following: traveltimes ranged from 1.3 to 31.9 ms, averaging 15.8 ms; survey distances ranged between 2.5 and 140 m, averaging 63.3 m; and veloc- ities ranged from about 1.2 to 6.19 km/s, averaging 3.82 km/s. By dividing the average velocity by the square root of the number of ray paths, we obtain the standard error of 0.08 km/s (this assumes a normal distribution). Significant deviations from the mean velocity (e.g. greater than two standard errors) can be interpreted as mine-induced phenomena, e.g. concen- trated stress, or induced fracturing (stress relief). While the mean value is in general agreement with the compressional velocity determined in regional seismic- ity studies (Stickney and Sprenke, 1992) the average velocity for this study is slightly lower, due to the pres- ence of openings and other low-velocity zones in the mine. Histograms of apparent velocities were used to iden- tify outliers and trends. For example, Fig. 7 shows the computed velocity histograms (0.25 km/s bins) for the 4900 level, 5100 level, level-to-level, and global surveys. In general, all of the velocity histograms are skewed toward lower values. The asymmetry and large standard deviation associated with the 4900 and 5 100 level surveys are attributed to both fracturing in the vicinity of the drifts as well as the presence of backfilled veins. The apparent lack of mining, i.e. no backfilled IO M.J. Friedel et al. /Journal ofApplied Geophysics 34 (1995) 1-21 I Global ? 300 IL E / 4900 Idye Fig. 7. Histograms of apparent velocities for 4900 level, 5 100 level, interburden (4900-5 100 level), and global data sets. stopes, between levels is indicated by an increase in the average velocity and corresponding decrease in stan- dard deviation. The bimodal distribution occurring at the 5 100 level is associated with sampling two distinct regions: the rockburst damaged region south of the fault, and the remaining portion of the pillar. 4.3. Seismic resolution The seismic resolution is related to the wavelength of the propagating energy. While seismic energy with shorter wavelengths (higher frequency) can provide greater resolution, it is more strongly attenuated in rock and limits the overall range of propagation. For trans- mission tomography, the resolution scale is propor- tional to the wavelength, which is theoretically related to the radius of the first Fresnel zone. Current theory suggests that the resolution of traveltime transmission tomography is about one wavelength (Williamson, 1991) . For a frequency of 5 12 Hz (the high-pass filter cutoff) and a mean velocity of 3.81 km/s, the average wavelength is about 7.44 m which provides an estimate of the minimum resolvable feature size. The use of first-arrival data leads to preferential emphasis on high-velocity regions. As a result, low- velocity regions are less well-resolved than fast regions ( Wielandt, 1987), This feature of over-sampling high- velocity regions presumably enhances the confidence in imaging regions of high concentrated stress. Con- versely, those interior low-velocity backfilled regions will likely be under sampled. The fact that source- receivers are placed along the drifts, however, ensures sampling of the stress relieved regions surrounding the drifts. 5. Tomographic imaging 5.1. Constraints Perhaps the most critical issue in transmission tomography is the mathematical nonuniqueness. The limited range of viewing angles in a tomographic sur- vey results mathematically in singular matrices, because the traveltime equations are not linearly inde- pendent. This can be true even when the number of known values (traveltimes) exceeds the number of unknown values (node velocities). Under these cir- cumstances, many different velocity models can satisfy the system of traveltime equations. To limit the range of possible solutions, the set of equations can be sup- plemented with available site information. For exam- ple, global constraints in the form of minimum and/or maximum velocity values are typically used. While a maximum value can be readily obtained from labora- tory stress-velocity tests, the use of a minimum veloc- ity value is somewhat subjective. In lieu of a-priori structural information, it becomes important to evaluate the robustness of the reconstruction with regard to the starting model. In this study, a uniform grid comprising 9261 3-D cells, or voxels (2 1 X 21 X 21 voxels, with dimensions of 6.4 mX5.2 m X4.3 m), was used for the tomo- graphic analysis (Fig. 8). The fact that the total number of computational nodes (10,648) was four times greater than the number of traveltimes (2365) resulted in a grossly underdetermined system of equations. To supplement these incomplete traveltime data, a maxi- mum velocity constraint was imposed. Since the host rock is primarily quartzite, the upper bound is assigned as the velocity through a single quartz crystal (6.0 km/ s) . Robustness was evaluated by performing the inver- sion numerous times, with various alternate starting velocity models. Initial models all had uniform veloc- ities, ranging from 4.76 to 6.0 km/s. Twelve starting velocities were selected between these bounds, with a regular interval pattern. The implementation of alter- nate starting models is thus equivalent to random selec- tion from uniform probability distribution (each equally likely). 5.2. Ray tracing The following tomographic reconstructions were obtained by simultaneous iterative adjustment of alter- M. J. Friedel et al. /Journal of Applied Geophysics 34 (1995) 1-21 11 0 2? Scale, m 50 Fig. 8. Grid used for tomographic inversion; 2 1 X 21 X 2 1 voxels (total of 926 1 voxels and 10,648 nodes). nate starting models subject to various constraints (described above) and ray tracing. The velocity recon- struction resulting from five straight-ray iterations was used as the starting point for the remaining 10 curved- ray iterations. The curved ray tracing used here is based on the method described by Urn and Thurber ( 1987). This approach involves perturbing a trial ray path until the traveltime is minimized. To facilitate the bending of rays outside the minimum survey volume, the com- putational domain is enlarged by 20%. In general, a satisfactory minimum global traveltime root-mean- square (RMS) residual (about 1.45 ms for all recon- structions) was reached following 15 iterations. RMS residuals were approximately 8.6% of the RMS trav- eltime, and residuals for the final models were approx- imately 30% of those for the starting models. Using a 486-66 MHz personal computer, this process required about five hours per reconstruction. Due to the com- putational time required, the total number of starting models was limited to twelve. 6. 3-D tomographic results 6.1. Global velocity b!istribution From the suite of tolmographic reconstructions, basic statistical properties, e.g. average, median, standard deviation, were comlputed for each node. Fig. 9a,b shows selected 2-D lslices through the 3-D average velocity distribution. By taking the average nodal val- ues, poorly-resolved local variations (those that are sensitive to the choice of initial model) are effectively filtered, resulting in a smoothed velocity field that is model-independent and on the average “best”. In gen- eral, conditions appear heterogeneous and anisotropic; the velocity ranges between 1.5 and 6.0 km/s display- ing global structure that dips roughly 60” to the south (coincident with the local geology). Fig. 9c combines the velocity and traced ray paths (for clarity, only 25% of the total number are shown). This figure illustrates the relative degree of spatial sampling, as well as the tendency of first-arriving energy to preferentially sam- ple high-velocity regions. It is interesting to note the relationship of microse- ismicity to local velocity. For example, the histogram, shown in Fig. 10, indicates a tendency for seismic events to occur in regions of intermediate velocity (stress). Of the 1184 seismic events, the majority occur in velocity regions of 4.78 km/s (standard deviation of 0.67 km/s). While this is roughly 25% greater than the average site velocity (3.81 km/s), it is about 22% less than the maximum site velocity. One possible explanation is that the high velocity-stress regions may be acting to preferentially close fractures, thereby pre- venting a release of local strain energy. 6.2. Vertical velocity distribution: 20,000 E Fig. 1 Ia shows the 20,000 E tomogram, a N-S ver- tical slice through the center of the 3-D average velocity M.J. Friedel et al. /Journal ofApplied Geophysics 34 (1995) I-21 1.49 km P 1.55 km . A 1 .49 1 km 1, 55 km B 1.49 1 V,,, km/s Depth C 0 Source/receiver location Fig. 9. Global average 13-D velocity distribution. (a) Vertical slices. (b) Horizontal slices. (c) Traced raypaths (only 25% are shown). model, together with microseismic activity. Velocities are low around the drifts, indicating stress relief, and they are high above the mined-out parts of the veins, indicating stress concentration. Thus there is good agreement between the reconstructed velocity distri- bution and the modeled stress distribution (despite the plane strain assumption). Other anomalously high- velocityregions (4.S6.Okm/s) occurabout 5 mabove and between the 5100 lateral and adjacent sandfilled stopes, at the end of the 93 crosscut, a few meters above but between the 4900 lateral and 93 crosscuts, as well as an increase in gradient below and away from the 4900 lateral and 93 crosscut. These high-velocity fea- tures, interpreted as h:igh stress, are also in agreement with the mechanical model described earlier. Note that the microseismic activity does not occur in these high 250 \"0 1.5 3 4.5 6.5 VELOCITY, km/s Fig. 10. Relationship of Ceismicity ( 1184 events) to local velocity; average and standard deviation are 4.78 and 0.67 km/s, respectively. velocity regions, but in those that reflect intermediate velocity (stress) levels. It is possible, however, that these high-velocity-high-stress regions may also release strain energy following the rotation of the stress field in response to nearby mining activities. Four anomalously low-velocity regions ( 1.5-2.0 km/s) appear in the 20,000 E tomogram; two at the intersection of crosscuts with the 4900 and 5 100 lateral, and at two sites of prior rockburst activity in the 93 crosscut and ramp (Blake and Cuvelier, 1990). The first two low-velocity regions are interpreted in terms of fracturing due to the crosscut curvature at the inter- section with the lateral haulage drift (see Fig. 3). The latter two sites are attributed to partial failure of the rockmass, which will be discussed later. While mine-related stress conditions appear to be identifiable through inspection of velocity distribu- tions, there are practical limitations on resolution that must be kept in mind when interpreting the velocity images. For example, by observing that the low-veloc- ity, stress-relieved envelope is not symmetric about the 4900 and 5 100 lateral drifts, it might be concluded that they are incorrectly located. In reality, however, it is not possible to completely resolve the stress-relieved zone due to asymmetric sampling, i.e. sampling only along one side. Further, it is not possible to directly resolve the sandfilled stopes, veins or fault. In some of the tomograms, however, these features can be located based on their apparent mechanical influence exerted within the remnant structure. The lack of direct reso- 14 M.J. Friedel et al. /Journal ofApplied Geophysics 34 (1995) 1-21 North 4 1 s 5.0 4.5 4.0 3.5 3.0 2.5 2.0 Depth VP, km/s Stress High Low Sandfill 0 Microseismic Event 0 , 25 Scale, m 5p North Depth 1 V,, km/s Stress 5.5 5.0 4.5 4.0 3.5 3.0 2.5 2.0 High Low Sandfill - Microseismic Event Scale, m Fig. 11 Selected vertical velocity slices from the average 3-D tomographic betv vee n velocity and stress (see Fig. 3). reconstruction. E. (b) 20,050 E. Note correspondence (a) 20 ~,OOO M.J. Friedel et al. /Journal of Applied Geophysics 34 (1995) 1-21 1.5 lution of internal low-velocity structures is a conse- quence of preferential sampling of fast rock, as explained earlier. Finally, note that the velocity distri- bution becomes uniform toward the edges of the tomo- gram. This reflects a region through which seismic energy did not pass; consequently, the uniform velocity reflects an average of the starting models used. 6.3. Vertical veloci9 distribution: 20,050 E The general agreement between velocity (tomo- gram, Fig. 1 lb) and stress (mechanical model, Fig. 3) is indicative of the apparent influence of mining on local stress conditions. While regions of anomalously high velocity ( > 4.5 km/s) also occur above the lower sandfilled stopes in this tomogram, the interburden dif- fers substantially from that in the previous tomogram (20,000 E) , particularly in and around the 4900 level. For example, the high-velocity anomaly in this region spans the largest lateral and vertical extent, ultimately transecting the 4900 level (between the 4900 lateral and sandfilled stopes ) . While moderately high veloci- ties may be expected along this section of the pillar defined by the 93-4900-95 drifts, it is the extremely high vetocities ( >5.5 km/s) existing immediately below it that are of potential concern. Another difference between the two vertical velocity distributions is that the structural features in this tomo- gram, such as the veins and fault, appear to exert some form of mechanical control on the concentration of stress. Evidence of this structural control is given by the appearance of these features at locations where high velocity regions ( > 5 km/s) become laterally discon- nected. Also, the microseismicity tends to be concen- trated within intermediate- to high-velocity zones in the region north of the SC fault and between the veins. The two other identifiable features appear as low-veloc- ity anomalies ( <2 km/s): the fractured rockmass related to the 4900 lateral and crosscut junction, and the zone of induced rockburst damage (below the upper right sandfilled stope). This latter observation suggests that seismic tomography may be used to assess the relative volume of damaged rock associated with rock- burst events. 6.4. Horizontal velocity distribution: 4900 Level Fig. 12a shows a plan view of both velocity and location of microseismic activity corresponding to the 4900 level at a depth of 1.49 km. A low-velocity region follows the general drift configuration, with high veloc- ities toward the center of the two pillars; west: 87-4900- 93, and east: 93-4900-95. This fits with the conventional engineering paradigm, where under- ground development results in maximum stress relief (through intense fracturing) nearest to the drift. The increased confinement toward the pillar center contrib- utes to the higher stresses and thus higher velocities. The observation that these high-velocity anomalies occur north of the geologic structure suggests that it may play a role in perturbing stress distribution and confinement. One plausible interpretation is that the west pillar is carrying higher stress than its counterpart, due in part to its inability to dissipate stresses over the structure because of the SC fault proximity. While this may be true, the total high-velocity volume is small compared to that immediately below the smaller east pillar. A comparison of the anomalous high- and low- velocity volumes is presented later using the combined concepts of threshold probability and isosurface. Based on the distribution of low velocities ( 1.5-3.0 km/s), the south side of the 4900 lateral drift appears to be the most intensely fractured of any drift, partic- ularly along the east pillar. At this location, the degree of fracturing extends to roughly three drift diameters. Of the three crosscuts, the 93 appears most fractured with an associated asymmetry; e.g. the east wall is fractured to a distance two times that of the west. More interesting, however, is the presence of a localized anomaly ( < 10 m radius) in each crosscut that tends to follow the strike of the vein and SC fault. These damage zones coincide with the location of prior rock- burst activity, where intensely rubblized rock was observed following a Richter 3.0 magnitude seismic event, and they are aligned with the inferred strike of the fault-slip seismic event. Based on visual inspection, there does not appear to be any particular correlation between current velocity and microseismicity. 6.5. Horizontal velocity distribution: 5100 Level Fig. 12b shows the tomographic plane, passing through the 5 100 level at a depth of 1.55 km., and microseismic activity. At this level, only the east pillar (5 l OO-92-ramp) was tomographically surveyed between drifts. The fact that this pillar had the lowest average velocity is reflected in part by the lack of veloc- M.J. Friedel et al. /Journal ofApplied Geophysics 34 (1995) 1-21 North t V,, km/s Stress 0 Microseismic Event North V,, km/s Stress 0 Microseismic Event Fig. 12. Selected horizontal velocity slices from the average 3-D tomographic reconstruction. (a) 4900 level. (b) 5100 level. High-velocity core illustrates effect of confinement on maintaining elevated stress. Intense fracturing as indicated by low velocities to roughly 3 drift diameters, e.g. along 49 lateral, ramp, and regions that experienced rockbursting (localized areas in each crosscut, just north of geologic structure). M.J. Friedel et al. /Journal of Applied Geophysics 34 (1995) l-21 11 ities above 5 km/s, as well as the largest region of velocities less than 2 km/s. In this pillar, the two small regions (diameter < .5 m) of maximum velocity appear to be separated due to geologic structure. The higher magnitude velocity anomalies appearing in the western pillar are the result of seismic transmission from sta- tions located in the bounding drifts (88-5100-92) to receiver locations on the 4900 level. In this case, the anomalies are likely to be poorly resolved in both mag- nitude and spatial extent due to limited ray-path cov- erage. In general, the low velocity response is similar to that of the 4900 level. The only exception is that the rock mass surrounding the ramp exhibits the greatest intensity and overall extent of fracturing in the remnant structure. It is interesting to note that the majority of seismic events in the east pillar have occurred in the region bounded by the ramp and SC fault. The remain- ing seismicity appears to be associated with secondary faulting outside the tomographic image plane. 7. Tomographic rislk assessment In displaying a tornogram that has the “best” (on the average) global velocity structure, the tomographic analysis becomes deterministic in the sense that it loses the ability to address uncertainty. In any risk assess- ment, deterministic measures typically result in over- estimating a particular outcome. For this reason, a probabilistic approach is implemented to better quan- tify the pending risk associated with anomalous regions determined using an iterative tomographic approach. The basic approach requires the investigating engineer to assign an appropriate velocity threshold. This value can be associated with either stressed or fractured rock. Next, the probability (or likelihood) that a node will exceed the threshold value is computed based on the previously derived nodal statistics. The final display is then contoured as a surface of constant probability, i.e. a probability isosurface. Note that these probability values assume a normal distribution exists, i.e. each nodal velocity distribution is well resolved and has infinite tails. Clearly, the dis- tribution is likely to be under-sampled and cannot have negative velocity values. While the nodal distributions could probably be inrproved by computing additional tomographic reconstructions (using additional and dif- ferent starting velocity models), the increased time involved renders it impractical. Also, the true distri- bution is better represented as truncated, i.e. having a local minimum and maximum; however, the small rel- ative error introduced does not outweigh the compu- tational convenience associated with the assumed normal distribution. Fig. 13a and b show two 3-D perspectives of anom- alous regions (defined by a 95% isoprobability sur- face) and mine related features. The anomalous regions, indicative of concentrated high (red) and low (blue) compressional stress, were identified by assign- ing threshold velocity values. Again, the low stresses follow drifts while the high-stress regions appear to be associated with (typically immediately above) mined features. For the case shown, the maximum (5.25 km/ s) and minimum (3.25 km/s) threshold limits were chosen to be 15% less than the theoretical value for intact quartzite and equal to the lower 25% (quartile) of the observed velocity measurements, respectively. It should be noted that these regions will either increase or decrease in volume depending on the assigned threshold values and/or magnitude of mapped surface probability. For example, when the threshold probabil- ity was decreased from 95 to 67%, the anomalous high- velocity volume regions increased by a factor of roughly four. This volume increase corresponding to probability decrease reflects the lack of starting model independence. To minimize uncertainty and thus max- imize robustness with regards to identification of anom- alous regions, only the 95% isoprobability surface should be mapped. Conversely, the ability to modify high and/or low threshold velocity values affords the latitude required to meet site specific objectives. 8. Conclusions Underground mining provides a unique environment that is amenable to active 3-D seismic tomography. In this case, the use of haulage drifts and crosscuts on multiple levels offers distinct advantages over conven- tional crosshole surveys. For example, the rock bolts installed along drifts provide a solid location to mount geophones and as strike points for introducing seismic energy using a sledge hammer. Further, the drilling and tapping of rock bolts minimize time, cost and logistics, as compared to hard rock drilling for conventional cros- shole studies. Finally, the practice of bolting through- 18 M.J. Friedel et al. /Journal @Applied Geophysics 34 (1995) 1-21 Fig. 13.3-D renderings of high (5.25 km/s threshold) and low (3.25 km/s threshold) velocity regions and related mine features. The high- and low-velocity regions indicate regions of anomalously high and low stress (95% probability), respectively. The silver vein appears as a translucent veil that connects the upper and lower backfilled portions (brown) ; drifts are shown in grey. 20 M.J. Fried& et al. /Journal of Applied Geophysics 34 (1995) l-21 out the mine permits flexibility in expanding or contracting future surveys. During this 3-D active tomographic survey, more than 2,400 seismic waveforms were recorded using multiple source-receiver geometries (24-channel com- mon shot gathers) between drifts on common levels, as well as between levels, using a 3-person team for three lo-hour shifts. The actual amount of time required for a tomographic study depends on a number of factors including the number of seismic records, manpower, number of recording channels, takeout spacing, ambi- ent noise, and other mine related concerns, e.g. drilling, blasting, mucking and hauling of the ore, and rockburst activity. Optimum tomographic survey conditions were achieved while recording during ‘ ‘off’ ’ production time, e.g. during a maintenance shift. From the suite of eleven tomographic reconstruc- tions, basic statistical properties, e.g. average, median, standard deviation, were computed for each node. By averaging the reconstructed velocities, model-depend- ent local variations were effectively filtered, resulting in a smoothed and model-independent velocity field. Globally, these tomograms appeared heterogeneous (velocity range: 1 S-6.0 km/s) and anisotropic (struc- ture was coincident with the dip of bedding). Locally, the velocity structure appears well correlated with mechanical models indicating the transfer of stress in direct response to mining. For example, anomalously high-stress regions indicated above and below sandfil- led stopes using plane strain models were also identifiable as high velocity regions in the tomographic reconstruction. Other localized concentrations of stress were also discemable, as was the extent of stress relief surrounding drifts and at locations of prior rockburst- ing. In general, low-velocity regions followed drift con- figurations extending roughly l/2 to three drift diameters into the remnant structure. Regions of anom- alous ground, i.e. intensely fractured and highly stressed, were identifiable using the concept of thresh- old probabilities and isosurfaces. The identification of anomalous regions using this approach provides a rational means whereby appropriate mitigation strate- gies can be devised and implemented based on the assessment of tomographic uncertainty. While mine-related stress conditions appear to be identifiable through inspection of velocity reconstruc- tions, several limitations on resolution are also evident. For example, it is not possible to completely resolve the stress-relieved zone surrounding outer drifts at this site due to asymmetric sampling, i.e. sampling only along one side. Similarly, it is not possible to directly resolve the sandfilled stopes because of the tendency for seismic energy to seek faster paths. While the seis- mic resolution (about 9 m) may be increased by increasing the low-cut filter threshold, the voxel dimen- sions must be made smaller. A reduction in voxel size, however, requires a corresponding increase in their total number and additional traveltime information to augment the system of equations. In contrast to the modeled stress distribution, the microseismicity did not demonstrate any consistent and observable correlation with the tomographic velocity distribution. One plausible reason may be that some small-scale events occur in response to tensile and not compressional failure. Other contributing reasons may be errors in source location and/or spatial velocity. To avoid these potential problems, a program is now under development that can simultaneously invert both active and passive seismic data as well as introduce curved ray tracing. Other efforts are directed at monitoring temporal variations of the Lucky Friday remnant struc- ture, as well as a larger structure at a depth of 2.2 km at the Homestake gold mine, Lead, SD. An investiga- tion of the statistical relationship between seismic velocity and in-situ stress measurements is also under- way. By combining these two types of information, it may be possible to map stress in three dimensions using a limited number of discrete measurements. Recogniz- ing that refracted shear waves are identifiable in the majority of seismic records, it is also reasonable to conduct the simultaneous inversion of both P- and S- wave travel for a more accurate 3-D distribution of elastic moduli. The use of dynamic moduli derived in this way could provide for more realistic mechanical models. Acknowledgements The authors thank the management of Hecla Mining Company for providing access to both the Lucky Friday mine and to the microseismic database. We also wish to express our gratitude to Terry McMahon and Jerry Hollis for drilling and tapping rock bolts and to Sean Killen for providing the 3-D LF rendering. M.J. Friedel el al. /Journal of Applied Geophysics 34 (1995) l-21 21 References Blake, W. and Cuvelier, D., 1990. Developing reinforcement requi- rements for rockburst conditions at Hecla’s Lucky Friday mine. Proc. 2nd lnt. Symp. Rockbursts and Seismicity in Mines ‘88. Balkema, Rotterdam, pp. 407-409. Dubinski, J. and Dworak, J., 1989. Recognition of the zones of seismic hazard in Po1ir.h coal mines by using a seismic method. In: S.J. Gibowicz (ENditor), Seismicity in Mines. Birkauser, Basel, pp. 609-617 [reprinted from Pure Appl. Geophys., 1291. Friedel, M.J., Tweeton, D.R., M.J. Jackson, Jessop, J.A. and Billing- ton, S., 1992. Mining applications of seismic tomography. Proc. Sot. Explor. Geophys. Annu. hit. Meet. New Orleans, LA, Octo- ber 25-28, 1992. Sot. Explor. Geophys., Tulsa, OK, pp. 58862. Friedel, M.J., Jackson, M.J., Tweeton, D.R. and Olson, MS., 1993. Application of seismic tomography for assessing yield pillar stress conditions. Proc. 12th Int. Conf. Ground Control Min., Morgantown, WV, Aulptst, pp. 292-300. Friedel, M.J., Jackson, M.J. and Williams, E.M., 1994. Tomographic imaging of coal pillar behavior: observations and implications. hit. J. Rock Mech. Min. Sci., in press. Ivansson, S., 1985. A study of methods for tomographic velocity estimation in the presence of low-velocity zones. Geophysics, 50: 969-988. Jackson, M.J., Tweeton, D.R. and Friedel, M.J., 1993. Approaches for optimizing the use of available information in crosshole seis- mic tomographic reconstruction. Proc. GeoTech ‘92, Denver, CO, August 31-September 3, pp. 13&143. Kormendi, A., Bodosky, T.. Hermann, L., Dianiska, L. and Kalman, T., 1986. Seismic measurements for safety in mines. Geophys. Prospect., 34: 1022-1037. Mason, I.M., 1981. Algebraic reconstruction of a 2-D velocity inhomogeneity in the High Hazels seam of Thoresby Colliery. Geophysics, 46: 298-308. Maxwell, SC., and Young, R.P., 1992.3D seismic velocity imaging using microseismic monitoring systems at Strathcona mine and mines Gaspe. Proc. 94th AGM. Can. Inst. Min. McCreary, R., McGaughey, J., Potvin, Y., Ecobichon, D., Hudyma, M., Kanduth, H. and Coulombe, A., 1992. Results from micro- seismic monitoring, conventional instrumentation, and tomog- raphy surveys in the creation and thinning of a burst-prone sill pillar. Pure Appl. Geophys., 139(3/4): 349-373. Scott, D.F., 1993. Geologic investigation near an underhand cut-and- fill stopc, Lucky Friday mine, Mullan, ID. BuMines RI, 9470,25 PP. Stickney, M.C. and Sprenke, K.F., 1992. A microseismic study of the Coeur d’ Alene Mining District, northern Idaho. In: J.R. Tillerson and W.R. Wawersik (Editors), Rock Mechanics. Bal- kema, Rotterdam, pp. 1087-1092. Telford, W.M., Geldart, L.P., Sheriff, R.E. and Keys, D.A., 1984. Applied Geophysics. Cambridge Univ. Press, Cambridge, 860 PP. Thill, R.E., 1989. Coal and rock properties for premine planning and mine design. BuMines IC, 8973, pp. 15-32. Urn, J. and Thurber, C.. 1987. A fast algorithm for two-point seismic ray tracing. Bull. Seismol. Sot. Am., 77( 3): 972-986. Whyatt, J.K., Williams, T.J. and Blake, W., 1993. Concentrations of rockburst activity and in situ stress at the Lucky Friday mine. In: R.P. Young (Editor), Rockbursts and Seismicity in Mines ‘93. Balkema, Rotterdam, pp. 135-139. Wielandt, E., 1987. On the validity of the ray approximation for interpreting delay times, In: G. Nolet (Editor), Seismic Tomog- raphy. Reidel, Dordrecht, pp. 85-98. Williamson, P.R., 1991. A guide to the limits of resolution imposed by scattering in ray tomography. Geophysics, 56( 2): 202-207. Young, R.P. and Hill, J.J., 1985. Seismic characterization of rock masses before and after mine blasting. Proc. 26th US Symp. Rock Mech., Rapid City, SD, pp. 1151-l 158. Young, R.P. and S.C. Maxwell, S.C., 1992. Seismic characterization of a highly stressed rock mass using tomographic imaging and induced seismicity. J. Geophys. Res., 97(B9): 12,361-12,373.