0

-

S
ψ

C
ψ

0

0

0

1

)

⁢

(

X

Y

Z

)

NED

where C represents cosine, S represents sine, ψ represents the azimuth, φ represents the gravity tool face, and θ represents the pitch angle. For a coordinate system employing inclination (rather than pitch angle), the relationship between the NED and original coordinate systems may be expressed, for example, as follows:

${(\begin{matrix} X \\ Y \\ Z \end{matrix})}_{NED} = {R_{NED - O} (\begin{matrix} X \\ Y \\ Z \end{matrix})}_{O} = (\begin{matrix} C_{ψ} & - S_{ψ} & 0 \\ S_{ψ} & C_{ψ} & 0 \\ 0 & 0 & 1 \end{matrix}) (\begin{matrix} S_{I} & 0 & - C_{I} \\ 0 & 1 & 0 \\ C_{I} & 0 & S_{I} \end{matrix}) (\begin{matrix} 1 & 0 & 0 \\ 0 & C_{ϕ} & - S_{ϕ} \\ 0 & S_{ϕ} & C_{ϕ} \end{matrix}) {(\begin{matrix} X \\ Y \\ Z \end{matrix})}_{O}$

${(\begin{matrix} X \\ Y \\ Z \end{matrix})}_{O} = {R_{O - N E D} (\begin{matrix} X \\ Y \\ Z \end{matrix})}_{NED} = (\begin{matrix} 1 & 0 & 0 \\ 0 & C_{ϕ} & S_{ϕ} \\ 0 & - S_{ϕ} & C_{ϕ} \end{matrix}) (\begin{matrix} S_{I} & 0 & - C_{I} \\ 0 & 1 & 0 \\ C_{I} & 0 & S_{I} \end{matrix}) (\begin{matrix} C_{ψ} & S_{ψ} & 0 \\ - S_{ψ} & C_{ψ} & 0 \\ 0 & 0 & 1 \end{matrix}) {(\begin{matrix} X \\ Y \\ Z \end{matrix})}_{NED}$

As disclosed in commonly assigned U.S. Pat. No. 12,123,297, triaxial magnetometer measurements may be modelled, for example, as follows:

${(\begin{matrix} B_{x} \\ B_{y} \\ B_{z} \end{matrix})}_{i} = (\begin{matrix} 1 & 0 & 0 \\ 0 & 1 & - γ_{s} ω_{s} \\ 0 & γ_{s} ω_{s} & 1 \end{matrix}) (\begin{matrix} 1 & 0 & 0 \\ 0 & 1 & - γ_{c} ω_{c} \\ 0 & γ_{c} ω_{c} & 1 \end{matrix}) (\begin{matrix} \end{matrix}) + (\begin{matrix} b_{x} \\ b_{y} \\ b_{z} \end{matrix}) + (\begin{matrix} ε_{B_{x}} \\ ε_{B_{y}} \\ ε_{B_{z}} \end{matrix})$

where B_x, B_y, and B_zrepresent the triaxial magnetometer measurements in the original coordinate system at survey location i, {circumflex over (B)}_x, {circumflex over (B)}_y, and {circumflex over (B)}_z, represent the true magnetic field vector (or true magnetic field measurements representative of reality), b_x, b_y, and b_zrepresent magnetometer bias, γ_cand γ_srepresent eddy current compensation terms for drill collar rotation γ_cand roll-stabilized sensor housing rotation γ_s, ω_cand ω_srepresent the rotation rates (angular frequency) of the drill collar and sensor housing, and ε_B_x, ε_B_y, and ε_B_zrepresent error terms.

U.S. Pat. No. 12,123,297 further disclosed that an interfering magnetic field may result from eddy currents in the rotating drill collar and/or sensor housing. Rotation of the collar 62 and/or housing 70 in the Earth's magnetic field may generate an eddy current in the collar 62 and/or housing 70 (owing to the Lorentz force created by the Earth's field penetrating rotating collar and/or housing). The interfering field is directed in the radial direction such that the measured field may deviate from the Earth's field. A methodology was disclosed to compensate for such eddy current interference.

