Continuous measurement-while-drilling surveying

[0002] This U.S. Patent Application is a continuation of U.S. patent application Ser. No. 09/790,591 filed Feb. 23, 2001, now U.S. Pat. No. 6,668,465 issued Dec. 30, 2003 entitled “Continuous Measurement-While-Drilling Surveying”.

[0003] The present invention relates to continuous measurement-while-drilling surveying methods and apparatuses.

[0004] To obtain hydrocarbons such as oil and gas, boreholes are drilled by rotating a drill bit attached at a drill string end. A large proportion of the current drilling activity involves directional drilling, i.e., drilling deviated and horizontal boreholes to increase the hydrocarbon production and/or to withdraw additional hydrocarbons from the earth's formations. Modem directional drilling systems generally employ a drill pipe having a drill bit at the bottom that is rotated by a drill motor or a “mud motor”. Pressurized drilling fluid or “mud” or “drilling mud” is pumped into the drill pipe to rotate the drill motor and to provide lubrication to various members of the drill string including the drill bit. The drill bit and drill motor form part of what is known as the bottom hole assembly (“BHA”). As required the drill pipe is rotated by a prime mover, such as a motor, to facilitate directional drilling and to drill vertical boreholes.

[0005] Measurement-While-Drilling (MWD) surveying for directional and horizontal drilling processes is performed to provide the orientation and the position of the BHA [Conti, 1999]. Azimuth, the inclination and the tool face angles determine the orientation of the BHA, while latitude, longitude and altitude determine the position of the BHA. The altitude directly determines the true vertical depth of the BHA. State of the art MWD surveying techniques are based on magnetic surveying which incorporates three-axis magnetometers and three-axis accelerometers arranged in three-mutually orthogonal directions. The three-axis accelerometers monitor the Earth gravity field to provide the inclination and the tool face angles. This information is combined with the magnetometer measurements of the Earth magnetic field to provide the azimuth [FIG. 1, Russel et al., 1979].

[0006] The magnetic surveying system determines the BHA orientation at certain survey stations with the assumption that the error which modifies the Earth's magnetic field vector at the surveying instrument is in the direction of the borehole [Russel et al., 1979]. This assumption is justified by installing these magnetometers inside a non-magnetic housing. Such housing system necessitates the use of non-magnetic drill collars around the surveying equipment at a cost approaching $30,000 per single installation [Rehm et al., 1989].

[0007] Although simple, the magnetic surveying system suffers from several inadequacies especially within the drilling environment [Thorogood, 1990; Thorogood et al., 1986]. The presence of downhole ore deposits deviate the measurements of the Earth magnetic field even with the non-magnetic drill collars surrounding the surveying instruments. In addition, magnetic surveying tools located in non-magnetic drill collars are subject to the influence of the other steel components of the drill string. It has been shown that drill string induced surveying error increases with inclination and as the borehole direction approaches the east-west azimuth. Drill string magnetic interference is particularly noticeable when inclination exceeds 30°. Although it has been reported that the effect of the drill string magnetic interference could be reduced (mitigated but never entirely eliminated) by running long lengths of non-magnetic materials above and below the survey instruments [Grindord et al., 1983], this solution could affect the cost benefits of the horizontal drilling technology. The effect of drill string magnetic interference as well as the presence of downhole ore deposits can neither be quantified nor compensated. In addition to these two effects, geomagnetic influences play an important role in the accuracy of the magnetic surveying system [Parkinson, 1983]. Geomagnetic influences are defined by the variation of the dip angle, the declination and the total magnetic field strength with respect to time. The dip angle is the angle between the direction of the Earth's magnetic field and the horizontal plane. The declination is the angle between the magnetic north and the true north. It was recorded that during any given day at a random location the standard deviations of the dip angle, the declination and the magnetic field strength were 0.3°, 0.9° and 0.3 μ-Tesla respectively [Parkinson, 1983]. The variation in the geomagnetic field is quite significant in relation to the performance capabilities of the magnetic surveying tools currently used. Therefore, geomagnetic effects must be taken into account when considering absolute survey accuracy and impose definite limitations on the accuracy levels that can ultimately be achieved.

[0008] Several investigations have been reported to improve the magnetic surveying accuracy. Shiells and Kerridge (2000) introduced the interpolated in-field referencing method in which absolute local geomagnetic field data is determined from spot measurements of the Earth's magnetic field. These measurements are taken at local site which is sufficiently close to the drilling site so that the measurement data is indicative of the Earth's magnetic field at the drilling site but is sufficiently remote from the drilling site so that the measurement data is unaffected by magnetic interference from the drilling site. This technique compensates for the geomagnetic influences but the measurements taken by the magnetometers still suffer from the drill string magnetic interference. In addition, measurements of the geomagnetic field in a local site might not be applicable in an area that has many drilling rigs. Moreover, such technique adds more complexity that may give rise to additional cost.

[0009] Hartmann (1998) introduced a method to improve the accuracy of borehole surveying. This method is based on determining the uncertainty in the magnetic field measurements by comparing the measured magnetic field with a corresponding known value at specific location. However, this method does not take into account two phenomena. The first phenomenon is the increase of the drill string magnetic interference at high inclinations. The second is the presence of ore deposits that can be found randomly downhole. These two phenomena add additional errors that are not taken into consideration when referencing the theoretical value of the magnetic field to the measured value at the known location.

[0010] Dipersio (1995) introduced a new method for the compensation of geomagnetic influences. This method depends on matching both the calculated magnetic field strength and the calculated dip angle to their corresponding nominal values at a particular geomagnetic location. This results in determining an error-free value of the axial component of the magnetic field which is directly used to determine the azimuth. However, at high inclinations, the drill string magnetic interference generates additional errors in the cross-axial directions.

[0011] Several other techniques have been described for magnetic MWD surveying [Nicholson, 1995; Engebretson, 1992; Helm, 1991; van Dongen et al., 1987; Trowsdale, 1982]. However, it is clear that the major weaknesses of the present directional sensing instruments stem from the use of magnetometers for monitoring the azimuth and from the hostile environment in which these devices work. The problem encountered with the use of magnetometers is the presence of massive amount of steel around the drilling rig. The abundance of ferromagnetic material necessitates the separation of the magnetometers by non-magnetic drill collars. Aside from the cost of utilizing non-magnetic drill collars, their use introduces a second problem. Since the non-magnetic drill collars impose an additional weight on the drill bit, the surveying tools are separated from the bearing assembly and the drill bit by about 30 feet [Conti, 1989; Rehm at al., 1989]. Elimination of the non-magnetic drill collars could reduce the distance between the instrument package and the drill bit by approximately 50%. The third problem associated with the use of magnetometers is their lack of reliability when used underground due to the deviation of the Earth's magnetic field from ore deposits.

[0012] It is known to use strapdown inertial navigation systems for borehole surveying using probes. Borehole surveying is performed after finishing the drilling process, or sections of the drilling process, for the purpose of quality assurance. U.S. Pat. No. 4,542,647 describes a borehole inertial guidance system which utilizes two ring laser gyroscopes having their sensitive axes normal to the borehole axis. In this case, azimuth is determined by a synthetic rotation signal provided by three-axis accelerometers. While this configuration may be suitable for borehole mapping, it fails to be useful for MWD surveying. There is considerable rotation about the tool-spin axis in a MWD process making the use of a synthetic signal to monitor this rotation and therefore azimuth useless.

[0013] In addition, ring laser gyros are unsuitable for MWD processes because they cannot withstand the harsh conditions encountered near the BHA during the drilling process. In particular, the vibration and shock forces associated with the progress of the drill bit impairs the effectiveness of ring laser gyros.

[0014] Therefore, there is a need in the art for a continuous measurement-while-drilling surveying method and apparatus which mitigates the difficulties of the prior art.

[0015] The present invention relates to the replacement of the MWD magnetic surveying system with the technology of inertial navigation. Inertial navigation systems (INS) determine the position, velocity and orientation of a moving body in three-dimensional space by integrating the measured components of the acceleration and the angular velocity of the body [Titteron et al., 1997]. The measurements of acceleration and angular velocity are provided by three-axis accelerometers and three-axis gyroscopes arranged in three mutually orthogonal directions.

[0016] Currently available 3-gyroscopes/3-accelerometers INS technology cannot provide high accuracy and small dimensions within the same setup and therefore at present, the direct installation of such system inside the BHA to replace the set of magnetometers is not feasible. Thus, the traditional high accuracy INS setup is not applicable downhole due to the limited space inside the BHA and the relatively large dimensions of the gyroscopes.

[0017] Therefore, herein disclosed is a method and apparatus for downhole strapdown INS utilizing only two novel gyroscopes. One gyroscope has its sensitive axis along the tool spin axis while the other has its sensitive axis normal to the tool spin axis and along the forward direction. This results in the removal of the costly non-magnetic drill collars because the gyroscopes are insensitive to magnetic fields. In one embodiment, the two gyroscopes are mounted inside the bearing assembly and are fiber optic gyroscopes (FOG). FOGs are suitable for the downhole drilling application because they do not contain moving parts, thus providing high reliability with no need for frequent calibration or maintenance [Kim, 1994; Merhav, 1993; Lefevre, 1993]. In addition, the FOGs exhibit low environmental sensitivity since they can withstand relatively high temperatures, shocks and vibrations [Noureldin et al., 1999; Noureldin et al., 2000]. Moreover, currently available FOGs are of small size, with a drift rate less than 1°/hr, long mean time between failure (60,000 hrs), no gravitational effects and excellent immunity to vibration and shock forces.

[0018] In one embodiment, the FOG that has its sensitive axis along the tool spin axis is prepared as torus to accommodate the flow of drilling fluid through the bearing assembly. The 3 accelerometers—2 FOG system is capable of providing continuous MWD surveying during the whole drilling process, up until the horizontal section of the well, instead of MWD station-based surveying. In addition, continuous MWD surveying is possible through an initial radical section of the well with 3 accelerometers and a single FOG which has its sensitive axis aligned with the tool spin axis. The initial radical section of the well may be that portion with an inclination less than about 45°, preferably less than about 30° and more preferably less than about 20°.

