Method for Obtaining Dielectric Constant and Conductivity Information on Subsoil Structure Using Continuous Complex Geotomography

1. Field of the Invention

The present invention relates to a method of and a system for deriving information regarding subsoil and rock structure. More particularly, the invention relates to a method and to a system for obtaining at least a two-dimensional plot of subsoil structure and the physical properties thereof, using complex geotomography. The method is used to obtain, for example, conductivity and dielectric constant values of the subsoil in the plane between two or more bore holes.

2. Description of the Related Art

Knowledge of the physical properties of the subsoil structure (stratigraphy) in an area is very useful to predict the condition and components of the soil. The information generated when interpreted is useful to determine the location of underground petroleum, oil shale, coal, mineral deposits, location of contaminants (such as hazardous chemicals, gasoline), location of fractures and faults, and the like. Having a multi-dimensional plot or graph of the subsoil and bedrock structure with sufficiently high resolution enables display of the relative locations of both small and large structures and layers.

Tomography was originally defined to denote an analytic method of X-ray photography in which cross-sectional planes are imaged without interference from adjacent planes. Since the first practical applications to medical diagnostics, tomography has evolved to include reconstruction using frequencies other than X-ray, such as those from acoustic or electromagnetic radiation.

Tomography was first used to reconstruct cross-sections in the earth by D. L. Lager and R. J. Lytle, "Computer Algorithms Useful for Determining the Subsurface Electrical Profile Via High Frequency Probing," Lawrence Livermore National Laboratory (LLNL), Livermore, Calif., Report No. UCRL-51748 (1975). Additional reports of the use of geotomography have applied medical X-ray reconstruction algorithms to earth and subsoil imaging. These include, for example: K. Dines and R. J. Lytle, "Computerized Geophysical Tomography," Proc. IEEE, Vol. 67, No. 7, pp. 1065-1073 (1979); D. L. Lager and R. J. Lytle, "Determining a Subsurface Electromagnetic Profile from High Frequency Measurements by Applying Reconstruction Techniques Algorithms," Radio Science, Vol. 12, No. 2, pp. 249-260 (1977); R. J. Lytle, et al., "Electromagnetic Cross-Borehole Survey of a Site Proposed for a Future Urban Transit Station," LLNL Report No. UCRL-52484 (1978).

In U.S. Pat. No. 4,214,226, K. Narasimhan et al. disclose a method for obtaining high resolution data regarding the subsoil structure between two spaced apart boreholes. Their method includes the steps of placing a number of signal sources at preselected distances from one another in one borehole, placing a number of signal detectors at preselected distances in at least a second borehole, activating the sources in a selected sequence, receiving at each of the detectors the signal which propagates to the detectors through the subsoil structure as a result of each activated source, and processing the signals received at each detector as a result of the activation of the sources. This method usually includes one-time signals at frequencies greater than 160 Hz up to about 1000 Hz to obtain the data adapted to produce subsoil structure information of incremental areas on the order of .lambda./3 for producing a tomograph of the substructure between the boreholes, where .lambda. is the wavelength of the processed signals. The reference tracks only explosive signals and seismic waves produced by the detonation of an explosive charge at different locations. Only a one-time signal is obtained from each source which is measured for amplitude and/or velocity. Once the destructive explosive charge has detonated, that region and the regions below the explosive charge can no longer be monitored on a repetitive basis. The patent does not teach the use of a radiation transmitter which produces a continuous signal which can be detected and measured for both amplitude and phase. On the other hand the borehole used in the present invention is not destroyed and can be used to monitor the subsoil again and again over a period of time.

Additional art of general interest include, for example:

"Controlled Source Geotomography Using a Circular Array," by G. L. Pavlis in Transactions, American Geophysical Union, Vol. 66, No. 18, pp. 303ff (1985);

"Application of Computed Tomography to Geological Investigation," by H. Sugawara et al. in Soil Mechanics and Foundation Engineering (Japan), Vol. 33, No. 7, Ser. No. 330, pp. 53-56 (1985);

