New Results from Strapdown Airborne Gravimetry Using Temperature Stabilisation

Abstract

In recent years, the use of a strapdown Inertial Measurement Unit (IMU) for airborne gravimetry has proven itself to be an accurate and resilient measurement system, improving the operational flexibility. The main concern is erroneous long-wavelength information in the resulting estimates, which is suspected to originate from uncompensated long-term drift of the accelerometers, probably originating from temperature variation. For this reason, iMAR navigation has designed a temperature stabilisation box, which allows for temperature stabilisation of their IMU systems. On a regional airborne gravity survey over the Kattegat Sea between Denmark and Sweden, such a temperature stabilised strapdown IMU was operated alongside a traditional spring-type platform-stabilised gravity system from ZLS. An analysis of the difference in gravity estimates at cross-over locations yielded a mean value of −0.3 mGal for the iMAR system with an indicated accuracy of 1.0 mGal. The temperature stabilisation unit therefore effectively limits the accelerometer drift and improves the long-wavelength information. However, a straightforward merging approach, adjusting the line-based mean values of the iMAR estimates to match that of the ZLS estimates, improved the accuracy to 0.8 mGal. This indicates that the long-wavelength information of the stabilised-platform system is still superior to that of the strapdown system.

1. Introduction

Measuring gravity from an airborne platform, commonly known as airborne gravimetry, has two fundamental challenges [1]: (1) Separating out the kinematic accelerations induced by aircraft movement and (2) determining the sensor orientation during aircraft dynamics. The advent of the Global Positioning System (GPS) in the 1990s provided a reliable and accurate solution to the first challenge, while the second challenge has traditionally been approached by mounting the gravimeter on a stabilised platform. Initially, marine gravimeters such as the LaCoste&Romberg (L&R) spring-type gravimeter [2], was modified for airborne purposes [3,4]. Later on, gravimeters have been designed for airborne applications [5,6,7], bringing along improvements in both resolution and accuracy of the gravity estimates. In recent years, quantum gravimeters have also been used for dynamic gravity observations [8,9]. However, since the price tag reflects both the market and production efforts, airborne-dedicated systems are usually more expensive.

Parallel to the platform systems, an alternative approach based on inertial technology has been evolving since early tests in the 1980s [10,11,12]. In this approach, the mechanical platform is omitted and the measurement system is mounted in a strapdown configuration. With no platform system, the orientation is solved for numerically using a triad of gyroscopes that measure aircraft rotation. Furthermore, in this configuration, the sensitive axis of the accelerometer will generally not be aligned with the direction of gravity, meaning that a triad of accelerometers are required to measure the full magnitude of the gravity vector. Although such a strapdown configuration poses more stringent requirements for the sensor performance, e.g., larger dynamic range, higher resolution and better scale factor stability, it does have a number of significant practical advantages over the platform system:

• An off-the-shelf Inertial Measurement Unit (IMU), designed for navigation purposes, can be used:
  • The instrument simultaneously provides a high-resolution navigation solution
  • As the market for navigation systems is much larger, the price is lower
• The mechanical platform is omitted:
  • As processing methods evolve, a better orientation can be obtained
  • No operation is required during flight
  • Smaller size
  • Less power consumption
  • Lower fail rate
  • Increased operational flexibility

The usefulness of the strapdown system has been repeatedly demonstrated [12,13,14,15] and a comparison test has shown that the accuracy is comparable to that of the platform systems [16]. In collaboration with the Technical University of Darmstadt, the Danish National Space Institute (DTU Space) has been operating such a strapdown IMU system along-side a platform-stabilised L&R gravimeter on a number of airborne campaigns since 2013. Our experience from these campaigns is that the strapdown system has a much better dynamic range, which is important in terms of operational effectiveness, robustness and flexibility. In 2016 DTU Space purchased the iNAT-RQH-4001 navigation-grade IMU from iMAR navigation to be used as a strapdown gravimeter. This gravimeter has confirmed our initial experiences on a number of campaigns [17,18], including the ability of draped flying, i.e., altitude following the terrain.

