Numerical evaluation of magnetic absolute measurements with arbitrarily distributed DI-fluxgate theodolite orientations

At geomagnetic observatories the absolute measurements are needed to determine the calibration parameters of the continuously recording vector magnetometer (variometer). Absolute measurements are indispensable for determining the vector of the geomagnetic field over long periods of time. A standard DI (declination, inclination) measuring scheme for absolute measurements establishes routines in magnetic observatories. The traditional measuring schema uses a fixed number of eight orientations (Jankowski et al., 1996).

We present a numerical method, allowing for the evaluation of an arbitrary number (minimum of five as there are five independent parameters) of telescope orientations. Our method provides D, I and Z base values and calculated error bars of them.

A general approach has significant advantages. Additional measurements may be seamlessly incorporated for higher accuracy. Individual erroneous readings are identified and can be discarded without invalidating the entire data set. A priori information can be incorporated. We expect the general method to also ease requirements for automated DI-flux measurements. The method can reveal certain properties of the DI theodolite which are not captured by the conventional method.

Based on the alternative evaluation method, a new faster and less error-prone measuring schema is presented. It avoids needing to calculate the magnetic meridian prior to the inclination measurements.

Measurements in the vicinity of the magnetic equator are possible with theodolites and without a zenith ocular.

The implementation of the method in MATLAB is available as source code at the GFZ Data Center (Brunke, 2017).

1 Introduction

Absolute measurements of the magnetic declination D and inclination I are taken by means of a non-magnetic theodolite with a fluxgate sensor mounted on its telescope that is nearly parallel to the optical axis. The reading of the magnetometer S becomes null when it points in a direction perpendicular to the field. A standard DI measuring scheme established routines in magnetic observatories. It uses eight such telescope orientations for absolute measurements. This standard DI scheme allows for a simple numeric evaluation and cancels out the influence of intrinsic instrument errors like sensor offset and misalignment angles between fluxgate sensor and telescope (hereafter referred to as instrument parameters).

The alternative method for DI measurements was motivated by a very practical reason. Zenith oculars are very seldom available. Without them, inclination measurements are not possible around the magnetic equator where the magnetic field is almost horizontal. We wanted to enable measurements with a DI-flux theodolite and without a zenith ocular. The need to evaluate DI measurements at arbitrary orientations was felt before we started our work on a viable solution. Peter Crosthwaite (Crosthwaite, 1994) and Anna Willer (née Nilsson) (Nilsson, 2010) gave us useful personal notes on their approaches. Our work is partly based on their notes but overcomes potential convergence problems faced in the previous work.

We present a numerical method, allowing for the evaluation of an arbitrary number (minimum of five as there are five independent parameters) of telescope orientations. No prescribed fixed measuring scheme is needed. We implement an instrument model to calculate the fluxgate reading S, dependent on field, instrument parameters and telescope orientation. The instrument model does not cover more instrument parameters than misalignment angles and offsets. Our method does not extend the physical description of a DI theodolite, but it allows for orientations of the ones used traditionally for D and I measurements. Inserting measured values for each telescope orientation gives one non-linear equation for each orientation. Eventually this results in an overdetermined system of non-linear equations. This system is solved in the sense of the least square solution using the Gauss–Newton method generalized to an overdetermined system. The accuracy of the resulting D and I base values depends both on the choice of the telescope directions and on the accuracy of the measurements. The method provides estimated variances for the resulting D and I base values. The calculated residuals provide a measure of quality for each individual measurement.

We took advantage of methods known from geophysical inversion theory such as assessing the residuals, using a priori information and objectively assessing the accuracy of the results (Schmucker, 1975; Tarantola, 1982). Exploitation of the method in practice also showed benefits for the routine work at observatories. We show, for example, that at a given rate of erroneous readings, a higher percentage of successful absolute measurements can be achieved. We apply our method to various data sets. We show the benefit to accuracy and reliability of routine absolute measurements at the Niemegk Observatory and the resulting base values.

2 The conventional D and I measurement

2.1 The measuring schema

For the D measurement, the telescope is put to a horizontal orientation (reading of vertical circle: 90^∘). Then the two orientations with the telescope pointing to the magnetic east and the west are determined by adjusting the magnetometer reading S to a small value. This is repeated with the telescope flipped over (vertical reading: 270^∘; sensor on the other side of the telescope). The four resulting readings of the horizontal circle (D readings ${\mathit{\phi }}_{\mathrm{1}},\mathrm{\dots },{\mathit{\phi }}_{\mathrm{4}}$) are used to calculate the direction of magnetic north (Eq. 1). Subsequently the horizontal circle is adjusted to the magnetic north, and the two orientations on the vertical circle are determined where the magnetometer reading S is null. This is finally repeated with the horizontal circle adjusted to the magnetic south. These last four measurements ${\mathit{\xi }}_{\mathrm{1}},\mathrm{\dots },{\mathit{\xi }}_{\mathrm{4}}$ serve to calculate the inclination I in Eq. (2).

This conventional DI scheme has two major advantages:

