1 Introduction

1.1 Dynamic gravimetry

In this article, airborne and shipborne gravimetry is called “dynamic gravimetry” to distinguish these techniques from satellite and static terrestrial gravimetry. Dynamic gravimetry is widely used for cost-efficient gravity determination at medium accuracy and spatial resolution. While point-wise terrestrial gravimetry enables very precise gravity measurements at the \(\upmu \)Gal level (\(1 \, \mathrm {Gal} = 10^{-2} \, \mathrm {m/s}^2 \approx 1\,\mathrm {mg}\)), dynamic gravimetry enables fast results at the 1 mGal level (Forsberg and Olesen 2010), which is especially useful at sea, great lakes and areas that are difficult to access. In comparison with satellite and static terrestrial gravimetry, dynamic gravimetry offers an intermediate spatial resolution (from a few hundred meters in shipborne to several kilometres in airborne gravimetry).

In most dynamic gravimetry campaigns, horizontally stabilised spring-type gravimeters have been installed, e.g. Nettleton et al. (1960), Brozena et al. (1997), Forsberg and Olesen (2010), Lu et al. (2017) and Ince et al. (2020). Since the 1990s, the use of inertial measurement units (IMUs) strapped-down to the moving vehicle became more and more convenient as it was shown that the accuracy is comparable to that of spring type gravimeters (Glennie et al. 2000; Becker et al. 2015; Johann et al. 2020; Yuan et al. 2020). Advantages of the strapdown technique are the lower power, space and maintenance requirements as well as reduced costs and especially a lower sensitivity to flight manoeuvres, turbulences and rough sea conditions (Förste et al. 2020). The influence of IMU sensor drifts can be significantly reduced with an appropriate thermal calibration (Becker et al. 2015; Becker 2016) or, even better, a thermally stabilised housing (Simav et al. 2020; Johann et al. 2020).

The output quantity of dynamic gravimetry is usually the actual gravity or the gravity disturbance \(\varvec{\delta g}^n\), which is defined as the deviation of the gravity from the normal gravity \(\varvec{\gamma } ^n\). The superscript n indicates that these quantities are expressed in the navigation frame (e.g. north, east, down). Gravity is obtained as the difference of kinematic acceleration \(\ddot{\varvec{r}}^n\) and the acceleration \(\varvec{C}_b^n \varvec{f}^b\), measured by the IMU and transformed to the navigation frame. Thereby, the basic equation

$$\begin{aligned} \varvec{\delta g}^n=\ddot{\varvec{r}}^n-\varvec{C}_b^n\varvec{f}^b+\varvec{\delta g}_{\mathrm {eot}}^n-\varvec{\gamma }^n \end{aligned}$$
(1)

follows, including the Eötvös correction \(\varvec{\delta g}_{\mathrm {eot}}^n\) that accounts for Coriolis and centrifugal accelerations (Wei and Schwarz 1998).

Heading-dependent errors have been observed in several strapdown campaigns, where Q-Flex servo accelerometers of type QA-2000 manufactured by Honeywell (formerly Sundstrand) were installed in the IMUs (Becker 2016; Johann et al. 2019). The maximum of the error was observed when moving approximately north or south, with opposite sign. When moving westwards or eastwards, no significant error was found. The heading-dependent error \(\varepsilon _\psi \) was empirically described as

$$\begin{aligned} \varepsilon _\psi =-c\,{\dot{r}}_N, \; \; \delta g_{D,\mathrm {cor}}=\delta g_{{D}}-\varepsilon _\psi , \end{aligned}$$
(2)

with \({\dot{r}}_N\) being the horizontal velocity in northern direction, \(\delta g_{D},\delta g_{D,\mathrm {cor}}\) being the uncorrected/corrected gravity disturbance in the local vertical direction (Johann et al. 2019). The quantity c was estimated empirically campaign-wise since it was noticed to be approximately constant during the individual campaigns, but varied significantly between campaigns far away on the globe, especially at highly different latitudes.

The cause for this effect remained unexplained so far. Systematic instrumental errors, possibly in combination with the Eötvös correction or modelling/approximation errors, were suspected at first. In recent campaigns, the used IMU (see Sect. 2) was installed with its X axis (see Fig. 3) pointing forwards or backwards with regard to the direction of motion. This difference in the IMU mounting attitude led to an opposite sign in the estimated correction factor c, indicating that the heading angle with respect to the physical sensor frame is most relevant for the correction, not the direction of motion itself. The assumption of an erroneous Eötvös modelling can be discarded to explain the effect since the Eötvös correction

$$\begin{aligned} \varvec{\delta g}_{\mathrm {eot}}^n=(2\varvec{\varOmega }_{ie}^n+\varvec{\varOmega }_{en}^n)\cdot \varvec{{\dot{r}}}^n, \end{aligned}$$
(3)

with \(\varvec{\varOmega }_{ie}^n,\varvec{\varOmega }_{en}^n\) being the skew-symmetric matrices of the Earth rotation and the transport rate, depends on the vehicle velocity \(\varvec{{\dot{r}}}^n\) in the navigation frame, not on the IMU axis directions. Instead, the sensor readings are suspected to be impaired by a locally homogenous and approximately constant effect that is attitude-dependent, e.g. a physical field like the Earth’s magnetic field.