Despite these results, the platform system has remained the preferred choice for geodetic purposes. This is because these previous studies have also demonstrated the presence of erroneous long-wavelength information in the resulting gravity estimates. This erroneous information has been proposed to originate from uncompensated long-term drift in the accelerometers, leaking into the gravity estimates. Considering that the accelerometers are sensitive to temperature variations, it has also been demonstrated that the gravity estimates can be improved by applying either precise temperature control [19,20] or off-line calibration methods [21]. In continuation of these efforts, iMAR navigation has designed an additional temperature stabilisation box, the iTempStab-AddOn, which allows temperature stabilisation of the entire iNAT system to a precision better than 0.1 $°$C at stable ambient room temperature. This paper will present results from an airborne gravity survey, using a prototype of this temperature stabilised strapdown gravity system.

Although temperature stabilisation potentially offers a convenient and effective way of improving the long-wavelength information in the gravity estimates, it also challenges some of the advantages of the IMU that have been outlined. First of all, we are no longer dealing with an off-the-shelf IMU, but rather a system modified for gravimetry. The additional temperature box brings along additional space, weight and power requirements as listed in

Table 1. Technical specifications of the iNAT-RQH-4001 IMU and the additional temperature stabilisation box (iTempStab-AddOn). Notice that the IMU fits inside the temperature box, which adds extra weight and power consumption

	iNAT-RQH-4001	iTempStab-AddOn
dimensions	19 × 13 × 30 cm	25 × 22 × 38 cm
weight	8 kg	+10 kg
power	10–36 V, <25 W	19–31 V, max. 175 W

. Furthermore, since it takes up to 24 h for the temperature to fully stabilise, the system must be constantly powered during the entire airborne campaign, limiting the operational flexibility. However, since heat is generated internally by the IMU, the non-temperature-stabilised system will also benefit from an initial 24 h saturation period and being constantly powered. It should also be noted that the L&R system has the same requirements in terms of temperature stabilisation.

2. Gravity Sensor Analysis

In order to investigate the performance of the iTempStab-AddOn before taking it into the field, we conducted a number of static recordings on the concrete floor of the basement at our institute. Although all six inertial sensors of the IMU are important in the processing, the vertical accelerometer will measure the majority of the gravity variation. The following analysis will focus on this sensor.

The raw observations are dominated by random noise, which was attenuated by applying a two-pass Butterworth filter with a time constant of 300 s. As the internal sensors generate heat, a temperature saturation period of approximately 24 h was observed, after which the dominant trend appears linear. Figure 1 shows the time series of an approximately four day static recording at a sampling rate of 2 Hz. The temperature variation for this recording is within 0.04 $°$C. The accelerometer observations have been low-pass filtered and corrected for a linear trend of −869 $\mathsf{\mu }$Gal/day. Also shown is the solid Earth tide, computed at the measurement location using the ETGTAB software [22]. The correlation coefficient between the two time series is 0.5.

An additional static recording of approximately 29 h at the full sampling rate of 300 Hz was collected for a sensor analysis in terms of the Allan variance method [23]. The Allan variance, or equivalently the Allan deviation, is a measure of the root mean square (RMS) random drift error as a function of averaging time. Since the dominant type of random process is dependent on the averaging time, i.e., the amount of filtering, the Allan deviation curve can be used to identify and quantify various noise terms existing in the data [24]. Specifically, on a log-log plot of Allan deviation versus averaging time, various noise terms are related to the slope of the curve as shown in Figure 2 and

Table 2. Error coefficients identified from the Allan deviation curve in Figure 2. Units are in terms of the standard gravity value $\mathrm{g}=9.80655$$\mathrm{m}/{\mathrm{s}}^2$.

