A spherical-stochastic methodology for microseismic event location

Computers & Geo~ciem'es Vol. 15, No. 7, pp. 1037-1052, 1989

0098-3004 89 $3.00 + 000 Pergamon Press plc

Printed in Great Britain

A SPHERICAL-STOCHASTIC M E T H O D O L O G Y FOR MICROSEISMIC EVENT LOCATION Energy Systems. Division 6213. Sandia National Laboratories. Albuquerque, NM 87185. U.S.A (Received 25 November 1987; revised 6 December 1988)

Abstract--The potential economic value of determining the geometry of hydraulic fractures for reservoir engineering is well established. Borehole seismic systems utilizing the polarization method of analysis have been employed successfully in numerous hydraulic fracture mapping experiments. The sensor for the borehole seismic system is a four-axis tool with four or seven geophones per axis that sense the magnitudes of the rectilinear components of seismic signals. From these data. the direction of incidence and distance of the seismic source to the sensor can be determined. This paper describes the methods for analyzing this angular data and for combining the data from more than one tool into a single estimate of the location of the origin of the seismic signal. Results from applying this methodology to the determination of fracture azimuth also are presented. Key rVords: Directional statistics, Fracture diagnostics, Fracture mapping. Microseismic event analysis,

Reservoir simulation, Spherical-stochastic analysis.

INTRODUCTION The National Petroleum Council estimates that in the U.S. there are 2.62 x 10t~m ~ (9.25 x 101eft;) of natural gas in the tight gas basins and 40% of the estimated recoverable tight gas lies in lenticular tight gas reservoirs (National Petroleum Council. 1980). Historically, tight gas fields have been uneconomical and inefficient to produce because of the low natural flow rates of the gas. Massive hydraulic fracturing has been used to increase production, but results have been disappointing. Improved results could be expected from the development of techniques for controlling the growth of the fracture geometry during the fracturing process. Current technology relies upon a mathematical model to predict the fracture geometry as a function of the formation properties and the treatment parameters. Accurate knowledge of formation properties is essential for input to the model (Pai and Garbis, 1983). In exercising the model, the treatment parameters are varied to determine their effect upon fracture geometry. Initially, a set of treatment parameters is selected that will maximize the economic return (i.e. the increase in revenues minus the treatment costs). Once the fracturing process has begun, the real-time data usually are limited (e.g. to fluid flow rate and possibly wellbore pressure and temperature); hence, the dimensions and rate of growth to the fracture cannot be estimated accurately, Without estimates of fracture size, values for the in-situ viscosity and leakoff rate (obtained by comparing fracture-model output with observed fracture geometry) cannot be updated: and, lacking these values, the fracture model cannot provide reliable information for use in controlling or altering the treatment to achieve the desired fracture geometry. A real-time, fracture-

diagnostics instrumentation system is needed to control treatment and to determine whether the ongoing treatment, in fact, is appropriate. If real-time information describing fracture geometry were available for feedback to the fractural model, then continual comparison of the predicted geometry from the model with the geometry sensed by the instrumentation system would make it possible to provide updated values for the in-situ viscosity and leakoff rate. Provided the fracture model is correct and convergence is obtained in this model-instrumentation loop, then fracture-treatment modifications can be tested to determine their effectiveness in controlling the growth of fracture geometry. Using the desired fracture geometry as a goal, iterations within the treatment parameter loop potentially could be useful for providing the information needed to guide and control the course of the treatment. It is possible that the iterative procedures will not lead to eventual agreement between observed and calculated fracture geometries. In this event, there is no modelderived rationale by which the ongoing treatment can be controlled or improved, and empirical methods must be used in determining the course of the treatment. Figure '1 is a flowdiagram that shows how a real-time instrumentation system might be used to optimize the fracture treatment. An integrated, real-time instrumentation system consisting of a number of diagnostic subsystems is being developed (Engi, 1983). A conceptual framework for the system is shown in Figure 2. To provide this real-time instrumentation system, four diagnostic subsystems are being developed currently. These subsystems, which are to be linked together by a microcomputer network, consist of (I) an internalfracture-pressure measurement system, (2) a fluidflow measurement system, (3) a borehole seismic sys-

1037

PREEIiAC SURFACE OIE|RVATIONE

ACTUAL FRACTURE GEOMETRY

CO~E MEASUmEMENT| I

COHE ERACTURES

,o___

t

I

DATA ACGUIEITION 1 ~ EYSTEM

MOOELING

l

Ile.=-.w.HIt

[~,.u.

FRACTURE MODEL

I*-r

REAL TOME

POITFRAC

'TREES t MEASUREMENTS

I,;~.*.~1,1 ME~;~;;ON \1 AND M D D E L ~ aonE~MENT

t ,#.E. •

f

I 'R'CtURE ,ALIDAT,O MODEL AND OTHER MESERVOIR 0ATA

NO F.m.plSANOi

I i

upOAte "MODEL

I SHOULD THIS RESEmVOlli i~1-. ] lie IIROOUCEO HOW?

I r::':'::°=°*':"

},']

'o,=;r,'::~ r i*¢'1,;:~;: ° I:i~,.Ire l, RCaI-Umc mslrumcriLLli~m ll~)v, di;IFr:lm

~

l

U

R

P SSS SIGNAL

SEP SIGNAL PROCESSING

PROCESSING

#WAm

FRACTURE PARAMETERE: l " FRACTURE ORIENTATION LuL J • MAJOR FRACTURE WING LENGTH LUIN • MINOR FRACTURE WING LENGTH H • FRACTURE HEIGHT W • FRACTURE WIOTN u " LmA~/LU,N 8 " FLUID LEAKOFF RATE

= H

4V(1 -/J) H

LMAo ~L LM,n (I -.)L (8, Lu•~ Lu, w H)N|Am

E

~ -y-

LuAj. LuI~lrAm

WIDTH EQUATION

(L. H. WlNiAn • FAro

10, LMi J . LugN , H, W)NIAI~• PAIr

Figure 2. Real-time rracture-diagnostics instrumentation system

(P)

/

r