• a.

  It allows for a simple numeric evaluation by just averaging four readings.

• b.

  The influence of the instrument parameters like sensor offset and the misalignment angles between the fluxgate sensor and the telescope cancel out.

A detailed description of the measurement and its evaluations can be found in (Jankowski and Sucksdorff, 1996; Kerridge, 1988; KringLauridsen, 1985; Matzka and Hansen, 2007).

In this paper we refer to the angle φ as the geographic direction. The angle φ is obtained from the actual reading of the horizontal circle, using a mark with known azimuth. Below we will use the symbols N_↑, N_↓, E_↑, E_↓, S_↑, S_↓, W_↑, W_↓, denoting the theodolite orientations according to the next cardinal point of the compass directions the telescope is pointing to. The arrow indicates whether the magnetometer sensor is situated above or below the telescope.

2.2 The evaluation

In the conventional scheme the calculations of declination D₁ and inclination I₁ are principally restricted to calculating mean values; see (Jankowski and Sucksdorff, 1996).

Declination and inclination usually vary during the measuring process. Let D₁ and I₁ be the values of declination and inclination at the time of the first magnetometer reading S₁. The successive horizontal readings ${\mathit{\phi }}_i,i=\mathrm{2}\mathrm{\dots }\mathrm{4}$ have to be corrected for the natural variation of the declination ΔD_i (reduction to the first measurement). The natural variation ΔD_i from the first measurement D₁ are known from the variometer. Additionally a correction ΔA_i is needed if the magnetometer reading S_i is not exactly zero (KringLauridsen, 1985). Assuming subsequently that D₁, I₁ and all measured angles are given in degrees, the corrected horizontal circle readings are

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{A}_{X_i}={\mathit{\phi }}_i-\mathrm{\Delta }{D}_i-\mathrm{\Delta }{A}_i\\ {\displaystyle }& {\displaystyle }={\mathit{\phi }}_i-\mathrm{\Delta }{D}_i-{\displaystyle \frac{\mathrm{360}{}^{\circ }}{\mathrm{2}\mathit{\pi }}}{\displaystyle \frac{S_i}{H}}\mathrm{with}\\ {\displaystyle }& {\displaystyle }{X}_i\in \mathit{\{}W\uparrow ,W\downarrow ,E\uparrow ,E\downarrow \mathit{\}};i=\mathrm{1}\mathrm{\dots }\mathrm{4}.\end{array}$$

The $X_i\in \mathit{\{}W\uparrow ,E\uparrow ,W\downarrow ,E\downarrow \mathit{\}}$ denotes the four possible horizontal telescope orientations perpendicular to the magnetic field. The telescope looks towards the indicated direction (W or E). The arrow (↑ or ↓) indicates whether the sensor is on the upper or lower side of the telescope. The exact orientations are found by adjusting the telescope horizontally to a small magnetometer readings S_i, ideally to zero. Eventually D₁ is determined as the mean value of $A_{X_i}$ plus 180^∘.

$$\begin{array}{}\text{(1)}& D_{\mathrm{1}}={\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}(A_{W\uparrow }+A_{W\downarrow }+A_{E\uparrow }+A_{E\downarrow })+\mathrm{180}{}^{\circ }.\end{array}$$

Using this declination value D₁, the telescope is now directed into the magnetic meridial plane. The inclination I₁ is determined using the four telescope orientations within the meridial plane, where magnetometer readings are null:

$$\begin{array}{}\text{(2)}& I_{\mathrm{1}}={\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}(V_{N\uparrow }+V_{S\downarrow }-V_{S\uparrow }-V_{N\downarrow })\pm \mathrm{90}{}^{\circ }.\end{array}$$

In Eq. (2) the + and − signs apply to the northern and southern hemispheres, respectively.

$$\begin{array}{ll}{\displaystyle }{V}_{X_i}& {\displaystyle }={\mathit{\xi }}_i-\mathrm{\Delta }{I}_i-\mathrm{\Delta }{V}_i;X_i\in \mathit{\{}N\uparrow ,S\downarrow ,S\uparrow ,N\downarrow \mathit{\}}\\ {\displaystyle }& {\displaystyle }={\mathit{\xi }}_i-\mathrm{\Delta }{I}_i-{\displaystyle \frac{\mathrm{360}{}^{\circ }}{\mathrm{2}\mathit{\pi }}}{\displaystyle \frac{S_i}{F}};i=\mathrm{5}\mathrm{\dots }\mathrm{8}\end{array}$$

Just as above, ΔI_i are corrections for field variations after the first measurement. The ΔV_i values account for non-zero magnetometer readings S_i. These I measurements, as well as the D measurements, can principally be taken in any order. A sequence is preferable, reaching each orientation by inverting the theodolite around just one axis.

Eight correct measurements are needed to apply Eqs. (1) and (2). If data from just a single orientation are corrupted, the entire measurement set is usually discarded, including the good data. Repeating this set of eight measurements several times at the observatories is good practice. We will show later that our method allows a single corrupted theodolite orientation to be discarded without needing to lose the rest.

