Method for Estimating Permeability of Earth Formations by Processing Stoneley Waves From an Acoustic Wellbore Logging Instrument

1. Field of the Invention

The invention is related to the field of acoustic wellbore logging. More specifically, the invention is related to methods for determining the permeability to fluid flow of earth formations by processing signals recorded by an acoustic wellbore logging instrument.

2. Description of the Related Art

Acoustic wellbore logging instruments are used to measure velocities of earth formations in one or more modes of acoustic energy propagation. Acoustic wellbore logging instruments are typically used in liquid-filled wellbores drilled through the earth formations. Velocity is generally determined using these instrument by measuring the time taken by an acoustic energy pulse to traverse a particular distance along the wall of the wellbore. The wall forms the interface between the liquid in the wellbore and the earth formations.

One form of acoustic energy pulses of particular interest to the invention is referred to as "Stoneley" waves. Stoneley waves, also known as tube waves, are essentially guided pressure pulses travelling in the wellbore. It had been determined in earlier research that a relationship is likely to exist between the transmission properties of Stoneley waves as they propagate along the wellbore, and the hydraulic characteristics of the formations along the wellbore wall. See for example, J. Rosenbaum, Synthetic Microseismograms: Logging in Porous Formations, Geophysics, vol. 39, pp. 14-32, Society of Exploration Geophysicists (1974). Determining formation permeability was not practical using the acoustic logging instruments available at the time the Rosenbaum reference was published because those instruments typically did not generate detectable Stoneley waves, and in the instances where they did so, separation of the Stoneley waves from the acoustic signals as detected was very difficult.

It was later shown to be possible to more reliably detect Stoneley waves generated by an acoustic well logging instrument. See for example, D. Williams et al, The Long Spaced Acoustic Logging Tool, paper "T", transactions of the 25th annual logging symposium, Society of Professional Well Log Analysts (1984). Verification of the theoretical relationship between Stoneley wave propagation and permeability proposed by Rosenbaum, supra, was then performed by laboratory experiment. See for example, K. Winkler et al, Permeability and Borehole Stoneley Waves: Comparison Between Experiment and Theory, Geophysics, vol. 54, pp. 66-75, Society of Exploration Geophysicists (1989). This theoretical model was also used to estimate permeability from acoustic well log data. See for example, D. R. Burns et al, Permeability Estimation from Full Waveform Acoustic Logging Data, The Log Analyst, vol. 29, pp. 112-122 Society of Professional Well Log Analysts, Houston, Tex. (1988).

In practice, the method used by Burns et al included estimation of the velocity and attenuation dispersion of Stoneley waves in permeable formations, and inverting the velocity and attenuation dispersion to obtain permeability using the relationship described by Rosenbaum, supra. The method described by Winkler et al was impractical and difficult to perform because of the number of unknown parameters, such as the wellbore diameter along the wave propagation path, and reflection coefficients at boundaries between different layers of the earth formations.

A simpler parametric model for estimating permeability from Stoneley wave propagation was developed by Tang et al. See, for example, X. M. Tang et al, Dynamic Permeability and Borehole Stoneley Waves: A Simplified Biot-Rosenbaum Model, Journal of the Acoustic Society of America, vol. 90, pp. 1632-1646 (1991). The model described by Tang et al separates propagation of Stoneley waves into an "elastic" (or non-permeability related) component, and a "permeable" component. The elastic component is well known in the art. The permeable component could be calculated using a model developed by Johnson et al. See for example, D. L. Johnson et al, Theory of Dynamic Permeability and Tortuosity in Fluid-Saturated Porous Media, Journal of Fluid Mechanics, vol. 176, pp. 379-400 (1987). Without more, the Tang et al model using the Johnson et al method for determining the permeable component has proven impractical because it is difficult to separate the permeable and elastic components of the Stoneley waves.

Another difficulty in determining formation permeability from Stoneley waves has been separation of reflected from directly transmitted Stoneley waves in the acoustic signals. Interference between transmitted and reflected waves obscures the true amplitude of the directly transmitted Stoneley wave.

The invention is a method for determining the permeability of earth formations penetrated by a wellbore from measurements made by an acoustic wellbore logging instrument. The method includes separating components of the acoustic measurements which represent Stoneley waves propagating directly through the earth formations. Stoneley waves are then synthesized using as parameters the elastic properties of the earth formations, and measurements of the wellbore diameter. The separated, directly propagated Stoneley wave components are then compared to the synthesized Stoneley waves to extract non-permeability related propagation effects from the Stoneley wave components measured by the instrument. The permeability is calculated from the difference between the permeability related propagation effects in the separated components of the Stoneley waves measured by the instrument, and the synthetic Stoneley waves.

In a preferred embodiment of the invention, the permeability is determined by inversion processing the Stoneley center frequency shift and travel time delay with respect to the permeability.

FIG. 1 shows an acoustic well logging instrument as it is used to acquire signals useful with the method of this invention.

FIG. 2 shows a graph of signals received at each one of the receivers on the acoustic logging instrument of FIG. 1.