This paper examines the influence of a magnetic field on a strapdown IMU with Q-Flex type servo accelerometers and investigates how to correct this effect. In Sect. 1.2, basics of the Earth’s magnetic field are introduced. Section 2 briefly presents the used strapdown IMU and the processing method, and Sect. 3 presents static experiments in a Helmholtz coil. A simple calibration approach is introduced in Sect. 4, which is applied to five shipborne and airborne gravimetry campaigns (Sect. 5).

1.2 The Earth’s magnetic field

Before analysing the influence of magnetic fields on measurements in strapdown gravimetry, basics on the magnetic field of the Earth are introduced in this section. The field originates from the geodynamo mechanism in the fluid outer core of the Earth (main field, about 95% of the total field), locally magnetised rocks in the Earth’s crust (crustal field) and atmospheric plus interplanetary electric currents and magnetic fields (external field) (Lanza and Meloni 2006; Clauser 2014; Oehler et al. 2018).

The magnetic field can be expressed in terms of the magnetic field strength \(\varvec{H}\) or the magnetic flux density \(\varvec{B}\). The latter considers the magnetisation measure of a material exposed to a magnetic field, the magnetic permeability \(\mu \). Similar to gravity acceleration, which is obtained as the gradient of the gravity potential, the magnetic flux density

$$\begin{aligned} \varvec{B}=\mu \varvec{H}=-\varvec{\nabla }V \end{aligned}$$
(4)

is obtained as the gradient of the scalar geomagnetic potential V. In the following, the magnetic flux density will be called “magnetic field” for better readability. While the Earth’s gravity field is assumed approximately constant over multiple decades, especially at the accuracy level of dynamic gravimetry, the Earth’s magnetic field varies significantly over long time periods (secular variation, more than 5 to 10 years) and short time periods (Lanza and Meloni 2006).

Fig. 1
figure 1

Horizontal magnetic field intensity (\(\upmu \mathrm {T}\)) on 1 January 2020 according to IGRF-13

Full size image

The geomagnetic potential including secular variation is modelled in the International Geomagnetic Reference Field (IGRF) based on satellite and terrestrial surveys as well as observatory data. The IGRF describes the long wavelengths mathematically using spherical harmonic functions, similar to the approach commonly applied for the Earth’s gravity field (see, for example, Torge and Müller 2012). The current IGRF-13 was released in 2019 by the International Association of Geomagnetism and Aeronomy (Alken et al. 2021). With the Gauß coefficients \(g_n^m,h_n^m\) given in the IGRF-13 up to degree and order 13, the geomagnetic potential can be approximated as

$$\begin{aligned} \begin{aligned} V(r,\theta ,\lambda ,t)&=a\sum _{n=1}^{13}\sum _{m=0}^{n}\left( \frac{a}{r}\right) ^{n+1}\cdot (g_n^m(t)\cos ({m\lambda })\\&\quad +h_n^m(t)\sin ({m\lambda }))\,P_n^m(\cos ({\theta })), \end{aligned} \end{aligned}$$
(5)

with \(a=6371.2\,\mathrm {km}\), r being the radial distance from the Earth’s centre of mass to the point of observation, \(\theta ,\lambda \) being the point’s colatitude and longitude, \(P_n^m\) being the Schmidt-normed Legendre polynomials. Note that the Gauß coefficients are dependent on the observation epoch t, since the IGRF considers the secular variation of the geomagnetic field. Therefore, the coefficient sets are given in 5-year intervals from 1900 to 2020, each at January 1. Interpolation is necessary between these epochs. For a time span of 5 years after the last epoch, the predicted rate of change of the coefficients is also given up to degree and order 8. Magnetic anomaly models of higher degree than 13 are not considered in the paper at hand, since the impact of higher degree coefficients on the results of dynamic gravimetry is assumed insignificant with Gauß coefficients far below \(1\,\upmu \mathrm {T}\) (Lowes 2000).

From the magnetic potential V, the vectorial flux density \(\varvec{B}=(B_{N}\;B_{E}\;B_{D} )^T\) in the navigation frame (north, east, down) is obtained with Eq. 4. Based on the vector \(\varvec{B}\), the horizontal component \(B_{H}\) of the magnetic field and the deviation of the local magnetic north from the geographic north, the declination \(\delta \), can be calculated (Lanza and Meloni 2006) as

$$\begin{aligned} B_{H}=\sqrt{B_{N}^2+B_{E}^2},\;\; \delta =\arctan \! 2{\left( \frac{B_{E}}{B_{N}}\right) }. \end{aligned}$$
(6)

While, according to the IGRF, the total magnetic field intensity varies locally between about 20 and 70 \(\upmu \mathrm{T}\), the horizontal magnetic field \(B_{H}\) varies between 0 and about 45 \(\upmu \mathrm {T}\). Approaching the polar regions, the horizontal intensity tends to decrease since the field direction is approximately vertical at the geomagnetic poles. The global horizontal field intensity map as plotted in Fig. 1 shows also longitude-dependent variations, e.g. a low intensity in the Southern Atlantic and a high intensity in Southeast Asia. The absolute declination is less than \(30^{\circ }\) in most areas of the world except for the polar regions due to the deviation between the geographic and magnetic poles (Fig. 2).

Fig. 2
figure 2

Magnetic field declination (\({}^{\circ }\)) on 1 January 2020 according to IGRF-13

Full size image

