Journal cover Journal topic
Geoscientific Instrumentation, Methods and Data Systems An interactive open-access journal of the European Geosciences Union
Journal topic
Geosci. Instrum. Method. Data Syst., 7, 1-9, 2018
https://doi.org/10.5194/gi-7-1-2018
Geosci. Instrum. Method. Data Syst., 7, 1-9, 2018
https://doi.org/10.5194/gi-7-1-2018

Research article 23 Jan 2018

Research article | 23 Jan 2018

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

Numerical evaluation
Heinz-Peter Brunke and Jürgen Matzka Heinz-Peter Brunke and Jürgen Matzka
• GFZ German Research Centre for Geosciences, Telegrafenberg, 14473 Potsdam, Germany
Abstract

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 (Brunke2017).

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 and Anna Willer (née Nilsson) (Nilsson2010) 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 .

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 D1 and inclination I1 are principally restricted to calculating mean values; see .

Declination and inclination usually vary during the measuring process. Let D1 and I1 be the values of declination and inclination at the time of the first magnetometer reading S1. 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 ΔDi (reduction to the first measurement). The natural variation ΔDi from the first measurement D1 are known from the variometer. Additionally a correction ΔAi is needed if the magnetometer reading Si is not exactly zero . Assuming subsequently that D1, I1 and all measured angles are given in degrees, the corrected horizontal circle readings are

$\begin{array}{ll}& {A}_{{X}_{i}}={\mathit{\phi }}_{i}-\mathrm{\Delta }{D}_{i}-\mathrm{\Delta }{A}_{i}\\ & ={\mathit{\phi }}_{i}-\mathrm{\Delta }{D}_{i}-\frac{\mathrm{360}{}^{\circ }}{\mathrm{2}\mathit{\pi }}\frac{{S}_{i}}{H}\phantom{\rule{1em}{0ex}}\mathrm{with}\\ & {X}_{i}\in \mathit{\left\{}W↑,W↓,E↑,E↓\mathit{\right\}};\phantom{\rule{1em}{0ex}}i=\mathrm{1}\phantom{\rule{0.25em}{0ex}}\mathrm{\dots }\phantom{\rule{0.25em}{0ex}}\mathrm{4}.\end{array}$

The ${X}_{i}\in \mathit{\left\{}W↑,E↑,W↓,E↓\mathit{\right\}}$ 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 Si, ideally to zero. Eventually D1 is determined as the mean value of ${A}_{{X}_{i}}$ plus 180.

$\begin{array}{}\text{(1)}& {D}_{\mathrm{1}}=\frac{\mathrm{1}}{\mathrm{4}}\left({A}_{W↑}+{A}_{W↓}+{A}_{E↑}+{A}_{E↓}\right)+\mathrm{180}{}^{\circ }.\end{array}$

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

$\begin{array}{}\text{(2)}& {I}_{\mathrm{1}}=\frac{\mathrm{1}}{\mathrm{4}}\left({V}_{N↑}+{V}_{S↓}-{V}_{S↑}-{V}_{N↓}\right)±\mathrm{90}{}^{\circ }.\end{array}$

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

$\begin{array}{ll}{V}_{{X}_{i}}& ={\mathit{\xi }}_{i}-\mathrm{\Delta }{I}_{i}-\mathrm{\Delta }{V}_{i};\phantom{\rule{1em}{0ex}}{X}_{i}\in \mathit{\left\{}N↑,S↓,S↑,N↓\mathit{\right\}}\\ & ={\mathit{\xi }}_{i}-\mathrm{\Delta }{I}_{i}-\frac{\mathrm{360}{}^{\circ }}{\mathrm{2}\mathit{\pi }}\frac{{S}_{i}}{F}\phantom{\rule{1em}{0ex}};\phantom{\rule{1em}{0ex}}i=\mathrm{5}\phantom{\rule{0.25em}{0ex}}\mathrm{\dots }\phantom{\rule{0.25em}{0ex}}\mathrm{8}\end{array}$

Just as above, ΔIi are corrections for field variations after the first measurement. The ΔVi values account for non-zero magnetometer readings Si. 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 1Graphic 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 $\left(\mathrm{90}{}^{\circ }-\mathit{\alpha }\right)$ with the telescope pointing towards north-north-east. Adjust S horizontally.

• 6.

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

• a.

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

• b.

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

• 7.

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

• 8.

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

• a.

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

• b.