FIG. 3 shows a graph of signals received by one of the receivers on the instrument of FIG. 1 when the instrument is moved through the wellbore.

FIG. 4 shows a representation of layers of earth formation used for a method of simulating Stoneley wave propagation through the earth formations.

FIG. 5 shows a diagram of propagation vectors used to update the Stoneley wave simulation as the logging instrument is moved through the wellbore.

FIG. 6 is a flow chart of the method for determining permeability according to the invention.

FIG. 1 shows an acoustic well logging instrument as it is used to acquire signals suitable for processing according to the method of this invention. The instrument 10 is inserted into a wellbore 2 drilled through earth formations 20. The instrument 10 can be inserted into and withdrawn from the wellbore 2 by means of an armored electrical cable 14 spooled from a winch (not shown) or any similar conveyance known in the art.

The wellbore 2 is typically filled with a liquid 4, which can be "drilling mud" or any similar fluid usually used for drilling or completion of wellbores. The instrument includes an acoustic transmitter 6 and one or more acoustic receivers, shown generally at 8. The transmitter 6 is periodically energized and emits acoustic energy, shown generally as waves 16, which travel through the liquid 4, along the wall of the wellbore 2, and eventually reach the receivers 8. The receivers 8 generate electrical signals in response to the detected acoustic energy. This is shown in the flow chart in FIG. 6 at 60. The instrument 10 typically includes a telemetry module, shown generally at 12, which converts the electrical signals from the receivers 8 into a form suitable for recording and transmission over the cable 14 to a surface recorder 18. It should be understood that the number of receivers 8 and transmitters 6 as shown in FIG. 1 is only an example and is not meant to limit the invention.

Several different types of transmitters 6 and receivers 8 are known in the art, some of which generate and detect primarily flexural energy, and others which generate and detect primarily compressional energy. The type of transmitter and receiver used is not meant to limit the invention. It is only necessary that the transmitters and receivers be capable of generating and detecting acoustic energy within a frequency range of about 1-2 KHz. Preferably, compressional transmitters and receivers which include 1-2 KHz within their bandwidths should be used, since the compressional mode of acoustic energy transmission is likely to generate detectable Stoneley waves.

It should also be understood that the method of the invention can be performed at any location where there is a suitable computer (not shown). The method of the invention can be performed in the surface recorder 18 if appropriately programmed, but this is not meant to limit the invention.

1. Separation of Directly Propagated Stoneley Waves From the Received Signals

The first step in the method of the invention is to separate various components of the signals detected by the receivers 8 so that Stoneley wave components having a particular travel path through the earth formations can be identified and processed. Typically, Stoneley waves will have substantially all of their energy within a frequency range of about 1-2 KHz. The signals from the receivers 8 should first be filtered to exclude energy at frequencies above this range. Filtering will remove other propagation mode components such as compressional and shear modes. The filtered signals will then consist substantially entirely of Stoneley wave components. This is shown on the flow chart in FIG. 6 at 62.

The next step in the process of separating the Stoneley wave components is to identify the Stoneley wave components which represent wave propagation along the particular travel path. In the example of FIG. 1, there are three likely travel paths for the Stoneley waves. Some of the Stoneley waves travel along the wellbore as shown at 16. Other Stoneley waves can travel downwardly from the transmitter 6, and reflect from a boundary 20B between layers of the earth formations 20 having different acoustic impedances. The reflected waves will generally be detected at later times by the receivers 8, and can mix with later portions of the directly propagated waves. Other Stoneley waves can travel upwardly and reflect from another boundary, such as shown at 20A, and be detected by the receivers 8.

Referring now to FIG. 2, the process by which the various travel path components of the Stoneley waves are separated from the filtered signals will be explained. The graph in FIG. 2 represents the filtered signals, shown as curves 22, 24, 26, 28 and 30, from each one of the receivers (8 in FIG. 1) on the instrument (10 in FIG. 1) which were measured at a single axial instrument position (depth) in the wellbore (2 in FIG. 1). As can be inferred from the individual curves 22-30, events in each curve corresponding to propagation of Stoneley waves along one or more transmission paths can be identified. These events occur in each signal at a time which varies with respect to the distance between the individual receiver and the transmitter. the amount of variation with respect to transmitter-to-receiver distance corresponds to the inverse of the propagation velocity ("slowness") of the Stoneley waves along the wellbore.

For example, in each of curves 22-30, a first event shown at 22A-30A, can be observed. The first events 22A-30A typically correspond to the direct transmission of Stoneley waves along the wellbore wall from the transmitter to the receivers, and therefore have the shortest overall propagation time to any receiver. The first events 22A-30A will have a time variation with respect to the transmitter-to-receiver distance, called the "moveout" and shown at 32, which corresponds to the Stoneley wave slowness of the earth formations through which the Stoneley waves directly propagate.

