Method of correcting axial and transverse error components in magnetometer readings during wellbore survey operations

Abstract

The present invention is directed to a method for compensating V2 for bias error on one or more magnetometers used for one or more borehole directional surveys, comprising determining the constant or first-order varying bias values which minimize the variance of computed local magnetic field vectors with respect to central values, which may be either independent estimates of the local magnetic field vector or the mean values of the computed local magnetic field vectors.

Description

BACKGROUND OF THE INVENTION

1. Field of the Invention:

Surveying of wellbore orientation is commonly performed by the use of instruments containing sets of three orthogonal accelerometers and magnetometers, which are inserted within the drillstring and used to measure the orientations of the local gravitational and magnetic field vectors. In order to measure the earth's magnetic field, which is used as a north reference from which wellbore azimuth may be computed, the instruments must be placed within a section of non-magnetic material extending between upper and lower ferromagnetic drillstring sections. These ferromagnetic portions of the drillstring tend to acquire magnetization as they are repeatedly strained in the earth's magnetic field during drilling operations. The nominally non-magnetic portion of the drillstring may also acquire some lesser magnetization as a result of imperfections. The result is that magnetometer measurements made by an instrument within a drillstring may measure not the undisturbed magnetic field, but the vector sum of the earth's field and an error field caused by drillstring magnetization. Since the tool is fixed with respect to the drillstring, the error field is fixed with respect to the tool's coordinate system and it appears as bias errors on the magnetometer measurements, which can lead to errors in the determination of wellbore azimuth and trajectory unless measures are taken to compensate for these bias errors.

2. Description of the Prior Art:

Since the greater part of the drillstring magnetization occurs in the ferromagnetic portions of the drillstring, which are displaced axially from the instrument, the bias error in the axial direction usually exceeds the transverse bias errors. Various methods have therefore been published which seek to determine axial magnetometer bias errors in a single directional survey, including U.S. Pat. Nos. 3,791,043 to Russell, 4,163,324 to Russell, U.S. Pat. No. Re. 33,708 to Roesler, 4,761,889 to Cobern, 4,819,336 to Russell, 4,999,920 to Russell, and 5,155,916 to Engebretson. All of these methods require the provision of an independent estimate of one or more components of the earth's magnetic field, and as a result all of them tend to lose accuracy in those attitudes in which the direction of the independent estimate is perpendicular to the drillstring and therefore contributes little or no axial information. In particular, all of these methods lose accuracy as the wellbore attitude approaches horizontal east-west. A number of methods have also been published which seek to determine magnetometer biases on all three axes, including U.S. Pat. Nos. 4,682,421 to van Dongen and 4,956,921 to Coles, and UK Pat. No. 2,256,492 to Nicolle. While certain of these methods can resolve transverse bias components without external estimates of the field, they all require an independent estimate of the earth's magnetic field in order to determine the axial bias component, and therefore they also tend to lose accuracy as the attitude approaches horizontal east-west. U.S. Pat. No. 4,709,486 to Walters discloses a method for determining axial bias errors without any external estimate, by the simultaneous use of transverse magnetometer data from a plurality of surveys. However the method fails to make use of the valuable information contained in the axial magnetometer measurements, since it does not require any correlation between the axial biases determined at the two locations. In U.S. Pat. No. 5,321,893, Engebretson discloses a method which may be used to determine magnetometer scale factor and bias errors from a plurality of surveys with or without requiring any external estimate of the earth's field. However, the method is inherently approximate since it requires the construction of a "measurement matrix", whose elements depend on the unknown borehole attitude and magnetic dip angle.

Additional objectives, features and advantages of the present invention will be apparent in the written description which follows.

SUMMARY OF THE INVENTION

It is therefore an objective of the present invention to provide a method for determining magnetometer biases during wellbore survey operations, which is capable of determining biases on up to three axes, with or without the use of an external estimate of the local magnetic field, and which is capable of providing an accurate result using data from a minimum number of surveys. Since drillstring magnetization may be acquired or lost during the course of the drilling operation, it is a further object of the present invention to provide a method for determining magnetometer biases which may vary between surveys in a predefined manner.

BRIEF DESCRIPTION OF THE DRAWINGS

The novel features believed characteristic of the invention are set forth in the appended claims. The invention itself, however, as well as a preferred mode of use, further objectives and advantages thereof, will best be understood by reference to the following detailed description of an illustrative embodiment when read in conjunction with the accompanying drawings, wherein:

FIG. 1 shows a typical drilling operation comprising a drilling rig, a drillstring including a survey instrument, and a fluid circulating system;

FIG. 2 shows a flowchart illustrating the typical acquisition of survey data at a number of survey stations;

FIG. 3 shows a typical tool-fixed coordinate system used by a magnetic survey instrument located within a drillstring;

FIG. 4 shows a typical relationship between magnetic and gravitational vectors in the Northern hemisphere;

FIG. 5 shows typical measured values of the local magnetic vector, with drillstring magnetization;

FIG. 6 shows typical corrected values of the local magnetic vector;

FIG. 7 shows a flowchart illustrating computation of constant bias errors on three axes by minimizing variance of the computed magnetic field vector about a mean value; and

FIG. 8 shows a bias error which is assumed to vary as a function of drillstring revolutions and approximate axial magnetic component.