Figure 1 gives a graphic presentation of all possible telescope orientations of a DI flux theodolite. The ordinate axis gives a reading of the vertical circle of the theodolite. Values of the abscissa are geographic headings (readings of the horizontal circle corrected using an azimuth mark). The four colours indicate the cardinal direction of the heading of the telescope. The faded colours indicate orientations requiring a zenith eye piece. The arrows indicate whether the sensor is above or below the telescope as in Eqs. (1) and (2). Dashed lines give orientations with null magnetometer reading S. These lines are specific to each observatory (indicated for each line) and depend on its local declination and inclination (Niemegk (NGK): $D_{\mathrm{NGK}}=\mathrm{3.6}{}^{\circ }$, $I_{\mathrm{NGK}}=\mathrm{67.5}{}^{\circ }$; Hyderbad (HYD): $D_{\mathrm{HYD}}=-\mathrm{0.4}{}^{\circ }$, $I_{\mathrm{HYD}}=\mathrm{24.3}{}^{\circ }$ and Tatuoca (TTB): $D_{\mathrm{TTB}}=-\mathrm{20.1}{}^{\circ }$, $I_{\mathrm{TTB}}=\mathrm{0.5}{}^{\circ }$). At a given observatory useful measurements are close to these lines. Blue and red marks indicate measuring orientations. Orientations from conventional DI measurement schemes are given on the left panel. D measurements are in red; I measurements are in blue. These orientations depend on the observatory. The right panel shows orientations of a new scheme presented in Sect. 2.3 as carried out at NGK.

Bildunterschrift

Figure 1 Graphic presentation of possible telescope attitudes of a DI flux theodolite. (a) Dashed lines give orientations with null magnetometer reading S. These lines are specific to the observatories of Niemegk (NGK), Hyderbad (HYD) and Tatuoca (TTB). Useful sensor orientations at a given observatory are close to this line. The marks represent measurements at these observatories according to the conventional scheme. Note the different φ values due to different local declination. (b) Measurements taken at the Niemegk Observatory using the new scheme.

2.3 A new scheme for absolute measurements, easy, precise and less error prone

The new scheme allows us to take DI measurements using a simplified and less error-prone method. Advantages with respect to the conventional methods are as follows:

• No need to calculate the magnetic meridian.

• The method tolerates imperfect levelling of the telescope.

• The measurement scheme facilitates adjusting the theodolite. Throughout the entire measurement procedure only the adjustment wheels for the horizontal circle are used to adjust the magnetometer reading to a small value.

• There is no need for telescope orientations along the magnetic meridian. This avoids orientations requiring a zenith eye piece at observatories close to the magnetic equator.

• We feel that, after some practice, the new method (including the optional measurements) needs slightly less time than the conventional, but even without gaining time, the resulting six measured orientations lead to a statistically firmer result than the four conventional ones.

• Repeating the new method with different tilt angles α (explained subsequently) allows systematical errors to be assessed.

Our scheme uses the procedure known from D measurements in the conventional DI scheme. The only difference is that the telescope is not necessarily oriented horizontally. A defined vertical angle ξ off the horizontal plane is used instead. The new scheme is just a repetition of the following two steps:

1. Set vertical circle to a certain value ξ.

2. For the given ξ value (vertical circle), adjust the magnetometer reading S to 0 (or a small value) using the horizontal circle (coarse and fine adjustment). Note the readings of both circles, magnetometer reading S and time of S measurement.

Subsequently we will we call the second step “adjust S horizontally”.

The first four orientations of the new scheme are identical to conventional D measurements. Using the procedure defined above, the new scheme reads as follows:

1. Set ξ to 90^∘ with the telescope pointing towards east. Adjust S horizontally.

2. Invert the telescope, set ξ to 270^∘, telescope pointing towards west. Adjust S horizontally.

3. Turn telescope to east, set ξ to 270^∘, telescope pointing towards east. Adjust S horizontally.

4. Invert the telescope, set ξ to 90^∘, telescope pointing towards east. Adjust S horizontally.

The next measurements are taken with the telescope tilted by a certain angle α and −α off the horizontal plane. The choice of α is not critical, but it must be lower than $\mathrm{90}{}^{\circ }-I$. Otherwise an orientation with S being null cannot be found. The angle α should not be too small in order to get orientations close enough to the magnetic N and S directions. At Niemegk the inclination is I≈67.5^∘, and $\mathit{\alpha }=\mathrm{20}{}^{\circ }$ is a good choice.

• 5.

  Set ξ to $(\mathrm{90}{}^{\circ }-\mathit{\alpha })$ with the telescope pointing towards north-north-east. Adjust S horizontally.

• 6.

  Invert telescope, set ξ to $(\mathrm{270}{}^{\circ }-\mathit{\alpha })$ with the telescope pointing towards south-south-west. Adjust S horizontally.

  • a.

    Optional: invert telescope, set ξ to $(\mathrm{90}{}^{\circ }-\mathit{\alpha })$ with the telescope pointing towards north-north-west. Adjust S horizontally.

  • b.

    Optional: invert telescope, set ξ to $(\mathrm{270}{}^{\circ }-\mathit{\alpha })$ with the telescope pointing towards south-south-east. Adjust S horizontally.