Second events 22B-30B can be observed in each of the curves 22-30. The second events 22B-30B correspond to Stoneley waves which reflect from an acoustic impedance boundary located axially below the transmitter (for the instrument configuration as shown in FIG. 1), and are subsequently detected by the receivers. The second events 22B-30B have a moveout shown at 34 which represents substantially the same slowness as the moveout 32 for the first events 22A-30A because the reflected Stoneley waves travel through the same portion of the formation to each receiver as do the directly transmitted Stoneley waves. The first events 22A-30A and the second events 22B-30B are shown as being substantially time separated for clarity of the illustration in FIG. 2, but it should be noted that if a boundary is located axially proximal to the transmitter, there could be substantial interference between the first events 22A-30A and the second events 22B-30B. A process for separating the first events from the second events will be further explained.

Third events, shown generally at 22C-30C, can be caused by Stoneley waves being reflected from the boundary located axially above the instrument (as shown in the configuration of FIG. 1). The moveout 36 of the third events 22C-30C is typically opposite in direction to moveout of the first and second events. It can be readily seen how interference between first, second and third events can make it difficult to determine which part of the Stoneley waves result from direct transmission of the energy through the earth formations from the transmitter to the receivers.

In the method of the invention, different components of the Stoneley waves resulting from wave propagation along different travel paths can be separated by selecting the signals from a single one of the receivers. This same wave mode separation can also be used with acoustic well logging instruments which only include one or two receivers, where an array graph such as shown in FIG. 2 would be impossible to construct. Referring now to FIG. 3, the signals from a single one of the receivers (8 in FIG. 1) are shown at 40-50 as being measured at a plurality of different axial (depth) positions of the instrument (10 in FIG. 1) within the wellbore (2 in FIG. 1). As a matter of convenience for the system designer, the axial (depth) separation between each one of the axial positions at which the signals were measured is about equal to the nominal sample interval between successive well log measurements. The nominal sample interval is typically 1/4 to 1/2 foot. The number of individual signals which can be used should correspond to a total depth span about equal to the distance between the receiver nearest to, and the receiver farthest away from, the source (6 in FIG. 1) (this distance being referred to as the "array aperture"). For a typical acoustic array well logging instrument this can include about 8-12 individual signals representing a receiver array aperture of about 31/2 feet.

As long as the Stoneley wave velocity does not substantially change in the formation located between the selected receiver and the source, the events such as shown at 40A-50A, which correspond to the directly propagated waves, will have substantially zero moveout as shown at 52, because the wave propagation path length does not substantially change as the instrument 10 is moved along the wellbore 2. Events corresponding to reflected waves, such as shown at 40B-50B for reflections from the boundary (such as 20A in FIG. 1) located above the instrument 10 will typically have a moveout 56 which corresponds to twice the slowness of the formation through which the Stoneley waves propagate. The moveout for reflected Stoneley waves is doubled relative to the moveout shown in the multiple-receiver signal graph in FIG. 2 because as the instrument 10 moves towards or away from a reflective boundary (such as 20A in FIG. 1), the total acoustic wave travel path between the source 6 and the particular receiver shortens at twice the rate of movement of the instrument 10. Events shown at 40C-50C correspond to Stoneley waves reflected from the boundary (20B in FIG. 1) below the instrument 10. The length of this wave travel path will increase at twice the rate of movement of the instrument, so the moveout of the reflected waves from the lower boundary 20B will typically represent twice the Stoneley wave slowness represented by the moveout shown in FIG. 2, and will have opposite sign with respect to the events 40B-50B representing the waves reflected from the upper boundary 20A.

The manner in which the Stoneley wave propagation modes can be separated from each other is explained as follows. A quantity p can represent the total number of individual wave propagation modes (such as the previously explained direct, upwardly reflected or downwardly reflected mode) present in the signals measured by an array (such as shown at 8 on the instrument 10 in FIG. 1) consisting of N equally spaced apart acoustic receivers. The receivers preferably have a uniform spacing between each of them represented by d. Each individual mode of propagation will have its own unique frequency spectrum, represented by h(.omega.). To estimate the Stoneley wave frequency spectrum for each one of the p propagation modes from signals which would be detected at the n-th receiver position, it is possible to "forward" or "reverse" propagate these wave modes to each individual receiver position in the array from any other individual receiver position in the array, by using the respective values of slowness for each one of the p Stoneley wave propagation modes. The slownesses can be determined by using the moveout, as shown in the graphs in FIGS. 2 and 3. Other method of determining slowness are known in the art. See for example U.S. Pat. No. 5,541,890 issued to Tang.

Propagation of one of the wave modes, h.sub.r (.omega.) from receiver position n to receiver position m can be calculated by the expression: h.sub.r (.omega.)Z.sub.r.sup.m-n, where Z=exp(i.omega.s,d). s.sub.r represents the slowness of the r-th wave propagation mode. m represents a receiver position index for the "propagated to" receiver position, and can be smaller than, equal to, or greater than the index number n of the receiver position of interest (the "propagated from" receiver position). It is then possible to equate the sum of the propagated frequency spectra to the frequency spectra of the measured signals at each receiver position. This results in the expression: ##EQU1## A typical array acoustic well logging instrument which generates signals useful in the method of this invention is described, for example, in a sales brochure entitled, Multipole Array Acoustilog, Western Atlas Logging Services, Houston, Tex. (1993). The array acoustic well logging instrument described in this reference includes a compressional transmitter and eight compressional receivers. Because the number of receivers in the receiver array, N, of this typical array acoustic logging instrument is greater than the number of propagation modes p, equation (1) can be solved by using a least-squares method for each frequency to obtain the estimated spectra (h) for each propagation mode at the receiver position n. Transforming the resulting frequency spectra into the time domain provides the signal waveform which would obtain for each one of the p wave propagation modes.