Error Term	Slope	Read-Off	Noise	Standard Dev.
quantisation noise	−1	$\sqrt{3}$	1325.91	53.21	$\mathsf{\mu }\mathrm{g}/\mathrm{Hz}$
velocity random walk	−1/2	1	145.11	14.58	$\mathsf{\mu }\mathrm{g}/\sqrt{\mathrm{Hz}}$
bias instability	0	minimum	3.36	0.58	$\mathsf{\mu }\mathrm{g}$
acceleration random walk	1/2	3	5.32	1.42	$\mathsf{\mu }\mathrm{g}/\sqrt{\mathrm{h}}$
drift rate ramp	1	$\sqrt{2}$	2.60	1.62	$\mathsf{\mu }\mathrm{g}/\mathrm{h}$

.

From Figure 2 it is evident that quantisation noise, represented by the part of the curve with slope −1, is the dominant noise source at short averaging times. From this long part of the curve, the slope transitions through −1/2 (velocity random walk), to 0 (bias instability), to +1/2 (acceleration random walk) and ends with +1, representing the long term drift of the error. The magnitude of each of these components can be read-off the corresponding straight lines at certain averaging times as listed in Table 2. For a more thorough description, the reader is referred to [24,25].

3. Survey Overview, Instrumentation and Operation

The Kattegat airborne gravity survey was carried out in October 2018, covering the Kattegat ocean area between Denmark and Sweden, including the surrounding coastal areas and islands. The purpose of the survey was to provide data in support of improved geoid models, augmenting the existing gravity data in the area, represented by marine observations from the 1960s.

The 22 flight lines are shown in Figure 3 and amounts to 5177 line kilometres with 7.5 km spacing in between the lines. In order to validate the results, four crossing lines were flown, resulting in 63 cross-over points. The average ground speed was 83 m/s and the average flight altitude was 930 m with respect to the WGS84 reference ellipsoid. Three of the flight lines, more specifically line Q and the crossing lines between X5/X6 and X7/X8, were flown at a higher altitude of 1240 m due to weather conditions.

The survey was carried out using a Beechcraft Super King Air 200. On-board this fixed-wing aircraft was the iNAT-RQH navigation-grade IMU from iMAR navigation along with the iTempStab temperature stabilisation box, mounted in a strapdown configuration. The complementary GNSS system consisted of a JAVAD DELTA receiver, connected to the aircraft GNSS antenna. Also on-board was a new ZLS Dynamic gravimeter (D-Type) provided by the Swedish Mapping, Cadastre and Land Registry Authority (Lantmäteriet). To the authors knowledge, this type of gravimeter was also flown for aerogravimetry for the first time.

The ZLS gravimeter platform is similar to the S-Type L&R gravimeter that DTU Space routinely operate. There are however some practical differences, one of those being that the D-type gravimeter does not have a clamping mechanism for the instrument beam. This means that the beam tends to drift during aircraft turns, resulting in substantial loss of data at the beginning of each straight line segment. In order to counteract the beam drift, we experimented with several operational procedures that counteract the drift by manually adjusting the spring tension of the instrument. After some trial-and-errors, an effective procedure was established, noting that such a procedure introduces additional operational complexity compared to beam clamping. It should also be noted that the D-Type gravimeter is designed for marine applications, which is an entirely different environment where sub-mGal accuracies are routinely obtained.

The iMAR IMU and temperature stabilisation box was constantly powered during the entire survey in order to maintain a constant temperature. The temperature reference point, about which the temperature is kept constant, is determined automatically during start-up based on the ambient room temperature. The ambient temperature is measured by a temperature sensor located inside one of the cables. However, due to some instrument communication issues, we had to re-boot the system before each flight, meaning that a new reference temperature may have been set before take-off. In the future, it will be possible for the operator to manually set the reference temperature [iMAR navigation, personal communication].

The temperature profiles for each of the six flight are shown in Figure 4, indicating that the temperature variation is within 0.4 $°$C after convergence. For some unknown reason, the initial temperature was around 44 $°$C for the last flight. The reason could be that the unit was left unpowered overnight and turned on in the morning a few hours before take-off. The effect of any residual temperature variation shown in Figure 4 was accounted for using a thermal correction curve as described in [27]. The slope of the correction curve is around 1 mGal/$°$C, meaning that the residual temperature effect is on the order of 0.5 mGal.