The '297 Patent further disclosed a method for compensating magnetometer bias (correcting for b_x, b_y, and b_z) that requires inclination measurements. In one disclosed method a measured wellbore inclination, an eddy current compensation term, and the magnetometer measurements are processed with a Kalman filter to remove the magnetometer bias. In another disclosed method, wellbore inclination measurements, gravity tool face measurements, and well azimuth are input into a downhole multi-station analysis model to compute the magnetometer bias. While the disclosed methods were found to provide excellent magnetometer bias compensation and to enable survey quality wellbore azimuth measurements, there is room for further improvement. For example, in near vertical wells (e.g., at wellbore inclinations of less than about two degrees, less than about three degrees, less than about five degrees, or less than about ten degrees), the bias compensation can be prone to error (since inclination measurements and gravity tool face measurements themselves are prone to error at low inclination). There is a need for an improved (or new) method of magnetometer bias compensation in near vertical wells.

The magnetometer bias may be thought of as a semi-constant offset parameter. For example, at the survey station (or location) i, if the true magnetic field is ( [Image Omitted] )^T, then the magnetometer bias offset may be expressed as follows:

$b = B_{i} -$

$where :$

$b = (\begin{matrix} b_{x} \\ b_{y} \\ b_{z} \end{matrix})$

$B_{i} = {(\begin{matrix} B_{x} \\ B_{y} \\ B_{z} \end{matrix})}_{i}$

$= (\begin{matrix} \end{matrix})$

One aspect of the disclosed embodiments was the realization that the magnetometer bias may be estimated by using magnetometer measurements and corresponding accelerometer measurements to compute an inclination independent magnetic field parameter, such as the magnetic dip, and then compare the computed parameter with a reference parameter value. By inclination independent it is meant that the wellbore inclination is not used (or otherwise needed) to compute the magnetic field parameter. For example, it was realized that the computed magnetic dip angle may be compared with a reference magnetic dip angle for a particular drilling location. It was then further realized that the magnetometer bias may be estimated via inversion by minimizing the difference (e.g., the root mean square difference) between a large number of measured (computed) magnetic dip angles and the reference magnetic dip angle.

FIG. 5 depicts a flowchart of one example method 100 for estimating magnetometer bias while drilling and optionally determining a dynamic wellbore azimuth. Method 100 includes rotating a drilling system in a wellbore (e.g., to drill) at 102. The drilling system may include a BHA having a drill bit and an RSS including a drill collar and a roll-stabilized control unit, e.g., as described above with respect to FIGS. 2 and 3. Triaxial magnetometer (magnetic field) measurements and triaxial accelerometer (gravitational field) measurements are made using corresponding sensors located in a roll-stabilized housing at 104. Rotation rates of the drill collar and/or the sensor housing may also be measured at 104. Magnetometer bias dependent magnetic dip angles are computed at 106 using the magnetometer and accelerometer measurements made at 104 and an unknown magnetometer bias. The computed magnetic dip angles are then compared with a reference dip angle for the drilling location at 108. A magnetometer bias (e.g., a bias vector including b_x, b_y, and b_z) is selected at 110 (e.g., using an inversion algorithm) that minimizes the difference (or error) between the computed magnetic dip angles and the reference dip angle. A wellbore azimuth may optionally be computed at 112 from the accelerometer measurements, the magnetometer measurements, and the computed magnetometer bias.

With continued reference to FIG. 5, the true magnetic field signal may be computed from the magnetic field measurements and the magnetometer bias (assuming no other magnetic interference in this particular example), for example, as follows:

$= B_{i} - b$

Assuming no accelerometer bias, the true gravitational field given is as follows:

$= G_{i}$

The magnetic dip depends on the true magnetic and gravitational fields (such that it depends on the accelerometer measurements, the magnetometer measurements, and the magnetometer bias) and may be expressed mathematically, for example, as follows:

$D_{i} (b) = \sin^{- 1} (\frac{}{❘ ❘ ❘ ❘}) = \sin^{- 1} (\frac{{G_{ι}}^{T} (B_{i} - b)}{❘ G_{i} ❘ ❘ B_{i} - b ❘})$

Note that D_idepends on the accelerometer measurements G_i, the magnetometer measurements B_iand the magnetometer bias b and is independent of the wellbore inclination (such that magnetic dip is an inclination independent parameter). It will be appreciated that magnetic dip is also independent of the wellbore azimuth, the toolface, the total magnetic field, and the total gravitational field. Differences between the computed magnetic dip angle and a reference magnetic dip angle D_refmay be computed, for example, as follows:

$e_{i} (b) = D_{i} (b) - D_{r e f}$

Those of ordinary skill in the art will readily appreciate that the magnetic dip angle (also referred to as the magnetic inclination) is the angle made with the horizontal by Earth's magnetic field lines with positive angular values indicating that the magnetic field is pointing into the Earth and negative values indicating that it is pointing out of the Earth. Those of ordinary skill will also readily appreciate that the magnetic dip angle strongly depends on the measurement location on the surface of the Earth (particularly on the latitude). It will still further appreciate that the magnetic dip angle at any location on Earth (e.g., as defined via the longitude and the latitude of the location) can be readily obtained using lookup tables or online calculators (e.g., from the U.S. National Centers for Environmental Information at https://www.ngdc.noaa.gov/geomag/calculators/magcalc.shtml). Such obtained magnetic dip angles are referred to herein as reference magnetic dip angles.

As described above, the foregoing equations may be inverted for the magnetometer bias that minimizes the errors e_i(such as a sum of the errors or a sum of the absolute value of the errors). In one particularly advantageous example, a magnetometer bias b=(b_xb_yb_z)^Tmay be selected that minimizes the root mean square error (difference) between the computed magnetic dip angle values and the reference value. The squared sum of the error ϵ may also be expressed, for example, as follows:

$ϵ (b) = e^{T} e$

Here, the error vector may be as defined, for example, as below:

$e = [\begin{matrix} e_{1} \\ ⋮ \\ e_{n} \end{matrix}] = [\begin{matrix} D_{1} - D_{ref} \\ ⋮ \\ D_{n} - D_{ref} \end{matrix}] = D (b) - D_{r e f}$

$where : D (b) = [\begin{matrix} D_{1} (b) \\ ⋮ \\ D_{n} (b) \end{matrix}], D_{ref} = D_{r e f} [\begin{matrix} 1 \\ ⋮ \\ 1 \end{matrix}]$

such that the magnetometer bias b may be selected by inversion as follows:

$b = \arg \min_{b} ϵ (b)$

FIG. 6 depicts a flowchart of another example method 120 for estimating magnetometer bias while drilling and optionally determining a dynamic wellbore azimuth. Method 120 is similar to method 100 in that it includes rotating a drilling system in a wellbore (e.g., to drill) at 122 and making triaxial magnetometer measurements and triaxial accelerometer measurements using corresponding sensors located in a roll-stabilized (or slowly rotating) housing at 124. Rotation rates of the drill collar and/or the sensor housing may also be measured at 124. An eddy current compensation is computed at 126, for example, from measured rotation rate(s) of the drill collar and/or the roll-stabilized housing. The accelerometer measurements, the magnetometer measurements, and the eddy current compensation are input into an inclination independent model that enables magnetometer bias dependent magnetic dip angles to be computed at 128. The inclination independent model may be further advantageously azimuth independent and gravity toolface independent. Moreover, the inclination independent model may be still further independent of the total magnetic and gravitational fields. The computed magnetometer bias dependent magnetic dip angles are compared with a reference dip angle for the drilling location at 130. A magnetometer bias (e.g., a bias vector including b_x, b_y, and b_z) is selected at 132 (e.g., using an inversion algorithm) that minimizes the difference (or error) between the computed magnetic dip angles and the reference dip angle. A wellbore azimuth may be optionally computed at 134 using the accelerometer measurements, the magnetometer measurements, the eddy current compensation term, and the selected magnetometer bias.

With continued reference to FIG. 6, the eddy current compensation term y may be estimated, for example, using the methods disclosed in the '297 Patent. For example, the eddy current compensation term y may be estimated from a change in angle X when the rotation rate changes. It will be appreciated that angle X is the angle between the gravity and magnetic field vectors in the y-z plane (the cross axial plane perpendicular to the axis of the BHA) and may be computed, for example, as follows:

$X = \tan^{- 1} (\frac{B_{z} G_{y} - B_{y} G_{z}}{B_{y} G_{y} + B_{z} G_{z}})$

