Method and Apparatus Employing Received Independent Magnetic Field Components of a Transmitted Alternating Magnetic Field for Determining Location

The present invention relates in general to method and apparatus for determining relative location of one object relative to another and, more particularly, to a method employing an alternating low frequency magnetic field transmitted from one location and received at another. The received signals are analyzed to determine the relative locations of the transmitting and receiving stations.

Heretofore, a number of prior art methods employing one or more low frequency magnetic field transmitters transmitting to one or more receiving stations have been employed for determining the relative location of one of the stations relative to the other. One such prior art method has been utilized for locating a transmitting magnetic dipole antenna operating at a frequency of about 1 kilohertz in a borehole by using a single axis receiver and triangulation methods from within a tunnel. However, the relatively high frequency of 1000 hertz does not permit direct application to cased boreholes and triangulation is not practical in some situations. This method is described in a publication titled, "Electromagnetic Technique for Locating Boreholes", published in 1978 by the U.S. Dept. of the Interior, Bureau of Mines, Report of Investigation 8302.

Another triangulation method for locating a subterranean drill bit transmitter relative to receiving stations on the surface is described in U.S. Pat. No. 3,828,867 issued Aug. 13, 1974. The problem with this method is that it requires a radiated field which is very difficult to do at low frequencies. If quasistatic fields are employed the location calculation is only correct if the dipole transmitting antenna is exactly vertical and at the surface. The dip angle measured is the compliment of the magnetic field dip, not the direction to the dipole. This gives a correct answer only for the trivial case of the transmitting dipole being at the surface and exactly vertical.

Another method for locating the transmitting station uses a two or three-axis transmitter which transmits to a two or three-axis receiver at an unknown location. The received signals are analyzed to give information on bearing but not range unless the transmitter moment is precisely known. This system is disclosed in U.S. Pat. No. 3,868,565 issued Feb. 25, 1975 and U.S. Pat. No. 3,983,474 issued Sept. 28, 1976. This method is not feasible for many applications because the antenna dimensions for a practical three-axis transmitter are too large and the power requirements are too large for most borehole applications.

Another method using a three-axis transmitting antenna and a three-axis receiving antenna is disclosed in U.S. Pat. No. 4,054,881 issued Oct. 18, 1977. Unlike the similar methods cited above, this method can operate open loop but the phase and amplitude of the transmitted signal must be known and communicated to the receiver by modulating the signal. This requires the use of relatively higher operating frequencies than are generally desirable in underground location applications and performance will be severly limited by unknown attenuation and inhomogeneities in the conducting earth. Another problem is that nine measured parameters and other known parameters are required to determine the relative position and orientation of the three-axis transmitting and three-axis receiving antennas and a three-axis transmitter is not practical in a borehole.

Another prior art system requires a relatively large linear scale dipole and quadrupole antenna loop operating with a three-axis receiving station located in a subterranean burrowing device. Such a system is disclosed in U.S. Pat. No. 3,529,682 issued Sept. 22, 1970. The problem with this system is that the large scale transmitting antenna loops are always cumbersome to deploy and impossible to use in many applications such as underground or where obstacles are present. The single three-axis receiver located in the burrowing device is subject to considerable electric noise induced by motion in the earth's field caused by vibration.

The location method described in U.S. Pat. No. 3,900,878 uses two meter diameter horizontal surface loop antennas to sense the magnetic field signals generated by an underground horizontal loop for mine rescue. The phase of the signals detected by portable and reference receivers are compared to determine the point of maximum phase lead which theoretically occurs directly over the underground transmitter. The accuracy will be limited by field scattering in an inhomogeneous medium at the relatively high frequency used for the communication function. Further limitations are an orthogonal relationship between the transmitting and receiving loops is required and range cannot be determined unless the transmitter moment is known and the signal amplitude is measured by a different method. Therefore, this method will provide the correct horizontal location, but not the depth, only if the transmitter moment is exactly vertical and the geology is horizontally stratified.

Other conventional underground survey methods, particularly in borehole work, typically employ gyroscopes or measurements of the earth's magnetic and gravity field at any given moment. Position is then determined by complicated approximations to give the integral over the path which generally leads to a significant cumulative error. Also, the gyroscopes, in particular, are subject to significant instrument drift and cannot withstand high levels of shock and vibration. The earth magnetic field measurements are strongly affected by local anomalies, such as drill pipes, railroad tracks, etc., which introduce errors. In underground mining, conventional optical survey methods are usually used but large cumulative errors typically result due to difficulties in transferring reference positions from the surface and, often, the inability to close surveys underground.

It is also known from the prior art to employ a low frequency solenoid transmitting antenna preferably with a ferromagnetic core which is particularly suited for use in subterranean locations such as boreholes, tunnels and the like and operating with receiving antennas located on the surface. Such a system is disclosed in U.S. Pat. No. 3,967,201 issued June 29, 1976 for communication, but does not address location.

There is a need to determine the relative position and/or orientation of remote objects under conditions (e.g., underground) where either the transmitter or receiver must be constrained to a single axis device. Also, there is a need to operate at frequencies as low as 1 Hz to minimize the magnetic field distortion due to various ferromagnetic and conductive anomolies resulting from man-made or geological structures. Further, there is a need to determine the position and/or orientation information with a minimum amount of data and without precise knowledge of the magnitude and phase of the transmitted signal or the attenuation information with a minimum amount of data and without precise knowledge of the magnitude and phase of the transmitted signal or the attenuation in the medium between it and the receiver and without the need to move some device in some cases.

The principal object of the present invention is the provision of an improved method and apparatus employing the transmission of alternating magnetic fields through the earth for determining the previously unknown position and orientation of a remote object. The restrictions placed on underground instrumentation by the environment and the confines of boreholes make the following cases the preferred embodiments:

1. Location of a single-axis transmitter, which could be run in a drill string behind the bit and operate reliably while actually drilling. Subclasses of this type of system are:

1.1 Use of three or more triaxial sensors at known locations and orientation to completely determine the transmitter's position and orientation without use of additional information.

1.2 Use of two triaxial sensors as in 1.1 which are adequate for complete determination when certain known ambiguous geometries are avoided.

1.3 Use of various arrays of one, two, and/or three-axis sensors, the total number of sensor elements required being determined to avoid ambiguities in the particular geometry involved. At least six elements (sensor axis) are required.

1.4 Use of less than the number of sensor elements required for unambiguous determination of transmitter position and orientation, plus enough additional information such as transmitter moment, inclination, azimuth and/or distance in the hole, to complete the solution.

2. Location of a triaxial sensor in a borehole with respect to two spaced single-axis transmitting loops on the surface or sea floor.

3. Location of a single-axis transmitter in a borehole with respect to a three-axis sensor by moving the transmitter through the borehole in known increments.

4. Location of a single-axis moving transmitter in a borehole with respect to a three-axis sensor, where the transmitter is attached to a steerable drill moving nominally toward the sensor, and the location measurements are used to guide the drill on a prescribed course.

In one feature of the present invention, alternating magnetic fields transmitted from a transmitter station or stations are received at a receiving station or stations to derive a set of signals containing a plurality of magnetic field components specifying the total complex vector of the a.c. magnetic vector field or fields being received. These received field components are quantified and employed either in a graphical method or in a method of successive approximations to determine the relative locations of the transmitting station or stations relative to the receiving stations.

In another feature of the present invention, the transmitting station employs a single axis magnetic dipole transmitting antenna and the receiving station uses a three-axis receiver in a known orientation and at a known location and is employed for determining the relative position of the transmitting station relative to the receiving station. This arrangement is particularly suited for locating a remote or buried object at the transmitting station.

In another feature of the present invention, the location of a remote or buried receiving station or stations is determined by locating the transmitting station at a known orientation and at a known location and measuring the resulting a.c. magnetic field at the receiving station or stations to derive a relative location of the receiving station or stations relative to the transmitting station.

In another feature of the present invention, one single axis transmitter of unknown magnetic moment or orientation is located by employing a pair of three-axis magnetic field receivers at a pair of known receiver stations. A third known receiver station can be added if the transmitter is expected to be near a locus that produces ambiguities when only two receiver stations are used. This method is particularly useful for borehole surveying, river crossing, drill guidance, pipeline tunnelling guidance, mine surveying and shaft location.