Referring once again to FIG. 2, the significance of the moveout 32 associated with the directly propagated waves is that this moveout will represent substantially the same velocity as the moveout 34 of the events 22B-30B present in the waves reflected from the boundary below the instrument. The reflected waves can be separated from the directly propagated waves by using the single receiver moveout technique as shown in FIG. 3, where the moveout value for the directly propagated waves will be substantially equal to zero. It should also be noted that if the acoustic logging instrument configuration is reversed from that shown in FIG. 1, where the transmitter is located above the receivers, then the moveout resulting from waves reflected from boundaries above and below the instrument will be reversed in direction from those moveouts as shown in FIGS. 2 and 3. The principle of wave separation is identical in either case. Separation of directly propagated wave components from reflected wave components is shown in the flow chart of FIG. 6 at 64.

2. Stoneley Wave Simulation by the Propagator Matrix Method

The next part of the method of the invention is to simulate Stoneley waves propagating within a wellbore having an irregular wall surface, and through formations having heterogenous acoustic propagation characteristics. The simulation method which will be explained herein allows the wellbore wall surface and the formation acoustic propagation characteristics to be selected and varied arbitrarily. Actual well log instrument measurements of the acoustic propagation characteristics and of the diameter of the wellbore can be used as inputs to this simulation method. The output of the simulation method is a data set representing Stoneley wave amplitude and phase variation due to wellbore diameter variations and due to changes in the formation acoustic propagation characteristics along the Stoneley wave propagation path.

Referring now to FIG. 4, the wellbore 2 wall and the surrounding earth formations 20 can be divided, or "discretized", into a "stack" of thin, individual layers, such as the ones shown at 20C and 20D. The wall of the wellbore 2 need not be smooth, and in fact can be of irregular shape, as represented by the different radius of each individual layer 20C, 20D. In FIG. 4, the radius of the wellbore 2 at the axial position of layers 20C and 20D is represented by R.sub.j and R.sub.i, respectively. Radii R.sub.j and R.sub.i are shown extending from the central axis 14A of the wellbore 2, outwardly to the portion of the wellbore wall corresponding to each layer 20C, 20D. The acoustic well logging instrument 10 is shown as being substantially centered in the wellbore 2.

Each layer used in the model shown in FIG. 4 represents a depth interval (an axial span along the wellbore 2), which as a matter of convenience can be set equal to the well logging instrument 10 data sample interval (typically 1/4 or 1/2 foot as previously explained). The wellbore 2 diameter (or radius) at each layer 20C, 20D, and the layer's elastic properties can be taken directly from measurements of acoustic log properties and a caliper log. The elastic properties include the compressional and the shear velocities, which can be measured using an acoustic well logging instrument such as one described in the Multipole Array Acoustilog reference, supra. The well log data needed to calculate the elastic properties should also include the bulk density of each formation layer. The effect on Stoneley wave propagation of the presence of the logging instrument 10 inside the wellbore 2 can be determined by using a method described in X. Tang et al, The Effects of a Logging Tool on the Stoneley Wave Propagation on Elastic and Porous Media, The Log Analyst, vol. 34 pp 46-56, Society of Professional Well Log Analysts, Houston, Tex. (1993).

Stoneley wave propagation is simulated using a propagator matrix method, wherein matrix elements describe the coupling of Stoneley waves across the layers, each layer having a different wellbore radius and different acoustic propagation characteristics. In the propagator matrix method, (b.sup.+ b.sup.-).sup.T can represent a vector containing both up-going (+) and down-going (-) Stoneley wave amplitude coefficients. The propagation of the amplitude coefficient vector from an axial (wellbore depth) position represented by z.sub.1, through some of the formation layers to another axial (wellbore depth) position, represented by z.sub.2, can be described by the equation: ##EQU2## where L represents the number of formation layers located between axial (depth) positions z.sub.1, and z.sub.2. The propagator matrix for the l-th layer can be derived using a method described in S. Gelinsky et al, Modeling of Stoneley Waves in a Heterogeneous Formation With an Irregular Borehole, 59th Annual Conference EAEG, Geneva Switzerland (1997). The propagator matrix can be described by the expression: ##EQU3## where z.sub.1 represents the axial (depth) position of the l-th earth formation layer and k.sub.l-1 and k.sub.l represent the Stoneley wavenumbers of the (l-1)-th and l-th formation layers, respectively. Q.sub.l and Q'.sub.l can be computed by an integral expression I(i, j), wherein:

Integral expression I(i, j) is based on an orthogonality relation. This relation is described for example, L. Tsang, Transient Acoustic Waves in a Fluid-Filled Borehole With Horizontal Bed Boundary Separating Two Solid Formations, Journal of the Acoustical Society of America, vol. 81, pp. 844-853 (1987). The orthogonality relation described in the Tsang reference connects wave-induced stresses and displacements across a boundary between adjacent earth formation layers i (shown as 20C in FIG. 5) and j (shown as 20D in FIG. 5). The method for calculating the integral expression I(i, j) will be further explained.

The radius of the well logging instrument 10 and the wellbore radius in the i-th and j-th layers can be represented by R.sub.t, R.sub.i and R.sub.j, respectively. .sigma..sub.z, .sigma..sub.rz, u.sub.z, and u.sub.r, respectively, represent the Stoneley wave-induced axial normal stress, the axial-to-radial shear stress, the axial displacement, and the radial displacement. The functional forms of these quantities are described, for example, in C. Cheng et al, Elastic Wave Propagation in a Fluid-Filled Borehole and Synthetic Acoustic Logs, Geophysics, vol. 87, pp. 5477-5484, Society of Exploration Geophysicists (1981). Assuming for purposes of this explanation that the Stoneley wave propagates from layer i (20C in FIG. 4) to layer j (20D in FIG. 4), the expression for calculating the integral expressions I(i, j) and I(j, i) should take into account two situations, the first being that R.sub.i >R.sub.j, and the second being that R.sub.i <R.sub.j.

For first the case, where R.sub.i >R.sub.j, the following expression applies: ##EQU4## For the second case, where R.sub.i <R.sub.j, the following expression applies: ##EQU5## where subscript f denotes that the quantity is defined as being located within the liquid (shown as 4 in FIG. 4) in the wellbore (2 in FIG. 4). The corresponding I(j, i) integrals can be obtained by exchanging the indices i and j within I(i, j). For example, integral expression I.sup.a (j, i) can be obtained in this way from I.sup.b (j, i), and vice versa. When the radius of the wellbore 2 is different from one layer to the next layer, a mass conservation method such as described in K. Tezuka et al, Modeling of Low-Frequency Stoneley Wave Propagation in an Irregular Borehole, Expanded Abstracts, 64th Annual Meeting, Society of Exploration Geophysicists (1994) can be used to account for the difference in radius on the propagation of the Stoneley waves. This method normalizes the axial fluid displacement for the individual layer having the larger wellbore radius. More specifically, for the case where R.sub.i >R.sub.j, the first integral in I.sup.a (i, j) is multiplied by the expression ›(R.sub.i.sup.2 -R.sub.t.sup.2)/(R.sub.j.sup.2 -R.sub.t.sup.2)!. Correspondingly, for the case where R.sub.j >R.sub.i, the first integral in I.sup.a (i, j) is multiplied by the expression ›(R.sub.j.sup.2 -R.sub.t.sup.2)/(R.sub.i.sup.2 -R.sub.t.sup.2)!. These expressions, or "normalization factors", represent the ratio of areas in the liquid 4 in the wellbore 2 between layer i 20C and layer j 20D.

The I(i, j) term in equation (4) can then be calculated using equation (5) or (6). For the case i=j=l-1 (representing the same individual layer), either equation (5) or (6) can be used to calculate the expected received signal W.sub.l-1. For the case where R.sub.i >R.sub.j, equation (6) can be used to calculate Q.sub.1 =I.sup.b (l, l-1), and equation (5) can be used to calculate Q.sub.1 '=I.sup.a (l, l-1). For the case where R.sub.i >R.sub.j, equation (6) can be used to calculate Q'.sub.1 =I.sup.b (l-1, l), and equation (5) can be used to calculate Q.sub.1 =I.sup.a (l, l-1). All these integral expressions have closed form analytical solutions, which can be calculated by using an appropriate computer program in order to provide fast simulation of Stoneley wave propagation across the layers of the earth formations.

A procedure to implement the simulation method just described to Stoneley wave propagation through actual earth formations will now be explained. Referring once again to FIG. 4, consider the acoustic instrument 10 as being positioned within the stack of formation layers. In the simulation method, the length (equivalent to the axial or wellbore depth extent) of the stack is set so that Stoneley wave reflections from the upper and lower ends of the stack will occur at times after end of the time window of acoustic signal measurement. For the typical logging instrument this time window is about 4 to 5 milliseconds (4,000 to 5,000 microseconds) from the time of actuation of the source 6. Since the Stoneley wave velocity through the layers for which the simulation is being run is known or can be determined from the measured acoustic signals, a stack length can be selected which will have a simulated reflected wave travel time calculated to cause the reflected Stoneley waves to arrive after the end of the recording time window. The stack length, therefore, is typically only a few feet long.

