Verification of the Jakarta Geoid Model from the Gravity Data of 2.5 km Resolution with Gravimetric Geoid

To use the Global Navigation Satellite System (GNSS) correctly, the height information should be transformed into orthometric height by subtracting geoid undulation from it. This orthometric height is commonly used for practical purposes. In 2015 geoid of Jakarta has been produced, and it has an accuracy of 0.076 m. In the year 2019, airborne gravimetry has been done for the entire Java Island. The area of DKI Province cannot be measured because there is inhibition from Airnav. For this reason, terrestrial gravimetric measurements are carried out in this region by adding points outside the previously measured area. To compute the geoid in the Jakarta region is needed the Global Geopotential Model (GGM). In this paper, the GMM used is gif48. The “remove and restore” method will be used in calculating the geoid in this Jakarta region. Besides that in this geoid calculation also uses Stokes kernel and FFT to speed up the calculation. The verification of the resulting geoid is carried with 11 points in DKI Jakarta Province. This verification produces a standard deviation of 0.116 m and a root mean square of 0.411 m.


Introduction
Based on Regulation of The Head of Indonesia Geospatial Information Agency No. 15 2013, the Indonesian Geospatial Reference System 2013 (SRGI2013) has been appointed as the only mapping reference system in Indonesia. In SRGI2013, the vertical datum of any Indonesian map, both base, and thematic maps, is geoid. The base map is a product of Basic Geospatial Information (IGD) dan the thematic map is a product of Thematic Geospatial Information (IGT). IGD and IGT are part of the geospatial information Agency. The geoid used in mapping should be modelled or computed based on terrestrial gravity measurements tied to the national gravity control network (JKG). The geoid is formed due to variations in the earth's gravity. An adequate number of gravity measurements with high precission data are needed to obtain an ideal geoid result [1].
In making a large-scale topographic map, the geoid datum should be precise or have a submeter precision number [2] [1]. One good example is the preprocessing step of LiDAR point cloud data which has elevation values measured by the Global Positioning System (GPS). The GPS (or GNSS in a broader context) elevation values are measured relative to the mathematical ellipsoid reference and then called ellipsoidal height (h) [3]. To give these values more scientifically and practically, they should be corrected by geoid undulation (N). These corrected ellipsoidal heights are known as orthometric heights (H). Relation among ellipsoidal heights, geoid undulation, and orthometric height was defined in Equation 1 and as illustrated in Figure 1.  [4]. This geoid determination method is often called geometric geoid or GNSS-Levelling. It is difficult to obtain reliable results to compute the high accuracy geoid from a regional gravimetric with sparse data, which can meet today's geodetic community [5]. Each point on the earth's surface can be expressed as an orthometric height (H) or a geodetic height (h). For mapping and engineering practice, orthometric height is used as a height reference [1]. However, measuring geodetic height is more accessible in the real world than measuring orthometric height by observing from a satellite orbiting above the reference ellipsoid. Therefore, geodetic heights in mapping, geoid undulation values are required, and specific calculations to obtain orthometric height. From Equation 1, it is clear that one can compute geoid undulation (shorten as undulation) by subtracting orthometric height from ellipsoidal height [4]. If these undulation values are plentiful and evenly distributed, the national scale geoid model could be calculated. This geoid determination method is also known as geometric geoid or GNSS-Levelling.
Besides the aforementioned geometric geoid determination, the geoid can also be modelled using the gravimetric method. In this method, the geoid was computed based on terrestrial gravity data and the global geopotential model (GGM). One of these GGM is the GIF 48 made by Ries et al. in the year 2011, which is based on GRACE and GOCE gravimetric satellite [6]. Based on previous research by Ramdani et al., GIF 48 is the most optimum GGM for the Jakarta region [7]. It thus will be used as the long-wavelength component of our regional geoid model. The GIF 48 is made by Ries et al. in the year 2011, which is based on GRACE and GOCE gravimetric satellite [8]. In 2015 The geoid of Jakarta had been calculated from terrestrial gravity data measured back in 2012. The data used in this study has a 1.5 km spatial resolution, with a 0.076 m standard deviation and 0.626 m RMS.
[9]. The precision value was computed by comparing the gravimetric geoid and the geometric geoid.
In 2019, Geospatial Information Agency conducted an extensive survey of the previous measurement area (Jakarta province). This extension covered area in coordinate 106˚-107.2˚ E and 6˚-6.5˚ S with a spatial resolution of 5 km for each measurement point. Further augmentation was done by the Research Division, BIG, with a spatial resolution of 2.5 km (see Figure 2). Those datasets were used as the input for Greater Jakarta Regional Geoid Model.  Figure 2. The distribution of 2019 terrestrial gravity measurement points.