• 7.

  Set ξ to $(\mathrm{90}{}^{\circ }+\mathit{\alpha })$ with the telescope pointing towards south-south-east. Adjust S horizontally.

• 8.

  Invert telescope, set ξ to $(\mathrm{270}{}^{\circ }+\mathit{\alpha })$ with the telescope pointing towards north-north-west. Adjust S horizontally.

  • a.

    Optional: invert telescope, set ξ to $(\mathrm{90}{}^{\circ }+\mathit{\alpha })$ with the telescope pointing towards south-south-west. Adjust S horizontally.

  • b.

    Optional: invert telescope, set ξ to $(\mathrm{270}{}^{\circ }+\mathit{\alpha })$ with the telescope pointing towards north-north-east. adjust S horizontally.

In steps 1 to 4 there is no need to set ξ exactly to values of 90 and 270^∘. Values close to these can also be used but have to be noted, of course. It is a matter of personal preference, because more effort is used when setting ξ to a given value than when reading and noting an exact arbitrary value. The latter is slightly more effort, but reading a value results in more precise data than when setting it.

The right panel of Fig. 1 shows the distribution of data points of the new scheme.

3 The general method for arbitrarily orientated telescope orientations

As explained above, the eight orientations of the standard DI scheme are used due to practical reasons. There is no principal reason to take the measurements at exactly the orientations as in the conventional scheme. Theoretically only five measurements are needed for the five unknown quantities: declination D₁, inclination I₁, the two angles of misalignment δ and ϵ and the magnetometer offset S_off.

The general method described below allows measurements at arbitrary telescope orientations to be evaluated. The number of measurements are arbitrary. Usually it should be at least five. Even fewer than five measurement can be sufficient if a priori information about the instrument parameters is available. The method has a lot of additional advantages. The quality of each single measurement can be assessed and outliers can be identified and discarded. As already mentioned, the measuring set of the eight orientations for the conventional scheme is usually repeated several times, so that 16 or even 24 orientations are available. Using all of them in one joint evaluation leads to a statistically firmer result and better allows us to identify outliers. Even though actual DI measurements at the Niemegk Observatory are known to be at a high accuracy level, we show that its accuracy and reliability could still be improved using the new method. Additionally the new method gives the observer more insight into sources of problems with magnetic cleanliness.

3.1 Theory and numerical approach

3.1.1 Instrument model

The DI-flux instrument model allows the magnetometer reading S to be calculated as a function of magnetic declination D, inclination I, total field strength F, readings of the theodolite orientation φ and ξ (horizontal circle and vertical circle), the misalignment angles of the sensor δ and ϵ and the magnetometer offset S_off. The angles δ and ϵ (misalignment in φ and ξ direction) are also called aberration errors. Crosthwaite (Crosthwaite, 1994) gave a formulation based on vector geometry. The following formulation and its derivation based on spherical trigonometry can be found in (Nilsson, 2010):

$$\begin{array}{ll}\text{(3)}& {\displaystyle }S& {\displaystyle }=f(D,I,\mathit{\delta },\mathit{ϵ},S_{\mathrm{off}},\mathit{\phi },\mathit{\xi }){\displaystyle }& {\displaystyle }=c\cdot F\cdot (-\sin (I)\cdot \cos (\mathit{\xi }+\mathit{ϵ})\\ {\displaystyle }& {\displaystyle }+\cos \left(I\right)\cdot \sin (\mathit{\xi }+\mathit{ϵ})\cdot \cos (D-\mathit{\phi })\\ {\displaystyle }& {\displaystyle }+\cos \left(I\right)\cdot \mathit{\delta }\cdot \sin (D-\mathit{\phi }))+S_{\mathrm{off}}.\end{array}$$

We used the latter, because it facilitated the calculation of the partial derivatives for $D,I,\mathit{\delta }$ and ϵ, which are needed in the Gauss–Newton method. The factor c is the scale factor of the fluxgate magnetometer. It can be c=1 or $c=-\mathrm{1}$, depending on the orientation of the sensor (c=1 if S>0 and the telescope pointing towards north). A slight deviation of $\left|c\right|$ from 1 can be neglected, as long as S is small (Marsal and Torta, 2007).

3.1.2 Accounting for field changes during the measurement by reduction of the first

For each measurement at time t_i the variations of the magnetic field after the first measurement t₁ have to be taken into account. This reduction of the first measurement is carried out by projecting the variation in the direction of the magnetometer sensor of the ith measurement. Changes in D and I after the first measurement are assumed small enough, justifying the assumption that the instrument function f is linear in D₁ and I₁:

$$\begin{array}{ll}\text{(4)}& {\displaystyle }{S}_i& {\displaystyle }=f(D_{\mathrm{1}}+\mathrm{\Delta }{D}_i,I_{\mathrm{1}}+\mathrm{\Delta }{I}_i,\mathit{\delta },\mathit{ϵ},S_{\mathrm{off}},{\mathit{\phi }}_i,{\mathit{\xi }}_i){\displaystyle }& {\displaystyle }\approx f(D_{\mathrm{1}},I_{\mathrm{1}},\mathit{\delta },\mathit{ϵ},S_{\mathrm{off}},{\mathit{\phi }}_i,{\mathit{\xi }}_i)+{\displaystyle \frac{\partial f}{\partial D_{\mathrm{1}}}}\mathrm{\Delta }{D}_i+{\displaystyle \frac{\partial f}{\partial I_{\mathrm{1}}}}\mathrm{\Delta }{I}_i\\ {\displaystyle }& {\displaystyle }\text{with}\mathrm{\Delta }{D}_i=D_i-D_{\mathrm{1}}\text{and}\mathrm{\Delta }{I}_i=I_i-I_{\mathrm{1}}.\end{array}$$

3.1.3 Declination, inclination and instrument parameters by solution of a system of non-linear conditional equations

Each measurement with the DI-flux theodolite consists of measuring the angles φ_i, ξ_i and the magnetometer reading S_i. It delivers one non-linear equation for the unknown quantities D₁, I₁ and the instrument parameters δ, ϵ and S_off. Let $\mathit{p}=(D_{\mathrm{1}},I_{\mathrm{1}},\mathit{\delta },\mathit{ϵ})$ be a vector of unknown parameters. The magnetometer offset S_off is not included. It will be treated separately as explained in Sect. 3.1.5. For a total number of N available measurements we get a system of N equations. Generally the number of equations exceeds the number of unknowns. Hence the resulting non-linear system is overdetermined. A solution of this system exists only in the sense of the least square solution minimizing the sum of squares of residuals r_i:

$$\begin{array}{ll}\text{(5)}& {\displaystyle }& {\displaystyle }{S}_{\mathrm{1}}-f(\mathit{p},S_{\mathrm{off}},{\mathit{\phi }}_{\mathrm{1}},{\mathit{\xi }}_{\mathrm{1}})=r_{\mathrm{1}}{\displaystyle }& {\displaystyle }\mathrm{\dots }\\ {\displaystyle }& {\displaystyle }{S}_N-f(\mathit{p},S_{\mathrm{off}},{\mathit{\phi }}_N,{\mathit{\xi }}_N)=r_N,\end{array}$$

with

$$S_i=S_{\mathrm{reading},i}-{\displaystyle \frac{\partial f}{\partial D_{\mathrm{1}}}}\mathrm{\Delta }{D}_i-{\displaystyle \frac{\partial f}{\partial I_{\mathrm{1}}}}\mathrm{\Delta }{I}_i.$$

The magnetometer readings reduced in the first reading, and ΔD_i , ΔI_i variations of D and I since the first reading. Additional conditional equations can be included by seamlessly adding additional equations. Such equations can be deduced from a priori information. A priori information stabilizes the solution, especially if the number of available measurements are small. If, for example, the collimation angles are known from former measurements (δ_apriory and ϵ_apriory), and if their accuracy can be estimated in terms of standard deviations ${\mathit{\sigma }}_{{\mathit{\delta }}_{\mathrm{apriory}}}$ and ${\mathit{\sigma }}_{{\mathit{ϵ}}_{\mathrm{apriory}}}$, these a priori information can be accounted for using the following additional equations:

$$\begin{array}{ll}\text{(6)}& {\displaystyle \frac{{\mathit{\sigma }}_S}{{\mathit{\sigma }}_{{\mathit{\delta }}_{\mathrm{apriory}}}}}(\mathit{\delta }-{\mathit{\delta }}_{\mathrm{apriory}})& {\displaystyle }=r_{N+\mathrm{1}}\text{and}{\displaystyle \frac{{\mathit{\sigma }}_S}{{\mathit{\sigma }}_{{\mathit{ϵ}}_{\mathrm{apriory}}}}}(\mathit{ϵ}-{\mathit{ϵ}}_{\mathrm{apriory}})& {\displaystyle }=r_{N+\mathrm{2}}.\end{array}$$

The factors $\frac{{\mathit{\sigma }}_S}{{\mathit{\sigma }}_{{\mathit{ϵ}}_{apriory}}}$ and $\frac{{\mathit{\sigma }}_S}{{\mathit{\sigma }}_{{\mathit{\delta }}_{\mathrm{apriory}}}}$ are indispensable for the correct weighting of the a priori information. An estimate of the deviation σ_S of the magnetometer readings S_i can be determined by empirically investigating the distribution of the residuals r_i.

3.1.4 The first estimate of D and I

In order to converge towards the proper solution, the Gauss–Newton method needs a good first estimate. Instrument parameters δ, ϵ and S_off are small and can always be assumed to be zero as a first estimate.

The orientation vectors P_i of the telescope at each measurement are perpendicular to the magnetic field (neglecting the fluxgate readings S_i). A good first estimate of declination and inclination can be calculated by fitting a plane to these vectors

