Determination of gravity anomalies in Java, Indonesia, from airborne gravity survey

The Geospatial Information Agency of Indonesia (BIG) recently carried out an airborne gravity survey mission to support a reliable Indonesian geoid model. The gravity observations covered all the main islands of Indonesia. This paper presents a state-of-the-art for gravity anomalies derivation using airborne gravity mission in Java, Indonesia. The common gravity corrections for deriving the scalar free-air gravity anomalies along the flight trajectory had been estimated using GNSS-derived positions. The corrected data were then filtered using the FIR method in which the cut-off frequency had been predetermined by considering aircraft altitude, geological setting, and instrument’s accuracy. To assess the airborne gravity results, we compared them with the upward continued terrestrial gravity measurements. In addition, we performed crossover analysis and adjusted the estimated biases to the airborne gravity measurements. The accuracy of adjusted airborne gravity anomaly was estimated to 3.37 mGal. In conclusion, the airborne gravity mission provided valuable data needed for further geodesy and geophysics applications


INTRODUCTION
The geoid represents the Earth's equipotential gravity surface which defines the definitive physical height reference in geodesy applications, i.e., orthometric height. This surface aligns with an ideal and undisturbed ocean surface, where tides, winds, and current influence were removed (Condi and Wunsch 2004;Amin et al. 2019). In addition, the geoid extends through the continents, creating a seamless surface omnipresent. It can be estimated using two common approaches, i.e., the geometric (Albayrak et al. 2020) and gravimetric (Vu et al. 2019) approach. Geoid computation using the geometric approach is a relatively straightforward method as it corporates the GNSS/leveling method to estimate the difference between well-known ellipsoid and geoid surface. On the other hand, gravimetric geoid modeling requires an extensive physical geodesy understanding, geological knowledge, and complex geoid computations. Furthermore, any geophysical influences, e.g., earth tides (Agnew 2015), ocean tides loading (Boy et al. 2003), atmospheric pressure loading (Boy et al. 2002), should be handled carefully from the data.
When using Stokes' formula to determine a gravimetric geoid, the gravity anomalies must be given continuously on the geoid surface around the world (Stokes 1849;De Witte 1966). However, since we typically measure gravity at or above the Earth surface, all gravity readings should be reduced to the geoid and satisfy a condition where there is no masses outside the geoid. The first bespoke condition, a continuous gravity anomaly surface, was a challenging task to undertake. Fortunately, thanks to the development of satellite gravimeters, such as The Gravity Recovery and Climate Experiment (GRACE), we are able to measure gravity throughout the Earth. The accuracy, however, may not satisfy engineering or scientific purposes that need a few Terr. Atmos. Ocean. Sci., Vol. 32, No. 5, Part II, 781-795, October 2021 centimeters accuracies of geoid model. Chen et al. (2018) reported that satellite gravimeter-based geoid accuracy ranges from 24 to 223 cm in Canada and Mexico. A service provided by the International Center for Global Earth Models (ICGEM) also corroborates this result that the accuracy of several satellite gravimeter-based geoid models in Australia, Brazil, Japan, and USA is spread from 40 to 300 cm. Moreover, the high spatial resolution of the geoid model is needed to mimic the short-scale geoid anomaly.
Numerous works have applied the so-called removecompute-restore (RCR) technique with a denser gravity data coverage in geoid modeling to enhance the geoid accuracy locally/regionally (Schwarz et al. 1990;Corchete et al. 2008;Yildiz et al. 2012). Using this approach, we first removed the long-wavelength and short-wavelength part from the gravity data observation using predicted gravity model from global/satellite gravity model and topography variation, respectively. Afterwards, the residual gravity anomalies were used to estimate the residual geoid model using proper Stokes' integration methods, e.g., Fast Fourier Transform (FFT;Schwarz et al. 1990) or Least Square Collocation (LSC; Tscherning 2015). Completing the RCR technique, both long and short-wavelength geoid signals were restored to the geoid model.
Several works have demonstrated that few centimeters geoid models is achievable most in lowland, flatter, and smaller areas where dense and uniform gravity data are available (Smith et al. 2013;Oršulić et al. 2020). However, in a rougher terrain condition and larger areas, e.g., in Indonesia, dense terrestrial gravity measurements may be not possible. If this is the case, measuring dense gravity networks covering the entire region will take more time and expenses. The airborne gravimeter allows us to traverse to places that are difficult or impossible to reach. This includes remote areas or mountainous areas with limited road access and areas mixed with land and ocean. Taking into consideration, the Geospatial Information Agency of Indonesia (BIG) started to deploy the airborne gravimeter covering Indonesia's main islands in 2008 (Pahlevi et al. 2015). During the first period of airborne gravity data acquisition (2008 -2010), BIG, together with Technical University of Denmark (DTU) completed airborne gravity measurements in Kalimantan, Sulawesi, and Papua islands. Eight years later (in 2018), BIG, with the help from National Chiao Tung University (NCTU), continued airborne gravity measurements in Sumatra island (Pahlevi et al. 2019). BIG finished airborne gravity measurements in Java in 2019 which covered all of Indonesia's major islands. LaCoste and Romberg (LCR) System II air/sea type gravimeter S-130 was used to perform gravity acquisitions in this study.
Although airborne gravimeter has several advantages in terms of time and cost efficiency, it possesses several challenges in data acquisition and processing. For example, one states an accuracy of 1 mGal of the observed gravity, but it may be worse than declared. The airborne gravity survey's accuracy is subject to several aspects, e.g., kinematic position of aircraft and error models of the measurement system (Hwang et al. 2007). Random errors further presented in the data caused by air turbulence and gravity sensor data containing systematic errors (Becker et al. 2015). The high-frequency random errors caused by air turbulences could be removed by using a proper low pass filter (Childers et al. 1999;Forsberg et al. 2000;Willberg et al. 2020). However, the resolvable wavelength might be longer than the sampling interval due to the filtering and averaging process. The systematic errors were first identified using the crossover analysis (Forsberg et al. 2000), and crossover adjustment was often used to correct the bias (Hwang et al. 2006(Hwang et al. , 2007. Hence, this paper presents the recently observed gravity anomalies using an airborne gravimeter with application to Java, Indonesia. We focused on the following issues: (1) software development to estimate gravity anomaly; (2) filter design and estimated resolvable wavelength; and (3) data accuracy assessment.
The following is how the paper is structured: section 2 reviews particulars aspects of utilizing airborne gravimeter to estimate gravity anomaly, section 3 presents the research findings and discussions, and lastly, section 4 provides the research summary and further considerations.