Methods
The basic equation to compute any geoid model is as written in Equation 2 [4].
The GGM-derived long-wavelength geoid component has an immense contribution to the whole geoid model, where φ dan λ are spherical coordinates; GM is gravitational constant (including earth atmospheric mass); a is the semi-major axis of reference ellipsoid; r is the distance from the computed point to the center of the earth; γ is normal gravity; Nmax is the maximum number of m degree and n order; C and S are harmonic coefficients, and P is the Legendre polynomial.
Meanwhile, to compute the local geoid model, Brun formulae will be applied (see  with is free-air gravity anomaly defined as the difference between gravity value on geoid surface (g) and normal gravity () as written in Equation 6.
The free-air gravity anomaly can be computed from terrestrial gravity measurement data reduced to geoid surface and corresponding normal gravity. This method is called boundary condition solution with Robin's problem. In this method, the geoid is the boundary itself; thus can be computed by using Stokes' integrals as written in Equation 7 [11]. with is the angular distance on the earth's surface. Gravity values on the earth's surface should be reduced to geoid surface except in areas of the ocean. There are two ways to calculate the undulation: Equation 7, Fast Fourier Transform (FFT), and least square collocation. The FFT algorithm needs regularized gridded data to optimize the computation's quality [12] Gridded data in the FFT algorithm were gathered using Trilinear Interpolation embedded in python's matplotlib plugin. The overall geoid model was calculated remove-restore method. The contributive of short and long-wavelength components were removed from the measured terrestrial gravity and restored on the final geoid. The step is written as in Equation 10 and Equation 11 [13].
where: N = total geoid undulation = long-wavelength geoid = residual geoid ℎ = short-wavelength geoid (terrain correction values) Δ = total free-air gravity anomaly Δ = long-wavelength gravity anomaly Δ = residual gravity anomaly Δ ℎ = short-wavelength gravity anomaly (terrain correction values) For this research purpose, we made use of about 707 gravity measurement data over the Greater Jakarta Area. Those data consisted of 2010 data which was measured by using Lacoste & Romberg gravimeter (1.5 km spatial resolution), and 2019 data, which was gathered by Scintrex CG-5 gravimeter (2.5 km spatial resolution).

Results and Discussion
Gravity measurement data were subtracted by the GGM gravity values as mentioned in the explanation about the remove-restore method and then interpolated into 1.5 x 1.5 km grids. The total interpolated area has the dimension of 50 x 25 km (see Figure 3), spanned from 6.4 º S -6.06 º S and 106.47 º E -107.15 º E. This step is known as the removing step of the remove-restore geoid computation method. As shown in Figure 3, the residual gravity data has a large value (75.079 mGal) in the dark blue area of North Jakarta and gradually decreases in the yellow area of South Jakarta (31.376 mGal).
This gridded gravity measurement data will be used as the primary input in an FFT processing algorithm embedded in spfour modules on Gravsoft [12]. In this spfour module, the long-wavelength geoid component from GGM will be added back to the aforementioned residual gravity data. This step is called the restore step of the remove-restore method. The final product of the Greater Jakarta local geoid model is depicted in Figure 4. As shown in Figure 4, the geoid has a low value (17.875 mGal) in the dark blue area of the Tangerang area and gradually increases in the yellow area of the Bekasi area (19.176 m).
A verification procedure will be conducted to analyze the quality of our local geoid model. As many as 11 points spanning all over the model's boundaries will be selected as the validation values. Those 11 points should have both ellipsoidal height (h) and orthometric height (H). The distribution of those samples depicted in Figure 5, while the data about those validation points is shown in Table 1.    . On the other hand, the RMS value of this research's result is somewhat smaller than that result above. The complete verification results are shown in Table 2. For this reason, terrestrial gravimetric measurements are carried out in this region by adding a point outside the previously measured area. This study was used to determine the accuracy of the results of terrestrial gravimeter measurements when using an interval of 2.5 km because, in previous studies, the measurement interval was larger.

Conclusion
The local geoid model of the Jakarta area based on 2.5 x 2.5 km gridded gravity measurement data has a standard deviation of 0.095 m and RMS of 0.365 m. These results are somewhat ambiguous compared to our previous local geoid model based on 1.5 x 1.5 km gravity data. The ambiguity itself was possibly caused by natural factors such as land subsidence and road condition, which affected the gravimeter readings. Another possible ambiguity source is instrumental factor such as GNSS and gravimeter measurement precision.