Optional: invert telescope, set ξ to $\left(\mathrm{270}{}^{\circ }+\mathit{\alpha }\right)$ 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 D1, inclination I1, the two angles of misalignment δ and ϵ and the magnetometer offset Soff.

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 Soff. The angles δ and ϵ (misalignment in φ and ξ direction) are also called aberration errors. Crosthwaite gave a formulation based on vector geometry. The following formulation and its derivation based on spherical trigonometry can be found in (Nilsson2010):

$\begin{array}{ll}\text{(3)}& S& =f\left(D,I,\mathit{\delta },\mathit{ϵ},{S}_{\mathrm{off}},\mathit{\phi },\mathit{\xi }\right)& =c\cdot F\cdot \left(-\mathrm{sin}\left(I\right)\cdot \mathrm{cos}\left(\mathit{\xi }+\mathit{ϵ}\right)\\ & \phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}+\mathrm{cos}\left(I\right)\cdot \mathrm{sin}\left(\mathit{\xi }+\mathit{ϵ}\right)\phantom{\rule{0.125em}{0ex}}\cdot \mathrm{cos}\left(D-\mathit{\phi }\right)\\ & \phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}+\mathrm{cos}\left(I\right)\cdot \phantom{\rule{2em}{0ex}}\mathit{\delta }\phantom{\rule{2em}{0ex}}\cdot \mathrm{sin}\left(D-\mathit{\phi }\right)\phantom{\rule{0.25em}{0ex}}\right)\phantom{\rule{0.25em}{0ex}}+{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 $|c|$ from 1 can be neglected, as long as S is small .

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

For each measurement at time ti the variations of the magnetic field after the first measurement t1 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 D1 and I1:

$\begin{array}{ll}\text{(4)}& {S}_{i}& =f\left({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}\right)& \approx f\left({D}_{\mathrm{1}},{I}_{\mathrm{1}},\mathit{\delta },\mathit{ϵ},{S}_{\mathrm{off}},{\mathit{\phi }}_{i},{\mathit{\xi }}_{i}\right)+\frac{\partial f}{\partial {D}_{\mathrm{1}}}\mathrm{\Delta }{D}_{i}+\frac{\partial f}{\partial {I}_{\mathrm{1}}}\mathrm{\Delta }{I}_{i}\\ & \phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}\phantom{\rule{1em}{0ex}}\text{with}\phantom{\rule{1em}{0ex}}\mathrm{\Delta }{D}_{i}={D}_{i}-{D}_{\mathrm{1}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\text{and}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\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 Si. It delivers one non-linear equation for the unknown quantities D1, I1 and the instrument parameters δ, ϵ and Soff. Let $\mathbit{p}=\left({D}_{\mathrm{1}},{I}_{\mathrm{1}},\mathit{\delta },\mathit{ϵ}\right)$ be a vector of unknown parameters. The magnetometer offset Soff 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 ri:

$\begin{array}{ll}\text{(5)}& & {S}_{\mathrm{1}}-f\left(\mathbit{p},{S}_{\mathrm{off}},{\mathit{\phi }}_{\mathrm{1}},{\mathit{\xi }}_{\mathrm{1}}\right)={r}_{\mathrm{1}}& \mathrm{\dots }\\ & {S}_{N}-f\left(\mathbit{p},{S}_{\mathrm{off}},{\mathit{\phi }}_{N},{\mathit{\xi }}_{N}\right)={r}_{N},\end{array}$

with

${S}_{i}={S}_{\mathrm{reading},i}-\frac{\partial f}{\partial {D}_{\mathrm{1}}}\mathrm{\Delta }{D}_{i}-\frac{\partial f}{\partial {I}_{\mathrm{1}}}\mathrm{\Delta }{I}_{i}.$

The magnetometer readings reduced in the first reading, and ΔDi , ΔIi 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)}& \frac{{\mathit{\sigma }}_{S}}{{\mathit{\sigma }}_{{\mathit{\delta }}_{\mathrm{apriory}}}}\left(\mathbit{\delta }-{\mathit{\delta }}_{\mathrm{apriory}}\right)& ={r}_{N+\mathrm{1}}\phantom{\rule{1em}{0ex}}\text{and}\frac{{\mathit{\sigma }}_{S}}{{\mathit{\sigma }}_{{\mathit{ϵ}}_{\mathrm{apriory}}}}\left(\mathbit{ϵ}-{\mathit{ϵ}}_{\mathrm{apriory}}\right)& ={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 Si can be determined by empirically investigating the distribution of the residuals ri.

### 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 Soff are small and can always be assumed to be zero as a first estimate.

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

$\begin{array}{}\text{(7)}& {\mathbit{P}}_{i}=\left(\begin{array}{c}\mathrm{cos}\left({\mathit{\phi }}_{i}\right)\mathrm{sin}\left({\mathit{\xi }}_{i}\right)\\ \mathrm{sin}\left({\mathit{\phi }}_{i}\right)\mathrm{sin}\left({\mathit{\xi }}_{i}\right)\\ -\mathrm{cos}\left({\mathit{\xi }}_{i}\right)\end{array}\right);\phantom{\rule{1em}{0ex}}i=\mathrm{1}\mathrm{\dots }N.\end{array}$

The normal vector $\stackrel{\mathrm{^}}{\mathbit{n}}$ on the least square plane fit through this set of points Pi gives a first but already pretty accurate approximation of the field direction. I is the angle between $\stackrel{\mathrm{^}}{\mathbit{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 . 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\left(\mathbit{p}+\mathrm{\Delta }\mathbit{p}\right)=S$, S being the magnetometer readings. Linearization of f around p with the Jacobian matrix J(p) leads to

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

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

$\mathbf{J}\left(\mathbit{p}\right)\mathrm{\Delta }\mathbit{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 fi(p) is the instrument model ${f}_{i}\left(\mathbit{p}\right)=f\left(D,I,\mathit{\delta },\mathit{ϵ},{S}_{\mathrm{off}},{\mathit{\phi }}_{i},{\mathit{\xi }}_{i}\right)$ 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 $|r|$. The normal equations reads

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

In Eq. (8) the instrument offset Soff does not show up. For the sake of simplicity and numerical stability we separated its calculation from the other parameters. Soff 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 pi. 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 }\mathbit{p}=\mathbf{G}\cdot r\phantom{\rule{0.25em}{0ex}}\mathrm{with}\phantom{\rule{0.25em}{0ex}}\mathbf{G}=\left({g}_{j,i}\right)={\left(\mathbf{J}\left(p{\right)}^{\mathrm{T}}\mathbf{J}\left(p\right)\right)}^{-\mathrm{1}}\mathbf{J}\left(p{\right)}^{\mathrm{T}},\end{array}$

or

$\begin{array}{}\text{(11)}& \mathrm{\Delta }{p}_{j}=\sum _{i}^{N}{g}_{j,i}{r}_{i}\phantom{\rule{1em}{0ex}}.\end{array}$

Let the ri by uncorrelated and normally distributed with the variance $V\left({r}_{i}\right)={\mathit{\sigma }}_{r}^{\mathrm{2}}$. Then the pi 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}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\mathrm{or}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}{\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) Dvar and HB+Hvar form a right-angled triangle with Habs 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 Habs, Dabs and Zabs instead of H1, D1 and Z1 in this framework. Baseline values are calculated by comparing Habs, Dabs and Zabs in the geographic coordinate system, as they result from the absolute measurement, to Hvar, Dvar and Zvar as “seen” by the variometer in its inherent coordinate system (Xvar, Yvar, Zvar).

We use the fact that scalar values Fabs and Habs are the same in the sensor system and in the geographic system. The base values HB, DB and ZB are calculated as the differences between Habs, Dabs and Zabs, stemming from the absolute measurements, and Hvar, Dvar and Zvar, deduced from variometer measurements. With ${H}_{\mathrm{abs}}=F\cdot \mathrm{cos}\left(I\right)$ it is (see Fig. 2)

$\begin{array}{ll}\text{(13)}& & \left({H}_{\mathrm{B}}+{X}_{\mathrm{var}}{\right)}^{\mathrm{2}}+{Y}_{\mathrm{var}}^{\mathrm{2}}={H}_{\mathrm{abs}}^{\mathrm{2}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\mathrm{or}& {H}_{\mathrm{B}}=\sqrt{{H}_{\mathrm{abs}}^{\mathrm{2}}-{Y}_{\mathrm{var}}^{\mathrm{2}}}-{X}_{\mathrm{var}},\end{array}$

$\begin{array}{ll}\text{(14)}& & {D}_{B}={D}_{\mathrm{abs}}-{D}_{\mathrm{var}}={D}_{\mathrm{abs}}-\mathrm{arcsin}\left(\frac{{Y}_{\mathrm{var}}}{{H}_{\mathrm{abs}}}\right)& \cong {D}_{\mathrm{abs}}-\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 \mathrm{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 ri are plotted as a scatter plot over the magnetometer readings Sreading,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 Si direction in Fig. 3b depends on the measuring method used and gives no information on data quality.

Figure 3Panel (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 ri over magnetometer readings Sreading,i, clearly showing an outlier.

Figure 4Panel (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 \mathrm{sin}\left(\mathrm{2}\cdot \mathit{\phi }+\mathit{\alpha }\right)$). 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 $\mathrm{sin}\left(\mathrm{2}\cdot \mathit{\phi }+\mathit{\alpha }\right)$ curve. No sin⁡(φ+α) dependency is observed.

Calculated error bars allow the reliability of each DI measurement to be assessed. An uncertainty of BH 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 ri 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 ri 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 $\mathrm{sin}\left(\mathrm{2}\cdot \mathit{\phi }+\mathit{\alpha }\right)$. 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 $\mathrm{sin}\left(\mathrm{2}\cdot \mathit{\phi }+\mathit{\alpha }\right)$ 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.

Code and data availability
Code and data availability.

The MATLAB source code implementing the method presented above is available from GFZ Data Services (Brunke2017). Test data sets are provided along with the source code.

Competing interests
Competing interests.

The authors declare that they have no conflict of interest.

Special issue statement
Special issue statement.

This article is part of the special issue “The Earth's magnetic field: measurements, data, and applications from ground observations (ANGEO/GI inter-journal SI)”. It is a result of the XVIIth IAGA Workshop on Geomagnetic Observatory Instruments, Data Acquisition and Processing, Dourbes, Belgium, 4–10 September 2016.

Acknowledgements
Acknowledgements.

We thank Anna Willer, Truls Lynne Hansen and Peter Crostwaite. They provided their experiences in numerically modelling a DI-flux theodolite and solving the needed non-linear system of equations.

We thank Christopher Turbitt for the editing. We thank two anonymous reviewers for an intense and fruitful reviewing effort. Kirsten Elger helped to publish the software as a digital object.

The article processing charges for this open-access
publication were covered by a Research
Centre of the Helmholtz Association.

Edited by: Christopher Turbitt
Reviewed by: two anonymous referees

References

Brunke, H.-P.: AlternateMagAbsolutes – A software package for Magnetic Absolute Measurements with arbitrarily distributed DI Theodolite Orientations V.1.0, GFZ Data Services, https://doi.org/10.5880/GFZ.2.3.2017.003, 2017. a, b

Crosthwaite, P.: Declination-Inclination Magnetometer: Theory and Practice, Australien Geological Survey Organisation (Geomagnetism Section) Geomagnetism Note 1994-11), 1994. a, b

Jankowski, J. and Sucksdorff, C.: IAGA Guide for Magnetic Measurements and Observatory Practice, IAGA, http://www.iaga-aiga.org/data/uploads/pdf/guides/iaga-guide-observatories.pdf, 1996. a, b

Kerridge, D.: The theory of the fluxgate theodolite, British Geological Survey Geomagnetic Research Group Report No. 88/14, 1988. a

KringLauridsen: Experiences with the DI-Fluxgate Magnetometer inclusive Theory of the Instrument and Comparison to other Methods, Danish Meteorological Institute, Geophysical Papers, 1985. a, b

Marsal, S. and Torta, J.: An evaluation of the uncertainty associated with the measurement of the geomagnetic field with a D/I fluxgate theod olite, Meas. Sci. Technol., 18, 2143–2156, 2007. a

Matzka, J. and Hansen, T. J.: On the Various Published Formulas to Determine Sensor Offset and Sensor Misalignment for the DI-flux, Publs. Inst. Geophys. Pol. Acad. Sc., C-99, 398, https://doi.org/10.1088/0957-0233/18/7/046, 2007. a

Nilsson, A. N., Hansen T. L., and Olsen, N. M. J.: Theory of the fluxgate variometer and absolute obserations with the fluxgate theodolite – including a new procedure for the theodolite observations, personal communication, 2010. a, b

Press, W., Flannery, B., Teukolsky, S., and Vetterling, W.: Numerical Recipes – The Art of Scientific Computing, Cambridge University Press, https://websites.pmc.ucsc.edu/~fnimmo/eart290c_17/NumericalRecipesinF77.pdf, 1986. a

Special issue