Gravity Readings
The zero-length spring principle is commonly used in any static platform relative gravimeter, and this also applies in the LCR air/sea type gravimeter (Valliant 1992). Since LCR air/sea type gravimeter is performed in a kinematic platform, any acceleration will affect gravity readings. The gravity reading eventually consists of gravity and any acceleration disturbances, and the effect reaches up to 100000 times larger than the gravimeter desired accuracy (Lacoste 1967). Neglecting the horizontal acceleration, one may write the characteristics of the zero-length spring as follows (Neumeyer et al. 2009): where g stands for the gravity, z p is the presented vertical acceleration on the data, B is the displacement of the mass test relative to its reference position, B o and B p are the first and second time derivative of B, while S refers to the spring tension. Assuming a linear gravimeter characteristic, b, f, k, and c are constants. The factor k is zero for the LCR air/sea type gravimeter as there is no restoring force of the spring.
The damping factor f is set to a really high value, but it remains constant. Eventually, the beam will reach its maximum velocity in a moment and b can also be ignored. We can further rewrite Eq. (1) as follows: Recalling that the measurement of gravity is carried out on a kinematic platform, the cross-coupling effect (CC) must be taken into account. The deflection of the gravimeter's lever arm caused by the horizontal acceleration adds additional vertical acceleration to the gravity readings. Neumeyer et al. (2009) reported that the cross coupling gravity effect spreads ±20 mGal. LCR air/sea type gravimeter implements five cross coupling monitors where each of the five monitors obtain each specific outputs. Outputs are typically extracted from the position and speed of the beam with the horizontal acceleration in cross and longitudinal axis direction. Finally, CC can be approximated as follows (Valliant 1992): with VE, VCC, AX, AL, and AX2 are the assigned variables for the five cross coupling monitors' output. C VE , C VCC , C AX , C AL , and C AX2 are the weight factors provided by the "Microg LaCoste". These weight factors were estimated using the cross correlation technique (Lacoste 1973). Figure 1 shows the cross-coupling gravity effects in airborne gravimetry data. By elaborating Eq.
(2) and CC, gravity reading can be calculated with: with Bv is the beam velocity, ST and ST c are the spring tension and its coefficient, respectively, and Sc is the gravity scaling factor.