$$\begin{array}{}\text{(7)}& {\mathit{P}}_i=\left(\begin{array}{c}\cos \left({\mathit{\phi }}_i\right)\sin \left({\mathit{\xi }}_i\right)\\ \sin \left({\mathit{\phi }}_i\right)\sin \left({\mathit{\xi }}_i\right)\\ -\cos \left({\mathit{\xi }}_i\right)\end{array}\right);i=\mathrm{1}\mathrm{\dots }N.\end{array}$$

The normal vector $\widehat{\mathit{n}}$ on the least square plane fit through this set of points P_i gives a first but already pretty accurate approximation of the field direction. I is the angle between $\widehat{\mathit{n}}$ and its projection to the horizontal plane. D is the angle between the projection to the horizontal plane and north.

3.1.5 Solution of the non-linear system with the Gauss–Newton method

The Gauss–Newton method is a straightforward extension of the well-known Newton method from one to more dimensions. It is also referred to as the Newton–Raphson method (Press et al., 1986). The method works by iteratively improving the first estimate using a gradient method. The method requires calculation of the partial derivatives of the instrument model for all unknowns. It is known to have quadratic converge, which is very fast. However, the convergence is only guaranteed for an appropriate first estimate.

We are looking for the improvement Δp to the estimated p so that $f(\mathit{p}+\mathrm{\Delta }\mathit{p})=S$, S being the magnetometer readings. Linearization of f around p with the Jacobian matrix J(p) leads to

$$\mathbf{J}\left(\mathit{p}\right)\mathrm{\Delta }\mathit{p}+f\left(\mathit{p}\right)=S.$$

With $r=S-f\left(\mathit{p}\right)$ as the residuals between the readings S and the results of the current model f(p) we get

$$\mathbf{J}\left(\mathit{p}\right)\mathrm{\Delta }\mathit{p}=r,$$

or by writing the Jacobian matrix J(p), Δp and r explicitly we obtain

$$\begin{array}{}\text{(8)}& \left(\begin{array}{cccc}\frac{\partial f_{\mathrm{1}}}{\partial D}& \frac{\partial f_{\mathrm{1}}}{\partial I}& \frac{\partial f_{\mathrm{1}}}{\partial \mathit{\delta }}& \frac{\partial f_{\mathrm{1}}}{\partial \mathit{ϵ}}\\ & \mathrm{\dots }& & \\ \frac{\partial f_N}{\partial D}& \frac{\partial f_N}{\partial I}& \frac{\partial f_N}{\partial \mathit{\delta }}& \frac{\partial f_N}{\partial \mathit{ϵ}}\end{array}\right)\left(\begin{array}{c}\mathrm{\Delta }D\\ \mathrm{\Delta }I\\ \mathrm{\Delta }\mathit{\delta }\\ \mathrm{\Delta }\mathit{ϵ}\end{array}\right)=\left(\begin{array}{c}{r}_{\mathrm{1}}\\ \mathrm{\dots }\\ r_N\end{array}\right),\end{array}$$

with N as the number of available measurements.

The function f_i(p) is the instrument model $f_i\left(\mathit{p}\right)=f(D,I,\mathit{\delta },\mathit{ϵ},S_{\mathrm{off}},{\mathit{\phi }}_i,{\mathit{\xi }}_i)$ of Eq. (3) with the measured values φ_i and ξ_i.

In our case Eq. (8) is overdetermined and must be solved in the sense of minimizing $\left|r\right|$. The normal equations reads

$$\begin{array}{}\text{(9)}& \mathrm{\Delta }\mathit{p}={\left(\mathbf{J}\left(p{)}^{\mathrm{T}}\mathbf{J}\right(p)\right)}^{-\mathrm{1}}\mathbf{J}(p{)}^{\mathrm{T}}r.\end{array}$$

In Eq. (8) the instrument offset S_off does not show up. For the sake of simplicity and numerical stability we separated its calculation from the other parameters. S_off is just an additive term in all equations as shown in Eq. (3). Accordingly it can be calculated after each iteration as the mean value of the residuals.

3.1.6 Calculating the uncertainty of $D,I,\mathit{\delta }$ and ϵ in terms of standard deviation σ

The linear relation between Δp and r given in Eq. (9) makes it easy to propagate the uncertainty of r to the results p_i. We assume that the uncertainty of the parameter vector p is equal to the uncertainty of the last improvement vector Δp within the Gauss–Newton iteration. If enough measurements are available, the uncertainty of r can be determined to investigate the distribution around its mean value. Otherwise the reading error of S has to be estimated. This approach may lack mathematical rigour, but we verified the result using numerical Monte Carlo experiments. We rewrite Eq. (9) with a single matrix G as

$$\begin{array}{}\text{(10)}& \mathrm{\Delta }\mathit{p}=\mathbf{G}\cdot r\mathrm{with}\mathbf{G}=\left(g_{j,i}\right)={\left(\mathbf{J}\left(p{)}^{\mathrm{T}}\mathbf{J}\right(p)\right)}^{-\mathrm{1}}\mathbf{J}(p{)}^{\mathrm{T}},\end{array}$$

or

