Experimental study on improving the accuracy of marine gravimetry by combining moving-base gravimeters with GNSS antenna array

The gravity field is one of the Earth’s basic physical fields. The geoid can be calculated and the tectonic activity underground can be inversed by gravity anomaly. With the development of various ship-borne gravimeters and navigation technology, including the Global Navigation Satellite System (GNSS) and Strapdown Inertial Navigation System (SINS), the precision of marine gravimetry has been significantly improved (achieve or better than 1mGal). Errors arising from calculations of the correction term have become the main source of gravity measurement errors. At present, the traditional approach is to deploy a GNSS antenna, connect the GNSS antenna to the gravimeter, record the real-time position through data acquisition software, and then use this position to calculate the gravity correction item afterward. Two errors are inevitable. (1) The GNSS antenna position error is large. The pseudorange point positioning method is generally used to obtain real-time GNSS antenna positions, and the positioning accuracy is poor compared with that of precise point positioning. (2) The position coordinates of the gravimeter contain systematic errors related to the ship’s attitude. In this paper, a joint experiment including GNSS antenna arrays and ship-borne gravimeters was designed to evaluate the measurement accuracy via repeat lines on the same ship. The experimental results show the following: (1) attitude accuracies of 0.0299° for the yaw angle, 0.0361° for the pitch angle, and 0.1671° for the roll angle can be obtained at baseline lengths of 25 and 4 m. (2) The GNSS antenna array has an obvious role in determining the point acceleration in the low-frequency band (0–0.01 Hz) and the point position and velocity in the high-frequency band (0.01–1 Hz). (3) The vertical position eccentricity causes an absolute error of 1 mGal and a relative error of 10-1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${10}^{-1}$$\end{document} mGal in gravity measurements and can be corrected by the GNSS antenna array method. (4) Using a GNSS antenna array can obviously improve the measurement accuracy of an instrument with a precision equaling or exceeding 1 mGal, but cannot obviously improve that to an instrument with poor precision (2 mGal or below).