In another feature of the present invention, the location of one single axis receiver is determined relative to the position of two three-axis transmitter stations by analyzing the received a.c. magnetic field and calculating the relative location of the transmitting and receiving stations. This method is particularly suited for location of a trapped miner, who can deploy a wire loop receiving antenna (of unknown size or axis orientation, due to the tunnel size restrictions, floor irregularites, rock falls, etc.). The receiving antenna is connected to a battery powered electronic system capable of receiving and measuring the received field signals and relaying the information to the surface.

In another feature of the present invention, a single axis transmitting station is moved to two or more related but unknown locations, i.e., lowered in known steps down a borehole, and a single three-axis receiving station is employed at a known location for analyzing the received signals and determining the relative location of the transmitting station relative to the receiving location.

In another feature of the present invention, a pair of single axis transmitters at transmitting stations of known location, orientation and relative moment are employed for transmitting low frequency magnetic fields to a three-axis receiving station of unknown location or orientation. The received signals are analyzed to determine the position of the receiving station relative to the known transmitting stations. This scheme is applicable to surveying of deep boreholes by deploying two separated transmitting loops on the surface, or sea floor, capable of producing the very large magnetic moment needed for deep penetration.

In another feature of the present invention, a single-axis transmitting station of known location, orientation and moment is employed for transmitting alternating magnetic field signals to a three-axis receiver station of known orientation and unknown location. The received signals are analyzed to determine the location of the receiving station relative to the transmitting station. This method is particularly useful for deep hole surveying and the orientation of the receiving antenna is conveniently obtained by including a tiltmeter and magnetic compass at the receiving station. This system is particularly useful by using all available real estate for one large transmitting loop for extremely deep holes or restricted access on the surface.

In another feature of the present invention, a two-axis transmitting station of known location and magnetic moment is employed for transmitting an alternating magnetic field to a three-axis receiving station of unknown location or orientation. The received signals are analyzed to determine the relative location of the three-axis receiver sensor, or knowing its location the orientation of the receiver sensor.

In another feature of the present invention, a quasi-static alternating magnetic field or fields is transmitted from essentially only a single axis transmitter movable along a course. A three-axis receiver receives the transmitted magnetic field and a steering signal is derived from the received fields for guiding the course of the transmitter relative to the receiver station.

In another feature of the present invention, quasi-static alternating magnetic fields are transmitted from a plurality of single axis colinear transmitter stations of known separation. The transmitted fields are received at essentially only a single three-axis receiver station of known location and known orientation relative to the transmitter axes. The location of the transmitting stations relative to the receiver station is derived from the received magnetic fields.

Other features and advantages of the present invention will become apparent upon a perusal of the following specification taken in connection with the accompanying drawings wherein:

FIG. 1 is a schematic line diagram including an enlarged detail portion depicting certain parameters associated with a transmitter and receiver location system employing features of the present invention,

FIG. 2 is a plot of received normalized vertical magnetic field intensity Hv vs distance in the Z direction of FIG. 1,

FIG. 3 is a plot similar to that of FIG. 2 depicting the horizontal magnetic field component,

FIG. 4 is a cartesian coordinate plot depicting a magnetic moment M and specifying certain parameters of the resulting magnetic field,

FIG. 5 is a schematic block diagram of a magnetic field transmitter (source) incorporating features of the present invention,

FIG. 6 is a schematic block diagram of a magnetic field receiver (sensor) incorporating features of the present invention,

FIGS. 7-14 (a-c) are schematic line diagrams depicting various arrangements of transmitter and receiver stations for practicing the present invention,

FIG. 15 is a computer program logic flow diagram depicting the logic of a portion of a program employed for practicing the present invention,

FIG. 16 is a computer program logic flow diagram depicting the logic for a second portion of a program used in practicing the present invention,

FIGS. 17A-17B is a computer program logic flow diagram depicting the logic for a gradient calculation subroutine of a program used in practicing the present invention, and

FIG. 18 is a diagram similar to that of FIG. 4 depicting a two pole flux source.

By making appropriate measurements of an electromagnetic field, it is possible to determine the source. In an essentially lossless environment, like the atmosphere, the transmitter location may be determined from coded or uncoded propagating electromagnetic waves. However, in a conducting medium, such as the earth or water, rapid exponential signal attenuation results from energy dissipated through currents generated in the medium. It has been shown in the aforecited U.S. Pat. No. 3,967,201 that low frequency magnetic fields are a very effective means for communicating in a lossy medium. If appropriate total complex vector measurements are made of these magnetic fields, the relative position and orientation of the transmitter and receiver can be calculated.

For the geometry illustrated in FIG. 1, the radial and tangential magnetic field strength (Hr and H0 in polar coordinates) components at the receiver, which result from a small magnetic dipole antenna, with a moment m, operating through a steel casing, with an attenuation A, located in the earth, with absorption parameter .GAMMA., are as follows:

In the location and navigation cases that are considered here, low frequencies F, (below 1 kHz and generally 1-100 Hz) are typically used so that for typical values of earth's conductivity, (0.001-0.1 S/m) the ranges, r, of interest (less than a few 1000 m) are significantly less than the earth's skin depth, .delta., so the absorption in the earth is negligible. Thus, the expressions reduce to those for simple quasi-static magnetic fields, or .GAMMA..sub.r =.GAMMA..sub..theta. =1.

The presence of an intervening metal structure, such as steel casing, can result in appreciable attenuation A. It can be shown that the effect is equivalent to that of a source with the same location but smaller amplitude than without the casing. Because of geometrical considerations, the relative location of the transmitter(s) and receiver(s) can be calculated without any knowledge of the attenuation, or even the dipole moment magnitude, if it is uniform over the paths of interest. Multiple frequencies may be used to resolve the attenuation in those cases where attenuation is nonuniform. By using relatively low operating frequencies, the attenuation can be minimized so that it does not seriously deteriorate system performance.

Field distortion, or scattering, caused by eddy currents induced in conductivity anomolies, such as metal structures, or nonhomogeneous or a stratified medium is an effect encountered as the operating frequency is increased. The location errors that will be produced by this effect are very difficult to predict because of the complexities and possible variations involved. However, these errors are generally small when the signal wavelengths are much greater than the size of the inhomogeneities in the local region. Location errors can also increase as frequency is reduced toward dc because of ferromagnetic concentration effects in nearby steel structures such as rails. However, many common structural steels appear diamagnetic in the presence of even low frequency alternating fields. It has been demonstrated that the location errors are very small, in at least the 2-10 Hz frequency range, for the actual conductors found in mine and drilling situations.

The use of magnetic field measurements to determine physical location is limited by noise introduced from electrical and mechanical sources. Atmospheric magnetic field noise caused by extraterrestial sources and worldwide thunderstorm activity is the basic, but not necessarily the largest component. A portion of this noise may be described by the simple empirical model H.sub.no.sup.2 .perspectiveto.10.sup.-11 /f(A.sup.2 /m.sup.2 Hz) which is the upper bound of noise levels reported by Maxwell and Stone in the IEEE Transactions on Antennas and Propagation, Volume AP-11, No. 3, pp. 339-343, May 1963. This model is valid for horizontal magnetic fields from about 2 to 1000 Hz, in midlatitudes for all seasons, and has been used for many years to optimize and predict the performance of underground communication and location systems. Although the model makes no allowance for the effects of local thunderstorms, if impulsive noise processing is incorporated in the receiver, the model still represents a reasonable worst case and becomes conservative most of the time. This is partially due to the fact, shown by Evans and Griffiths in the IEEE Transactions on Communications, Volume COM-22, No. 4, April 1974, that local thunderstom activity raises the impulsiveness as well as the mean level of the noise, rendering it more susceptible to processing.

To be more accurate, the magnetospheric or "micropulsation" noise that is generally predominate only below 2 Hz, should be added to the atmospheric noise model. However, magnetospheric noise, even though it rises very rapidly below 2 Hz, is only of concern in very special cases since the optimum frequency for most location problems is generally above this value. Also, the low-frequency end of the "1/f' model should be modified to incorporate the Schumann resonances. As illustrated by Evans and Griffiths in the reference cited above, the first three modal peaks at 8, 14, and 20 Hz would be somewhat above the "1/f' model, while the troughs at 5, 11, and 17 Hz would be well below it. It should also be noted that the vertical component of noise is generally considerably smaller than the horizontal component, typically by 10 to 20 dB. Thus, the application of frequency selection, impulse noise processing, and possibly long baseline cancellation, might lower the "1/f' model by 20 dB.