DETAILED DESCRIPTION OF THE INVENTION

FIG. 1 illustrates a rig engaged in drilling operations; the equipment includes a derrick 1, drawworks 2, cable 3, crown block 4, traveling block 5, and hook 6, supporting a drillstring which includes a swivel joint 7, kelly 8, drillpipe 9, drill collars 10, and drill bit 11. Pumps 12 circulate drilling fluid through a standpipe 13 and flexible hose 14, down through the hollow drillstring and back to the surface through the annular space 15 between the drillstring and the borehole wall 16. During the course of drilling a borehole for oil or gas exploration, it is advantageous to measure from time to time the orientation of the borehole in order to determine its trajectory. This can be accomplished by the use of a survey tool 17 located within the drill collars 10, for measuring the direction and magnitude of the local gravitational and magnetic fields with respect to a tool-fixed coordinate system. It is customary to take one survey each time the drilling operation is interrupted to add a new section to the drillstring; however, surveys can be taken at any time as indicated in the flowchart of FIG. 2. Referring again to FIG. 1, the measured data are transmitted to the surface by modulating a valve (not shown) placed in the flow passage within or adjacent to survey tool 17, causing pressure pulses to propagate in the mud column up the drillstring, where they are detected by a pressure transducer 18 placed in the standpipe 13 and communicated to data processing system 24 which may be located on the rig floor or in a logging trailer or other work area, which is approximately programmed to (1) to interpret the pressure pulses (2) eliminate the influence of magnetic field bias error components and (3) calculate one or more conventional wellbore orientation indicators. Data processing system 24 may be programmed in accordance with the present invention.

The borehole inclination can be determined by use of the gravitational measurements alone, while the borehole azimuth is determined from the gravitational and magnetic measurements; since the azimuth uses the direction of the local magnetic field as a north reference, it is necessary for the survey tool 17 to be placed in non-magnetic portions 19 and 20 of the drillstring situated between upper and lower ferromagnetic sections 21 and 22. Magnetization of the upper and lower ferromagnetic sections 21 and 22, as well as imperfections in the non-magnetic materials comprising the survey tool 17 and the non-magnetic collars 19 and 20 can produce a magnetic error field, which is fixed in the tool's frame of reference and which therefore appears as bias errors affecting the magnetic measurements. The present invention is directed to determining these errors in order to compensate for their presence and thus to provide more accurate measurements of borehole azimuth.

The invention will first be described as it pertains to solving for constant bias errors along each axis. It is conventional to define the tool-fixed coordinates as x,y and z, the z-coordinate being aligned with the drillstring axis as illustrated in FIG. 3. The instrument measures three components Gx, Gy and Gz of the gravitational vector G, and three components Bx, By and Bz of the magnetic field vector B. The influence of drillstring magnetization may therefore be described as bias errors .epsilon..sub.x, .epsilon..sub.y, and .epsilon..sub.z affecting the measured components Bx.sub.m, By.sub.m and Bz.sub.m of the earth's field, thus:

Bx.sub.m =Bx+.epsilon..sub.x (1) By.sub.m =By+.epsilon..sub.y (2) Bz.sub.m =Bz+.epsilon..sub.z (3)

The procedure of the present invention requires initial values for .epsilon..sub.x, .epsilon..sub.y, and .epsilon..sub.z. Initial estimates for the small transverse biases .epsilon..sub.x and .epsilon..sub.y will normally be zero or the values obtained from a previous computation, while the estimate for the axial bias .epsilon..sub.z may either be zero, the value obtained from a previous computation, or a value obtained by an alternative method such as that of Cobern U.S. Pat. No. 4,761,889, incorporated herein by reference, or other alternative prior art methods, all of which can be considered to be "independent" of the immediate calculations. Further computations are performed in an iterative loop using corrected magnetic components Bx.sub.c, By.sub.c, and Bz.sub.c, defined by:

Bx.sub.c =Bx.sub.m -.epsilon..sub.x (4) By.sub.c =By.sub.m -.epsilon..sub.y (5) Bz.sub.c =Bz.sub.m -.epsilon..sub.z (6) in which .epsilon..sub.x, .epsilon..sub.y, and .epsilon..sub.z represent the current estimates of the bias errors. It is apparent from FIG. 4 that the instrument can be used to determine the vertical and horizontal components of the magnetic field vector, using the vector dot product relationship: Bv.sub.c =(Bx.sub.c *Gx+By.sub.c *Gy+Bz.sub.c *Gz)/.vertline.G.vertline.(7) and the vector sum relationship: Bh.sub.c =(Bx.sub.c.sup.2 +By.sub.c.sup.2 +Bz.sub.c.sup.2 -Bv.sub.c.sup.2).sup.0.5 (8)

Values for Bv.sub.c and Bh.sub.c as obtained from a number of surveys N will typically exhibit some scatter as shown in FIG. 5, reflecting the varying direction of the residual uncompensated error components. The invention seeks to minimize this scatter by minimizing the variance .SIGMA.L.sup.2 of the computed components Bv.sub.c and Bh.sub.c with respect to central values Bv.sub.0 and Bh.sub.0. The objective is therefore to minimize the variance function:

.SIGMA.L.sup.2 =.SIGMA.[(Bh.sub.c -Bh.sub.0).sup.2 +(Bv.sub.c -Bv.sub.0).sup.2 ] (9) where .SIGMA. indicates summation over all N surveys.

This can be accomplished by differentiating Equation 9 with respect to the unknowns .epsilon..sub.x, .epsilon..sub.y, and .epsilon..sub.z, and setting the results to zero:

F(.epsilon..sub.x, .epsilon..sub.y, .epsilon..sub.z)=.SIGMA.[(Bh.sub.c -Bh.sub.0)*(.differential.Bh.sub.c /.differential..epsilon..sub.x -.differential.Bh.sub.0 /.differential..epsilon..sub.x)+(Bv.sub.c -Bv.sub.0)*(.differential.Bv.sub.c /.differential..epsilon..sub.x -.differential.Bv.sub.0 /.differential..epsilon..sub.x)]=0(10) G(.epsilon..sub.x, .epsilon..sub.y, .epsilon..sub.z)=.SIGMA.[(Bh.sub.c -Bh.sub.0)*(.differential.Bh.sub.c /.differential..epsilon..sub.y -.differential.Bh.sub.0 /.differential..epsilon..sub.y)+(Bv.sub.c -Bv.sub.0)*(.differential.Bv.sub.c /.differential..epsilon..sub.y -.differential.Bv.sub.0 /.differential..epsilon..sub.y)]=0(11) H(.epsilon..sub.x, .epsilon..sub.y, .epsilon..sub.z)=.SIGMA.[(Bh.sub.c -Bh.sub.0)*(.differential.Bh.sub.c /.differential..epsilon..sub.z -.differential.Bh.sub.0 /.differential..epsilon..sub.z)+(Bv.sub.c -Bv.sub.0)*(.differential.Bv.sub.c /.differential..epsilon..sub.z -.differential.Bv.sub.0 /.differential..epsilon..sub.z)]=0(12)

The unknown error components .epsilon..sub.x, .epsilon..sub.y, and .epsilon..sub.z are small perturbations, and an approximate solution can therefore be found by linearizing the Equations 10 through 12. The linearized form of these equations is:

a.sub.1 *.epsilon..sub.x +b.sub.1 *.epsilon..sub.y +c.sub.1 *.epsilon..sub.z =g.sub.1 (13) a.sub.2 *.epsilon..sub.x +b.sub.2 *.epsilon..sub.y c.sub.2 *.epsilon..sub.z =g.sub.2 (14) a.sub.3 *.epsilon..sub.x +b.sub.3 *.epsilon..sub.y +c.sub.3 *.epsilon..sub.z =g.sub.3 (15) where: a.sub.1 =.differential./.differential..epsilon..sub.x [F(.epsilon..sub.x ', .epsilon..sub.y ', .epsilon..sub.z ')] (16) b.sub.1 =.differential./.differential..epsilon..sub.y ]F(.epsilon..sub.x ', .epsilon..sub.y ')].epsilon..sub.z ')] (17) c.sub.1 =.differential./.differential..epsilon..sub.z [F(.epsilon..sub.x ', .epsilon..sub.y ', .epsilon..sub.z ')] (18) a.sub.2 =.differential./ .differential..epsilon..sub.x [G(.epsilon..sub.x ', .epsilon..sub.y ', .epsilon..sub.z ')] (19) b.sub.2 =.differential./.differential..epsilon..sub.y [G(.epsilon..sub.x ', .epsilon..sub.y ', .epsilon..sub.z ')] (20) c.sub.2 =.differential./.differential..epsilon..sub.z [G(.epsilon..sub.x ', .epsilon..sub.y ', .epsilon..sub.z ')] (21) a.sub.3 =.differential./.differential..epsilon..sub.x [H(.epsilon..sub.x ', .epsilon..sub.y ', .epsilon..sub.z ')] (22) b.sub.3 =.differential./.differential..epsilon..sub.y [H(.epsilon..sub.x ', .epsilon..sub.y ', .epsilon..sub.z ')] (23) c.sub.3 =.differential./.differential..epsilon..sub.z [H(.epsilon..sub.x ', .epsilon..sub.y ', .epsilon..sub.z ')] (24) and: g.sub.1 =a.sub.1 *.epsilon..sub.x '+b.sub.1 *.epsilon..sub.y 'c.sub.1 *.epsilon..sub.z 'F(.epsilon..sub.x ', .epsilon..sub.y ', .epsilon..sub.z ') (25) g.sub.2 =a.sub.2 *.epsilon..sub.x '+b.sub.2 *.epsilon..sub.y '+c.sub.2 *.epsilon..sub.z '-G(.epsilon..sub.x ', .epsilon..sub.y ', .epsilon..sub.z ') (26) g.sub.3 =a.sub.3 *.epsilon..sub.x '+b.sub.3 *.epsilon..sub.y '+c.sub.3 *.epsilon..sub.z '-H(.epsilon..sub.x ', .epsilon..sub.y ', .epsilon..sub.z ') (27) and the primed error terms .epsilon..sub.x ', .epsilon..sub.y ', and .epsilon..sub.z ' represent the previous estimates of these values. The linearized Equations 13 through 15 can be solved by standard methods for the unknowns .epsilon..sub.x, .epsilon..sub.y, and .epsilon..sub.z, and these solutions can then be refined by further iteration in the loop described by Equations 4 through 15, if necessary. A suitable stopping criterion (or "convergence" criterion) may be used to determine whether further iterations are needed, such as the change in computed values of .epsilon..sub.x, .epsilon..sub.y, and .epsilon..sub.z falling below a predefined minimum limit, or a specified number of iterations having been performed. Upon completion of the computation, the corrected magnetic vectors are more closely grouped, as illustrated in FIG. 6.

The central values Bv.sub.0 and Bh.sub.0 may be taken to be either independent constants obtained from an external source, or the mean values .SIGMA.Bv.sub.c /N and .SIGMA.Bh.sub.c /N, with the latter values preferred in cases where at least some of the surveys are in attitudes approaching horizontal east-west. The various derivatives can be evaluated by use of the chain rule and differentiation of Equations 10 through 12. FIG. 7 illustrates the solution procedure in flowchart form, with the coefficients in equations 16 through 27 expanded for the case where Bv.sub.0 and Bh.sub.0 are the mean values of the computed vectors.

After determining the biases .epsilon..sub.x, .epsilon..sub.y, and .epsilon..sub.z, the compensated magnetic field components may be found by use of Equations 7 and 8, and compensated values of the borehole azimuth may be found by the use of well-known equations.

The method of obtaining coefficients a.sub.1 through c.sub.3 will be illustrated through the following examples in which coefficient a.sub.1 is computed, as defined by Equation 16. First, in the case where Bv.sub.0 and Bh.sub.0 are independent constants, Equation 10 can be simplified to:

F(.epsilon..sub.x, .epsilon..sub.y, .epsilon..sub.z)=.SIGMA.[(Bh.sub.c -Bh.sub.0)*.differential.Bh.sub.c /.differential..epsilon..sub.x +(Bv.sub.c -Bv.sub.0)*.differential.Bv.sub.c /.differential..epsilon..sub.x ] (28)

Taking the derivative with respect to .epsilon..sub.x,

a.sub.1 =.SIGMA.[(.differential.Bh.sub.c /.differential..epsilon..sub.x).sup.2 +(Bh.sub.c -Bh.sub.0)*.differential..sup.2 Bh.sub.c /.differential..epsilon..sub.x.sup.2 +(.differential.Bv.sub.c /.differential..epsilon..sub.x).sup.2 +(Bv.sub.c -Bv.sub.0)*.differential..sup.2 Bv.sub.c /.differential..epsilon..sub.x.sup.2 ] (29)

This expression can be simplified by use of the following relationships, which follow from Equations 4, 7 and 8:

.differential.Bx.sub.c /.differential..epsilon..sub.x =-1 (30) .differential.Bv.sub.c /.differential..epsilon..sub.x =-.differential.Bv.sub.c /.differential.Bx.sub.c =-G.sub.x /.vertline.G.vertline. (31) .differential.Bh.sub.c /.differential..epsilon..sub.x =-(Bx.sub.c +Bv.sub.c *Bv.sub.c /.differential..epsilon..sub.x)/Bh.sub.c (32) .differential..sup.2 Bv.sub.c /.differential..epsilon..sub.x.sup.2 =0(33) .differential..sup.2 Bh.sub.c /.differential..epsilon..sub.x.sup.2 =[1-(.differential.Bv.sub.c /.differential..epsilon..sub.x).sup.2 -(.differential.Bh.sub.c /.differential..epsilon..sub.x).sup.2 ]/Bh.sub.c( 34)

Substituting these expressions in Equation 29, yields:

a.sub.1 =.SIGMA.(1-Bh.sub.0 *.differential..sup.2 Bh.sub.c /.differential..epsilon..sub.x.sup.2)=N-Bh.sub.0 *.SIGMA.(.differential..sup.2 Bh.sub.c /.differential..epsilon..sub.x.sup.2) (35)

The coefficients become more complicated in the case where the central values Bv.sub.0 and Bh.sub.0 are not constants, but are instead defined as the means of the Bv.sub.c and Bh.sub.c values. In this case,

F(.epsilon..sub.x, .epsilon..sub.y, .epsilon..sub.z)=.SIGMA.[(Bh.sub.c -.SIGMA.Bh.sub.c /N)*{.differential.Bh.sub.c /.differential..epsilon..sub.x -.SIGMA.(.differential.Bh.sub.c /.differential..epsilon..sub.x)/N}]+.SIGMA.[(Bv.sub.c -.SIGMA.Bv.sub.c /N)*{.differential.Bv.sub.c /.differential..epsilon..sub.x -.SIGMA.(.differential.Bv.sub.c /.differential..epsilon..sub.x)/N}](36) which can be simplified to: F(.epsilon..sub.x, .epsilon..sub.y, .epsilon..sub.z)=.SIGMA.(Bh.sub.c *.differential.Bh.sub.c /.differential..epsilon..sub.x)-.SIGMA.Bh.sub.c *.SIGMA.(.differential.Bh.sub.c /.differential..epsilon..sub.x)/N+.SIGMA.(Bv.sub.c *.differential.Bv.sub.c /.differential..epsilon..sub.x)-.SIGMA.Bv.sub.c *.SIGMA.(.differential.Bv.sub.c /.differential..epsilon..sub.x)/N(37) Differentiating with respect to .epsilon..sub.x, a.sub.1 =.SIGMA.(.differential.Bh.sub.c /.differential..epsilon..sub.x).sup.2 +.SIGMA.(Bh.sub.c *.differential..sup.2 Bh.sub.c /.differential..epsilon..sub.x.sup.2)-[.SIGMA.(.differential.Bh.sub.c /.differential..epsilon..sub.x)].sup.2 /N-.SIGMA.Bh.sub.c *.SIGMA.(.differential..sup.2 Bh.sub.c /.differential..epsilon..sub.x.sup.2)/N+.SIGMA.(.differential.Bv.sub.c /.differential..epsilon..sub.x).sup.2 +.SIGMA.(Bv.sub.c *.differential..sup.2 Bv.sub.c /.differential..epsilon..sub.x.sup.2)-[.SIGMA.(.differential.Bv.sub.c /.differential..epsilon..sub.x)].sup.2 /N-.SIGMA.Bv.sub.c *.SIGMA.(.differential..sup.2 Bv.sub.c /.differential..epsilon..sub.x.sup.2)/N (38) and simplifying using Equations 30 through 34: a.sub.1 =N-[.SIGMA.(.differential.Bh.sub.c /.differential..epsilon..sub.x)].sup.2 /N-[.SIGMA.(.differential.Bv.sub.c /.differential..epsilon..sub.x)].sup.2 /N-.SIGMA.Bh.sub.c *.SIGMA.(.differential..sup.2 Bh.sub.c /.differential..epsilon..sub.x.sup.2)/N (39)

Similar calculations may be used to find the remaining coefficients. The results are given in FIG. 7, wherein each coefficient has been multiplied by N to simplify further.

It is commonly believed that drillstring magnetization results from repeated mechanical strains applied to ferromagnetic portions of the drillstring in the presence of the earth's magnetic field, during the drilling operation. It is therefore possible for the drillstring magnetization to change slowly during the course of the drilling operation. Such a change might reasonably be approximated by a linear change over time, or by an amount proportional to strain history as measured by accumulated

drillstring revolutions, or by a function including the approximate or uncompensated magnetic field along the magnetometer axis, or by functions of these and other quantities which can be determined in advance of the computation. The present invention therefore includes a further method for computing magnetic bias errors which are first-order functions of known or measured parameters. For example, FIG. 8 illustrates a case wherein the axial magnetometer bias is assumed to vary as a linear function of the accumulated product of drillstring rotations and the axial component of the magnetic vector.

To solve for bias errors which are assumed to vary as linear functions of known parameters P.sub.x, P.sub.y, and P.sub.z, the bias errors are represented as:

.epsilon..sub.x =m.sub..epsilon.x *P.sub.x +.epsilon..sub.x0(40) .epsilon..sub.y =m.sub..epsilon.y *P.sub.y +.epsilon..sub.y0(41) .epsilon..sub.z =m.sub..epsilon.z *P.sub.z +.epsilon..sub.z0(42) in which the six unknowns are the slopes m.sub..epsilon.x, m.sub..epsilon.y, and m.sub..epsilon.z, and the intercepts .epsilon..sub.x0, .epsilon..sub.y0, and .epsilon..sub.z0. Differentiation must now be performed with respect to these unknowns, thus the iterative portion of the computation must now solve the following equations: F(.epsilon..sub.x, .epsilon..sub.y, .epsilon..sub.z)=.SIGMA.[(Bh.sub.c -Bh.sub.0)*(.differential.Bh.sub.c /.differential..epsilon..sub.x -.differential.Bh.sub.0 /.differential..epsilon..sub.x)+(Bv.sub.c -Bv.sub.0)*(.differential.Bv.sub.c /.differential..epsilon..sub.x -.differential.Bv.sub.0 /.differential..epsilon..sub.x)]=0(10) G(.epsilon..sub.x, .epsilon..sub.y, .epsilon..sub.z)=.SIGMA.[(Bh.sub.c -Bh.sub.0)*(.differential.Bh.sub.c /.differential..epsilon..sub.y -.differential.Bh.sub.0 /.differential..epsilon..sub.y)+(Bv.sub.c -Bv.sub.0)*(.differential.Bv.sub.c /.differential..epsilon..sub.y -.differential.Bv.sub.0 /.differential..epsilon..sub.y)]=0(11) H(.epsilon..sub.x, .epsilon..sub.y, .epsilon..sub.z)=.SIGMA.[(Bh.sub.c -Bh.sub.0)*(.differential.Bh.sub.c /.differential..epsilon..sub.z -.differential.Bh.sub.0 /.differential..epsilon..sub.z)+(Bv.sub.c -Bv.sub.0)*(.differential.Bv.sub.c /.differential..epsilon..sub.z -.differential.Bv.sub.0 /.differential..epsilon..sub.z)]=0(12) I(.epsilon..sub.x, .epsilon..sub.y, .epsilon..sub.z)=.SIGMA.[(Bh.sub.c -Bh.sub.0)*P.sub.x *(.differential.Bh.sub.c /.differential..epsilon..sub.x -.differential.Bh.sub.0 /.differential..epsilon..sub.x)+(Bv.sub.c -Bv.sub.0)*P.sub.x *(.differential.Bv.sub.c /.differential..epsilon..sub.x -.differential.Bv.sub.0 /.differential..epsilon..sub.x)]=0(43) J(.epsilon..sub.x, .epsilon..sub.y, .epsilon..sub.z)=.SIGMA.[(Bh.sub.c -Bh.sub.0)*P.sub.y *(.differential.Bh.sub.c /.differential..epsilon..sub.y -.differential.Bh.sub.0 /.differential..epsilon..sub.y)+(Bv.sub.c -Bv.sub.0)*P.sub.y *(.differential.Bv.sub.c /.differential..epsilon..sub.y -.differential.Bv.sub.0 /.differential..epsilon..sub.y)]=0(44) K(.epsilon..sub.x, .epsilon..sub.y, .epsilon..sub.z)=.SIGMA.[(Bh.sub.c -Bh.sub.0)*P.sub.z *(.differential.Bh.sub.c /.differential..epsilon..sub.z -.differential.Bh.sub.0 /.differential..epsilon..sub.z)+(Bv.sub.c -Bv.sub.0)*P.sub.z *(.differential.Bv.sub.c /.differential..epsilon..sub.z -.differential.Bv.sub.0 /.differential..epsilon..sub.z)]=0(45) where Equations 10 through 12 are obtained by differentiation of Equation 9 with respect to the unknown intercepts .epsilon..sub.x0, .epsilon..sub.y0, and .epsilon..sub.z0, and Equations 43 through 45 are obtained by differentiation with respect to the unknown slopes m.sub..epsilon.x, m.sub..epsilon.y, and m.sub..epsilon.z. Improved estimates of the six unknowns can be found by solving the six linearized equations: a.sub.1 *e.sub.x0 +b.sub.1 *e.sub.y0 +c.sub.1 *e.sub.z0 +d.sub.1 *m.sub..epsilon.x +e.sub.1 *m.sub..epsilon.y +f.sub.1 *m.sub..epsilon.z =g.sub.1 (46) a.sub.2 *e.sub.x0 +b.sub.2 *e.sub.y0 +c.sub.2 *e.sub.z0 +d.sub.2 *m.sub..epsilon.x +e.sub.2 *m.sub..epsilon.y +f.sub.2 *m.sub..epsilon.z =g.sub.2 (47) a.sub.3 *e.sub.x0 +b.sub.3 *e.sub.y0 +c.sub.3 *e.sub.z0 +d.sub.3 *m.sub..epsilon.x +e.sub.3 *m.sub..epsilon.y +f.sub.3 *m.sub..epsilon.z =g.sub.3 (48) a.sub.4 *e.sub.x0 +b.sub.4 *e.sub.y0 +c.sub.4 *e.sub.z0 +d.sub.4 *m.sub..epsilon.x +e.sub.4 *m.sub..epsilon.y +f.sub.4 *m.sub..epsilon.z =g.sub.4 (49) a.sub.5 *e.sub.x0 +b.sub.5 *e.sub.y0 +c.sub.5 *e.sub.z0 +d.sub.5 *m.sub..epsilon.x +e.sub.5 *m.sub..epsilon.y +f.sub.5 *m.sub..epsilon.z =g.sub.5 (50) a.sub.6 *e.sub.x0 +b.sub.6 *e.sub.y0 +c.sub.6 *e.sub.z0 +d.sub.6 *m.sub..epsilon.x +e.sub.6 *m.sub..epsilon.y +f.sub.6 *m.sub..epsilon.z =g.sub.6 (51) in which the coefficients are defined as follows: a.sub.1 =.differential./.differential..epsilon..sub.x0 [F(.epsilon..sub.x.sup.', .epsilon..sub.y.sup.', .epsilon..sub.z.sup.')](52) b.sub.1 =.differential./.differential..epsilon..sub.y0 [F(.epsilon..sub.x.sup.', .epsilon..sub.y.sup.', .epsilon..sub.z.sup.')](53) c.sub.1 =.differential./.differential..epsilon..sub.z0 [F(.epsilon..sub.x.sup.', .epsilon..sub.y.sup.', .epsilon..sub.z.sup.')](54) d.sub.1 =.differential./.differential.m.sub..epsilon.x [F(.epsilon..sub.x.sup.', .epsilon..sub.y.sup.', .epsilon..sub.z.sup.')](55) e.sub.1 =.differential./.differential.m.sub..epsilon.y [F(.epsilon..sub.x.sup.', .epsilon..sub.y.sup.', .epsilon..sub.z.sup.')](56) f.sub.1 =.differential./.differential.m.sub..epsilon.z [F(.epsilon..sub.x.sup.', .epsilon..sub.y.sup.', .epsilon..sub.z.sup.')](57) and coefficients a.sub.2 through f.sub.6 are obtained similarly, by differentiation of functions G(.epsilon..sub.x, .epsilon..sub.y, .epsilon..sub.z) through K(.epsilon..sub.x, .epsilon..sub.y, .epsilon..sub.z), and g.sub.1 =a.sub.1 *.epsilon..sub.x0.sup.' +b.sub.1 *.epsilon..sub.y0.sup.' +c.sub.1 *.epsilon..sub.z0.sup.' +d.sub.1 *m.sub..epsilon.x.sup.' +e.sub.1 *m.sub..epsilon.y.sup.' +f.sub.1 *m.sub..epsilon.z.sup.' -F(.epsilon..sub.x.sup.', .epsilon..sub.y.sup.', .epsilon..sub.z.sup.')(58) g.sub.2 =a.sub.2 *.epsilon..sub.x0.sup.' +b.sub.2 *.epsilon..sub.y0.sup.' +c.sub.2 *.epsilon..sub.z0.sup.' +d.sub.2 *m.sub..epsilon.x.sup.' +e.sub.2 *m.sub..epsilon.y.sup.' +f.sub.2 *m.sub..epsilon.z.sup.' -G(.epsilon..sub.x.sup.', .epsilon..sub.y.sup.', .epsilon..sub.z.sup.')(59) g.sub.3 =a.sub.3 *.epsilon..sub.x0.sup.' +b.sub.3 *.epsilon..sub.y0.sup.' +c.sub.3 *.epsilon..sub.z0.sup.' +d.sub.3 *m.sub..epsilon.x.sup.' +e.sub.3 *m.sub..epsilon.y.sup.' +f.sub.3 *m.sub..epsilon.z.sup.' -H(.epsilon..sub.x.sup.', .epsilon..sub.y.sup.', .epsilon..sub.z.sup.')(60) g.sub.4 =a.sub.4 *.epsilon..sub.x0.sup.' +b.sub.4 *.epsilon..sub.y0.sup.' +c.sub.4 *.epsilon..sub.z0.sup.' +d.sub.4 *m.sub..epsilon.x.sup.' +e.sub.4 *m.sub..epsilon.y.sup.' +f.sub.4 *m.sub..epsilon.z.sup.' -I(.epsilon..sub.x.sup.', .epsilon..sub.y.sup.', .epsilon..sub.z.sup.')(61) g.sub.5 =a.sub.5 *.epsilon..sub.x0.sup.' +b.sub.5 *.epsilon..sub.y0.sup.' +c.sub.5 *.epsilon..sub.z0.sup.' +d.sub.5 *m.sub..epsilon.x.sup.' +e.sub.5 *m.sub..epsilon.y.sup.' +f.sub.5 *m.sub..epsilon.z.sup.' -J(.epsilon..sub.x.sup.', .epsilon..sub.y.sup.', .epsilon..sub.z.sup.')(62) g.sub.6 =a.sub.6 *.epsilon..sub.x0.sup.' +b.sub.6 *.epsilon..sub.y0.sup.' +c.sub.6 *.epsilon..sub.z0.sup.' +d.sub.6 *m.sub..epsilon.x.sup.' +e.sub.6 *m.sub..epsilon.y.sup.' +f.sub.6 *m.sub..epsilon.z.sup.' -K(.epsilon..sub.x.sup.', .epsilon..sub.y.sup.', .epsilon..sub.z.sup.')(63)

After determining the various slopes and intercepts by solving the set of linearized equations, new estimates of the biases can be made at each survey station using equations 28 through 30, and these improved estimates can be used in further iteration as required.

Although the invention has been described in detail in forms wherein the unknown biases are either all constant or all varying as functions of known parameters, it will be understood by those skilled in the art that the invention may also be used to solve for biases of mixed character; for example, some biases constant and others varying, some biases zero and others varying, or some biases zero and others constant. This can be understood by reference to Equations 46 through 51, which represent the general case of six linearized equations with six unknowns. Reduction of one bias from variable to constant reduces the number of unknowns by one since the previously unknown slope is now known to be zero, and it reduces the number of equations by one because the row of coefficients obtained by differentiation with respect to the slope may be eliminated. Similarly, reduction of one bias from constant to zero will also reduce the number of unknowns and equations by one. In particular, if the three slopes are set to zero, the six Equations 46 through 51 reduce to the three Equations 13 through 15 which can be solved for constant biases.

Preferably, the techniques of the present invention are utilized in wellbore surveying operations. As is conventional, the measurements are transmitted to the earth's surface utilizing measurement-while-drilling data transmission techniques. These data may be entered into a computing device which is preprogrammed in accordance with the method of the present invention. The program will include as inputs the x-axis, y-axis and z-axis components of the local magnetic and gravitational field vectors at each survey station. The calculations are performed in accordance with the description provided hereinabove, and the program provides as an output for each survey station the wellbore's azimuth and inclination which are so useful in defining the position of a wellbore.

While the invention has been shown in only certain forms, it is not thus limited but is susceptible to various changes and modifications without departing from the spirit thereof.

Claims

1. A method of surveying in a wellbore, comprising the method steps of:

(a) providing a drillstring;

(b) providing a sensor array including magnetic field sensors and gravitational field sensors;

(c) locating said sensor array within said drillstring;

(d) locating said drillstring in said wellbore;

(e) utilizing said drillstring to extend said wellbore;

(f) at a number of survey stations within said wellbore, utilizing said sensor array to make measurements of (a) three components of gravitational field, and (b) three components of magnetic field, with each of said three components of magnetic field possibly including a magnetic field bias error component;

(g) communicating said measurements to a processor which is programmed to execute a set of program instructions, including:

(i) instruction means for selecting for said number of survey stations within said wellbore, an initial estimate for said magnetic field bias error component for particular ones of said three components of magnetic field;

(ii) instruction means for computing, for each survey station within said wellbore, corrected values for said three components of magnetic field utilizing said initial estimate for said magnetic field bias error component for particular ones of said three components of magnetic field;

(iii) instruction means for computing, for each survey station within said wellbore a corrected magnetic field vector utilizing said corrected values for said three components of magnetic field;

(iv) instruction means for computing an improved estimate for said magnetic field bias error component for particular ones of said three components of magnetic field;

(h) establishing a preselected convergence criterion for said corrected magnetic field vector;

(i) iteratively executing said set of program instructions in order to generate, for each survey station within said wellbore, a corrected magnetic field vector which satisfies said preselected convergence criterion and which minimizes the variance of said corrected magnetic field vector with respect to a predetermined value; and thereby

(j) providing substantially corrected magnetic field measurements for each survey station within said wellbore.

2. A method of surveying in a wellbore according to claim 1, further comprising:

(k) utilizing said instruction means for selecting and supplying, for each survey station within said wellbore, an initial estimate for said magnetic field bias error component for particular ones of said three components of magnetic field which is at least partially based upon a magnetic field bias error component determined from at least one preceding survey station.

3. A method of surveying in a wellbore according to claim 1, further comprising:

(k) utilizing said measurements of said three components of gravitational field and said substantially corrected magnetic field measurements to calculate at least one wellbore orientation indicator.

4. A method of surveying in a wellbore according to claim 3, wherein said at least one wellbore orientation indicator includes wellbore azimuth.

5. A method of surveying in a wellbore according to claim 1, further comprising:

(k) deriving said predetermined value from an independent estimate of earth's magnetic field.

6. A method of surveying in a wellbore according to claim 1, further comprising:

(k) deriving said predetermined value from a plurality of previous calculations of corrected magnetic field vector.

7. A method of surveying in a wellbore according to claim 1, further comprising:

(k) deriving said predetermined value from an averaging of a plurality of previous calculations of corrected magnetic field vectors.

8. A method of surveying in a wellbore according to claim 1:

wherein said step of establishing a preselected convergence criterion comprises establishing a preselected minimum variance for said corrected magnetic field vector for each survey station;

wherein said step of iteratively executing comprises:

iteratively executing said set of program instructions in order to generate, for each survey station within said wellbore, a corrected magnetic field vector which satisfies said preselected minimum variance and which minimizes the variance of said corrected magnetic field vector with respect to a central value.

9. A method of surveying in a wellbore according to claim 1:

wherein said step of establishing a preselected convergence criterion comprises establishing a preselected number of iterations for calculation of said corrected magnetic field vectors for each survey station;

wherein said step of iteratively executing comprises:

iteratively executing said set of program instructions in order to generate, for each survey station within said wellbore, a corrected magnetic field vector in a manner which satisfies said preselected number of iterations.

10. A method of surveying in a wellbore according to claim 1, further comprising:

determining at least one linear relationship between a parameter and said magnetic field bias error component for at least one of said three components of magnetic field; and

including in said set of program instructions means for developing an improved estimate for said magnetic field bias error component which utilizes said at least one linear relationship.

11. A method for determining bias errors affecting magnetometers used in borehole surveying, comprising the steps of:

a) measuring three components of the gravitational field and three components of the local magnetic field at a number of survey stations;

b) selecting initial estimates for said bias errors;

c) computing corrected magnetic field components by use of said estimated bias errors;

d) using said corrected magnetic field components to compute a corrected magnetic field vector for each of said survey stations;

e) computing improved estimates of said magnetometer bias errors which minimize variance of said magnetic vectors with respect to a central value; and

