Contribution to solving the orientation problem for an automatic magnetic observatory

The problem of the absolute calibration of a vectorial (tri-axial) magnetometer is addressed with the objective that the apparatus, once calibrated, gives afterwards, for a few years, the absolute values of the three components of the geomagnetic field (say the Northern geographical component, Eastern component and vertical component) with an accuracy on the order of 1 nT. The calibration procedure comes down to measure the orientation in space of the three physical axes of the sensor or, in other words, the entries of the transfer matrix from the local geographical axes to these physical axes. Absolute calibration follows indeed an internal calibration which provides accurate values of the three scale factors corresponding to the three axes – and in addition their relative angles. The absolute calibration can be achieved through classical absolute measurements made with an independent equipment. It is shown – after an error analysis which is not trivial – that, while it is not possible to get the axes absolute orientations with a high accuracy, the assigned objective (absolute values of the Northern geographical component, Eastern component and vertical component, with an accuracy of the order of 1 nT) is nevertheless reachable; this is because in the time interval of interest the field to measure is not far from the field prevailing during the calibration process.


Introduction
The geomagnetic field is continuously measured in a network of magnetic observatories, which, however, has significant gaps in the remote areas and over the oceans.This uneven distribution is linked to the fact that currently it is not possible to operate fully automated observatories which do not require manual operation of any instrument.Already, some fifty years ago, Alldregde planned an automatic standard magnetic observatory (ASMO) - (Alldregde, 1960;Alldregde and Saldukas, 1964), i.e. a device providing at each time the absolute values of -say -the Northern geographical component, Eastern component and vertical component of the geomagnetic field, without extra independent absolute measurements, at least for a long enough timespan.This idea has remained in the geomagnetism community and, over the last years, attempts have been made to automate a DItheodolite (van Loo and Rasson, 2006;Rasson and Gonsette, 2011), a proton vector magnetometer (Auster et al., 2006), or to build a device which can perform discrete absolute measurements automatically (Auster et al., 2007).The idea to use absolute measurements for solving orientation problems has been used by Schott and Leroy (2001) when developing a DIDD magnetometer.
In a former paper (Gravrand et al., 2001), the question was addressed of the internal calibration of a vectorial (or tri-axial) magnetometer, such as the 4 He pumped magnetometer built by the Laboratoire d'Electronique et de Technique de l'Information (LETI) of the French Commissariat à l'Energie Atomique (CEA), or the Oersted space mission fluxgate magnetometer, built by the Institute for automatization of the Danish Technical University (DTU).The problem of the internal calibration is to determine the six ( 4 He magnetometer) or nine (fluxgate magnetometer) intrinsic parameters: the three scale values (one for each modulation coil physical axis in the 4 He case, or each fluxgate sensor axis in the DTU case), the three angles between the three physical axis in both cases, and the three offset values in the case of the fluxgate variometer.The calibration problem can be solved by rotating in space the triaxial magnetometer while making simultaneously absolute intensity measurements with a proton or optical pumping magnetometer (Olsen et al., 2003), and it is essentially linear (Gravrand et al., 2001).The calibration algorithm was developed for ground operation but can be -and has been -extended to in flight calibration.But this internal calibration is not enough if we want to install somewhere a genuine automatic magnetic observatory.In the following u 1 , u 2 , u 3 are the unit vectors corresponding to these north east down directions.
It was stated in Gravrand et al. (2001) that such an absolute calibration should not be too difficult in a place where standard absolute measurements can be performed.This is the question (of quite practical interest) that we address in the present paper.We show that the above statement is both valid and invalid, depending on the objective.We also show that only the determination of the scale factors provided by the internal calibration process are of interest for the absolute calibration (the accurate determination of all angles between physical axes is not necessary however provided by the internal calibration; see also Appendix).