Power line transient noise is a major factor in most industrial applications. For lack of available measurements, this has been modeled by assuming that noise outside the immediate vicinity of 60 Hz is due to an ensemble of random-sized steps in load current occurring at random times. This source drives a magnetic line dipole, representing the infinitely long, narrow loop of circulating current in the offending power line. The resultant magnetic field varies inversely as the square of distance from the line, and has the spectrum of a delta function at .omega..sub.o =2.pi..times.60 plus a noise continuum of spectral density H.sup.2 (.omega.)=K(.omega..sup.2 +.omega..sub.o.sup.2)/(.omega..sup.2 -.omega..sub.o.sup.2). To be complete a similar term should be added for each harmonic.

The remaining significant noise component, that due to random vibration of the receiving antenna (in rotation) in the earth's field, must be considered but is not felt to be a serious problem. Although no attempt has been made to model it, it has been observed and successfully combatted in the field. Even with a very crude system, measurements within a few hundred feet of an operating oil well drill rig showed total noise levels (atmospheric plus power line plus motional) comparable to the "1/f" model in the 1 to 30 Hz region when the antenna was shielded from the wind. Proper design of a shock-mounted, wind shielded magnetic cored search coil array, or simply a very large horizontal loop for vertical field sensing, should avoid motional noise contamination due to wind and vibration. The problem is much more difficult when dealing with measurements from moving platforms such as aircraft. However, an adaptive processing technique has been devised which can be used to compensate for small rotations in such situations.

The location errors are a function of the signal-to-noise ratio (SNR) of the magnetic field strength measurements and the field distortions caused by inhomogeneous conditions. In general, the preferred frequency for location in such that the earth's skin depth is on the order of or greater than the operating path length at the maximum operating frequency to minimize attenuation and distortion effects due to the conducting medium. As the operating frequency is increased significantly above this point, the attenuation in the conducting earth increases rapidly. Although the most common noise conditions tend to decrease with increasing frequency, the rate is lower than the signal reduction. Since the location accuracy is roughly proportional to the SNR, it will degrade in some proportion of the frequency increase for a given transmitting and receiving system (e.g., if a SNR=40 dB results in a 1% location error at f, the SNR might be .ltoreq.20 dB at 10.times.f, which would result in an error >10%).

For communication through the earth, the optimum frequency is essentially that which maximizes the ratio of received signal-to-noise density for the given (or assumed) profile of transmitter moment vs. frequency, propagation loss, and noise model as shown in the aforecited Pat. No. 3,967,201. For location, accuracy is constrained by the same signal-to-noise ratio as for communication, although a much higher value (approximately 40 dB compared to approximately 10 dB) is required. For homogenous earth, therefore, the same optimum frequency will result. In the real world, inhomogeneities and scattering obstacles could modify this choice in a very complicated manner. For instance, in the presence of inhomogeneous formations and nonmagnetic obstacles, the optimum frequency will always be reduced to a value where the unpredictable scattering is comparable to the additive noise. On the other hand, if magnetic obstacles are present, particularly long, thin rods, the frequency may have to be increased to reduce scattering.

For location, the optimum frequency will depend on the magnitude of the transmitting moment as well as on the noise spectrum, since scattering "noise" is multiplicative rather than additive. The optimum frequency for scattering alone may not occur at zero frequency. In practice, however, a frequency of a few hertz is sufficient to greatly reduce ferromagnetic concentration effects, at least in the mines where this technique has been employed.

The specific operating frequency is chosen to minimize the effects of the ambient noise, the conducting medium and inhomogeneities in the medium.

Any signal acquisition method may be used which minimizes noise effects and produces an accurate measurement of the signed signal amplitude. The exact signal phase angle, as distinct from polarity, is often determined to achieve the advantages of synchronous detection but it is not necessary for the calculation of location. One example of a cross-correlation method is described in the following but there are many others.

Suitable processing of noisy data should be employed in order to achieve satisfactory location accuracies for practical applications.

The data processing steps for the problem can be divided into two general categories:

1. from samples of the magnetic field sensor outputs, determine estimates of the amplitude and phase of the sinusoidal field at each of the sensors;

2. given these estimates of the amplitude and phase of each component of the magnetic field, estimate the position, orientation, and source strength of the transmitting antenna.

This discussion considers these problems and presents the derivation of algorithms that will provide the desired estimates.

For purposes of this derivation only one signal will be considered; the process is repeated identically for the remaining five or more signals to be processed.

The signal to be processed is provided as a sequence of amplitude samples taken at times t.sub.n and can be represented in the form

where

v(t.sub.n)=the sample value at time t.sub.n

A sin (.omega.t.sub.n +.phi.)=value of the useful signal at t.sub.n

n(t.sub.n)=additive noise value at t.sub.n

.phi.=phase of desired signal with respect to arbitrary time epoch established by the sampling

.omega.=signal frequency--presumed known.

From these amplitude samples, an estimate for the two parameters A and .phi. can be obtained. This can be done by fitting a sinusoid to the samples, and determining the sinusoid parameters that lead to a minimum squared error fit. (This process is equivalent to synchronous detection followed by an integrating filter). The error of the fit for each sample is expressed by

and the total squared error over the sample set of N samples is ##EQU2##

It is the expression that is to be minimized via the selection of A and .phi.. At this point some manipulation of the terms in this expression will be helpful. Employing the identity:

the second term becomes ##EQU3## Note that this has removed the variable to be estimated from within the summations. Next, the third term can be rewritten as ##EQU4##

Observe that the second term of Eq. (7) does not depend on the measured data and will be identically zero, independent of .phi., if (1) the samples are taken uniformly in time, and (2) the sampling rate is a multiple of four times the signal frequency, and (3) N is selected so that samples are taken over an integer number of signal cycles. Since the signal frequency .omega. is known, these constraints on the signal sampling are easy to meet.

A second point must be observed about Eq. (7): if it is determined that a particular sample is badly contaminated by noise, then it may do more harm than good to include it in the determination of A and .phi.. The only effect on all other terms of the expression for .epsilon..sup.2 except this last one is that the edited out sample is deleted from the indicated summations (in effect V(t.sub.n) is set to zero for that t.sub.n) and the value of N is reduced to reflect the deleted samples. In this case of arbitrarily deleted samples, the last term (the .SIGMA. cos 2 () term) will no longer be identically zero because it will be mssing the edited out terms. This can be corrected by adding these terms to the deficient summation, at the same time subtracting them from the total expression. Thus, for the case where K data samples are deleted, the expression becomes ##EQU5## where t.sub.n are those times when the data samples are deleted. ##EQU6## Written in this form, the two parameters we are to determine, A and .phi., are not contained in any of the summations involving the measured data. Therefore, to simplify the remaining discussion, define the following ##EQU7## then

and

and setting these expressions each to zero yields ##EQU8## and

Substituting the expression for A into this latter equation leads to ##EQU9## Multiplying these expressions and collecting terms leads to the equation

or ##EQU10## Note that this editing process should only be used when the signal is much smaller, i.e., 10% of the noise, in the predetection bandwidth. When applied to noise superposed on larger signals, a bias can result which distorts the apparent signal amplitude. However, this is not a problem in practice since noise editing is not generally needed when signals are, for example, larger than 10% of the noise.

A few comments on how the data can be collected for use in these expressions: Since the samples are uniformly spaced in time, the values for sin .omega.t.sub.n, cos .omega.t.sub.n, sin 2.omega.t.sub.n and cos 2.omega.t.sub.n can be pre-computed and stored in a table for a priori chosen values of t.sub.n. (For example, if 48 samples per cycle are taken, only 12 values need be stored, and the table entered at different points for the sin, cos, and double frequency terms.)

Two data collection approaches are possible (1) Collect and sum into an array of storage registers the samples from corresponding points on each successive sample, and in a second array the number of samples deleted from each sampling point. When the total number of samples (N) have been collected, multiply each sum by its corresponding trig value and sum to obtain the terms a, b, c, d. Alternatively, (2) multiply each sample by the appropriate trig values as the sample is taken, and form the sums a, b, c, d directly. The trade-off is that (1) requires more memory while (2) requires more computing time per sample.

The two general computational categories are:

1. Determination of the locus of a moving object using one single-axis transmitter and one three-axis receiver.

2. Determination of a "fix" of the location and orientation of a stationary object, or the average fix of a slowly moving object during the observation time, using one single-axis transmitter and multiple three-axis receivers on two single-axis transmitters and one three-axis receiver. There are various sub-categories in which only a partial solution is of interest.

The algorithms used for the calculations are described below. Although the examples discussed are for the cases of a single-axis transmitter as the unknown, the same algorithms are useful by exchanging moving for stationary elements or exchanging transmitters for receivers.