$$\begin{array}{}\text{(11)}& \mathrm{\Delta }{p}_j=\sum _i^N{g}_{j,i}{r}_i.\end{array}$$

Let the r_i by uncorrelated and normally distributed with the variance $V\left(r_i\right)={\mathit{\sigma }}_r^{\mathrm{2}}$. Then the p_i are also normally distributed and their variance can be calculated as

$$\begin{array}{}\text{(12)}& V\left(\mathrm{\Delta }{p}_i\right)=\sum _i^N{g}_{j,i}^{\mathrm{2}}{\mathit{\sigma }}_r^{\mathrm{2}}\mathrm{or}{\mathit{\sigma }}_{p_{\mathrm{1}}}={\mathit{\sigma }}_r\sqrt{\sum _i^N{g}_{j,i}^{\mathrm{2}}}.\end{array}$$

3.1.7 Baseline determination

Figure 2 (a) D_var and H_B+H_var form a right-angled triangle with H_abs as hypotenuse. (b) The variometer sensor is not perfectly aligned with the principal geographic directions. Therefore $D_{\mathrm{abs}}=D_B+D_{\mathrm{var}}$.

Absolute measurements are needed to determine the vector magnetometer offsets and their long-term drift. The dynamic range of the instrument is often adapted to the natural variations of the field, which is far smaller than the constant part. Hence the vector magnetometer is referred to as a “variometer”. A slight drift due to mechanic and electronic instabilities can never be excluded. A very small rotational movement, e.g. in the instrument or its pillar due to temperature, can result in a measurable effect. The variometer orientation in the geographical reference frame must be taken into account. Absolute measurements are indispensable for determining the vector of the geomagnetic field in the geographic reference frame over long periods of time. The wording “absolute” traces back to Gauss and signifies that absolute measurements are inherently adjusting for first-order instrument inexactness (e.g. misalignment of telescope and magnetometer sensor).

The offsets are called “base values”. The stability of the base values is indicative for the stability of the variometer, for the accuracy of the absolute measurements, and for unwanted small scale local magnetic field changes. Thus they are a measure of the quality of an observatory.

Baseline formulas: traditionally baselines refer to H, D and Z (not to the field components X, Y and Z). We assume that the variometer is set up with its X sensor almost pointing towards magnetic north, the Z sensor downwards and consequently the Y sensor eastward. Small misdirections are covered by the base values. For a better comprehensibility, we will subsequently use the notation H_abs, D_abs and Z_abs instead of H₁, D₁ and Z₁ in this framework. Baseline values are calculated by comparing H_abs, D_abs and Z_abs in the geographic coordinate system, as they result from the absolute measurement, to H_var, D_var and Z_var as “seen” by the variometer in its inherent coordinate system (X_var, Y_var, Z_var).

We use the fact that scalar values F_abs and H_abs are the same in the sensor system and in the geographic system. The base values H_B, D_B and Z_B are calculated as the differences between H_abs, D_abs and Z_abs, stemming from the absolute measurements, and H_var, D_var and Z_var, deduced from variometer measurements. With $H_{\mathrm{abs}}=F\cdot \cos \left(I\right)$ it is (see Fig. 2)

$$\begin{array}{ll}\text{(13)}& {\displaystyle }& {\displaystyle }(H_{\mathrm{B}}+X_{\mathrm{var}}{)}^{\mathrm{2}}+Y_{\mathrm{var}}^{\mathrm{2}}=H_{\mathrm{abs}}^{\mathrm{2}}\mathrm{or}{\displaystyle }& {\displaystyle }{H}_{\mathrm{B}}=\sqrt{H_{\mathrm{abs}}^{\mathrm{2}}-Y_{\mathrm{var}}^{\mathrm{2}}}-X_{\mathrm{var}},\end{array}$$

$$\begin{array}{ll}\text{(14)}& {\displaystyle }& {\displaystyle }{D}_B=D_{\mathrm{abs}}-D_{\mathrm{var}}=D_{\mathrm{abs}}-\arcsin \left({\displaystyle \frac{Y_{\mathrm{var}}}{H_{\mathrm{abs}}}}\right){\displaystyle }& {\displaystyle }\cong D_{\mathrm{abs}}-{\displaystyle \frac{\mathrm{360}{}^{\circ }}{\mathrm{2}\mathit{\pi }\cdot H_{\mathrm{abs}}}}{Y}_{\mathrm{var}}\end{array}$$

and finally

$$\begin{array}{}\text{(15)}& Z_B=Z_{\mathrm{abs}}-Z_{\mathrm{var}}=Z_{\mathrm{abs}}-F\cdot \sin \left(I_{\mathrm{abs}}\right).\end{array}$$

4 Results

4.1 Application of the method for outlier detection in traditional data sets

