Abstract
Airborne scalar gravimetry, a kinematic survey technology, is one of the most efficient techniques to acquire the gravity data in the areas where it is neither practical nor possible to make terrestrial measurement. Recently, studies have shown that the precision reaches sub-mGal. Besides the improvement of the instruments, the data processing and gravity anomaly solution algorithms are evolving consistently. This paper investigates an approach based on developing a state space model and using Kalman filtering and smoothing to determine gravity anomaly. The state space model is developed based on the gravimeter measurement equations and measurement errors equations for the specific airborne gravimetry system. The proposed method is implemented to the airborne data collected in southeast China from GT-1A with a stabilized physical platform. The solution results based on the presented model using Kalman filtering and smoothing confirm that the approach was able to implement the solution and acquire the gravity anomaly information, and the comparison between the proposed method and the traditional method indicates that a remarkable improvement in the solution precision is achieved when the proposed method is used.
1. Introduction
Airborne gravimetry is a kinematic survey technology for obtaining near-earth air gravity field information. Its theory is based on Newton’s second law of motion expressed in either inertial, Earth, or local-lever coordinate frames. The main objective of airborne scalar gravimetry is to determine the vertical component of gravity disturbance vector at flight altitude from inertial turn rate and acceleration data along with the global navigation satellite system (GNSS) position or velocity solution [1].
Despite the fact that several reformatory systems were devised and experimented, there are two kinds of airborne scalar gravimetry known as stabilized platform system and strapdown system [2]. The former one is conducted using gravimeters on a Schuler-tunes stabilized platform, and it was successfully implemented for scalar airborne gravity survey in the nineties. Although the system could be less portable because of the platform, this can be mitigated with complex design on dedicated instrumentation and the vertical orientation of gravimeters can be kept stably. The accuracy of 3–5 mGal was achieved with 10 km resolution in the initial period [3–5]. Along with the improvement in instruments and processing algorithms, sub-mGal precision can be achieved for airborne scalar survey in the southeast China [6]. The other one is conducted using integrated navigation systems which are known as GNSS/INS, and it is considered to a simple design due to the absence of the stabilized platform. This system was presented in the last eighties, though the successful test was reported in the nineties due to the availability of precision, velocity, or altitude provided by the phase and Doppler data from GNSS [7–9]. According to the experiments by Gerlach et al. and Li [10, 11], the accuracy of about 3 mGal can be achieved due to low speed of aircraft and smoothing by low-pass filters.
In addition to the instrument accuracy, the data processing algorithms are important and necessary for gravity anomaly solution. Here, the solution is regarded as the calculation process to convert the raw airborne gravimeter data into a set of free-air gravity anomaly. There are two distinct approaches used in gravity anomaly solution. In the first approach, the process starts with subtracting the sensed gravity accelerations from GNSS/INS data in the same coordinate frame as well as other effects determined by specific equations. Then, the results are smoothed by low-pass filter for suppressing the noises distributed in the high-frequency spectrum [12, 13]. Hence, the design of the filter and its parameters (order and cut-off frequency) are important to the gravity anomaly solution. The approach based on second-order Butterworth filter, which is combined with exponential filter, was designed to obtain gravity anomaly information [14, 15]. However, this approach does not provide full separation between gravity anomaly signal in low-frequency part and the noise because the boundary between gravity signal and noise in spectrum is not distinct. It has been proved that the performance of FIR (Finite Impulse Response) low-pass filter designed by optimal algorithm is superior to that of the FIR filters designed by other algorithms in gravity anomaly solution [16]. Moreover, some literatures presented the approach based on FIR low-pass filter designed by window function to determine the gravity anomaly [17–19]. Kinds of window functions were proposed such as Rectangular and Hamming to take into account performance in obtaining gravity anomaly. An interesting approach presented by Zhang and Xu is based on zero-phase filter to ensure the consistency of the gravity anomaly and measurement time [20, 21]. But the approach based on FIR filter leads to the loss of measurement data due to the boundary effect.
On the other hand, since the Kalman filtering is used to estimate the random errors of inertial instruments and the gravity disturbance is often considered as a kind of measurement noises in GNSS/INS postprocessing, the approach of gravity anomaly solution based on developing the state space model for Kalman filtering and smoothing to implement the solution was developed. This approach integrates relevant measurements and quantities by Kalman filter, and the desired gravity disturbance is modeled as a stochastic process and it is added to the state vector as estimated states. The approach has been implemented in a few researches with new processing strategies. Stepanov and Koshaev presented the accuracy analysis of using various GNSS measurements and the altitude data was suggested for solving the filtering and smoothing problem [22]. Considering the solution precision depending on the prior stochastic model of gravity, modeling parameters for stochastic process of gravity were identified by Stepano et al. [23]. The joint problem of the parameters identification and gravity anomaly solution was formulated as a nonlinear adaptive filtering problem, and these were confirmed by simulation results. Besides, Bolotin et al. presented the general statements and the approach of measurement error compensation for GT-1A airborne gravimetry system [24]. Then, a simplified gravity anomaly model described as linear state equations was developed by Bolotin [25]. An innovative method reported by Skaloud et al. combines the satellite observations and those of strapdown inertial system, and the solution was conducted using dynamic network with no additional low-pass filtering [26]. Actually, this method is regarded as a simple solution method due to absence of low-pass filtering to the Kalman filter/smoother outputs, as well as differential process for vertical kinematic acceleration. Last but not the least, all available information of the whole system can be combined for optimal estimation.
In this paper, the problem of gravity anomaly solution using Kalman filtering and smoothing is solved by the development of specific state space model. It should be noted that the Kalman filtering and smoothing used in gravity anomaly solution has been discussed before, for example, in Zheng et al.’s work, in which adaptive filtering was used to take into account the uncertainty of the altitude measurement [27]. Here, the presented state space model is developed for specific scalar airborne gravimetry system called GT-1A, and the state vector involves the variables of the gravimeter calibration.
The objective of this paper is to investigate the effectiveness and accuracy of the presented approach based on the state space model using Kalman filtering and smoothing in airborne gravity anomaly solution. First, the basic equations of airborne scalar gravimeter are presented. Then, the state space model development based on combining gravimeter measurement equations with the basic equations is given, and the altitude measurement equation is used to obtain the measurement model for filtering. Finally, practical surveying data provided by GT-1A are processed using Kalman filtering and smoothing based on the presented approach and the results are discussed and compared with those of the FIR filtering approach.
2. Basic Equations of Airborne Gravity Anomaly Solution
According to Newton’s second law of motion, the vertical dynamic equation for the sensitive mass (SM) of an airborne gravimeter can be expressed as follows [25]:where is the vertical acceleration of the gravimeter, which is considered to be the second derivative of the altitude; f3 is the vertical component of the specific force working on the SM; is the Eotvos correction; is the normal gravity and free air correction, in which is the geographical latitude and is the altitude above the reference ellipsoid; is the desired free-air gravity anomaly.
The Eotvos correction can compensate the error caused by the centrifugal force produced by Earth’s rotation and aircraft motion. It can be determined as follows [25]:where is the velocity of Earth’s rotation; and are the eastern and northern velocities, respectively; and are the curvature radii in the eastern and northern directions, respectively. The normal gravity and free-air correction can be determined by Helmert’s equation [25]:
In equation (3), m/s2, , , m/s2, and s−1.
Equation (1) can be transferred as
It can be seen that in equation (4) and can be determined using the position and velocity, and is the gravimeter reading along the survey line processed by essential corrections and error compensations. Therefore, the gravity anomaly solution can be conducted by solving equation (4) using the gravimeter and navigation measurements.
Note that the positioning observations (latitude, longitude, and altitude) of SM are provided by postprocessing differential phase measurements. The positioning solution errors are neglected due to the fact that the solution errors level is at several centimeters, and the accuracy of normal gravity and free-air correction and the Eotvos correction is up to 0.1 mGal in such positioning accuracy [28].
3. The State Space Model Equations
There are two steps of Kalman filtering implemented in gravity anomaly solution: system propagation and measurement update. During the propagation step, gravimeter measurement equations are executed to predict the altitude, vertical velocity, and calibration parameters using gravimeter, horizontal accelerometer, and gyro measurements. Then, the altitude data from GNSS are used to update the innovation for correcting the predictions, which is called estimation in measurement update step. In this study, the propagation process is based on the development of system state equations, and the measurement update is based on that of altitude measurement equations.
3.1. Gravimeter Measurement Equations
The equation of gravimeter can be expressed as [25]where is the observation of gravimeter; is projection component of specific force working on the gravimeter sensitive axis; is gravimeter measurement error; is the gravimeter scale factor; is the onboard delay time.
Assume that the horizontal accelerators can be set as perpendicular to each other and one points to the east and the other points to the north. The readings of the horizontal accelerators are and , respectively. Their measurement equations can be expressed aswhere and are the actual horizontal accelerations in the east and north; and are the horizontal accelerometers errors.
Considering the compensations of error when the gravimeter was installed on the stable platform, as well as compensations for the flight deviation from horizontal, can be expressed aswhere and are the angles of the stable platform deviating from horizontal; and are the angular deviation from horizontal when gravimeter was installed on platform, and they are constant during flight but vary from flight to flight.
Therefore, the equation of gravimeter involving sensor error correction can be obtained from equation (5) by replacing true value of equation (7).where , , , and are regarded as gravimeter calibration parameters. In general, they are modeled as random constant and estimated by Kalman filtering in gravity anomaly solution.
According to equation (8), the equation of the specific force component working on the SM projection on the gravimeter sensitivity axis can be determined by the readings of gravimeter, the readings of horizontal accelerators, the angles of the stable platform deviating from the horizontal, and the estimated calibration parameters.
3.2. Altitude Measurement Equations
Actually, the altitude and vertical velocity from GNSS can be used and the performance is better in gravity anomaly solution when the altitude is used [22]. Therefore, the measurement equation is presented by using the altitude data from GNSS.
Assume that the altitude of antenna of GNSS receiver can be expressed aswhere is the measurement value of antenna altitude; is the true value of antenna altitude; is the measurement error. Because of the attitude changing during flight, the attitude correction is expressed as [7]where and are the pitching angle and rolling angle, respectively, which are measured during flight; , , and are the distances between the antenna and the SM of gravimeter in the carrier coordinates. Therefore, the altitude correction is
Then, the altitude of SM of gravimeter can be expressed aswhere is the gravimeter altitude provided by altitude correction; is the true value of altitude of gravimeter; is the measurement error, which is assumed as zero-mean white noise and its intensity does not equal zero.
3.3. State Space Model Development
Let us start from developing the system state equation. The following equation is obtained by substituting equation (4) from equation (8):where notations are used asand in equations (13) and (14), is the sum of the various vertical forces on gravimeter; is the sum of the various noises, and it is considered as system noise. Then, the differential equations can be obtained aswhere is the vertical velocity. According to the discussion about the calibration parameters (, , , and ), they are modeled as random constants, the intensity of which is zero; then
To obtain the gravity anomaly by Kalman filtering, the gravity anomaly considered as a stochastic process along the flight line can be modeled by a filter and expressed as [25]
Combining equations (16) and (17), they can be described by the following linear dynamic model. In the model, is state vector, is dynamics matrix, expresses the system noise distribution matrix, is system noise vector, and is an additional term.
The dimensions in equation (18) are specified as
The measurement model contains the functional relationship between the observations and the state vector. Generally, GNSS altitude and/or velocity data can be used as measurements in Kalman filtering, and the altitude will be used in this paper. Equation (12) is used to develop the measurement model, which is described aswhere is the measurement vector; is the measurement noise vector; is the measurement matrix, which is specified as
Moreover, it is evident that the desired is the component of state vector. Then, can be estimated by forward Kalman filtering, and the results are smoothed using backward RTS smoother. Note that the calibration parameters can also be estimated to correct the readings of gravimeter.
4. Data and Processing
4.1. Surveying Data
The surveying data used in the experiment was supplied by the airborne scalar gravimetry system GT-1A imported by China Aero Geophysical Survey and Remote Sensing Center for Natural Resources. GT-1A is the system with three-axis stabilized platform. Three accelerators are installed on the platform and they are perpendicular in pairs. The two horizontal accelerators are used for measuring the accelerations during flight. Schuler tuning platform is used as the stabilized platform, and the vertical orientation of the gravimeter (vertical accelerator) is maintained by the feedback from accelerators and gyros.
The survey was conducted over a sea-land transitional region located in Jiangsu province of southeast China. The survey unit was mounted inside Cessna 208 with the x-axis pointing to the front, y-axis pointing to the right-wing, and z-axis pointing down in directions and switched on more than 30 minutes in both precalibration phase and postcalibration phase at the airport. The distances between the SM of gravimeter and GNSS antenna were precisely measured before precalibration. There are two repeated flight survey lines (lines 1 and 2), the lengths of which are nearly 80 km, and survey lines 1 and 2 are both in the latitudinal direction. During the flight, the two lines flew at nearly a constant altitude of 400 m above the sea level, though the actual altitude fluctuated within 10 m, which is shown in Figure 1. The flight was conducted at an average ground speed of 60 m/s and at night for better meteorological conditions prevailed for reducing noises by air flowing.
4.2. Processing and Results
In this section, the airborne gravity anomaly solution based on the developed state space model using Kalman filtering and smoothing (KFS) is conducted. When the GNSS/INS postprocessing is implemented, the Kalman filtering formulas are used to estimate the states and to update the estimating covariance matrix using measurement matrix, measurement noise covariance matrix, filtering gain vector, measurement innovation vector, and predicted error covariance matrix in every measurement epoch. Then, the estimated states are used to predict the states of the next epoch. Finally, the estimated states of Kalman filtering are smoothed using RTS smoother. Note that the gravity anomaly solution obtained by GT-1A software GTGRAV was used as the solution reference.
Besides, the gravity anomaly solution is also conducted using FIR filter for comparison with that of KFS. According to the basic equation (1), the direct gravimetric measurements minus vertical kinematic accelerations derived by differential calculation from GNSS/INS as well as other deterministic effects, such as Eotvos and normal gravity correction mentioned in equations (2) and (3). The result is the vertical component of the gravity with the random errors and noises, and it is called original gravity anomaly in this section. The FIR processing is based on a frequency analysis of the original gravity anomaly. The PSD of the original gravity anomaly and reference standard on survey lines 1 and 2 are shown in Figures 2 and 3. It can be seen from the figures that the power of the gravity anomaly signal is mainly distributed in the low-frequency range. In the range from 0 to 0.004 Hz, the PSD of the original gravity anomaly is comparable with that of the reference. Then, the FIR filter is designed and implemented using Hamming function with 600 orders and 0.004 Hz cut-off frequency. It should be noted that the option of FIR properties is a trade-off between the solution accuracy and resolution [17, 18], and zero-phase algorithm is implemented in the FIR solution according to Zhang and Xu [20, 21].
The gravity anomaly solution based on the developed state space model using KFS (black) is shown in Figures 4 and 5, and the solution obtained by FIR (blue) is depicted in the figures. Note that, in the results of FIR, the parts affected by the boundary effect have been removed.
As can be seen from the figures, the gravity anomalies along line 1 and line 2 obtained by the developed model using KFS are consistent with the reference. Especially in the flat parts of reference, the gravity anomaly curve of KFS is almost coincident, which indicates that the result from the presented model using KFS is valid and effective.
Moreover, comparing the results obtained by FIR with those of KFS, the deviation between the FIR filtering result and the reference result was obvious. Especially at the initial part on line 1 and the end part on line 2, the fluctuations and the deviation with the reference are relatively severe. In contrast, the trend of KFS is more fit, and the noises are limited due to the process in time domain with calibration compensation considered in the solution.
The error statistics for the approaches based on KFS and FIR are given in Table 1, respectively. According to the table, the result of FIR showed less accurate solution than the approach based on KFS. The accuracies of the solution results provided by FIR are almost 50% (line 1) and 40% (line 2) worse than those of KFS, and the mean error of the solution results showed the same trend. Moreover, both mean error and RMSE using KFS are less than 1 mGal, which is generally considered as acceptable accuracy in geophysical prospecting.
5. Conclusions
An approach to determine airborne gravity anomaly in airborne scalar gravimetry is proposed in this paper. The approach is developed based on the airborne gravimeter equations for the state space modeling. The gravimeter equations are generated by several measurement errors, corrections, and calibrations. The developed state space model with the measurement equation using Kalman filtering and smoothing is conducted for gravity anomaly solution.
The practical surveying data were processed using proposed approach, and the results showed the accuracy of 0.79 mGal and 0.72 mGal along survey line 1 and line 2, respectively. Moreover, the approach was compared with the traditional one based processing in frequency domain. The accuracy of the solution obtained using the proposed approach is two times better than that of FIR. In conclusion, the proposed approach is effective for gravity anomaly solution in airborne scalar gravimetry, and the accuracy is acceptable in geophysical prospecting.
Data Availability
The measurement data used to support the findings of this study were supplied by China Aero Geophysical Survey and Remote Sensing Center for Natural Resources (AGRS) under license and so cannot be made freely available. Request for access to these data should be made by contacting Rui Li, 510231865@qq.com.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This research was supported by the Natural Science Foundation of Henan Province (Grant no. 202300410101) and the Initiation Project for Ph.D. of Henan Institute of Technology (Grant no. KQ1813).