Airborne Scalar Gravimetry Principle
The simplified equation for the airborne scalar gravimetry involves several additional corrections which leads to: with subscript i and 0 in g r d are the number of gravity readings and base gravity reading, respectively, the term a v p refers to the aircraft's vertical acceleration gravity effects (positive value at zenith direction), g T d is the platform's misalignment correction, and g 0 is reference gravity value at the apron or other specified point which is typically done using absolute gravimeter, e.g., FG5 or A10-type absolute gravimeter.
The motion of the aircraft is described using a rotating Earth coordinate system. Hence, the aircraft's vertical acceleration follows the following equation (Harlan 1968): where the first term on the right-hand side is the aircraft's acceleration with r and t are the position vector and time, respectively. The mentioned aircraft position can be precisely estimated using the Precise Point Positioning (PPP) or kinematic differential method on Global Navigation Satellite System (GNSS) technology (Novák et al. 2003;Hwang et al. 2007;Xincun et al. 2014). The third term refers to the acceleration induced due to the Earth's rotation speed (~) variations. According to the length of day (LOD) variations provided by the International Earth Rotation Service (IERS), the LOD variations vary between -1 to 4 milliseconds (see Fig. 2). Assuming the LOD is set to 5 milliseconds, the gravity effect due to ~ variation reaches 3 × 10 -4 mGal. This effect is considered small relative to the accuracy of airborne gravimetry. Thus, the ~ variation induced gravity effects can be neglected. The remaining second and fourth term of its vertical component refer to the Eötvös effects. Eötvös effect is the observed gravity change when the measurement is done in kinematic platform. Since aircraft moves relatively to the rotating Earth, both centrifugal and Coriolis acceleration occur and impact the gravity observation. When ground speed is considered, Harlan (1968) defines the Eötvös correction as follows: with v g is the ground speed, a is the semimajor axis of the reference ellipsoid, h is the altitude of the aircraft, z is the corresponding latitude, a is the aircraft heading relative to the North direction, and Ẽ is the Earth angular velocity.
The variable e in Eq. (7) presented as follows: The platform's misalignment gravity effects give substantial impact on gravity reading, in which can reach up to several hundreds of mGal on given trajectory (Neumeyer et al. 2009). To correct this platform's misalignment, one may apply the direct method by Peters and Brozena (1995): with X acc and L acc are the acceleration derived from the accelerometers in cross and longitudinal directions. The additional X accGNSS and L accGNSS are the same as the mentioned acceleration but are derived from the GNSS. g T d is the median of the misalignment effect.

Time Synchronization of Airborne Gravity Data
Multiple hardware was involved in the airborne gravity surveys, i.e., gravimeters and GNSS. The lag time from each corresponding output leads to a systematic error in airborne gravity data processing, and the application of the corrections mentioned in Eq. (5) may worsening the data. To correct the lag time, the maximum similarity of two data series as a function of one relative to the other needs to be estimated first using a cross-correlation technique (Olesen 2002;Neumeyer et al. 2009). For continuous functions of two data series u(t) and v(t), the cross-correlation is as follows (Papoulis 1962): and x is known as lag or displacement. Here, we followed the procedure mentioned by Olesen (2002) to decimate the 1 Hz data to 10 Hz before applying the cross-correlation technique. Figure 3 displays one example of cross-correlation application between gravity readings and calculated vertical accelerations at the selected observation period. The lag time between the two datasets was estimated to 1.7 seconds, and we used the estimated lag time to correct the corresponding flight's time.

