Geoid of South East Sulawesi from airborne gravity using Hotine approach

In the past, geoid was computed from gravity anomaly data using Stokes or Molodensky approaches. Obtaining gravity anomaly data is difficult because it needs some reductions of gravity from surface of the earth to the geoid using orthometric height from spirit level measurement. In the modern era, gravity anomaly data may be replaced by gravity disturbance data. It only required gravity and GNSS (Global Navigation Satellite System) measurement. This research aimed to determine geoid using Hotine’s approach. Disturbance data were generated from archived free air anomaly of airborne gravimetry in Sulawesi area. South East Sulawesi province was selected as a case study area. In this study, gravity observation was calculated at an altitude of 4000 m above the reference ellipsoid. Gravity estimation at the same height aims to increase the precision of the downward continuation process to the geoid. Hotine integral is calculated above the geoid, so that the gravity disturbance data is downwarded to the geoid. The geoid undulation is graded from north to south. Geoid from airborne gravity around Pegunugan Mekongga in the northern part of Southeast Sulawesi Province has the largest geoid undulation which reaches 63 m, while the geoid in the southern part of Buton Island reaches 52 m. Geoid validation of airborne gravity at 13 test points produces a standard deviation of ± 0.050 m. The standard deviation is much smaller than the results of geoid testing from airborne data in North Sulawesi, Central Sulawesi and Southeast Sulawesi. This fact indicates that the Hotine approach has the potential to produce a precise geoid if used in geoid-based airborne gravity calculations.


Introduction
Geoid is an equipotential gravity field which coincides with mean sea level. Visually, on sea region, geoid can be imagined as the average sea level in a undisturbed state. On land region, geoids are more imaginary or abstract. The position of the geoid below the earth's surface is represented by the height value or geoid undulation above the reference ellipsoid. Practically, geoid is applied to convert the geodetic height to the orthometric height. Geodetic height of GNSS refers to the reference ellipsoid while orthometric height of spirit leveling refers to the mean sea level.
Globally, geoid can be calculated from the EGM2008 harmonic coefficient. In Indonesia, geoids are graded from -40 m in Banda Aceh Province to +80 m in Papua. Since the computatior of EGM2008 coefficients in Indonesia is mostly based on data from satellites, the geoid accuracy of EGM2008 is still in the decimeter fraction. Geoid accuracy can be improved by combining terrestrial, airborne and shipborne gravity data.
Due to the large geographical area and the problem of accessibility, terrestrial gravity measurements are still inadequate to make precise geoids. Marine gravity measurements have also not been carried out because of Indonesia's vast ocean area. The serious effort to make a precision geoid began with airborne gravity measurements in Indonesia, which were first carried out in 2008 to 2011. The measurements were carried out by BIG in collaboration with DTU and NGA covering the regions of Kalimantan, Sulawesi and Papua.
There are not many publications about geoid in eastern Indonesia using airborne gravity data. Pahlevi et al (2015) calculated geoids using the Molodensky approach in the Sulawesi region. Geoid validation of airborne gravity with geoid geometries on levelling lines in three provinces. The quality of geoid can be seen from the root mean squares error (RMSE) value and standard deviation (SD). The RMSE value represents the accuracy of the geoid, while SD represents the precision of geoid. Practically, RMSE can be considered as offset or datum shift which shows the global relationship IOP Conf. Series: Earth and Environmental Science 731 (2021) 012014 IOP Publishing doi:10.1088/1755-1315/731/1/012014 2 between the reference ellipsoid and the geoid. Meanwhile sd represents the consistency of the geoid gravimetric deviation from the geometric geoid. In North Sulawesi, RMSE and SD are ± 1,061 m and ± 0.352 m. In Central Sulawesi, RMSE and SD are ± 0.521 m and ± 0.167 m. In South Sulawesi, RMSE and SD are ± 0.686 m and ± 0.168 m. Gravity testing was not carried out in the Southeast Sulawesi province on the grounds that the geoid validation path in the province had not yet been completed.
Sulawesi gravimetric Geoid calculation with free anomaly air data from airborne using Molodensky approach. This approach requires elevation when observing in orthometric height systems. The orthometric height for the reduction of free air is obtained from the reduction in geodetic height from GPS observations on the plane to the geoid undulation of EGM2008. The procedure actually has the potential to reduce altitude accuracy because there will be geoid error contributions. To avoid the height conversion process, the geoid can be calculated using the Hotine approach using gravity disturbance data. With the availability of geoid geometric points in the province of Southeast Sulawesi, it is necessary to try to calculate the geoid in the region. This study aims to calculate the geoid in the Southeast Sulawesi region with the Hotine approach