For any selected set of simulation calculations according to this method, the source 6 would be positioned in the wellbore 2 adjacent to one particular layer, such as shown at 20G and referred to hereinafter as the "source layer". hsr and h can denote, respectively, the product matrix from the upper boundary of the source layer 20G to a layer 20H in which a centrally positioned receiver 8A in the array of receivers 8 is located (this layer being referred to as the "receiver layer"), and the product matrix to the upper end 20F of the stack. hb can denote the product matrix from the lower boundary of the source layer 20G to the lower end 20E of the stack of layers. The synthetic propagation of a Stoneley wave emitted from the source 6 positioned in the wellbore 2 adjacent to the source layer 20G is completely described by the product matrices hsr, h, and hb. This can be better understood from the following analysis.

With respect to the source layer 20G, the Stoneley wave amplitude coefficient vectors at the upper boundary and at the lower boundary of each layer (in this case the source layer 20G) are related to excitation caused by the source 6 according to the expression: ##EQU6## where k.sub.s represents the source layer 20G Stoneley wavenumber and d.sub.s represents the distance between the source 6 (the transmitter) and the lower boundary of the source layer 20G. S(.omega.) represents the Fourier spectrum of the energy output of the source 6. The amplitude coefficient vector at the upper boundary of the source layer 20G, from equation (7), can be propagated to the upper end 20F of the stack of layers by using the expression for h, and the source layer 20G lower boundary amplitude coefficient vector can be propagated to the lower end 20E of the stack using the expression for hb. Using boundary conditions which may be specified for the two ends 20E, 20F of the stack, plus equation (7) is sufficient to determine the two Stoneley wave amplitude coefficient vectors in equation (7). Then by propagating the amplitude coefficient vectors from the source 6 position to the central receiver 8A position using hsr provides the propagated Stoneley wave signal which would obtain at the central receiver 8A position.

To simulate the logging process of an acoustic well logging instrument in an actual wellbore, the procedure just described can be repeated for a plurality of consecutive source 6 positions. This can be described as a solution to the following problem: given hsr, h, and hb, how can these matrices be updated to obtain corresponding matrices hsr',h', and hb' when the source layer is moved into the next layer within the earth formation. The updating can be performed by calculating the propagator matrix (from equation (3)) adjacent to the top end of the current model stack (see FIG. 4), and denote the resulting matrix by p. The new matrices can be updated according to the expressions:

where p.sub.s-, p.sub.s+, p.sub.r and p.sub.b represent the propagator matrix at the source layer 20G, the layer on top of the source layer, the receiver layer 20H, and the lowermost layer of the stack, respectively. These matrices have already been calculated for previous depth positions of the logging instrument. The updating scheme is illustrated in FIG. 5. Once hsr', h', and hb' have been determined, the Stoneley wave propagation for the new source position is completely determined.

Since the simulated Stoneley waves will include both directly propagated and reflected wave components, the simulated Stoneley waves can also be gathered into common source gathers, such as shown in FIG. 2 and into common offset gathers, such as shown in FIG. 3, for the purpose of applying the wave separation method of equation (1). The simulated Stoneley waves can then be split into directly propagated, up-going and down-going simulated Stoneley waves, just as was done for the measured Stoneley waves. The directly propagated synthetic Stoneley waves are preferably used in the permeability analysis which will be explained below. Synthesis of Stoneley waves using the propagator matrix method is shown in the flow chart in FIG. 6 at 66.

3. Permeability Indication

Indication of permeable earth formations can be useful even without quantitative determination of the permeability. Permeability indication can be performed by comparing the measured Stoneley wave data derived from the wave separation processing described previously herein with the numerical simulation of the Stoneley waves just described. By simulating the elastic wave propagation effects of the actual earth formations penetrated by the wellbore, any effects on the Stoneley waves which are not related to the permeability can be accounted for. Then, the differences between the measured and simulated Stoneley waves will mainly reflect the permeability effects.

First, a "reference" formation layer, in which the permeability is known (which can also include zones known to be impermeable) is selected. Then a synthetic waveform "log", representing a synthetic Stoneley waveform for each axial position in the wellbore, can be calculated using the expression: ##EQU7## where M represents the synthetic Stoneley wave spectral response (or transfer function) directly along the wave propagation path. The response of the measured Stoneley wave signals, as previously explained, can be obtained through wave separation. W represents the spectrum of the measured Stoneley wave signals. The spectral ratio in equation (9) represents a source (excitation) spectrum, obtained by deconvolving the propagation transfer function from the measured wave at the reference position z.sub.r. Then, convolving the source spectrum with the transfer function at any depth (axial) position z will generate the synthetic Stoneley wave at that depth.

For a given wave spectrum (synthetic or measured), the wave center frequency and the variance of the spectrum can be calculated using the expressions:

A center time can be computed from the waveform of W using the expression:

The frequency shift and travel-time delay between the measured wave and synthetic wave can be defined by the expressions:

Alternatively, the time delay, .DELTA.T.sub.c, can be calculated using a cross-correlation method, where the delay is determined as the value for time shift of the measured wave required to maximize the degree of correspondence (correlation) between the measured wave and the synthetic wave.