Introduction
The gravity field is one of the Earth's most important physical fields; and all of humankind lives within and is influenced by the gravity field. It is well known that the ocean covers 70% of the Earth's surface. With the application of satellite altimetry, the marine gravity field can be inverted efficiently. However, the data obtained through satellite altimetry are sometimes insufficient for establishing a highly accurate geoid model because of their low spatial resolution and low accuracy (Sandwell and Smith 1997;Ebbing et al. 2013;Zhu et al. 2019). Therefore, highprecision dynamic gravimetry is important and necessary .
Marine gravimetry was developed after terrestrial gravimetry. Vening Meinesz (1929) first measured the gravity field at sea in 1923 with a precision of 4-5 mGal; he used two pendulums with the same amplitude but opposite phases to correct the influence of horizontal acceleration. He also investigated the relationship between gravity anomalies and tectonic activity. Subsequently, many instruments appeared using pendulums (Hu et al. 2017). In the 1950s, Lacoste (1967) and Lacoste and Bowles (1982) developed the first gravimeter with a stable platform utilizing a zero-length spring as the sensing unit. Ouyang et al. (2013) conducted repeated line and crossing point tests on 5 sea-air gravimeters carried by the same Y-8 fixed-wing aircraft, and the results showed that the Russian GT-1A airborne gravimeter had the best comprehensive performance and accuracy. Cai et al. (2017) analyzed the error characteristics of the strapdown gravimeter SAG-ZW. Yuan et al. (2020) selective analyzed the error characteristics of different types of gravimeters with different principles. Yuan finds that GT-2M gets the best performance in the rough sea condition compares with the other five gravimeters (SAG-2M, CHZ-II, SGA-WZ, ZL11, and L&R). The precision of 5 gravimeters (except L&R) is less than 0.5 mGal. With the development of associated hardware and data processing methods, the precision of gravity measurements continues to increase (Ishihara et al. 2018;Hwang et al. 2014). After obtaining specific force observations, two fundamental steps must be performed to obtain accurate gravity values. The first step is various error corrections, such as Eötvös corrections, vertical acceleration corrections, and free air corrections. The accuracy of the correction directly affects the gravity observation. Accordingly, this paper will focus on methods for applying these corrections. The second step is low-pass filtering, which is important because gravity data are composed of low-frequency signals, and thus, the noise amplitude is generally several times that of the true gravity signal. A finiteimpulse response (FIR) filter, a common low-pass filter, is often used to restore the gravity signal (Liang 2012;Ouyang 2013;Saha et al. 2012;Hehl 1992). However, this paper will not focus on the low-pass filtering of gravity data.
Eötvös derived the formula for the Eötvös correction in 1919, after which Harlan (1968) provided another rigorous formula in which the Global Positioning System (GPS)-derived east-west or north-south velocity can be used directly. The free air correction and normal gravity correction are related to the normal elliptic parameters (physical geodesy), for which there are strict formulas. Zheleznyal and Mikhailov (2012) found that the accuracy of the GNSS-derived height is insufficient when employed for the free air correction. Thus, they deployed a GNSS receiver on a fixed base near the route of a gravity survey and applied the 1996 Earth Gravitational Model (EGM-96) to transform ellipsoidal heights into orthometric heights; accordingly, a higher precision was obtained. However, this method is not practical when the region to be surveyed is far from land.
Most of the above-mentioned corrections are related to or derived from position data. At present, GNSS systems are the most commonly used approach for measuring position data, and the coordinates of the antenna can be calculated conveniently. However, the coordinates of instruments in the laboratory cannot be obtained directly, and taking the antenna position as the instrument position can produce errors in the correction calculations (Huang 1995). Sun et al. (2003) analyzed the eccentricity correction method in the airborne gravity survey with attitude sensor and concluded that the attitude angle accuracy should be better than 0.3° to obtain sufficient eccentricity correction accuracy. The eccentricity correction method of gravity measurement has been mentioned all the time, but few people have analyzed the specific influence of eccentricity correction on different gravimeters and different corrections in detail. Alternatively, the position of a gravimeter installed indoors can be calculated if we can obtain the attitude of the vessel, the position of the GNSS antenna and the baseline between the antenna and instrument; this process is called an eccentricity correction. The attitude of the vessel can be measured by a SINS or determined by several GNSS antennas (no fewer than 3) mounted on the surveying vessel (Lavrov et al. 2017;Varble et al. 2020). This paper will compare the accuracy and error corrections of two different methods to compute the position: using the positions of GNSS antennas and using the position of a gravimeter after an eccentricity correction.

Methodology for measuring and correcting the gravity field
The theory employed for surveying the gravity field is Newton's second law, and the equation is: where δ g is the gravity anomaly, which is located on the geoid; f U is the specific force in the vertical direction output by the gravimeter; v U is the vertical acceleration of the gravimeter at the moment when obtaining f U ; δ a E is the Eötvös correction, which is related to the velocity of (1) the gravimeter; and γ is the normal gravity value, which is located on the geoid.
However, in relative gravity measurements, the following formula is more commonly used: where g b is the gravity anomaly at the gravity benchmark; f 0 U is the mean specific force in the vertical direction at the gravity benchmark; δ a H is the horizontal correction; δ a F is the free air correction, which is related to the orthometric heights; δ a K is the drift correction (Jekeli 2012;Wang et al. 2018;Dodson et al. 1997).
The formula of the Eötvös correction is expressed as follows (Harlan 1968;Sun 2004): where v ( m/s ) is the velocity of the gravimeter; v E is the velocity on the east-west component; ω is the Earth's angular velocity of rotation; f is the first oblateness of the ellipsoid; a is the semimajor axis of the ellipsoid.
The velocity (v) and acceleration ( v ) are calculated by taking the derivative of the position using the following equations: In the corrections mentioned above, v U , δ a H , δ a E , δ a F , δ a K and γ are related to the instrument position. There is a large body of literature that studies their formulas, which will not be introduced here. Their differences before and after the correction are compared in "Experimental design and data processing" section.

