New Algorithm for Gravity Vector Estimation from Airborne Data Using Spherical Scaling Functions

The paper presents an approach to determination of the gravity disturbance vector from airborne gravimetry measurements at an aircraft’s flight path. A navigation-grade inertial navigation system (INS) and the carrier-phase differential mode of GNSS are assumed. To improve observability of the gravity horizontal components, which are observed in combination with the INS systematic errors, we use a spatial model of gravity. We parameterize the disturbing potential in the observation area using the spherical scaling functions. The unknown coefficients of the parameterization and the INS systematic errors are estimated simultaneously using the Kalman filter. Due to ill-conditioning of the estimation problem, the information form of the Kalman filter and regularization are used. The numerical results obtained from simulated data processing shows that the approach based on spatial modeling is capable to improve accuracy of the gravity horizontal component determination comparing to a typical modeling of gravity in the time domain.


Introduction
Airborne gravimetry is capable to provide Earth's gravity observations of high accuracy and spatial resolution. A typical airborne gravimetry system includes an inertial navigation system (INS) (strapdown or a gyro-stabilized platform) and global navigation satellite system (GNSS) receivers operating in the carrier-phase differential mode. Airborne gravimetry is well established for determining the vertical component of the gravity disturbance vector (scalar gravimetry based on platformed INS) and is capable to provide all three gravity vector components simultaneously from airborne measurements (vector gravimetry) (Kwon and Jekeli 2001;Becker et al. 2016 in vector gravimetry is based on INS/GNSS integration via Kalman filtering. A high accuracy (navigation grade) strapdown INS is assumed. The main difficulty is in determining the gravity horizontal components (deflections of the vertical) as these are observed in combination with the INS systematic errors (the inertial sensor biases, the attitude errors, etc.), which have a similar spectral content. An additional information (e.g., an a priori gravity model, external gravity data) is required for separation. The common approach to airborne vector gravimetry is based on modeling gravity in the time domain during Kalman filtering, e.g., assuming each gravity component to be a stationary stochastic process (Becker et al. 2016;Jensen et al. 2019). As an alternative, a gravity model may not be used explicitly, and the gravity components can be determined from the Kalman filter residuals (Kwon and Jekeli 2001). However, remains of the systematic errors may be present in the gravity estimates of both approaches, and additional corrections should be applied.
In this study, we develop an approach that takes into account spatial correlation of gravity at adjacent survey lines to improve observability of the gravity horizontal components. The approach is based on spatial modeling of the disturbing potential in the observation area using the spherical scaling functions (SSF) of a certain resolution level, which are harmonic and spatially localized (decrease with distance from their origins) (Freeden and Michel 2004). The SSFs and other spherical radial basis functions are widely used in local and regional gravity field modeling and are mainly applied for processing satellite gravity data and combining gravity data, e.g., Lieb et al. (2016) and Klees et al. (2008). We use the SSFs to airborne vector gravimetry problem solving. The developed estimation algorithm is based on simultaneous estimation of the INS systematic errors and the scaling coefficients (the coefficients of the parameterization of the disturbing potential) via Kalman filtering given measurements at the aircraft's flight path. The information form of the Kalman filter (Kailath et al. 2000) and regularization are applied due to the ill-conditioning of the estimation problem. The approach was partially presented in Bolotin and Vyazmin (2016).
The aim of this work is to investigate the performance of the developed approach by processing simulated airborne data and by comparing with the common approach based on modeling the gravity vector components in the time domain. The numerical results obtained with the developed approach are optimistic and show that an increase in the accuracy of the gravity horizontal component determination can be achieved.

Airborne Vector Gravimetry Equations
The gravity disturbance vector determination is based on INS/GNSS integration. The basic equation is the equation of motion of the INS proof mass M expressed in the navigation (local-level) frame M x as Here v x is the 3 1 vector of the aircraft velocity relative to the Earth expressed in M x frame coordinates, ıg x is the gravity disturbance vector, x is the reference (normal) gravity vector, f z is the specific force expressed in coordinates of the body frame (denoted by M z). Further, L zx is the transformation matrix from the navigation frame (M x) to the body frame which satisfies the kinematic equation (Farrell 2008), ! x x is the (absolute) angular velocity of the navigation frame with respect to the Earth-centered inertial frame, expressed in M x, ! x Á is the Earth absolute angular velocity vector expressed is a skew-symmetric matrix formed by the components of the vector ! x Á in such a way that where means the cross product.
From Eq.
(1), we arrive at the INS mechanization equations (Farrell 2008) by replacing true values in Eq. (1) with measurements from the INS sensors (measurements of the body-frame angular velocity and the specific force) and with GNSS position and velocity, and omitting the gravity disturbance vector term. Here, we assume GNSS positioning and velocity solutions to be obtained using the carrierphase differential mode. In this case, inaccuracies in GNSS positions can be neglected (Vyazmin and Bolotin 2019). Denote by v 0 x the velocity solution obtained from solving the INS mechanization equations and by L 0 zx the transformation matrix from the navigation frame to the computed body frame.
By subtracting Eq. (1) from the INS mechanization equations, the INS error dynamics equations can be obtained as follows (Farrell 2008; Vyazmin and Bolotin 2019) where˛x is the attitude error, ıv x is the velocity error defined as the sum of the total velocity error v x D v 0 x v x and the term O v x˛x . Note that the equations do not contain the INS positioning solution as we assume that the INS proof mass position is perfectly determined from GNSS.
The terms !, f are the gyroscope and accelerometer errors, respectively, for which the following models are assumed in the study: where b ! , b f are the gyroscope and accelerometer biases, respectively, the vectors q ! , q f , " ! , " f are the random errors. The error models in Eq. (5) represent the sensor bias drifts. Denote by v GNS S x the GNSS velocity of the aircraft relative to the Earth and by y x the observation computed as The observation equation is written as (Vyazmin and Bolotin 2019) where e v is the GNSS velocity error.