[0019] Therefore, in one aspect, the invention comprises a continuous measurement-while-drilling surveying apparatus for surveying the drilling progress of a bottom hole assembly (“BHA”) comprising:

[0020] (a) a housing adapted to be part of the BHA, said housing having a tool-spin axis;

[0021] (b) a first fiber-optic gyroscope mounted within the housing for generating a first angular rotation signal representing angular rotation of the BHA about the tool-spin axis;

[0022] (c) a second fiber-optic gyroscope mounted within the housing for generating a second angular rotation signal representing angular rotation of the BHA about an axis normal to the tool-spin axis;

[0023] (d) accelerometer means for generating three acceleration signals representing the components of acceleration of the BHA along three mutually orthogonal axes;

[0024] (e) first processing means responsive to the acceleration signals for determining the angle of the BHA away from the vertical and for generating a third angular rotation signal representing rotation of the BHA about an axis normal to the sensitive axes of the first and second gyroscopes;

[0025] (f) second processing means responsive to the first, second and third angular rotation signals and the acceleration signals for transforming signals representing movement of the BHA in a BHA coordinate system to a earth local-level coordinate system;

[0026] (g) third processing means operatively connected to the second processing means for determining the orientation of the BHA, determining the velocity changes of the BHA, updating the velocity components of the BHA and updating the position components of the BHA.

[0027] The first gyroscope may be toroidal. The first, second and third processing means may be a general purpose computer programmed with appropriate software, firmware, a microcontroller, a microprocessor or a plurality of microprocessors, a digital signal processor or other hardware or combination of hardware and software known to those skilled in the art.

[0028] In another aspect of the invention, the invention is a method of continuous MWD surveying of a wellbore which includes, on completion, a vertical section, a radical section and a horizontal section, using a BHA comprising a first gyroscope having its sensitive axis aligned with the tool spin axis, a second gyroscope having its sensitive axis on an axis normal to the tool spin axis and three accelerometers arranged in three mutually orthogonal directions, said method comprising the steps of:

[0029] (a) receiving the rotation signals from the first gyroscope and the second gyroscope and the acceleration signals from the three accelerometers;

[0030] (b) establishing the desired azimuth direction;

[0031] (c) determining the pitch angle from the accelerometer signals;

[0032] (d) determining the rotation rate about an axis normal to the sensitive axes of the two gyroscopes by determining the time rate of change of the pitch angle;

[0033] (e) determining the effect of the earth rotation and local-level (l-frame) change of orientation to the angular changes monitored along the three axes and removing these two effects from the rotation rate signals;

[0034] (f) transforming the rotation signals from a BHA coordinate frame to the earth local-level coordinate frame;

[0035] (g) determining the pitch, roll and azimuth of the BHA;

[0036] (h) determining the time rate of change of the velocity of the BHA expressed in the l-frame and update the velocity components of the BHA; and

[0037] (i) updating the position components of the BHA.

[0038] In one embodiment, this method is used for an initial portion of the radical section of the wellbore. The initial portion of the radical section of the well may be that portion with an inclination less than about 45°, preferably less than about 30° and more preferably less than about 20°. After the initial portion, the method is altered by switching the axes orientation of the two gyroscopes. Once the wellbore reaches the horizontal section and rotary mode drilling is required, the surveying method may then switch to station-based surveying.

[0039] FIG. 1 is a block diagram illustrating the functionality of a prior art magnetic surveying system.

[0040] FIG. 2 is a sketch representing the inertial frame.

[0041] FIG. 3 is a sketch representing the Earth-fixed frame.

[0042] FIG. 4 is a sketch representing the local level frame.

[0043] FIG. 5 is a sketch showing the body frame axes arranged inside the bearing assembly of a BHA.

[0044] FIG. 6 shows the change of l-frame orientation along the Earth's surface.

[0045] FIG. 7 shows the arrangement of the fiber optic gyroscopes inside a bearing assembly.

[0046] FIG. 8 is a bottom view for the bearing assembly showing the flow of mud through the toroidal FOG.

[0047] FIG. 9 shows the initial arrangement of the FOGs and the accelerometers with the body frame axes inside the bearing assembly.

[0048] FIG. 10 shows the setup while surveying the radical section until 20° inclination.

[0049] FIG. 11 shows the setup after achieving 20° inclination.

[0050] FIG. 12 shows the new axis orientation after 20° inclination.

[0051] FIG. 13 shows the relationship between the pitch and the inclination angles before and after changing the axes orientation.

[0052] FIG. 14 is a block diagram describing the functionality of the surveying procedure until 20° inclination.

[0053] FIG. 14A is a block diagram describing the functionality of the surveying procedure until 20° inclination using only the toroidal FOG.

[0054] FIG. 15 is a block diagram describing the functionality of the surveying procedure after 20° inclination.

[0055] FIG. 16 shows the arrangement of the surveying instruments inside the bottom hole assembly while drilling the horizontal section of the well.

[0056] FIG. 17 is a block diagram showing the station-based surveying procedure.

[0057] FIG. 18 is a block diagram showing the technique of aided inertial navigation.

[0058] The present invention provides for a continuous measurement-while-drilling surveying method and apparatus. All terms not defined herein have their common art-recognized meanings.

[0059] Strapdown Intertial Navigation

[0060] There are two common INS categories [Titteron et al., 1997; Britting, 1971]. The first category is the gimbaled INS that maintains a stationary platform for the accelerometers by torque motors. These systems use the angular rate measurements of the gyroscopes attached to the platform to control the motors. The second category is the strapdown INS (SINS) that mathematically transforms the output of the accelerometers attached to the body into a locally-leveled coordinate system before performing integration. These systems use the output of the gyroscopes attached to the body to continuously convert from body coordinate to locally leveled coordinate.

[0061] Three reference frames are commonly used in the technology of inertial navigation [Schwarz, 1999; Salychev, 1998]. These are the inertial frame, the earth-fixed frame and the local level frame. The inertial frame (i-frame) is a non-rotating and non-accelerating frame. The origin of the inertial frame is at the center of mass of the Earth, its Z-axis is parallel to the spin axis of the Earth, its X-axis points toward the mean equinoctial colure and its Y-axis completes a right handed orthogonal frame (see FIG. 2). The Earth-fixed frame (e-frame) has its origin at the center of mass of the Earth with its Z-axis parallel to the spin axis of the Earth, X-axis pointing towards the mean meridian of Greenwich and Y-axis completing a right handed orthogonal frame (see FIG. 3). The local level frame (l-frame) has its origin at the origin of the sensor frame with its Z-axis orthogonal to the reference ellipsoid of the Earth and pointing upward, Y-axis pointing toward the geodetic north and X-axis completing a right handed orthogonal frame (see FIG. 4). When the inertial sensors (gyroscopes and accelerometers) are installed inside a moving object, they are oriented in three mutually orthogonal directions known as the body frame (b-frame). The axes of the b-frame point toward the forward (Y), transverse (X) and vertical (Z) directions (see FIG. 5).

[0062] The measurements of linear acceleration and angular velocities are taken in the b-frame. These measurements should be transformed into the corresponding reference frame within which they are processed to provide the position vector, the velocity vector and the attitude components. The transformation from the b-frame to the l-frame is performed with the rotation matrix R_b^l, which is expressed as a function of the attitude components (azimuth ψ, roll (tool face) φ and pitch (inclination) θ) as follows: $\begin{array}{cc}{R}_b^l=\left(\begin{array}{ccc}\cos \psi \cos \phi -\sin \mathrm{\psi sin}\theta \sin \phi & -\sin \mathrm{\psi cos}\theta & \cos \psi \sin \phi +\sin \mathrm{\psi sin}\theta \sin \phi \\ \sin \psi \cos \phi +\cos \mathrm{\psi sin}\theta \sin \phi & \cos \psi \cos \theta & \sin \psi \sin \phi -\cos \mathrm{\psi sin}\theta \cos \phi \\ -\cos \theta \sin \phi & \sin \theta & \cos \theta \cos \phi \end{array}\right)& \left(1\right)\end{array}$

[0063] The position of the BHA as described in the l-frame is expressed in terms of the curvilinear coordinates (latitude φ, longitude λ and altitude h). The velocity is alternatively expressed by three components along the East direction (V^e), north direction (Vⁿ) and vertical direction (V^u). Therefore, the time rate of change of the position components is related to the velocity components as follows [Schwarz et al., 1990]: $\begin{array}{cc}{\stackrel{.}{r}}^l=\left(\begin{array}{c}\stackrel{.}{\varphi }\\ \stackrel{.}{\lambda }\\ \stackrel{.}{h}\end{array}\right)=\left(\begin{array}{ccc}0& \frac{1}{M+h}& 0\\ \frac{1}{\left(N+h\right)\cos \varphi }& 0& 0\\ 0& 0& 1\end{array}\right)\left(\begin{array}{c}{V}^e\\ V^n\\ V^u\end{array}\right)=D^{-1}{V}^l& \left(2\right)\end{array}$

[0064] where

[0065] N :is the prime vertical radius of curvature of the Earth's surface (East-West);

[0066] M :is the corresponding meridian radius of curvature (North-South); $D^{-1}\text{:}\mathrm{is}a3\text{-}\mathrm{by}\text{-}3\mathrm{matrix}\mathrm{defined}\mathrm{as}\left(\begin{array}{ccc}0& \frac{1}{M+h}& 0\\ \frac{1}{\left(N+h\right)\cos \varphi }& 0& 0\\ 0& 0& 1\end{array}\right).$

