Solving the orientation problem for an automatic magnetic observatory

Solving the orientation problem for an automatic magnetic observatory A. Khokhlov, J. L. Le Mouël, and M. Mandea International Institute of Earthquake Prediction Theory and Mathematical Geophysics 79, b2, Warshavskoe shosse, 113556 Moscow, Russia Institut de Physique du Globe de Paris – UMR7154, 1, Rue Jussieu, 75005 Paris, France Centre National d’Etudes Spatiales, 2, Place Maurice Quentin, 75001 Paris, France


Introduction
The geomagnetic field is continously 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, Alldreage planned an automatic standard magnetic observatory (ASMO) - (Alldregde, 1960;Alldregde and Saldukas, 1964) Eastern component and vertical component of the geomagnetic field (without extra independent absolute measurements, at least for a long enough time span).This idea has remained in the geomagnetism community and, over the last years, attempts have been made to automate a DI-theodolite (van Loo and Rasson, 2006), a proton vector magnetometer (Auster et al., 2006) or to build a device which can perform discrete absolute measurements automatically (Auster et al., 2007).
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' électronique et de technique de l'information (LETI) of the French Commissariat à l' énergie atomique (CEA), or the Oersted space mission fluxgate magnetometer, built by the Institute for automatisation 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 beenextended 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 Introduction

Conclusions References
Tables Figures

Back Close
Full 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 Alldredge's statement, to obtain one minute absolute values of the field components X , Y , Z in, say, the North-East-Verical Down frame, fitting Intermagnet standard (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 magnetometervariometer, 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 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, ... 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 .. K of the components of the magnetic vector V (t k ) Introduction

Conclusions References
Tables Figures

Back Close
Full 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 magnetometervariometer 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 again that we note B(t 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 independant 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 magnetometer-variometer.We have:

Conclusions References
Tables Figures

Back Close
Full 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 obtained straightforwardly in function of the u i through the B k : where F and D are the matrices of coefficients (components) f j k and d j k .
The magnetometer-variometer provides the values and (V 1 (t), V 2 (t), V 3 (t)) of the physical componenets of V : The problem is solved in the error-free case; we have obtained the geographical componenets (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:

Conclusions References
Tables Figures

Back Close
Full The values of the components along the physical axes {e i }, f j k are not either errorfree.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 (Gravrand et al., 2001;L éger et al., 2002) 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 (there is no difficulty nor methodological change when adding errors on f .It now immediately comes: where ω 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

Conclusions References
Tables Figures

Back Close
Full 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 time span of say a week (see Sect. 4 for numerical values).
The matrix F whose lines are close to one another is a priori poorly conditionned; its inverse F −1 may have large eigen-values, and a strong amplification of error εb might affect the directions of e i (i.e. an error on e i might be amplified 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 erroneus 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) 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 dirrectly 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. (9) 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 (Eq.3).The balls are drawn in dark grey 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 variometermagnetometer, in different locations at the Earth's surface.

Observatory data
We use the hourly means of the three components of the field recorded during the year 1999, as available on the INTERMAGNET 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 farther apart the B k are, the better the configuration is for calibration.Then, at a given observatory, the larger the Introduction

Conclusions References
Tables Figures

Back Close
Full magnetic activity, the larger the possibility 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 indexes (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, specificially the day 6 September 1999, a quiet day belonging to Q CLF .Figure 2 presents two illustrations of the walk of the vector B(X , Y , Z) during this day (see figure caption).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. ( 9) -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 grey ball of Fig. 1); note that the set of vectors V is partly simulated.We compute, for a given calibration triplet B k , the mean value of |∆V | = |V − V | (Eq.9) and its maximum value taken on the set of (B 1 , B 2 , B 3 ) and the set of V .Those quantities |∆V | av and |∆V | max are in this way computed for each of the B k triplets B 1 = B(t), B 2 = B(t + t 0 ), B 3 = B(t + 2 t 0 ), for the (1441 − 2 t 0 ) values of t, as defined above.We are left with the populations |∆B| av (t) and |∆B| max (t).We order them using a parameter η which grossly characteristizes the quality of the configuration B k , i.e. the aperture of this triplet.We choose:

GID Introduction Conclusions References
Tables Figures

Back Close
Full < > is for the mixt product.Figure 3 represents the distribution of |∆B| av (t) and |∆B| max (t) versus η.The parameter η is not discriminant enough to rank univoqually the |∆B|(t) distribution; many points of the plot have the same abscissa.Nevertheless, it appears clearly that the calibration errors |∆B| av (t) and |∆B| max (t) decrease when η increases.

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, specificially 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 ).(the most 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 presentation of the calibration error, more practical, which gives the time-spans during which this error is Introduction

Conclusions References
Tables Figures

Back Close
Full smaller than a given threshold of αnT.We choose again values of year 1999, consider the triplets B 1 = B(t), B 2 = B(t + t 0 ), B 3 = B(t + 2 t 0 ), and compute the corresponding calibration errors |∆B| av (t) as explained in the last but one paragraph.Figures 7 to 11 present the results for the four observatories; values of the parameters are given in the captions.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 quiestet days of December at its top.
In Figs. 7 to 10 the considered value 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 upper panels, t 0 = 6.0 h for the lower 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 and δb = 50 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; otherwhise, a colored tiret, red or blue, is drawn.Continous red or blue time intervals are such that, for any first measurement with t in this interval, leads to a calibration error smalles than 2 nT on V (t).
Looking at graphs of Figs. 7 to 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  measurements.The values provided by the magnetometer-variometer should be exact in absolute value within a few nT (it is said 5 nT in 99 % of cases).From the results of Sect.4, it appears that, after the calibration performed as in Sects. 2 and 3, the magnetometer-variometer -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 checkings.We conclude with three remarks, of different nature, which could not be developed in this paper.First, 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.Second, 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 third remark, that we already touched upon in Sect. 5 is, although relevant to he problem at hand, more general.Long term stability generally is 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 full understanding of the Earth's magnetic field will come from improvements in measuring it and separating its different componenets with a better spatial and temporal resolution.The upcoming ESA Swarm mission will provide the best-ever survey of the Figures , i.e. a device providing at each time the absolute values of -say -the Northern geographical component, Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Here ε 1 and r j k are O(1).In matrix form: D = D + εb R. 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.
Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | For each vector V (t) (1441 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) Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | is an exemple of the effect of changing t 0 .6 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 (generaly weekly) man-made absolute Discussion Paper | Discussion Paper | Discussion Paper |