4. Data Processing

The gravity estimates are derived using software developed at DTU Space. A gravimeter is essentially an accelerometer measuring a combination of gravitational attraction, $\mathbf{g}$, and acceleration due to movement, $\ddot{\mathbf{r}}$, as

$$\mathbf{f}=\ddot{\mathbf{r}}-\mathbf{g},$$

(1)

where $\mathbf{f}$ is commonly known as specific force and dots denote the derive with respect to time. The ZLS and iMAR gravimeters therefore provide specific force observations, $\mathbf{f}$, whereas the GNSS system provides position estimates, $\mathbf{r}$, meaning that the gravity signal is contained within the difference between those observations.

4.1. The GNSS Data

Position estimates were derived from GNSS observations using NovAtel’s Waypoint software suite. The final satellite ephemerides for GPS and GLONASS are processed by the Center for Orbit Determination in Europe (CODE) and made available by the International GNSS Service (IGS) with a latency of 12–18 days, after which a Precise Point Positioning (PPP) solution was produced of aircraft position and associated error covariance matrix. The position estimates are introduced into the data processing as described in the following two subsections.

4.2. The ZLS Data

The processing methodology for the ZLS observations is essentially a direct approach to the expression in Equation (1) as

$$g=f_z-\ddot{h}+\delta g_{\mathrm{E}\mathsf{ö}\mathrm{tv}\mathsf{ö}\mathrm{s}}+\delta g_{\mathrm{tilt}}+\left(g_{\mathrm{tie}}-g_{\mathrm{base}}\right),$$

(2)

where the specific force, $f_z$, along the vertical axis is readily available from the instrument output and the kinematic accelerations, $\ddot{h}$, are derived from GNSS height estimates using a double-difference approach. The additional terms are the Eötvös correction [28], a tilt correction and the introduction of a tie value, by subtracting the base reading at the apron and adding the known tie value.

The tilt correction, $\delta g_{\mathrm{tilt}}$, accounts for the platform being off-level, meaning that the vertical accelerometer is not fully aligned with the gravity plumb line. The gravimeter itself provides specific force observations along a vertical axis, $f_z$, and two additional accelerometers, mounted on top of the gravimeter, provide specific force observations along the two horizontal directions, $f_x$ and $f_y$. The observations along the horizontal directions can be used to estimate the platform tilt as

$${\varphi }_x=\frac{f_x-q_x}{g}\mathrm{and}{\varphi }_y=\frac{f_y-q_y}{g},$$

(3)

where $q_x$ and $q_y$ are horizontal accelerations derived from the GNSS position estimates. For small tilt angles, the tilt correction can be approximated as a linear combination of specific force observations as ([29], Equation 2.11)

$$\delta g_{\mathrm{tilt}}=\left(1-cos{\varphi }_{x}cos{\varphi }_y\right)f_z+sin{\varphi }_x{f}_x+sin{\varphi }_y{f}_y,$$

(4)

noting that data will be discarded for large tilt angles.

Finally, the gravity estimates are filtered using a six-fold cascaded second-order Butterworth filter with a half-power point of 170 s, leading to a half-wavelength spatial resolution of approximately 7 km. This filter is applied both forward and backward in time in order to avoid the introduction of time lags.

4.3. The iMAR Data

The data processing of the iMAR strapdown IMU has been previously described in detail [18]. This approach is commonly applied for integrated IMU/GNSS navigation systems using a loosely-coupled closed-loop Kalman filter framework. Basically, one defines a state vector, $\mathbf{x}$, governing the essential modes of the system. In this case, these are the physical variables

$$\mathbf{x}=\left[{\psi }_{nb}^n,{\mathbf{v}}^n,\mathbf{p},{\mathbf{b}}_a,{\mathbf{b}}_g,\delta {\mathbf{g}}^n,\delta {\dot{\mathbf{g}}}^n,\delta {\ddot{\mathbf{g}}}^n\right],$$

(5)