A theory known in the art for Stoneley wave propagation in a permeable borehole indicates that permeability has two direct effects on the Stoneley wave, these being an increase in wave attenuation and decrease in the wave speed. The first effect is characterized by the frequency shift, while the second is characterized by the travel-time delay. Therefore, these two quantities can be used to indicate formation permeability. For example, if the frequency shift corresponds to the travel-time delay within a particular depth interval, it is very likely that these effects are indicative of formation permeability. A formation exhibiting these Stoneley wave properties can thus be characterized as a permeable formation. Indication of a permeable formation by comparing synthetic Stoneley waves to measured Stoneley waves is shown in the flow chart of FIG. 6 at 68.

4. Permeability Inversion Processing

Quantitative estimation of permeability can be performed using a fast model theory such as described in X. Tang et al, Dynamic Permeability and Borehole Stoneley Waves: a Simplified Biot-Rosenbaum Model Journal of the Acoustical Society of America, vol. 90, pp. 1632-1646 (1991) and an inversion process such as one described in R. Chunduru et al, 2-D Resistivity Inversion Using Spline Paramerterization and Simulated Annealing, Geophysics, vol. 61, pp. 151-161, Society of Exploration Geophysicists (1996). The Stoneley wave attenuation due to the fluid in the wellbore and due to the intrinsic attenuation of the earth formations is also estimated, because the intrinsic attenuation also affects the Stoneley wave attenuation and the Stoneley wave travel time. As compared to the effects caused by permeability, the intrinsic attenuation has a different frequency dependency and a different effect on wave attenuation and travel time. This enables these factors to be simultaneously estimated along with the permeability. Including intrinsic attenuation in the inversion procedure allows accounting for the non-permeability related effects on Stoneley wave propagation through the wellbore and earth formations.

Within a permeable earth formation, the Stoneley wavenumber k can be determined by the expression: ##EQU8## where K.sub.n (n=0, 1) represents a modified Bessel function of the second kind of order n, R represents the wellbore radius in the permeable formation and R.sub.t represents the well logging instrument radius. .mu., .rho..sub.0 and .nu..sub.0 represent, respectively, the viscosity of the fluid in pore spaces of the permeable formation, the density of this fluid, and its Stoneley wave velocity. .phi. represents the fractional volume of pore space of the permeable formation (called "porosity"). .gamma. represents a correction factor which is related to the permeable formation's shear wave velocity. .zeta. represents a correction factor which is related to the compressibility of the solid rock portion of the permeable formation (referred to in the art as the "matrix"). The correction factor can be calculated if the compressional velocity and the shear wave velocity, the porosity, and the density of the earth formation are known (see the Tang et al 1991 reference, supra). .kappa.(.omega.) represents the dynamic permeability (see for example, D. Johnson et al, Theory of Dynamic Permeability and Tortuosity in Fluid-Saturated Porous Media, Journal of Fluid Mechanics, vol. 176, pp. 379-400 (1987)), which can be calculated by the expression: ##EQU9## where .kappa..sub.0 represents the quantitative formation permeability. Wavenumber k.sub.0 in equation (14) corresponds to an equivalent elastic formation whose properties are similar to those of the fluid saturated earth formation. The X. Tang et al 1993 reference describes a method for calculating the wavenumber k.sub.0 in the permeable formation when the well logging instrument 10 is present in the wellbore 2.

The effect of the intrinsic attenuation can be accounted for by a transformation described in Cheng et al, Determination of In-Situ Attenuation from Full Waveform Acoustic Logs, Journal of Geophysical Research, vol. 87, pp. 5477-5484 (1982). This transformation can be expressed as: ##EQU10## where f.sub.0 represents an arbitrary reference frequency, which can be the Stoneley wave center frequency. Q represents the quality factor; 1/Q represents the attenuation, which is mostly caused by the liquid 4 in the wellbore 2, but is partially caused by the permeable formation itself.

With this theoretical model, the theoretical center frequency and theoretical frequency variance of the Stoneley waves can be calculated using the expressions:

where P represents the excitation reduction due to the formation permeability, which can be calculated as P=k.sub.0 /k (See X. Tang et al, The Effects of a Logging Tool on the Stoneley Waves Propagation in Elastic and Porous Formations, The Log Analyst, vol. 34, pp. 46-56, Society of Professional Well Log Analysts, Houston, Tex., 1993). d in this expression represents the Stoneley wave travel distance and e.sup.ikd represents the amplitude decay caused by permeability and by attenuation along the wave propagation path. In equation (16), the theoretical frequency shift is defined as:

The theoretical travel-time delay can be computed using the expression: ##EQU11## Equations (17) and (18) provide the overall frequency shift and travel-time delay of Stoneley waves in a permeable and attenuating wellbore relative to an impermeable and nonattenuating wellbore. Although the travel-time delay is calculated in the frequency domain, it can be shown that this delay is equivalent to the time delay computed directly in the time domain using the method of equations (11) and (12).

By comparing the measured frequency shift and travel-time delay given by equation (12), with the theoretical predictions given by equations (17) and (18), an inverse problem can be formulated, with the intended result of minimizing the following objective function:

where .sigma..sup.2.sub.syn and .sigma..sup.2.sub.th can be computed from equations (10) and (16), respectively. The

