Anti-collision method for drilling wells

Methods for drilling a new well in a field having a plurality of existing cased wells using magnetic ranging while drilling are provided. In accordance with one embodiment, a method of drilling a new well in a field having an existing cased well includes drilling the new well using a bottom hole assembly (BHA) having a drill collar having by an insulated gap, generating a current on the BHA while drilling the new well, such that some of the current passes through a surrounding formation and travels along a casing of the existing cased well, measuring from the BHA a magnetic field caused by the current traveling along the casing of the existing cased well, and adjusting a trajectory of the BHA to avoid a collision between the new well and the existing cased well based on measurements of the magnetic field.

BACKGROUND OF THE INVENTION

The present invention relates generally to well drilling operations and, more particularly, to well drilling operations using magnetic ranging while drilling to avoid collisions with existing cased wells.

With conventional drilling practices, the uncertainties in a well's position increase as the depth of the well increases. These uncertainties are usually represented as ellipsoids that are centered on the location of the well as determined by Measurement While Drilling (MWD) or wireline survey data. An ellipsoid corresponds to a certain probability density corresponding to whether the well bore is actually located within the ellipsoid. The uncertainties in the well position arise from the limited accuracy of the well bore direction, inclination, and depth measurements which may be obtained from MWD and/or wireline surveys, as documented extensively. For example, MWD inclination measurements are typically accurate to no better than 0.1°, while MWD directional measurements are typically accurate to no better than 1°. Moreover, MWD survey points may be acquired only once every 90 feet in practice. Thus, under-sampling may significantly increase the actual errors in the well position.

An additional source of survey error arises because the directional measurement is based on the magnetic field, which requires correction for variations in the Earth's magnetic field, and which can also be strongly perturbed by nearby casing. If the casings are very close to the well path, then the MWD directional measurement may not even be useful. Under such conditions, a gyro may be used to provide the directional information. The gyro may be run with the MWD tool, or it may be run on wireline with periodic descents inside the drill pipe to the bottom hole assembly (BHA). Finally, an accurate MWD depth measurement is difficult to achieve, with depth errors of 1/1000 common.

Further complications may arise in older fields with existing wells. In older fields, the survey information on existing wells may be very low quality, survey data may have been lost, or the wells may have been drilled without running a MWD or wireline survey.

Wells associated with a typical offshore platform are drilled vertically for a considerable depth before they are deviated to reach distant portions of the reservoir. These vertical sections typically range from several hundred feet to a few thousand feet before they reach the kick-off point (KOP) where directional drilling begins. Because offshore production platforms are very expensive and have as many wells as possible given the limited surface area of the platform, well heads are packed as closely as possible. The distances between well heads, and therefore the number of wells, are limited primarily by the uncertainty in well positions and the risk of accidentally drilling into a cased well. Since an existing cased well and the drill bit could be located anywhere inside the respective ellipsoids of uncertainty, well heads are spaced a distance apart so that any two ellipsoids cannot overlap.

Existing platforms may have filled many or all of the available slots (i.e., locations for well heads) based on factors derived from MWD direction and inclination technology. In order to tap additional oil or gas resources, new wells may be drilled. Unless there is a reliable method to avoid drilling into an existing well, another platform may have to be built. However, if one could thread new wells among the existing wells without risk of collision, then a new platform may not be needed.

SUMMARY

In accordance with one embodiment of the invention, a method of drilling a new well in a field having an existing cased well includes drilling the new well using a bottom hole assembly (BHA) having a drill collar having by an insulated gap, generating a current on the BHA while drilling the new well, such that some of the current passes through a surrounding formation and travels along a casing of the existing cased well, measuring from the BHA a magnetic field caused by the current traveling along the casing of the existing cased well, and adjusting a trajectory of the BHA to avoid a collision between the new well and the existing cased well based on measurements of the magnetic field. The relative position of the new well to the existing well may be estimated based on measurements of the magnetic field. An alarm may be triggered if an apparent distance between the new well and the existing cased well approaches less than a threshold distance.

DETAILED DESCRIPTION OF SPECIFIC EMBODIMENTS

FIG. 1is a schematic10illustrating the spacing of two proximate wells at an offshore platform. A first well12and a second well14have wellheads16and18, respectively, extending from a platform area20. The initial placement of the first well12and the second well14is based on a well head separation Xd, the determination of which is discussed below. Based on potential survey errors associated with drilling and the casing diameter Xc, as the first well12and second well14extend to a depth D, ellipsoids of uncertainty22increase correspondingly until reaching a kick-off point (KOP)24. Each ellipsoid of uncertainty22corresponds respectively to a certain probability density corresponding to whether the well bore is actually located within the ellipsoid. As apparent in the schematic10, the final ellipsoids of uncertainty22at the KOP24are represented as E1 and E2. Upon reaching the KOP24, the first well12and the second well14deviate for directional drilling.

Well head separation Xd for the first well12and the second well14may be based on a relationship known as oriented safety factor (OSF). To ensure no collision occurs, the final ellipsoids of uncertainty22at the depth D may not overlap. The OSF may be defined according to the following equation:

In equation (1) above, Xdrepresents the well head separation, Xcrepresents the casing diameter, and E1and E2represent the radii of the ellipsoids at the depth D. The larger the oriented safety factor, the less likely that two wells will collide. Typically, one wants OSF>1.5 for a sufficient safety factor to avoid a collision.