In cases where a transmitter of fixed orientation is moved linearly, or if multiple sensors are placed on a common axis, the a.c. field components observed as a function of the linear displacement of the moving (or multiple) item will always be a linear combination of three basic forms, known as the Anderson functions, regardless of the direction of motion or orientation. Two of these functions, normalized to unit moment and unit distance of closest approach, are illustrated in FIGS. 2 and 3.

The process of computing the transmitter locus from data then reduces to correlating data from each axis with these three functions, for various trials of vector moment and closest approach distance and time (position), until maximum correlations are found. The mechanisms of this process require an initial estimate of all parameters, then iteration of trials varying each of the parameters for a best fit.

A relatively simple solution, employing graphical techniques, can be used with satisfactory results when the transmitter orientation is fairly well known (e.g., hanging in a vertical borehole) and conditions are reasonably uniform. Although the applications of the graphical method are fairly limited, it illustrates the location principles and problems involved. For the cases where the receiver is a leveled three axis vector magnetometer as suggested in FIG. 1, the vertical and horizontal components actually measured are related to the field expression in polar corrdinates by ##EQU11## where Z and D are the vertical and horizontal distance, respectively, from the transmitter to receiver and .phi. is the angle of the horizontal field with respect to the receiver coordinates.

These expressions have characteristic patterns which are plotted in FIGS. 2 and 3, normalized to the maximum vertical field at Z=0 which is defined as the transmitter antenna center elevation equal to the receiver elevation. This illustration assumes that the transmitting antenna is moved to generate the data, but similar results would be achieved if the receiver is moved although alignment accuracy will be more critical.

The direction, or bearing azimuth (.phi.), of the transmitter locus from the receiver is given directly by the direction of the horizontal field

This is best evaluated at the horizontal field maximum which occurs when the transmitter is located at Z=.+-.D/2. Quadrant resolution must be determined by observation of the relative polarity of the three field components with respect to the antenna phase.

The distance D, or range, of the transmitter locus from the receiver can be determined in several ways. The value could be obtained from the derivative of the horizontal field at its peak, as inferred above, but the peak is too broad to permit any more than a rough estimate in practice. In those cases where the transmitter can be moved a sufficient distance in the borehole, the range can be determined from the field zeroes, that occur at Z.sub.A, Z.sub.B, Z.sub.0 using the relationships ##EQU12##

An alternate procedure involves calculation of the horizontal field slope at Z=0 when the vertical field zeroes cannot be obtained. In those instances where the borehole does not reach the receiver elevation, it is useful to extrapolate or to use a more involved calculation procedure.

These equations (12-15) are strictly applicable only when the borehole axis, or the path the transmitter follows, is perfectly parallel to the receiver Z axis although reasonable results can be achieved with up to several degrees of misalignment. Also, the simple graphical methods require very good signal-to-noise conditions and an absence of attenuation changes (e.g., due to casing discontinuities) over the transmitter's path to achieve reasonable accuracy.

An iterative calculation procedure can be used to improve the accuracy by correlating data points to find a best fit with the theoretical curve and to automate the entire data reduction process. This method can also be used to solve for the locus of a linear bore of unknown inclination if the transmitter moment and direction remains constant, and a curved bore by assuming it consists of segments of a circular arc. Inclinometers can be used when it is not possible to control the vertical orientation of the transmitter or receiver.

If instead of moving the transmitter in known increments, the received information is used for steering to force the transmitter to remain nearly colinear and nearly on track (as in drilling application), then the locus can be determined from the data after the fact. In this case, the receiver must be located ahead of the transmitter, within the core having an angle of 54.7.degree. with respect to the transmitter axis, to avoid unresolvable ambiguities for steering.

The steering direction to an intercept point at the receiver location can be determined from the ratio of the measured fields, appropriate for the geometry, as suggested by Equation 14. Although some combinations of offset and orientation will cancel at a given range, the process will converge as drilling progresses. If the receiver is offset from the planned locus, but still within the core of unambiguous operation, or a specific steering radius is required, the steering direction to stay on track can be determined by a combination of the received fields and an independent measurement of transmitter inclination and/or azimuth depending on the geometry.

Transmitter location can be determined from calculations using measurements of the total complex vector field from at least two stationary receiver (or transmitter) stations. Although the computations are very involved, particularly when dealing with low quality data, they can readily be done by using integrated circuit microprocessors in equipment that is designed for easy field use.

The method is based on the fact that the vector properties of the magnetic field generated by a dipole transmitter can be used to locate the transmitter directly, without requiring any sort of search pattern. Assuming a rectangular coordinate frame, and a transmitter within that frame located at x, y, z, and generating a magnetic dipole moment with components m.sub.x, m.sub.y, m.sub.z, then the components of the magnetic field vector at any point 1 in the space can be written

where the functions f.sub.1x, f.sub.1y, f.sub.1z are determined by the source and propagation models chosen.

Measuring (receiving, detecting and quantifying) the magnetic field in two distinct locations or from two sources generates two such sets of equations, resulting in six simultaneous equations in the six unknowns x, y, z, m.sub.x, m.sub.y, m.sub.z. In principle, this set could then be solved to yield not only the position of the transmitting antenna but its orientation and signal strength as well. Because the functions f are highly nonlinear, the six simultaneous equations cannot be solved directly. A number of iterative solution techniques can be used to solve this set of equations that are suitable for use by a programmable calculator. In each case the computation begins with an initial estimate of the solution, which need not be accurate, and then converges toward the correct solution of the given set of equations.

The primary advantage of this approach is that it gives both the location and orientation of the transmitter with a minimum number of measurements. More than two measurements may be desirable under poor conditions in order to minimize uncertainty by doing a best fit analysis or to eliminate ambiguities for certain loci of transmitter position.

One method of calculation is described in the following which is a form of sucessive approximations called the Gradient Projection Method.

As inferred in Equation 1 and 2, a dipole source model of the magnetic field transmitting antenna is often used for computation in applications where the distances to the sensor are large compared to the dimensions of the transmitting antenna. Referring now to FIG. 4, given a magnetic source specified by the six parameters of dipole moment (m) azimuth (.alpha.), elevation (.beta.) and position (x,y,z), the components of alternating magnetic field amplitude (H) produced at n locations P.sub.i (x.sub.i, y.sub.i, z.sub.i), expressed in rectilinear coordinates, are

where

Therefore, for the case of n.gtoreq.2 sensors there will be 3n equations of the form

where

i=1, 2, 3, . . . n

j=x, y, z

Note that f.sub.ij, while all of the same form, will be parametized by x.sub.i, y.sub.i, z.sub.i. If a fit error defined as the sum square error of each component is established at each sensor then ##EQU13## where H.sub.ij is the measured value at the i.sup.th sensor in the j.sup.th direction and f.sub.ij is the predicted value for a given set of six source parameters. (A weighting function can be added to the fit error by multiplying each term of the summation by an arbitrarily selected constant A.sub.ij.) This can be minimized by simultaneously setting to zero

and solving for m, .alpha., .beta., x, y, z. This takes the form ##EQU14## where the V is the variable and similarly for the other 5 variables. If the trial solution for the estimated source location is correct b.sub.v will=0 otherwise it will not. This yields a set of six equations whose solution is known to exist but which may not necessarily be unique. Because they are also highly non-linear, a straight-forward solution of this set of six equations is not possible and an iterative solution is necessary.

The partial derivatives of f.sub.ij are of the form

The Gradient Projection method of sucessive approximations is applied directly to the expression for .epsilon.. Constraints, based on a priori knowledge of where the solution cannot lie, may be incorporated to limit the field of search. The solution begins with a trial solution at the point X=(m, .alpha., .beta., x, y, z). At this point both the value of .epsilon.(m, .alpha., .beta., x, y, z) and each of its first partial derivatives ##EQU16## Unless these first partials are all equal to zero, this implies that changing the variable m . . . z of the trial solution will produce a smaller value of .epsilon.; closer to the desired result of .epsilon.=minimum.

The question remains in which direction, and by how much the value of each variable should be changed. The six first partials, when considered as the components of a vector in six dimensional space, constitute the "gradient" vector. This is the vector that is perpendicular to the surface of .epsilon.=.epsilon..sub.x at the point X, and hence points in the direction of maximum rate of change of the function .epsilon. for a change in its independent variables m, .alpha., .beta., x, y, z. Therefore, an iterative procedure that will force the independent variables to the values corresponding to .epsilon.=minimum, is ##EQU17## where the gain parameter k permits control of the convergence and stability of the algorithm. Convergence of the procedure can be detected either by noting that each of the gradient vector components is acceptably close to zero, or by noting that the value of .epsilon. is changing by only an infinitesimal amount on each iteration indicating its arrival at its asymptotic value.