Local Gravity Model and Estimation Algorithm
In order to use a Kalman filter framework, we introduce a model of the gravity disturbance vector. We use a local spatial model by parameterizing the disturbing potential at the points at the aircraft's flight path by the SSFs (Freeden and Michel 2004;Klees et al. 2008) where a i , i D 1; : : : ; N , are the unknown parameters (the scaling coefficients) to be estimated from measurements, z i is a knot of a regular grid on a reference sphere˝R of radius R, r is a geocentric position vector, jrj R, j is the resolution level, which is defined by the spatial resolution of measurements. The number of the scaling coefficients N depends on the density of the regular grid and on the size of the modeling area (see Fig. 1) as the SSFs are spatially localized (decrease with distance from r).
The SSF is harmonic outside the sphere˝R and expressed via the Legendre polynomials P n as (Freeden and Michel 2004) where r D jrj, V z i , V r denote the unit vectors, b n j is the Legendre coefficient, which determines the spectral behaviour of the SSF. We use the Abel-Poisson SSF, which is defined by b j ;n D b n j , b j D e 1=2 j . The Abel-Poisson SSF has a closed-form expression (Freeden and Michel 2004), so there is no necessity to calculate the sum in Eq. (8). Figure 2 shows the spectral content of the Abel-Poisson SSFs of different resolution levels j . A greater value of j corresponds to a larger bandwidth of the SSF in the spherical harmonic domain and therefore to a higher spatial resolution of the gravity model Eq. (7).
The gravity disturbance vector can be represented via the scaling coefficients by differentiating Eq.
The Kalman Filter in the Information Form Let us assume that the random errors q ! , q f , " ! , " f , e v in Eqs. (4)-(6) are zero-mean white noises with known variances. Transforming the state-space model Eqs.
(2)-(6), (9)-(10) to the discretetime form, we can further use a Kalman filtering framework. The system's state vector x k at the time moment t k , k D 0; : : : ; K, where K is the total number of observations, includes N C 12 unknown variables ıv x ;˛x; b ! ; b f ; a 1 ; : : : ; a N : Using Kalman filtering, an optimal (in the sense of the minimum mean squared error) estimate of the state vector is obtained.
Since the state vector x k at each time moment t k includes the scaling coefficients for the entire observation area, the estimation problem is ill-conditioned. For this reason, we use the Kalman filter in the information form (Kailath et al. 2000;Park and Kailath 1995). The filter provides the information matrix Q k D P 1 k and the estimate of the information vector u k D P 1 k x k at each time moment t k given the initial values: Q 0 D 0, u 0 D 0. Here P k denotes the covariance matrix of the state vector estimate error. The estimates of the scaling coefficients are obtained at the last time moment t K from the state vector estimate O x K , which is obtained by solving the equation Here, the information vector estimate O u K and the information matrix Q K are provided by the filter at the last time moment t K . As the information matrix Q K is ill-conditioned (due to extension of the observation area to reduce edge effects as shown in Fig. 1, and due to the inverse problem of determining the disturbing potential from the gravity disturbance vector), regularization is required. Using Tikhonov regularization with a unit matrix I , we finally obtain the state vector estimate and the covariance matrix of the estimate error, respectively, as where > 0 is the regularization parameter.