The practical problem
Suppose we want to install an automatic observatory in some new place, say a remote island in the Pacific ocean.What is required is, to refresh Alldregde's statement, to obtain one minute absolute values of the field components X, Y , Z in, say, the north east vertical down frame, fitting INTER-MAGNET standards (see http:/www.intermagnet.org),without needing an observer to visit the place in the few years following the installation.
One first builds a pillar (the permanent pillar) in a location propitious to install the 4 He magnetometer, and an auxiliary pillar a few meters apart.The calibration process can start.The observer determines the differences X, Y , Z between the absolute values of X, Y , Z at the two pillars.This is classical observatory work, not negligible, but which can be completed in a few days using a DI-flux theodolite and a proton magnetometer; modern devices for determining azimuths are welcome.The magnetometer-variometer, as we call it, can now be installed on the permanent pillar (in fact after a non magnetic house has been built around it; we do not develop here these practical aspects).By construction, the unit vectors e 1 , e 2 , e 3 of the physical axes, or coil axes, of the apparatus are nearly orthogonal, and its installation on the pillar is generally made in such a way that e 1 is close to u 1 , e 2 close to u 2 and e 3 close to u 3 ; although this is by no way a necessary condition.The observer makes at the auxiliary pillar a series of absolute measurements of the magnetic field at time moments t 1 , t 2 , ...t k , ... and corrections X, Y , Z are applied to get the corresponding absolute values on the permanent pillar.At the same time moments t k , the magnetometer-variometer to be calibrated provides the values ..K of the components of the magnetic vector V (t k ) along its physical axes {e i } whose orientations with respect to {u 1 , u 2 , u 3 } are not exactly known.
The observer, with his equipment, now leaves the place.The magnetometer-variometer in place continues to provide the values V 1 (t), V 2 (t), V 3 (t) of the (contravariant) components of V along its physical axes.The problem to solve is the following: how, relying on the set of absolute measurements made previously at times t 1 , t 2 , ...t k , to compute the geographical components (X 1 (t), X 2 (t), X 3 (t)) of V (t) at any following time t (in fact depending on the sampling rate), and estimate the error on those computed values?
This error, as we will see it, is a direct function of the errors on the absolute measurements made at the times t 1 , t 2 , ... We call it the calibration error.Let us say that we note B(t k ) = B k the absolute measurements at time t k and V (t) the measurements provided by the magnetometer-variometer.

The principle of the calibration
To compute the e i vectors in the u i frame, we go through the B k .Obviously one B k is not enough; but as a linear operator in R 3 is uniquely defined by its action on three linearly independent vectors, we take three of them, that we note B 1 , B 2 , B 3 , to present the algorithm of the calibration.In practice, several triplets among the K measurements available, if K > 3, are used.
Let d j k be the geographical components of of B k (k = 1, 2, 3) (as measured by the observer), and f j k the values of the (contravariant) components of B k along the physical axes e 1 , e 2 , e 3 provided at the same times by the magnetometervariometer.We have where the hat symbol is to stress the error-free nature of the corresponding quantities.The solution of the calibration is trivial, the e i being straightforwardly obtained in function of the u i through the B k : where F and D are the matrices of coefficients (components) f j k and dj k .The magnetometer-variometer provides the values (V 1 (t), V 2 (t), V 3 (t)) of the physical components of V : The problem is solved in the error-free case; we have obtained the geographical components (X 1 , X 2 , X 3 ) of V , at each measurement time.But the problem is in dealing with errors.

General statements
The absolute measurements of B 1 , B 2 , B 3 are affected by errors which can be viewed as errors on the geographical components d j k of B k vectors (Eq.1).They could be discussed at some length; but it is a very well-known topic.For the purpose of the present study, we suppose the magnitude of the errors on B k to be ε in relative value and randomly distributed in direction.In other words: Here ε 1 and r j k are O(1).In matrix form: R is the matrix of the r j k , and b a characteristic value of the B intensity at the station and the epoch of interest (see infra).
The values of the components along the physical axes {e i }, f j k are also not error-free.Nevertheless, in the case of the 4 He magnetometer that we have especially in mind, these f j k are measured with a very high accuracy, better than 0.1 nT (Léger et al., 1992;Gravrand et al., 2001) after the internal calibration has been performed.To simplify the writing, we consider the values f j k as error-free, i.e. f j k = f j k .Indeed, the f j k , the components of B k along the physical axes of the variometer-magnetometer, are determined with a high accuracy: 0.1 nT, i.e. a relative error of a few 10 −6 .The poor conditioning of F matrix comes from the B 1 , B 2 , B 3 of the calibration being only δbω apart, with the values considered in this paper (|ω i | = 0.1, δb ≈ 50 nT), i.e. δ ≈ 10 −3 .Therefore, the measurement error of f does not significantly distort the inverse matrix F −1 .It is in fact possible to develop the theory taking into account errors on f (supposed larger than above).It makes the writing larger and heavier.We have prefered to present the simplified version, essentially relevant, in this paper.
It now comes immediately: where εbω i , i = 1, 2, 3 denote the error along the i-th direction, and Multiplying Eq. ( 4) by F −1 (recall that F −1 = F−1 ) and using Eq. ( 2): In other words, when computing the physical unit vectors e i using the "measured" transformation matrix C = F −1 D (instead of Ĉ = F −1 D), an error is made which depends on F −1 .The difficulty to be expected is rather obvious.We go from the orthogonal frame u i to the nearly tri-orthogonal frame e i through the B k frame.But the three vectors B 1 , B 2 , B 3 have directions close to one another (remember that they are measurements made at the station during a timespan of say a week; see Sect. 4 for numerical values).The matrix F whose lines are close to one another is a priori poorly conditioned; its inverse F −1 may have large eigenvalues, and a strong amplification of error εb might affect the directions of e i , and the error on V (t) might be much larger than the error εb on B k (see Appendix).But, in fact, the practical conditions of the calibration process (B k ) and of the following measurements of the current magnetic field V (t) by the magnetometer-variometer discard such error amplification as shown later.We now build a simple algorithm allowing a statistical modeling and providing realistic error estimation, sufficient for the present study.