Methodology for calculating the position of the gravimeter
Generally, the coordinates of GNSS antennas are geographic coordinates comprising the latitude ( ϕ ), longitude ( ) and ellipsoidal height (h). These coordinates are converted into a spatial Cartesian coordinate system under an Earth-fixed reference system using the following equations: (2) where (x, y, z) are the coordinates of the Earth-fixed reference system; N is the curvature radius of the unitary circle of the Earth at the coordinate point; e is the Earth's first eccentricity: ,n sin β sin γ −y 13,n sin β cos γ +z 13,n cos β x 13,n cos γ −y 13,n sin γ where α is the roll angle, β is the pitch angle, γ is the yaw angle and x ij , y ij , z ij n is the baseline vector from antenna i to antenna j in the navigation coordinate system. Subscript e, n, b represents the Earth-fixed coordinate system, the navigation coordinate system, and the body coordinate system, respectively. The detailed derivation process can be found in the literature (Cai et al. 2018;Liu and Ou 2003;Aleshechkin 2011;Zhang and Hao 2017). After obtaining the attitude angle, the position of a gravimeter installed indoors can be easily calculated using the following formula: where x, y, z e are the coordinates in the Earth-fixed coordinate system, the subscript g indicates the position of the gravimeter, the subscript a indicates the position of the antenna; x ag , y ag , z ag T b is the baseline vector from the antenna to the gravimeter; C b n T is the rotation matrix from the body coordinate system to the navigation coordinate system, and C n e T is the rotation matrix from the navigation coordinate system to the Earth-fixed reference system.

Error of the attitude angle
First, we take the derivative with respect to the (x, y, z) of Eq. (6): where d u means differentiate with respect to variable u, in the actual situation, it also means the error of the variable; s 12 (m) is the distance from antenna 1 to antenna 2, S is the distance from antenna 3 to baseline s 12 equal to x 13,n cosγ − y 13,n sinγ , and Q is an auxiliary quantity that is quite small, where Q = −x 13,n sinβsinγ −y 13,n sinβcosγ +z 13,n cosβ x 13,n cosγ −y 13,n sinγ . In the real world, Q, d α , d β , d β andz 12,n ≪ x 12,n ≈ y 12,n < s . ( ≈ means that x 12,n is of the same order of magnitude as y 12,n ). Thus, Eq. (10) is simplified to the following: (10) x 2 12,n +y 2 12,n � d z 12 ,n −x 12,n z 12,n d x 12,n −y 12,n z 12,n d y 12,n � x 2 12,n +y 2 12,n +z 2 12,n �� x 2 12,n +y 2 12,n x 12,n +d 2 y 12,n s 12 .
Hence, the accuracies of the two horizontal angles α and β are related to the vertical positional accuracy of the antenna and the distance of the baseline. Similarly, the yaw accuracy is related to the horizontal positional accuracy of the antenna and the equivalent distance of the baseline. In reality, the distance of the baseline on the ship can range from 10 to 80 m or more. Using the GNSS differential positioning method, we can obtain a horizontal positional accuracy on the millimeter or centimeter scale and a vertical positional accuracy on the centimeter scale (Dong et al. 2020). Then, we can forecast accuracies of approximately 0.01° for the yaw angle, 0.1° for the pitch angle and 0.2° for the roll angle.