"The Use of Computer-Assisted Tomography to Determine Spatial Distribution of Soil Water Content," by J. M. Hainsworth, Australian Journal of Soil Research, Vol. 21, No. 4, pp. 435-443 (1983);

"Geophysical Diffraction Tomography," by A. J. Devaney in a special issue on Electromagnetic Methods in Applied Geophysics, edited by A. Q. Howard, Jr., IEEE Transactions on Geoscience and Remote Sensing, GE-22:1 (1984);

"Future Exploration Techniques for the Coal Mining Industry," by A. Richmond, Queensland Government Mining Journal, Vol. 86, No. 1005, pp. 321-327 (1985);

"Two-Dimensional Velocity Anomaly Reconstruction by Seismic Tomography," N. J. Wattrus, Doctoral Thesis at the University of Minnesota (Duluth) 1984 [Available from Dissertation Abstracts, University Microfilms, Ann Arbor, Mich.];

"Resolution and Error of the Back Projection Technique Algorithm for Geophysical Tomography," by N. R. Burkhard, URCL 52984 (1980)(available from U.S. Nat. Tech. Info. Serv., Springfield, Va.;

"Exploration Technology Development Program of the Board of Industrial Leadership and Development: Summary of Research 1984-1985 by V. G. Milne (ed) and R. B. Barlow (ed) Ontario Geological Survey Miscellaneous Paper 125, Ontario Geological Survey, Ontario, Canada, pp. 37-44 (1985);

"Seismic Tomography of Boreholes," by G. A. McMechan, The Geophysical Journal of the Royal Astronomical Society, Vol. 74, No. 2, pp. 601-612 (1983); and

"A Stable Geotomography Technique for Refractive Media," by R. D. Radcliff, IEEE Transactions Geoscience Remote Sensing (NY), Vol. 22, No. 6, pp. 698-699 (1984).

The reconstruction data found in the above patents and most of the articles are based on the assumption that the energy travels in a straight line from the source of the receiver in the same manner as X-rays do in medical applications (ray-optic). This assumption is not correct, however, for frequencies low enough to penetrate the earth. Previous attempts have been made to try to correct this assumption using geometrical optics, back propogation, and ray tracing methods, however, these corrections have only limited success.

Further, the previous approaches have neglected the effect of wave diffusion. Previous geotomography applications have made measurements of wave amplitude (electromagnetic and acoustic applications) to reconstruct an image of attenuation-per-unit length, or have measured wave velocity (acoustic). These attenuation-per-unit length geotomography reconstructions, although related in a complex way to electric or elastic parameters of the earth, cannot be used to accurately derive such parameters.

It is therefore very desirable to have a non-destructive method for obtaining high resolution information of subsoil structure between two or more boreholes by measuring the amplitude and phase of a continuously transmitted signal, wherein the signal can be received, measured, processed and correlated to useful physical property measurements of the subsoil, particularly dielectric constant and conductivity in the case of electromagnetic signals, and elastic moduli values in the case of seismic signals.

The present invention relates to an improved method for obtaining high resolution information regarding the properties of subsoil structure in the plane between two or more spaced apart boreholes, which method comprises:

placing at least one electromagnetic radiation transmitter in at least one borehole;

placing at least one electromagnetic radiation receiver in at least one of the other boreholes;

activating the at least one of the electromagnetic radiation transmitter using a continuous frequency of 100 Hz or greater;

receiving continuously at the at least one electromagnetic radiation receiver, the electromagnetic radiation signal which originates from the transmitter of the electromagnetic radiation;

processing the reduced electromagnetic radiation signal continuously received at least at one receiver as a result of the transmission of the at least one transmitter; and

measuring both the amplitude and phase of the received electromagnetic signal.

It is the primary object of the present invention to provide a method for obtaining useful data related to subsoil structure.

Another object of this invention is to provide a system for obtaining useful data related to subsoil structure from which a multi-dimensional graph of the subsoil between two or more boreholes can be produced.

Another object of the present invention is to provide a non-destructive monitoring technique to obtain useful data related to subsoil and rock structure and physical properties.

Yet another object of this invention is to provide a method of deriving conductivity and dielectric constant data of subsoil structure.

A further object of the present invention is to provide a method of deriving data related to dynamic sheer modulus and dynamic Young's modulus of subsoil and rock structure.

A further object of the invention is to provide a method of deriving data of subsoil to assist in determining location, volume, and/or concentration of coal deposits, oil shale deposits, petroleum deposits, organic chemical contaminants, inorganic chemical contaminants, moisture, salinity, fractures, joints and the like.

FIG. 1, one part, found in the lower half of the drawing, is a combination of two aspects of the invention in a cross-sectional view of the subsoil structure between two boreholes and is used to explain the basic principles of the invention. This shows an underground gasoline tank leaking into the space between the two boreholes. The second part found in the top half of the drawing is a schematic diagram showing the computer controlled mechanical winches for controlled raising and lowering of both the transmitter and the receiver, and also the electrical and computer components to reduce the electromagnetic radiation data to useful form.

FIG. 2 is a combined table diagram which shows the expected dielectric constant and conductivity values for the patches in the plane between the boreholes of FIG. 1 where gasoline has leaked into the subsoil.

FIG. 3 is a diagram of the scanned area shown in FIG. 2 illustrating at least some of the dimensional relations used in reducing the measured data.

FIGS. 4A-4B is a flow chart of a representative system adaptable to a computer for determining the geometry of the scanned area.

FIGS. 5A-5B is a flow chart of a representative system adaptable to a computer for calculating the basic matrices necessary for solution using among other inputs the background information.

FIGS. 6A-6B is a flow chart of the solution using the measured frequency and phase.

In a report prepared for the U.S. Nuclear Regulatory Commission (NRC), entitled "Indirect Rock Mass Investigations for Optimizing Borehole Drilling Programs: Wave Diffusion Geotomography," Vol. 4, A. Q. Howard, Jr. et al. discuss several theoretical approaches to ascertain the geologic properties and internal geometry of a soil and rock mass. [See NRC Report No. NUREG/CR-3143, Vol. 4, October 1983, L under Contract No. NRC 04-78-269.] Geotomography and wave front reconstruction (holography) were considered two promising techniques for site investigation and characterization. Algorithms were presented for both ray-optic geotometry and diffusion geotomography. At one point the measurement of both wave amplitude and phase, treating the inversion using complex variables with the VECTOR algorithm in theory is discussed. However, the ability to solve the problem, or to determine whether or not the problem was capable of being solved, was limited by the unavailability of an adequate computer program.

In a Final Report entitled "High Resolution Geotechnical Exploration Using Wave-Diffusion Geotomography," by Charles E. Glass, et al., a high resolution method to remotely image rock mass features important to the construction of deep basing facilities is disclosed. [See Small Business Innovation Research (SBIR) Program Phase I, Contract No. F04704-85-C-0121, supported by: USAF Ballistic Missile Office, Norton Air Force Base, Calif., Report submitted March 1986.] The reference discloses that wave-diffusion geotomography furnishes precise images of subsurface rock mass electrical properties, properties which relate directly to the subsurface distribution of ground water.

The complex patterns of both electrical conductivity and dielectric constant can be reconstructed precisely using wave-diffusion geotomography. A number of wave-diffusion algorithms can be implemented using iterative approaches, however, conventional projection approaches are unstable and do not converge. One iterative algorithm which is useful in the present invention is the one published by V. Pan and J. Reif, entitled "Efficient Parallel Solution of Linear Systems," Proceedings, 17th Annual ACM Symposium on the Theory of Computing, Providence, R.I., pp. 143-152, May 1985. This reference and particularly the algorithm disclosed is hereby incorporated by reference.

A portion of the text related to the Pan Reif algorithm described in the aforementioned reference is reproduced hereinbelow: ##SPC1##

With reference to the lower portion of FIG. 1 as described above, the lower portion shows a cross-sectional view of two spaced apart boreholes 11 and 12 having a vertical plane 13 between them. Land surface 14, is between boreholes 11 and 12 is over plane 13, and is the surface on which rests the computer controlled winch 15 for electromagnetic radiation transmitter 16 on cable 17 and for computer controlled winch 18 for electromagnetic radiaton receiver 19 on cable 20.

Multiple electromagnetic radiation transmitters and receivers at different locations have been demonstrated by computer in the present invention to produce the desired conductivity and dielectric constant values. The invention contemplates the use of multiple transmitters 16, 16A, 16B, etc. and multiple receivers 19, 19A, 19B, etc., normally with the number of transmitters 16 being less than the number of receivers 19. Here, multiple transmitters and multiple receivers also means multiple transmitter and multiple receiver locations.

The boreholes 11 and 12 are shown being essentially parallel. However, accurate data are also obtained if the boreholes are skewed, angled, helical and the like, so long as the distance between transmitter 16 and receiver 19 are accurately known. It is to be understood that electromagnetic radiation (or radiation) in this invention include, radiation having wavelengths from about 100 meters to 0.1 meter, and also includes seismic radiation wavelengths in the same range. A preferred range for the wavelength is between about 10 meters and 1 meter. It is noted that for the purpose of this invention that "electromagnetic radiation" or "radiation" also includes wavelengths in the seismic and sonic ranges.

VECTOR-Geotomography Design Requirements--Using Volume Eddy Current Technique of Reconstruction (VECTOR)-geotomography, all transmitters (16, 16A, 16B, etc.) send a signal simultaneously and in phase so that increases in unique receiver data must be acquired by increasing the number of receivers (19, 19A, 19B, etc.) in the receiver borehole 12. The unique information that can be acquired depends upon the probing frequency being used. For probing frequencies of 100 MHz (wavelengths on the order of 1 m for our models) very little new information can be acquired by increasing the receiver density beyond approximately four or five receivers per meter. A preferred number of receivers is about twice the number of patches. Length of receiver borehole 12 and signal attenuation will combine to limit the number of patches 21 possible using the VECTOR algorithm. A maximum depth for the boreholes is less than about 7 miles. A more preferred depth is less than about 5 miles. An especially preferred depth is about 3 km or less. For example, a 100.times.100 patch network will require 20,000 receivers for a redundancy factor (ratio of receivers to patches) of only 2. At a probing frequency of 100 MHz, sampling five times per wavelength requires 4 km of borehole and the ability to measure the scattered field along the entire 4 km. Receiver density can be increased, but even if one could measure every 10 mm along the borehole, 100 m of borehole would be required. It is doubtful that measurements at 10 mm intervals will provide crisp new data even with a probing frequency of 1 GHz. Hence, a rule of thumb for the maximum number of patches for a given VECTOR-geotomography field situation is: ##EQU1## where NP is the number of patches in the reconstruction, BL is the length of receiver borehole in which effective measurements can be made, and .lambda. is the probing wavelength. A borehole 200 meters deep would permit a VECTOR-geotomography reconstruction of approximately 500 patches (22.times.22) at 100 MHz. For boreholes 100 m apart, one could achieve resolutions of between 2 to 5 m depending upon placement of the reconstruction grid.

Transmitter borehole 11 is designed to the same requirements as a linear antenna array, because in the VECTOR-geotomography algorithm all transmitters must effectively work simultaneously. This arises, as discussed below, because the x vector reconstructed in the VECTOR algorithm is not invariant with respect to transmitters; it contains the dependence on .sup..phi. 2 (see below). Consequently, the number of transmitters 16 (16A, etc) and their spatial placement in transmitter borehole 11 determine the width of the effective antenna beam. Signal strength must be great enough to permit the scattered field to irradiate the measurement borehole. Signal power on the order of 100 W should be adequate for most field situations. This is quite a modest requirement considering the systems available today which are capable of generating several KW of power. Best coverage will be obtained when receivers and transmitters are aligned symmetrically with respect to the reconstruction region, but this is not a strict requirement. Receivers must, however, be capable of measuring both signal amplitude and phase.

TENSOR-Geotomograhy Design Requirements--The added iterations required during reconstruction using the TENSOR algorithm are partially offset by advantages in field design. Unlike VECTOR-geotomography, Transmitter-Enhanced Numerical Solution of Reconstruction (TENSOR)-geotomography permits each transmitter 16, or array of transmitters (16, 16A, 16B, etc.) to operate individually much as the typical ray-optic approach does. Hence, under the same conditions as specified above for VECTOR-geotomography, instead of 20,000 receivers (19, 19A, 19B, etc.), TENSOR-geotomography can achieve equal redundancy of 2 using only 150 of both transmitters and receivers. For the same 200 m borehole and frequency of 100 MHz, TENSOR-geotomography needs only 30 m of receiver borehole 12 rather than 4 km. Resolution of TENSOR-geotomography is then less dependent on borehole geometry. Transmitter 16 and receiver 19 design, transmitter power requirements, and receiver measurements are essentially the same as for VECTOR-geotomography.

In the upper portion of FIG. 1 is a schematic electrical diagram describing computer controlled winches 15 and 18 to raise and lower transmitters 16, 16A, etc., and receivers 19, 19A, etc., respectively. In the schematic a signal generator 22 producing the electromagnetic radiation is connected to a power divider 23. In power divider 23, a portion of the signal is conveyed to the network analyzer 24 as a measure of the transmitted electromagnetic radiation. There are also commercial instruments available which incorporate the functions of signal generator 22, power divider 23 and network analyzer 24, and are purchased as network analyzers. The signal then passes through the RF power amplifier 25 and then to the coupler 26. A small portion of the signal is shunted through power meter 27 for storage in computer 28.

The electromagnetic radiation reaches the transmitters 16, 16A via cable 17. Cable 17 is protected by ferrite beads to minimize uphole conduction of transmitter radiation. The electromagnetic radiation signal may be pulsed or continuous and is transmitted from one or more of the transmitters 16, 16A, etc. and is received by one or more receivers 19, 19A, etc. The signal is visually monitored for at least the amplitude and phase of the wave on optional digital scope 29, or visually monitored on the screen of computer 28, and measured in network analyzer 24 for both amplitude and phase. Network analyzer 24 is connected to the signal generator 22 wherein the signal is divided in power divider 29 and also the network analyzer 27. Computer 28 has an input 30 on which data can be manually entered into the computer. The output 31 of the data processing is expected to be the series of values for dielectric constant and conductivity such as are found in the 36 patches as shown in FIG. 2.

Background information regarding the base (background) condition of the soil may be obtained in a number of ways. One method is analyzing samples of the core of the borehole as each is being drilled. Another is to have a third borehole at the edge of the area to be monitored between the active boreholes. The background information (e.g. conductivity, dielectric constant, etc.) is a standard to which is compared the unknown information derived in the present invention.

In a preferred embodiment the following equipment is used to generate the amplitude and phase data regarding the transmitted and received electromagnetic radiation.

Also in FIG. 1 in the lower portion is shown a 6.times.6 patch grid plane 13 having square patches 21 each having a 2-meter side. The surface area has a gasoline tank 40 which is suspected to be leaking into the subsoil. This grid plane is interrogated using the system described hereinabove at about 300 MHz. The output from the data processing using the computer program is expected to generate the values for dielectric constant and conductivity for each of the patches (1-36) shown in FIG. 2. Using this generated information, the concentration of gasoline and water in each patch plane can be reconstructed. The plume areas 41, 42 and 43 are then easily created using a graphics computer program. Area 41 is expected to have the largest concentration of gasoline with the available soil porosity (35%) being filled with about 95% gasoline/2% water. Area 42 is expected to have about 50% gasoline/40% water mix. Area 43 is expected to have about 25% gasoline/32% water mix. The remaining patches including patch 21, are found to have the measured background wherein the porosity of the soil has about a 24% water content and no gasoline.

In this example, either the dielectric constant or the conductivity can be used to derive the gasoline/water mix. In more complex soil environments, both dielectric constant and conductivity used together will prove very useful to discern soil structure and enable reconstruction of contaminant distributions. The movement of the gasoline plume areas 41, 42 and 43 from underground storage tank 40 into the subsoil can be easily monitored over time using this method. Similarly, any type of underground tank containing organic or inorganic chemicals (i.e., contaminants from any impoundment) can be monitored for leaking and extent of contamination by adaptation of this method (i.e., plume monitoring a one-time or multiple monitoring over a period of time).

VECTOR-Geotomography Solution-The solution to the VECTOR-geotomography problem requires the assistance of computer technology principally because of the size of the matrices necessary to obtain the solution. In order to better understand the flow charts included herein as FIGS. 4-6, the diagram FIG. 3 should first be discussed. In FIG. 3 an array of transmitters T1 through T5 (transmitters 16A, 16B, etc. in FIG. 1) is located in bore hole 11 while a corresponding line of receivers R1 through R8 (receivers 19A, 19B, etc. in FIG. 1) is located in a second borehole 12. Between these two bore holes is the area to be scanned. This area may well be the area located under the gasoline storage tank 40 as illustrated in FIG. 1.

It is appropriate to divide the scanned area into patches. In a gross sense, a patch, when represented on a cathode ray tube or other display device, can be compared to a pixel. For convenience sake, the patches hereinafter are numbered in rows starting with the top row (See FIG. 3). The dimension of each patch will be referred to as .DELTA. making a square patch of area .DELTA..sup.2. In the mathematics necessary for the solution of the problem, it is appropriate to convert the square patches into an equivalent circle. The equivalent circle has a radius a equal to 0.564.DELTA.. The distance from the center of each patch XP, YP to the receiver will be referred to with the uppercase R and the subscript of the patch number and the receiver number. Similarly, the distance from an appropriate transmitter to a patch number will be referred to as uppercase T with the subscript of the patch number followed by the subscript of the transmitter number. In referring to the patch numbers in the matrix algebra that follows, the conventional Row i, Column j convention will be used. A particular patch may also be referred to as an uppercase P with the patch number being a subscripted lowercase p. The subscript denoting each receiver will be lowercase n and the subscript denoting the particular transmitter will be a lower case t. The center of each patch will be referred to by its conventional cartesian coordinates as XP.sub.p, YP.sub.p. The geometrical location of each patch will be referred to by its coordinates, with the lower left hand corner of the patch pattern being referred to as X.sub.1, Y.sub.1. Other subscripts and the like will be readily recognized from a study of the flow charts appended hereto.

Referring to FIG. 4A and 4B, the patch geometry and transmitter information is computed. It is necessary to establish parameters which is done in step 100 shown in FIG. 4A. Specifically, the number of rows NROW and number of columns NCOL in the patch geometry along with the lower left hand corner of the image and the individual patch sizes should be entered. From this the calculation done in step 102 is to determine the patch size. In FIG. 3, it is readily apparent that this would result in 36 being the patch size in the sample. Next a loop 104 performs the necessary steps to locate the XP component of the center of each individual patch. This is followed by a similar loop 106 that computes the YP component of the center of each patch. At this point the receiver and transmitter locations are entered at step 108 so that the positioning in the various bore holes can be computed in the loops 110 and 112.

Having this information the dimensons R.sub.pn, T.sub.pt can be computed. Of course the physical location X.sub.1, Y.sub.1 must also be known in relation to the two bore holes 11 and 12. Since computation of R.sub.pn and T.sub.pn is a relaively straightforward geometric problem, the flow chart does not include a method to calculate these figures.

FIGS. 5A and 5B computes basic matrix functions necessary for solution. It is to be understood that the parameters computed in the geometry program just described are carried through into this particular step and the subsequent step. In addition at step 112, the frequency to be utilized is read in MHz along with the background conductivity and background dielectric constant. These two figures can be determined experimentally from the adjacent soil. One possible solution, previously noted, is to take core samples during the drilling of the bore hole to determine the conductivity and the dielectric constant of the soil in the vicinity of the area to be sampled. It should be remembered that one of the uses of this particular method is to determine the conductivity and dielectric constant between two spatially separated boreholes. In particular, the method is useful where it is not practical to drill a borehole (e.g., under or through a contaminant). At step 114 the Hankel function of the first kind of order 0 with a complex argument is represented as H.sub.0. This function H.sub.0 (k.sub.1 R.sub.pn) is computed for each patch to each receiver. In block 114, H.sub.0 is the Hankel function of the first kind of order 0 with a complex argument while k.sub.1 is a complex background wave number and R.sub.pn is the patch to receiver distance.

Once this A.sub.pn matrix is computed the patch interaction matrix G is calculated. This is accomplished in the loop 116. In general, the calculations in loop 116 calculate the interaction between the patches. For clarification sake, the Hankel function appears again in this loop along with the Bessel function which is denoted in equations as J.sub.1. This Bessel function is of the first kind with a complex argument where k.sub.1 is complex in both the Bessel and Hankel equations. Argument a, as noted above, is equal to the equivalent patch circle, i.e., 0.546.DELTA.. Because the arguments are necessarily complex, this technique to determine useful information is called complex geotomography.

Upon completion of the computation of the interaction matrix G, the incident field must be computed at each patch F1. This is accomplished in loop 118 again using the Hankel function H.sub.0. In this case the distance is from the transmitter T.sub.t to the patch P.sub.p and is represented by the term T.sub.pt.

At this point, the program titled VECTOR, and illustrated in the flow charts in FIGS. 6A-6B may be used to solve the conductivity and dielectric constant problems necessary to arrive at the composition of the various patches. Certain constants should be entered either at this time or obtained from the earlier programs illustrated in FIGS. 4A-4B and 5A-5B. Specifically, the previously mentioned equivalent patch circle, .omega., .mu..sub.0, .epsilon..sub.0 which is equal to 8.85.times.10.sup.-12, and k.sub.1 which is equal to ##EQU2## The transmitted and received fields are also entered at this point wherein .PHI..sub.1 (x.sub.j) is the transmitted field and .PHI..sub.2 (x.sub.j) is the received field. This is represented by Y(I) as indicated in FIG. 6A at 120. The vector Y(I) is modified by the factors shown in the right hand side of the equation shown in block 120.

Next the Hankel function matrix A, computed in the previous step, is converted to a square matrix as indicated in step 122 and this square matrix is combined with the modified Y vector that results from the operation in 120 at step 124. To this point, the computation followed well established lines. However, in operation block 126, the Pan-Reif algorithm for inversion of this rather large matrix is utilized. The Pan-Reif inversion process is noted above. The resulting matrix is multiplied with the vector A.sup.T Y to form a vector V(I) which is utilized in the step 128 to compute the total field for each patch Pp. It should be noted that the variable X.sub.n is the nth vector from the V(I) VCOMP vector computed at step 127.

At this point, sufficient information has been obtained to compute in box 130, the dielectric constant .epsilon..sub.2 and the conductivity, in box 132, represented by .sigma..sub.2.

Having computed the conductivity and the dielectric constant for each patch, the characteristics of that patch can be accurately described as indicated above.

VECTOR-GEOTOMOGRAPHY--VECTOR-geotomography is formulated from the integral equation describing the received electromagnetic (EM) field given as

where .phi..sub.1 (x) is the transmitted EM field, g(x,x') is the two-dimensional Green's function corresponding to the point source equation, and j.sub.v (x') is the volume current source within the anomalous region given by ##EQU3## Here k(x') corresponds to the wavenumber and the integral is carried out over the anomalous area A. The idea is to solve Equatin (1) for the anomalous wavenumber and subsequently acquire the dielectric constant .epsilon..sub.r (x') and conductivity .sigma.(x') from k.sub.2 of Equation (2).

Equation (1) can be solved by discretizing the interborehole interrogation area into a series of M square patches having dimension .DELTA. and center of mass .xi..sub.m. By defining a characteristic function .rho..sub.m (x) as ##EQU4## and assuming that the anomalous wavenumber k.sub.2 and field .phi..sub.2 (x) are constant within each patch, Equation (1) becomes ##EQU5## If the square patches are approximated by circles having a radius a, such that the circle area equals .DELTA..sup.2 (a=0.564.DELTA.), the integral in equation (3) can be solved analytically and equation (3) becomes ##EQU6## where R.sub.m is the distance from patch m to the point x. In the borehole, where the field is measured, x=(X, yi)=xj and we have, after rearranging, ##EQU7## By setting the left-hand side of equation (5) (a function of borehole position j only) to Y, the H.sub.o (R.sub.jm k.sub.1) term to A.sub.jm, and the remaining term (containing the unknown wavenumber k.sub.2m) to X.sub.m, we have the equation ##EQU8## which is identical in form to the ray-optic approach. As in ray-optic geotomography, we wish to insure that the system is overdetermined, i.e., that the set of N measurements (j=1, 2, 3, . . . N) is such that N>M. An interpretation of X.sub.m for low frequencies is that X.sub.m .about.i.omega..mu.J.sub.a (.xi..sub.m) where J.sub.a (.xi..sub.m) is an anomalous eddy current in patch .xi..sub.m. Hence, the inverse problem physically becomes one of determining the set of eddy currents ##EQU9## given the set of measurements ##EQU10## Physically this is far more satisfying than ray-optic geotomography because reflection, refraction, and wave diffusion (or diffraction, depending upon frequency) are accommodated. In addition, from the reconstructed anomalous eddy currents, the rock mass dielectric constant (a function of the bulk porosity) and conductivity (a function of the style of pore or fracture interconnection) can be computed as ##EQU11##

In the following section titled "VECTOR Elastic Geotomography" is a discussion of the use of algorithms to obtain useful elastic properties. Preferably, the dynamic elastic modulus and dynamic Young's modulus are obtained. This information is very useful in the determination of a geological fault or fracture using complex geotomography.

The original SBIR solicitation sought imaging techniques capable of yielding information on rock mass elastic properties. VECTOR-geotomography, however, is derived for the reconstruction of scattered electromagnetic wavefields. Nevertheless, we will show in this appendix that the VECTOR, TENSOR, and STAR algorithms can be applied in a straightforward manner to imaging systems employing seismic waves. The elastic wave equation satisfies the following: ##EQU12## where F(x,y,z) is the source density distribution, .phi. is the compressional wave potential, .psi. is the shear wave potential, k.sub..alpha. is the compressional wavenumber given by k.sub..alpha. =.omega./V.sub.p =.omega.[(.lambda.+2.mu.)/.rho.].sup.-1/2, k.sub..beta. is the shear wavenumber given by k.sub..beta. =.omega.[.mu./.rho.].sup.-1/2, .lambda. and .mu. are Lame's constants, .omega. is the angular frequency, and .rho. is the density. Using the Green's function approach we can derive the same equation for the elastic and electromagnetic case, that is equation (1) as ##EQU13##

If now we assume that the region of propagation can be subdivided into cells small enough so that k.sub..alpha.2, k.sub..beta.2, .phi..sub.2, and .psi..sub.2 are constant within each, we have ##EQU14##

Please refer to equation A.3 for definitions of .rho.m and .xi.m. Changing to circular cells of radius "a" to evaluate the integral in equation B.3 one obtains ##EQU15##

Letting the following equalities hold ##EQU16## one arrives at an equation equivalent to equation (6).

If we assume a viscoelastic medium with solid friction damping C (small compared to the elastic moduli) we can compute the following dynamic elastic parameters given the density .rho. ##EQU17##

The correspondence between wavelength and frequency is given by ##EQU18##

While only a few embodiments of the present invention have been shown and described herein, it will be apparent to those skilled in the art that various modifications and changes can be made in the method of obtaining high resolution information regarding the physical properties of subsoil structure in the plane between two or more spaced-apart boreholes without departing from the spirit and scope of the present invention. All such modifications and changes coming within the scope of the appended claims are intended to be covered thereby.