A simple algorithm
Let us now consider the vector V at time t.From Eqs. (1), (3) and (5): The B k are the true values, the B k the erroneous absolute measurements of the field at times t k .V (t) is the true value of V (t) at time t and V (t) the erroneous measurement of V (t) www.geosci-instrum-method-data-syst.net/2/1/2013/Geosci.Instrum.Method.Data Syst., 2, 1-9, 2013 provided by the magnetometer-variometer due to the error on the determination of the physical axes directions e i .The calibration error on V (t) appears directly as a linear form of the measurement errors on the B k , without explicit reference to the frames u i and e i : Of course we do not know the true values of B 1 , B 2 , B 3 , but to estimate the error | V | = |V − V |, it is enough to replace in Eq. ( 8) the quantities (B i − B i ) by their estimates.For that, for a given triplet B k and a given vector V we resort to a statistical estimate.We consider that the measurement errors B i −B i are randomly and uniformly distributed in a ball of center B i and radius εb.The small balls are drawn in dark gray in Fig. 1.Note that the B i and B i can be interchanged, in this error estimation.

Numerical results
In this section we estimate the calibration error using the simple algorithm presented above, for different configurations of V (t) and B k , k = 1, 2, 3.These vectors are essentially taken or simulated from observatory records.In other words, we estimate the calibration error which would affect the vector data provided by the variometer-magnetometer, in different locations at the Earth's surface.

Observatory data
We use one-minute values of the three components of the field recorded during the year 1999, as available on the IN-TERMAGNET CDROM 1999, from four observatories: a high-latitude observatory, Resolute-Bay (RES), an equatorial observatory, Bangui (BNG), and two middle-latitude observatories, one in the Northern and one in the Southern hemispheres, Chambon la Forêt (CLF) and Hermanus (HER).Their coordinates are given in Table 1.
From the theory developed above, it is obvious that the more diverse in direction the B k are, the better the configuration is for calibration.Then, at a given observatory, the larger the magnetic activity, the larger the probability for the triplet B k to be open.To evidence this effect, we select, in 1999, for each observatory O i , two subsets of 60 days each, Q i containing the five quietest days of each month of 1999, and D i containing the five most disturbed days.The day selection is made using the Kp indices (Mayaud, 1980).

Effect of the (B 1 , B 2 , B 3 ) configuration
To study this effect, we take a full day of one-minute values of X, Y , Z from, for example, CLF, specifically the day 6 September 1999, a quiet day belonging to Q CLF .We form at each minute t the triplet B 1 = B(t), B 2 = B(t + t 0 ), B 3 = B(t + 2 t 0 ).And, along the lines indicated supra, we associate to each of these triplets a set of vectors (B 1 , B 2 , B 3 ), B i being in the ball of center B i and radius εb (Fig. 1).
Note that we have (1441 − 2 t 0 ) triplets B k (t 0 in minutes).We then compute the calibration error -through formula Eq. ( 8) -affecting a set of vectors V = W + v, W being the mean value of the (recorded) field for day 6 September 1999, and v a vector uniformly distributed in a ball of center W and radius δb (the big light gray ball of Fig. 1); note that the set of vectors V is partly simulated.We compute, for a given calibration triplet B k (t), (Eq.9) the value |V − V | for all represent the measurement error εb; value δb is the upper bound for all ; b is a typical value of the intensity of the geomagnetic field at the site and epoch of measurements.We make this computation for all B k , t = 1, 2, ...1441 − 2t 0 .We order the collection | V | av and | V | max using a parameter η which grossly characterizes the quality of the configuration B k , i.e. the aperture of this triplet.We choose <> is for the vector triple product.Figure 3