The global magnetic field variations are not distributed homogeneously. In some regions, the magnetic field intensity rate of change in terms of secular variation amounts to several tens of nT per year, and the declination varies up to several arc minutes per year. Short-term variations of external origin effect the field additionally, e.g. diurnal effects in the order of 10...30 nT depending on latitude and solar activity; solar storms induce magnetic field variations at the Earth’s surface up to 1 \(\upmu \mathrm {T}\) (Leitgeb 1990; Lanza and Meloni 2006; Clauser 2014).

The magnetic field is affected by local influences that can be of natural origin (e.g. magnetised rocks) as well as of artificial origin (e.g. railways, mains current). The local field can be disturbed by magnetic materials close to the observing point (e.g. ferromagnetic construction material in floors, walls and ceilings in buildings) (de Vries et al. 2009) and electrical devices (Tenforde 1995; Bachmann et al. 2004). Especially in a dynamic environment, the gravimeter is surrounded by several electrical devices and cables that are located all across the vehicle body. Those lead to inhomogeneous changes of the electromagnetic field within the vehicle as the current flow is not constant during the survey. Increasing the distance to disturbing sources in a building or vehicle strongly reduces the disturbance of the magnetic field.

2 Measurement system and processing method

The base of the inertial measurement system used as a gravity meter is an iNAV-RQH-1003 by iMAR (Fig. 3) with some specific modifications, which are under manufacturer confidentiality. The main sensors in this IMU, each three Honeywell QA-2000 accelerometers (Honeywell International Inc. 2005) and Honeywell GG1320AN ring-laser gyroscopes, as well as a specific data digitisation, a temperature sensor and a clock form the so-called inertial sensor assembly (ISA). The ISA is mechanically linked to the IMU chassis via shock mounts damping vibrations and mechanical shocks which act on the IMU and damp the dithering of the gyros. Further technical details on this device, e.g. sensor specifications, can be found in iMAR Navigation (2012) and Becker (2016). Optionally, the used IMU can be enclosed in a specific housing, the iTempStab-AddOn by iMAR. Thereby, major temperature-dependent accelerometer drifts can be avoided due to a highly precise thermal stabilisation.

Fig. 3
figure 3

iMAR iNAV-RQH-1003 without iTempStab-AddOn with IMU axis directions X, Y and Z (Förste et al. (2020), adapted)

Full size image
Fig. 4
figure 4

Q-Flex accelerometer construction (reprinted by permission from Springer Nature: Lawrence (1993)) with accelerometer axis directions \({x_a,y_a,z_a}\) in green

Full size image

Since the accelerometer readings will be analysed in the following sections, it is useful to introduce the working principle of the QA-2000 accelerometer. It is a pendulous accelerometer based on the “Q-Flex” design (Fig. 4): The U-shaped pendulum and the hinge are made of a single piece of stable, insulating quartz. If the sensor is accelerated in its sensitive direction \({x_a}\), the partly metallised pendulum is pulled out of its zero position. The displacement is measured capacitively and determines the signal that is sent to a servo drive. The back swing of the pendulum is initialised using forcer coils placed around a permanent magnet. The force is proportional to the magnetic field (Lawrence 1993). Even though the Q-Flex design was initially patented in the 1970s (Jacobs 1972), Q-Flex accelerometers (like the Honeywell QA series) are still used without fundamental design changes as a reference for navigation grade accelerometers nowadays (Touboul et al. 2016).

Apparently, an exterior magnetic field might impact the accelerometer reading since it is based on a magnetic feedback system. Additionally, other parts of the IMU, e.g. the gyroscopes and the electronics, might be affected by an exterior magnetic field. Sufficient magnetic shielding is not available for the used off-the-shelf IMU and might be hard to realise if the temperature should be stabilised at the same time. In particular, for other IMUs that use magnetic sensors for attitude determination, shielding would not be an option.

For each accelerometer, a separate accelerometer coordinate system \({x_a,y_a,z_a}\) (see Fig. 4) is defined, where the axis directions depend on the specific mounting of each accelerometer inside of the IMU. The three accelerometers are mounted in the iNAV-RQH-1003 in a way that their sensitive directions \({x_a}\) are up, right, back. The IMU coordinate axes X, Y, Z illustrated in Fig. 3 are obtained by permutation and will be exclusively used in the following analysis.

Campaign data were processed applying the direct method of strapdown gravimetry where Eq. (1) is used to determine the gravity disturbance directly in the acceleration domain. That is why the direct method has also been called “accelerometry approach” (Kwon and Jekeli 2001; Ayres-Sampaio et al. 2015). In this study, the commercial software NovAtel Waypoint Inertial Explorer generated position solutions in the Precise Point Positioning (PPP) mode based on GPS and GLONASS code and phase observations plus precise satellite orbit and clock products from the Center for Orbit Determination in Europe (CODE). Inertial Explorer also performed the GNSS/IMU integration in order to calculate the sensor attitude. The kinematic acceleration \(\ddot{\varvec{r}}^n\) was obtained by a double numerical differentiation of the PPP position solution. In the shipborne campaigns, the kinematic acceleration was neglected, since the disturbing vertical movement is eliminated mostly by the more rigorous low-pass filtering (Krasnov and Sokolov 2015; Lu et al. 2019). Gravity was computed using a MATLAB software developed at the chair of Physical and Satellite Geodesy at the Technical University of Darmstadt (PSGD). The detailed processing scheme, including several corrections, is described in Johann et al. (2019).

3 Static experiments using a Helmholtz coil