Methodolog? for microseismlce',ent location tem (BSS). and (4) a surface-electric-potential (SEP) measurement system. This paper is de,.oted to a description of the BSS and the analysis of the data acquired from a BSS. DESCRIPTION OF THE BOREHOLE SEISMIC SYSTEM The borehole seismic system (BSS) that is being developed by Sandia National Laboratories is capable of detecting, recording, and analyzing microseismic signals generated during the hydraulic fracturing of an oil or gas reservoir (Thorne and Morris. 1988). The sensor for the BSS is a four-axis tool, on each axis of which are mounted four-seven geophones. The tool is contained in a cylindrical housing, which protects the instruments from downhole pressures and which is clamped tightly against the wellbore wall to provide good acoustic coupling between the geologic formation and the gcophones. The depth of the tool below the surface is determined by monitoring the length of cable rcclcd out while lowering the instrument into the borchole. The directional orientation of the tool inside the borcholc can be determined by ,tnalyzing the signals detected by the tool from calibration shots detonated tit known localities with respect It) the horchole. TIll': %'EI,OCrI'Y FA('I'OR |'ach seismic signal consists of two elastic pulses tra~,eling at different speeds outward from the source v, ith approximately sphcric:il wavefronts. For the faster P-wave (primary, coniprcssional), the particle nit)tions tire in the direction of advance of the wave. For the slower S-wave (secondary, shear), the particle motions are orthogomd to the direction of advance of the wave (see Fig. 3). The wave velocities, % and v , , of the P-wave and S-wave, respectively, are determined from calibration shots. The ~,elocity factor then is calculated from: Ii

(%

v, % -

(1) v,)

Alter l; has been determined, the distance, d, from the origin of a microseismic event to the tool can be estimated from: d

SHEAR WAVE \

=

~;(t,

-

SEISMIC SOURCE DIRECTION The direction of a seismic source can be described by the spherical coordinates 0 and • (see Fig. 4). During the period between arrival of the P-wave and arrival of the S-wave, a tool placed to intercept the wave will resolve the seismic signal into its rectilinear components: X,

=

k.4, sin~, cos O,

Y,

=

k.4, sin ~, sin0,

Z, =

kA, cos~,

(3)

~here k

= proportionality constant × (dimensionless)

A, = signal amplitude (V) at time slice i. From Equation (31. the expressions for 0. and ~, are: O, =

tan i ( y , ! x , ) +_ n

(4)

~ll =

tan i[()t," + }]")):!:,] + ~.

(5)

The ambiguity of ~ radians in the values for 0 and (i) arises from the inability of the tool to distinguish between the comprcssional phase of a wave incident from one direction and the dilutational phase of a wave incident from the opposite direction. Because of this ambiguity, 0 and (ll are described more accurately as axes (rather than directions). When data from more than one calibration shot are avaihtble, the ambiguity can be resolved. If no noise or distortion is present in the signal then till the O, (and (I),) tire equal and the direction of the seismic source can be determined from Equations (4) and (5). ANALYSIS OF DATA FRONt A SINGLE TOOL The seismic events resulting from the hydraulic fracturing process are thought to be caused by shear failure induced by the localized high-pore-pressure zone surrounding the fracture, by the presence of inhomogeneities near the fracture, and by the fracture itself (Dobecki, 1983). The locations of these seismic events are believed to cluster in a tight band near the z

(2)

t~).

SEISMIC SOURCE

1039

PARTICLE MOTION

I

Figure 3. Seismic wave types.

Figure 4. Rectilinear and spherical coordinate systems.

1040

E~;GI

D

FORMATION ANISOTROPHY AND - ' ~ TOOL ' " ~ ERRORS

':,72;7""'.\

NOISE

~

l.cou,.,c,... U,...o...o..o..,o.I [.,,..,..o i-i...oo°.,o-.,o. I -I

ACOUSTIC WAVE INCIDENT ON GEOPHONES

I I I I/ TRANSFORMATIONTO I I " - - I COMMON COOROINATE ~ I/ SYSTEM |li~ It It')_

' ESTIMATE OF LOCATION (~, y, x)

/

A

i

J

I I I I I

~ ANALOG SIGNALS | - i PROOUCED BY ~| GEOPHONES | (X. V, Z)

i

, .o.......o..0., IX, ,, ¥) [~ t

L:

I

)\

I

I ANALOG SIGNAL I

~ ~ a I / ~ ' ~ POLARIZATION PLOT W

IlL ~

O'omzeo (X~, V~,Z)}

I |,

I i |

Figure 5. Error introduction into seismic signal.

hydraulic face and arc used to map the fracture (AIbright and Pearson. 1980). Errors arise from many sources; for example, the tool itself can introduce noise, extraneous seismic waves present in the I\)rmation cvn mix with the source signal, .'misotropics and inhomogeneitics in the formation can distort the signal, and signal processing can introduce amplitude changes and phase shifts between signal components. Figure 5 illustrates some of the points at which errors can be introduced. Because of the error and distortion, 0 and 4) are regarded as random variables defined within a limited range and are handled probabilistically, It is tempting to treat the data in the conventional linear fashion to determine the mean and variance for 4) and 0; however, this would be inappropriate because the mean and variance are affected by a change in the zero direction. A different method of analysis is required. The method used here is based upon directional statistics (Mardia, 1972). Each data point is considered as representing a vector, ~, extending from the origin (0,0.0) centered on the tool. to the point (X,. Y,, and ,7..,). Each vector is described by its direction cosines l,. m,, and n,. which are the cosines of the angles ~t,, fl,. 2, respectively, as shown in Figure 6. The mean direction is defined as the direction of the resultant, p, which is the vector sum of the p,. Clearly. the direction cosines of p are given by: /o

=

cos20 =

mo = cos~o no

=

=

COS 70 =

5"x,,p ~r,/p [Z,,"p

written

4).

=

a.:t~,n [(~.v,): + ( ~ r , ) : l ' :

The mean resuhant, l',, is p.

detined

=

as:

I, ~I,,.

(8)

If all the fi, have the same direction then I'(, = 1. otherwise ¢,, < I. Using the translbrtnation Po = exp[-(Var(4))):/2], the variance of 4) becomes: Var(4)) =

- 2 In(p,)).

(9)

Expressions for the 1st and 2nd moments of the azimuth can be obtained by suppressing the 7., components and treating the data as two dimensional. The results are: 0o =

arctan

Var(0)

=

[ZY,/~X,]

(10)

-2In(r0)

(11)

where

ro

[(Ix,): + ( [ r , ) : ] ' : =

~(X,:

+

(~2)

Y,:)J :

,%

(6)

S

where

p = Noting that ;'

-

(7)

y z,

magnitude of p. 4). the mean value of 4) can be

Figure 6. Direction angles for A.

;y

1041

Methodology for microseismic event location POL-~RIZATION PLOTS

independent random variables are given by:

Polarization plots may' be used to assist in the analysis o f d a t a (Gal'perin, 1984). If)(, is plotted vs Y. the line from the origin to the point (X,, Y, ) defines the angle ~,. For the vertical angle the abscissa and ordinate are given by: abscissa =

e(x)

=

y[~:E(x,)]

Var(x)

=

'~.[:t~Var(.s) ].

The expression for :t< that yields the minimum value for Var(x) whereas observing the constraint is: I

(X," + Y,:)t: cos (0, - 0,,) % =

ordinate

Var(%) V Var(.r, ).

in analyzing data. it is essential to remember that the geophone data specify axes. not directions. Successful data analysis requires that the ambiguity in direction be eliminated temporarily. This can be done by assuming that all seismic signals arrive at the geophone from the positive-x and the positive-z direction. that is set Z, = abs (Z,) and X, = absl.'t;). The mean directions iO and O) then can be computed, but they arc ambiguous by _+.rt radians. In thrcc-dimensional space, there v, ill bc four directions that arc consistent with the two possible values of 0. and the two possible values of O,,. Additiomtl information is needed to resolve the ambiguity, if data for the same microscism are available from two (or more) tools. then the true direction can be resolved. ('OMBININ(; I)ATA F'RO.M MORE TIIAN ONE TOOl. Each tool Ioc:ttcs the origin of the seismic signal at point (p, 0, O),. (The subscriptj will be used to denote the tool. as opposed to the subscript i, which denotes the time-slice.) Because of signal errors, the points (p. 0,0), will not be the same. To combine the measure from each tool into a single estimate of the location of the seismic source, it is necessary to transform to a common Cartesian coordinate system. The transformation equations are:

z,

p, cosO, + q

=

The derivation of Equation (16) is given in Engi (1987). The equations for)' and : are identical to those f o r .V.

For a two-tool system, the expression for Var(x) is a family of ellipses. The axes of each ellipse are proportional to l,rVar(x~ ) and I Var(x:), respectively. In three-dimensional space, with :g and ztz as the horizontal axes and Var(x) as the vertical axis, the filmily of ellipses is represented as an elliptical cone (Fig. 7). The constraint is represented by a vertical plane through the points (I,0,0) and (0,l,0). The intersection of the plane and cone identifies all ellipses satisfying the constraint. The lowest point [minimum Var(x)] on the curve of intersection is the desired solution. Similar equations for y and : can be written. Because :t, is a function of the variance, the Ist and 2nd moments of the coordinates (x,.v,z)~ are required to determine the optical location estimates. Using Equatton (13), and assuming independence of p, 0, and q~. the lirst moment can be expressed as: E(x,)

=

E(,g) E(sinO~) E(cos0,).

(17)

The second moment can be expressed as: Var(.V~) =

E(t,~)E(sin:O~)E(cos:O~) - [E(p,)E(sinO,)E(cosO~)]:.

(18)

The problem can be stated concisely as follows: Determine:

% = p, s i n ~ cos 0~ + a, p, sinO, sin0, + b,

(16) I

= Z,.

.v, =

(15)

E(h u) =

(13)

~ Wf(~)d¢

(19)

-I

where (a~. b,. q ) are the coordinates of tool j in the common Cartesian system. This transformation assumes that the zero direction for all the ~j is vertical and that the zero direction for all the 0j is the same geographical azimuth. It is assumed that the coordinates (x,y,:) of the true location of the seismic source are linear (convex) combinations of the corresponding measures from each tool, having the form: x =

~:t,vj, with the constraint ~:tj =

1.

(14)

The constraint assures that x is a weighted average of the %s. For the discrete situation, the 1st and 2nd moments of a function that is a linear combination of

and E(W:)

=

*f" hu-'f(,~)d,~

(20)

-at

where f(O

=

the probability density function for

q' ~ [pj, sin0j, cos0 r sin%, cosOj], j ~ [1,2] ( ~ [p,, 0 , Oj].j ~ [I,2]. The general solutions of Equations (19) and (20) cannot be obtained. A closed form solution can be obtained if it is assumed that p,. O,. and 0, each are

1(~-12

D. Exo! Var(x)

ELLIPTICAL CONE

al2

+

Var(x)/Var(x+)

a2~

=1

Var(x)/Var(x 2)

JRVE OF INTERSECTION E L L I P T I C A L CONE AND VERTICAL PLANE

M I N I M U M Var(x)

VERTICAL PLANE ~rI + a' 2 ~ 1

(~, o, o) J,~

oi

":',~L(o,1, o)

~

.0. 2

Figure 7. Geometrical representation of selection of minimum Vat(x).

distributed triangularly, Because the triangular distribution delivers a probability density function that is linear in form, solutions to Equations (19) and (20) can be obtained. For a triangular distribution the probability distribution function is defined as

are known. ["or 0 and O, the mean and variance can Ix: calculated as shown. The mode for 0 (or ~) will bc taken as the value of 0, (t~,) for the data point having the maximum signal strength. The maximum and minimum can be calculated from:

C, f(.:)

= 2(.:. - 51.~)[(.:. - ¢,)(;. = O.

- ,~.)] z. il',~ ~< ~.

. . . .i

~,

/.") ¢ -l-

-

t-'~,,, -

~h +

3(,;~ -

(,,):]l'}

[24 Var(~) (22)

otherwise ,;~

(21) v.here ,;,,. ,;~,. and ,;, are the minimum, mode. and maximum, respectively, of the triangular distribution. The values for the minimum and maximum can be determined provided the mean. variance, and mode

=

3~. -

(~h-

~.1.

(23)

The derivations of Equations (22) and (23) are given in Engi (1987). To determine the minimum, mode. and maximum o f d . the minimum, mode, and maximum oft, and t,

.'Methodolog.'+ for microseismac event location (the arrival times of the S-reave and P-wa,,e. respectively} first must be determined. In general, both t, and t~ can be specified as falling within a range of values spanning tmo consecutive cycles of the digitizing frequency v.ith the range of t, occurring later than the range of t r. It is assumed that the probability of occurrence of t. is uniform ~,ithin the range of t.. that is a uniform distribution: and similarh for tp. Because the digitizing frequency is constant for the entire seismic signal, the magnitude of the range of t, is equal to that of t+,. Because t, and t~ are independent+ continuous, uniformly distributed random variables, the probability densit) function of their difference is given b,,:

Ii

1

;'

z[(O, - f~)cos20~

2 (0 - Ohlcos20~ - (0~ - 0~)cos20 ]

-

(30)

E(cos-" 0) =

I

~[(0

-

O.)cos2tL

- (0 - 0,)cos20,, - (0~ - 0.}cos20,] (31) where 7 =

At/b

/v), 1. - lp ~< A ~< h. - h r (Region II)

K(h,

A - I~,). h,

O,

2[(0 - 0,)(0 - 0~)(0~ -

0,,)] -~.

(32)

othcrv,isc

i~

=

I,

--

1.

~t,~h,

Ip

]l, <- t l, ~ K

hp

I

=

(h r, - l~)(h, - L ) The derivation ol+g(A) is given in Engi (1987). For the situation under consideration, (h r - /+,) -- (h, - 1,) and Region II vanishes: :rod g(A) is a triangular distribution. The minimum, mode, and maximum of (t, tr,) for this triangular distribution arc: Mhfimurn (t, - tp)

=

l, - h r

M o d e ( t , - It,)

=

h, - hi, =

M~,ximum (t, - t+,) =

l, - 1~,

h, - Ip.

(25)

The corresponding values of d then are computed from Equation (2). With the minimum, mode, and maximum values for the triangular distributions now known, Equations (19) and (20) can be integrated to obtain the required expressions for E ( ~ ) and E(Ud z). The results are as follows:

E(d)

= ~(+1. + ,I,, +

~'[,l~(d,

E(,IZ) =

d,)

(26)

- d+) - d,~(d, - dh)

3 -

=

d, (d,. -

=

d~)]

(27)

7 [ ( 0 , - 0,,)sin0b - (0, - 0,,)sin0,, -

(0~ - 0+)sin0,]

(28)

7[(0, - 0.)cos0~ - (0,. - 0Dcos0. -

(24)

lip <~ A <~ tt, - Ir (Region III)

here

E(cosO)

=

K(hp - I, + L%).1. - /h, ~< A ~< l, - lp (Region I)

g(A)

E(sin0)

E(sin:0)

1043

(0,. -

0,,)cos0,]

(29)

The subscript j hits been suppressed. The equations for • are identical to those for 0. The derivations of Equations (26)-(32) are given in Engi (1987). Two BASIC programs have been written ( I B M BASI('. Version 3.2) to perform these computations on an IBM PC AT. The first. SPIISTAT.BAS. takes the digitalized amplitude data from the geophones as input and calculates the average and the variance fi~r the ~,nglcs 0 :rod O. It also computes the nfinimum. mode, and maxinwm t\~r the assunled triangular distribution for both angles. The second program. [X)CATE.BAS, takes the nfininlum, mode. and m a y imunl ft,- d. 0, and qJ. respcctb+ely, as input par+,meters and computes the rectilinear coordinates .+,,)., ;.ind 2. This program also will combine the angular tk, ta from two tools into :L single estimate of the seismic source location. The program listings are in Appendices I and 2 respectively. A sample problem is provided in Appendix 3. APPI+ICATION OF TIlE METIIOI)OI.OG~'

The event location gi,,en by the (x..v. :) coordinates can be interpreted as an ellipsoid centered at the estimated location of the seismic source and with its ~,xes parallel to the axes of the common coordinate system. The length of each major axis is proportional to the standard de~,iation of the location estimate for that axis. The volume of the ellipsoid is related directly to the uncertainty of loc:,tion. The location uncertainty for each individual tool will be greater than the uncertainty for the minimum-error event location in which the data from both tools have been combined into a single estimate. Figure 8 illustrates the location method as applied to a single microseismic event resulting from a stimulation in the Fluvial B sandstone of the Muhiwell Experiment Site. Rifle. Colorado (Thorne and Morris. 1988). In Figure 8. the horizontal locations deter-

D. l::.~ot

1044

lO0

liillllilllll tl6VOlllll'i'~ 6v~OallO mlv~ VION , l • I i vw.-]

)sO ~

I

i

lSO

•

]

60 1

{

O;

0-

-60

.,w.-I ~ l

-)00 L -1SO

-260

-~00

-t~O

-)00

0

-~0

$0

lO0

-iO0

*tSO

t -iO0

i

-50

WEST - EAST

0

SO

)00

ISO

(it)

W E S T - E A S T (ft)

Figure 8 Horizontal projection of microseismic event location.

Figure 10, Horizontal projection of microseismic event It)cation for stimulation experiment in Fluvial B sandstone,