Error of the eccentricity correction
For convenience, an accuracy analysis is carried out in the navigation coordinate system. Equation (7) can be rewritten as follows: The error of the eccentricity correction clearly consists of two parts: the GNSS antenna position error and the attitude angle error. The emphasis here is on the latter. We take the derivative of Eq. (12) with respect to the attitude angle: ig. 1 Errors of the eccentricity correction varying with the attitude angle caused by the attitude angle error. a, d and g show the east-west position errors of the eccentricity correction; b, e, and g show the north-south position errors of the eccentricity correction; and c, f, and i show the elevation position errors of the eccentricity correction. The baseline is taken as (10, 20, 10), and the attitude angle error is set to 0.1° Each error is related to nine different quantities: three attitude angles, their errors, and three baseline coordinates. Thus, to simplify the above, we employ the baseline (10 m, 20 m, 10 m) and analyze three different yaw angles (0°, 45°, 90°) to determine the error distribution characteristics.
As shown in Fig. 1, where the error of each of the three attitude angles is 0.1° and the baseline is (10 m, 20 m, 10 m), the direction of the maximum horizontal position error rotates with the yaw angle, while the vertical position error is independent of the yaw angle. The maximum errors are listed in Table 1.

Gravity and GNSS data
The data employed in this study are from the 2018 Xiangyanghong No. 06 gravity comparison survey in the South China Sea. Three GNSS antennas were mounted on the exterior of the ship, and three gravimeters, namely, GT-2M (No. 39), SAG-2M and L&R Air Sea System II (No. S129), were installed on the ship. Their coordinates are shown in Table 2.
It should be noted that the SAG-2M gravimeter is a strapdown gravimeter with a triaxial accelerometer and a triaxial gyroscope that can be used as a SINS. Its corresponding data processing software also uses inertia's integrated navigation algorithm to calculate gravity. SAG-2M gravimeter adopts a high-precision fiber optic gyroscope and a high-precision quartz accelerometer, and the attitude resolution accuracy is better than 0.002° through inertial /GNSS combined navigation calculation (Table 2). Figure 2 demonstrates that the GNSS antennas are positioned higher than the gravimeters. Given this allocation of positions, the value of S in Eq. (11) is not long enough to obtain a higher accuracy of roll angle precision. L&R is installed at a lower height than the other gravimeters. The receiver corresponding to Antenna 1 is a NovAtel PwrPak7. Antennas 2 and 3 correspond to a Trimble R9 receiver.
There are four repeated ship lines: a north-south line named S1, an east-west line named H1, and two northwest-southeast lines named X1 and X2. The information of these repeated ship lines is provided in Table 3. The location on the map is shown in Fig. 3. In Fig. 4, thirty-one crossing points were selected for subsequent accuracy assessment.
In the following text, X1-2 denotes the second repeated line of X1.
The roll angle reflects the sea conditions for each route. In Fig. 5, the conditions of the X1 lines are the worst; line S1-1 is also bad. In contrast, the sea conditions for all X2 lines are good. All H1 lines are not bad (Fig. 6, Table 4).
In Table 4 and Fig. 6, the average velocity and velocity change of each survey line are shown.

