Improving Gravity Estimation Accuracy for the GT-2A Airborne Gravimeter Using Spline-Based Gravity Models

In this article, we propose a technique for processing airborne gravimetry measurements at repeat drape lines. To model gravity along a line, we use cubic B-splines. The unknown parameters of the model (the spline coefficients) are to be estimated from the measurements. In order to take into account dependance of gravity on position, the same set of the spline coefficients is assumed for gravity modeling at each repeat line. Gravity estimation is performed with the Kalman filter. The spline coefficients are included in the state vector of the filter as unknown constants. The developed estimation algorithm was tested using measurements of the GT-2A gravimeter collected at repeat drape lines. The gravity estimates obtained with the use of B-splines are more accurate compared with those obtained with the standard technique, in which gravity is modeled as a stochastic process in the time domain.


Introduction
Airborne scalar gravimetry is widely used to obtain high accuracy information about the Earth gravity in hard-toreach regions (shelves, mountains, polar areas, etc.). A typical platformed airborne gravimetry system consists of a gravity sensor (a single-axis accelerometer) mounted on a gyro-stabilized platform and two (or more) global navigation satellite system (GNSS) receivers (one on board the aircraft, others on the ground). Postprocessing of airborne gravimetry measurements includes three main stages such as determination of GNSS positioning and velocity solutions, estimation of the platform tilt angles via the inertial navigation system (INS) and GNSS integration, and estimation of gravity using the Kalman filtering technique or low-pass filtering. Both filtering techniques are based on representing gravity in the time domain (e.g., (Lu et al. 2017) or modeling gravity as a stochastic process (Bolotin and Golovan 2013; Kwon and Jekeli 2002;Stepanov et al. 2015)).
In this research, we focus on the last stage of postprocessing (i.e., gravity estimation) given airborne measurements at repeat drape lines. Using a stable platform gravimetry system (e.g., the GT-2A (Berzhitzky et al. 2002)) in a drape flight is challenging, since such a system is more suitable for benign dynamics (Studinger et al. 2008). For the GT-2A, underestimation of the scale factor of the gravity sensor can be observed during postprocessing of drape lines. For this reason, specific calibration lines are flown, at which the scale factor is reliably estimated. However, accuracy of gravity estimation at repeat drape lines is still lower than in the case of a flight at a constant height. To our opinion, the standard approach to gravity estimation based on modeling gravity in the time domain is not well suitable for processing repeat lines, as gravity modeled along a line, strictly speaking, differs from gravity modeled along a repeat line. Actual gravity is, of course, the same at repeat lines assuming that they were flown at the same altitude. Therefore, in the case of a repeat line flight, it is preferable to model gravity in the spatial domain rather than in the time domain.
In this work, we propose an approach based on spatial modeling of gravity to process repeat drape lines flown over the same ground track. We parameterize the gravity disturbance along a line using B-splines of order three. To take into account dependance of gravity on position, we use the same set of unknown parameters of the model (the spline coefficients) for gravity modeling at each repeat line. The model is easily described in the state space assuming that the spline coefficients are constant over the time during the aircraft's flight. The standard stochastic estimation problem is posed, and solved via Kalman filtering. The state vector for the filter includes the spline coefficients and the parameters of the gravimeter systematic errors (the calibration parameters).
The developed approach was applied to processing airborne measurements collected with the GT-2A system during a repeat drape line flight in South Africa in 2009. The obtained gravity estimates are compared with those provided by the standard approach based on modeling gravity in the time domain, and with the gravity values derived from the global model EGM2008. An increase in accuracy in gravity estimation using the new approach is shown.