Data Simulation
The developed approach was applied to processing simulated airborne data. Error-free measurements of the inertial sensors and error-free GNSS data (positions and velocities) were simulated using the software from Bogdanov and Golovan (2017). The simulated aircraft's flight path contained 4 lines flown in the south-north direction (Fig. 1) at a constant height of 1,600 m above the reference ellipsoid. The initial geodetic latitude is 0.3 ı . The line spacing is about 9.5 km, the length of each line is 100 km. The ground speed at each line is constant and equals 100 m/s. The attitude angles are constant at each line (the pitch and roll angles equal 0 ı ). The gravity disturbance vector at the aircraft's flight path is synthesized using EGM2008 up to d/o 1,500 (equivalent to minimum wavelength of about 26 km).
The inertial sensor errors with characteristics shown in Table 1 were simulated and added to the error-free data. The gyroscope errors were obtained from measurements of the calibrated FOG gyroscopes by Optolink recorded during a static test at constant temperature. A small constant bias The accelerometer measurement error was simulated as a sum of a zero-mean white noise and a bias drift modeled as a random walk process (see Table 1). The constant bias of 2 mGal was added to the horizontal accelerometer measurements. The bias of the vertical accelerometer was assumed to be determined and removed by comparing preand post-flight static measurements with a terrestrial gravity value.
The simulated INS measurement frequency is 100 Hz. The GNSS velocity error e v is modeled as the time derivative of a zero-mean white noise at 20 Hz, the standard deviation (STD) of e v is 2-3 cm/s.

Data Processing
The data processing procedure consisted of four steps. First, the INS mechanization equations were solved using data from the entire flight. The gravity vector x was associated with the reference gravity field represented as the sum of the normal gravity field and a longwavelength part of the anomalous field provided by a global gravity model (EGM2008 up to d/o 100 was used). The residual gravity disturbance vector (i.e., the gravity disturbance vector with a long-wavelength part subtracted) to be estimated varied from 15 to 15 mGal for the horizontal components and from 20 to 10 mGal for the vertical component.
Second, the spatial model of the residual disturbing potential was designed using Eq. (7). The resolution level of the Abel-Poisson SSFs was chosen equal to 8. The 10 km 4.5 km latitude-longitude regular grid of the scaling coefficient knots was used (Fig. 1). The total number of the scaling coefficients N equals 154. The dimension of the state vector in the Kalman filter is 166.
Third, the scaling coefficients of the gravity model and the INS systematic error parameters were estimated using the Kalman filter in the information form and regularization Eq. (11) applied at the last time moment t K . The value of the regularization parameter was set equal to s 1 10 14 , where s 1 is the maximum singular value of the information matrix Q K , which was found to be suitable.
Finally, the scaling coefficient estimates obtained at the last time moment t K were used to build the along-path estimates of the three components of the residual gravity disturbance vector using Eq. (9), and the long-wavelength part of the gravity disturbance vector was readded to the along-path estimates (EGM2008 up to d/o 100). The final

Numerical Results
The statistics for the errors of the alongline gravity estimates are given in Table 2. The standard deviation varies from 0.6 to 1.9 mGal for the east component and from 0.8 to 1.3 mGal for the north component. The average STDs are 1.2 and 1.0 mGal for the east and north component, respectively. The mean values are less than 1.7 mGal for the east and 0.7 mGal for the north component. A slightly lower accuracy of the east component, in our opinion, is due to the south-north direction of the lines and the small total number of the lines. It should be noted here that in the case of a smaller number of processed lines, the estimation accuracy for the east component is significantly worse than for the north component. The accuracy of the along-line estimates of the vertical component is better than 0.5 mGal (1 ), the mean value is not greater than 1.5 mGal.  Fig. 3 and Table 3. The accuracy of the horizontal component estimation is significantly lower than in the case of the spatial gravity modeling. Large biases and drifts can be observed in the gravity estimates. The average STDs are 3.1 and 2.8 mGal for the east and north component, respectively.
For the vertical component, the estimation accuracy (1 ) is 0.3 mGal that is better than in the case of the developed approach. The better accuracy is probably due to a better approximation of the band-limited gravity signal using the Fourier expansion than using the non-bandlimited SSF. The mean values of the vertical component estimate errors are approximately the same for the two approaches.

Discussion and Conclusions
An approach based on spatial modeling of gravity using the spherical scaling functions was developed for the airborne vector gravimetry problem solving. The approach assumes the use of a navigation-grade INS and the carrier-phase differential mode of GNSS. From the results of simulated data processing, the new approach provides significantly more accurate estimates of the gravity horizontal components in comparison with the standard approach based on modeling gravity in the time domain. The results, however, may be too optimistic as simplified models were used for simulating the GNSS errors and the inertial sensor errors (only the bias drift and noise were simulated). Further, a long-wavelength part of the gravity vector was assumed precisely known from a low-resolution global gravity model. Despite this, the developed approach, in our opinion, can be put into practice in the near future. There are, however, a number of issues, the study of which is of interest and probably will allow to improve accuracy of the gravity estimates. Namely, it should be investigated in more detail how the accuracy of gravity estimation depends on inaccuracies in the inertial sensor calibration results (e.g., on the bias level or a long-term drift in the horizontal accelerometers). At the same time, possible ways to improve the gravity estimation accuracy should be also studied, e.g., using a greater number of parallel survey lines, cross-lines, or several overlapping flights. In addition, more careful design of the spatial gravity model (e.g., selecting the type of the SSF, regularization parameter value, etc.) may also improve the resulting estimates.