[0067] The acceleration of the bearing assembly is measured in three mutually orthogonal directions in the b-frame f^b=(𝑓_x𝑓_y𝑓_z)^T. These measurements are transformed into the l-frame using the rotation matrix R_b^l(f^l=R_b^lf^b). However, the acceleration components expressed in the l-frame, f^l, can not directly give the velocity components of the bearing assembly due to three reasons. The first reason is due to the rotation rate of the Earth about its spin axis (ω^e=15°/hr) which is interpreted in the l-frame as ω_ie^l=(0 ω^ecosφ ω^esin φ)^T. The second reason is the change of orientation of the l-frame with respect to the Earth as the bearing assembly penetrates the downhole formation. The l-frame change of orientation along one of the meridians can be shown on FIG. 6. This change of orientation is due to the definition of the local North and vertical directions. The North direction (N) is tangent to the meridian all the time while the vertical direction (V^le) is normal to the Earth's surface. This effect is interpreted as rotation rate vector expressed as ω_el^l=(−{dot over (φ)} {dot over (λ)} cos φ {dot over (λ)} sin φ)^Twhich can be rewritten as ${\omega }_{el}^l={\left(\begin{array}{ccc}\frac{-V^n}{M+h}& \frac{V^e}{N+h}& \frac{V^e\tan \varphi }{N+h}\end{array}\right)}^T$

[0068] (see equation 2). The third reason is the gravity field of the Earth which is given as g^l=(0 0−γ)^T, where γ is obtained from the well-known normal gravity model given as follows [Schwarz et al., 1999]:

[0069] γ=a₁(1+a²sin²φ+a₃sin⁴φ)+(a₄+a₅sin²φ)h+a₆h²  (3)

[0070] where

[0071] Then, the time rate of change of the velocity components V^lof the bearing assembly while penetrating downhole can be expressed as follows [Schwarz et al., 1999]:

[0072] {dot over (V)}^l=R^l_b𝑓^b−(2Ω^l_ie=Ω^l_el)V^lg^l  (4)

[0073] where Ω_ie^land Ω_el^lare the skew-symmetric matrices corresponding to ω_ie^land ω_el^lrespectively and they are expressed as follows: $\begin{array}{cc}{\Omega }_{\mathrm{ie}}^l=\left(\begin{array}{ccc}0& -{\omega }^e\sin \varphi & {\omega }^e\cos \varphi \\ {\omega }^e\sin \varphi & 0& 0\\ -{\omega }^e\cos \varphi & 0& 0\end{array}\right)& \left(5\right)\\ {\Omega }_{\mathrm{el}}^l=\left(\begin{array}{ccc}0& \frac{-V^e\tan \varphi }{N+h}& \frac{V^e}{N+h}\\ \frac{V^e\tan \varphi }{N+h}& 0& \frac{V^n}{M+h}\\ \frac{-V^e}{N+h}& \frac{-V^n}{M+h}& 0\end{array}\right)& \left(6\right)\end{array}$

[0074] The rotation matrix R_b^lcan be obtained by solving the following differential equations [Schwarz et al., 1999]:

[0075] {dot over (R)}_b^l=R_b^lΩ_lb^b=R_b^l(Ω_ib^b−Ω_il^b)  (7)

[0076] where Ω_ib^bis the skew-symmetric matrix of the measurements of angular velocities provided by the gyroscopes. Ω_ib^bcorresponds to the angular velocity vector ω_ib^b=(ω_xω_yω_z)^Tand is given as: $\begin{array}{cc}{\Omega }_{\mathrm{ib}}^b=\left(\begin{array}{ccc}0& -{\omega }_z& {\omega }_y\\ {\omega }_z& 0& -{\omega }_x\\ -{\omega }_y& {\omega }_x& 0\end{array}\right)& \left(8\right)\end{array}$

[0077] Since the gyroscopes measure both the Earth rotation and the change in orientation of the l-frame in addition to the angular velocities of the bearing assembly, the term Ω_il^bin equation 7 is subtracted from Ω_ib^bto remove these two effects. The term Ω_il^bconsists of two parts. The first part is Ω_ie^bwhich accounts for the Earth rotation rate and the second part is Ω_el^bwhich accounts for the orientation change of the l-frame. Therefore, Ω_il^bcan be written as follows:

[0078] Ω_il^b=Ω_ie^b+Ω_el^b  (9)

[0079] Since ω_ie^b=R_l^bω_ie^land ω_el^bω_el^l, then ω_il^b=R_l^b(ω_ie^l+ω_el^l) and can be given as follows: $\begin{array}{cc}{\omega }_{\mathrm{il}}^b=R_l^b\left[\left(\begin{array}{c}0\\ {\omega }^e\cos \varphi \\ {\omega }^e\sin \varphi \end{array}\right)+\left(\begin{array}{c}\frac{-V^n}{M+h}\\ \frac{V^e}{N+h}\\ \frac{V^e\tan \varphi }{N+h}\end{array}\right)\right]=R_l^b\left(\begin{array}{c}\frac{-V^n}{M+h}\\ \frac{V^e}{N+h}+{\omega }^e\cos \varphi \\ \frac{V^e\tan \varphi }{N+h}+{\omega }^e\sin \varphi \end{array}\right)& \left(10\right)\end{array}$

[0080] Consequently, the corresponding skew-symmetric matrix Ω_il^bcan be determined.

[0081] When equations 2, 4 and 9 are combined together, they form what is known as the mechanization equations in the l-frame and they are usually given together as follows [Schwarz et al., 1990; Schwarz et al., 1999]: $\begin{array}{cc}\left(\begin{array}{c}{\stackrel{.}{r}}^l\\ {\stackrel{.}{V}}^l\\ {\stackrel{.}{R}}_b^l\end{array}\right)=\left(\begin{array}{c}{D}^{-1}{V}^l\\ R_b^l{f}^b-\left(2{\Omega }_{\mathrm{ie}}^l+{\Omega }_{\mathrm{el}}^l\right)V^l+g^l\\ R_b^l\left({\Omega }_{\mathrm{ib}}^b-{\Omega }_{\mathrm{il}}^b\right)\end{array}\right)& \left(11\right)\end{array}$