Histograms of the calibration error
We now present some reciprocal numerical experiments, closer to the real situation to be met, using again minute data of day 6 September 1999.This time we choose a single absolute measurement triplet B k , k = 1, 2, 3, picked up in the observatory records, specifically at t = 0300, 0600, 1500 on 6 September 1999, and retain as current vectors V (t) all the one-minute values recorded at CLF over the 1999 year (instead of the simulated vectors in the ball of center W ). Again a set of triplets (B 1 , B 2 , B 3 ) is associated with (B 1 , B 2 , B 3 ), uniformly distributed in a ball of radius εb centered respectively at (B 1 , B 2 , B 3 ).realistic estimate), for εb = 1 nT, is most of the time smaller than 2 nT (Fig. 5).In Fig. 6 | V | max (t) values are simply ranked versus time t.An examination of this figure in regard of the magnetic situation shows that, as expected, the largest values of | V | max (t) are associated with magnetic storms.V (t), during these events, leaves the ball of centre W and radius δb (Fig. 1; δb = 50 nT).Note in passing that it is not important, in general, to know with a high accuracy the absolute value of V (t) at each minute of a magnetic storm.

Time tables
We now give, for each of our four observatories, a different, more practical presentation of the calibration error, which gives the timespans during which this error is smaller than a given threshold of αnT.We choose again values of year 1999, consider the triplets   | B| av (t) as explained in Sect.5.2.Figures 7-11 present the results for the four observatories.All graphs -that we call time tables -are to be read in the following way: the upper sub-panel (in blue) is for the 60 most disturbed days of the year, the lower one (in red) for the 60 quietest days.In each of the sub-panels the days are ranked as follows: the five quietest days of January, according to their calendar date are at the bottom of the lower sub-panel, the five quietest days of December at its top.The same for the upper sub-panel.In Figs.7-10 the value of α is 2 nT.In Fig. 11, relative to Bangui, a time table for α = 4 nT is also presented.Everywhere t 0 = 7.5 h for the left panels, t 0 = 6.0 h for the right panels.Of course, t < 0900 for t 0 = 7.5 h (24 − 2 × 7.5), and t < 1200 for t 0 = 6.0 h.In all the computations εb = 1 nT.We plot a characteristic function which is equal to zero at time t (white) if the triplet (B 1 = B(t), B 2 = B(t + t 0 ), B 3 = B(t + 2 t 0 )) leads to a | B| error >2 nT; otherwise, a colored tiret, red or blue, is drawn.Continuous red or blue time intervals are such that, for any first measurement with t in this interval leads to a calibration error smaller than 2 nT (or 4 nT) on V (t).
Looking at graphs of Figs.7-10, we remark, as expected, that it is easier to get time intervals with < 2 nT for the disturbed days than for the quiet days, and for a high latitude observatory (RES) than for an equatorial one (BNG).Figure 11 is an example of the effect of changing α.