Data and equipment
In general, airborne gravity on Sulawesi Island was recorded in the north to south, east to west, and northwest to southeast directions, as shown in Figure 1. For the Southeast Sulawesi region, airborne gravity measurements were carried out on 11 main routes in the northwest to southeast and they were controlled by 1 cross over from the one perpendicular to the main line, as in Figure 2. The distance between the flight paths was about 18 km or the equivalent of 5 arc minutes. The flight path intervals were similar to the spatial resolution of EGM2008 with n = 2160. Geoid calculation using Hotine approach requires gravity disturbance data, Global Geopotential Model from EGM2008, and Digital Elevation Model from SRTM30 and GEBCO. Gravity disturbance data obtained from the conversion of Free Air Anomaly. BIG and DTU stored airborne gravity data in the form of free air anomalies that were archived with geodetic coordinates and data acquisition time on the platform. The Geoid was calculated in Gravsoft software of DTU by modifying the Stokes function script into the Hotine function.

Quality Control of airborne gravity data
Measurements above Southeast Sulawesi were carried out on 26 September 2008 to 6 October 2008. In case of calculation using raw data, the measurement precision can be seen from the standard deviation of measurements and tilt at the time of measurement. For this research, lack of standard deviation and tilt data were overcome by checking the consistency of main line and cross over data. Difference of more than 60 mgal were found in the cross over lines on 6 October 2008 at the border of Southeast Sulawesi and Central Sulawesi. The cross over data was discarded to avoid calculation errors. The cross over measured 29 September 2008 found differences in the main lines and cross over data that were graded from 2 to 10 mgal. Milligals are unit of gravity equivalent to 10-5 m.s -2 .
The speed of aircraft at measurement must be stable to avoid contributing aircraft acceleration to gravity data records. They varied from 224 to 326 km per hour (kph) which were calculated from aircraft position data and measurement time. The distance between measurement points per second ranged from 68 m to 91 m. The changing speed difference produced an acceleration of 0.287 to 0.323 m.s -2 or equivalent to 287 to 323 mgal. In this study, measurement points whose accelerating greater than 0.020 m.s -2 were removed, so that the amount of observation data was reduced from 63,619 to 9,177 points. The distance between the measurement points ranged from 71 m to 4405 m.