mined from BSS data taken in wells MWX-2 and MWX-3 are shown by the triangles nearest to the respective observation wells. The ellipse drawn about each of these two locations represents the horizontal projection of the error ellipsoid for that location. The third triangle and the ellipse drawn about it represent the two-tool minimum-error horizontal location ;ind the horizontal projection of the minimum-error ellip. sold. The line through the stimulation well, MWX-I, is given by the spherical average of the three angles determined by these three locations and is at an azimuth of N70°W. Figure 9 gives the projection of these same three points into the vertical plane determined by this line. Here, the positions of the bottoms of the BSS tools in MWX-2 and MWX-3 are indicated by the square symbols on the projections of the wellbores onto this plane. Depth is measured from the surface at well MWX.I, and the top and bottom of the Fluvial B zone are indicated by horizontal lines. The lower triangle corresponds to the location determined from MWX-2 and the ellipse about it to the projection of the error ellipsoid for that location onto the vertical plane. The upper triangle and ellipse correspond to the location determined from MWX-3; the third

triangle and its associated ellipse correspond to the minimum-error location and the minimum-error ellipsoid. Note that in both projections the ellipse about the minimum-error position is smaller than the ellipse about either single-tool location. Figure 10 shows the locations given to 29 microseismic events which were located from both tools. All of these locations lie between parallel lines 25ft (7.6m) on either side of this line, consistent with the scatter in locations determined for the perlbration (vek)city factor determining) shots in well MWX-I. Thus, reasonable confidence can be had in the N68°W fracture azimuth as determined from these locations, even though the spherical statistics on the set of angles measured from well MWX-I to event locations yield a 21.3 ° sd. Moreover, the N68°W azimuth is in excellent agreement with previous measurements at this site (Teufel and others, 1984).

*~ObO

~llwx.l I

I

ll*;I-i

lwl-~

-6700

-ITSO

4O00

I

-$11~

-

0

I'%-

-6100

-$O60 -llO

'

"|~

Li

o1~

I *100

SPHERICAL

I

i

40

MEAN

~)O

OIRECTION

I tOQ

IIO

(it)

Figure 9. Vertical projection ol'microseismic event location.

SUMMARY

In the procedures shown here, the techniques for analyzing angular data have been applied to seismic data acquired from orthogonally oriented geophones. A simple weighting factor has been introduced to weight most heavily data with the largest signal-tonoise ratio. The limitations of the geophone in distinguishing directions (as opposed to axes) have been accounted for, but not entirely overcome. Additional information, beyond just the seismic data, is needed to identify a final direction. The distortions arising from random noise and from phase shifts between signals from different geophones are treated as random error. The difference in velocity between compressional waves and shear whves is used to obtain a distance estimate that is then combined with the directional estimate to form a completed location estimate. Data from more than one tool can be combined into a single location estimate with an associated uncertainly than that achieved by a single tool.

Methodology for microseismic e'.ent location NOMENCLATURE

REFERENCES Albnght. J N.. and Pearson. C. F., 1980, Location of h>draulic fractures using microseismic techniques: So<. Pet. Eng AIME. Paper SPE9509. 10p. DobeckL T. L.. 1983. Hydraulic fracture orientation using passl'.e borehote seismics: So¢. Pet. Eng. AIME, Paper SPEI2110, 6p. Engi. D . 1983. Integrated real-time fracture diagnostics instrumentation system: Western Gas Sands Subprogram Review Technical Proceedings. D O E METC 84-3, p, 195-199. Engi. D.. 1~87. A numerical signal processing technique for borehole fracture diagnostics: SAND86-1865 (Albuquerque: Sandia National Laboratories}. p. AI-A2. BI-B2, CI-C4. D I - D 6 Gal'pcrin. E. I.. 19~..I., The polarization method of seismic exploration: Klav.er Academic Publishers. Hingham. Ma,,,,achusetts. p. 5, 12-14.20-21. Mardia. K. V. 1972. Statiqics ofdirectional d~lta: Academic Prcs'~, Inc.. Orlando. Florida. p. 19-28. 74. National Petroleum Council. 1980, Unconventional gas ,~ource - tight gas reser,.oirs: Library of Congress Card No SO-82488. p. 5-11. P~,. \'..and Garbis. S.. 1083. Tight sands require detailed planmng: Oil and Gas Jour.. v. 81, no, 30. p. 135-139. Teul~'l, L. W.. Ilart. C. M.. Saltier. A. R.. and Cl,.rk, J . A . . 1'4,"¢4. Determination of hydraulic fracture azimuth by gcoph~,.ical, geological and orientated core methods at the multtv.ell site. Rifle. Colorado: So¢. Pet. Eng. AIME. Paper Sl'|!1.1226. I I p. Thorne, B J.. ;.lnd Morris, I I E., 1988. Passive sc~isnlic nlt'qlilormg o f h,~dr,tulic fr,tclurc expcrimcnt~, ;it the rnultiv, cll cx~.rlmcnt site: SANI)X~-1284 (Albuquerque: Sandia Nat,real I,abor:ttortcs), p, 14.17. 32-34

Progrant

1190 1200 1210 1220 1230 1240 1250 1200 1270 1280 1290 1300 1310 1320 1330 1340 1350 1360

Listing

,4(t~ d E(.,c)

= Signal amplitude = Distance from seismic source to tool = The expected value of .x. defined as EI.~I

f(:c) I.m.n r t

= The density function of x = Direction cosines = Radial distance, r: = ~: = Time = \'elocit~ factor = \'clocit'.

=

l)

t' \'ar(r~