Data processing flow
This test involves the processing of three kinds of data: gravity data, GNSS data and SINS data. Gravimeters provide specific force observations for gravity measurements. SINS instruments can provide position, velocity, acceleration, and attitude information for gravity measurement corrections, but the accuracy of pure inertial navigation data drifts with time, so external reference information must be provided. Data based on GNSS antenna arrays can also provide auxiliary information for gravity measurement corrections, and the accuracy does not drift over time. Differential positioning can obtain the relative position with a high precision (millimeter-level accuracy for the horizontal direction and centimeter-level accuracy for the vertical direction) (Dong et al. 2020) and thus can be used to determine the attitude. Using precise point positioning (PPP) can yield decimeter-level and even centimeter-level absolute dynamic position accuracies (Li et al. 2013(Li et al. , 2010Heroux and Kouba 2001); hence, PPP can provide the location reference for SINS. Alternatively, the difference between epochs for computing the speed, Doppler observations and carrier phase observations can provide higher-precision velocity and acceleration information (Xu et al. 2015). In this paper, the Moving-Base mode of rtklib (http:// gpspp. sakura. ne. jp/ rtklib/ rtklib. htm) is used to process the multi-antenna GNSS data to get the baseline. The internal coincidence accuracy of the processing results can reach the level of millimeters. The attitude angle is calculated by the baselines and the attitudes calculated by the two methods (SINS and GNSS antenna arrays) are compared. The absolute position is calculated by the PPP method of MG-APP (Xiao et al. 2020), which was able to achieve accuracy on the centimeter to millimeter scale. Finally, a unified 200 s FIR filter was applied to the corrected gravity data (Fig. 7).

Comparison of the multi-antenna attitude with the inertial navigation system attitude
The GNSS antenna array attitude angles and SINS attitude angles are compared, and error statistics are obtained. In Table 5, the precision of the roll angle is worse than the precision of the yaw and pitch angles because the value of S in Eq. (11) is not long enough, as mentioned in "Gravity and GNSS data" section (Figs. 8,9,10). Attitude determination requires simultaneous observation of at least three GNSS measurements, and the absence of one will cause the attitude to fail to be calculated. During the measurement of the X1-2 survey line, the receiver corresponding to Antenna 2 was not turned on, so the multi-antenna attitude determination could not be carried out. The attitude angle of SINS was used in the subsequent processing of the X1-2 survey line.

Fourier time-frequency analysis of the specific force
In specific force observations, which are inevitably accompanied by high-frequency noise, the amplitude of noise is often hundreds of times that of the gravity signal, and the noise signal occupies a wide frequency band. In practical processing, a reasonable filter should be designed according to the frequency band of the gravity signal to extract the weak gravity signal. Accordingly, this paper compares the post-processed GNSS antenna array data with the GNSS observation data to recalculate whether the gravity measurement corrections are improved in the corresponding frequency band of gravity signals.
Spectrum analysis of the specific force observations from the entire voyage (Fig. 11, the small windows are enlarged views of the low-frequency band. The small windows are enlarged views of the low-frequency band) reveals that the effective range of the gravity signal is below 0.005 Hz, and there is a wide frequency band without a signal from 0.005 to 0.01 Hz. Therefore, the cutoff frequency of the low-pass filter can be selected within this frequency band, and thus, the focus of the gravity correction calculation should be within 0.005 Hz (or even 0.003 Hz).

Speed comparison before and after correction
The horizontal velocity error is the main influencing factor in the calculation of the Eötvös correction. In gravimetry, the low-frequency domain should be of interest. As Fig. 12 shows, the horizontal velocity is greatly improved in the time domain; in the frequency domain, waves cause the ship to sway at a high frequency (0.08-0.22 Hz), and the gravimeters are less affected than the GNSS antennas because the former are installed at a lower height. While the speeds of these instruments are very different in the high-frequency band, they do not differ much in the lowfrequency band. Therefore, the Eötvös correction error caused by velocity eccentricities will be small.

Eötvös correction
As mentioned in "Speed comparison before and after correction" section, the Eötvös corrections are scarcely different when calculated by the velocities before and after in the low-frequency band, but there is a large difference at high frequencies (Fig. 13).

Horizontal acceleration correction
For gravity measurements, specific force observations in the vertical direction must be acquired. Platform gravimeters need a stable platform to ensure the sensitivity of Fig. 6 The speed of the ship: a for H1, b for S1, c for X1, and d for X2. The speed is somewhere between 9 and 11 knots their internal unit pointing in the vertical direction, while strapdown gravimeters need to measure the attitude of the instrument in real time to calculate the projection of the triaxial acceleration in the vertical direction. However, certain errors are always encountered in the control of a stabilized platform, and strapdown gravimeter attitude calculations also contain errors. These errors will detrimentally impact the vertical acceleration measurement accuracy and result in an unnecessary horizontal acceleration component; therefore, a horizontal acceleration correction is needed. Moreover, the correction intensity rises as the platform angle and horizontal acceleration increase. In Fig. 14   acceleration changes obviously after the eccentricity correction and behaves differently at the lower frequencies where the gravity signal is located.

Vertical acceleration correction
In marine gravimetry, the vertical acceleration is generally considered to have little influence on the gravity signal because of the specific measurement environment. However, the vertical acceleration is a necessary correction in airborne gravimetry (Schwarz and Li 1997;Krasnov et al. 2014). In this paper, the vertical acceleration signal of a marine gravity survey is analyzed, and the necessity of its correction and the influence of the eccentricity correction on the vertical acceleration are analyzed. Figure 15 reveals that the change in the vertical acceleration after the eccentricity correction is not as great as the change in the horizontal acceleration in the time domain, and the same is true in the frequency domain. However, the horizontal and vertical accelerations both exhibit considerable changes in the low-frequency band.
The main change is in the vertical acceleration calculated by the antenna position, as it has a larger amplitude in the low-frequency band than the vertical acceleration corrected by the eccentricity. However, the vertical acceleration signal in the low-frequency band should be smaller in ocean gravity measurements, so the vertical acceleration should not be corrected during postprocessing. Furthermore, because the vertical Fig. 8 Difference in the yaw angle: a for S1-1 and b for X1-1 Fig. 9 Difference in the pitch angle: a for S1-2 and b for X1-1 acceleration after the eccentricity correction obtained by postprocessing is different in the low-frequency band from the vertical acceleration before the correction (at the same antenna position), in airborne gravimetry, it is necessary to carry out the eccentricity correction before the vertical acceleration correction.

Free air correction
The free air correction adjusts the normal ellipsoid vertical gravity gradient. The free air correction values are approximately 1 mGal for 3 m. As shown in Fig. 16, the eccentricity causes an elevation deviation of approximately 4.3 m, which will result in an absolute error of

Precision statistics of repeat lines and crossover points
In this section, we calculate the maximum (Max.), minimum (Min.), average (Mean.), and root mean square  (RMS.) of the measurement error for each repeated line. The statistical results are shown in the tables below. The notation H1-1-2 indicates the difference between the first and second lines of H1.
As shown in Tables 6 and 7, the accuracy of gravimeters GT and SAG is improved because of the low measurement noise of the instruments. However, in Table 8, the improvement in the measurement accuracy is not

Conclusion
Concerning the proposed postprocessing method for gravity calculations using the GNSS antenna array attitude positioning method, the following conclusions are obtained: 1. The determination of an instrument's position and attitude using a GNSS antenna array is an efficient, low-cost and high-precision method. In marine gravimetry, the advantages of the ship's large size and long baseline can be utilized to obtain a more accurate attitude, and the marine measurement environment is free of shielding effects; hence, the measurement stability can be guaranteed. When the effective baseline length is 25 m and 4 m, the attitude angle accuracies obtained in this test are 0.0299° for the yaw angle, 0.0361° for the pitch angle and 0.1671° for the roll angle. 2. The eccentricity correction has a great influence on the overall calculation of the horizontal velocity, but it has little influence on the low-frequency part corresponding to the gravity signal and thus has little influence on the calculation of the Eötvös correction in gravity measurements.  3. The horizontal acceleration after implementing the eccentricity correction exhibits obvious changes in both the low-frequency and the high-frequency parts, especially under poor sea conditions. Therefore, the eccentricity correction should be carried out first when the horizontal acceleration correction is considered in gravity measurement. The vertical acceleration does not change obviously in the highfrequency band, but displays considerable variation in the low-frequency band. Hence, the eccentricity correction should be carried out first when the vertical acceleration is corrected in airborne gravity data, but can be ignored in marine gravity surveying. 4. The absolute error of the free air correction caused by the position eccentricity is large (at the mGal level) and thus cannot be ignored in marine and airborne gravity measurements. If the correction is carried out with a constant value, an error on the scale of 10 −1 mGal will remain, which should be considered in gravity measurements. Measured data show that the relative accuracy of gravity data can be steadily improved by 0.02 mGal after applying the eccentricity correction in this experiment. 5. The method of using a GNSS antenna array can obviously improve the measurement accuracy of instruments with a precision equal to or exceeding 1 mGal, but this approach is not effective for instruments with poor precision (2 mGal or below).