Computation of geoid
Geoid calculations with Stokes and Molodensky Approaches apply gravity anomaly. Gravity anomaly is the difference between gravity observation and normal gravity as shown by equation (1) as follows: 'g gobs J 0 (1) Gravity normal is computed from geodetic latitude by equation (2) as shown below: where a is semi-major axis, b is semi-minor axis, γ e is normal gravity at equator, γ p is normal gravity at pole, φ is geodetic latitude at observation point.
Geoid calculations require anomaly data between gravity in geoid and gravity in the reference ellipsoid, so the gravity value of the observation at the measurement point must be changed to the gravity observed on the geoid by reducing free air with equation (3) as follows: ( 3) where f is flattening of the ellipsoid, m is the comparison of centrifugal force at equator and gravity at the equator, h is the geodetic height above reference ellipsoid. Combination of equation (1) and (3) (2) gives free air anomaly formula as shown by equation (4) below: 'gFA gobs J 0 FA (4) In airborne measurements, gravity observations are measured on a platform, whereas normal gravity in ellipsoids is calculated from the geodetic latitude of GNSS measurements on the plane. Free air anomaly is obtained from anomaly gravity reduction by free air correction using orthometric altitude data from sea level. BIG and DTU estimate the orthometric height of the geodetic height reduced by geoid undulation from EGM2008. In this study, gravity observation on the platform was reconstructed from FAA values and normal gravity values as shown by equation (5) as follow: Normal gravity on the platform is estimated using geodetic latitude and height data with the following equation (6) Gravity disturbance on the platform are then calculated using combination of actual gravity and normal gravity on the same level as described by equation (7) below: The value is then amplified with the downward continuation procedure up to the geoid.
This research applied R-C-R (Remove-Compute-Restore) method [9]. The first step of R-C-R methods is reducing the long wave and short wave components in the airborne gravity data. The long wave component is represented by EGM2008 n=360 meanwhile short wave component was not reduced explicitly. Residual gravity disturbance (δgRES) were only obtained by reducing airborne gravity disturbance data (δgP) with gravity disturbance from EGM2008 (δgEGM2008) without reducing the topographic effect (δgtopo), as shown below: Residual gravity disturbances (δgRES) at measurement points are then interpolated to obtain continuous gridded residual gravity disturbance data (δgGRID) over calculation area. the COMPUTE part of R-C-R technique is to apply gridded residual gravity disturbance for computing residual geoid (NRES) as follows: where R is radii of the earth and H is Hotine's function. Hotine's function can be computed form spherical distance between evaluation point and gridded data as shown by equation (10) below: (10)

Result and discussion
Reconstruction of the free air anomaly data with geodetic position measurement data produces gravity data on the platform and gravity disturbance on the geoid. In this study, gravity observation was calculated at an altitude of 4000 m above the reference ellipsoid. Gravity estimation at the same height aims to increase the precision of the downward continuation process to the geoid. Airborne gravity measurements are only carried out along the land in the northwest to southeast direction. The oblique flyway configuration results in the measurement data not forming a rectangle as expected in geoid calculations. The results of gravity interpolation in ocean areas that are not measured by airborne gravimetry are of low quality, as shown in Figure 3. On flights over the Mekongga Mountains, gravity observations reach the highest values up to 976970 mgal. The combination of gravity observations and normal gravity at 4000 m elevation produces gravity disturbance at 4000 m elevation. In this study, the Hotine integral is calculated above the geoid, so that the gravity disturbance data is downward down to the geoid. Downward continuation precision is improved by filling in the unmeasured area with secondary data in the form of gravity disturbance from EGM2008 n = 2190 upward to an elevation of 4000 meters above the reference ellipsoid. EGM2008 data is only used in areas with a radius of 20 km from airborne gravity data. The results of the downward continuation to the geoid can be seen in Figure 4. The largest gravity disturbance value was in the geoid in the Mekongga Mountains reaching 175 mgal.
Referring to the results of cross-over checking, the precision of airborne gravity data rangedbetween 2 mgal to 10 mgal. The distance between the measurement lanes was 18 km or equivalent to the resolution of EGM2008 n = 2160. This condition was suspected to result in the geoid accuracy of the airborne being no better than the EGM2008 data. Checking the measurement line proved that the airborne gravity data resembled EGM2008 data, as shown by  The geoid from EGM2008 n = 2190 and the geoid calculated from the airborne data had almost the same configuration, as in Figure 6 and Figure 7. The geoid undulations were graded from north to south. Geoid from airborne gravity around Pegunugan Mekongga in the northern part of Southeast Sulawesi Province had the largest geoid undulation which reached 63 m, while the geoid in the southern part of Buton Island reached 52 m.
For geoid validation, BIG had measured the geodetic height at 14 Vertical Network points around Kendari City in Southeast Sulawesi. The geoid undulation of the airborne data at the test point ranged from 57.970 m to 59.540 m, while the difference between gravimetric and geometric geoids ranged from -0.312 m to -0.718 m, as shown in Table 1. The RMSE and the Standard deviation of the geoid from airborne were greater than the geoid from EGM2008. GNSS. If the outlier is removed, the geoid deviation standard from airborne gravity decreases to 0.050 m, while the geoid from EGM2008 n = 2190 becomes 0.088 m.