In the absence of eddy currents, angle X is essentially constant at a survey station (or location). Since angle X has been found to change with changing collar rotation rate (e.g., increase with increasing rotation rate), this dependency may be used to estimate the eddy current compensation term γ (and to estimate changes in the eddy current compensation term with a changing rotation rate of the collar). For example, the eddy current compensation term γ may be expressed as follows:

$γ = α \frac{1 - \sqrt{1 - 4 \frac{ω^{2}}{α^{2}}}}{2 ω^{2}}$

where ω represents the collar rotation rate and α=(γ²ω²+1)/γ.

While the preceding equation provides a suitable solution for the eddy current compensation term γ, a simplified solution may be obtained by recognizing that the eddy current compensation term may be approximated as ∂X_ω/∂X_ω, for example, as follows when ω²<<α²(e.g., when ω²/α²approaches zero):

$γ = α \frac{1 - (1 - 2 \frac{ω^{2}}{α^{2}})}{2 ω^{2}} = \frac{1}{α} = \frac{\partial X_{ω}}{\partial ω}$

This approximation provides for a more robust computation of y and advantageously has an error of less than 1 percent for most drilling conditions (e.g., collar rotation rates of less than about 300 rpm).

With still further reference to FIG. 6, the magnetometer bias may be determined at 128, 130, 132 in a manner similar to that described above with respect to method 100 (FIG. 5). In the presence of eddy current interference, the true magnetic field may be expressed, for example, as follows:

$= R_{eddy} (B_{i} - b) - w_{B_{i}}$

where, as before, [Image Omitted] represents the true magnetic field vector, B_irepresent the magnetometer measurements, b represent the magnetometer bias vector, and w_B_iis random noise with zero mean. R_eddyrepresents a correction matrix that includes the computed eddy current compensation term γ and may be given, for example, as follows:

$R_{e d d y} = (\begin{matrix} 1 & 0 & 0 \\ 0 & f_{a} & f_{b} \\ 0 & - f_{b} & f_{a} \end{matrix})$

where, in some embodiments, ƒ_aand ƒ_bmay be mathematical relations of the eddy current compensation term γ such that ƒ_acompensates the amplitude of the eddy current interference and ƒ_bcompensates the phase of the eddy current interference. Assuming that there is only negligible bias in the corresponding accelerometer measurements:

$= G_{i} - w_{G_{i}}$

where [Image Omitted] represents the true gravitational field vector, G_irepresent the accelerometer measurements, and w_G_iis random noise with zero mean. As before, the magnetic dip angle may be expressed based on the true magnetic field and the true gravitational field (not including the random noises), for example, as follows:

$D_{i} = f_{i} = f (+ w_{B_{i}}, + w_{G_{i}}) = f (R_{eddy} (B_{i} - b), G_{i})$

$where :$

$f (B, G) \equiv \sin^{- 1} (\frac{G^{T} B}{❘ G ❘ ❘ B ❘})$

such that:

$D_{i} = \sin^{- 1} (\frac{G_{i}^{T} (R_{e d d y} (B_{i} - b))}{❘ G_{i} ❘ ❘ R_{e d d y} (B_{i} - b) ❘})$

In example embodiments a nonlinear optimization may be carried out using Gauss-Newton. However, it will be appreciated that the random noises w_B_iand w_G_iof magnetometer and accelerometer measurements can potentially invalidate local linearity requirements of the equation if the noise levels are high. A Monte Carlo sampling was carried out to understand the effect of random noises on the computed dip angles using the foregoing equation. In the example simulation, wellbore inclination and wellbore azimuth values were sampled over a uniform distribution from 0 to 360 degrees, the magnetic dip angle varied from −60 to 60 degrees, the total magnetic field varied from 40000 to 60000 nT and the total gravitational field was 1000 mG. Random Gaussian noises were assumed with amplitudes of 50 nT for the magnetometers and 10 mG for the accelerometers. It will be appreciated that these assumed noises are larger than what is generally observed in practice.

The dip angle error dD resulting from the magnetometer and accelerometer random noises w_B_iand w_G_imay be defined as given below:

$d D = f (+ w_{B_{i}}, + w_{G_{i}}) - D_{r e f}$

It will be appreciated that the original dip angle equation may be further simplified, for example, as given below:

$f_{i} = f (,) + w_{D_{i}} = D_{ref}$

where w_Drepresents a random dip angle noise (with both accelerometer measurement noise and magnetometer measurement noise contributions) with zero mean such that:

$D_{i} = f (,) = D_{r e f} - w_{D_{i}}$

In the foregoing equation, the function ƒ contains data samples from the magnetometers and accelerometers, which are a series of fixed value measurements. With continued reference to FIG. 6, the foregoing equation may be expressed with respect to those measurements, for example, as follows:

$f (R_{e d d y} (B_{i} - b), G_{i}) = D_{r e f} - w_{D_{i}}$

Since, B_iand G_iare fixed valued magnetometer and accelerometer measurements and since R_eddyis a three-by-three matrix that is defined by the previously determined eddy current compensation term (at 126), the only remaining variable is the magnetometer bias b. This may be expressed mathematically in a simpler form, for example, as follows:

$f^{'} (b) = D_{r e f} - w_{D_{i}}$

Optimizing b to minimize the squared error from the true dip angle yields the following matrix equation:

$(e = {(F (b) - D)}^{T} (F (b) - D)$

$where :$

$F (b) = (\begin{matrix} f_{1}^{'} (b) \\ ⋮ \\ f_{n}^{'} (b) \end{matrix}), D = D_{ref} (\begin{matrix} 1 \\ ⋮ \\ 1 \end{matrix}), w_{D_{i}} = (\begin{matrix} w_{D_{i}} \\ ⋮ \\ w_{D_{n}} \end{matrix})$

The magnetometer bias vector b which minimizes e may be determined by an iteration process, such as the Gauss Newton method, for example, as defined below:

$b_{s + 1} = b_{s} + {(J_{F}^{T} J_{F})}^{- 1} J_{F}^{T} (D - F (b_{s}))$

where J_Fis the Jacobian matrix of F such that:

$J_{F} = \frac{\partial F (b)}{\partial b} = (\begin{matrix} \frac{\partial f_{1}^{'} (b)}{\partial b} \\ ⋮ \\ \frac{\partial f_{n}^{'} (b)}{\partial b} \end{matrix}) = (\begin{matrix} \frac{\partial f_{1}^{'} (b)}{\partial b_{x}} & \frac{\partial f_{1}^{'} (b)}{\partial b_{y}} & \frac{\partial f_{1}^{'} (b)}{\partial b_{z}} \\ ⋮ & ⋮ & ⋮ \\ \frac{\partial f_{n}^{'} (b)}{\partial b_{x}} & \frac{\partial f_{n}^{'} (b)}{\partial b_{y}} & \frac{\partial f_{n}^{'} (b)}{\partial b_{z}} \end{matrix})$

An analytical expression that defines ∂ƒ′(b)/∂b is given below:

$\frac{\partial f}{\partial b} = \frac{1}{A} {(\begin{matrix} r_{x} r_{y} (G_{y} f_{a} - G_{z} f_{b}) + r_{x} r_{z} (G_{y} f_{b} + G_{z} f_{a}) - G_{x} (r_{y}^{2} + r_{z}^{2}) \\ G_{x} r_{x} r_{y} + r_{y} r_{z} (G_{y} f_{b} + G_{z} f_{a}) - (r_{x}^{2} + r_{z}^{2}) (G_{y} f_{a} - G_{z} f_{b}) \\ G_{x} r_{x} r_{y} + r_{y} r_{z} (G_{y} f_{a} - G_{z} f_{b}) - (r_{x}^{2} + r_{y}^{2}) (G_{y} f_{b} + G_{z} f_{a}) \end{matrix})}^{T}$

$where :$

$r_{x} = B_{x} - b_{x}$

$r_{y} = B_{y} - b_{y}$

$r_{z} = B_{z} - b_{z}$

$A = (r_{x}^{2} + r_{y}^{2} + r_{z}^{2}) (r_{x}^{2} + r_{y}^{2} + r_{Z}^{2}) (G_{x}^{2} + G_{y}^{2} + G_{Z}^{2}) - {(G_{x} r_{x} + G_{y} r_{y} + G_{Z} r_{Z})}^{2}$

Note that in the foregoing equation r