Basic Equation and Gravimeter Measurements
The gravity sensor of the GT-2A is mounted on a gyrostabilized platform, which makes the instrument sensitivity axis to be vertical. The gravity disturbance ıg is determined from the fundamental equation of airborne scalar gravimetry, which is written as: where V 3 is the vertical velocity of the gravimeter proof mass in the ENU coordinate frame, g 0 is the normal gravity at the flight height, g Etv is the Eötvös correction term, f 3 is the vertical projection of the specific force.
Raw measurements of the gravimeter f 0 3 can be described by the equation Here k 3 is the scale factor error, is the time delay of gravimeter measurements (due to physical delay and averaging of gravity sensor measurements performed in the internal software of the gravimeter), ıf hor is the horizontal acceleration influence on f 0 3 due to imperfect placement of the gravimeter on the platform, i.e., more explicitly, ıf hor D k 1 f 0 1 C k 2 f 0 2 , where f 0 1 , f 0 2 are the platform accelerometer measurements, k 1 , k 2 are unknown factors that can differ from one flight to another. Further, ıf tilt is the horizontal acceleration due to the platform tilt error (it is determined from INS-GNSS integration), ıf drif t is the sum of the gravity sensor bias and drift. Assuming linear drift, the term ıf drif t can be computed given pre-and post-flight measurements and terrestrial gravity value. The term q f is the noise.
The frequency of the GT-2A raw measurements is 18.75 Hz, and the time delay is about 2 s.
As the terms ıf tilt and ıf drif t can be determined and compensated for using standard procedures, we will omit them further on.
State-Space Modeling Rewrite Eq. (1) replacing unknown variables on the right-hand side of the equation by measurements and omitting the gravity disturbance term: where g 0 and g Etv are known from the GNSS position and velocity solutions (assuming that the lever arm effect is compensated for).

Denote by V 3 the vertical velocity error defined as (Bolotin and Golovan 2013)
Then we obtain the error dynamics equation of airborne scalar gravimetry by subtracting Eq. (2) from Eq. (1) and using the above gravimeter measurement model where k 1 ; k 2 are the residual platform misalignment errors in the sum with errors of placement of the gravity sensor on the platform, and k 3 is the scale factor error of the gravity sensor.
Representing the GNSS velocity solution as V GNSS 3 D V 3 C e vel , where e vel is the GNSS velocity random error, and introducing the velocity observation as y D V GNSS 3 V 0 3 , we obtain the observation equation Equations (4) and (5) contain six unknown variables: V 3 ; ıg; ; k 1 ; k 2 ; k 3 : The variables q f and e vel in Eqs. (4) and (5) are assumed to be zero-mean white noises with known variances. We adopt the following models for the gravimeter systematic error parameters (calibration parameters): where q ; q k 1 ; q k 2 ; q k 3 are zero-mean white noises with known (sufficiently small) variances. By introducing an equation for ıg (a gravity model), we obtain the system of equations in the state-space form. The gravity estimation algorithm includes the following steps: 1. Determining V 0 3 from Eq. (2). 2. Estimation of gravity and the gravimeter calibration parameters via Kalman filtering and smoothing given the state-space system Eqs. (4)-(6) and the equation describing a gravity model.

Modeling Gravity Along Repeat Lines Using B-Splines
To arrive at estimation algorithm, we introduce a spatial gravity model assuming a repeat line flight. We use Bsplines of order three for gravity modeling, which form a basis in the space of cubic splines in R 1 and each function of which has finite support ( Fig. 1). Recall that the B-spline, by definition, is a fourfold convolution of a rectangular function with itself (De Boor 2001). We assume that the repeat lines are flown along the same ground track at the same altitude and have equal length L. Let x D x.t/ 2 R 1 denote the distance travelled by the aircraft (more precisely, the gravity sensor proof mass) along the first line. Then x varies from 0 to L. For simplicity, assume that the subsequent repeat lines coincide with the first line. Thus, x varies from 0 to L or from L to 0 at the repeat lines, depending on the direction of aircraft's flight.
We represent the gravity disturbance at a point of a line as follows where B i .x/ is the cubic B-spline, c i is the i -th spline coefficient to be estimated, N is the total number of the spline coefficients. Note that the same spline coefficients c i , i D 0; : : : ; N 1, are used for describing gravity at each repeat line. Note also that the B-spline B i .x/ in fact is the function of .x x S i /= x, where x S 0 ; : : : ; x S N 1 are the knots located at the first line, x is the length of the knot step (measured along the line). The value of x is to be chosen.
Assuming that the spline coefficients are constant over the time during the flight, we can easily transform the gravity model Eq. (7) into the state-space form by writing the obvious equations P c i D 0; i D 0; : : : ; N 1: Therefore, we obtain the complete state-space system described by Eqs. (4)-(6) and (8) given the gravity model (7). Using Kalman filtering, the optimal (in the minimum mean squared error sense) estimates of the spline coefficients, c i , i D 0; : : : ; N 1, the velocity error V 3 and the calibration parameters, ; k 1 ; k 2 ; k 3 , are obtained. The total dimension of the state vector is N C 5. At the initial iteration of the filter, the a priori variance of the state vector, including the variance of the spline coefficients, should be chosen.