Discussion and conclusion
Stability in absolute values, and particularly long-term stability -say up to a few years -used to be the most difficult requirement to fulfill in magnetic observatories.Let us adopt the standards of the INTERMAGNET program, which are up to now essentially intended to classical observatories with regular (generally weekly) man-made absolute measurements.The one minute magnetic field values provided by the magnetometer-variometer should be characterized by a resolution of 0.1 nT and a long-term stability of 5 nT yr −1 .
From the results of Sect.4, it appears that, after the calibration performed as in Sects. 2 and 3, the magnetometervariometer -as already said, we have especially in mind the LETI (CEA) apparatus -can function as an automatic observatory, fitting INTERMAGNET standards for a time-span of one to a few years, depending on the amplitude of the secular variation.A special study is necessary in the case of the highest latitude observatories.The necessity of a visiting the station every other year or so to renew the calibration is not so hard a constraint; in any case, such visits should be necessary for other purposes and checking of instruments and Fig. 7. Time table for RES data, with t 0 = 7.5 h (left panel) and t 0 = 6 h (right panel) periods for the B 1 (t).Note the distinction between subsets D i and Q i : blue lines for disturbed days, red lines for quiet days.Yearly variation reads from down to top (starting with the first five days from January, up to the last five days of December, separately for D i and Q i families).environmental conditions.We conclude with four remarks, of different nature, which could not be developed in this paper.
First, at each new calibration, a step in the series provided by the observatory will happen; but it is only of the magnitude discussed above, i.e. small.Smoothing the small steps may be considered; in the absence of extra information, no better technique than a linear correction between the two calibrations, distant by one year or so, exists.
Second, we have to stress that we only discussed the effect of the inaccuracy of the absolute measurements of B k on the values V given subsequently by the magnetometer supposed to remain identical to itself, in particular geometrically invariant.The 4 He magnetometer is built in such a way as to ensure this stability.We do not discuss either the important question of the stability of the pillar.To our knowledge, there are no available data to make any good estimation of the stability of the pillars.Therefore, it is important to build the best possible pillar and retain an adequate magnetometer.
Third, we stress again that we only made an excursion in the (calibration) error space, using the simple algorithm described in Sect. 4. A full exploration of this space would be a heavier task; in the Appendix we give a glimpse of it.
The fourth remark, which we already touched upon in Sect. 5 is, although relevant to the problem at hand, more general.Long-term stability is generally required for the study of long time scale phenomena (secular variation of the main field, solar cycle related variations, seasonal variations).For this kind of studies, what is relevant is not the absolute accuracy of one-minute values, but of some means (annual, monthly, daily, hourly); and, briefly speaking, averaging reduces the error.
A fuller understanding of the Earth's magnetic field will come from improvements in measuring it and separating its different components with a better spatial and temporal resolution.The upcoming ESA Swarm mission will provide the best-ever survey of the geomagnetic field and its temporal evolution.This constellation will benefit from a new generation of instruments, as each satellite will carry two 4 He magnetometers; these Absolute Scalar Magnetometers (ASM) are the nominal instruments for measuring the magnetic field intensity, but it is planned to operate them in vector mode, as demonstrator.

Fig. 1 .Fig. 2 .Fig. 1 .
Fig. 1.Calibration triplet B 1 , B 2 , B 3 and the geographical North-East-Down frame {u i }.The unit vectors e 1 , e 2 , e 3 (defining the physical axes) are nearly orthogonal, and each e i is close to corresponding u i .The large gray ball represents the variation of vector V (the one to be measured after calibration); small balls radii represent the measurement error εb; value δb is the upper bound for all |B i − B k |, i, k = 1, 2, 3 ; b is a typical value of the intensity of the geomagnetic field at the site and epoch of measurements.
Figure 2 presents two illustrations of the path of the vector B(X, Y, Z) during this day (see figure caption).
represents the distribution of | B| av (t) and | B| max (t) versus η.The parameter η is not discriminant enough to rank unequivocally the | B|(t) distribution; many points of the plot have the same abscissa.Nevertheless, it clearly appears that the calibration errors | B| av (t) and | B| max (t) decrease when η increases.
For each vector V (t) (1440 × 365 of them) we compute the average and maximum values of | V | over the set of (B 1 , B 2 , B 3 ).The histograms of the set of | V | av (t) and | V | max (t) are shown in Figs. 4 and 5 for εb = 0.75, 1, 2 nT.It appears that | V | av (t) (the most

Fig. 6 .
Fig.6.Sequential observations of the maximal possible calibration error for CLF data, during the year 1999; the value εb is supposed to be 1 nT.The largest values are associated with known magnetic storms of the year 1999 as those of 13 January, 18 February, or 20 October.

Fig. 10 .
Fig. 10.Time table for HER data.Same legend as for Fig. 7.

Fig. 11 .
Fig.11.Time table for BNG data, with t 0 = 6 h, α = 2 nT (left panel) and α = 4 nT (right panel) periods for the B 1 (t).Note the distinction between subsets D i and Q i : blue lines for disturbed days, red lines for quiet days.Yearly variation reads from top to down (starting with the first five days from January, up to the last five days of December, separately for D i and Q i sets).

Table 1 .
The magnetic observatories used in this study.