State space constraints can be easily incorporated into the gradient projection method: A constraint prohibiting the solution from lying in some area of the state space can be specified in the form of an inequality like

and similarly for additional constraints. Examples might be m.ltoreq.5000 or 300.ltoreq.x.ltoreq.500. These constraints are incorporated into the algorithm as follows: At each step, each constraint inequality is evaluated from the trial solution at that point by answering the question of whether

For those constraint inequalities that are satisfied, no further action is required. If, however, a constraint inequality is not satisfied, this indicates that the trial solution has moved into a prohibited area.

The procedure then is as follows:

1. Reverse the last iteration to return to an allowed area of state space.

2. Using a very much reduced value of the gain constant k, (say 1%) slowly approach the constraint surface as close as possible.

3. Once the value is as close as possible (in theory right at the point where the constraint equality is satisfied--i.e., lying right on the constraint surface) calculate the projection of the gradient vector tangent to the constraint surface. The next adjustment of the trial solution variables will be determined by this projected gradient rather than the gradient itself.

The projection of the gradient of the function

.epsilon.(m, .alpha., .beta., x, y, z)

that is tangent to the surface

f.sub.1 (m, .alpha., .beta., x, y, z)=0

at the point

X(m, .alpha., .beta., x, y, z)

is found as follows:

Denote the gradient of the constraint surface 1 by the column vector N.sub.1 ##EQU18## and the six by six matrix Q by

where I is the six by six identity matrix and the matrix N is composed of those column vectors N.sub.i corresponding to the constraint surfaces we have hit. (In most cases the matrix N will have only one column corresponding to only one constraint surface.)

The desired gradient projection vector is then given by ##EQU19##

The aforedescribed computational steps are most conveniently performed by a programmed computer.

Referring now to FIG. 18, there is shown an alternative method using two magnetic poles (plus and minus) instead of a dipole transmitter antenna model for computation of location. This method reduces errors at a short range. The model of FIG. 18 and the method of using the two magnetic poles instead of the dipole is described as follows:

Assume a two-pole flux source located symmetrically about the point x, y, z such that the displacement of the + pole from the point is .DELTA.x, .DELTA.y, .DELTA.z. See FIG. 18.

The magnetic flux source is of magnitude .phi. (webers) and is assumed to radiate isotropically. The pole separation is 2D, a known quantity, and leads to the expression

The magnetic flux density at any point P.sub.i (x.sub.i, y.sub.i, z.sub.i,) due to each isotropic pole of this transmitter source will be

where

and the total magnetic flux density at the point will be the vector sum of that due to each pole of the transmitter.

Since the magnetic flux vector from each pole is oriented radially with respect to the pole location, the components of each flux vector at point P.sub.i are given by ##EQU20## and therefore the net x component of magnetic flux density is ##EQU21## with similar expressions for the y and z components. Using these three expressions, and given a set of values for .phi., x, y, z, .DELTA.x, .DELTA.y, .DELTA.z, it is possible to compute the three-component outputs of a magnetic flux receiver sensor 22 placed at point P.sub.i with its three-component axes aligned with our original coordinate axis system. Furthermore, using two or more such receiver sensors 22, resulting in six equations, and the expressions relating .DELTA.x, .DELTA.y, .DELTA.z, to D as the seventh equation, if it is possible to solve the system of equations for the variables .phi., x, y, z, .DELTA.x, .DELTA.y, .DELTA.z, given the six receiver sensor readings. This, however, is not possible in closed form due to the complexity of the six receiver sensor equations but may be done by using an iterative technique such as discussed above.

One particular example of the computer programming steps required to implement the algorithm for calculating the location and orientation of an oscillating magnetic dipole field source is illustrated in the Software Flow Diagram, FIGS. 15, 16, and 17, and the computer program listing in Fortran language for use in a model PDP-11 computer is found in Appendix A. The Initialization Routine permits the entry of a number of values which define the system geometry and control parameters as follows:

1. The receiving system configuration

(a) Number of sensors

(b) Location of each sensor (X, Y, Z)

(This case assumes the orientation of all receivers is parallel to the reference coordinate system.)

2. A best estimate for the transmitter, or signal source, which provides a starting point for the calculation but need not be accurate.

(a) Location of the transmitter center (X, Y, Z)

(b) Dipole vector azimuth (A.sub.z or .alpha.,)

(c) Dipole vector elevation (El or .beta.,)

(d) Magnetic dipole moment (m).

3. Calculation control parameters which relate to the desired accuracy of the result

(a) Minimum signal-to-noise ratio (SNR) to be reached by the processed measured data before the computation of source location and orientation is attempted.

(b) An error threshold above which the routine computing location and orientation will continue to iterate to achieve an answer of sufficient accuracy.

(Time limits can be added to (a) and (b) in the event excessive noise or other error sources prevent a satisfactory solution.)

Once the program is initialized, the operator is able to command a Verify Routine which allows the review and modification of previously entered data. The program may also be restarted or terminated at this point.

Upon entry of a Continue command, the program initializes information for a Data Acquisition Routine. This information includes the number of sensors and previously determined information such as system gains, number of samples per signal period and the number of signal periods to be measured. The Data Acquisition Routine then accumulates data as instructed and processes it to compute signal amplitude, phase and signal-to-noise ratio. (Note that manual entry of data is an alternative procedure.) If the SNR is below the specified level, additional data is accumulated and added to the previous data. When the SNR threshold is reached, the process of computing the transmitting source location and orientation can begin.

The location and orientation calculation is an iterative process where one or more of the six source parameters (X, Y, Z, Az, El, m) are adjusted on each iteration. An output of the computation, called EPSILON, is an error function which is a measure of how well a set of predicted measurements computed from the estimated source parameters fit the measured data. This process is repeated until the error function is reduced to the specified level. When the iteration process is completed, the computed "best fit" source location and orientation is displayed. A tolerance on each source parameter can also be determined as a function of the final error value and geometry. At this point, the program re-enters a command mode where the present data may be processed with revised limits or a new set of data may be taken and processed.

The computation subroutine called on each iteration is called Gradient Calculation, or GRDCAL. The purpose of this program is to calculate an error value as a function of the measured data and a data estimate and to determine the slope of the error with respect to each of the six source parameters X, Y, Z, Az, El and M. Changes to the source parameters can then be computed to improve the source parameter estimates and thus reduce the error function. (All references to source parameters below are to be interpreted as "the latest best estimate".)

When this subroutine is called, it first computes the sin and cos functions of the source azimuth and elevation from the source data (either first best estimate or last iteration) and frequently used products of these trigonometric functions. The X, Y and Z components of the field vector are computed from the source amplitude, azimuth and elevation and the 6 error slope accumulators are cleared.

A "DOLOOP" is then entered which is executed once for each sensor in the system.

First the X, Y, and Z components of Range are computed for the sensor from the source, then the Range (R) and other frequently used values are computed. Computed field values for each axis are then determined (FX, FY, FZ). These values represent the fields that would be measured if the estimate of the source parameters were in fact the actual parameters.

The first derivatives of each field component with respect to each of the six source parameters are computed, followed by computing the differences in measured fields and computed fields for each axis (HXMFX, HYMFY, HZMFZ). The contributions of that sensor to the six error slopes are determined and added to the error slope accumulators.

The error value for that sensor is determined and stored.

Once the above process is repeated for each sensor, Delta values can be computed for each source parameter. This process is both computational and analytical.

First the total field error is computed by summing the field error of each sensor. This value is then compared to the field error from the last iteration. If the new error is larger, the new error is replaced by the old error times a factor where the factor is less than 1.0 to prevent error oscillation and thus solution oscillation. The error slopes for elevation and azimuth are then weighted to make them comparable to the slopes for X, Y, and Z in terms of relative contribution to the error.

Five of the error slopes (magnitude is not included) are compared to determine the highest slope (absolute magnitude). Each of the five slopes are compared to the maximum slope.

The delta values for those slopes less than the maximum slope by some factor are set to 0. The delta values for the remaining slopes are, in this particular example, set as follows:

where the error (E) is raised to a power P. This has the effect of reducing the change as the error approaches 0 Since E is less than 1.0.