["~

xf(xldx

~:

= The ~ariance of ~. defined as: Var{.~:) =

i [x - E(~)]" t i ~ ~ d . x

X.Y. Z = Signal amplitude components in the r. v. _directions x..r. : = Cartesian coordinates :t,/I. ;' = Direction angles 2 = Langrangian multiplier e.~.tl = Spherical coordinates q., = A generalized coordimtte = A generalized coordinate a = Standard deviation ~ = Weighting factor Suh~xril, ts

a h c t /

Minimum - Mudc - Maximum = Time slice idcntiticr, t = I. 2 . . . "= l'ool identiticr. / = I. 2 . . . . v = hliti:tl ',~11tie = Prml,lry wave = SIle~t r ',~,a'~ c ~|Cdll

0

p s

APPENDIX

I000 I010 I020 1030 1040 1050 I060 1070 1080 1090 I I00 III0 1120 1130 1140 II 50 1160 1170 I180

10,15

!

of S P I I S T A

T. B..I S

REM SPIISTAT REM This program computes the directional coordinates of a microseism. REM The input data are the digitized amplitude components of the signals RKM generated by orthogonally oriented geophones. REM The variables used are: REM For the computation of the statisticsof TIIETA : REM X(1),Y(1),Z(1)-the amplitude components or the seismic signal REM for the l-th time slice. REM R(1)- SQR(X(1)A2 ÷ Y(1)^2) REM S U M R - the sum of the R(1)'s. REM THETA(1)- the angle defined by ATN(Y(1) / X(1)) REM SUMXSTARthe sum of the x components in the expanded domain. REM SUMYSTARthe sum of the y components in the expanded domain. REM RBARHthe mean resultant in the x-y plane. REM MODETHETAthe THI:TA(I) associated with the m a x i m u m R(1) REM VARTHETA,, the variance of T H E T A . REM AVETSTARthe average of T H E T A in the expanded domain. REM AVETHETA- the average of THETA. REM MAXTHETAthe m a x i m u m of the triangular distribution of T H E T A . REM MINTHETA- the minimum of the triangular distribution of TI-IETA. REM For the computation of the statistics of PHI : REM RHO(I),,SQR(R(I)'2 + Z(I)'2) REM PHI(l)- the angle between RBARH and RHO(1). REM S U M X - the sum of the X(1). REM SUMY- the sum of the Y(1). REM SUMZ- the sum of the Z(I). REM SUMRHO,, the sum of the RHO(I). REM RBARV- the mean resultant in 3-dimensional space. REM MODEPHI- the PHI(I) associated with the maximum RHO(I). REM VARPHI- the variance of PHI. REM AVEPHI- the average of PHI. REM MAXPHI- the maximum of the triangular distribution of PHI. REM MINPHI- the minimum of the triangular distribution of PHI. REM D X and D Y - d u m m y variables. CLEAR:CLS KEY OFF PI-3.141593:K- 180/Pl

n

1040

D. E.~Gi 1370 REM 1380 REM Determine the type of data to be input. 1390 REM 1400 PRINT "This program will accept any o[ three types of input" 1410 PRINT " I. X and Y data input : THETA output." 1420 PRINT " 2, X,Y, and Z input : THETA and PHI output." 1430 PRINT " 3, R and Z input ; PHI output." 1440 PRINT " Choose I, 2. or 3." 1450 INPUT OPT 1460 IF OPT3 THEN CLS:GOTO 1400 1470 INPUT "How many datapoints";N 1480 INPUT "Do you want a hard copy printout (Y OR N)';B$ 1490 DIM X(NI,Y(N).Z(N),RHO~N),R(NLTHETA(N),PHI(N) 1500 SCREEN 3 1510 REM 1520 REM Depending on the input options selected, route program execution 1530 REM to the appropriate subroutines for data input and computation. 1540 REM 1550 IF OPT=I THEN GOSUB 1770 1560 IF OPT=2 THEN GOSUB 1880 1570 IF OPT=3 THEN GOSUB 2000 1580 1F B$='Y* AND OPT=I THEN GOSUB 3000 ISg0 IF B$='Y" AND OPT=2 THEN GOSUB 3130 1600 IF B$='Y" AND OPT-3 THEN GOSUB 3280 1610 CH2$='Y" 1620 LOCATE I,I:PRINT "Do you wish to alter some data and rerun.'? (Y or N)" 1630 LOCATE 1.63:INPUT CH25 1640 WHILE CH2$-'Y" 1650 GOSUB 3410 1660 IF OPT-I TtlEN GOSUB 2110 1670 IF OPT-2 TItEN GOSUB 2110 1680 IF OPT=3 TtlEN GOSUB 2570 1690 IF B$-'Y" AND OPT-I TtlEN GOSUB 3000 700 IF B$-'Y" AND OPT-2 TIIEN GOSUB 3130 710 IF B$-"Y" AND OPT-3 TIIEN GOSUB 3280 720 LOCATE I.I:PRINT "Do you wish to alter some data and rerun? (Y or N)" 730 LOCATE 1.63:lNPUT CII25 740 WEND 750 LOCATE I,I:CLS:KEY ON 760 END 770 REM 780 RI'M This subroutine accepts X,Y input data. 790 REM 800 FOR I-1 TO N 810 PRINT "For datapoint number "1 1820 INPUT "X-'.X(I) 1830 INPUT "Y-".Y(I) 1840 NEXT I 1850 GOSUB 3410 1860 GOSUB 2110 1870 RETURN 1880 REM 1890 REM This subroutine accepts X.Y,Z input data. 1900 REM 1910 FOR I=I TO N 1920 PRINT "For datapoint number "I 1930 INPUT " X - ' , X ( I ) 1940 INPUT "Y=',Y(I) 1950 INPUT "Z-'.Z(I) 1960 NEXT I 1970 GOSUB 3410 1980 GOSUB 2110 1990 RETURN 2000 REM 2010 REM This subroutine accepts R,Z input data. 2020 REM 2030 FOR 1=1 TO N 2040 PRINT "For datapolnt number "1 20,50 INPUT "R,.',R(I) 2060 INPUT "Z-',Z(I) 2070 NEXT 1 2080 GOSUB 3410 2090 GOSUB 2570 2100 RETURN 2110 REM 2120 REM This subroutine computes the mean and variance of THETA. 2130 REM From the values of the mean and variance, it then computes 2140 REM the minimum, mode, and maximum for the triangular distribution 2150 REM of THETA. 2160 REM 2170 SUMR-0:SUMXSTAR=0:SUMYSTAR.,0 2180 FOR i - I TO N 2190 R(1)-SQR(X(I)* X(I)÷Y(I)* Y(I)) 2200 IF R(I)>MAXR THEN MAXIT-I 2210 IF R(I):~MAXR THEN MAXR-R(I)

Methodology for microseismic e~ent location 2220 2230 2240

SUMR-SUMR+R(1) THETA(I)-ATN(Y(I)/X(I)) SUMXSTAR-SUMXSTAR+R(I)*COS(?*THETA([)) 2250 SUMYSTAR-SUMYSTA R+R(1)*SIN(2eTHETA(I)) 2260 NEXT I 2"70 R BA RH-SQR(SUMXSTA R*SUMXSTA R+SU MYSTA R*SUMYSTA R)/SU M R 2280 MODETHETA-ATN(Y(MAXIT)/X(MAXIT)) 2-'-90 VA RTHETA--LOG(RBARH)/2 2300 A VETSTA R-ATN(SUM YSTA R/SUM XSTA R) 2310 IF SUMXSTAR>0 AND SUMYSTAR>0 THEN AVETHETA-AVETSTAR/2 2320 IF SUMXSTAR>0 AND SUMYSTAR<0 THEN AVETHETA-AVETSTAR/2 2330 IF SUMXSTAR<0 AND SUMYSTAR>0 THEN AVETHETA~(AVETSTAR÷PI)/2 23,10 IF SUMXSTAR<0 AND SUMYSTAR<0 THEN AVETHETA.~(AVETSTAR-PI)/2 2350 IF MODETHETA-AVETHETA>P[/2 THEN MODETHETA=MODETHETA-PI 2360 IF MODETHETA-AVETHETA<-PI/2 THEN MODETHETA=MODETHETA+PI 2370 MA XTHETA=.5*(3*A V ETHETA- MODETHETA+SQR(3"(8* VA RTH ETA-(MODETH ETA AVETHETA)'2))) 2380 MINTHETA-3*AVETHETA- MODETHETA-MAXTHETA 2390 CLS 2400 LOCATE 2,1 2410 PRINT TAB(12) "AVERAGE" TAB(22) "VARIANCE" TAB(32) "MINI MUM" TAB(42) "MODE" TAB(52) "MAXIMUM" 7'420 PRINT "THETA "; 2430 PRINT USING "+a~at.a~ ";AVETHETA*K.SQR(VARTttETA)*K.MINTHETA*K. MODETHETA*K.MAX] HETA*K ?440 LINE (I.IJ7)-(201.347)..8 2450 LINE ([.247)-(201.247) 2460 LINE (101.147)-(101.347) 2470 LOCATE 18.23:PRINT'X" 2480 LOCATE 11.12:PRINT "Y" 2490 FOR i-2 TO N 2500 LINE (100*X(I-I)/MAXR+I01.-100"Y(I-I)/MAXR~-247)-(100*X(I)/MAXR+I01.100*Y(I)/MAXR+247) 2510 NEXT I 2520 A- 100*CO~ A V ETII ETA):B- 100*SIN(A V ETiIETA ) 2530 LINE (101.247)-(A+I01.-B÷247) 25.10 CIRCLE (A÷I01.-B~-247).3 2550 IF OPT-2 TIIEN GOSU[/ 2570 2560 RETURN 2570 REM 2580 REM This subroutine computes the mean and variance of Pill. 2500 REM From the values of the mean and variance, it then computes 2600 REM the minimum, mode. and maximum for the triangular distribution 2610 REM of Pill. 2620 REM 2630 FOR [-I TO N 2640 SU MX-0:SUMYo0:SUMZ=0:SUMR HO-0 2650 R HO(I)-SQR(R(I )e R(1)+Z(I)e Z([)) 2660 IF RHO(I)>MAXRtlO THEN MAXIP-I 2670 IF RHO(I)>MAXRIiO THEN MAXRIlO-RHO([) 2680 SUM R HO-SUMR tlO+R HO(i) 2690 PHI(I)oABS(ATN(R(I)*COS(THETA(I)-AVETHETA)/Z(I))) 2700 DX-X(I):DV-Y(I) 2710 IF X(I)<0 TtlEN DX--DX 2720 IF X([)<0 TIIEN DV--DY 2730 SUMX-SUMX+DX 2740 SUMY-SUMV+DV 2750 SUMZ-SUMZ+ABS(Z(I)) 2760 NEXT I 2770 RBA RVoSQR(SUMX*SUM X÷SUMY%UMY+SUMZ*SUMZ)/SU M R HO 2780 MODEPHI-ATN(R(MA XIP)*COS(THETA(MA XIP)-A V ETH ETA }/A BS(Z(MA XIP))) 2790 VARPHI--LOG(RBARV)/8 2800 A V EPHIoATN(SQR(SUMX*SUMX,SUMV*SUMY)/SUMZ) 2810 MAXPHI-.5*(3*AVEPHI-MODEPHI÷SQR(3e(g*VARPHI-(MODEPHI-AVEPHI)'2))) 2820 MINPHI-3"AV EPHI-MODEPHI-MA XPHI 2830 LOCATE 2.1 2840 IF OPT-3 THEN PRINT TAB(12) "AVERAGE" TAB(22) "VARIANCE" TAB(32) "MINIMUM" TAB(42) "MODE" TAB(52) "MAXIMUM" 2850 IF OPT-2 THEN PRINT:PRINT 2860 PRINT "PHI "' 2870 PRINT USING'÷###.#ak ";AVEPHI K.SQR(VARPHI) K.MINPH[ K.MODEPHI K. MAXPHI*K 2880 LINE (300,147)-(500,347),,B 2890 LINE (300,247)-(500,247) 2900 LINE (400,147)-(400,347) 2910 L O C A T E [8,57:PRINT "R(proj.)" 2920 L O C A T E I 1,45:PRINT "Z" 2930 F O R I-2 T O N 2940 LINE (100"R(I- [)*COS(THETA(I- I) - A V E T H E T A ) / M A X RHO+400,- 1000Z([I)/MAXRHO÷247)-( IOO*R(1)*COS(THETA(1)-A V E T H E T A )/MAXRHO+400,- 100*Z([) MAXRHO+247) 2950 NEXT I 2960 A- 100%IN(AVEPHI):B= 100"COS(AV EPHI) 2970 LINE (400,247)-(A÷400.-B+2,17)

1047

D. ENGI CIRCLE (A+400,-B÷247).3 RETURN REM REM This subroutine prints out the X,Y input data and the REM statistics of THETA. REM LPRINT " n'," X'," Y','THETA" FOR I-1 TO N LPRINT I,X(1),Y(I),INT(K* 100*THETA(1))/100 NEXT ! LPRINT;LPRINT:LPRINT TAB(10) "Average" TAB(20) "Std. Dev." TAB(30) "Minimum" TAB(40) "Mode" TAB(S0) "Maximum" 3090 LPRINT "Theta'; 3100 LPRINT TAB(10) INT(100*K"AVETHETA)/100 TAB(20) INT(100"K*SQR(VARTHETA))/100 TABOO) INT(100"K*MINTHETA)/100 TA B(40) INT( 100"K'MODETHETA)/100 TAB(50)INT( 100*K*MAXTHETA)/100 3110 LPRINT CHRSII2) 3120 RETURN 3130 REM 3140 REM This subroutine prints out the X,Y,Z input data and the 3150 REM statistics of THETA and PHI. 3160 REM 3170 LPRINT " n" TAB(7) "X" TAB(20) "Y" TAB(33) "Z" TAB(44) "THETA" TAB(59) "PHI" 3180 FOR I-I TO N 3190 LPRINT I TAB(5) X(1) TAB(18) Y(I) TAB(31) Z(D TAB(44) INT(K*I00*THETA(1))/100 TAB(ST) INT(K*I00*PHI(I))/I00 3200 NEXT I 3210 LPRINT:LPRINT:LPRINT TAB(10) "Average" TAB(20) "Std. Dev." TAB(30) "Minimum" TAB(40) "Mode" TAB(50) "Maximum" 3220 LPRINT "Theta'; 3230 LPRINT TAB(10) INT(100*K'AVETHETA)/100 TAB(20) INT(100*K"SQR(VARTHETA))/100 TAB(30) INT(100*K*MINTHETA)/100 TAB(40) INT(100*K*MODETHETA)/100 TAB(50)INT(100*K*MAXTHETA)/100 3240 LPRINT "Phi'; 3250 LPRINT TAB(10) INT(100eK*AVEPHI)/100 TAB(20) INT(100*K*SQR(VARPHI))/100 TAB(30) INT(100"K*MINPHI)/100 TAB(40) INT(100*K*MODEPHI)/100 TAB(50)INT(100*K*MAXPHI)/100 3260 LPRINT CHR$(12) 3270 RETURN 3280 REM 3290 REM This subroutine prints out the R and Z input data and the 3300 REM statistics of Pill. 3310 REM 3320 LPRINT " n"," RIIO","Z'," PHI" 3330 FOR I-I TO N 3340 LPRINT I,RHO(I),Z(I),INT(K*IOO*PHI(I))/IO0 3350 NEXT 1 3360 LPRINT;LPRINT:LPRINT T~B(10) "Average" TAB(20) "$td. Dev." TAB(30) "Minimum" TAB(40) "Mode" TAB(50) "Maximum" 3370 LPRINT "Phi"; 3380 LPRINT TABOO) INT(100*K*AVEPHI)/100 TAB(20) INT(IOO*K*SQR(VARPHI))/IO0 TAB(30) INT(100*K*MINPHI)/100 TA B(40) INT( 100" K*MODEPHI)/100 TA B(50) INT( 100*K*MAXPHI)/100 3390 LPRINT CHR$(12) 3400 RETURN 3410 REM 3420 REM This subroutine prints the input data to the screen 3430 REM and allows alteration of that data. 3440 REM 3450 CLS 3460 PRINT SPC(5) "n','X(n)','Y(n)',"Z(n)",'RHO(n)" 3470 FOR I-I TO N 3480 PRINT SPC(5) I,X(1),Y(1) 3490 NEXT I 3500 IF OPT=I GOTO 3600 3510 FOR I,,I TO N 3520 LOCATE I+1,41 3530 PRINT Z(1) 3540 NEXT I 3550 IF OPT-2 GOTO 3600 3560 FOR l-I TO N 3570 LOCATE 1+1,55 3580 PRINT RHO(i) 3590 NEXT I 3600 PRINT "The data you entered are listed above" 3610 PRINT "To change any value follow this example:" 3620 PRINT X(7)I2.g9" 3630 PRINT "After completing all changes type:" 3640 PRINT " CONT" 3650 STOP 3660 CLS:RETURN 2980 2990 3000 3010 3020 3030 3040 3050 3060 3070 3080

"

Methodology for micros~ismic event location

APPENDIX 2

Program Listing of LOCA TE.BAS 1000 REM 1010 REM 1020 REM 1030 REM 1040 REM 1050 REM 1060 REM 1070 REM 1080 REM 1090 REM 1100 REM I 1l0 REM 1120 REM 1130 REM 1140 REM I 150 REM 1160 REM 1170 REM 1180 REM 1190 REM 1200 REM 1210 REM 1220 REM 1230 REM 1240 REM 1250 REM 1260 REM 1270 REM 1280 REM 1290 REM 1300 REM 1310 REM 1320 REM 1330 REM 1340 REM =350 REM 1360 REM 1370 REM 1380 REM 1390 REM 1400 REM 1410 REM 1420 REM 1430 REM 1440 REM 1450 REM 1460 REM 1470 REM 1480 REM 1490 REM 1500 REM 1510 REM 1520 REM 1530 REM 1540 REM 1550 REM 1560 REM 1570 REM 1580 REM lJ90 REM 1600 REM 1610 REM 1620 REM 1630 REM 1640 REM 1650 REM 1660 REM 1670 REM 1680 REM 1690 REM 1700 REM 710 REM 720 REM 730 REM 740 REM 750 REM 760 REM 1770 REM 1780 REM 1790 REM 1800 REM

This program computes the (minimum variance) location o f a microseismic event generated during hydraulic fracturing. The location is calculated as a function of the location parameters meuured by individual borehole seismic systems. Location parameters are statistically described by input and the minimum varianc~ location is calculated then output. This program was written by Dennis Engi at Sandia National Laboratories, Albuquerque, NM. (1985) Definition of program indices.., ITOOL - Borehole seismic system : I refers to tool 1 2 refers to tool 2 IDIM

= Cartesian dimension : I refers to X 2 refers to Y 3 refers to Z

ISCOR - Spherical coordinates : I refers to RHO 2 refers to PHI 3 refers to T H E T A I T R A N - Trigonometrically transformed spherical coordinates : I refers to RHO 2 refers to sin(PHI) 3 refers to sin(THETA) 4 refers to cos(PHI) 5 refers to cos(THETA) IPAR

1 2 3 4 5

Parameters for the triangular density function : refers to minimum value (A) refers to most likely value (B) refers to maximum value (C) refers to (B-A) refers to (C-A) 6 refers to (C-B) 7 refers to 2/((B-A)(C-A)(C-B))

IMOM = Statistical moments of the density functions : I refers to Ist moment 2 refers to 2nd moment Definition of program variables... PAR(IPAR) = Parameters for the triangular density functions. M O M E N T ( I T O O L , I M O M , I T R A N) : Statistical moments ( I M O M ) o f the transformed spherical coordinates ( I T R A N ) per tool (ITOOL). CCOOR D(ITOOL,IDIM,MOM) : Statistical moments of the Cartesian coordinates per tool (ITOOL) per dimension (IDIM). LSTIM(ITOOL,IDIM) : Event locations with respect to the common reference point per tool (ITOOL) per dimension (IDIM). VAR(ITOOL,IDIM) : Statistical variance of the event location per tool (ITOOL) per dimension (IDIM). TVAR(IDIM) : Total statistical variance o f the event location per dimension (IDIM). EVNTLOC(IDIM) : Minimum variance estimator o f the event location per dimension (IDIM). Definition o f program parameters... TLOC(ITOOL,IDIM) : Location o f each tool (ITOOL) per Cartesian dimension (IDIM) with respect to a common reference point.

1049

D. E~Ga

1050

1810 1820 1830 1840 1850 1860 1870 1880 1890 1900 1910 1920 1930 1940 1950 1960 1970 1980 1990 2000 2010 2020 2030 2040 2050 -'060 2070 2080 2090 2100 2110 2120 2130 2140 2150 2160 2170 2180 2190 2200 2210 2220 2230 2240 2250 2260 2270 2280 2290 2300 2310 2320 2330 2340 2350 2360 2370 2380 2390 2400 2410

2420 2430 2440

2450 2460 2470 2480 2490 2500 2510 2520 2530 2540 2550 2560 2570 2580 2590 2600 2610

2620

REM SCOORD(ITOOL,ISCOR,IPAR) : REM Triangular probability density function REM parameters (IFAR) per spherical coordinate REM (ISCOR) per tool (ITOOL). REM REM REM Dimension the program arrays. REM DIM SCOORD(2,3,3),TLOC(2.3I,CCOORD(2,3,2),LSTIM(2,3),EVNTLOC(3) D I M MOMENT(2,2,5),PAR(7),TVAR(3),VAR(2,3),VARLOC(3) CLS REM REM Compute the degrees-lo-radians conversion factor. REM CONV = 3.1416/180 REM REM Input the n u m b e r o f borehole seismic tools v,'hich REM measured location parameters for this microseism. REM INPUT " NUMBER OF TOOLS='.NTOOLS REM REM Read the tool locations. REM FOR ITOOL = I TO NTOOLS F O R 1DIM = I T O 3 READ TLOC(ITOOL.IDIM) NEXT IDIM,ITOOL DATA -]50,-100.150 DATA -100.80.100 REM REM Input the seismic source location parameters. REM FOR I T O O L = I TO N T O O L S PRINT" T O O l , N U M B E R : *:ITOOL FOR ISCOR - I "lO 3 VSCOR$ = "RIIO" IF ISCOR - 2 T t l E N VSCOR$ - "Pill" IF ISCOR - 3 T I I E N VSCOR$ - " T I I E T A " PRINT " C O O R D I N A T E : ";VSCOR$ FOR IPAR - I T O 3 IF IPAR - I T I I E N I N P U T " M I N I M U M - ";SCOORI)(ITOOI..ISCOR,IP/~,R) 11' I I ' A R . , 2 1 l I E N I N P U T " MODE o ";SCOORD(ITOOI.,IS('OR.IPAR) IF IPAR - 3 "I'IIEN I N P U T " M A X I M U M .. ";SCOORD(I I O O I . , I S C O R . I P A R ) NEXT IPAR,ISCOR CI.S NEXT ITOOL RI!M REM Define the f u n c t i o n a l relationships for the first a n d REM second m o m e n t s of the t r i g o n o m e t r i c a l l y t r a n s f o r m e d , REM t r i a n g u l a r l y d i s t r i b u t e d , spherical r a n d o m variables. REM DEE FNESIN(A,B,C, BMA,CMA,CMB,FAC)-EAC'~((CMA*SIN(13))-(CMI3%IN(A)) (BMA*SIN(C))) D E F FNVSIN(A,B,C,BMA,CMA,CMB,FAC)-(FAC/8)*((4/FAC)( C M A * C O S ( 2 " B))+(CM B"COS(2* A))+(BM A*COS(2"C))) D E F FNECOS(A,B,C,BMA,CMA,CMB,FAC)-FAC*((CMA*COS(BI)-(CMB*COS(A))(BMA'COS(C))) D E F FNVCOS(A,B,C,BMA,CMA,CMB,FAC)-I-FNVSIN(A,B,C,BMA,CMA,CMB.EAC) REM FOR ITOOL=I TO NTOOLS F O R I S C O R - I TO 3 REM REM Initialize the p a r a m e t e r vector for spherical REM c o o r d i n a t e ISCOR per tool I T O O L REM PA R( I ) - S C O O R D ( I T O O L , I S C O R , I ) PA R ( 2 ) - S C O O R D ( I T O O L , I S C O R , 2 ) PA R ( 3 ) - S C O O R D ( I T O O L , I S C O R , 3 ) IF P A R ( I ) > - P A R ( 2 ) O R P A R ( 2 ) : , - P A R ( 3 ) T H E N G O T O 3320 PAR(4)=PAR(2)-PAR( I ) PA R(5),.PA R(3)- PAR(I ) PA R ( 6 ) - P A R ( 3 ) - P A R ( 2 ) IF ISCOR>I T H E N G O S U B 3280 PA R(7)..2/(PA R(4}" PA R(5)* PA R(6)) REM REM Compute the first a n d second m o m e n t s o f the REM t r i g o n o m e t r i c a l l y t r a n s f o r m e d spherical coordinates. REM IF ISCOR>I G O T O 2620 M O M E N T ( I T O O L . I , I ) - ( P A R( I )+PA R(2)*PA R(3))/3 MOM-MOMENT(ITOOL.I.I ) VARi.(PAR(3)*PAR(5)-PAR(I)*PAR(4)-PAR(2)*PAR(6))/18 MOMENT(ITOOL,2, I)-VARI~MOM*MOM) G O T O 2680 ITRAN-ISCOR

Methodology for microseismic event location 2630 MOMENT(ITOOL. I.ITRA N)-FNESIN(PAR( I ).PAR(2),PA R(3).PA R(4).PA R(5).PA R(6), PAR(7)) 2640 MOMENT(ITOOL.2.ITRAN)-FNVSIN(PAR( I ).PA R(2).PAR(3).PA R(4).PA R(5),PA R(6). PAR(7)) 2650 ITRAN.,.ITRAN÷2 2660 MOMENT(ITOOL, 1.ITRAN)- FNECOS(PAR( I ),PAR(2),PAR(3).PA R(4).PA R(5).PA R(6), PAR(7)) 2670 MOMENT(ITOOL,2.1TR AN)- FNVCOS(PA R( I ).PA R(2).PA R(3).PA R(4).PA R(5).PA R(6), PAR(7)) 2680 NEXT ISCOR 2690 NEXT ITOOL 2700 REM 2710 REM Compute the first and second moments in Cartesian coordinates. 2720 REM 2730 FOR ITOOL - I TO NTOOLS 2740 FOR IMOM - I TO 2 2750 CCOOR D(ITOOL.1,1MOM),, MOM ENT(ITOOL.I MOM. 1)* MOM ENT(ITOOL.IMOM.2)* MOMENT(ITOOL.IMOM.5) 2760 CCOORD(ITOOL.2.1MOM),,MOMENT(ITOOL,IMOM.I )*MOMENT(ITOOL,IMOM.2)* MOMENT(ITOOL,IMOM.3) 2770 CCOORD(ITOOL,3.1MOM)-MOMENT(ITOOL,IMOM,I )"MOMENT(ITOOL,IMOM.4) 2780 NEXT IMOM 2790 NEXT ITOOL 2800 FOR IDIM,,I TO 3 2810 TVAR(IDIM)-0 2820 FOR ITOOL - I TO 2 2830 V AR(1TOOL,IDIM)-(CCOORD(ITOOL,IDIM,2)-(CCOORD(ITOOL.IDIM, I ))*2) 28.10 IF VAR(ITOOL,IDIM) q: 0 THEN VAR(ITOOL.IDIM) - 0 2850 TVAR(IDIM)-TVAR(IDIM)÷VAR(ITOOL,IDIM) 2860 NEXT 1TOOL 2870 NEXT IDIM 2880 FOR IDIM- I TO 3 2890 FOR ITOOL - I TO NTOOLS 2900 LSTIM(ITOOL.IDIM)-CCOORD(ITOOL.IDIM.I)÷TLOC(ITOOL.IDIM) 2910 NEXT ITOOL 2920 NEXT IDIM 2930 FOR IDIM - I TO 3 2940 EVNTLOC(IDIM)-.VA R(2,1DIM)*LSTIM( I .IDIM)÷VAR(I.IDIM)*LSTIM(2.1DIM) 2950 EVNTLOC(IDIM)-EVNTLOC(IDIM)/TVA R(IDIM) 2960 NEXT IDIM 2970 CLS 2980 P R I N T " " 2990 PRINT " Event locations and errors for each tool" 3000 PRINT" " 3010 PRINT "TOOL DIMENSION LOCATION ERROR" 3020 FOR ITOOL - I TO NTOOLS 3030 FOR IDIM - 1 TO 3 3040 VDIM$-'X" 3050 IF IDIM,,2 THEN VDIM$-'Y" 3060 IF IDIM,,3 THEN VDIM$-'Z" 3070 PRINT ITOOL,VDIM$;" "; 3080 PRINT USING " ~##~t.m';LSTIM(ITOOL,IDIM)0SQR(VAR(ITOOL,IDIM)) 3090 NEXT IDIM,ITOOL 3100 IF NTOOLS-I THEN GOTO 3360 3110 P R I N T ' " 3120 PRINT " " 3130 PRINT" Minimum-error event location" 3140 PRINT " " 3150 PRINT " DIMENSION LOCATION ERROR" 3160 FOR IDIM - I T O 3 3170 VARLOC(IDIM)-(VAR(I .IDIM)*VAR(2.1DIM))/TVAR(IDIM) 3180 VDIM$-'X" 3190 IF IDIM - 2 THEN VDIM$-'V" 3200 IF IDIM - 3 THEN VDIM$..'Z" 3210 P R I N T " ";VDIM$;" "; 3220 PRINT USING " ##~'~m.~';EVNTLOC(IDIM),SQR(VARLOC(IDIM)) 3230 NEXT IDIM 3240 GOTO 3360 3250 REM 3260 REM Transform from degrees to radians 3270 REM 3280 FOR IPAR - I TO 6 3290 P A R ( I P A R ) . CONV*PAR(IPAR) 3300 NEXT IPAR 3310 RETURN 3320 CLS 3330 PRINT " INVALID PARAMETER SET" 3340 PRINT" " 3350 PRINT " STRICT INEQUALITIES ARE REQUIRED : MINIMUM ,c MODE .c MAXIMUM" 3360 END

CAGEO 1 5 : 7 - B

1051

1052

D.E.~G!

APPENDIX 3 A S a m p l e Problem

The BSS tool is assumed to be located at the origin of a coordinate system in which the seismic source is located at x = 750. y = 433, ." = 500. In addition, the amplitude o f the signal induced at the tool is assumed to be A(t)

= 9 sin [800ru] V.

Since the coordinates above correspond to p = d = 100, ~ = 60. 0 - 30, the rectilinear components o f the induced voltage are 3A(t):4

X(t)

=

Y(t)

= x3A(t)4

Z(t)

= A(t) 2

which in tabular form gives

n

t(msec)

X(V)

Y(V)

X(V)

1 2 3 4 5 6 7 8 9 I0

0.20 0.45 0.70 0.95 1.20 1.45 1.70 1.95 2.20 2.45

3.3 6.1 6.6 4,6 0.8 -3.3 -6.1 -6.6 -4.6 -0.8

1.9 3.5 3.8 2.7 0.5 -1,9 -3.5 --3.8 --2.7 -0.5

2.2 4.1 4.4 3.1 0.6 -2.2 -4.1 -4.4 -3.1 -0.6

rounded to a single decimal place. Using these values as inputs to SPHSTAT, BAS, the following output is obtained.

Theta Phi

AVE

STDEV

MIN

MODE

MAX

30.09 59.79

0.44 0.60

29. I I 57.04

29.93 59.98

31.23 62.35

Line 2080 of LOCATE.BAS is edited to read "2080 DATA 0,0,0" to indicate that the common coordinate system is identical to that o f tool I, Assuming that the velocity factor and time-delay data yield a minimum, mode, and maximum of 950, I000, and I050, respectively,for the distance ("RHO") to the seismic event, the entry data for L O C A T E . B A S are as follows.

Parameter

Value

pmin pmode pmax • min •mode emax 0rain 0mode 0max

950 I000 1050 57.04 59.98 62,35 29. l 29.93 31.23

and the output is

TOOL

DIMENSION

LOCATION

STDEV

l 1 I

X Y Z

747.5 433.2 503. I

20.4 10.9 17.4

which may be compared to the original coordinates for the seismic source.