where n denotes the navigation reference frame aligned with the north, east and down directions and b denotes the aircraft body frame aligned with the forward, starboard and through-the-floor directions. The physical variables are the attitude, ${\mathit{\psi }}_{nb}^n$, between those two reference frames, the Earth-referenced velocity along the n-frame directions, ${\mathbf{v}}^n$, the geodetic position, $\mathbf{p}$, in terms of latitude, longitude and ellipsoidal height, the accelerometer sensor biases, ${\mathbf{b}}_a$, the gyroscope sensor biases, ${\mathbf{b}}_g$ and the gravity disturbance vector, $\delta {\mathbf{g}}^n$, along with its first two derivatives, $\delta {\dot{\mathbf{g}}}^n$ and $\delta {\ddot{\mathbf{g}}}^n$.

In an error state implementation, a linear dynamic system model governing the temporal evolution of the errors on those variables is then defined as

$$\delta \dot{\mathbf{x}}\left(t\right)={\mathbf{F}}^n\left(t\right)\delta \mathbf{x}\left(t\right)+\mathbf{G}\left(t\right){\mathbf{w}}_s\left(t\right),$$

(6)

where $\delta \mathbf{x}$ is a vector representing the error on the state variables defined in Equation (5) and ${\mathbf{w}}_s$ is a vector of noise terms that serve as driving input to the system. The system matrix, ${\mathbf{F}}^n$, determines how the errors evolve by themselves and is defined as