Transfer Function of the Optimal Filter
The Fourier transform of a cubic B-spline can be written as (De Boor 2001): where T D x=V is the time step of the knots, V is the average speed of the aircraft at the lines, j D p 1. By varying the time step of the knots T , one can obtain B-splines with different frequency bandwidths. Since then, it is desirable to determine how the spatial resolution of the gravity estimate provided by the Kalman filter depends on T .
For this, let us derive the expression for the transfer function of the filter. Consider the state-space equations Eqs. (4)-(8) in the discrete-time form denoting by t the measurement time step. Neglect for simplicity the gravimeter calibration parameters and the gravity sensor noise q f in the equations. Let the measurement time step t and the spline knot step T satisfy the relationship T D m t for some integer m 2 and assume no a priori information on the variances of the spline coefficients (i.e., the variance of each coefficient tends to infinity). As above, assume that the GNSS velocity error that appears in Eq. (5) is a white noise. Then it can be shown that the optimal filter transfer function which maps the acceleration observations to the Fig. 2 Transfer functions of the optimal filter in the spline-based approach and the standard approach based on the Gauss-Markov model of order 2 for gravity, and the Butterworth filter (two-pass) of order 2 gravity disturbance estimate can be written as (Unser et al. 1993) H .z/ D B The obtained transfer function is close to the Butterworth filter (two-pass) of order 2 (see Fig. 2). Figure 2 also shows the transfer function of the optimal filter in the standard approach to gravity estimation, which is based on stochastic modeling of gravity (e.g., using a Gauss-Markov process of order 2 with infinite correlation time, that is, a process whose second time derivative is a white noise). In this case, the transfer function of the filter is close to the Butterworth filter (two-pass) of order 4. Thus, with the use of the above expression for the transfer function of the optimal filter (for the spline-based approach), the relationship between the filter cutoff wavelength g and the length of the spline knot step x (or the time step of the knots T ) can be derived as g 3 x; x D V T :

Survey Flight Overview
The developed approach based on spatial gravity modeling using B-splines was applied to airborne gravimetry data collected within a flight over the Vrederfort Dome impact crater (120 km southwest of Johannesburg, South Africa) on March 21, 2009. The GT-2A gravimetry system used in this survey was  (Olson 2010). Two dual-frequency GPS receivers operating at 10 Hz were used during the flight (the rover on board the aircraft and the base placed on the ground). The flight consisted of six repeat drape lines (with the ID labels R01-R06) and two calibration lines (with the ID labels R07-R08). The latter were flown in order to obtain more accurate estimate of the gravimeter scale factor error k 3 . All the lines were flown in the north-south and southnorth direction over the same ground track (Fig. 3). The drape lines were flown at the heights from 1,400 to 1,600 m above the WGS-84 ellipsoid and accurately followed the terrain (Fig. 4). The calibration lines were flown at the height of 3,000 m (Fig. 5).
The length of each line is 85 km. Average flight speed is 57 m/s with standard deviation (STD) 1.3 m/s. The flight duration is 3 h 50 min. Accuracy of Estimating the Scale Factor Error k 3 The maximum vertical acceleration at the repeat drape lines (with the Eötvös effect and the normal gravity subtracted) equals 0.1 m/s 2 (Fig. 6). Hence, for gravity estimation with the accuracy of 0.5 mGal the scale factor error k 3 should be estimated with the accuracy of 0.00005 or better.   Table 1 shows the calibration parameter estimates obtained with the developed spline-based approach from processing the drape lines and the calibration lines. The corresponding standard deviations of the estimate errors provided by the Kalman filter are shown in Table 2. Equal estimation accuracy of 0.00001 for k 3 is achieved at the drape lines and at the calibration lines. This accuracy is sufficient for the gravity estimation error due to k 3 to be at the level of 0.5 mGal. Moreover, since the two estimates of the scale factor error k 3 are close enough to each other (the  Fig. 7 The along-track gravity disturbance estimates obtained with the two approaches for the cutoff wavelength of (a) 100 s and (b) 70 s, and gravity from EGM2008; (c) GPS height (averaged over six repeat drape lines) difference equals 0.00005), it is sufficient to process only the repeat drape lines to obtain a reliable estimate of k 3 , and there is no necessity to use the calibration lines. It should be noted here that the standard approach, which is used in the GTGRAV software and based on modeling gravity in the time domain, provides a reliable estimate of k 3 only from processing the calibration lines.