f) iteratively repeating steps c) through e) using said improved estimates, until a defined convergence criterion is met.

12. The method of claim 11 wherein said bias errors comprise a single constant bias error along the magnetometer axis parallel to the drillstring.

13. The method of claim 11 wherein said bias errors comprise constant bias errors along all three magnetometer axes.

14. The method of claim 11 wherein said bias errors comprise a single bias error which varies as a defined function of known quantities, along the magnetometer axis parallel to the drillstring.

15. The method of claim 11 wherein said bias errors comprise bias errors which vary as defined functions of known quantities, along all three magnetometer axes.

16. The method of claim 11 wherein said bias errors comprise two constant bias errors along the transverse magnetometer axes, and a single bias error which varies as a defined function of known quantities along the magnetometer axis parallel to the drillstring.

17. The method of claim 11 wherein said central value is obtained from an independent estimate of the earth's magnetic field.

18. The method of claim 11 wherein said central value is obtained from the mean values of said corrected magnetic field vectors.

19. The method of claim 11 wherein said convergence criterion consists of determining whether the change in bias error between successive iteration exceeds a predefined minimum limit.

20. The method of claim 11 wherein said convergence criterion consists of performing a predefined number of iterations.

21. The method of claim 11 wherein borehole azimuth at one or more of said survey stations is determined from said components of the gravitational field and said corrected magnetic field components.