$${\mathbf{F}}^n=\left[\begin{array}{cccccccc}{\mathbf{F}}_{11}^n& {\mathbf{F}}_{12}^n& {\mathbf{F}}_{13}^n& {\mathbf{0}}_3& {\mathbf{C}}_b^n& {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_3\\ {\mathbf{F}}_{21}^n& {\mathbf{F}}_{22}^n& {\mathbf{F}}_{23}^n& {\mathbf{C}}_b^n& {\mathbf{0}}_3& {\mathbf{I}}_3& {\mathbf{0}}_3& {\mathbf{0}}_3\\ {\mathbf{0}}_3& {\mathbf{F}}_{32}^n& {\mathbf{F}}_{33}^n& {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_3\\ {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_3\\ {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_3\\ {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{I}}_3& {\mathbf{0}}_3\\ {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{I}}_3\\ {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_3& -{\beta }^3& -3{\beta }^2& -3\beta \end{array}\right],$$

(7)

with the upper left elements listed in ([30], Equations 14.64–14.71) ${\mathbf{C}}_b^n$ is the transformation matrix from the b-frame to the n-frame, ${\mathbf{I}}_3$ is the $3\times 3$ identity matrix and ${\mathbf{0}}_3$ is a $3\times 3$ matrix of zeros. The term $\beta$ is a $3\times 3$ diagonal matrix containing correlation parameters defining the temporal evolution of the gravity disturbance vector. In this case, we will assume that the gravity field is constant in time, but that the covariance between two points in space can be described by a third-order Gauss-Markov model with standard deviation, ${\sigma }_{\mathrm{GM}3}$, and correlation parameter, ${\beta }_{\mathrm{GM}3}$. As the correlation parameter is defined in terms of distance, it can be converted into the time domain using the along-track velocity, ${\left|v_{\mathrm{hor}}\right|}^2=v_N^2+v_E^2$, as

$$\mathit{\beta }=\left[\begin{array}{ccc}{\beta }_N& 0& 0\\ 0& {\beta }_E& 0\\ 0& 0& {\beta }_D\end{array}\right]=\left|v_{\mathrm{hor}}\right|\left[\begin{array}{ccc}{\beta }_{N,\mathrm{GM}3}& 0& 0\\ 0& {\beta }_{E,\mathrm{GM}3}& 0\\ 0& 0& {\beta }_{D,\mathrm{GM}3}\end{array}\right],$$

(8)

with $\left(N,E,D\right)$ denoting the north, east and down directions, respectively. The two parameters of the stochastic model govern the amplitude and resolution of the resulting gravity estimates and are turned during processing as described in Appendix A. The choice of parameters is however ultimately a subjective choice. For the results presented here, a standard deviation of 15 mGal and a correlation parameter of 1/(8km) was used, indicating a correlation length of approximately 23 km [31].

The GNSS position estimates and external gravity observations are introduced in the processing as measurement updates, which is described in [18]. Finally, the estimates are smoothed by combining estimates processed forward and backward in time using an implementation of the Rauch-Tung-Striebel smoother ([32], Chapter 5). This type of smoothing makes it difficult to determine the spatial resolution in terms of wavelengths.

4.4. Merging of ZLS and iMAR Gravity Estimates

In order to examine whether the iMAR gravity estimates will still benefit from the more well-tested long-wavelength information present in the ZLS gravity estimates, two merging approaches are tested:

• Bias Adjustment: The iMAR gravity estimates are interpolated onto the time stamps of the ZLS estimates of the particular flight line. The mean value of the ZLS estimates and the interpolated iMAR estimates are determined. The iMAR mean value is subtracted from the iMAR estimates and the ZLS mean value is added instead
• Bias&Trend Adjustment: The iMAR gravity estimates are interpolated onto the time stamps of the ZLS estimates of the particular flight line. A straight line is fitted to the ZLS estimates and the interpolated iMAR estimates. The iMAR straight line fit is is subtracted from the iMAR estimates and the ZLS straight line is added instead

Thus, not all iMAR gravity estimates along the line are used for estimating bias and trend. The correction is however applied to all iMAR estimates along the line. It should be noted that ZLS estimates do not exist for all lines and are sparse for others. This means that it can be difficult to obtain a representative estimate of bias and/or trend for some flight lines. Lines with no ZLS data are not adjusted.

5. Results

Gravity estimates along the flight lines are derived differently for the two systems as described in the previous section. The ZLS gravity estimates are derived directly as the difference between accelerations, low-pass filtered with a half-power point of 170 s, while the iMAR gravity estimates are modelled as a third-order Gauss-Markov model with a standard deviation of 15 mGal and correlation parameter of 1/(8 km), corresponding to a correlation length of 23 km. In Figure 5, these two datasets are interpolated over the survey region using least-squares collocation. The covariance matrix for the interpolation was formed using a third-order Gauss-Markov model having a standard deviation of 10 mGal and a correlation parameter of 1/(15km), indicating a correlation length of 44 km.

In

Table 3. Statistics from the difference between gravity estimates at cross-over locations.

	ZLS	iMAR	iMAR	iMAR
		No Adjustment	Bias Adjusted	Bias&Trend Adjusted
no. of crossings	12	63	63	63
mean	0.3	−0.1	−0.3	−0.4	mGal
standard dev.	2.6	1.5	1.1	1.9	mGal
minimum	−3.0	−4.2	−4.0	−9.4	mGal
maximum	6.2	3.5	1.8	4.5	mGal
RMS	2.5	1.5	1.1	2.0	mGal
RMSE	1.8	1.0	0.8	1.4	mGal

is shown some statistics from the differences in gravity estimates at cross-over locations. The Root-Mean-Square-Error (RMSE) is obtained as RMS/$\sqrt{2}$ and reflects the distribution of a single stochastic variable, based on the distribution of the difference between two stochastic variables, assuming that these have identical standard deviation.

The difference between the ZLS estimates and the non-adjusted iMAR estimates are shown to the left in Figure 6. The figure to the right shows the difference between the bias adjusted and interpolated iMAR gravity estimates and the Earth Gravitational Model 2008 (EGM2008) [33]. The EGM2008 model is the most complete global gravity model available, incorporating all types of gravity observations, and is presented in Appendix A.1. Such global models are typically formulated in terms of spherical harmonic coefficients up to some maximum degree and order, characterising their spatial resolution. The maximum degree of the EGM2008 global model is 2190, indicating a half-wavelength spatial resolution of around 12 km [34]. For reference, the span of the figures shown here are around 400 km along the latitudinal direction and 200 km along the longitudinal direction.

Table 4. Mean, standard deviation, minimum and maximum of gravity estimates and global models, along with the difference between airborne estimates and external gravity information.

	Mean	Standard Dev.	Minimum	Maximum
iMAR (no adj.)	13.6	12.6	−11.0	52.6	mGal
iMAR (bias adj.)	13.7	12.8	−11.0	52.6	mGal
iMAR (bias & trend adj.)	13.6	12.9	−12.3	52.6	mGal
ZLS	13.6	11.9	−11.1	40.6	mGal
EGM2008	16.1	13.1	−9.5	52.2	mGal
GO_CONS_GCF_2_DIR_R6	15.5	11.7	−11.4	48.7	mGal
No adjustment
iMAR - ZLS	0.5	1.6	−6.9	14.5	mGal
iMAR - EGM2008	−2.5	2.2	−17.6	8.3	mGal
iMAR - GO_CONS_GCF_2_DIR_R6	−2.0	6.5	−19.4	16.5	mGal
iMAR - DTU17	−1.8	1.9	−10.1	9.1	mGal
iMAR - NKG15	1.4	1.3	−3.1	8.1	mGal
Bias adjustment
iMAR - ZLS	0.0	1.3	−5.5	6.2	mGal
iMAR - EGM2008	−3.0	2.1	−17.6	7.8	mGal
iMAR - GO_CONS_GCF_2_DIR_R6	−2.4	6.4	−20.5	16.5	mGal
iMAR - DTU17	−1.2	1.8	−10.2	10.2	mGal
iMAR - NKG15	0.9	1.3	−3.4	7.8	mGal
Bias & trend adjustment
iMAR - ZLS	0.0	1.2	−5.4	6.4	mGal
iMAR - EGM2008	−3.2	2.5	−17.0	7.4	mGal
iMAR - GO_CONS_GCF_2_DIR_R6	−2.6	6.6	−24.2	16.5	mGal
iMAR - DTU17	−1.0	2.2	−10.0	9.6	mGal
iMAR - NKG15	0.7	2.1	−10.3	11.5	mGal
ZLS platform estimates
ZLS - EGM2008	−3.1	2.2	−16.1	6.5	mGal
ZLS - GO_CONS_GCF_2_DIR_R6	−2.0	6.2	−19.3	6.5	mGal
ZLS - DTU17	−1.4	2.0	−10.3	4.0	mGal
ZLS - NKG15	0.9	1.5	−4.0	6.6	mGal

lists the mean, standard deviation, minimum and maximum values of both the iMAR and ZLS line gravity estimates, i.e., not the interpolated estimates. The same statistical parameters are shown for EGM2008 and ESA’s Release 6 GOCE gravity field model (GO_CONS_GCF_2_DIR_R6), both evaluated along the flight lines. The GO_CONS_GCF_2_DIR_R6 model is based on GRACE and GOCE satellite measurements only, meaning that this model represents the long wavelength information of the gravity field [35]. The maximum degree is 300, indicating a half-wavelength spatial resolution of around 85 km.

Also shown are statistical values for the differences between the various datasets. Additional to the global models are shown comparisons with the DTU17 model based on satellite altimetry [36] and to the gravity database of the Nordic Geodetic Commission (NKG15). These datasets represent gravity anomalies at the surface of the Earth and are upward continued to flight altitude as described in Appendix A.2 and Appendix A.3 respectively.

In addition to the datasets presented in Table 4, a bathymetric/topographic model was formed as described in Appendix A.4. Although the gravitational response from this model is very small for most of the survey area, three areas with significant response can be identified. With respect to Figure 3, these are (1) Hallandsåsen in southern Sweden, (2) the area north of Gothenburg in southern Bohuslän, Sweden and (3) Djursland in Denmark. Figure 7 shows the topographic response computed along the flight lines in these areas and the iMAR gravity disturbance estimates. In order to compare the two, one has to take into account the inherent filtering of the airborne estimates. This is done here by applying a six-fold cascaded second-order Butterworth filter in the time domain using various filter lengths. By forming the difference between the filtered topographic response and the iMAR estimates, the RMS difference can be computed as a similarity measure between the two. Except for the line X1/X2 over Hallandsåsen, this leads to local minima using filters with a half power point ranging from 80 to 150 s, indicating a half-wavelength spatial resolution of 3 to 6 km.

6. Discussion

The gravity estimates shown in Figure 5 show a consistent map of the gravity field with a large degree of across-track correlation. The iMAR estimates show a continuous field, illustrating its robustness against aircraft dynamics, whereas a large portion of the ZLS gravity estimates are discarded. Since it is difficult to interpret the spatial resolution of the iMAR gravity estimates in terms of signal wavelength, there is no direct way of comparing the resolution of the two datasets. However, the iMAR gravity map does appear more detailed and a comparison with filtered topographic response indicates a half wavelength spatial resolution in the range of 3 to 6 km.

In terms of cross-over statistics in Table 3, both the ZLS and iMAR datasets show small mean values of 0.3 and −0.1 mGal, respectively. The accuracy in terms of RMSE is 1.8 and 1.0 mGal, respectively. This indicates that the iTempStab-AddOn is significantly limiting the accelerometer drift, thereby improving the accuracy of the long wavelength information contained in the resulting gravity estimates. However, the accuracy seems to improve when adjusting the iMAR line estimates as a bias with respect to the ZLS line estimates. There is therefore indication that the long-wavelength information of the ZLS estimates is still superior to that of the iMAR estimates. Whether this remaining error originates from sensor drift, processing artefacts or some other source is unknown. Adjusting the iMAR gravity estimates using both a bias and trend from the ZLS estimates seems to worsen the accuracy in terms of cross-over statistics. This is probably because ZLS estimates on some lines are sparse, making it difficult to obtain a stable estimate of the trend.

The presence of a bias between the iMAR and ZLS gravity estimates is also indicated in Table 4, since the difference between the two datasets has a mean value of 0.5 mGal. Inspecting the gravity estimates along individual flight lines shows excellent agreement for most lines, but a bias can be identified for five particular lines. Removing these five lines reduces the mean difference to 0.1 mGal. The reason for this is however unclear.

The difference between EGM2008 and the adjusted iMAR gravity estimates is shown in Figure 6. First of all, it is evident that the airborne gravity data is more detailed. Second, regarding the more long-wavelength part, the figure shows that gravity in EGM2008 is generally too high over the ocean and too low in the southern-most part of the survey area, along the northern coast of Zealand. Such erroneous long-wavelength information will induce a tilt in the resulting geoid with an amplitude of up to 10 cm.

The long-wavelength disagreement between the airborne data and EGM2008 is also evident when inspecting the mean differences in Table 4. The mean difference improves with respect to GO_CONS_GCF_2_DIR_R6, although this model has a lower resolution, indicating that long-wavelength information of the satellite-only model is more reliable. The mean difference of −2.4 mGal with respect to GO_CONS_GCF_2_DIR_R6 is probably an effect of the large difference in resolution between the two datasets. The survey area is not large enough to be properly represented by the satellite model, having a half-wavelength resolution of more than 80 km. This theory is supported by the comparison with the satellite altimetry model DTU17 and the NKG15 terrestrial/marine gravity observations, which leads to mean differences of −1.2 and 0.9 mGal, respectively.

7. Conclusions

The cross-over statistics yields an accuracy of 1.0 mGal for the iMAR strapdown system and 1.8 mGal for the ZLS platform system. The half-wavelength resolution of the ZLS estimates is 7 km and likely better for the iMAR estimates, although no straightforward method of determining the resolution from the Kalman filter parameters is known.

The cross-over statistics also yields a mean value of -0.3 mGal for the iMAR gravity estimates which, together with the derived accuracy, indicates that the iTempStab-AddOn is significantly improving the long-wavelength information of the gravity estimates. However, a straightforward merging approach using the ZLS mean values improves the iMAR gravity estimates, indicating that the long wavelength information contained in the ZLS estimates is superior and that some erroneous component still remains.