In the scope of this work, the influence of the Earth’s magnetic field on the IMU introduced in Sect. 2 is examined. An artificial magnetic field with varying field intensity and field direction imitating the Earth’s field was generated using a Helmholtz coil, which consists of two narrow coils with the same direction of electric current. Inside of the coils, a homogeneous magnetic field is induced in the direction of the coil axis. The field intensity is dependent on the electric current. A 3-D magnetic field can be generated by combining a set of three perpendicular Helmholtz coils.

For the experiments, a 3-D Helmholtz coil constructed and owned by iMAR Navigation has been used (Fig. 5). The Helmholtz coil is realised using a perpendicularly arranged triple of double square coils with the square dimensions 585 mm (blue), 610 mm (red), 635 mm (yellow). Each of the square coils consists of 84 windings. This setup is capable for neutralising the Earth’s magnetic field as well as possibly relevant disturbing fields in the laboratory and generating an overlaying magnetic field exceeding \(\pm 70\;\upmu \mathrm {T}\).

Fig. 5
figure 5

iMAR iNAV-RQH-1003 encased in the iTempStab-AddOn within the 3-D Helmholtz coil at iMAR Navigation facilities with accelerometer axes X, Y, Z in attitude 1 and artificial magnetic north N, east E, down D

Full size image

3.1 Homogeneity of the generated magnetic field

The interpretation of the experiments that will be described in Sect. 3.2 is based on the assumption that the magnetic field inside of the Helmholtz coil is homogeneous. This assumption was verified by placing a solid-state 3-D magnetometer of the type iMAR iTAHS at 60 representative points in a 3-D raster within 10 cm around the central point of the coil with a point spacing of 25 mm. At all positions, 13 different magnetic fields were generated in all three IMU axis directions and the opposite directions at field intensities of 32 and 65 \(\upmu \mathrm {T}\) (\(2\cdot 6=12\) fields) plus a field with an intensity close to 0 \(\upmu \mathrm {T}\) (“zero field”) at which the exterior magnetic field is neutralised. The measurement repeatability of the iTAHS magnetometer is 10 nT according to iMAR Navigation (2016).

The actual magnetic fields were generated with deviations to the above-mentioned target fields up to 2 \(\upmu \mathrm {T}\) at the Helmholtz coil’s central point. This should not be problematic for the evaluation since the actual field intensities and directions have been used in the analysis and were approximately constant.

Figure 6 visualises the root mean square (RMS) of the magnetic field norm residuals at the individual positions. A slight increase in the RMS can be observed with increasing distance to the coil centre as well as a systematic effect in conjunction with the different horizontal layers. The overall RMS is 1.1 \(\upmu \)T.

Fig. 6
figure 6

Magnetic field RMS inside the Helmholtz coil. Each point represents a grid position and is computed with all 13 magnetic field settings. The horizontal layers indicate the vertical distance (mm) from the coil centre to the grid position

Full size image

In the static experiments, the IMU is placed in the Helmholtz coil in such a way that the centre of observation of the IMU is located in the central point of the Helmholtz coil in order to minimise errors caused by inhomogeneity of the magnetic field. Hence, all accelerometers are placed within 0.1 m around the coil centre and are expected to be accurate to about 1 \(\upmu \mathrm {T}\) assuming that the field inside the Helmholtz coil is not heavily reshaped by the IMU itself.

3.2 Methods

In order to evaluate the influence of the thermal stabilising housing, the iTempStab-AddOn (abbreviated as “iTempStab”), on the magnetic field and to verify the repeatability of the results, the static experiment was divided into two parts: First, the IMU without the iTempStab was placed in the Helmholtz coil in several attitudes in order to evaluate the response of the three accelerometers to various magnetic fields. Second, the IMU was encased in the iTempStab and was exposed to various magnetic fields, with the Z-axis pointing vertically down as shown in Fig. 5. This was also the standard orientation used for the various campaigns discussed in Sect. 5. In Sects. 3.2.1 and 3.2.2, both parts of the experiment are described in detail. The settings are summed up in Table 1. Before the start of the measurements, the IMU was turned on for several hours in order to minimise temperature-dependent accelerometer drifts.

Table 1 Overview about the static experiments
Full size table

The first and second axes of the 3-D Helmholtz coil are defining the horizontal plane, representing artificial magnetic north and east components, respectively, and the third axis is pointing down vertically. The IMU was placed inside of the Helmholtz coil in four different attitude settings as indicated in Table 2 and Fig. 7. In the standard setting, attitude 1, the IMU axes X, Y and Z are pointing towards the artificial north, east and down.

Table 2 IMU and vertical accelerometer axis directions in four attitude settings with respect to the artificial north, east and down axes of the 3-D Helmholtz coil
Full size table
Fig. 7
figure 7

Mounting and axis directions X, Y, Z of the IMU (black) and axis directions \({x_a,y_a,z_a}\) of the vertical accelerometers (green) in the four attitude settings with respect to artificial magnetic north, east, down. The circles illustrate the shape of the accelerometer housing

Full size image

In post-processing, zero field measurements were used to remove linear trends (smaller than 1 mGal/h) from the low-pass filtered accelerometer readings between these updates. All results are related to the zero field updates by setting the gravity in the corresponding epochs to zero. For every measurement under a constant magnetic field, the mean and its standard deviation are computed.

3.2.1 Part 1 without iTempStab