[0082] The inputs to these mechanization equations are the gyroscope and accelerometer measurements while the outputs are the curvilinear coordinates, three velocity components and three attitude components. The l-frame is chosen to describe the drill bit motion for three reasons. The first is that the axes of the l-frame are aligned towards the local East, North and vertical (normal to the Earth's surface and pointing upward). Using this orientation, the pitch, the roll and the azimuth angles describing the drill bit attitude (inclination, tool face and azimuth) can be obtained directly as outputs of the mechanization equations solved in the l-frame. The second reason for choosing the l-frame is that the geodetic coordinates expressing the drill bit latitude, longitude and depth can be obtained directly by solving the mechanization equations described in this frame. Finally, the computational errors in the navigation parameters on the horizontal plane (North-East plane) are bounded due to Schuller's effect. This effect stipulates that the inertial system errors of the horizontal plane components are coupled together producing what is called a Schuller loop [Mohammed, 1999; Titterton et al., 1997]. Therefore these errors oscillate with frequency called the Schuller frequency (1/5000 Hz).

[0083] Apparatus

[0084] At present, the FOG size determines the overall size of the inertial measurement unit. As mentioned earlier, present technological level of FOG-based inertial measurement units do not permit direct installation of 3-FOGs inside the BHA. Therefore, the present invention comprises a novel structure of inertial measurement unit with two specially designed FOGs in order to accommodate space and precision requirements. These two FOGs monitor two rotations rate signals along the Y and Z directions of the b-frame axes. The third rotation rate signal may be generated using the time rate of change of the inclination angle calculated with the three-axis accelerometer measurements. The two FOGs and the three-axis accelerometers may be mounted anywhere inside the drill pipe and the technique provided in the following sections can process their measurements to determine the surveying data.

[0085] In one embodiment, the MWD inertial surveying sensors are fitted inside the bearing assembly close behind the drill bit as is shown in FIG. 7. The system comprises a first toroidal FOG (10), a second FOG (12), three accelerometers (14, 16, 18), a microprocessor (20), wireless communication transmitter (22) and a power source such as a battery. The microprocessor (20) receives the output signals from the two FOGs (10, 12) and the three accelerometers (14, 16, 18) and determines the inertial navigation components as described herein. These components are then transmitted to the surface for review by a drilling engineer. Alternatively, some or all of the processing may take place at the surface if the output signals of the inertial sensors are transmitted to the surface. Standard transmission or telemetry apparatus and methods may be used, as is well known in the art of MWD.

[0086] The toroidal FOG (10) is installed with its sensitive axis along the tool spin axis. It's shape allows the flow of mud through the drill pipe as is clearly shown in FIG. 8. The dimensions of the toroidal FOG (10) are dictated by the dimensions of the bearing assembly. The second FOG (12) is installed so that its sensitive axis is normal to the tool spin axis and towards the forward direction (Y-axis) Since its sensitive axis is normal to the tool spin axis, this FOG is denoted as the normal FOG (12). Both FOGs should preferably have drift rate less than 1 deg./hr and angle random walk of less than 0.1 deg./hr/√{square root over (Hz)} in 100 Hz bandwidth. The bias drift is defined as the deviation in the measured rotation rate at constant temperature. The angle random walk is the broadband random noise component at the FOG output that results either from the shot noise or the thermal noise in the photodetector.

[0087] The three-axis accelerometers (14, 16, 18) are arranged in three mutually orthogonal directions with one of them having its sensitive axis parallel to the tool spin axis.

[0088] Standard wireless communication between the system inside the bearing assembly and the surface MWD processor is provided [Skillingstad, 2000]. A power source such as a battery package provides the necessary power supply to the two FOGs, the accelerometers, the microprocessor and the wireless communication system.

[0089] Once the vertical section of the well is established, the horizontal drilling process involves three main tasks:

[0090] 1. Establishing the desired azimuth direction while the drill pipe is still in the vertical direction.

[0091] 2. Drilling the radical section of the well using steering mode of operation.

[0092] 3. Drilling the horizontal section of the well with using rotary mode of operation.

[0093] The Establishment of the Desired Azimuth Direction

[0094] In order to establish accurately the azimuth direction, precise value of the initial azimuth should be known. This initial azimuth value is determined when the whole system is completely stationary. Therefore, the two FOGs monitor the Earth rotation, while the accelerometers monitor the Earth gravity. The Earth rotation vector monitored by the gyroscopes, ω_ie^b, and the Earth rotation vector expressed in the l-frame, ω_ie^l, are related through the rotation matrix R_l^b=(R_b^l)^Twhich transforms the l-frame to the b-frame. This relationship can be written as follows: $\begin{array}{cc}\begin{array}{c}{\omega }_{\mathrm{ie}}^b={\left(\begin{array}{ccc}{\omega }_x& {\omega }_y& {\omega }_z\end{array}\right)}^T=R_l^b{\omega }_{\mathrm{ie}}^l\\ =\left(\begin{array}{ccc}\cos \psi \cos \phi -\sin \mathrm{\psi sin}\theta \sin \phi & \sin \psi \cos \phi +\cos \mathrm{\psi sin}\theta \sin \phi & -\cos \theta \sin \phi \\ -\sin \psi \cos \phi & \cos \psi \cos \theta & \sin \theta \\ \cos \psi \sin \phi +\sin \mathrm{\psi sin}\theta \cos \phi & \sin \psi \sin \phi -\cos \mathrm{\psi sin}\theta \cos \phi & \cos \theta \cos \phi \end{array}\right)\left(\begin{array}{c}0\\ {\omega }^e\cos \varphi \\ {\omega }^e\sin \varphi \end{array}\right)\end{array}& \left(12\right)\end{array}$

[0095] The body frame axes are aligned as shown on FIG. 9 with the normal FOG having its sensitive axis along the Y-axis and the toroidal FOG having its sensitive axis along the Z-axis (the tool spin axis). Using equation 12, the rotation rate measurement provided by the normal FOG, ω_y, can be given as function of the Earth rotation rate ω^eas follows:

[0096] ω_y=(cos ψ cos θ)ω^ecos φ+(sin θ)ω^esin φ  (13)

[0097] Consequently, the azimuth can be determined as follows: $\begin{array}{cc}\cos \psi =\frac{\frac{{\omega }_y}{{\omega }^e\cos \varphi }-\sin \theta \tan \varphi }{\cos \theta }& \left(14\right)\end{array}$

[0098] Since the bottom hole assembly is almost in the vertical direction, the inclination angle (i.e. the pitch angle, θ) is very close to zero. Therefore, we can simplify equation 14 by considering sin θ≈0 and cos θ≈1. The initial azimuth can then be given as: $\begin{array}{cc}\cos \psi =\frac{{\omega }_y}{{\omega }^e\cos \varphi }& \left(15\right)\end{array}$

[0099] It should be noted that equation 15 has a singularity when the latitude angle φ equals to ±π/2. This situation corresponds to drilling at the Earth poles that can not exist.

[0100] Therefore the initial azimuth can be determined using equation 15 by substituting the angular velocity output of the normal FOG, ω_y, after removing the offset and the constant drift values. The residual errors in the measurements provided by this FOG, δω_y, affect the calculations of cos ψ and consequently the accuracy of the initial azimuth. The error in cos ψ due to these errors δω_ycan be obtained by differentiating equation 15: $\begin{array}{cc}\delta \cos \psi =\frac{\delta {\omega }_y}{{\omega }^e\cos \varphi }& \left(16\right)\end{array}$

[0101] It is crucial to remove the error δ cos ψ from the cos ψ term calculated from equation 15 before proceeding to calculate the corresponding value of the azimuth ψ. This is necessary because of the inversely proportional relationship between the error in the azimuth δψ and the azimuth angle itself, ψ, which can be obtained from equation 16 as follows: $\begin{array}{cc}\delta \psi =\frac{\left(\frac{\delta {\omega }_y}{{\omega }^e\cos \varphi }\right)}{\sin \psi }& \left(17\right)\end{array}$

[0102] Since sin ψ≦1, the term (1/sin ψ) in equation 17 tends to magnify the error δω_y. When ψ is close to π/2, i.e. approaching the East direction, the term 1/sin ψ becomes closer to unity and has minor effect on δψ. On the other hand, when ψ is close to 0, i.e. approaching the North direction, the term 1/sin ψ magnifies the random error δψ. In order to avoid this problem, the error term δ cos ψ is determined first using equation 16 and removed from the values obtained for cos ψ using equation 15. Once an accurate estimation of cos ψ is determined, the azimuth can be computed independently from the drilling direction.

[0103] Once the initial azimuth of the bottom hole assembly is determined, the drill pipe can be rotated to achieve the desired azimuth direction. Although the normal FOG is responsible for the determination of the initial value of the azimuth, the toroidal FOG is responsible for the establishment of the desired azimuth direction.

[0104] With the BHA oriented along the vertical direction, the drill pipe is rotated along the tool spin axis until the desired azimuth direction is established. Rotation along the tool spin axis is preferably the only motion, either angular or translational, performed during this operation. Under these circumstances, the toroidal FOG monitors the rotation rate along the tool spin axis in addition to the component of Earth rotation along the vertical direction (ω^esin φ). Therefore, the time rate of change of the azimuth, {dot over (ψ)}, can be written as follows:

[0105] {dot over (ψ)}=ω_z−ω^esin φ  (18)

[0106] where ω_zis the measurement delivered by the toroidal FOG.

[0107] Equation 18 is solved numerically in real-time to provide continuous monitoring of the azimuth angle while rotating the drill pipe. This first-order differential equation can be solved numerically using the Euler method [Yakowitz et al, 1989] after substituting {dot over (ψ)} by (ψ_k+1−ψ_k)/(Δt), where ψ_k+1is the azimuth value at time t_k+1, ψ_kis the azimuth value at time t_kand Δt=t_k+1−t_k. Then, the azimuth can be continuously monitored using the following equation:

[0108] ψ_k+1=ψ_k+(ω_z−ω^esin φ)Δt  (19)

[0109] As the drill pipe continuously rotates about its tool spin axis, the azimuth angle is continuously determined from equation 19. Once the desired azimuth direction is established, the rotation of the drill pipe is stopped and building the radical section of the well is initiated.

[0110] Surveying the Radical Section of the Well

[0111] During the building of the radical section of the well, the bottom hole assembly starts the drilling process from complete vertical direction and continues toward a complete horizontal direction. At the beginning, while the BHA is still in the vertical direction, the body frame axes are chosen as shown on FIG. 9. The X, Y, and Z axes point toward transverse, forward and vertical directions respectively. Therefore, the Z-axis points toward the tool spin axis. With this choice of axes orientation, we assume that the motion of the BHA is along the Y-axis, and sliding along the Z-axis (FIG. 10). Hence, the toroidal FOG monitors the changes in the azimuth angle ψ while the normal FOG monitors the changes in the roll angle φ. However, as the wellbore inclination increases, the Z-axis approaches the horizontal and the Y-axis approaches the vertical. The Y-axis is no longer in the forward direction and the Z-axis is no longer in the vertical direction. In addition, the toroidal FOG is no longer in the horizontal plane to successfully monitor the azimuth. Therefore, at a predetermined inclination such as 20°, a new orientation of the axes may be set to keep Y along the forward direction and Z along the vertical direction. This can be achieved by virtually choosing new axes system X¹, Y¹, and Z¹, as is illustrated in FIG. 12. X¹points toward the transverse direction exactly like X while Y¹points toward the forward direction (i.e. along the tool spin axis) and Z¹points toward the vertical direction. This axes orientation is maintained until the end of the drilling process without any physical change in the mounting of the two FOGs.

[0112] In a preferred embodiment, there are two surveying methodologies for the radical section of the well. The first is applied during an initial portion of the radical section. As used herein, an “initial portion” is that portion of the well where the inclination is less than about 45°, preferably less than about 30° and more preferably less than about 20° inclination. The body frame is defined by the X, Y and Z axes as is shown in FIG. 10. During this initial portion of the well, the toroidal FOG (with its sensitive axis along the Z-axis) is responsible for monitoring the changes in the azimuth angle while the normal FOG (with its sensitive axis along the Y-axis) is responsible for monitoring the changes in the roll angle. Since the toroidal FOG monitors the rotation along the tool spin axis, the tool face angle is directly obtained from the toroidal FOG output.

[0113] After the initial portion, the tasks of the two FOGs are switched. In addition, the roll and the tool face angles become the same and are monitored by the toroidal FOG while the normal FOG monitors the changes in the azimuth angle.

[0114] It should also be noted that the relationship between the pitch angle and the inclination angle is not the same for the two cases during and after the initial portion of the well. Before 20° inclination the pitch angle is equal to the inclination angle while after 20° inclination there is π/2 radians phase shift between them. This is because the pitch angle is defined as the angle between the forward direction (Y-axis) of the moving body and the horizontal plane (FIG. 13). Since the definition of the Y-axis is changed at 20° inclination, the definition of the pitch angle is changed accordingly and the π/2 radians phase shift is created between the inclination and the pitch angle.

[0115] The pitch angle is determined from the three-axis accelerometers. The determination of the inclination angle is based on the assumption that the three-axis accelerometers monitor the gravity components in three mutually orthogonal directions and their outputs can be integrated together to determine the pitch. This assumption is valid as long as the penetration rate through the downhole formation is relatively small. The average angle-building rate is about 4°/100 ft. for long-radius application, 12°/100 ft. for medium-radius applications and 4°/ft. for short-radius applications. In each of these three categories of horizontal drilling, the above assumption can be applied and the pitch can be determined from the measurements of the accelerometers.

[0116] Surveying the Radical Section Before 20° Inclination

[0117] The surveying technique explained in this sub section provides continuous surveying based on the measurements delivered by the three-axis accelerometers, the toroidal FOG and the normal FOG. The toroidal FOG determines the rotation rate along the Z-axis, ω_z, while the normal FOG determines the rotation rate along the Y-axis, ω_y. The rotation rate along the X-axis, ω_x, can be determined from the time rate of change of the pitch angle. Therefore, the process starts with receiving the accelerometer measurements that can be combined to determine the pitch, which, until 20° inclination, is the same as the inclination angle. Based on the assumption of very small penetration rate through the downhole formation, we assumed that the three accelerometers are basically affected by the gravity field of the Earth only. Therefore, the accelerometer measurements taken at the body frame, g^bare related to the gravity field of the Earth expressed in the l-frame, g^l, through the rotation matrix R_l^bas follows: $\begin{array}{cc}{g}^b=\left(\begin{array}{c}{f}_x\\ f_y\\ f_z\end{array}\right)=R_l^b{g}^l=\left(\begin{array}{ccc}\cos \psi \cos \phi -\sin \mathrm{\psi sin}\theta \sin \phi & \sin \psi \cos \phi +\cos \mathrm{\psi sin}\theta \sin \phi & -\cos \theta \sin \phi \\ -\sin \psi \cos \phi & \cos \psi \cos \theta & \sin \theta \\ \cos \psi \sin \phi +\sin \mathrm{\psi sin}\theta \cos \phi & \sin \psi \sin \phi -\cos \mathrm{\psi sin}\theta \cos \phi & \cos \theta \cos \phi \end{array}\right)\left(\begin{array}{c}0\\ 0\\ -\gamma \end{array}\right)& \left(20\right)\end{array}$

[0118] where 𝑓_x, 𝑓_yand 𝑓_zare the measurements provided by the three mutually orthogonal accelerometers.

[0119] γ is the value of gravity obtained from the normal gravity model given in equation 3.

[0120] From equation 20, the accelerometer measurements 𝑓_x, 𝑓_yand 𝑓_zcan be given in terms of γ as follows:

[0121] 𝑓_x=(cos θ sin φ)γ  (21)

[0122] 𝑓_y=−(sin θ)γ  (22)

[0123] 𝑓_z=−(cos θ cos φ)γ  (23)

[0124] where θ is the pitch (inclination) angle and φ is the roll angle.

[0125] It is clear that the pitch angle can be obtained directly from equation 22 as long as γ is determined from the normal gravity model (equation 3). The determination of γ depends on the available information for the latitude φ and the altitude h. However, since the normal gravity model is an approximate model, and in order to prevent the propagation of errors, equations 21 through 23 are manipulated to remove the gravity term and to directly give the pitch angle. From equations 21 and 23, we can write

[0126] √{square root over (𝑓_x²+𝑓_z²)}=(cos θ)γ  (24)

[0127] Then the pitch angle can be obtained by dividing equation 22 by equation 24. $\begin{array}{cc}\theta =\arctan \left(\sqrt{\frac{f_y^2}{f_x^2+f_z^2}}\right)& \left(25\right)\end{array}$

[0128] Equation 25 gives the pitch angle directly in a way that is completely independent of the gravity field of the Earth. Then the rotation rate along the X-axis, ω_x, is determined using the time rate of change of the pitch angle as follows: $\begin{array}{cc}{\omega }_x=\frac{{\theta }_{k+1}-{\theta }_k}{t_{k+1}-t_k}& \left(26\right)\end{array}$

[0129] where θ_kand θ_k+1are the values of the pitch angle at time t_kand t_k+1respectively.

[0130] The surveying procedure delivers the navigation data (the position and the orientation) for the bottom hole assembly by solving the set of first order differential equations given in equation 11. The first step in solving these equations is the parameterization of the rotation matrix R_b^l[Schwarz et al., 1999; Salychev, 1998; Titterton et al., 1997]. One method for this purpose is the quaternion approach. According to Euler's Theorem [Salychev, 1998], the rotation of a rigid body (represented by the b-frame) with respect to a reference frame (represented by the l-frame) can be expressed in terms of the rotation angle Θ about a fixed axis and the direction cosine of the rotation axis that defines the rotation direction. Thus, quaternion parameters Q=(q₁q₂q₃q₄)^Tare introduced to describe the rotation of the b-frame with respect to the l-frame and they are expressed as follow: $\begin{array}{cc}Q=\left(\begin{array}{c}{q}_1\\ q_2\\ q_3\\ q_4\end{array}\right)=\left(\begin{array}{c}\left({\Theta }_x/\Theta \right)\sin \left(\Theta /2\right)\\ \left({\Theta }_y/\Theta \right)\sin \left(\Theta /2\right)\\ \left({\Theta }_z/\Theta \right)\sin \left(\Theta /2\right)\\ \cos \left(\Theta /2\right)\end{array}\right)& \left(27\right)\end{array}$

[0131] where Θ={square root over (Θ_x²+Θ_y²+Θ_z²)} is the rotation angle and $\frac{{\Theta }_x}{\Theta },\frac{{\Theta }_y}{\Theta },\frac{{\Theta }_z}{\Theta }$

[0132] are the three direction cosines of the rotation axis with respect to the l-frame.

[0133] The definition of the quaternion parameters described in equation (27) implies that the four quaternion components are not independent because q₁²+q₂²+q₃²+q₄²=1. This means that only three independent quaternion components are sufficient to describe the rigid body rotation. The time rate of change of the quaternion is described by the following first-order differential equation: $\begin{array}{cc}\stackrel{.}{Q}=\frac{1}{2}\Omega \left(\omega \right)Q& \left(28\right)\end{array}$

[0134] where Ω(ω) is a skew-symmetric matrix given as follows: $\begin{array}{cc}\Omega \left(\omega \right)=\left(\begin{array}{cccc}0& {\omega }_z& -{\omega }_y& {\omega }_x\\ -{\omega }_z& 0& {\omega }_x& {\omega }_y\\ {\omega }_y& -{\omega }_x& 0& {\omega }_z\\ -{\omega }_x& -{\omega }_y& -{\omega }_z& 0\end{array}\right)& \left(29\right)\end{array}$

[0135] where ω_x, ω_y, ω_zare the angular velocities of body rotation which are determined by equation 26 for ω_xand are monitored by the normal FOG and the toroidal FOG for ω_yand ω_zrespectively, after compensating for the Earth rotation and the l-frame change of orientation.

[0136] To solve the set of first order differential equations given in equation 28, Euler method can be used to determine the quaternion parameters Q_k+1at time t_k+1based on the values of the quaternion parameters Q_kat time t_kas follows: $\begin{array}{cc}{Q}_{k+1}=Q_k+\left(\frac{1}{2}\Omega \left({\omega }_k\right)Q_k\right)\Delta t& \left(30\right)\end{array}$

[0137] where Δt=t_k+1−t_k.

[0138] Once the quaternion parameters are determined at certain time, the rotation matrix R_b^lcan be obtained using the following direct relationship: $\begin{array}{cc}\begin{array}{c}{R}_b^l=\left(\begin{array}{ccc}{R}_{11}& R_{12}& R_{13}\\ R_{21}& R_{22}& R_{23}\\ R_{31}& R_{32}& R_{33}\end{array}\right)\\ =\left(\begin{array}{ccc}{q}_1^2-q_2^2-q_3^2+q_4^2& 2\left(q_1{q}_2-q_3{q}_4\right)& 2\left(q_1{q}_3+q_2{q}_4\right)\\ 2\left(q_1{q}_2+q_3{q}_4\right)& -q_1^2+q_2^2-q_3^2+q_4^2& 2\left(q_2{q}_3-q_1{q}_4\right)\\ 2\left(q_1{q}_3-q_2{q}_4\right)& 2\left(q_2{q}_3+q_1{q}_4\right)& -q_1^2-q_2^2+q_3^2+q_4^2\end{array}\right)\end{array}& \left(31\right)\end{array}$

[0139] On the other hand, the quaternion parameters can be obtained from the rotation matrix R_b^las follows: $\begin{array}{cc}\left(\begin{array}{c}{q}_1\\ q_2\\ q_3\\ q_4\end{array}\right)=\left(\begin{array}{c}0.25\left(R_{32}-R_{23}\right)/q_4\\ 0.25\left(R_{13}-R_{31}\right)/q_4\\ 0.25\left(R_{21}-R_{12}\right)/q_4\\ 0.5\sqrt{1+R_{11}+R_{22}+R_{33}}\end{array}\right)& \left(32\right)\end{array}$

[0140] The quaternion parameters are introduced for the parameterization of the rotation matrix R_b^lfor three reasons. The first reason is that only four differential equations are solved numerically instead of six differential equations if we manipulate the rotation matrix R_b^ldirectly. The second reason is that the quaternion solution avoids the singularity problem that might exist with some other solution methods. The third reason is the computational simplicity introduced by the quaternion.

[0141] One embodiment of the method may be described as follows:

[0142] 1. Receive the measurements from the two FOGs and the three accelerometers. Some modem FOGs deliver angular increments at their output instead of rotation rate. Similarly, the accelerometers deliver velocity increments instead of specific force (acceleration) measurements. The angular changes provided by the two FOGs should be compensated for the drift rate d as follows:

[0143] Δθ_y=Δ{circumflex over (θ)}_y−d_yΔt  (33)

[0144] Δθ_z=Δ{circumflex over (θ)}_z−d_zΔt

[0145] where

[0146] Δ{circumflex over (θ)}_yis the raw measurement of the angular increments provided by the normal FOG;

[0147] Δθ_yis the measurement of the angular increments provided by the normal FOG after compensating for the drift rate d_y;

[0148] Δ{circumflex over (θ)}_zis the raw measurement of the angular increments provided by the toroidal FOG;

[0149] Δθ_zis the measurement of the angular increments provided by the toroidal FOG after compensating for the drift rate d_z;

[0150] Δt is the difference between the current time sample t_k+1and the previous time sample t_k.

[0151] Similarly, the velocity increments provided by the accelerometers should be compensated for the bias terms b as follows:

[0152] Δν_x=Δ{circumflex over (ν)}_x−b_xΔt

[0153] Δν_y=Δ{circumflex over (ν)}_y−b_yΔt

[0154] Δν_z=Δ{circumflex over (ν)}_z−b_zΔt  (34)

[0155] where

[0156] Δ{circumflex over (ν)}_xis the raw measurement of the velocity increments provided by the accelerometer installed along the X-axis;

[0157] Δv_xis the measurement of the velocity increments provided by the accelerometer installed along the X-axis after compensating for the bias b_x;

[0158] Δ{circumflex over (ν)}_yis the raw measurement of the velocity increments provided by the accelerometer installed along the Y-axis;

[0159] Δv_yis the measurement of the velocity increments provided by the accelerometer installed along the Y-axis after compensating for the bias b_y;

[0160] Δ{circumflex over (ν)}_zis the raw measurement of the velocity increments provided by the accelerometer installed along the Z-axis;

[0161] Δv_zis the measurement of the velocity increments provided by the accelerometer installed along the Z-axis after compensating for the bias b_z;

[0162] Δt is the difference between the current time sample t_k+1and the previous time sample t_k.

[0163] 2. Determine the pitch angle from the accelerometer measurements as given in equation 25. This equation is slightly modified to suit the velocity increment output Δν from the accelerometers. This modification is based on the direct relationship between the specific force measurement f, representing the acceleration, and the velocity increments Δν given as 𝑓=Δν/Δt. Thus, the pitch angle, θ, is determined as follows: $\begin{array}{cc}\theta =\arctan \left(\sqrt{\frac{\Delta v_y^2}{\Delta v_x^2+\Delta v_z^2}}\right)& \left(35\right)\end{array}$

[0164] Then, the angular changes along X-direction, Δθ_x, can be obtained from equation 26 as follows

[0165] Δθ_x=ω_xΔt=[θ(t_k+1)−θ(t_k)]  (36)

[0166] 3. Determine the effect of the Earth rotation and the l-frame change of orientation in the angular changes monitored along the three axes. The term responsible for the Earth rotation and the l-frame change of orientation is ω_il^band is given as follows (see equation 10): $\begin{array}{cc}{\omega }_{\mathrm{il}}^b=R_l^b\left(t_k\right)\left(\begin{array}{c}\frac{-V^n\left(t_k\right)}{M+h}\\ \frac{V^e\left(t_k\right)}{N+h}+{\omega }^e\cos \varphi \\ \frac{V^e\left(t_k\right)\tan \varphi }{N+h}+{\omega }^e\sin \varphi \end{array}\right)& \left(37\right)\end{array}$

[0167] Then the angular changes corresponding to ω_il^bcan be determined as:

[0168] θ_il^b=ω_il^bΔt  (38)

[0169] Consequently, the vector of measured angular increments θ_ib^b=(Δθ_xΔθ_yΔθ_z)^Tis compensated for θ_il^bto determine the actual bottom hole assembly angular changes θ_lb^bas follows:

[0170] θ_lb^b=θ_ib^b−θ_il^b=(δθ_xδθ_yδθ_z)^T  (39)

[0171] 4. Update the quaternion. The time rate of change of the quaternion as given in equation 28 corresponds to the time rate of change of the rotation matrix R_b^las given in equation 7. The quaternion update follows equation 30. The time interval Δt is multiplied by the elements of the skew-symmetric matrix of angular velocities Ω, so that ω_xΔt=δθ_x, ω_yΔt=δθ_yand ω_zΔt=δθ_z. The quaternion update equation is then given as follows: $\begin{array}{cc}\left(\begin{array}{c}{q}_1\left(t_{k+1}\right)\\ q_2\left(t_{k+1}\right)\\ q_3\left(t_{k+1}\right)\\ q_4\left(t_{k+1}\right)\end{array}\right)=\left(\begin{array}{c}{q}_1\left(t_k\right)\\ q_2\left(t_k\right)\\ q_3\left(t_k\right)\\ q_4\left(t_k\right)\end{array}\right)+\frac{1}{2}\left(\begin{array}{cccc}0& \delta {\theta }_z& -\delta {\theta }_y& \delta {\theta }_x\\ -\delta {\theta }_z& 0& \delta {\theta }_x& \delta {\theta }_y\\ \delta {\theta }_y& -\delta {\theta }_x& 0& \delta {\theta }_z\\ -\delta {\theta }_x& -\delta {\theta }_y& -\delta {\theta }_z& 0\end{array}\right)\left(\begin{array}{c}{q}_1\left(t_k\right)\\ q_2\left(t_k\right)\\ q_3\left(t_k\right)\\ q_4\left(t_k\right)\end{array}\right)& \left(40\right)\end{array}$

[0172] where the initial value of the quaternion Q(t₀)=(q₁(t₀) q₂(t₀) q₃(t₀) q₄(t₀))^Tis chosen to be equal to (0 0 0 1)^T.

[0173] Consequently, the rotation matrix R_b^lcan be determined from its direct relationship to the quaternion Q=(q₁q₂q₃q₄)^Tas given by equation 31. From the R_b^ldefinition given in equation 1, the pitch (inclination) is determined as θ=arc sin(R₃₂), the roll is determined as φ=arc tan(−R₃₁/R₃₃) and the azimuth (tool face) is determined as ψ=arc tan(−R₁₂/R₂₂)

[0174] 5. Determine the changes in velocity of the bottom hole assembly using equation 4 that describes the time rate of change of the velocity of the bottom hole assembly expressed in the l-frame. This equation can be rewritten as follows: $\begin{array}{cc}\begin{array}{c}\frac{\Delta V^l}{\Delta t}=R_b^l{f}^b-\left(2{\Omega }_{\mathrm{ie}}^l+{\Omega }_{\mathrm{el}}^l\right)V^l+g^l\\ \mathrm{or}\\ \Delta V^l=R_b^l{f}^b\Delta t-\left(2{\Omega }_{\mathrm{ie}}^l+{\Omega }_{\mathrm{el}}^l\right)V^l\Delta t+g^l\Delta t\end{array}& \left(41\right)\end{array}$

[0175] Since the specific force measurement f^b=(𝑓_x𝑓_y𝑓_z)^Tis related to the velocity increment measurements Δν^b=(Δν_xΔν_yΔν_z)^Twith f^b=Δν^b/Δt, equation 41 can be given as:

[0176] ΔV^l(t_k+1)=R_b^lΔν^b(t_k+1)−(2Ω_ie^l+Ω_el^l)V^l(t_k)Δt+g^lΔt  (42)

[0177] where Ω_ie^land Ω_el^lare given in equations 5 and 6 respectively, while g^lis the gravity field vector and is defined as shown in equation 3.

[0178] The first term on the right-hand side of equation 42 is the measured velocity increments transformed to the l-frame. The second term is the Coriolis correction that compensates for the Earth rotation and the l-frame change of orientation. The third tem is the gravity correction.

[0179] The velocity components at the current time epoch t_k+1are computed using the modified Euler formula as follows: $\begin{array}{cc}{V}^l\left(t_{k+1}\right)=V^l\left(t_k\right)+\frac{1}{2}\left(\Delta V^l\left(t_k\right)+\Delta V^l\left(t_{k+1}\right)\right)& \left(43\right)\end{array}$

[0180] where V^l=(V^eVⁿV^u) are the velocity components as described in the l-frame.

[0181] 6. Determine the altitude of the bottom hole assembly, h, (i.e. the true vertical depth TVD) using the vertical component of the velocity vector, V^u. Since {dot over (h)}=V^uaccording to equation 2, the TVD can be computed at time instant t_k+1using the modified Euler formula as follows: $\begin{array}{cc}h\left(t_{k+1}\right)=h\left(t_k\right)+\frac{1}{2}\left(V^u\left(t_{k+1}\right)+V^u\left(t_k\right)\right)\Delta t& \left(44\right)\end{array}$

[0182] Once the TVD is determined, the latitude φ, and the longitude λ, can be computed. According to equation 2, φ and λ can be computed at time instant (t_k+1) using the modified Euler formula as follows: $\begin{array}{cc}\varphi \left(t_{k+1}\right)=\varphi \left(t_k\right)+\frac{1}{2}\frac{\left(V^n\left(t_{k+1}\right)+V^n\left(t_k\right)\right)}{M+h}\Delta t& \left(45\right)\\ \lambda \left(t_{k+1}\right)=\lambda \left(t_k\right)+\frac{1}{2}\frac{\left(V^e\left(t_{k+1}\right)+V^e\left(t_k\right)\right)}{\left(N+h\right)\cos \varphi }\Delta t& \left(46\right)\end{array}$

[0183] It should be noted that the computation of the longitude, λ, will suffer from singularities if this procedure is performed closer to the poles where the latitude φ becomes equal to π/2.

[0184] The six steps described above are repeated for each new set of measurements received from the two FOGs and the three accelerometers until 20° inclination. A block diagram describing the functionality of this procedure is shown on FIG. 14.

[0185] In an alternative embodiment of the invention, surveying may be accomplished with only the toroidal FOG and the three accelerometers until 20°inclination. A block diagram describing the functionality of this procedure is shown in FIG. 14A. In this case, the toroidal FOG may monitor both azimuth and roll angles because of the near verticality of the BHA. As the pitch angle increases, the unity between azimuth and roll decreases and the ability of the toroidal FOG to monitor roll angle diminishes.

[0186] Surveying the Radical Section After 20° Inclination

[0187] As shown earlier on FIG. 13, change of axes orientation may take place at 20° inclination in one embodiment of the invention. The Y¹-axis is now oriented along the tool spin axis while Z¹-axis is oriented along the orthogonal direction. The X¹-axis remains collinear to the old X-axis. The initial values for the navigation parameters (position components, velocity components and attitude components) for surveying this part of the radical section of the well can be obtained from the last values determined just before changing the axes orientation. It should be reiterated that the location of the two FOGs and the accelerometers would not physically change due to the change of axes orientation. Once the procedure for the first part of the radical section of the well (before 20° inclination) indicates an inclination of 20°, the output of the toroidal FOG is interpreted as Δ{circumflex over (θ)}_y, since its sensitive axis is reoriented along Y¹-axis. Similarly, the output of the normal FOG is interpreted as Δ{circumflex over (θ)}_z, since its sensitive axis is reoriented along the Z¹-axis (see FIG. 12). In addition, the output of the accelerometer installed along the tool spin axis (Y¹-axis) is interpreted as Δ{circumflex over (ν)}_yinstead of Δ{circumflex over (ν)}_zwhile the output of the accelerometer installed along the orthogonal direction (Z¹-axis) is interpreted as Δ{circumflex over (ν)}_zinstead of Δ{circumflex over (ν)}_y. The third accelerometer installed along the X¹-axis has its output interpreted as Δ{circumflex over (ν)}_xthe same as in the first part of the drilling of the radical section of the well.

[0188] The navigation procedure for this part of the radical section of the well is exactly the same as the procedure described in the previous section except for three changes. The first change is the different interpretation of the measurements provided by the two FOGs and the three accelerometers. The second change is that the tool face angle becomes no longer equivalent to the azimuth. Instead, it becomes equivalent to the roll. The third change is that the inclination is no longer equivalent to the pitch (θ) and is given as π/2−θ. The navigation procedure for surveying after 20° inclination is shown schematically in the block diagram shown as FIG. 15.

[0189] Surveying the Horizontal Section of the Well

[0190] Surveying the horizontal section of the well is different from the radical section due to the rotary mode of drilling during which the whole drill pipe starts to continuously rotate about its spin axis. Therefore, continuous surveying can no longer be performed and station-based surveying procedure should be prepared. The station-based surveying is based on stopping the drilling operation at certain survey stations. Since the surveying setup is completely stationary at each station, the two FOGs are affected by the Earth rotation rate and the three accelerometers are affected by the Earth gravity field. The arrangement of the surveying instruments inside the bottom hole assembly is shown on FIG. 16. The relationship between the Earth rotation rate vector monitored by the FOGs and the same vector defined at the l-frame is given using the rotation matrix R_l^bas follows: $\begin{array}{cc}{\omega }_{\mathrm{ie}}^b=\left(\begin{array}{c}{\omega }_x\\ {\omega }_y\\ {\omega }_z\end{array}\right)=R_l^b{\omega }_{\mathrm{ie}}^l=(\begin{array}{ccc}\cos \psi \cos \phi -\sin \psi \sin \theta \sin \phi & \sin \psi \cos \phi +\cos \psi \sin \theta \sin \phi & -\cos \theta \sin \phi \\ -\sin \psi \cos \theta & \cos \psi \cos \theta & \sin \theta \\ \cos \psi \sin \phi +\sin \psi \sin \theta \cos \phi & \sin \psi \sin \phi -\cos \psi \sin \theta \cos \phi & \cos \theta \cos \phi \end{array})\left(\begin{array}{c}0\\ {\omega }^e\cos \varphi \\ {\omega }^e\sin \varphi \end{array}\right)& \left(47\right)\end{array}$

[0191] The angular velocity ω_yis monitored by the toroidal FOG and can be obtained from the corresponding angular increment output as follows: $\begin{array}{cc}{\omega }_y=\frac{\Delta {\theta }_y}{\Delta t}& \left(48\right)\end{array}$

[0192] The component of the Earth rotation rate monitored by the toroidal FOG can be obtained from equation (47) and written as:

[0193] ω_y=(cos ψ cos θ)ω^ecos φ+(sin θ)ω^esin φ  (49)

[0194] This equation is similar to equation 13 for the establishment of the desired azimuth direction. The only difference is that the angular velocity measurement, ω_y, is provided by the toroidal FOG instead of the normal FOG. With simple mathematical manipulation, equation 49 can be rewritten as: $\begin{array}{cc}\cos \psi =\frac{\frac{{\omega }_y}{{\omega }^e\cos \varphi }-\sin \theta \tan \varphi }{\cos \theta }& \left(50\right)\end{array}$

[0195] Since the drilling process for this section of the well is along the horizontal plane, we can assume that the pitch angle is equal to zero and equation 50 can be simplified. However, general solution is presented in this section to suit a special case that might exist while drilling the radical section of the well. This case involves drilling some parts of the radical section with rotary mode of drilling instead of utilizing the steering mode. In this special case the pitch angle cannot be assumed equal to zero.

[0196] To determine the azimuth angle from equation 50, the pitch angle, θ, should be computed. Moreover, complete station based surveying requires the determination of the inclination angle (π/2−θ) and the tool face angle (the roll angle, φ). This information is determined by manipulating the accelerometer measurements. The accelerometer measurements (𝑓_x𝑓_y𝑓_z)^Tare related to the Earth gravity filed through the rotation matrix R_l^b. This relationship can be written as: $\begin{array}{cc}{g}^b=\left(\begin{array}{c}{f}_x\\ f_y\\ f_z\end{array}\right)=R_l^b{g}^l=(\begin{array}{ccc}\cos \psi \cos \phi -\sin \psi \sin \theta \sin \phi & \sin \psi \cos \phi +\cos \psi \sin \theta \sin \phi & -\cos \theta \sin \phi \\ -\sin \psi \cos \theta & \cos \psi \cos \theta & \sin \theta \\ \cos \psi \sin \phi +\sin \psi \sin \theta \cos \phi & \sin \psi \sin \phi -\cos \psi \sin \theta \cos \phi & \cos \theta \cos \phi \end{array})\left(\begin{array}{c}0\\ 0\\ -\gamma \end{array}\right)& \left(51\right)\end{array}$

[0197] The expressions for 𝑓_x, 𝑓_yand 𝑓_zcan be expressed in terms of the pitch and the roll angles by manipulating equation 51 as follows:

[0198] 𝑓_x=(cos θ sin φ)γ  (52)

[0199] 𝑓_y=−(sin θ)γ  (53)

[0200] 𝑓_z=−(cos θ cos φ)γ  (54)

[0201] From equations (52-54) the pitch and the roll angles can be expresses as: $\begin{array}{cc}\theta =\arctan \left(\sqrt{\frac{f_y^2}{f_x^2+f_z^2}}\right)& \left(55\right)\\ \phi =\arctan \left(-\frac{f_x}{f_z}\right)& \left(56\right)\end{array}$

[0202] The specific force measurements 𝑓_x, 𝑓_yand 𝑓_zare obtained from the velocity increment outputs of the accelerometers Δν_x, Δν_yand Δν_zas follows: $\begin{array}{cc}{f}_x=\frac{\Delta v_x}{\Delta t},f_y=\frac{\Delta v_y}{\Delta t},f_z=\frac{\Delta v_z}{\Delta t}& \left(57\right)\end{array}$

[0203] Therefore, station-based surveying during the horizontal section of the well delivers only the orientation of the bottom hole assembly using equations 50, 55 and 56 to provide the azimuth, the pitch and the roll respectively. The inclination is then equal to π/2−θ and the tool face is exactly the roll.

[0204] Error Analysis for the Station-based Surveying

[0205] Although the measurements from the toroidal FOG and the accelerometers are compensated for all bias terms, there are still some residual errors that affect the accuracy of the attitude angles. These residual errors include the cross-axis sensitivity, the bias calibration error, the bias temperature shift, and the scale factor calibration error. Some typical values are assigned to these errors and they are defined by the manufacturer.

[0206] The error components in the specific force measurements can be obtained by differentiating equations 52, 53 and 54.

[0207] δ𝑓_x=(−sin θ sin φδθ+cos θ cos φδφ)γ  (58)

[0208] δ𝑓_y=−(cos θδθ)γ  (59)

[0209] δ𝑓_z=(sin θ cos φδθ+cos θ sin φδφ)γ  (60)

[0210] The error in the pitch angle, δθ, can be obtained from equation 59 and is expressed as follows: $\begin{array}{cc}\delta \theta =-\frac{\left(\frac{\delta f_y}{\gamma }\right)}{\cos \theta }& \left(61\right)\end{array}$

[0211] Using equations 58 and 60, (δ𝑓_x)²+(δ𝑓_z)²can be given as follows:

[0212] (δθ_x)²+(δθ_z)²=(sin²θ(δθ)²+cos²θ(δφ)²)γ²  (62)

[0213] Equation (62) can be manipulated with the substitution of the error in the pitch angle from equation (61) to give the error in the roll angle δφ as follows: $\begin{array}{cc}\delta \phi =\frac{\left(\frac{\sqrt{\delta f_x^2+\delta f_z^2-\delta f_y^2{\tan }^2\theta }}{\gamma }\right)}{\cos \theta }& \left(63\right)\end{array}$

[0214] The errors in both the pitch and the roll are determined from equations 61 and 63 and are used to correct the values of the pitch and the roll angles obtained by equations 55 and 56. The correct value of the pitch angle is then used in equation 50 to determine the azimuth. However, any residual errors existing in the toroidal FOG affect the azimuth accuracy. The azimuth error can be determined from equation 49. If this equation is differentiated with the assumption that the pitch angle is free from error, the toroidal FOG error can be expressed as:

[0215] δω_y=δ(cos ψ)cos θ(ω^ecos φ)  (64)

[0216] Then the error term δ(cos ψ) can be written as $\begin{array}{cc}\delta \left(\cos \psi \right)=\frac{\left(\frac{\delta {\omega }_y}{\left({\omega }^e\cos \varphi \right)}\right)}{\cos \theta }& \left(65\right)\end{array}$

[0217] As discussed above regarding the establishment of the desired azimuth direction, the error in the cos ψ term should be determined using equation 65 and completely removed before the value of the azimuth angle is computed. This is to avoid the magnification of the azimuth error if calculated directly as follows: $\begin{array}{cc}\delta \psi =-\frac{\left(\frac{\delta {\omega }_y}{\left({\omega }^e\cos \varphi \right)}\right)}{\cos \theta }·\frac{1}{\sin \psi }& \left(66\right)\end{array}$

[0218] It can be inferred from equation 66 that the term 1/sin ψ largely magnifies the azimuth error, δψ, as the bottom hole assembly approaches the North-South direction (i.e. when ψ≈0). In addition, singularity problems exist in the calculation of δψ if ψ=0. Therefore, to avoid all these problems, the error term δ(cos ψ) is computed prior to determining the azimuth.

[0219] The whole station-based surveying procedure can be shown as a block diagram in FIG. 17. It should be mentioned that this procedure gives coarse computation of the attitude angles. Precise computation should be applied to improve the estimation accuracy. The technique for precise computation is described in the following sections.

[0220] Error Analysis of the Motion of the Bottom Hole Assembly

[0221] In practical systems employing the technology of strapdown inertial navigation the accuracy is limited by the errors originating in the raw measurements as well as by the imperfections in the various components which are combined to build the system. These error sources are categorized as inertial sensor errors and computational errors [Stephenson et al., 1992; Schwarz et al., 1990]. Any lack of precision in the measurements of an inertial navigation system is passed from one estimate to the next with the overall uncertainty in the precision of the estimated quantity drifting with time. This explains why the accuracy of the INS deteriorates in the long term. For MWD surveying, the drilling process continues for several hours. Therefore, error models are required for analysis and estimation of the various error sources associated with the inertial system proposed in the previous section.

[0222] Errors in dynamic systems are variable in time and therefore described by differential equations. Linearization of the non-linear dynamic system is the most common approach for deriving a set of linear differential equations which define the error states of the dynamic system. The error states expressed in the l-frame for the mechanization equations given in the previous section are defined as follows:

[0223] e^l=(δφ,δλ,δh,δV_e,δV_n,δV_u,δθ,δφ,δψ,δω_x,δω_y,δω_z,δ𝑓_x,δ𝑓_y,δ𝑓_z)  (67)

[0224] where

[0225] δφ,δλ,δh are the coordinate errors;

[0226] δV^e,δVⁿ,δV^uare the velocity errors;

[0227] δθ,δφ,δψ are the attitude errors;

[0228] δω_x, δω_y,δω_zare the errors in the angular velocity measurements;

[0229] δ𝑓_x,δ𝑓_y,δ𝑓_zare the errors in the specific force measurements.

[0230] Attitude Errors

[0231] The attitude errors are defined as the misalignment due to the orthogonal transformation between the b-frame and the l-frame. These errors (δθ,δφ and δψ) depend on two main sources. The first source are the errors in the measurement of angular velocities provided by the FOGs. The second source are the errors generated due to the Earth rotation and the change of orientation of the l-frame. The expression describing the time rate of change of the attitude errors is given as follows:

[0232] {dot over (ε)}^l=−Ω_il^lε^l−δω_il^l+R_b^ld  (68)

[0233] where

[0234] ε^l=(δθ δφ δψ)^Tis the attitude error;

[0235] Ω_il^l=Ω_ie^l+Ω_el^lis the skew symmetric matrix for the rotation of the l-frame with respect to the i-frame;

[0236] δω_il^l=δω_ie^l+δω_el^lis the angular velocity error vector corresponding to ω_il^l;

[0237] d=(δω_xδω_yδω_z)^Tis the vector of angular velocities measurement errors (FOG drift rates).

[0238] The second term in the right hand side of equation 68, δω_il^l, depends implicitly on the velocity errors δV^l=(δV^eδVⁿδV^u)^Tand the coordinate errors δr^l=(δφ δλ δh)^T. This angular velocity error consists of two terms, δω_ie^land δω_el^l, which are given as: $\begin{array}{cc}\delta {\omega }_{\mathrm{ie}}^l=\left(\begin{array}{ccc}0& 0& 0\\ -{\omega }^e\sin \varphi & 0& 0\\ {\omega }^e\cos \varphi & 0& 0\end{array}\right)\left(\begin{array}{c}\delta \varphi \\ \delta \lambda \\ \delta h\end{array}\right)& \left(69\right)\\ \delta {\omega }_{\mathrm{el}}^l=\left(\begin{array}{ccc}0& 0& \frac{V^n}{{\left(M+h\right)}^2}\\ 0& 0& \frac{-V^e}{{\left(N+h\right)}^2}\\ \frac{V^e}{\left(N+h\right){\cos }^2\varphi }& 0& \frac{-V^e\tan \varphi }{{\left(N+h\right)}^2}\end{array}\right)\left(\begin{array}{c}\delta \varphi \\ \delta \lambda \\ \delta h\end{array}\right)+& \left(70\right)\\ \left(\begin{array}{ccc}0& \frac{-1}{M+h}& 0\\ \frac{1}{N+h}& 0& 0\\ \frac{\tan \varphi }{N+h}& 0& 0\end{array}\right)\left(\begin{array}{c}\delta V^e\\ \delta V^n\\ \delta V^u\end{array}\right)& \end{array}$

[0239] With the substitution of the expressions of δω_ie^l(equation 69) and δω_el^l(equation 70), the attitude error equation (equation 68) can be rewritten as: $\begin{array}{cc}\begin{array}{c}\left(\begin{array}{c}\delta \stackrel{.}{\theta }\\ \delta \stackrel{.}{\phi }\\ \delta \stackrel{.}{\psi }\end{array}\right)=\left(\begin{array}{ccc}0& {\omega }^e\sin \varphi +\frac{V^e\tan \varphi }{N+h}& -{\omega }^e\cos \varphi -\frac{V^e}{N+h}\\ -{\omega }^e\sin \varphi -\frac{V^e\tan \varphi }{N+h}& 0& -\frac{V^n}{M+h}\\ {\omega }^e\cos \varphi +\frac{V^e}{N+h}& \frac{V^n}{M+h}& 0\end{array}\right)\\ \left(\begin{array}{c}\delta \theta \\ \delta \phi \\ \delta \psi \end{array}\right)+\left(\begin{array}{ccc}0& 0& \frac{-V^n}{{\left(M+h\right)}^2}\\ {\omega }^e\sin \varphi & 0& \frac{V^e}{{\left(N+h\right)}^2}\\ -{\omega }^e\cos \varphi -\frac{V^e}{\left(N+h\right){\cos }^2\varphi }& 0& \frac{V^e\tan \varphi }{{\left(N+h\right)}^2}\end{array}\right)\left(\begin{array}{c}\delta \varphi \\ \delta \lambda \\ \delta \lambda \end{array}\right)+\\ \left(\begin{array}{ccc}0& \frac{1}{M+h}& 0\\ \frac{-1}{N+h}& 0& 0\\ \frac{-\tan \varphi }{N+h}& 0& 0\end{array}\right)\left(\begin{array}{c}\delta V^e\\ \delta V^n\\ \delta V^u\end{array}\right)+R_b^{l}d\end{array}& \left(71\right)\end{array}$

[0240] Coordinate Errors

[0241] Since the coordinates of the bottom hole assembly depend only on the velocity components V^e, Vⁿand V^u, the coordinate errors δφ, δλ and δh are given directly in terms of velocity errors δV^e,δVⁿand δV^u. The expression for the coordinate errors is obtained by differentiating both sides of equation 2 that describes the relationship between the time rate of change of the coordinate components and the velocity components. Since this relationship is given by the matrix D⁻¹with its elements depending on the coordinate components φ and h, the coordinate errors δφ, δλ and δh are written as follows: $\begin{array}{cc}\begin{array}{c}\delta {\stackrel{.}{r}}^l=\left(\begin{array}{ccc}0& \frac{1}{M+h}& 0\\ \frac{1}{\left(N+h\right)\cos \varphi }& 0& 0\\ 0& 0& 1\end{array}\right)\left(\begin{array}{c}\delta V^e\\ \delta V^n\\ \delta V^u\end{array}\right)+\\ \left(\begin{array}{ccc}0& 0& \frac{-V^n}{{\left(M+h\right)}^2}\\ \frac{V^e\tan \varphi }{\left(N+h\right)\cos \varphi }& 0& \frac{-V^e}{{\left(N+h\right)}^2\cos \varphi }\\ 0& 0& 0\end{array}\right)\left(\begin{array}{c}\delta \varphi \\ \delta \lambda \\ \delta h\end{array}\right)\end{array}& \left(72\right)\end{array}$

[0242] Velocity Errors

[0243] The complete implementation of the l-frame mechanization equations requires the computation of the three velocity components V^e, Vⁿand V^uusing equation 4. The corresponding error equations are obtained by differentiating equation 4 and are given as: $\begin{array}{cc}\left(\begin{array}{c}\delta {\stackrel{.}{V}}^e\\ \delta {\stackrel{.}{V}}^e\\ \delta {\stackrel{.}{V}}^u\end{array}\right)=\left(\begin{array}{ccc}\left(2{\omega }^e\sin \varphi \right)V^u+\left(2{\omega }^e\cos \varphi \right)V^n+\frac{V^n{V}^e}{\left(N+h\right){\cos }^2\varphi }& 0& 0\\ \left(-2{\omega }^e\cos \varphi \right)V^e-\frac{{\left(V^e\right)}^2}{\left(N+h\right){\cos }^2\varphi }& 0& 0\\ \left(-2{\omega }^e\sin \varphi \right)V^e& 0& \frac{2\gamma }{R}\end{array}\right)\left(\begin{array}{c}\delta \varphi \\ \delta \lambda \\ \delta h\end{array}\right)+\left(\begin{array}{ccc}0& f_u& -f_n\\ -f_u& 0& f_e\\ f_n& -f_e& 0\end{array}\right)\left(\begin{array}{c}\delta \theta \\ \delta \phi \\ \delta \psi \end{array}\right)+(\begin{array}{ccc}\frac{-V_u}{N+h}+\frac{V_n}{N+h}\tan \varphi & 2{\omega }^e\sin \varphi +\frac{V^e\tan \varphi }{\left(N+h\right)}& -2{\omega }^e\cos \varphi -\frac{V^e}{N+h}\\ -2{\omega }^e\sin \varphi -\frac{2V^e\tan \varphi }{N+h}& -\frac{V^u}{M+h}& \frac{-V^n}{M+h}\\ 2{\omega }^e\cos \varphi +\frac{2V^e}{N+h}& \frac{2V^n}{M+h}& 0\end{array})\left(\begin{array}{c}\delta V^e\\ \delta V^n\\ \delta V^u\end{array}\right)+R_b^l\left(\begin{array}{c}\delta f_x\\ \delta f_y\\ \delta f_z\end{array}\right)& \left(73\right)\end{array}$