Setting the .DELTA.V's to zero for relatively small slopes thus prevents large changes in those parameters since they contribute only a small portion of the total error. However, using a factor in the slope comparisons accomodates inaccuracies in weighting the error slopes for azimuth and elevation relative to those for X, Y, and Z, since they vary as a function of geometry.

Each .DELTA.V is then multiplies by its respective gain and divided by the number of nonzero deltas. Dividing by the number of nonzero deltas simply reduces the gain when more than one change is made, preventing large overshoots.

Changes to the source magnitude are determined independent of the other five parameters and are simply proportional to the error slope for magnitude.

In FIG. 17 of the computer logic flow diagram, lines of the program listing of Appendix A are related to the various logic functions and the respective algorithm is also identified by equation number and related to the respective logic function.

A specific example of a solution of the field equations of Eqs. (17) for the case of FIG. 7 using the computer program of Appendix A is given in Apprendix B. In this example, the program of Appendix A and FIGS. 15-17 is entered at the point of the functional box labeled "input control parameters" of FIG. 15 at which point the real data from the sensors, etc., is inputted manually and the program bypasses that portion of the flow diagram intervening between that functional point and "CALL GRDCAL" of FIG. 16 so that, the program listing picks up again at "CALL GRDCAL".

In the example shown, the actual position orientation and moment of the unknown transmitter are given by:

x=-10.000

y=5.000

z=36.000

A=4.000

E=3.000

M=1,000,000

The fields entered as data are those calculated from these parameters.

The initial trial parameters given the computer are:

x=+10

y=+10

z=+40

A=+10

E=-3

M=900,000

As can be seen from the listing of successive approximations, after several tens of iterations, the solution is rapidly convergent to the correct answer.

In the present invention, as thus far described, the field equations of Eqs. (16, 17, etc.) have been solved for relative position (x,y,z) coordinates of the transmitter relative to the receiver. It is to be understood that these equations may also be solved for orientation of the receiver axes or transmitter axes, if unknown. There are important cases where orientation will be the unknown and the method and apparatus described herein is useful for determining the unknown orientation parameter.

The location equipment of the present invention includes one or more transmitters 11 of the type shown in FIG. 5 which serve as the source of the oscillating magnetic field H, and one or more receivers 12 shown in FIG. 6 which serve to measure the resulting magnetic fields.

The transmitter 11 consists of a magnetic field antenna 16 along with the electronics and power source to drive it. The antenna 16 may be an air core loop 16a, e.g., in planar form or a ferromagnetic core loop 16, e.g., in solenoid form as described in the aforecited U.S. Pat. No. 3,967,201. These loops each produce a magnetic field pattern equivalent to those for an infinitesimal dipole at distances that are large compared to the antenna dimensions. A measure of the transmitter source field strength and field direction is given by the equivalent magnetic dipole moment in terms of Ampere.multidot.Turn.multidot.Meter.sup.2 and a vector direction normal to the plane of the loop. In some cases, a slightly different field model is useful such as two opposite monopoles separated by a distance, d, when working relatively close to a long solenoid (see FIG. 18). The primary advantage of the ferromagnetic core antenna 16 is that it produces a relatively large dipole moment in a compact solenoid configuration that is ideal for borehole use, and it is the only practical way to generate sufficiently large signals in the majority of underground applications.

The transmitter 11 includes an oscillator 13 which, in a preferred embodiment, includes a crystal for precisely controlling the operating frequency to parts per million per degree centigrade. The oscillator output frequency is preferably within the range of 1 to 1000 Hz and its output is fed to a high efficiency transmitter drive 14, the output of which is coupled via coupler 15 into transmitting antenna 16. The inductive antenna is driven with an alternating current by the high efficiency electronics 14 that may be of the direct switching, resonant capacitor storage or other type depending on the application. The efficiency can be optimized by tailoring the drive voltage and current, the winding loss, the ferromagnetic core loss, and the other electronic loss to suit a given application.

One or more antennas 16 may be coupled to the output of the coupler 15. The antennas 16 may have their axes oriented in a certain relationship, i.e., orthogonal, depending upon the particular circumstances of the location system being employed.

In some embodiments, it is important to know more information about the transmitter antenna 16. Therefore, certain accessory navigation components, i.e., tiltmeters, magnetic compass, etc., and communication functions may be employed with the transmitter to provide that such accessory information is either superimposed upon the transmitted magnetic fields as by modulation or the accessory information may be fed from the respective function generators via a separate communication link to the receivers 12 or to a central location where the location computations are to be performed.

In one embodiment, the accessory navigation of uplinked communication functions include a two or three-axis tiltmeter, the outputs of which are fed into a multiplexer and A to D converter 17 and thence into the input of a modulator 18 which serves to modulate the phase of the transmitted oscillating magnetic field. Other monitoring and logging functions such as temperature sensing, static magnetic (earth) field direction, etc., may have their inputs fed into the multiplexer and A to D converter 17 and similarly modulated onto the transmitted magnetic field. (This information can be communicated directly when a wireline is used to drive the antenna.) The information may be digitized and encoded in a pseudo noise format for efficient communication under poor signal to noise conditions.

Batteries 19 may be employed at the transmitter station for suplying power to the various power consuming units or, as an alternative, power may be supplied from a central source via a wire or cable. The transmitting antenna 16 may also be operated by wireline from some remote location.

In addition, an auxiliary down-link command or communication function may be performed by means of a receiver 21 coupled to the antenna 16 via coupler 15 and the output of the receiver being fed to a decoder 20 which decodes the information received by the antenna. The decoder outputs the decoded information to the appropriate circuit functions.

Referring now to FIG. 6, there is shown the receiver 12. The receiver 12 includes the basic magnetic field sensor 22 such as a one to three-axis magnetic field vector sensor 22. Such a sensor may comprise, for example, a magnetometer of some type such as a flux gate, superconducting etc., or a search coil, i.e., loops of conductors. In addition, the magnetic field sensor 22 may include a one to three-axis tiltmeter coupled thereto and tiltable with the field sensor 22 so as to derive the orientation of the field sensor relative to the earth's gravitational field.

The output of the magnetic field sensor 22 is fed to a wide band filter and amplifier 23 for filtering and amplification and the output of the amplifier is fed to a multiplexer and analog-to-digital converter 24 for quantifying the tiltmeter readings, if present. The transmitter 11 and receiver 12 operating frequency stability is controlled to permit coherent (amplitude and phase) detection under noisy signal conditions.

The digitized and multiplexed outputs from the multiplexer 24 are fed to the input of a data transmission system such as a wireline, radio frequency transmitter or the like 26 for transmission to the data reduction equipment including a digital processor 27. In some applications, more than one set of magnetic field sensors 22 are employed in which case the outputs from the other sensors 22 are fed into the digital processor 27 via a digital multiplexer 28 which multiplexes all information derived from the various sensors. In a typical example, the digital processor 27 includes a microprocessor, APU, ROM, RAM, etc. A keyboard and display 29 is interactively coupled with the digital processor 27 and a printer or graphic plotter 31 is interactively coupled with the keyboard and display 29 for printing out and plotting data derived from the keyboard 29 and digital processor 27.

In the preferred form, the weak fields measured by the sensor 22 undergo some filtering, amplification (with automatic gain control to track varying conditions) and digitizing right at the remote detector location prior to communication with the main portion of the receiver or data reduction equipment. After transmission, the various channels of information are multiplexed into a digital processor 27 such as a microprocessor based computing system. The important signal processing functions of impulse noise editing, filtering the signals to suitable quality, noise evaluation, coherent detection to determine phase and amplitude, and decoding any information transmitted are then performed digitally from software programs prior to doing the location calculations. The location calculations are performed after the receiver (or transmitter) location and alignment reference information is entered with the keyboard 29 or alternatively in fixed installations, with a programmable read only memory (PROM). The results are presented on the display 29, printer 31 or graphics plotter 31.

In some applications it may be desirable to use a microprocessor right at the remote sensor(s) for doing all of the signal processing prior to transmission to the central computer 27 (or for doing the entire location calculation) and transmitting only the results when communication channel capacity is limited. The basic receiver 12 can also be implemented by using conventional hardware techniques for the filtering and coherent detection functions, etc., and the location calculations can be done using a programmable hand calculator.

In applications that only require a single sensor, such as for steering or location by the moving transmitter method, the data reduction functions are often included in the sensor data unit. The computed steering or location information may then be transmitted to a remote display 32.