In the first part of the experiment, the IMU was placed in the Helmholtz coil in all four attitude settings. For all settings, magnetic field directions were generated in all three IMU axis directions and in the opposite directions, just like in the homogeneity verification. The field intensity was 65 \(\upmu \mathrm {T}\), which is approximately the maximum total intensity of the Earth’s magnetic field. At least before and after each attitude change, zero fields were generated. Hence, the IMU was exposed to six nonzero magnetic field directions per attitude plus in total eight zero fields. Each field was set for three minutes allowing averaged accelerometer readings at an accuracy level of some tenths of mGal.

3.2.2 Part 2 with iTempStab

In the second part of the experiments, which was conducted about three months later, the IMU encased in the iTempStab was tested in the standard attitude setting solely. In order to obtain more detailed and more accurate results in this setting, more fields have been generated and each magnetic field was held for five minutes allowing more precise averaged accelerometer readings. The same vertical fields like in the first part of the experiments were generated. Horizontal fields were generated in \(45^{\circ }\) intervals at field intensities of 10, 20 and 32 \(\upmu \mathrm {T}\) for all directions plus 65 \(\upmu \mathrm {T}\) for four of the eight directions. In addition, the zero field has been measured seven times leading to a total of 39 field measurements. The horizontal and vertical fields will be analysed independently.

3.3 Results

This paper focuses on the reaction of the vertically oriented accelerometer (downwards positive) to the magnetic field since the vertical accelerometer is assumed to be most relevant for dynamic gravimetry. Furthermore, already small instabilities in the attitude angles at the order of few arc-seconds lead to errors in the observed horizontal acceleration of several mGal (e.g. \(1''\rightarrow 4.8\,\mathrm {mGal}\)), which impedes the evaluation of the horizontal accelerometer readings. Hence, in accordance with Table 2, for the four attitude settings, the analysis in this section is related to the IMU axes \(Z, -\,Z, Y\) and X, respectively.

First of all, it should be noted that the magnetic field significantly affected the accelerometer readings in all attitude settings: The maximal deviations compared to the zero field were between 3 and 6 mGal.

3.3.1 Part 1 without iTempStab

In the first part of the experiment, the IMU was placed in the Helmholtz coil in the four attitude settings shown in Table 2 without the iTempStab. At first, the effects of the horizontal magnetic field component are analysed. Figure 8 visualises the mean deviation from the zero field readings in the four declination angles of the horizontal magnetic field. The vertical field was approximately zero. The 65 \(\upmu \mathrm {T}\) field strength caused errors in a sinusoidal pattern in all attitude settings with varying amplitudes and phase shifts. Note that the readings of attitude 4 were impaired by a too short warm-up period before the start of the measurements.

Fig. 8
figure 8

Means and standard deviations of the down component as a function of the horizontal magnetic field direction (without iTempStab; blue: zero field; green: field intensity of 65 \(\upmu \mathrm {T}\))

Full size image
Fig. 9
figure 9

Most relevant horizontal magnetic field directions in the four attitude settings. The red (blue) arrows indicate the horizontal field direction of greatest positive (negative) influence on the down sensor reading

Full size image

In the attitudes 1 and 2, the same accelerometer (Z) has been analysed after a rotation around the X-axis. The turnaround of the sensitive axis \({x_a}\) of the accelerometer was taken into account by changing the sign of its readings in Fig. 8b. Both plots can be approximated by a cosine function with an amplitude of about 5 mGal, but there might be a small phase shift of a few degrees. The error of the down component induced by a horizontal magnetic field appears to be determined by its direction relative to the accelerometer axis directions \({y_a,z_a}\). Depending on this input direction, the measured absolute gravity reading always increases or decreases by a specific amount, no matter if the sensitive axis of the sensor is pointing downwards or upwards (Fig. 9). This hypothesis is supported by the attitudes 3 and 4. In attitude 3, the Y accelerometer inside of the IMU is oriented exactly like the \(-\,Z\) accelerometer in attitude 2, showing a similar behaviour with differences regarding the sign of the phase shift and the amplitude (about 4.5 mGal in attitude 2 and about 3 mGal in attitude 3). The X accelerometer in attitude 4 is rotated around its sensitive axis by \(90^{\circ }\) in comparison with accelerometer Z in attitude 1, which might explain the phase shift of about \(90^{\circ }\).

Fig. 10
figure 10

Means and standard deviations of the down component as a function of the vertical magnetic field direction in the four attitude settings (without iTempStab)

Full size image
Fig. 11
figure 11

Most relevant vertical magnetic field directions in the four attitude settings. The red (blue) arrows indicate the horizontal field direction of greatest positive (negative) influence on the down sensor reading

Full size image

Second, the effects of the vertical magnetic field component are analysed. While the absolute values of the mean gravity values are different for the different attitude settings, they are approximately equal for the upwards and downwards directions regarding the individual attitude settings. The sign is dependent on the vertical field direction (Fig. 10). Like for the horizontal fields, the reaction of the accelerometer to a vertical field seems to depend on the direction of the magnetic field with respect to the accelerometer axis directions. Since the field input direction is parallel to the accelerometer’s sensitive axis \({x_a}\), the reaction of the accelerometer to a vertical magnetic field depends on the mounting direction of the sensitive axis (Fig. 11). The amplitude of the error is different for both attitudes (1 and 2) of accelerometer Z.