first and second terms are so scaled that the frequency shift and travel-time delay have almost equal importance in the objective function. The third term in the objective function is a penalty function related to the differences in variance between the synthetic Stoneley waves' spectral distribution (calculated using equation (10)) and the measured Stoneley waves' spectral distribution (calculated using equation(18)). The coefficient .alpha. is typically very small, so that this term in equation (19) has almost zero contribution in the inversion if the values of the parameters vary only within a range of values likely to be encountered in real earth formations. If the values of .kappa..sub.0 and 1/Q are very large, the third term of equation (19) becomes large also, which drives the inversion parameters back into the expected range of values.

The inversion process can use any appropriate minimization algorithm, including well known local algorithms such as conjugate gradient and Levenberger-Marquardt. The inversion can also use global minimization algorithms such as simulated annealing. See for example L. Ingber, Very Fast Simulated Reannealing, Mathematical Computing, vol. 12, pp. 967-973 (1989) and R. Chunduru et al., 2-D Resistivity Inversion Using Spline Parameterization and Simulated Annealing, Geophysics, vol. 61, pp. 151-161, Society of Exploration Geophysicists (1996). The local minimization algorithms require assigning initial values for the model parameters. The initial values can be selected using simulation results prom the previously processed axial position (wellbore depth). The global algorithms do not need initial values, but may require somewhat longer computational time to complete the inversion. Calculating permeability by inversion processing is shown in the flow chart in FIG. 6 at 70.

5. Improving Axial Resolution

In the previously described inversion process, the permeability effect which is calculated typically represents the average value within an axial span comprising a distance from the source (6 in FIG. 4) to the centrally located receiver (8A in FIG. 4) in the receiver array, typically about 10 to 13 ft. The features in the processed Stoneley wave data, however, usually have an axial resolution of about the aperture of the receiver array (8 in FIG. 1), typically about 3 to 4 ft. With the process to be described below, the permeability can be estimated to an axial resolution of about equal to the receiver array aperture.

The wave travel distance d (representing the distance between the source 6 to the center receiver 8A in the array 8) can be divided into two parts:

where d' represents a distance equal to the array aperture, this being measured from the center of the receiver array 8 to the source 6. d.sub.0 represents the remaining distance linking to the source 6. Assuming that the theoretical travel time delay and frequency shift over d.sub.0 are known from previous estimates, then the theoretical travel time delay and frequency shift over d.sub.0 can be subtracted from the measured time delay and frequency shift using the expression:

Then d' can be used in substitution of d in Equations (17) and (19) for calculating the theoretical frequency shift and time delay. Finally, the calculated quantities and the modified data can be used in the inversion formulation ›from equation (20)! to estimate the permeability over the reduced distance d'.

After the permeability over the reduced distance d' is estimated, the theoretical frequency shift and time delay over the reduced distance d' can be calculated. These quantities are combined with those of the previous instrument axial positions (depths) to provide a new average of the quantities over the distance d.sub.0, for estimating permeability in the next depth interval. This scheme is repeated over the entire processing interval.

6. Fluid Parameter Calibration

By performing an asymptotic analysis of the model theory given in equation (14), it has been determined that the permeability-induced time delay and frequency shift are affected by the following parameter combination (in addition to the formation porosity, and wellbore and instrument radii): ##EQU12## where B.sub.0 represents the pore fluid bulk modulus or "incompressibility". This asymptotic analysis is valid for most permeability values of typically earth formations and for the Stoneley wave frequency range used in this invention, of about 1 kHz. This parameter combination indicates that formation pore fluid parameters such as fluid viscosity, density, and acoustic travel speed, must be determined to calculate the permeability quantitatively. The pore fluid parameters, however, may be difficult to determine if multiple fluids are present in the pore spaces of the formations and the fractional composition of these multiple fluids is unknown. In the invention, a calibration method is used to determine a useful value for these parameters.

The calibration method consists of selecting a few (a minimum of two) depth intervals whose permeability values have been estimated from other measurements (including, but not limited to such measurements as gas flow in core samples, formation fluid pressure tests, nuclear magnetic resonance (NMR) well logs). By choosing one of these depths as reference depths for equation (9), synthetic Stoneley wave data can be computed for the remaining depths. This process (i.e., reference selection and modeling) is rotated for all the depths. Then, by comparing the synthetic data with the measured data, we can estimate .mu.B.sub.0.sup.1/2 by minimizing the following misfit function: ##EQU13## where n represents the total number of depths chosen and the pair i, j indicates using depth position i as the reference to compare with depth j. In the minimization, the known permeability values at the calibration depths are used as a known parameter in the parameter combination described above. Only .mu.B.sub.0.sup.1/2 is varied as the fitting parameter. The value of this parameter that minimizes the value of the misfit function in equation (23) is determined to be the desired fluid parameter value for the depth zone of interest. This fluid parameter value can be used to estimate a continuous permeability over the entire zone of interest.

Those skilled in the art will devise other embodiments of this invention which do not depart from the spirit of the invention as disclosed herein. Accordingly, the invention should be limited in scope only by the attached claims.