Discussion
The gravity disturbance was estimated with the new approach based on spatial modeling of gravity using B-splines and with the standard approach based on modeling gravity in the time domain. Figure 7 shows the along-track gravity disturbance estimates obtained for the filter cutoff wavelength (inverse cutoff frequency) of 100 and 70 s, which are equivalent to the spatial resolution g =2 of 3 and 2 km, respectively. From Eq. (9), the corresponding values of the length of the spline knot step x are 2 and 1.4 km. These values resulted in N D 46 and N D 64 spline coefficients in the gravity model Eq. (7), respectively. Figure 8 shows how the estimates of the spline coefficients evolve during processing of the repeat drape lines. The estimate changes occur due to the measurement update step of the Kalman filter. Since B-splines have finite support in the time domain, the estimate changes are local. As the spline Mean value of the difference of the gravity estimates obtained with the two approaches (100 s cutoff wavelength) is 0.05 mGal, the STD is 1.13 mGal. It can be seen from Fig. 7 that the gravity estimates provided by the new algorithm are more detailed and better correspond to the height profile in the interval from 60 to 70 km.
In the absence of ground gravity data, we compare the gravity estimates with the global gravity model EGM2008, which is sufficiently accurate in the area of interest, as dense gravity datasets were available at the model construction stage. Table 3 summarizes the statistics for comparing the gravity estimates and the gravity values derived from EGM2008. It can be seen that both approaches (spatial modeling and modeling in the time domain) show the same result for the cutoff wavelength of 180 s. For 320 s (equivalent to the nominal spatial resolution of EGM2008, which equals 9 km), the new approach showed greater discrepancy than the standard one. The reason is likely in short-wavelength errors present in the gravity estimate provided by the new approach, as the optimal filter has less steep slope than in the case of the standard approach (see Fig. 2).
At last, the gravity estimates were computed by processing a different number of repeat drape lines (from 1 to 6) to analyze how each new repeat line affects the gravity estimate. The estimates were obtained with the two approaches for the cutoff wavelength of 100 s and compared to EGM2008. The results are shown in Fig. 9. The dependance on the number of the processed lines is different for the two approaches. For the new approach, the discrepancy with EGM2008 grows with the number of processed lines, which can be explained by the more and more detailed along-track estimate of gravity. For the standard approach, on the contrary, the along-track gravity estimate obtained by averaging over the lines becomes less detailed with the number of processed lines, and tends to gravity from EGM2008.  9 Standard deviations of the difference of the gravity estimates obtained for the cutoff wavelength of 100 s and gravity from EGM200

Conclusions
In this research, an approach to gravity estimation at repeat drape lines was developed using B-splines for spatial modeling of gravity. This approach was compared with the standard one based on modeling gravity in the time domain by processing the GT-2A gravimeter data. The spatial gravity modeling showed an increase in accuracy of estimating the gravity sensor scale factor significantly. The achieved accuracy was sufficient for gravity estimation along the repeat drape lines with the accuracy of 0.5 mGal. Moreover, the calibration lines, which are necessary for accurate estimation of the scale factor when using the standard approach, are not required in the case of the developed approach. The reliable estimate of the scale factor was obtained from processing the repeat drape lines. The gravity estimate provided by the new approach is in good agreement with EGM2008 (for the filter cutoff wavelength close to the minimum wavelength represented by the global model). For the shorter cutoff wavelengths, the new approach provides more detailed gravity estimates than the standard one.