To conclude, the above-mentioned results suggest that for the tested IMU, the direction and magnitude of the horizontal and the vertical magnetic field relative to the axis directions of the vertical accelerometer cause a reading error of the type of a scale factor. Hence, the orientation of the sensor-triad inside the IMU determines the error characteristics. All three accelerometers show similar reactions to changes in the magnetic field. However, the amplitudes differ for the individual sensors. For vertical fields, in addition, the amplitude of the response of each sensor is different depending on the direction of the vertical magnetic field relative to the sensitive axis of the accelerometer.

3.3.2 Part 2 with iTempStab

In the second part of the experiment, the IMU was encased in the iTempStab and was exposed to multiple magnetic fields of various intensities and directions in attitude 1. Under thermal stabilisation, the effect of temperature changes on the accelerometer readings is expected to be minimised. Figure 12 shows the mean gravity values for the horizontal and vertical magnetic fields. The calibration function that is also included in the figure will be introduced in Sect. 4.

Fig. 12
figure 12

Means and standard deviations of the down component as a function of the horizontal or vertical magnetic field direction (attitude 1, accelerometer Z, with iTempStab). The solid lines illustrate the calibration function at the tested magnetic field intensities

Full size image

The measurements with a field intensity of 65 \(\upmu \mathrm {T}\) can be compared to part 1 of the experiment, attitude 1. All differences at the corresponding magnetic field angles are smaller than 0.4 mGal, with a mean of 0.01 mGal and an RMS of 0.19 mGal. This indicates, firstly, that the results are repeatable and, secondly, that the iTempStab housing does not reshape the local magnetic field in a way that has a significant influence on the vertical accelerometer readings. The absolute mean gravity increases with increasing magnetic field intensity.

4 Calibration approach

This section proposes a basic magnetic field calibration approach with the purpose of correcting the vertical accelerometer readings of the IMU for the Earth’s magnetic field during a dynamic gravimetry campaign, using the results of the static experiments described in the previous section. Since the differences between the measurements with and without the iTempStab were small, it is assumed that the calibration will work with sufficient accuracy in both cases.

While the heading variation of the vehicle regularly covers the full circle of \(360^{\circ }\) during dynamic gravimetry surveys, the vehicle roll and pitch angles are usually small except during aircraft manoeuvres. Hence, the IMU attitude typically is similar to the attitude 1 of the static experiments. Furthermore, the vertical accelerometer (IMU Z axis) is most important for the determination of the vertical gravity. The results of the second part of the static experiments were obtained in attitude 1 under various magnetic fields (Fig. 12a) and were used as input data for the calibration. The influence of vertical magnetic fields is neglected in the calibration, since it is assumed approximately constant during a typical gravimetry campaign; a constant bias is irrelevant for relative gravimetry. The function

$$\begin{aligned} \varDelta g_{\mathrm {mag},D}=(c_1 B_{H}+c_2 B_{H}^2+\cdots +c_k B_{H}^k)\cdot \cos {(\alpha )} \end{aligned}$$
(7)

is used as a basis for the calibration function, with the corrected vertical gravity disturbance

$$\begin{aligned} \delta g_{\mathrm {cor},{D}}=\delta g_{\mathrm {unc},{D}}+\varDelta g_{\mathrm {mag},{D}} \end{aligned}$$
(8)

being the sum of the uncorrected gravity disturbance \(\delta g_{\mathrm {unc},{D}}\) as computed with Eq. (1) and the correction \(\varDelta g_{\mathrm {mag},{D}}\). The correction is approximated as the cosine of an angle \(\alpha \) with an amplitude given by a polynomial of the order k depending on the horizontal magnetic field intensity \(B_{H}\). The polynomial coefficients \(c_i\) do not contain a constant \(c_0\) since the correction is defined to be zero at a zero field. The horizontal angle \(\alpha \) is the angle between local magnetic north and the horizontal direction at the IMU where a horizontal magnetic field causes the maximal shift of the vertical accelerometer readings and depends on the vehicle heading \(\psi \), the local magnetic field declination \(\delta \), the angle \(\beta \) between the vehicle front axis (F) and the IMU front axis (X) and the angle \(\kappa \) between the IMU front axis (X) and the direction of the maximal susceptibility to a horizontal magnetic field of the vertical accelerometer. The angles required to determine

$$\begin{aligned} \alpha =\psi -\delta +\beta +\kappa \end{aligned}$$
(9)

are illustrated in Fig. 13. The polynomial coefficients \(c_i\) and the angle \(\kappa \) are estimated in a least squares adjustment (Gauß–Markov model).

Fig. 13
figure 13

Horizontal angles relevant for the calibration approach with the horizontal axes pointing towards geographic north \(N_\mathrm {geo}\), magnetic north \(N_\mathrm {mag}\), the vehicle front axis F, the IMU front axis X and the IMU direction \(\varDelta g_\mathrm {max}\) of maximal susceptibility to a horizontal magnetic field

Full size image

Higher-order coefficients of the calibration functions turned out not to be significant. Hence, the error of the respective vertical accelerometer due to a horizontal magnetic field can be described as being proportional to the horizontal field intensity. Accordingly, the calibration function

$$\begin{aligned} \varDelta g_{\mathrm {mag},D}=c_1 B_{H} \cos {(\psi -\delta +\beta +\kappa )} \end{aligned}$$
(10)