Either the transmitter(s) or the receiver(s) may be the unknown (u) with the other being the known (k). Because of its rugged nature in one particular form, one single axis transmitting antenna 16 is usually used at the unknown location (u), e.g., in a borehole or at a drill bit. The receiver 12 is usually located for direct operator access since it is easier to very precisely align the relatively small receiver sensors and to limit noise induced by mechanical motion of the sensor 22 in the earth's field. However, the roles can be reversed, through appropriate system design, to suit the needs of specific applications.

In general, an array of some appropriate combination of transmitter and receiver elements 11 and 12 is used to obtain enough information for calculating the relative orientation and position of the array elements. For example, the two three-axis receivers 12, of known orientation and location (k), may be used to determine the orientation and location of an unknown single axis transmitter T(u). Other combinations are possible as described below. In all methods, except for steering only, the best accuracy will be achieved when the span between the known array elements (e.g., receiver or transmitter stations) is on the order of the distance between the known array and the unknown. However, it is usually possible to obtain very satisfactory results with much shorter spans, although higher signal to noise ratios and better reference alignment will be required in such cases.

In the basic application, one or more three-axis sensors 22 are used to measure total vector field and are positioned and aligned with respect to some surveyed reference location. This may be accomplished manually by using a transit or automatically by using techniques described below. Two or more sensors 22 may also be pre-mounted on a rigid beam, when the practical spans are adequate, to maintain accurate relative alignment or simplify field set-up. This also facilitates the use of two or more sensors 22 as the unknowns R(u) in inaccessable locations.

Various features can be incorporated to simplify or automate the alignment process or to aid the location calculation process. For example, sensitive two or three-axis tiltmeters can be used to eliminate the need for manual leveling, and the information used to calculate vertical orientation with respect to any set of reference coordinates. In those limited cases where local magnetic anomalies are very small, the earth's magnetic field can be used to determine azimuth for automatic correction of rotation. In both of these cases it is still desirable to accurately measure position when the receiver 12 is the known reference R(k). There are many situations when it is not practical to determine receiver position or orientation by conventional means, such as when it is necessary to deploy receivers at the bottom of a waterway. In these cases, two or more reference transmitters T1(k), T2(k) can be installed at a convenient and precisely surveyed location or locations which will permit automatic calculation of receiver position and orientation. The receiver(s) can then be located as the unknown then used to determine the position of another, unknown transmitter. In general, the orientation of the transmitting antenna's major axis (loop axis) and the antenna location can be determined from the field measurements. However, the rotation of the antenna about its major axis (commonly referred to as the tool face angle in drilling) cannot be directly determined. Except in near vertical holes, rotation may be determined when necessary by adding a two or three-axis tiltmeter to measure orientation with respect to vertical and transmitting this information to the receiver 12 by means described above. Integral tiltmeters can also be used as a cnvenient accessory in certain applications such as aiding the calculation of the unknown transmitter location under very noisy conditions or for determining vertical profile with a single sensor or in aligning the transmitter when it is used as a known reference T(k). Co-located two and three-axis transmitting dipoles are an alternative method for determining rotation as well as being useful in other applications described below.

Location system accuracy is dependent on a number of mechanical and electrical factors. The mechanical factors involved primarily have to do with the alignment and location survey precision of the basic reference used whether it's the receiver(s) or transmitter(s) in the particular application. Mechanical factors can generally be controlled to achieve sufficiently accurate measurements, in practical applications, unless initial hardware cost is a major consideration in which case performance tradeoffs may be necessary. The electrical factors involved are much more diverse and uncontrollable. The effects of attenuation, scattering in a highly inhomogenous conducting medium, distortion by magnetic materials and mechanically induced noise have already been mentioned. Magnetic noise from extraterrestrial sources, worldwide thunderstorm activity and electric power line transients is also one of the major functions which limit system performance. The effects of these factors are minimized through suitable choice of operating frequency, signal processing methods and control or measurement of the reference set-up.

In general, six independent items of information (signal measurements) are required to define the six unknowns of location (3 degrees of freedom) orientation (2 degrees of freedom), and magnitude (1 degree of freedom) of the transmitter or receiver. The absolute phase of the alternating transmitter magnetic moment, although it may be unknown to the receiver, is not an additional unknown in the solution of the location problem as long as the relative phases of all received magnetic field components are measured.

Of the various cases described below, some are capable of the complete solution for six unknowns. Others, as noted, require either specification of some parameters (e.g. moment, one component of transmitter orientation, one component of transmitter position) or some constraint, such as the assumption that the axis orientation of a moving transmitter, although unknown, is invariant or follows a known constraint such as a circular arc. Some cases may be overspecified, that is, more measurements are obtained than there are degrees of freedom. In these cases, a best fit to the known constraints should be applied to make use of all available data, since the measurements are known to be noisy. Even in cases that are not overspecified (measurements match unknowns), due to the inevitable errors in data, a best fit approach is employed rather than seeking an exact solution, which may not exist.

Referring now to FIGS. 7-9, there is shown a class of cases utilizing one single axis transmitter and plural three-axis receivers. This is an important case, as it is a method useful in borehole surveying, river crossing drill guidance, pipeline tunneling guidance, mine surveying, and shaft location. The solution is generally completely specified by the data, meaning, that the six measurements of field (3 axes.times.2 locations) allow simultaneous solution of six equations (three for each receiver) of the type of Eqs. (17) for all six unknowns by the field calculation method described herein.

Only a single axis (magnetic dipole) transmitter is required. This is particularly advantageous in borehole applications, where geometry is limited to one large (axial) and two small (radial) dimensions. Under this constraint, a ferromagnetic-cored solenoid 16 is capable of generating a much larger magnetic moment in the axial direction than in the radial directions, and a very much larger axial moment than any nonmagnetic-cored coil 16a, for the same power input.

The magnitude of the transmitter magnetic moment need not be known a priori to complete the solution, which includes computation of the moment from the data. This is important for use in cased holes, where the effective moment "seen" outside the casing by the receiver 12 is an attenuated version of the moment actually produced by the transmitter 11 inside the casing. This attenuation is a function of the diameter, thickness, electrical conductivity, and magnetic permeability of the casing, and of the frequency and position of the transmitter inside it. In some applications, such as a transmitter fixed inside a drill head which will be beyond the casing in operation, the attenuation can be calibrated by measurement and then considered invariant. In other cases, such as surveying an already cased hole, the attenuation is not only unpredictable but may vary along the length of the hole due to multiple casings, different materials used in different sections, corrosion, etc.

The orientation of the transmitter magnetic moment need not be known a priori to complete the solution, which includes computation of this orientation from the data. This is important in cased holes, where the effective moment orientation may be distorted by non-uniform attenuation, or by the position of the transmitter within the casing. It also provides a by-product of direct measurement of the drill orientation in casing-free operation.

Computation for this general case requires knowledge of the (relative) position of the two receiver sensors 22, and the orientation of their axes with respect to their baseline (the line connecting the centers of the sensors) and one other independent direction, such as the vertical, which together define a sensor coordinate base. The position and orientation of the transmitter are then calculated in this sensor (receiver) coordinate base; to the extent that the sensor coordinate base origin and orientation are known in some other reference coordinates, to within the additional accuracy of the location system. This accuracy is in turn determined by the magnitudes of the moment and the local noise level; the time available for measurement; the geometry (spacing between sensors 22 and from each sensor 22 to the source 16); and the accuracy with which the gain and orientation (within the sensor coordinate base) of each sensor axis are known.

As can be seen by reference to FIG. 7, the location of the transmitter T may be unknown (u) and computed, given the known location (k) of the receivers R.sub.1 and R.sub.2.

Alternatively, in FIG. 8, the location of the receivers R.sub.1 and R.sub.2 may be computed from knowledge of the known position (k) of the transmitter T and measurements of the magnetic field components derived from the receivers R.sub.1 and R.sub.2 and the coordinate system defined by R.sub.1 and R.sub.2. The reciprocal arrangement of FIG. 8 is also useful, i.e., the location of the T(u) can be computed from data derived from the known location of the receivers R.sub.1 and R.sub.2.

This latter case is depicted in FIG. 9 as well, and is useful for blind surveying where only approximate location data is acceptable.

There are certain loci of ambiguity, such as when the transmitter axis lies (1) in a plane bisecting the baseline between the sensor 22 or (2) tangent to any circle whose axis is the baseline between sensors 22, where there is insufficient information to obtain a solution. In those cases where it is not possible to control the geometry to avoid the problem, the required additional information can be obtained in several ways. Independent measurements of one or more components of transmitter orientation and some other constraint (e.g., moment or depth in the hole) can be made in some application such as drilling. Alternatively, a third sensor can be used when the transmitter must be kept extremely simple such as for trapped miner location. Similar techniques are also useful when the signals are contaminated by extremely high noise levels.