The new method data errors, otherwise undetected, can be identified and discarded, because data from each measured orientation can be assessed. Figure 3a shows H base values obtained with the new method (red) in comparison with H base values obtained conventionally (blue). Two measurements within about half a year of the observatory routine could be improved. Typing errors in the respective data sheets could be identified and corrected. They were detected because they produced residuals off the normal range (outliers). The corrected base values are flagged with red circles in Fig. 3a. Identification of outliers is illustrated in Fig. 3b. The residuals r_i are plotted as a scatter plot over the magnetometer readings S_reading,i as introduced in Eq. (5). The outlier can easily be detected. Outlier identification is more successful when more measurements are taken. If only a restricted number of measurements are available, a priori information should be used. This can be instrument parameters known from other measurements. After identifying an erroneous measurement, it is often possible to find the according typo in the data sheet. If this is not the case, the according measurement can simply be discarded. The spread in S_i direction in Fig. 3b depends on the measuring method used and gives no information on data quality.

Figure 3 Panel (a) shows conventionally measured baseline values of the H component in NGK in blue. The red points show results obtained with the new method. Calculated error bars allow their quality to be assessed. Red circles show baseline data which could be improved, because typos in the data sheets could be identified and corrected. Panel (b) gives a scatter plot of residuals r_i over magnetometer readings S_reading,i, clearly showing an outlier.

Figure 4 Panel (a) shows unsystematically increased residuals of a certain measurement. This indicates magnetic cleanliness problems with the operator. Panel (b) shows residuals which have a sinusoidal dependance of the pointing direction of the alidade ($\sim \sin (\mathrm{2}\cdot \mathit{\phi }+\mathit{\alpha })$). This indicates a systematic problem with the instrument. Measurements with the sensor above the telescope are coloured in blue; else they are in red. The green dots show the differences with respect to the fitted $\sin (\mathrm{2}\cdot \mathit{\phi }+\mathit{\alpha })$ curve. No sin(φ+α) dependency is observed.

Calculated error bars allow the reliability of each DI measurement to be assessed. An uncertainty of B_H less then 0.5nT can be expected. This is the case for the major part of the red data points shown in Fig. 3a. Furthermore the distribution of the residuals r_i allows the quality of a measurement to be assessed. Figure 3b gives a good example. Undisturbed measurements are characterized by residuals ranging between about −1nT and 1nT with a nice Gaussian distribution. In Fig. 3b the deviation of r_i would be in the expected range if the outlier was discarded. Figure 3a shows that measurements taken in the first half of April are slightly less precise than others. In this case it is worthwhile taking a closer look at the distribution of the residuals. Figure 4a shows residuals r observed in two routine measurements on 5 April (each measurement taken at 16 orientations). The observed residuals are clearly bigger than normal. The same behaviour can be observed on 12 April (Fig. 3a). We assume a problem in the magnetic cleanliness of the operating staff in the first half of April only.

A systematic error of a specific theodolite (Zeiss Theo 020, 817992) is revealed in Fig. 4b. In this figure the residuals are plotted over the reading of the horizontal circle φ. They are proportional to the sine of twice the horizontal reading $\sin (\mathrm{2}\cdot \mathit{\phi }+\mathit{\alpha })$. In Fig. 4b measurements with the sensor above the telescope are plotted in blue. The rest are plotted in red. The green dots show the differences with respect to the fitted $\sin (\mathrm{2}\cdot \mathit{\phi }+\mathit{\alpha })$ curve. They are plotted to check whether there is a sin(φ+α) dependency, which would indicate a mechanical problem with the theodolite. We found exactly the same behaviour when re-evaluating data taken 1 year earlier with the same theodolite. We could produce a similar effect by sticking magnetically soft material to the alidade of a good theodolite.

4.1.1 Avoiding problematic orientations close to the magnet equator

In areas of the globe with a small magnetic inclination, i.e. typically within 2000 km to the north and south of the magnetic equator, the conventional DI-flux procedure involves vertical circle readings at steep telescope orientations. This is not possible without zenith oculars mounted on the theodolite. The new scheme circumvents this problem if the tilt angle α (see former section) is sufficiently smaller than the inclination.

5 Discussion

We have been testing the new method for more than 1 year at the Niemegk Observatory. A first test was applying the new numeric evaluation to data produced with the conventional scheme. In case of faultless data, we got exactly the same result but this time with error bars. We also evaluated partly corrupted data, which could not be treated conventionally. The method can still be used if one or even two out of eight measurements have to be discarded. Investigation of the residuals often allowed us to identify and correct typos in the data sheet. A data set with the vertical reading set to 90 Gon instead of 100 Gon, produced by a person who was used to a degree scale, could be evaluated without problems. Furthermore we found that the major advantage in observatory routine is that several data sets measured on the same day can be evaluated at once as a single data set. At observatories it is good practice to repeat measurements several times. Evaluating them at once leads to a statistically firmer result and facilitates the identification of outliers. Investigating the calculated error bar allows measuring conditions to be assessed. As shown in Fig. 3a we could identify a problem in magnetic cleanliness in the first half of April. Figure 4a shows the respective residuals. Measuring in more than the four principal directions (see Fig. 4b) reveals properties of the DI theodolite, which otherwise remain hidden. This gives a new perspective to assess a DI-flux theodolite. A simpler and less error-prone measuring scheme was developed.