was applied in this study. For the examined IMU, the dependency of the horizontal field intensity was estimated as \(c_1\approx (85.0\pm 1.8)\mathrm {\frac{\upmu Gal}{\upmu T}}\), and the direction of the maximal susceptibility of the sensor as \(\kappa \approx (7.0^{\circ } \pm 1.1^{\circ })\). The calibrated function is illustrated in Fig. 12a for several horizontal field intensities.

5 Verification in shipborne and airborne campaigns

The presented simple magnetic field calibration function was tested in several dynamic gravimetry campaigns. First, a shipborne campaign at Lake Müritz in 2019 is presented. Then, an overview on statistics of four additional campaigns is given referring to recent publications that introduced these campaigns.

5.1 Example case: Müritz 2019

In November 2019, shipborne gravimetric test cruises have been conducted by the GeoForschungsZentrum Potsdam (GFZ) in cooperation with PSGD at Lake Müritz, the second largest lake in Germany. The iMAR iNAV-RQH-1003, encased in the iTempStab, was installed at the indoor passenger area of the excursion boat “Klink” (Fig. 14). The NovAtel GNSS receiver integrated in the IMU recorded multi-GNSS and multi-frequency signals of a GNSS antenna installed on top of the roof at the bridge. Two measurement cruises were conducted with durations of about 7 and 9 h. The trajectories were planned with the focus on three approximately straight lines that have been repeated multiple times. The GFZ performed terrestrial gravity measurements in order to establish the absolute gravity reference. Further statistics on the campaign can be found in the last column of Table 3.

Fig. 14
figure 14

Excursion boat “Klink”, owned by “Weiße Flotte Müritz”, at Waren/Müritz

Full size image

The quality of the results is assessed using the gravity disturbance residuals at line crossings (cross-over points) and at repeated lines with a maximal horizontal line separation of 200 m. If a line has been repeated multiple times, each pair is regarded separately in order to obtain comparable precision values for repeated lines and cross-over points in this paper. Assuming equal accuracy over all trajectories, the precision of the gravity disturbance can be given in terms of the RMS error (RMSE), which is obtained as the RMS divided by \(\sqrt{2}\). The RMS is computed using all cross-over residuals or line pair differences at once. The influence of sensor drifts can be reduced performing a cross-over adjustment where a bias is removed from each cruise/flight (stint-wise) or individually from each (straight) trajectory line (line-wise). Similarly, in order to focus on short wavelength residuals in the repeated lines analysis, “zero-mean” data can be regarded where the mean is removed from each line.

Fig. 15
figure 15

Vertical gravity disturbance (mGal) Lake Müritz 2019 (map data: Google, Europa Technologies, GeoBasis-DE/BKG)

Full size image

The obtained down component of the gravity disturbance along the trajectory is illustrated in Fig. 15. Without adjustment and heading-dependent correction, a cross-over RMSE of 1.46 mGal was obtained for the 67 cross-over points. Applying the heading-dependent correction based on Eq. (2), where the campaign constant is determined empirically, the RMSE is strongly reduced to 0.27 mGal. The same RMSE is obtained using the calibration function (Eq. (10)) based on the experiments with the Helmholtz coil. After a line-wise cross-over adjustment with 48 remaining cross-over points, the RMSE was further reduced to 0.10 mGal (0.11 mGal) with the empirical correction (calibration), the best values that have been obtained in PSGD’s campaigns so far. Since the heading is approximately constant at the lines, there are no significant differences between the line-wise adjusted RMSE when a heading-dependent correction is applied. The repeated lines analysis confirms the cross-over precision with RMSE of 1.19 mGal before and 0.26 mGal (0.28 mGal) after having applied the heading-dependent corrections, with the calibration correction being marginally worse than the empirical correction. After a line-wise bias removal, the RMSE is lowered to 0.13 mGal.

Table 3 Strapdown gravimetry campaign overview (none: no heading-dependent calibration; empirical: empirical heading-dependent correction based on Eq. (2); calibration: horizontal magnetic field calibration according to Eq. (10); MY results based on Johann et al. (2019), ODW and BTS results based on Johann et al. (2020))
Full size table

5.2 Results of other campaigns

In order to verify the developed magnetic field calibration approach, it has been applied to various further airborne and shipborne strapdown gravimetry campaigns with the IMU that was also used in the static experiments and at Lake Müritz (“MRZ2019”). The campaigns are briefly introduced below; for detailed campaign descriptions and evaluations, the reader is referred to Johann et al. (2019) for the Malaysia 2014 campaign and Johann et al. (2020) for the other campaigns. The campaign statistics are shown in Table 3.

  • “MY2014”: Using a medium-size aircraft, a Beechcraft King Air 350, an airborne gravimetry campaign was undertaken above the South China Sea northwest of Borneo in 2014.

  • “ODW2017”: In 2017, several test flights were conducted at a German region with high gravity variations, the Upper Rhine Graben and the low mountain range Odenwald near Darmstadt. The quality of the results was limited by vibrations of the GNSS antenna and the instable flight with the small motor glider of type Grob G 109B.

  • “BTS2017” and “BTS2018”: In the scope of the FAMOS project (Finalising Surveys for the Baltic Motorways of the Sea), shipborne gravimetry campaigns were undertaken with the survey, wreck-search and research vessel DENEB in German and adjacent Danish, Swedish and Polish territories in the Baltic Sea in 2017 and 2018. In 2018, PSGD’s thermally stabilised housing (iTempStab) was used for the first time. Some cruises from port to port were exceptionally long taking up to 60 h. Hence, the influence of long-term sensor drifts was increased compared to the other campaigns with stints of less than 8 h.