Free-air Gravity Anomaly Reduction
The free-air gravity anomaly is principally given by the following equation: where c is the normal gravity at specified reference ellipsoid and latitude, g FAC d is the free-air correction, and g atm d is the atmospheric correction.
Normal gravity is also known as the latitude correction. Its correction varies as a function of latitude caused by the centrifugal acceleration. To correct it, we applied a standard gravity normal calculation of the Somigliana-Pizetti closed formula as given by (Moritz 1980 where e c is the equatorial normal gravity value, e is the first numerical eccentricity of the specified ellipsoid, and z is the latitude of observations. The term k stands for the normal gravity constant, which can be described as: with p c is the polar normal gravity value and b is the semimajor axis of the reference ellipsoid. Free-air correction is made to partially downward or upward-continue the gravity observation to a reference datum or geoid. It is common to use a linear vertical free-air gradient of 3.086 mGal m -1 to approximate such correction. However, to compensate for the ellipsoid representation of earth, a second-order correction was used instead. This correction was derived by Featherstone (1995) as: with f is the flat earth and H is the orthometric heigh. The last mentioned m can be calculated as follows: where GM is the product of the gravitational constant and the mass of the earth. The difference between the conventional free-air correction and second-order free-air correction in Java, Indonesia reached 0.5 mGal. Note that this difference was computed in the height of 4200 m which is the mean altitude of the airborne gravity survey in Java, Indonesia. The fourth term on the right-hand side in Eq. (11) arose due to the systematic error that existed during the determination of ellipsoidal parameters. The GM was computed using satellite-based geodetic data, causing the earth's atmospheric masses' inclusion in the estimated parameter (Featherstone and Dentith 1997). Therefore, an additional atmospheric correction must be added to the observed gravity (Rapp and Pavlis 1990). Equation (16) is the given approximation of the atmospheric correction in mGal outputs by Yang (2013).

Time-Series Filtering
It has been mentioned earlier that airborne gravimetry is subject to high-frequency noises. In fact, the high-frequency spectrum is sometimes not associated with noise only in a   dynamic airborne gravimetry system. In addition, the aircraft motion in which the aircraft pitches up and followed by pitches down resulting in speed variation accordingly also contributes lower frequency noises (Li and Jekeli 2004). To account for these problems, Zhao et al. (2015) implemented the Empirical Mode Decomposition (EMD) in their data processing. Meanwhile, others used the wavelength approaches to improve the resolution of estimated gravity anomalies (Li and Jekeli 2004). Although these approaches do not require assumptions about the nature of the signal, they are relatively more complex than the single direct cutoff method. In addition, the low-pass filtering is comparable with a relatively slight RMS difference of 0.3 mGal when Zhao et al. (2015) compared it with the EMD method.
Presently, we are implementing the single cut-off method using Finite Impulse Response (FIR) to tackle highfrequency noises. The advantages of utilizing FIR include the following: it is relatively easy to design and requires less computational time. FIR can be mathematically described as a convolution of an input signal x with a filter impulseresponse h yielding (Oppenheim et al. 1999): with N is the filter order. Let h d (k) be an ideal frequency selective filter's unit sample response, the frequency response can be given by: By evaluating the inverse Fourier transform to the desired frequency response ( ) H d~, we obtained the infinite duration of the sample response h d (k). Therefore, it must be truncated at some extent that can be performed using some suitably window functions w(k), yielding h(k) = h d (k) × w(k). The different types of window functions are the Rectangular, Hanning, Hamming, Blackman, Bartlett, and Kaiser.

Software Implementation
We developed a Graphical User Interface (GUI) software called AiNG-ITB (Airborne Gravity Data Processing-Institut Teknologi Bandung) to apply the mentioned corrections and reductions for airborne gravimetry data processing. AiNG-ITB is written entirely in Matlab® language and allows the users to obtain free-air gravity anomaly at flight altitude. Moreover, it features the ability to process the data in a single file or batch mode. The main window of AiNG-ITB is displayed in Fig. 4.

Data Inputs and Outputs
The required input files for AiNG-ITB are the raw output of the LCR air/sea gravimeter (*.dat file) and the output from the GrafNav software (*.txt file). The raw gravity reading data consists of all the needed variables for data processing, such as the beam position, VE, VCC, AX, AL, and AX2. The position-dependent gravity reductions (i.e., normal gravity, atmospheric, and free-air correction) are estimated using the estimated positions from the GrafNav's output. Further, the aircraft's dynamic gravity effects are also estimated using the mentioned output.
There are two main outputs for AiNG-ITB: (1) processed data file and (2) summary report. The processed data file (with a file extension of *.out) lists the gravity anomalies (with one additional output of gravity disturbances) at their corresponding location. It also comes with some data processing summary, i.e., the standard deviation of tilt and their differences with the chosen Global Geopotential Model (EGM2008). A Portable Document Format (PDF) file comes along with the processed data file. It functions as a quick user's guide for describing the range of data acquisition/processing performance, in which several figures, such as time synchronization, tilt correction, and the comparison with GGM are shown in this report.

Things to Consider when using AiNG-ITB
Prior to using AiNG-ITB, the following should be defined in the data processing inputs: (1) Reference or absolute gravity value at the reference point. It is commonly obtained using the absolute gravimeter of FG5 or A10. Otherwise, it can also be obtained using a relative gravimeter. However, an absolute-type gravimeter is preferred to ensure the highest accuracy of the gravity reference point.
(2) Offset between the GNSS antenna height and gravimeter. The GNSS antenna and the gravimeter are not at the same level. Eventually, the free-air mass lies between them and needs to be corrected. Considering the two meters of height difference between antenna and gravimeter, the free-air gravity gives gravity effects of approximately 0.6 mGal and leads to a systematic error in further data analysis. (3) Cut-off frequency. It has been described earlier that a kinematic platform (in this case using the aircraft) experiences substantial noises in higher frequency spectrum. Therefore, a proper cut-off frequency must be defined by considering several aspects (will be described in section 3.2). We implemented two window functions for the FIR implementation in the software. They were the Hamming and Blackman windows which are user selectable.

Identifying the Optimum Filter of Airborne Gravimetry Use
The optimum filter, in term of cut-off frequency, in the airborne gravimetry data processing holds a crucial aspect in the data processing. Some errors which occur when selecting the cut-off frequency can cause loss of expected anomaly gravity signal or even allowing high-frequency noise in the filtered data. The shortest observable anomaly wavelength defines the selection of an ideal cut-off frequency from the aircraft altitude. Additionally, it depends on the hardware and environment conditions (Wei and Schwarz 1998;Childers et al. 1999).
According to Newton's universal gravitation law, the force acting between two objects is proportional to the masses and the inverse distance between their masses' centers. Therefore, the smallest observable anomaly wavelength in airborne gravimetry will eventually attenuate as the altitude increased. In line with the bespoke Newton's universal gravitation law, the study area's geological setting also contributes to the corresponding wavelength. Assuming a low contrast density in the study area, we may obtain relatively low gravity responses. Thus, a longer observable wavelength is obtained in such conditions. The last factor to be considered in identifying the observable wavelength is the accuracy of the gravimeter.
We followed the procedure proposed by Childers et al. (1999) to estimate the appropriate cut-off frequency used for the recently deployed airborne survey. First, the smallest observable gravity effects needed to be set. Accommodating the noise inherent in the measurement, any geological setting that lies beneath the earth's surface generates at least three mGal of gravity effects. Three mGal was chosen according to the three-sigma-rule of thumb for the claim of 1 mGal of LCR sea/air type gravimeter's accuracy. Afterwards, the expected contrast density that lies beneath the earth's surface needed to be defined. Indonesia's geological structure, particularly in Java, is formed by various Geological formations, i.e., dominated by Cenozoic formation and covered with volcanic and quaternary formation (Darman and Sidi 2000). The expected contrast density in Java was set to 1100 kg m -3 . It followed the low and highdensity zone ranging from 1700 to 2800 kg m -3 from previous research (Tiede et al. 2005;Febriani 2014). Finally, a rearranged formulation of simple sphere gravity effects was used to determine the dimension of buried geological structures. It is defined as follows [see details in Childers et al. (1999)]: with v is the expected contrast density, Δg min is the smallest observed amplitude of anomaly as characterized by the instrument's accuracy, z is the mean aircraft altitude, and R is the dimension of buried geological structures. An iterative procedure was applied to estimate the R until the desired contrast density value obtained. Afterward, the geological wavelength was estimated as a function of distance by g m = 1.54(z + R). In terms of the frequency domain, the 'Fourier' wavelength was defined by f m = 3.1(z + R). With application to the recent Indonesia gravity survey, the corresponding wavelength was estimated using the mentioned approach. The corresponding 'Fourier' wavelength using an average aircraft altitude of 4200 m was estimated to be 14.9 km. Using an average aircraft speed of 70 m s -1 , the proposed cut-off frequency utilized in this airborne gravimetry survey was set to 0.0047 Hz. Figure 5 shows the calculated free-air gravity anomaly at the selected flight trajectory. The unfiltered free-air gravity anomaly varied ±10000 mGal. At epoch 13000 to 13750, we saw that the GNSS-derived kinematic acceleration could not entirely reduce air turbulences' noises and led to a relatively high gravity variation compared to other periods. This was also mentioned by Lu et al. (2017) that gravity effects due to turbulence may not be entirely reduced using GNSS-derived kinematic acceleration correction. In the frequency domain, we observed that noises were dominating the unfiltered free-air gravity anomaly time-series data. Its amplitudes reached up to 2800 mGal extending from 0.05 to 0.5 Hz. Furthermore, the application of FIR tackled approximately 98% of observed gravity amplitude, and the filtered free-air gravity anomaly ranged between 30 and 160 mGal at the selected trajectory.

Derived Free-air Gravity Anomaly
We used the modeled anomaly from Earth Gravitational Model 2008 (EGM2008) by Pavlis et al. (2012) to assess our observed airborne gravity anomaly. Figure 5 shows the overall trend of the observed gravity follows the EGM2008's trend. The standard deviation of their differences was estimated to 12.6 mGal. We used EGM2008 to evaluate the long-wavelength component of the observed gravity anomaly, as the short-wavelength component might not be appropriate for the Indonesia because of the lack of terrestrial gravity observations. Sarsito et al. (2020) reported that gravity anomalies in Indonesia's land areas are mostly based on a forward modeling approach using a terrain model. Eventually, we expect our gravity survey will improve the short-wavelength gravity component over the study area. Figure 6 shows free-air gravity anomalies along all tracks at their corresponding altitude. They vary from -220 to 250 mGal. Compared to the topographical variations, most of the high gravity anomalies are associated with mountainous landscapes. However, the high gravity anomaly seems not to be correlated with the height increment at the south of Jakarta (i.e., 106.5°E and 7.2°S or south coast of Java). It suggests that the density at the vicinity point of interest is relatively larger than surrounding, i.e., due to the existence of ultramafic materials that were produced as a result of the submergence of the Indian oceanic plate under the Eurasian continental plate (Gunawan et al. 2019;Widiyantoro et al. 2020). Several areas were not covered by airborne gravimetry due to flight permit limitation, such as Jakarta. Terrestrial gravity observations will be used to fill the lack of measurements for further gravity application, i.e., geoid modeling.

Comparison with Terrestrial Gravity Data
Quality assessment of airborne gravity survey was begun by comparing the terrestrial gravity with airborne grav-ity data. BIG collected the existing terrestrial gravity data at several major cities in Java, Indonesia (see Fig. 7). Gravity measurements were made using Scintrex CG-5/CG-6 gravimeter, with a distance between observation points of approximately 5 km. In Jakarta, denser gravity acquisitions were able to collect with approximately 2 km spacing. The comparison was performed by first removing the long-wavelength gravity anomaly component in terrestrial gravity measurements using the GGM. Subsequently, a residual surface was computed by applying a least-squares collocation (LSC; Tscherning 2015) on the residual gravity anomaly data. The residual surface was then upward continued to the mean flight altitude using 2D-FFT as (Verdun et al. 2003;Hwang et al. 2007): where F z and F 0 are respectively the gridded upward continue surface and reference surface in spectra, u and v are spatial frequencies. Finally, the long-wavelength component of GGM was restored to the residual surface to obtain the upward continued gravity anomaly surface.
Since terrestrial gravity data are only presented in several major cities, the comparison was performed at three selected cities, namely, Yogyakarta (Southern part of Java), Semarang (Northern part of Java), and Surabaya (Eastern part of Java). In the bottom of Fig. 7 shows the comparison of airborne and upward continued gravity data. The comparison revealed that the differences ranged ±40 mGal, with an average difference of 5.4 mGal and a standard deviation of 14.2 mGal. This average difference highlights possible errors in both airborne and terrestrial gravity measurements. Assuming there was no systematic errors in the data, we expected to get a constant average difference for each zone. However, this is not the case of this research. We obtained different average residues for each zone, i.e., the average residual for Jogjakarta, Semarang, and Surabaya of -1.8, 10.3, and 9.8 mGal, respectively. Hence, the gravity reference points that were used need to be evaluated.
We also considered the upward continuation computation errors. The classical Poisson integral works under several assumptions, i.e., a uniform assumption of gravity gradient and rock density (Cruz and Laskowski 1984). Given the complexity of Indonesia's geological setting, this assumption may not work well in Indonesia. Supriyadi et al. (2020) reported that Semarang had the gravity gradient varied from 1.5 to 3 mGal m -1 in July 2019. That case becomes more complex and the gravity gradient also varies as a function of time due to hydrological variations. Concerning the rock density, the Yogyakarta region is dominated by three materials (i.e., unconsolidated materials, tertiary sediments, and igneous rocks) with a density ranging from 2.0 to 2.8 kg m -3 (Nurwihastuti et al. 2014).

Crossover Analysis and Adjustment
We also use crossover difference analysis to evaluate the accuracy of airborne gravity data. The crossover analysis is performed by comparing gravity values at the crossover point of two intersecting airborne gravity trajectory lines. Since no gravity observation is likely done in the exact locations, interpolation is needed to be done at crossover points. Crossover differences are defined by subtracting the North-South lines' interpolated data with the East-West lines' interpolated data to make it consistent.
A total of 113 crossover points was collated for crossover analysis. Figure 8 displays the crossover differences ranging from -31.44 to 32.85 mGal. These results are considered high if compared to the other research (Hwang et al. 2007). However, we considered these differences due to random errors as its histogram follows the normal distribution (Fig. 9). Random errors possibly occurred due to the air turbulence. If we considered an additional 'Morlet' shape noise present in the middle of observation time due to air turbulence, it led to an additional long-wavelength signal existing in the pre-filtered data. We used the low-pass filter to tackle high-frequency noises. Hence, the mentioned additional long-wavelength signal still presented in the filtered data. The accuracy of GNSS-derived positions possibly caused these random errors. It affected the estimated gravity correction used for airborne gravimetry, e.g., vertical acceleration and Eötvös correction. Hwang et al. (2007) mentioned that the required horizontal position accuracy for East-West direction is about 4 cm for the Taiwan case. Only a single GNSS station was used as a reference for each airborne gravity survey in this research. The corresponding baseline reached up to 250 km. It may affect GNSS-derived positions' accuracy since GNSS accuracy will be poorer as a distance function (Xu 2007;Hofmann-Wellenhof et al. 2008;Bramanto et al. 2019).
We tried to reduce any possible systematic errors using a crossover adjustment. A detailed explanation of crossover adjustment can be seen in Hwang et al. (2006). Each airborne gravity survey lasts for a few hours, and the reported drift of the LCR sea/air type gravimeter was relatively low compared to the survey period (3 mGal month -1 ). Hence, the effect of gravimeter drift was excluded, leaving a parameter of bias during the computation. We set three North-South tracks, i.e., lines 3, 9, and 73 to be fixed for the computation. These tracks have relatively low turbulence during the flight. Moreover, these tracks also give relatively slight variances for the tilt corrections. The application of crossover adjustment is presented in Table 1 and Fig. 10. Prior to the application of crossover adjustment, the standard deviation of crossover differences was estimated to 12.54 mGal. We also observed that four crossover differences were larger than ±12 mGal after applying crossover adjustment. These significant discrepancies were probably caused by sudden turbulence and poor GNSS accuracy at these corresponding tracks. Besides these outliers, the application of crossover adjustment reduced the crossover variance difference of approximately 62% from its original.

SUMMARY AND OUTLOOK
This paper presents an airborne gravimetry-derived free-air gravity anomaly from the recent airborne gravity survey in Java, Indonesia. The resolvable 'Fourier' wavelength     10. Same as Fig. 9, but after the application of crossover adjustment.
of airborne gravity survey was estimated to be 14.9 km. By using this information, we expected our airborne gravity survey to capture the corresponding 7.4 km 'geologic' wavelength. It equals to an object with a radius of 0.6 km based on the estimated R in Eq. (19). Although the smallest observable wavelength was defined by three factors (i.e., aircraft altitude, contrast density, and instrument's accuracy), we considered using a single value for the cutoff frequency rather than the adaptive value. It ensures the consistency for obtained resolvable wavelength at all tracks. The standard deviation of differences at crossover points for adjusted gravity anomaly was 4.76 mGal. Hence, along-track accuracy was calculated to 4.76/ 2 = 3.37 mGal according to the law of error propagation. The accuracy information of derived gravity anomalies was necessary for numerical aspects of regional geoid modeling. It helped us reviewing and calculating the formal estimates of the representation error for regional geoid modeling in future milestones. Concerning the geoid computation, we suggest that an extensive investigation is needed for the terrestrial gravity observations. We found that systematic errors possibly exist in terrestrial gravity observations based on the comparison between upward continued terrestrial gravity and airborne gravity observations.
We highlight some topics of future improvements for airborne gravity data processing. First, the present software applied a straightforward low-pass filter of FIR to the data. In ideal occasions, comparisons of other low-pass filter options give similar results (Hammada 1996;Olesen 2002). However, air turbulence creates sudden jumps at time-series data and lead to an additional long-wavelength gravity signal. Therefore, various options to filter on noisy airborne gravity data is an interesting subject for research. For example, the iterative Gaussian filter (Hwang et al. 2006), the second-order Butterworth filter (Olesen 2002), EMD (Zhao et al. 2015), and wavelets (Li and Jekeli 2004) methods. Moreover, the accuracy of GNSS-based positions should be analyzed further. PPP might be a good candidate since it can achieve few cm positioning accuracies. In addition, it is not related to the baseline length. A research by Lu et al. (2017) mentioned that the velocity and acceleration estimated from raw GNSS Doppler observations having better performance than GNSS carrier phase observations to correct airborne gravity disturbances. If any airborne gravity survey is conducted in the future, we suggest attaching several GNSS antennas on the aircraft and the Inertial Measurement Unit (IMU). These instruments can be used to calculate the flight-state of an aircraft (velocity, acceleration, heading, pitch, and roll) better. Eventually, the corrections, subject to the aircraft movement, can be estimated precisely.
To summarize, the recent airborne gravity survey is considered as a successful mission. It complements the need for evenly distributed gravity data in Java, Indonesia. The results can be further used in any geodesy or geophys-ics applications, including the geoid or geological structural modeling.