By way of example, suppose the first well12and the second well14are vertical for a depth D=500 m, and that the casings on both wells will be 30 inches in diameter, such that Xc=0.76 m. Also, assume that the ellipsoids of uncertainty22are solely determined by the accuracy of the measurement while drilling (MWD) inclination measurement (α=2·10−3radians, ˜0.1°, and that the accuracy is the same for any new well as for existing cased wells. Hence, at 1500 ft, E1=E2=α·D=0.9 m, and a new well must be separated from existing wells by Xd=Xc+OSF·√{square root over ((E1)2+(E2)2)}{square root over ((E1)2+(E2)2)}=0.76 m+1.5·√{square root over (2)}·(0.9 m)≈2.8 m.

Note that the slot spacing may be primarily determined by the accuracy of the MWD tool. If the MWD measurements are less accurate, or if the wells must go to greater depths, or if a greater safety margin is desired, the distance between slots may generally be increased. Using the techniques disclosed herein, however, a driller may plan and subsequently drill within the ellipsoids of uncertainty22that may be determined based on MWD tool capabilities. Thus, the slot spacing may be reduced, as discussed below.

FIG. 2illustrates a schematic view26of existing wells from an offshore platform. In the schematic view26, an offshore platform28includes a plurality of wells30. After penetrating a seabed32, the wells30remain in a largely parallel configuration34through a depth D. Upon reaching a kick-off point (KOP)36, the wells30deviate into directional wells38.

FIG. 3depicts an exemplary well slot pattern40for drilling additional wells amid the plurality of wells30ofFIG. 2. Within a platform perimeter42, each existing well44is represented by a circle and each proposed well46is represented by a star. The existing wells44have been drilled with a well head spacing Xd of 2.8 meters (m). Given the limited space within the platform perimeter42, this spacing provides a maximum number of existing wells30when the ellipsoids of uncertainty22have a 2.8 meter diameter at the depth D of the kick-off point (KOP)36where the wells30deviate.

Using a technique discussed below, the ellipsoids of uncertainty22may be reduced to 2.0 meters in diameter at the depth D. Accordingly, an additional thirty-seven proposed wells46may be drilled within the platform perimeter42amid the existing wells44, more than doubling the total number of wells30on the offshore platform28. To accommodate the new well heads, a second floor may be added to the offshore platform28, above or below the initial floor. This configuration could save the cost of building an additional offshore platform when additional wells are desired.

Turning toFIG. 4, a well placement schematic48illustrates a placement of a new well50amid four existing wells52,54,56, and58on the offshore platform28when well head spacing of 2.0 meters (m) for new wells may be achieved. For the purposes of the discussion, the new well50and the existing wells52,54,56, and58are assumed to be vertical for the first few hundred meters before diverging at different angles. The well head of the new well50is located at (x,y,z)=(0,0, zp), and the well heads of the existing wells52,54,56, and58are located at (x,y,z)=(2,0,zp), (0,2,zp), (−2,0,zp), (0,−2,zp), respectively, where the floor of the offshore platform28is at zpand the z-direction is vertical. The well head spacing Xd between the four existing cased wells is 2.8 m, consistent with the example inFIG. 3.

FIG. 5provides a schematic64of a bottom hole assembly (BHA)66for drilling amid the four existing wells52,54,56, and58ofFIG. 4. The BHA66is aligned vertically on the z-axis68, drilling downward with a drill bit70coupled to a rotary steerable system (RSS)72for setting the direction of the drill bit70. The BHA66further includes an electric current driving tool74, which may be a component of a measurement while drilling (MWD) tool or a standalone tool, such as Schlumberger's E-Pulse or E-Pulse Express tool. The electric current driving tool74provides an electric current76to an outer drill collar78of the BHA66. The outer drill collar78is separated from the rest of the BHA66by an insulated gap80in the drill collar, over which electric current may not pass.

As discussed above, the electric current driving tool74may provide the electric current76to the outer drill collar78. The current76produced by the electric current driving tool74may, for example, have a frequency between about 1 Hz and about 100 Hz, and may have an amplitude of around17amps. Beginning along the outer drill collar78of the BHA66, the current76may subsequently enter the formation surrounding the BHA66. The portion of the current76that enters the surrounding formation is depicted as an electric current82.

The casing on existing wells52,54,56, and58provides very low resistance to electricity as compared to the surrounding formation. As a result, a substantial portion of the current82will pass along the casing of the existing wells52,54,56, and58. For purposes of simplification, the current82is depicted as flowing toward the casing of the existing well52, but it should be noted that the current82will be divided among the existing wells52,54,56, and58. The portion of the current82which travels along the casing of the existing well52is illustrated as current84. The current84travels along the casing of the existing well52before re-entering the formation as a current86toward the BHA66. When the current86reaches the BHA66, the resulting current is depicted as a current88, which completes the circuit at the electric current driving tool74.

The movement of the current84along the casing of the existing well52creates an azimuthal magnetic field90centered on the casing of the existing well52. A magnetometer tool92having a three-axis magnetometer94may detect both the magnitude and the direction of the magnetic field90along three axes. The magnitude and direction of the magnetic field90may provide measurements for estimating the direction and distance from the BHA66to the existing well52according to techniques discussed below.

The BHA66may include a variety of tools and configurations. For example, the RSS72may be a PowerDrive RSS. Circulating drilling mud may power the PowerDrive RSS cartridge. Because the PowerDrive RSS has a magnetometer at 126 inches behind the bit, the magnetometer tool92may form a part of the PowerDrive RSS. Such a configuration could be used to measure the induced magnetic field90generated by the current84on the casing of the existing well52. To do so, the control cartridge of the PowerDrive RSS could be maintained in geostationary mode while it is measuring the induced magnetic field90.

Above the RSS72, the BHA66may include a SlimPulse MWD tool. Because the SlimPulse MWD tool has a magnetometer located at 254 inches from the bit, the magnetometer tool92may alternatively or additionally form a part of the SlimPulse MWD tool. The SlimPulse tool is battery powered, so it can acquire data with the mud pumps on or off. After the induced magnetic field90has been measured, the data may be transmitted to the surface by the MWD pulser.

Alternatively, another MWD tool, such as a PowerPulse tool, may replace the SlimPulse tool. It is also possible to replace the PowerDrive RSS by an Exceed RSS or simply by a mud motor with a steerable assembly. A special purpose tool including both the magnetometer tool92and the electric current driving tool74may be used in place of the SlimPulse MWD tool, and the E-Pulse tool used to send data to the surface via electromagnetic (EM) waves. Moreover, if continuous steering data and instantaneous feedback to the steerable system are desired, a wired drill pipe may be used for telemetry.

Continuing to viewFIG. 5, the generation of the magnetic field90may be further described. The electric current76generated by the electric current driving tool74may be given by I(z,t)=(z)·cos(2πt+φ), where t represents time, f represents frequency, and φ represents phase. Hereafter, the time t and frequency f dependence is suppressed in the formulas, but should be understood. The electric current76on the BHA66, I(z), decreases with distance (z) from the insulated gap80as it flows from the BHA66into the surrounding formation. For example, between the insulated gap80and the drill bit70, the current76decreases in a nearly linear manner as I/(z)≈I/(0) (1+z/L), where L is the distance from the insulated gap80to the tip of the drill bit70, and where z<0 below the insulated gap80.

As discussed above, most of the current76that enters the surrounding formation also flows onto the casing of the existing wells52,54,56, and58to return to the BHA66above the insulated gap80. In the foregoing description, the current84, which may represent a return current moving along any ithexisting well casing may be denoted as Ii. Further, L may be assumed to be larger than the inter-well spacing for simplicity in the mathematical analysis, but the technique described herein does not depend on this assumption.

Turning toFIG. 6, a schematic96depicts geometry underlying the calculation of magnetic field90at the BHA66which, in a general case, arises due to the current84on an ithwell casing98. The magnetometer94may be located in the center of the BHA66may be understood to be located at {right arrow over (r)}m=(xm,ym,zm); the ithwell casing98may be understood to be located at {right arrow over (r)}i=(xi,yi,zi), and a vector pointing from the ithwell casing98to the BHA66may be {right arrow over (S)}i={right arrow over (r)}m−{right arrow over (r)}i. For simplicity, the BHA66and the ithwell casing98may be assumed to be parallel and aligned in the z-direction. Hence, the distance from the BHA to the ithcasing may be represented by Si2=(xm−xi)2+(ym−yi)2. Because it should be understood that the quantities are evaluated at the same depth, the explicit z dependence may be neglecting in the equations that follow.

The induced magnetic field90measured at the magnetometer94due to the current Ii on the ithwell casing98may be described according to the following equation:

It should be appreciated that equation (2) represents an expression for induced magnetic field from a long line of constant current. Under the assumption that L□ Si, this is a reasonable approximation.

Further, a total induced magnetic field90at the magnetometer94may be represented by a sum of the induced magnetic fields from all nearby casings (not depicted) according to the following equations:

It should be noted that equations (3) and (4) lack a Bz component. Due to the assumption that the BHA66and the existing wells52,54,56, and58all extend in the z-direction, the induced azimuthal magnetic field90which forms on the casing of the existing wells52,54,56, and58accordingly includes components in only the x- and y-directions.

The sum of the currents on all of the casing of the existing wells52,54,56, and58must not exceed the current76on the BHA66, as represented by the relationship

I≥∑i=1n⁢Ii.
The current84on any casing of the existing wells52,54,56, and58depends on the position of the well relative to the BHA66, the resistivities of both the formation and the cement surrounding the casing of the existing wells52,54,56, and58, and on the presence of other nearby casings. The current84and resulting induced magnetic field90for each of the existing wells52,54,56, and58may be obtained from a full 3-D numerical model, but simpler approaches may yield sufficient results.

With the assumption that L□ Si, the current distributions on adjacent casings may be approximated with a simple formula describing the conductance between two long, parallel cylinders. If two parallel conductors have a diameter D and are separated by the distance Si, then the conductance per unit length between them is given by the following relationship:

Equation (5) above applies for a homogeneous formation with a conductivity σ. The current Ii on the casing of the ithwell98is therefore proportional to Giaccording to the following equation:

In equation (6), the sum considers a total of n adjacent casings. Distant casings have a small effect and can be neglected for this analysis. Also, a small fraction of the current76of the BHA66will return though the borehole and shallow formation, but this minor effect may be neglected. However, the effects may be considered in a more rigorous analysis.

It should be noted that {right arrow over (B)}(xm,ym) is not a vector magnetic field in the normal sense. Rather, it represents the induced magnetic field90at the location of the magnetometer94inside the drill collar of the BHA66when the magnetometer94is located at coordinates (xm,ym). The current76on the BHA66itself does not produce a magnetic field inside the BHA66, but it does produce a strong magnetic field outside the BHA66. This external field due to the current76on the BHA66is not included in the expression for {right arrow over (B)}(xm,ym) for the reasons stated above, but the external magnetic field would be included in any expression for the magnetic field outside of the BHA66. Also, the expression for {right arrow over (B)}(xm,ym) includes any changes in any casing current84as the BHA66changes position.

Some specific examples of {right arrow over (B)}(xm,ym) are now given. The four existing wells52,54,56, and58surrounding the BHA66may be located at (x1,y1)=(2,0), (x2,y2)=(0,2), (x3,y3)=(−2,0), and (x4,y4)=(0,−2), while the BHA66is located at (xm,ym). Unless explicitly indicated otherwise, all distances are in meters. The current76generated at the insulated gap80of the BHA66may be I(0)≈17 amp, where the insulated gap80is defined at z=0. The diameter D of the BHA66and of the casing on the existing wells52,54,56, and58may be D=0.18 m, the length L of the BHA66below the insulated gap80may be L=15 m, the drill bit70may be located at z=−15 m, and the magnetometer94may be located at zm=−9 m. With the assumption that the current76decays linearly from the BHA66, the current on the BHA66at the location of the magnetometer94is I(−9)≈=(1−9/15) amp≈7 amp. The sum of the currents on the four adjacent casings of the existing wells52,54,56, and58is thus

If the BHA66is located at (xm,ym)=(0,0), as depicted in the well placement schematic48ofFIG. 4, then all four casings of the existing wells52,54,56, and58will have the same currents and, as the distances from the BHA66to the four casings of the existing wells52,54,56, and58are identical, the induced magnetic fields from the four casings of the existing wells52,54,56, and58will cancel. Hence, the magnetic field at the magnetometer will be {right arrow over (B)}(0,0)=0. If the BHA is closer to any ithcasing98, representing one of the existing wells52,54,56, and58, then the distance Siwill decrease, the conductance Gi will increase, and the current Ii will correspondingly increase. As a result, the induced magnetic field90, or Bi(xm,ym), due to the current84on the casing of the ithwell98will increase due to the increase in the current84and the factor Si−1in equation (4). Meanwhile, the induced magnetic fields from the casings of the other existing wells52,54,56, or58will decrease.

FIGS. 7 and 8plot the induced magnetic field90amplitude Bt (xm,ym)=|{right arrow over (B)}(xm,ym) as a function of the magnetometer94position (xm,ym) over the ranges xmε[−2.6,2.6] and ymε[−2.6,2.6]. Turning first toFIG. 7, a 3-D plot100clearly indicates the locations of casings of the four existing wells52,54,56, and58. The 3-D plot100illustrates the amplitude Bt102for the magnetic field90over the ranges xmε[−2.6,2.6] and ymε[−2.6,2.6]. A numeral104indicates the y-direction and a numeral106indicates the x-direction, such that point108is located at (x,y)=(2.6,2.6), point110is located at (x,y)=(2.6,−2.6), and point112is located at (x,y)=(−2.6,−2.6). A numeral114indicates the location of the BHA66at the center of the 3-D plot100. Four spikes in amplitude Bt102denoted by numerals116,118,120, and122indicate respectively a location of the existing wells52,54,56, and58.

FIG. 8similarly represents the induced magnetic field90amplitude Btin the form of a contour plot124. The contour plot124illustrates magnetic field90amplitude Btin microTesla (μT) using distinct hatching, as indicated in the legend126. An ordinate128illustrates the y-direction and an abscissa130illustrates the x-direction, such that point132is located at (x,y)=(2.6,2.6), point134is located at (x,y)=(2.6,−2.6), point136is located at (x,y)=(−2.6,−2.6), and point138is located at (x,y)=(−2.6,2.6). The center of the contour plot124indicates a location140of the BHA66. Four spikes in amplitude Bt denoted by numerals142,144,146, and148indicate respectively a location of the existing wells52,54,56, and58.

Turning toFIG. 9, an expanded view150of the contour plot124ofFIG. 8represents the induced magnetic field90amplitude Btover the ranges xmε[−1,1] and ymε[−1,1]. The expanded view150illustrates magnetic field90amplitude Btin microTesla (μT) using distinct hatching, as indicated in the legend152. An ordinate point154illustrates the y-direction and an abscissa156illustrates the x-direction, such that158is located at (x,y)=(1,1), point160is located at (x,y)=(1,−1), point162is located at (x,y)=(−1,−1), and point164is located at (x,y)=(−1,1). The center of the contour plot166indicates a location140of the BHA66. Though the four spikes in amplitude Bt denoted by numerals142,144,146, and148ofFIG. 8are not visible in the plot150ofFIG. 9, the very steep gradient patterns in the induced magnetic field amplitude Bt168,170,172, and174indicate respectively that the casings of the existing wells52,54,56, and58are nearby.

A simple alarm may be triggered if the induced magnetic field amplitude Btexceeds a certain value which indicates that the casing is too close to the BHA66. The alarm may indicate a potential collision between the drill bit70and a casing of one of the existing wells52,54,56, or58if the drilling continues unchanged. A driller controlling the BHA66may be prompted to stop and evaluate the situation upon the triggering of the alarm.

As indicated byFIGS. 7-9, the induced magnetic field amplitude Btis quite large if the BHA66is more than 1 m from the origin in the center of each plot. If the induced magnetic field90amplitude exceeds 150 nanoTesla (nT), then the BHA66is more than 1 m from the origin in the center of each plot. Because the value exceeds the minimum resolution of conventional MWD magnetometers, approximately 10 nanoTesla (nT), and because magnetometers with a resolution of 1 nanoTesla (nT) or smaller are available, the presently described technique may be performed using existing magnetometer technology.

The position of the BHA66relative to the casings of the existing wells52,54,56, and58may further be determined by measuring the induced magnetic field90components Bx(xm,ym) and By(xm,ym). Note that resolving the Bx-By components of the induced magnetic field90requires an independent measurement of the BHA66orientation, i.e. x-y, or North and East. Under normal conditions, the orientation is provided by a measurement of the Earth's magnetic field using the magnetometer94when the current76on the BHA66is not active. However, nearby steel casings of the existing wells52,54,56, or58may perturb the Earth's magnetic field and thus degrade the directional measurement, reducing the accuracy with which one may resolve the x-y directions.

Accordingly, an MWD gyro in the BHA66may additionally or alternatively be used to determine the direction, or a wireline gyro may be periodically run in the drill string attached to the BHA66to determine the x-y directions. The MWD gyro or the wireline gyro could be employed to calibrate the effect of the casings on the Earth's magnetic field or to directly determine orientation with respect to North. If the existing wells52,54,56, and58and the BHA66are slightly inclined, then a gravity tool face may be used to determine the x-y directions. In the foregoing discussion, it may be assumed that the x-y directions have been determined according to the above-described manners or any other appropriate manner.

FIGS. 10 and 11illustrate respectively the magnetic field components Bx(xm,ym) and By(xm,ym) over the region xmε[−1,1] and ymε[−1,1]. Turning first toFIG. 10, a 3-D plot176illustrates the magnetic field component Bx(xm,ym) over the region xmε[−1,1] and ymε[−1,1]. A legend178indicates magnetic field strength in microTesla (μT), which is illustrated along the height180of the 3-D plot176. A numeral182indicates the y-direction and a numeral184indicates the x-direction, such that a point186is located at (x,y)=(1,1), a point188is located at (x,y)=(1,−1), and a point190is located at (x,y)=(−1,−1). A numeral192marks the location of the BHA66in the center of the 3-D plot176.

Turning next toFIG. 11, a similar 3-D plot194illustrates the magnetic field component By(xm,ym) over the region xmε[−1,1] and ymε[1,1]. A legend196indicates magnetic field strength in microTesla (μT), which is illustrated along the height198of the 3-D plot194. A numeral200indicates the y-direction and a numeral202indicates the x-direction, such that a point204is located at (x,y)=(1,1), a point206is located at (x,y)=(1,−1), and a point208is located at (x,y)=(−1,−1). A numeral210marks the location of the BHA66in the center of the 3-D plot194.

FromFIGS. 10 and 11, it should be noted that there is additional information in the amplitudes and phases of the component data, which may be distinguished from the total induced magnetic field90amplitude. The total induced magnetic field90amplitude may be described according to the following equation:
Bt(xm,ym)=√{square root over (Bx(xm,ym)2+By(xm,ym)2)}{square root over (Bx(xm,ym)2+By(xm,ym)2)}  (7).

FIG. 12provides a schematic212which depicts a situation where the BHA66is located more closely to the casing of the existing well52than to any other of the existing wells54,56, or58. The magnetometer94within the BHA66measures the Bx and By components of the magnetic field90which surrounds the casing of the existing well52. In the schematic212ofFIG. 12, the x-axis is denoted by numeral60and the y-axis is denoted by the numeral62. A drift trajectory214shows a path, along which the BHA66slowly drifts from its original position at the origin due to slight errors in the MWD inclination measurements in the BHA66.

The situation depicted in schematic212ofFIG. 12may illustrate a manner of obtaining additional information from the individual magnetic field90components Bx(xm,ym) and By(xm,ym). Because the casing of the existing well52has the largest current84, the induced magnetic field90from this casing will be stronger than that of any other of the existing wells54,56, or58. Moreover, because the current84flows in the +z direction, both components of magnetic field90will be negative, such that Bx<0 and By<0.

Both the phases and amplitudes of Bx and By may provide additional information about the location of the BHA66with respect to the casings of the existing wells52,54,56, and58. For the purposes of plotting the magnetic field90components, it may be assumed that the magnetometer94in the BHA66moves along the drift trajectory214, represented by a line defined by y=m·x+b=0.2x. This may occur if the MWD inclination measurement of the BHA66is slightly erroneous, such that the vertical well trajectory drifts away from vertical with increasing depth. For a specific example, suppose that the new well drilled by the BHA66drifts 0.25 m in the x-direction and 0.05 m in the y-direction for every 10 m increase in depth. Such drift corresponds to an angle of about 1.4° deviation from vertical.

FIG. 13provides a schematic216which depicts geometry for estimating the direction and distance from the BHA66to the closest existing well52. The magnetometer94within the BHA66measures the Bx and By components of the magnetic field90which surrounds the casing of the existing well52. In the schematic216ofFIG. 13, the x-axis is denoted by numeral60and the y-axis is denoted by the numeral62.

By neglecting the effect of casings of the other existing wells54,56, and58, an apparent distance (Sa) and an apparent direction (γa) from the magnetometer94at the BHA66to the nearby casing of existing well52may be estimated. As illustrated in the schematic216, the BHA66is located at {right arrow over (r)}m=(xm,ym) and the casing of the existing well52is located at {right arrow over (r)}1=(x1,y1). Accordingly, an apparent direction to the casing can be derived from the induced magnetic field90components according to the following equation:

If the existing well52were the only casing, the above result would be exact, since the azimuthal magnetic field90is perpendicular to a radial vector which is directed from a line current to the observation point. As derived from the geometry depicted in the schematic216, the true direction (γ) from the BHA to the casing may be represented according to the following equation:

Turning next toFIG. 14, a plot218illustrates a change in angle over distance when the drift trajectory214is defined by y=0.2×. An ordinate220represents the direction in degrees and an abscissa222represents distance in meters (m). A curve224illustrates a change in apparent direction (γa) over distance from 0.5 m to 2.6 m, while a curve226illustrates a change in true direction (γ) over the distance from 0.5 to 2.6 m.

In the example shown by the plot218, the apparent direction (γa) is within 10° of the true direction (γ) over the range xmε[0.5, 2.6]. The difference results by neglecting the casings of the other existing wells54,56, and58, particularly the existing well54located at (x2,y2)=(0,2). Nonetheless, the apparent direction (γa) is sufficient information to steer the BHA66back toward the origin and away from the casing of the existing well52at (x1,y1)=(2,0).

FIG. 15is a plot228illustrating lines of constant apparent angle γa(xm,ym) for the area surrounding the casing of the existing well52at (x1,y1)=(2,0). An ordinate230indicates the y-coordinate value over a range of ymε[−1,1] and an abscissa232indicates the x-coordinate value over a range of xmε[0.5,2.6]. Each of the lines illustrated in the plot228shows a constant apparent angle γa(xm,ym) as a multiple of 10. Every third line is labeled accordingly. The plot228ofFIG. 15shows that the error in the apparent direction γa(xm,ym) reduces as the BHA66approaches this casing of the existing well52.

FIG. 16is a plot234illustrating the corresponding contour lines for the induced magnetic field90amplitude Bt(xm,ym) surrounding the casing of the existing well52at (x1,y1)=(2,0). An ordinate236indicates the y-coordinate value over a range of ymε[−1,1] and an abscissa238indicates the x-coordinate value over a range of xmε[−0.5,2.6]. Each contour line indicates an increase in magnetic field90amplitude Bt(xm,ym) in increments of 0.2 microTesla (μT) as the BHA66approaches this casing of the existing well52.

As indicated by the plot234, the magnetic field90amplitude Bt(xm,ym) lines are approximately circular near the casing of the existing well52, so that it is possible to invert for the approximate distance to the casing of the existing well52with the total induced magnetic field90. A first order approximation is given by

Sa=μ0⁢IC2⁢π⁢⁢B⁢⁢t,
where ICrepresents an estimate of the current84on the casing of the existing well52. The simplest approach is to allocate ¼thof the total current76(IZ) to the casing of the existing well52, namely IC=I(z)/4. The factor of ¼ is chosen because the BHA66is surrounded by the four casings of the existing wells52,54,56, and58.

Turning toFIG. 17, a plot240illustrates a change in distance from the BHA66to the casing of the existing well52when the drift trajectory214is defined by y=0.2x. An ordinate242represents the distance from the BHA66to the casing of the existing well52in meters (m) and an abscissa244represents distance in the x-direction in meters (m). A curve246illustrates a change in apparent distance (Sa) over distance in the x-direction from 0.5 m to 2.6 m, while a curve248illustrates a change in true distance (S) over distance in the x-direction from 0.5 m to 2.6 m. Further denoted in the plot240is a threshold distance250, which may trigger an alarm indicating that the BHA66is too close to another well.

The true distance (S) between the BHA66and the casing of the existing well52at (x1,y1)=(2,0) may be represented as S1=√{square root over ((x1−xm)2+(y1−ym)2)}{square root over ((x1−xm)2+(y1−ym)2)}. As mentioned above, the plot240illustrates the true distance in curve248and the apparent distance (Sa) in curve246for the same drift trajectory214, y=0.2x. The apparent distance (Sa) is an overestimate for x<1.4m because the other three casings of the existing wells54,56, and58reduce the magnetic field90amplitude around the origin. The apparent distance (Sa) is an underestimate for x>1.4 m as the BHA66approaches the casing of the existing well52at (x1,y1)=(2,0) because the current84on the casing will be greater than ¼thof the total current.

FIG. 18is a flowchart254for employing the apparent distance (Sa) for avoiding a collision with one of the existing wells52,54,56, or58. The flowchart254begins with step256, in which drilling begins in a field having at least one existing well such as the existing wells52,54,56, or58. In step258, magnetic ranging while drilling may be periodically or consistently employed generating the current76on the BHA66using the electric current driving tool74. The current76will enter the surrounding formation as the current82and run along the casing of one of the existing wells52,54,56, of58as the current84, which induces the azimuthal magnetic field90. In step260, the components of the magnetic field90, Bx and By, may be measured from the magnetometer94in the BHA66.

Step262involves estimating the apparent distance (Sa) and apparent direction (γa) using the first order approximation described above. As indicated by a decision block264, if the apparent distance (Sa) drops below the predetermined threshold distance250, then the process turns to step266. An alarm may alert the driller that the drill bit70of the BHA66is approaching a well casing, allowing the driller to take evasive action by steering in the direction opposite the apparent direction (γa). For example, if the threshold distance250is set at Sa=1 m, then the driller would be alerted at an alarm trigger distance252of x=1.2 m, which corresponds to a true distance of S1=0.8 m. Of course, the threshold distance250could be set to be a larger apparent distance (Sa). For example, if the threshold distance250were instead Sa=2 m, then the closest true distance would be S1=1. Returning to decision block264, if the apparent distance (Sa) remains above the threshold distance250, the process returns to step258to continue drilling.

As noted, the collision-avoidance solution above represents a first order solution for locating the BHA66with respect to the casings of the existing wells52,54,56, and58. The accuracy could be further improved by accounting for the current84on the casings of the existing wells54,56, and58in the inversion process, starting from the first order result. In addition, the currents84could be adjusted to reflect the relative distances from the BHA66to the casings of the existing wells52,54,56, and58. The apparent distance calculation may be improved by including an estimate of the conductance Gibetween the BHA66and any ithcasing. The conductance Giincreases as the distance between the BHA66and the ithcasing decreases. Accordingly, the current on the casing, Ii, increases. This effect may be included in the inversion by replacing the approximation for current84IC=I(z)/4 with an approximation that includes estimates for the conductances Gifor each existing well52,54,56, and58.

Alternatively, the first order solution may be practiced in other ways. For example, the apparent direction γa(xm·ym) may be plotted as inFIG. 15, and the total field amplitude Bt(xm,ym) may be plotted as inFIG. 16. The comparison of the two plots may provide a better estimate of the BHA66location, since only the (x,y) points where both conditions are satisfied are possible locations for the BHA66. A related approach using least squares will be described below.

Summarizing, the first order inversion process, which assumes a single well, involves estimating the apparent angle from the BHA to the cased well as

γa=tan-1⁡(-BxBy)
and estimating the apparent distance to the cased well according to the following equation:

In equation (10) above, the current ICis chosen depending on the situation. If there is only one cased well nearby, then a reasonable choice is IC≡I(0)(1+zm/L), where I(0) represents the current76generated at the insulated gap and where the magnetometer94is located at zm. If there are four casings nearby, as occurs when the BHA66is surrounded by the existing wells52,54,56, and58, then IC≡/(0)(1+zm/L) 4 is a reasonable choice. When the apparent distance Sadrops below a threshold value, the driller may be warned via an alarm of an impending collision with a cased well. The apparent angle γapoints toward the casing, and so the driller can avoid the collision by steering the drill bit in the opposite direction.

Using inversion and assuming a single cased well may apply to any arbitrary arrangement of cased wells. One may avoid a collision following the procedure described above. Knowing the location of the cased well is not required, as such information is not needed for Saor γa. It is not even necessary to know that there are any cased wells in the immediate vicinity, as the threshold alarm may indicate the proximity of a nearby cased well. Further, while the process has been illustrated with parallel wells, it may also be employed with non-parallel wells.

The above analyses assumed that the location of a casing of the existing wells52,54,56, or58may be unknown. If the positions of the existing wells52,54,56, and58are known, such data, in combination with measurements of the magnetic field90, may be used to locate the BHA66. The foregoing technique for locating the BHA66amid the existing wells52,54,56, and58involves calculating a theoretical magnetic field distribution and comparing the theoretical values to actual measurements of the magnetic field90. A least squares analysis may be employed for estimating the position of the BHA66.

The theoretical magnetic field that is measured at the magnetometer is denoted by {right arrow over (B)}(xm,ym)=(xm,ym){circumflex over (x)}+By(xm,ym)ŷ; where (xm,ym) refers to the position of the magnetometer94in the BHA66. For the purposes of illustrating the concept, simplifying assumptions about the theoretical model for {right arrow over (B)}(xm,ym) are employed. First, the BHA66and the casings of the existing wells52,54,56, and58are parallel or nearly parallel. Second, the positions of the existing wells52,54,56, and58are known. Third, resistivity of the surrounding formation is homogenous. Fourth, the current84on a casing of the existing wells52,54,56, or58may be calculated using the theoretical conductance between the BHA66and the casing. With a more sophisticated analysis, the above assumptions may be relaxed accordingly, but the underlying principles of the method will remain the same.

The present embodiment may explained by returning to view the geometry illustrated inFIGS. 4 and 5. From the geometry of theFIGS. 4 and 5, a resulting theoretical field {right arrow over (B)}(xm,ym) is plotted inFIGS. 7-11. The position of the BHA66may be assumed not well known, owing to accumulated errors in the standard MWD direction and inclination measurements. The actual measurement of the induced magnetic field90observed by the magnetometer94in the BHA66may be denoted as {right arrow over (β)}(x,y)=βx(x,y){circumflex over (x)}+βy (x,y)ŷ. Also, the actual position of the magnetometer94may be denoted as (x,y), which is treated as unknown. An objective of the present embodiment is to estimate (x,y) by comparing the actual magnetometer94measurement {right arrow over (β)}(x,y) to the theoretical model {right arrow over (β)}(xm,ym).

One approach for comparing measured or experimental values to theoretical values is to employ a least squares method, whereby the differences between the measured and theoretical values are minimized The quantity Q to be minimized may be defined according to the following relationship:

In equation (11) above, the actual position of the BHA66, (x,y), is an unknown quantity. Moreover, xmε[−2.6,2.6] and ymε[−2.6,2.6] are variables. To estimate the actual position of the BHA66, the objective is to minimize Q(xm,ym) on the xm-ymplane.

FIG. 19illustrates a 2-D plot268of the function Q(xm,ym) when the BHA66is at the origin, so that the true position of the magnetometer94in the BHA66is (x,y)=(0,0) and the measured values of the magnetic field90from the existing wells52,54,56, and58are βx(0,0)=0 and βy(0,0)=0. An ordinate270represents a range of ymε[−2.6,2.6] in the y-direction and abscissa272represents a range of xmε[−2.6,2.6] in the x-direction. The 2-D plot268for Q(xm,ym) includes contour lines274in increments of 20 nanoTesla (nT). The largest value plotted is 100 nanoTesla (nT). The location of the casings of the existing wells52,54,56, and58in the plot268are marked accordingly. The contour line closest to the origin is a minimum of Q(xm,ym), which has a value less than 20 nT within this area. If the magnetometer94is accurate to 20 nanoTesla (nT) and reads a value less than or equal to 20 nT, then the BHA66must be within ±0.5 m of the origin where the theoretical value for the magnetic field is zero. The more accurate the measurement, the better to estimate the actual location of the BHA66. Defining the magnetometer94accuracy as σBallows for the definition of a unit-less quantity ξ(xm,ym) as follows:
ξ(xm,ym)=Q(xm,ym)/σB(12).

FIGS. 20-24offer similar 2-D plots of the function Q(xm,ym) for different positions of the BHA66following the drift trajectory214of y=0.2x. Turning first toFIG. 20, a plot276of the function Q(xm, ym, ym) indicates the true position of the BHA66at (x,y)=(0.5,0.1). An ordinate278represents a range of ymε[−2.6,2.6] in the y-direction and abscissa280represents a range of xmε[−2.6, 2.6] in the x-direction. The location of the casings of the existing wells52,54,56, and58in the plot276are marked accordingly. The 2-D plot276for Q(xm,ym, ym) includes contour lines282in increments of 20 nanoTesla (nT). The largest value for a contour line is 100 nanoTesla (nT). The smallest value for a contour line is 20 nT, and it lies to the right of the origin, centered near (x,y)=(0.5,0.1). The area within this contour line indicates that the measured magnetic field is within 20 nT of the theoretical value for the magnetic field. This contour line2indicates that the BHA66is within the contour line centered on (x,y)=(0.5,0.1). However, it should be noted there are also two areas to the left of the origin that are also minima284of Q(xm,ym).

FIG. 21depicts a plot286of the function Q(xm,ym) where the true position of the BHA66is at (x,y)=(1.0,0.2). An ordinate288represents a range of ymε[−2.6,2.6] in the y-direction and abscissa290represents a range of xmε[−2.6,2.6] in the x-direction. The location of the casings of the existing wells52,54,56, and58in the plot286are marked accordingly. The 2-D plot286for Q(xm,ym) further includes contour lines292in increments of 20 nanoTesla (nT). The largest value plotted is 100 nanoTesla (nT).

As apparent in the plot286ofFIG. 21, there are three minima294,296, and298of Q(xm,ym). The minimum294to the right of the origin at (x,y)=(1.0,0.2) represents the true position of the BHA66, and is located to within ±0.05 m for measurement accuracy of 20 nanoTesla (nT). However, the two minima296and298to the left of the origin at (x′,y′)=(−0.90,1.15) and (x″,y″)=(−0.60,−1.10), respectively, are false positions or ghost images.

FIG. 22depicts a plot300of the function Q(xm,ym) where the true position of the BHA66at (x,y)=(1.5,0.3). An ordinate302represents a range of ymε[−2.6, 2.6] in the y-direction and abscissa304represents a range of xmε[−2.6, 2.6] in the x-direction. The location of the casings of the existing wells52,54,56, and58in the plot300are marked accordingly. The plot300for Q(xm,ym) includes contour lines306in increments of 20 nanoTesla (nT). The largest value plotted is 200 nanoTesla (nT).

As apparent in the plot300ofFIG. 22, there are four minima308,310,312, and314of Q(xm,ym). The minimum308to the right of the origin at (x,y)=(1.5,0.3) represents the true position of the BHA66, and is located to within ±0.05 m for measurement accuracy of 20 nanoTesla (nT). As in the plot286ofFIG. 22, the remaining minima310,312, and314are ghost images.

FIGS. 23 and 24illustrate plots of the function Q(xm,ym) when the BHA66is located at (x,y)=(2.0,0.4) and (x,y)=(2.5,0.5), respectively. Turning first toFIG. 23, the true position of the BHA66is (x,y)=(2.0,0.4). An ordinate318represents a range of ymε[−2.6,2.6] in the y-direction and abscissa320represents a range of xmε[−2.6, 2.6] in the x-direction. The locations of the casings of the existing wells52,54,56, and58in the plot316are marked accordingly. The 2-D plot316for Q(xm,ym) includes contour lines322in increments of 20 nanoTesla (nT). The largest value plotted is 200 nanoTesla (nT).

As apparent in the plot316ofFIG. 23, there are four minima324,326,328, and330of Q(xm,ym). The minimum324to the right of the origin at (x,y)=(2.0,0.4) represents the true position of the BHA66. However, the remaining minima326,328, and330are ghost images. Thus, a single measurement at one depth would not provide sufficient data to ascertain which minimum corresponds to the position of the BHA66and which minima are ghost images.

Similarly,FIG. 24depicts a plot322where the true position of the BHA66is at (x,y)=(2.5,0.5). An ordinate334represents a range of ymε[−2.6,2.6] in the y-direction and abscissa336represents a range of xmε[−2.6,2.6] in the x-direction. The locations of the casings of the existing wells52,54,56, and58in the plot332are marked accordingly. The 2-D plot332for Q(xm,ym) includes contour lines338in increments of 20 nanoTesla (nT). The largest value plotted is 200 nanoTesla (nT).

As apparent in the plot332ofFIG. 24, there are four minima340,342,344,346of Q(xm,ym). The minimum340to the right of the origin at (x,y)=(2.5,0.5) represents the true position of the BHA66. However, the remaining minima342,344,346are ghost images. Thus, a single measurement at one depth would not provide sufficient data to ascertain which minimum corresponds to the position of the BHA66and which minima are ghost images.

To distinguish the true location of the BHA66from the false positions or ghost images which may arise, a sequence of measurements may be obtained at different depths which may indicate the true position of the BHA66over the ghost images. Turning toFIG. 25, a plan view348shows the minima of Q(xm,ym) for BHA66at various depths. A legend350indicates the true position of the BHA66and three ghost images. An ordinate352represents a range of ymε[−3,3] in the y-direction and abscissa354represents a range of xmε[−3,3] in the x-direction. In the plan view348, the minima of Q(xm,ym) are plotted for increments of Δx=0.25 m, Δy=0.05 m for every 10 m increase in BHA66depth.

The initial position356of the BHA66is at the origin, (x,y)=(0,0), a logical starting point at the surface to drill another well amid the existing wells52,54,56, and58. Since the initial position356of the BHA66is known, the sequence of measurements versus depth may be used to differentiate the true trajectory358from the ghost trajectories360,362, and364. At the first measured depth (10 m), the minima of Q(xm,ym) which are plotted are labeled “1.” Among the points labeled “1”, the point labeled “1” in the true trajectory358may be more probably understood to be the true location of the BHA66than the first ghost trajectory360or the second ghost trajectory362because the step-out is smaller. Moreover, the step-out should be appreciated to be more consistent with an expected deviation from the BHA66drilling tendencies or MWD direction and inclination errors.

As the well is drilled, the true trajectory358follows a relatively straight line with relatively consistent increments in the position on the x-y plane. Meanwhile, the first ghost trajectory360and the second ghost trajectory362are curved and their increments are more erratic. Furthermore, the third ghost trajectory364does not even appear until the sixth depth measurement is made, and thus may clearly be eliminated as a ghost image. An interpreter could differentiate the true trajectory358from the ghost trajectories360,362, and364based on a plot such as the plot348.

FIGS. 26-28illustrate how additional information may clarify the interpretation and further distinguish the true trajectory from ghost trajectories which may arise. Turning first toFIG. 26, a plot366denotes the computed apparent direction

γa=tan-1⁡(-B⁢⁢xBy)
to the casing of the nearest well, existing well52, for the true trajectory358. In the plot366, a numeral368denotes the y-axis and a numeral370denotes the x-axis. Directional arrows372indicate the apparent direction (γa) to the nearest casing for each point along the true trajectory358and an arrow374indicates the movement of the true trajectory358. As illustrated in the plot366, the apparent positions and directions show a high degree of consistency with the casing of the existing well52located at (x1,y1)=(2,0). All of the directional arrows372point toward the casing at (x1,y1)=(2,0), beginning with the point labeled “1.”

FIG. 27depicts a plot376denoting the computed apparent direction

γa=tan-1⁡(-B⁢⁢xBy)
for each point of the ghost trajectory360. In the plot376, the numeral368denotes the y-axis and the numeral370denotes the x-axis. Arrows378indicate the movement of the ghost trajectory360and directional arrows372indicate the apparent direction (γa) to the nearest casing for each point along the ghost trajectory360.

As illustrated in the plot376, the apparent positions and directions for the ghost trajectory360are not as consistent as those associated with the true trajectory358. The inconsistencies are especially notable near the origin. For example, the first point, labeled “1,” is located to the left of the origin to (x,y)=(−0.55,0.60), and hence is thus further from the casing of the existing well52at (x1,y1)=(2,0) than the casing of the existing well54at (x2,y2)=(0,2). However, the directional arrow for point “1” points toward the casing of the existing well52. Thus, point “1” is clearly shown not to represent a part of the true trajectory358. Not until the sixth point in the ghost trajectory360does the directional arrow point toward the nearest casing, located at (x2,y2)=(0,2).

Similar conclusions may be drawn fromFIG. 28, which depicts a plot382denoting the computed apparent direction

γa=tan-1⁡(-B⁢⁢xBy)
for each point of the ghost trajectory362. In the plot382, the numeral368denotes the y-axis and the numeral370denotes the x-axis. Arrows384indicate the movement of the ghost trajectory362and directional arrows386indicate the apparent direction (γa) to the nearest casing for each point along the ghost trajectory362. As similarly illustrated in the plot376ofFIG. 27, in the plot382ofFIG. 28, the apparent positions and directions for the ghost trajectory362are not as consistent as those associated with the true trajectory358.

The data presented inFIGS. 25-28may greatly enhance the ability to avoid a collision with one of the existing wells52,54,56, or58. However, even without such data, a driller may be able simply to steer the BHA66away from a well casing. Suppose a driller were to make a decision as to which way to steer the BHA66based solely on the data illustrated in the plot316ofFIG. 23. The true position is (x,y)=(2.0,0.4), as indicated by the minimum324, and the ghost images are at (x′,y′)=(0.05,2.45), (x″,y″)=(0.05,−1.65), and (x′″,y ′″)=(−1.9,0.4), as indicated by the minima326,328, and330. Suppose an alarm based on the apparent distance has alerted the driller to an impending collision, but the driller does not have the historical sequence of measurements to tell him which minima of the plot316are ghosts. For all four possible positions indicated by the minima324,326,328, and330, the apparent direction remains the same, γa=−1.69 radians or −97°. Thus, the driller would know to steer at 83°, thus avoiding a collision with the casing, despite not knowing which minimum represents the true position and which minima represent ghost images.

FIG. 29is a flowchart388representing a general embodiment of the same approach which may be applied for other well configurations with any number of cased wells surrounding the BHA66. The principle remains the same, but the geometry may be different. In a first step390, the locations of cased wells versus depth are defined as {right arrow over (r)}i=(xi,yi,zi) for i={1, 2, 3, . . . , n} where {right arrow over (r)}irepresents the assumed location of the ithcased well and n represents the total number of nearby cased wells. The {{right arrow over (r)}i} will remain fixed throughout the procedure. The diameter of each cased wells is similarly defined as Di.

In step392, for a given depth zm, a location for the magnetometer94may be assumed as {right arrow over (r)}m=(xm,ym,zm), where xmand ymwill incremented over a range of values. In a subsequent step394, the conductance G, between the BHA66and each cased well may be computed according to the relationship

Gi=πσcosh-1⁡(Si/D),⁢⁢whereSi=(xm-xi)2+(ym-yi)2.
Similarly, the conductance may also be computed between each pair of cased wells. In both cases, the computations should take into account formation resistivity, cement resistivity, and bedding.

Turning next to step396of the flowchart388, the current84on each casing, Ii, may be computed for the assumed position of the BHA66, {right arrow over (r)}m. In step398, the magnetic field90at the magnetometer94for the assumed BHA66position {right arrow over (r)}mmay be computed according to the relationship

B->⁡(xm,ym,zm)=∑i=1n⁢B->⁢i⁡(xm,ym,zm)=∑i=1n⁢μ0⁢Ii⁡(zm)2⁢π⁢⁢Si2⁢n^×(r->m-r->i),
where {circumflex over (n)} represents a unit vector in the direction of the ithwell.

In step400, the induced magnetic field90may be measured with the three-axis magnetometer94to obtain the quantities {right arrow over (β)}(x,y,z)=βx(x,y,z){circumflex over (x)}+βy (x,y,z)ŷ+βz(x,y,z){circumflex over (z)}, where {right arrow over (r)}=(x,y,z) represents the actual position of the BHA66which is to be determined. Having obtained the magnetic field90measurements, in step402, the quantity

Q⁡(xm,ym,zm)=[β⁢⁢x⁡(x,y,z)-B⁢⁢x⁡(xm,ym,zm)]2+[β⁢⁢y⁡(x,y,z)-B⁢⁢y⁡(xm,ym,zm)]2+[β⁢⁢z⁡(x,y,z)-B⁢⁢z⁡(xm,ym,zm)]2
may be computed for the assumed location for the BHA66, {right arrow over (r)}m.

Continuing with step404of the flowchart388ofFIG. 29, the value for xmmay be incremented by Δx. Unless the maximum value for xmhas been reached, the process returns to the second step392. However, if the maximum value for xmhas been reached, the process continues to a ninth step406. In step406, the value for ymmay be incremented by Δy. Unless the maximum value for ymhas been reached, the process next returns to the second step392. However, if the maximum value for ymhas been reached, the process continues to a tenth step408.

Tenth step408involves locating the minima of Q(xm,ym,zm) for the given depth zm. In step410, a direction to the nearest casing for each minimum value of Q(xm,ym,zm) may be computed. Once computed, the apparent direction may be plotted on a plan view, such that

Continuing to drill in step412, measurement data may be obtained at a new depth zm+Δz. In step414which follows, the process returns to second step392to perform steps392-410with data obtained at the new depth. Finally, in step416, the position of the BHA66may be determined from the minima plotted in step410. Using both the positional information and the directional information, the true trajectory of the BHA66may be differentiated from the ghost trajectories of the minima

The approaches described above rely entirely on magnetic ranging data to resolve ambiguities that arise in estimating the actual position, (x,y), of the BHA66containing the magnetometer94when the objective function Q(xm,ym) has multiple minima Another approach may be to use the survey data to supplement the ascertainment of the actual position of the BHA66from the many ghost positions which may be represented by the minima in Q(xm,ym). As discussed above, when wells are tightly clustered, as in the example discussed above involving the existing wells52,54,56, and58, available survey data may not provide sufficient precision for drilling to continue within a desired margin of error. Nevertheless, the survey data may still contain additional information to resolve some ambiguities that may arise in the inversion of the ranging data.

The uncertainty in the position of a well bore resulting from survey errors can be described by a Gaussian probability distribution of the following form:

In equation (13) above, (x′,y′,z′) represents the well bore location obtained from the survey data, and σx, σy, and σzrepresent the standard deviations derived from measurement errors. It should be noted that the coordinate system, (x,y,z), is chosen such that there is null covariance between any two directions. Thus, the coordinate system to achieve such a result generally defines z along the wellbore, x in the vertical plane containing the wellbore, and y perpendicular to the x-z plane. As such, the coordinate system tends to decouple measured depth (“along hole”) errors, inclination errors, and azimuth errors.

An ellipsoid of uncertainty22(as depicted inFIG. 1) may be defined such that there is a given probability that the actual well falls inside the ellipsoid. Such an ellipsoid of uncertainty22may be centered on the location indicated by the survey data, (x′,y′,z′), may have semi-axes kσx, kσy, and kσz, and may be described according to the following equation:

By way of example, there is a 20% probability that the well lies within the ellipsoid defined by equation 14 when k=1. Similarly, there is an 86% probability that the well lies within the ellipsoid defined by equation 14 when k=2.

For the case of nearly parallel, vertical wells, the “along hole” errors correspond to σz, while the inclination and direction errors may combine to affect σxand σy. Because the relative angle between the BHA66and a cased well is small, an error in depth does not translate to a significant error in the x or y directions, in which there may be a risk of a collision. Hence, the probability distribution may be reduced to two dimensions (x,y) at any given depth z. Although not necessarily true in general, it may also be assumed that σx=σy=σ for simplicity. The probability density function at a given depth z may be defined by the following equation:

The three dimensional ellipsoid may reduce to a two dimensional circle, as defined by the following equation:
(x−x′)2+(y−y′)2=(kσ)2(16)

For such a special case, the probability is given by 1−exp(−0.5 k2). Thus, there is a 39% probability that the well lies within the circle defined by equation (16) when k=1, and a 95% probability that the well lies within the ellipsoid defined by k=2.45.

FIG. 30Aillustrates the situation described above with a well placement schematic418. The well placement schematic418depicts the predicted location of the BHA66relative to an ithcased well98. The numeral60represents the x-axis, while the numeral62represents the y-axis. The survey data predicts the BHA66location to be {right arrow over (r)}′=(x′,y′), with a one sigma circle420of radius σ centered on {right arrow over (r)}′. The survey data for the ithcased well98indicates that it is located at {right arrow over (ri′)} and hence the two surveys predict that the separation between the BHA66and the ithcased well98is {right arrow over (Si′)}={right arrow over (r′)}−{right arrow over (ri′)}. If the only uncertainty came from the BHA66survey, but the position of the ithcased well was known exactly, then one would need |{right arrow over (Si′)}|≧2.456 for a 5% probability of collision with the cased well. However, the above equation is true only with perfect knowledge of the location of the ithcased well98. Equations to here . . . .

In reality, the position of the ithcased well98is also described by a Gaussian probability distribution with an uncertainty, σi, associated with it. Hence, the actual condition for a 5% probability of a collision may be described according to the following equation:
|{right arrow over (Si′)}|≧2.445√{square root over (σ2+σi2)}  (17).

The uncertainty of the ithcased well98may be accounted for in the Gaussian probability distribution with the following equations:

σ->σ~=σ2+σi2;(18)F⁡(x,y)=1(2⁢π)⁢σ~2⁢exp⁢{-(x-x′)22⁢σ~2-(y-y′)22⁢σ~2}.(19)
Equation (18) combines the standard deviation for the BHA66with the standard deviation for a cased well to obtain an effective standard deviation {tilde over (σ)}. Equation (19) expands the width of the Gaussian probability distribution to include the uncertainties from the surveys of the cased wells. In equation (19), the most likely position for the BHA66is still the survey result, {right arrow over (r′)}.

FIG. 30Bdepicts the actual position of the BHA66in a well placement schematic422. In the well placement schematic422, the numeral60represents the x-axis, while the numeral62represents the y-axis. The BHA66is actually located at {right arrow over (r)} which, according to the Gaussian probability distribution, has a 39% probability of being in the one sigma circle420centered on {right arrow over (r′)}. The true location for the ithcased well98is {right arrow over (ri)}, and the true separation between the BHA66and the ithcased well98is {right arrow over (Si)}={right arrow over (r)}−{right arrow over (ri)}. However, to proceed with the analysis it may be assumed that the ithcased well98is actually located at a point of maximum probability424, such that {right arrow over (ri)}={right arrow over (r′)}. While this assumption is not true in general, the uncertainty in the separation between the BHA66and the ithcased well has been accounted for by equations (18) and (19). Alternatively, a Gaussian probability distribution function for each cased well can be used with that for the BHA66. However, this alternative approach only adds to the mathematical complexity. The simpler approach using equations (18) and (19) adequately illustrates the principle.

FIGS. 31 and 32depict two views of a Gaussian probability function for the magnetic ranging illustrated inFIG. 21. Considering that it is desirable to resolve magnetic ranging ambiguities using the survey data while including the uncertainties in the survey data, a Gaussian probability function as given by equations (18) and (19) may be combined with the magnetic ranging illustrated inFIG. 21. RecallingFIG. 21, there are three possible locations for the BHA66derived from the quantity Q(xm,ym). One location is the true position at {right arrow over (r)}=(1.0,0.2), while the other two locations are ghosts.

Turning toFIG. 31, a 3-D probability density plot426illustrates probability428from 0 to 1 in increments of 0.1 for the locations of the existing wells52,54,56, and58and the BHA66. A numeral430indicates the y-direction over the range ymε[−2.6,2.6] and a numeral432indicates the x-direction over the range xmε[−2.6,2.6], such that a point434is located at (x,y)=(2.6,2.6), a point436is located at (x,y)=(2.6,−2.6), and a point438is located at (x,y)=(−2.6,−2.6). The locations of the existing wells52,54,56, and58are represented by a probability of 1, as such data is assumed to be known. The casing diameters for the existing wells52,54,58, and58are shown inFIG. 31, while the Gaussian probability density is shown for the BHA66. A peak amplitude440of the probability density distribution of the location of the BHA66is normalized to 1, representing survey data which may be available predicting the BHA66location as {right arrow over (r′)}=(1.5,0.5) with an uncertainty of {tilde over (σ)}=1.

FIG. 32depicts a probability density plot442corresponding to the 3-D probability density function plot426ofFIG. 31. The probability density plot442similarly illustrates the location of a one sigma circle444, which indicates a high probability of the location of the BHA66. The x-axis60indicates the x-direction over a range xmε[−2.6,2.6] and the y-axis62indicates the y-direction over a range ymε[−2.6,2.6]. The probability density plot442further indicates the location of the existing wells52,54,56, and58. The one sigma circle444encircles the casing of the existing well52located at {right arrow over (r1)}=(2,0), indicating a high probability of a collision between the BHA66and the existing well52. Because the probability density data is provided by survey data alone, the new well being drilled by the BHA66could not be drilled with certainty if only survey data were available.

The survey data can be combined with the magnetic ranging information to improve the knowledge of the BHA66location. The probability distribution can be modified to include the magnetic ranging data by weighting the Gaussian probability density by ξ(x,y) as indicated by the following relationship:

FIG. 33depicts a plot446illustrating the weighted probability density function H (xm,ym), for {right arrow over (r′)}=(1.5,0.5) and {tilde over (σ)}=1, when the true BHA position is at {right arrow over (r)}=(1.0,0.2). An ordinate448represents a range of ymε[−2.6, 2.6] in the y-direction and abscissa450represents a range of xmε[−2.6,2.6] in the x-direction. The location of the casings of the existing wells52,54,56, and58in the plot446are marked accordingly. Weighted probability density function contour lines452indicate three maxima454,456, or458. However, as apparent in the plot446, the maxima454vastly outweighs the other two maxima456and458. Thus the maxima454clearly represents the true location of the BHA66, while the remaining locations456and458are clearly ghost images.

FIG. 34represents a flowchart460illustrating a process for employing the weighted probability density function of equation (20) to estimate the location of the BHA66when the locations of the existing wells52,54,56, and58are known. In a first step462, the locations of cased existing wells52,54,56, and58versus depth may be defined as {right arrow over (ri)}=(xi,yi,zi) for i={1, 2, 3, . . . , n}, where {right arrow over (ri)} represents the assumed location of the ithcased well98and n represents the total number of nearby cased wells. The {{right arrow over (ri)}} will remain fixed throughout the procedure. The diameter of each cased well is similarly defined as Di. In step464, the new well is drilled using the BHA66down to a depth zm.

In step466, MWD survey data may be used to obtain the probability distribution function

F⁡(x,y)=12⁢π⁢σ~2⁢exp⁢{-(x-x′)22⁢σ~2-(y-y′)22⁢σ~2}
at the given depth zm, where {right arrow over (r′)}=(x′,y′,zm) represents the most likely position of the BHA66determined by the survey data, where {tilde over (σ)}=√{square root over (σ2+σi2)}, σ represents the standard deviation in the x-y plane for the BHA66, and σirepresents the standard deviation for survey data for the cased wells. Step468, which follows, involves assuming a location for the magnetometer94in the BHA66, {right arrow over (r)}m=xm, ym, zm, for the given depth zm. As discussed further in the flowchart460ofFIG. 34, xmand ymwill be incremented over a range of values.

With further reference to the flowchart460ofFIG. 34, in step470, the conductance Gibetween the BHA66and each cased well may be computed according to the relationship

Gi=πσcosh-1⁡(Si/D),⁢⁢whereSi=(xm-xi)2+(ym-yi)2.
Similarly, the conductance may also be computed between each pair of cased wells. In both cases, the computations should take into account formation resistivity, cement resistivity, and bedding. In step472, the current84on each casing, Ii, may be computed for the assumed position of the BHA66, {right arrow over (rm)}.

In step474, the magnetic field90at the magnetometer94for the assumed BHA66position {right arrow over (rm)} may be computed according to the relationship

B->⁡(xm,ym,zm)=∑i=1n⁢B->⁢i⁡(xm,ym,zm)=∑i=1n⁢μ0⁢Ii⁡(zm)2⁢π⁢⁢Si2⁢n^×(rm->-ri->),
where {circumflex over (n)} represents a unit vector in the direction of the ithwell98. In step476, the induced magnetic field90may be measured with the three-axis magnetometer94to obtain the quantities {right arrow over (β)}(x,y,z)=βx(x,y,z){circumflex over (x)}+βy (x,y,z)ŷ+βz(x,y,z){circumflex over (z)}, where {right arrow over (r)}=(x,y,z) represents the actual position of the BHA66which is to be determined The standard deviation in the measured magnetic field components is σB. Having obtained the magnetic field90measurements in step476, in step478, the quantity

ξ⁡(xm,ym,zm)=Q⁡(xm,ym,zm)/σB=[β⁢⁢x⁡(x,y,z)-Bx⁡(xm,ym,zm)]2+[β⁢⁢y⁡(x,y,z)-By⁡(xm,ym,zm)]2+[β⁢⁢z⁡(x,y,z)-Bz⁡(xm,ym,zm)]2σB
may be computed for the assumed location for the BHA66, {right arrow over (r)}m.

Continuing to step480of the flowchart460ofFIG. 34, the value for xmmay be incremented by Δx. Unless the maximum value for xmhas been reached, the process returns to the fourth step468. However, if the maximum value for xmhas been reached, the process continues to an eleventh step482. In step482, the value for ymmay be incremented by Δy. Unless the maximum value for ymhas been reached, the process next returns to the fourth step468. However, if the maximum value for ymhas been reached, the process continues to a twelfth step484.

In step484, the Gaussian probability density function F(xm,ym) is divided by ξ(xm,ym) to obtain the weighted probability distribution

H⁡(xm,ym)=F⁡(xm,ym)ξ⁡(xm,ym).
Using the weighted probability distribution H(xm,ym) calculated in step484, in step486, the minima of H(xm,ym) may be located for the given depth zmwhich corresponds to the most probable location for the BHA66. Continuing to drill in step488, measurement data may be obtained at a new depth zm+Δz, before returning to the fourth step468to perform steps468-486with data obtained at the new depth. From the data obtained in the flowchart460, the position of the BHA66may be estimated by locating the true position as distinguished from any ghost images which may arise.

Another approach to finding the ‘best estimate’ for the location of the BHA66is to use a method described in U.S. Pat. No. 6,736,221, assigned to Schlumberger Technology Corporation, incorporated by reference herein [NOTE: we may not be able to incorporate this patent by reference; the cited patent incorporates matter by reference in the background. I do not know whether the incorporated matter is essential.]. This technique requires covariance matrices for the positions calculated from the ranging and survey data. The covariance matrices can be evaluated by standard methods.

The previous example was based in part on the assumption that the uncertainty in the positions of the existing wells52,54,56, and58may simply be included in the uncertainty for the BHA66location, and that the locations of the existing wells52,54,56, and58may be assumed to be at the most probable locations provided by the survey data for the cased wells, i.e. {right arrow over (ri)}={right arrow over (ri′)}. In the manner described above, the magnetic ranging data used to compute Q(xm,ym) derived from a model in which the cased well locations are assumed to be known. In a more general case, however, this assumption may be substituted by describing the locations of the cased wells using Gaussian probability distributions.

For example, the ithcased well98may have a Gaussian probability distribution of the form represented by the following equation:

In equation (21) above, σirepresents the standard deviation, and {right arrow over (r′i)}=(x′i,y′i) represents the survey position of the ithcased well98, which corresponds to the most probable location of the ithcased well98. For simplicity, the probability distributions are assumed to be symmetric, i.e. σix=σiy=σi.

FIGS. 35A and 35Bmay illustrate the geometry used in estimating the location of the BHA66using equation (21). Turning first toFIG. 35A, a well placement schematic490depicts the predicted location of the BHA66relative to the ithcased well98. The numeral60represents the x-axis, while the numeral62represents the y-axis. The survey data for the cased well98indicates that r′iis the most likely location for it, which is surrounded by a one sigma circle492. Likewise, survey data for the BHA66indicates that {right arrow over (r′)} is its most likely location of the BHA66, which is surrounded by a one sigma circle494. The relative displacement between the BHA66and the cased well is thus S=

In contrast,FIG. 35Bdepicts a well placement schematic496represents the actual location of the BHA66and the actual location of the ithcased well98. The numeral60represents the x-axis, while the numeral62represents the y-axis. As indicated in the well placement schematic496, the ithcased well98is actually at a different location, {right arrow over (ri)}=(xi,yi), and the BHA66is actually at a different location {right arrow over (r)}=(x,y). The relative displacement between the BHA66and the ithcased well98is {right arrow over (Si)}={right arrow over (r)}−{right arrow over (ri)}. Because the magnetic field90will be different for the two cases, i.e. when the cased well is at {right arrow over (r′i)} or {right arrow over (ri)}, the procedure is more complex.

The Monte Carlo method provides one method for combining two or more probability distributions with magnetic ranging in order to avoid a collision between the BHA66and a cased well, and to improve the knowledge of the relative positions of the BHA66and any cased wells, such as the existing wells52,54,56, or58. The Monte Carlo method is a well known computational process where random numbers and a large number of calculations are performed to model a physical process. Modern computers are capable of performing large numbers of calculations rapidly. To apply the Monte Carlo method to this particular problem, a set of values is chosen for the locations of the n nearby cased wells (i.e., for {{right arrow over (r1)}, {right arrow over (r2)}, {right arrow over (r3)}, . . . , {right arrow over (rn)}, }). The procedure described by the steps of the flowchart460ofFIG. 34from step462to step486may then be executed. The magnetic field90may be calculated for various possible positions of the BHA66given the set of values for {{right arrow over (r1)}, {right arrow over (r2)}, {right arrow over (r3)}, . . . , {right arrow over (rn)}, }. The quantity ξ(xm,ym) may be calculated and used to weight the probability distribution for the BHA66. The result, Hi(xm,ym), may be recorded or stored (the subscript “1” indicates that this is the first calculation). Then a different set of values for {{right arrow over (r1)}, {right arrow over (r2)}, {right arrow over (r3)}, . . . , {right arrow over (rn)}, } may be chosen, and the procedure described by the steps of the flowchart460ofFIG. 34from step462to step486may then be executed again. The result, H2(xm,ym), may be recorded or stored. The process may be repeated many times, but with the proviso that the probability distributions Fi(x′i,y ′i) are honored by the values chosen for {{right arrow over (r1)}, {right arrow over (r2)}, {right arrow over (r3)}, . . . , {right arrow over (rn)}, }

For example, 68% of the random values chosen for the location of the ithcased well98, located at {{right arrow over (ri)}}, should fall within the circle of radius σithat is centered on the point {right arrow over (r′i)}. After a sufficiently large number of calculations (p) are performed to achieve statistical accuracy, the quantity described according to the following equation is calculated:

The results of the equation above may be plotted in a manner similar to that shown by the plot446ofFIG. 33. The greatest of the maxima of H(xm,ym) corresponds to the best estimate for the location of the BHA66amongst the n cased wells, and takes both the probability distributions and the magnetic ranging data into account. It should be appreciated that the same techniques used for determining the position of the BHA66relative to the n cased wells may also be used to determine the position of the n cased wells relative to the BHA66. Thus, using MWD direction and inclination measurements from the BHA66, combined with the above-described methods of determining apparent distance and direction to the n cased wells, the position of the n cased wells may be similarly determined.

FIG. 36illustrates the procedure discussed above with a flowchart498. In a first step500, a for-do loop from 1 to p may be initialized by setting j=1 and choosing a value for p to achieve statistical accuracy. In step502, a set of random values for the locations of the n cased wells, {{right arrow over (r1)}, {right arrow over (r2)}, {right arrow over (r3)}, . . . , {right arrow over (rn)}, }, may be chosen, such that the random values honor the probability distributions {F1(x′1,y′1), F2(x′2,y′2), F3(x′3,y′3), . . . , Fn(x′n,y′n)}. Step504involves executing the procedure described by the steps of the flowchart460ofFIG. 34from step462to step486.

Having obtained a result for Hj(xm,ym) in step504, the result Hj(xm,ym) may be recorded and stored in a subsequent step506. In step508, the variable j may be incremented by 1. If j=p then the process continues to step510. Otherwise, the process returns to step502. In step510, the quantity

H⁡(xm,ym)=1p⁢∑j=1p⁢Hj⁡(xm,ym)
may be calculated, and in step512, the greatest of the maxima of

H⁡(xm,ym)=∏j=1p⁢Hj⁡(xm,ym)
may be ascertained. As discussed above, the greatest of the maxima of

H⁡(xm,ym)=1p⁢∑j=1p⁢Hj⁡(xm,ym)
represents a most probable position of the BHA66relative to the n cased wells.

Another application is determining the location of a cased well that has inaccurate survey data or no survey data. For example, old cased wells may have been surveyed with old and less accurate equipment, or the well surveys may have been lost, or the wells may not have been surveyed at all. When drilling a new well in the proximity of such an existing well, magnetic ranging while drilling and the MWD survey data from the well being drilled can be used to establish the cased well's location. Magnetic ranging can determine the relative displacement {right arrow over (S)}={right arrow over (r′)}−{right arrow over (rc)} of the cased well to the well being drilled. The MWD measurements provide data for the well being drilled, i.e. {right arrow over (r′)}—the survey position. Hence, the location of the cased well {right arrow over (rc)} is determined from {right arrow over (rc)}={right arrow over (r′)}−{right arrow over (S)}.

While these methods have been demonstrated for wells that are essentially parallel, this has been done only to simplify the equations and to provide a clear understanding of the technique. The condition of parallel wells is not essential for these methods to be applied. In particular, techniques for using magnetic ranging while drilling as applied to non-parallel wells are described in described in Published Application No. US 2007/016426 A1, Provisional Application No. 60/822,598, application Ser. No. 11/833,032, and application Ser. No. 11/781,704, each of which is assigned to Schlumberger Technology Corporation and incorporated herein by reference.

Moreover, the probability distribution functions for the well position may be three-dimensional, using arbitrary orientations of the ellipsoids for the cased wells and for the well being drilled. The probability distributions need not be Gaussian, although these are commonly used for describing oil and gas wells. Additionally, as discussed above, the above description illustratively discusses vertical wells only to simplify the mathematical analysis. When the wells are vertical, magnetic fields 90 which are induced on around the casings of the existing wells52,54,56, and58lie in the x-y plane, while the electric currents on the BHA66and casings of the existing wells52,54,56, and58flow in the ±z-direction. However, it is not necessary in general for the existing wells52,54,56, and58to be vertical or exactly parallel. The magnetic fields induced on a non-vertical well that is not parallel to the BHA can be modeled using the techniques described in the patent applications referenced above.