Table 3 includes the precision indicators of all five campaigns: the cross-over RMSE and the repeated lines RMSE if available. Line-wise adjustments were computed if at least two cross-over points were available at the individual flight lines, stint-wise adjustments were computed if at least three stints/flights were conducted. Line data were used for the evaluation; sections of minor accuracy or harsh sea conditions have not been removed from the results. For this reason and also due to minor changes in the algorithm, the precision indicators slightly differ from previous publications.

In all campaigns, the precision of the gravity data was significantly improved applying a heading-dependent correction. Without adjustment, the empirical heading-dependent correction performs on a par or slightly better than the heading-dependent magnetic field calibration. The biggest improvements using the magnetic field calibration were observed at the MY2014 and MRZ2019 campaigns, where the non-adjusted cross-over RMSE was lowered by 54 and 82%, respectively. The other campaigns showed improvements between 7 and 18%. Applying the calibration, the cross-over RMSE was between 0.27 and 2.94 mGal before adjustment and between 0.11 and 1.80 mGal after adjustment. The best precision was reached at the shipborne campaigns using the iTempStab. In addition to MRZ 2019, repeated line data were collected at the MY2014 campaign, with an RMSE of 1.13 mGal (0.71 mGal for zero-mean data).

6 Discussion and outlook

Both heading-dependent corrections significantly improved the precision. However, in most campaigns, the empirical heading-dependent correction performed slightly better than the presented magnetic field based heading-dependent calibration approach. Possible reasons are listed below:

  • Since the constant in the empirical correction is user-defined, non-heading dependent systematic or arbitrary errors might be removed by this correction, leading to a too optimistic precision estimation.

  • Deviations of the Earth’s true magnetic field from the IGRF model and disturbing fields in the vehicles might impair the results. To overcome this, readings of a magnetometer that might be installed close to the ISA inside or outside of the IMU might be used instead of the IGRF data as input for the calibration function. However, in the analysed campaigns, the IGRF model seemed to be sufficient to eliminate at least the bulk of the magnetic effect.

  • The accuracy of the Helmholtz coil, the homogeneity of its field and, thus, the accuracy of the performed calibration are limited. Calibration performance might be optimised using a larger Helmholtz coil giving an increased volume of homogeneous magnetic field.

  • The basic calibration model might be further improved. By extending the model with the effect of a vertical magnetic field and by analysing the X and Y accelerometers in detail, slight improvements of the results are conceivable for trajectories where roll and pitch angles deviate from zero significantly, especially during airborne manoeuvres.

The greatest improvement applying heading-dependent corrections was observed at the MRZ2019 campaign, where the heading of most lines was close to the north–south or opposite direction (mean deviation only about \(25^{\circ }\), see Table 3) and where a lot of cross-over points were obtained on almost opposite tracks, which maximises the heading-dependent effect. The second largest improvement at the MY2014 campaign is suspected to be due to the high mean horizontal magnetic field intensity near Borneo of 40.6 \(\upmu \mathrm {T}\), which is close to the global maximum according to the IGRF (see Fig. 1) and approximately twice the value of the other campaigns with horizontal intensities between 17 and 21 \(\upmu \mathrm {T}\).

The presented findings point out that the sensors of further strapdown gravimeters might be susceptible to environmental influences like magnetic fields. If other gravimeter types, especially but not exclusively gravimeters with Honeywell Q-Flex accelerometers, are also affected by a magnetic field, basic calibration approaches like the one presented in this paper might yield further accuracy gains in strapdown gravimetry.

In simplified versions of the static experiment in the 3-D Helmholtz coil introduced in Sect. 3.2.1, two other IMUs of the same series showed dependencies on the magnetic field in the same order of magnitude, but a device-specific calibration seems to be advisable.

Regarding other strapdown IMUs, it is recommended to reappraise their results by analysing if a systematic heading-dependent behaviour can be observed, e.g. by visualising cross-over point gravity residuals as a function of the corresponding line heading angles. If a systematic behaviour is detected, a calibration in a 3-D Helmholtz coil would be desirable, but an empirical heading-dependent correction based on Eq. 10 might be sufficient to remove the bulk of the erroneous effect by estimating the amplitude \(c_1\) and, possibly, the angle \(\kappa \). The remaining parameters in Eq. (10) can be measured, calculated or taken from the IGRF.

It is expected that the presented compensation will be considered soon in state-of-the-art strapdown gravimeters like the iMAR iCORUS (iMAR Navigation 2019). The discussed results might also have a relevant impact on certain inertial navigation systems (INS) used for navigation, especially during dead reckoning or periods of free inertial navigation.

7 Conclusions

So far, in the literature, little attention has been paid to the influence of the Earth’s magnetic field on strapdown gravimetry. In the scope of the paper at hand, it was shown in static experiments that the accelerometer readings of the evaluated IMU were influenced by a magnetic field even only within an intensity at the order of the Earth’s magnetic field. A basic calibration approach for the compensation of the magnetic field impact on the vertical accelerometer was introduced, based on the horizontal magnetic field direction and intensity. The approach was verified in diverse airborne and shipborne strapdown gravimetry campaigns with completely different vehicles, line velocities, measurement regions and horizontal mounting directions of the evaluated IMU. It was demonstrated that the calibration improved the precision of diverse airborne and shipborne strapdown gravimetry campaigns between 7 and 82%.