Referring now to FIG. 10, there is shown a class of cases which can be considered the reciprocal of the case for FIGS. 7-9 and, more particularly, a class of cases wherein one single axis receiver is utilized with a pair of three-axis transmitters. This example is particularly useful in mining and relies upon the fact that the position of the first transmitter T.sub.1 is known relative to the position of the second transmitter (T.sub.2) to define a coordinate system for determining the relative location of a receiver at an unknown (U) location. Conversely, if the position of the receiver is known (K) the position of the transmitters relative to the receiver may be calculated, again assuming that the position of one transmitter relative to the other is known such as the transmitter being co-located on a single beam or structure, the axes of the first and second transmitters being known relative to each other.

By reciprocity (exchange of transmitter and receiver), all the characteristics denoted in the case of FIGS. 7-9 above apply. This present case has potential application to the location of a trapped miner, who can deploy a wire loop receiving antenna (of unknown size or axis orientation, due to tunnel size restrictions, floor irregularities, rock falls, etc.) connected to a battery powered electronic system capable of receiving and measuring the six transmitted signals, and relaying the information to the surface. In this embodiment, the six necessary transmissions could be made in time sequence, and the magnitude and crude phase of each measurement relayed. Alternate embodiments include transmission of multiple frequencies, with modulation to provide a phase reference; transmission of independent pseudorandom-coded modulation; and computation of position before relaying to minimize the amount of relayed data.

Referring now to FIG. 11, there is shown a class of cases employing a pair of distinctively encoded (e.g. time shared, frequency separated or otherwise encoded) single-axis transmitters and one three-axis receiver. In this class of cases, the relative position of the two transmitters T.sub.1 and T.sub.2 is known to define a reference set of coordinates. The location of the transmitter stations may be known to derive the unknown location of the three axes receiver or, in the alternative, the position of the three axes receiver may be known to determine the location of the single axis transmitters, T.sub.1 and T.sub.2.

The case of FIG. 11 is particularly applicable to surveying deep boreholes by deploying two separated transmitting loops 16a on the surface or sea floor, capable of producing the very large moments needed for deep penetration. (The source locations, moments, and orientations will be determined by measuring their fields at multiple known locations, as in the cases of FIGS. 7-9.) This method usually results in complete specification of the problem in six unknowns, e.g., both source loops 16a and the sensor 22 location all co-planar, with unknown attenuation. Thus, the location and orientation of a three-axis receiver sensor 22 in a borehole may be determined by transmission of received or processed data from the receiver 12 to the surface either by wireline or by telemetry. In this case, as in all cases in which a multi-axis sensor (or transmitter) is located in a borehole, penetration of casing will invalidate the computations unless accurate values of axial and radial attenuation are known, since these values will generally be different.

Alternatively, if the location of the receiver sensor 22 is known and the relative positions and orientations of the two transmitters T.sub.1 and T.sub.2 one with respect to the other is known, the location of the transmitters T.sub.1 and T.sub.2 can be calculated from the field measurements derived from the three axes receiver sensor 22.

As in the cases of FIGS. 7-9, if the geometry happens to contain a locus of ambiguity, additional information such as the orientation of the unknown (u) and the constraints can be obtained in order to achieve a solution for location.

Referring now to FIG. 12, there is shown a case employing a single, distinctively encoded, two-axis transmitter and a single three-axis receiver for determining relative location of the transmitter station relative to the receiver station. In this embodiment, the orientation of the transmitter is known and the moment of the transmitter is known. If the location of the two axis transmitter is known, the unknown position of the three-axis receiver (R) can be calculated from the field measurements derived from the receiver. Alternatively, if the receiver location is known, the location of the two-axis transmitter can be calculated from field signals derived from the receiver. In this embodiment, the two axes of the transmitter are preferably nearly perpendicular and the transmitting antennas 16 are co-located. This permits use of a single transmitting antenna structure, although it requires at least two large dimensions (e.g., crossed iron-cored solenoids), precluding its use in a borehole if substantial magnetic moment is required. It should also be noted that penetration of casing of unknown attenuation, or any other loss of calibration of received signals relative to the transmitter moment, will invalidate the position calculation in this case.

Referring now to FIGS. 13 and 14, there is shown a class of cases employing only one single-axis transmitter and one three-axis receiver.

In FIG. 13, a case is shown wherein the single-axis transmitter T is moved to two or more unknown locations, e.g., lowered in known steps down a borehole and a single three-axis receiver sensor 22 is positioned at a known location.

The case of FIG. 13 is useful for locating cased vertical shafts relative to drifts. An iron-cored solenoid transmitter antenna 16, driven by wireline from the surface, is lowered in measured steps down the shaft, and three axes field measurements are made by a fixed receiver sensor 22 in the drift for each transmitter location, T, T', etc. The casing attenuation and, therefore, the effective transmitter moment, are unknown. Thus, even assuming the transmitter antenna orientation to be known (vertical), no single location of the transmitter is completely specified, since there are four unknowns (3 in position; 1 in moment) and only three measurements. To provide the additional information needed to specify a solution, it is assumed that the transmitter moment is constant, vertical, and moved in known steps on a vertical axis. Theoretically, two such locations are sufficient. In practice, there may be anomalies (of unknown location) in the casing, such as holes, pipes, and pumps. Therefore, data are collected at many steps, and those points that obviously violated the assumptions are discarded. A "best fit" (overspecified) solution is then computed from the remainder.

In alternate embodiments of this method, utilizing more general computation approaches, constraints will be relaxed. For instance, an arbitrary but invariant transmitter moment and direction, moved in known steps along an axis of arbitrary direction, is solvable by the iterative "best fit" computation process above described relative to FIG. 13.

Referring now to FIG. 14a, there is shown a special case for the class of FIGS. 13 and 14 wherein a single axis transmitter T of known location, T(k) orientation and moment is employed with a three-axis receiver sensor 22 of know orientation. This method is particularly suitable for deep hole surveying, similar to the embodiment of FIG. 13, but when the receiver sensor orientation is known, such as by incorporation of tiltmeters and a magnetic compass at the receiver location (which may use the same magnetic sensor as the location system) in the downhole tool. In this case, only a single transmitting antenna 16 is required allowing the dedication of all available real estate to one large loop 16a for extremely deep holes, or restricted access on the surface. The position of the receiver is calculated from the magnetic field components received at the receiver station and the sensor orientation information.

In FIG. 14b, a single-axis transmitter (T) of unknown (u) location is employed in a drilling assembly and a single three-axis receiver sensor 22 is positioned within a core (eg 45.degree. half angle) ahead of the drill at a known (k) location (e.g., on the surface or in a tunnel or borehole) and generally in the plane of the hole being drilled. With a calibrated transmitter moment and continuous measure of its inclination (by wireline on modulation of the navigation signal) or length of the drill string in the hole, a predetermined vertical profile can be followed by calculating the depth and steering accordingly. An alternative is to use the inclination and length of drill string only for vertical steering. Horizontal steering can be achieved by using a magnetic compass or the horizontal magnetic field components received at the sensor station. If the sensor is at the final target location, the received fields can be used to steer the drill straight toward an aligned sensor or in an arc toward an inclined sensor (mechanically or electrically rotated) without knowledge of moment, inclination or pipe length. The locus of the final borehole can be calculated by the iterative "best fit" process using data collected during drilling.

An alternative is shown in FIG. 14c in which the single three-axis receiver sensor 22 is at the unknown (u) location in the drill assembly and a single-axis transmitter is positioned at the known (k) location. This configuration has advantages when it is desired to use the sensor as a magnetic compass and a receiver but vibration-induced noise will usually preclude its use while actually drilling.

With regard to the above examples of FIGS. 7-14, more than one transmitter frequency may be employed to resolve propagation attenuation problems. When employing more than one transmitting antenna or transmitting axis, the transmitters may be switched in an unambiguous sequence (e.g. off, A, B) to provide reference information to the receiver sensor 22.

With regard to the receiver stations, the receiver sensor 22 may be adapted for sensing the static magnetic (earth) field to define receiver sensor orientation. This information may be utilized in the calculations deriving the position of the unknown transmitter or receiver stations. Also with regard to the receiver station, the receiver sensors 22 may be mounted on a common beam to control relative alignment and spacing therebetween. ##SPC1## ##SPC2## ##SPC3##