Gravimetric geoid modeling in the northern region of Peninsular Malaysia (NGM17) using KTH method

In this study, a new geoid model for the northern region of Peninsular Malaysia (NGM17) was computed using an alternative method known as the Least Squares Modification of Stokes formula (LSMS) with Additive Corrections (AC) or commonly called the KTH method. The NGM17 geoid was derived from the recent terrestrial gravity data provided by Department of Surveying and Mapping Malaysia (DSMM), the most recent global digital elevation model ALOS World 3D (AW3D-30) GDEM, global geopotential model (GGM) derived from three satellite gravity missions, marine gravity anomalies extracted from DTU 10 Global Gravity Field and WGM2012 Earth’s gravity anomalies. The gravimetric geoid model derived in this study (NGM17) as well as the geoid obtained from DSMM were then evaluated against the GNSS-levelling data. The statistical analysis obtained shows that NGM17 gives slightly better accuracy with the mean error of NGM17 and DSMM geoid model were 0.2568m and 1.1648m respectively, and the RMSE of ±0.2686m and ±1.1656m respectively.


Introduction
Various approaches have been proposed by geodesists to achieve the goal of determining precise geoid model with each different approach is based on its own technique and philosophy. One of the techniques is Remove-Compute-Restore (RCR) that was successfully applied in a number of geoid computation tasks [1], [2]. Recently, a version of the KTH method has been introduced and has attracted becoming a popular method among the geodetic community. The method was successfully implemented in Iran [3], Tanzania [4], Central Turkey [5] and Sweden [6]. One of the characteristics which made KTH method differs compared to other methods is that the KTH approach matches the error of truncation, gravity anomaly and the Global Geopotential Model (GGM) [7].
In geoid computation, a number of data sources are needed and Digital Elevation Model (DEM) being one. It is used in the computation of Bouger gravity correction and the downward continuation effects. The Shuttle Radar Topography Mission (SRTM) DEM is the most well-known GDEMs and usually used in geoid computation whether in R-C-R method as well as in KTH method. The most recent release of the GDEM is ALOS World 3D (AW3D-30), released by the Japan Aerospace Exploration Agency (JAXA). According to [8], this GDEM generated based on the optical stereo matching technique. Based on preliminary validation using 5,121 independent checkpoints distributed around the world, the accuracy claimed of AW3D-30 is 4.40 m [9].
The main aim of the study reported in this paper is to model a geoid model over northern region of Peninsular Malaysia by utilizing AW3D-30 GDEM and the geoid is to be computed using the KTH method. Among the motivations of carrying out this work is that this being the first study that used the AW3D-30 to determine geoid model in Peninsular Malaysia. The focus area for geoid computation in this study is northwestern of Peninsular Malaysia (Figure 1).

Gravity data
The land gravity data were measured at 6047 points as shown in Figure 2. In addition to land gravity data, WGM2012 surface free-air anomalies were used for the area which not cover by the land gravity data. The WGM2012 data are freely available and details about the gravity anomalies can be found in the paper by [10]. Next, the DTU10 Global Gravity Field is utilized to cover free-air gravity anomalies at ocean area. The DTU10 data was derived from satellite altimeter by using double re-tracking technique in analyzing the waveform [11]. The data was limited to ocean area only since the quality of altimetry data is not standardized for land area. Therefore, the gravity anomalies on land area were removed.

Figure 2.
Location of land gravity data recent release of the GDEM is ALOS World 3D (AW3D-30), released by the Japan Aerospace Exploration Agency (JAXA). According to [8], this GDEM generated based on the optical stereo matching technique. Based on preliminary validation using 5,121 independent checkpoints distributed around the world, the accuracy claimed of AW3D-30 is 4.40 m [9].
The main aim of the study reported in this paper is to model a geoid model over northern region of Peninsular Malaysia by utilizing AW3D-30 GDEM and the geoid is to be computed using the KTH method. Among the motivations of carrying out this work is that this being the first study that used the AW3D-30 to determine geoid model in Peninsular Malaysia. The focus area for geoid computation in this study is northwestern of Peninsular Malaysia (Figure 1).

Gravity data
The land gravity data were measured at 6047 points as shown in Figure 2. In addition to land gravity data, WGM2012 surface free-air anomalies were used for the area which not cover by the land gravity data. The WGM2012 data are freely available and details about the gravity anomalies can be found in the paper by [10]. Next, the DTU10 Global Gravity Field is utilized to cover free-air gravity anomalies at ocean area. The DTU10 data was derived from satellite altimeter by using double re-tracking technique in analyzing the waveform [11]. The data was limited to ocean area only since the quality of altimetry data is not standardized for land area. Therefore, the gravity anomalies on land area were removed.

Figure 2.
Location of land gravity data recent release of the GDEM is ALOS World 3D (AW3D-30), released by the Japan Aerospace Exploration Agency (JAXA). According to [8], this GDEM generated based on the optical stereo matching technique. Based on preliminary validation using 5,121 independent checkpoints distributed around the world, the accuracy claimed of AW3D-30 is 4.40 m [9].
The main aim of the study reported in this paper is to model a geoid model over northern region of Peninsular Malaysia by utilizing AW3D-30 GDEM and the geoid is to be computed using the KTH method. Among the motivations of carrying out this work is that this being the first study that used the AW3D-30 to determine geoid model in Peninsular Malaysia. The focus area for geoid computation in this study is northwestern of Peninsular Malaysia (Figure 1).

Gravity data
The land gravity data were measured at 6047 points as shown in Figure 2. In addition to land gravity data, WGM2012 surface free-air anomalies were used for the area which not cover by the land gravity data. The WGM2012 data are freely available and details about the gravity anomalies can be found in the paper by [10]. Next, the DTU10 Global Gravity Field is utilized to cover free-air gravity anomalies at ocean area. The DTU10 data was derived from satellite altimeter by using double re-tracking technique in analyzing the waveform [11]. The data was limited to ocean area only since the quality of altimetry data is not standardized for land area. Therefore, the gravity anomalies on land area were removed.

Global Geopotential Model (GGMs)
The use of GGM data is needed in geoid computation scheme using the modified Stokes formula to provide long-wavelength contribution of gravity anomalies. In practice the GGMs from satellite-only (Goce, Grace and Lageos) missions were used and details about the GGMs used in this study are shown in Table 1. A total of six (6) GGMs were used in this study and the evaluation was conducted to choose the best fit GGMs for the study area. The evaluation is carried out by comparing the GGMs derived geoid height (N GGMs ) with GNSS/levelling derive geoid height (N GNSS ) based on the model:

Global Digital Elevation Model (GDEM)
Another source of data which become requirement in geoid computation is the digital elevation model. In this study the DEM chosen to be used is AW3D-30 GDEM and the elevations from this DEM (H DEM ) were compared with local mean sea level (H msl ) to evaluate the accuracy of DEM. The interpolation of H DEM was implemented using the Krigging approach which is available in SURFER package (http://www.goldensoftware.com/). The GDEM has also been compared with SRTM accuracy and the DEM model can be downloaded via USGS website (http://earthexplorer.usgs.gov).

GNSS/Leveling Data
The GNSS/leveling data are the importance source of data used to analyze the accuracy of GDEM as well as the GGM. The verification is required for the purpose of choosing the best model to be adopted for geoid determination. The same data set (i.e., GNSS/Levelling) will be used once again to assess the accuracy of model resulted from the computation. A total of 38 GNSS/levelling data is used in this study. The GNSS observations were conducted on the Standard Benchmark (SBM) and Benchmark (BM) (see Figure 3) using static method. By using static method, the data were processed using TOPCON Tool software. In the processing, two MyRTKnet stations were used as the control stations. Figure 3. Distribution of validation point

Determination of gravimetric geoid using KTH method
The first step of implementing geoid computation using KTH method is the reduction of gravity anomalies. The surface gravity anomalies (Δg) were computed as;

Global Geopotential Model (GGMs)
The use of GGM data is needed in geoid computation scheme using the modified Stokes formula to provide long-wavelength contribution of gravity anomalies. In practice the GGMs from satellite-only (Goce, Grace and Lageos) missions were used and details about the GGMs used in this study are shown in Table 1. A total of six (6) GGMs were used in this study and the evaluation was conducted to choose the best fit GGMs for the study area. The evaluation is carried out by comparing the GGMs derived geoid height (N GGMs ) with GNSS/levelling derive geoid height (N GNSS ) based on the model:

Global Digital Elevation Model (GDEM)
Another source of data which become requirement in geoid computation is the digital elevation model. In this study the DEM chosen to be used is AW3D-30 GDEM and the elevations from this DEM (H DEM ) were compared with local mean sea level (H msl ) to evaluate the accuracy of DEM. The interpolation of H DEM was implemented using the Krigging approach which is available in SURFER package (http://www.goldensoftware.com/). The GDEM has also been compared with SRTM accuracy and the DEM model can be downloaded via USGS website (http://earthexplorer.usgs.gov).

GNSS/Leveling Data
The GNSS/leveling data are the importance source of data used to analyze the accuracy of GDEM as well as the GGM. The verification is required for the purpose of choosing the best model to be adopted for geoid determination. The same data set (i.e., GNSS/Levelling) will be used once again to assess the accuracy of model resulted from the computation. A total of 38 GNSS/levelling data is used in this study. The GNSS observations were conducted on the Standard Benchmark (SBM) and Benchmark (BM) (see Figure 3) using static method. By using static method, the data were processed using TOPCON Tool software. In the processing, two MyRTKnet stations were used as the control stations. Figure 3. Distribution of validation point

Determination of gravimetric geoid using KTH method
The first step of implementing geoid computation using KTH method is the reduction of gravity anomalies. The surface gravity anomalies (Δg) were computed as;

Global Geopotential Model (GGMs)
The use of GGM data is needed in geoid computation scheme using the modified Stokes formula to provide long-wavelength contribution of gravity anomalies. In practice the GGMs from satellite-only (Goce, Grace and Lageos) missions were used and details about the GGMs used in this study are shown in Table 1. A total of six (6) GGMs were used in this study and the evaluation was conducted to choose the best fit GGMs for the study area. The evaluation is carried out by comparing the GGMs derived geoid height (N GGMs ) with GNSS/levelling derive geoid height (N GNSS ) based on the model:

Global Digital Elevation Model (GDEM)
Another source of data which become requirement in geoid computation is the digital elevation model. In this study the DEM chosen to be used is AW3D-30 GDEM and the elevations from this DEM (H DEM ) were compared with local mean sea level (H msl ) to evaluate the accuracy of DEM. The interpolation of H DEM was implemented using the Krigging approach which is available in SURFER package (http://www.goldensoftware.com/). The GDEM has also been compared with SRTM accuracy and the DEM model can be downloaded via USGS website (http://earthexplorer.usgs.gov).

GNSS/Leveling Data
The GNSS/leveling data are the importance source of data used to analyze the accuracy of GDEM as well as the GGM. The verification is required for the purpose of choosing the best model to be adopted for geoid determination. The same data set (i.e., GNSS/Levelling) will be used once again to assess the accuracy of model resulted from the computation. A total of 38 GNSS/levelling data is used in this study. The GNSS observations were conducted on the Standard Benchmark (SBM) and Benchmark (BM) (see Figure 3) using static method. By using static method, the data were processed using TOPCON Tool software. In the processing, two MyRTKnet stations were used as the control stations. Figure 3. Distribution of validation point

Determination of gravimetric geoid using KTH method
The first step of implementing geoid computation using KTH method is the reduction of gravity anomalies. The surface gravity anomalies (Δg) were computed as; where g is the terrestrial gravity data and γ is the normal gravity. In this study, gravity anomalies adopted is based on Molodensky's solution as discussed in [18] where the normal gravity at the telluroid (let say at point γ p ) can be determined as; where γ is the normal gravity determined using the formula: γ = 978032.67715 (1+a 2 sin 2 φ+a 4 sin 4 φ+a 6 sin 6 φ+a 8 sin 8 φ) where a 2 =0.0052790414, a 4 =0.0000232718, a 6 =0.0000001262, and a 8 =0.0000000007 and φ is the latitude.
Meanwhile, the H N can be computed based on the formula: Here, the N-ζ is the geoid-to-quasi-geoid separation and determined using the formula [19]; where Δg B is the simple Bouguer gravity anomaly.
In order to archive outlier-free data, the surface gravity anomalies were analyzed using the crossvalidation approach [20]. In KTH method, the full gravity anomaly without any reduction will be used in geoid determination and it has become a challenge in interpolating unreduced gravity anomalies. Therefore, a special interpolation technique (see detail in [7]) has been used to grid the gravity data. According to [6], this interpolation technique reduces the discretization errors in the KTH technique. As shown in Figures 1 and 2, the study area is mostly surrounded by marine region and there is no data at the vicinity of Sumatra Island. These gaps were filled with WGM2012 surface gravity anomaly and DTU 10 Global Gravity Field to cover Sumatra Island and the marine region respectively. By using the method of Kriging, the final grid of surface gravity anomaly derived from combination of terrestrial gravity, WGM2012 and DTU 10 Global Gravity Field, was constructed.
The second step in geoid determination using KTH method is accounting the error degree variance using the white noise and the reciprocal distance model as suggested by [6]. The third step is computation of the modification parameters. As mentioned previously, the s n and b n are two sets of arbitrary modification parameters needed to be determined accurately. In general, the behaviour of these parameters is influenced by a set of conditions, namely choice of cap size, correlation length, and standard deviations of the terrestrial gravity. The conditions in determination of geoid for Peninsular Malaysia have been studied in detail [21]. As suggested by the study [21], the optimum condition parameter values for the cap size, correlation length, and standard deviations are 3.0º, 0.4º and 5.0mGal, respectively.
Computation of approximate geoid model is the fourth step in KTH method. In this study, a 1'x1' grid of the surface gravity anomalies from terrestrial gravity, DTU marine gravity and WGM2012 has been combined with the best fit GGM to produce approximate geoid height model for the northern region of Peninsular Malaysia as shown in Figure 4. At the final stage, the four different additive corrections (see Figure 4) have been computed separately and added to approximate geoid to yield final geoid as shown in Figure 5.

Assessment of GGM
In order to choose the best fit GGMs six (6) GGMs were evaluated to identify the best fit GGMs for the study area. The evaluation was implemented by comparing the geoid height form GGMs (N GGMs ) with the geometrical geoid derived from GNSS/levelling (N GNSS ) computed at 38 points based on the following: The geoid height accuracy is evaluated based on the computed Mean Error (ME) and Root Mean The result of statistical analysis conducted in comparing the geoid height from GGMs with geometrical geoid height derived from GNSS /levelling is presented in Table 2. The finding was go_cons_gcf_2_dir_r3 model performs well compared to other GGMs. As presented in Table 2, the ME values and RMSE for go_cons_gcf_2_dir_r3 were ±~0.2804m and 0.2473m, respectively. Meanwhile, go_cons_gcf_2_dir_r4 exhibited the largest error where the RMSE and mean error value for the GGMs model were calculated as ± 0.4263m and 0.4244m respectively.

Assessment of GGM
In order to choose the best fit GGMs six (6) GGMs were evaluated to identify the best fit GGMs for the study area. The evaluation was implemented by comparing the geoid height form GGMs (N GGMs ) with the geometrical geoid derived from GNSS/levelling (N GNSS ) computed at 38 points based on the following: The geoid height accuracy is evaluated based on the computed Mean Error (ME) and Root Mean The result of statistical analysis conducted in comparing the geoid height from GGMs with geometrical geoid height derived from GNSS /levelling is presented in Table 2. The finding was go_cons_gcf_2_dir_r3 model performs well compared to other GGMs. As presented in Table 2, the ME values and RMSE for go_cons_gcf_2_dir_r3 were ±~0.2804m and 0.2473m, respectively. Meanwhile, go_cons_gcf_2_dir_r4 exhibited the largest error where the RMSE and mean error value for the GGMs model were calculated as ± 0.4263m and 0.4244m respectively.

Assessment of GGM
In order to choose the best fit GGMs six (6) GGMs were evaluated to identify the best fit GGMs for the study area. The evaluation was implemented by comparing the geoid height form GGMs (N GGMs ) with the geometrical geoid derived from GNSS/levelling (N GNSS ) computed at 38 points based on the following: The geoid height accuracy is evaluated based on the computed Mean Error (ME) and Root Mean The result of statistical analysis conducted in comparing the geoid height from GGMs with geometrical geoid height derived from GNSS /levelling is presented in Table 2. The finding was go_cons_gcf_2_dir_r3 model performs well compared to other GGMs. As presented in Table 2, the ME values and RMSE for go_cons_gcf_2_dir_r3 were ±~0.2804m and 0.2473m, respectively. Meanwhile, go_cons_gcf_2_dir_r4 exhibited the largest error where the RMSE and mean error value for the GGMs model were calculated as ± 0.4263m and 0.4244m respectively.

Assessment of AW3D-30 and SRTM GDEMs
The accuracy of height from AW3D-30 and SRTM GDEMs were evaluated using mean sea level height (H msl ) at the 38 Benchmark around Perlis and the result is shown in Table 3. In general, there is small difference between SRTM and AW3D-30 model. Therefore, both models can be said as in good agreement with the H msl with AW3D-30 model has slightly strong correlation value of 0.9957. Based on this statistical analysis, AW3D-30 model showed more accurate than SRTM where the ME value between the H msl and the GDEM was 2.748m. The estimated accuracy for the AW3D-30 is ±3.048m. Meanwhile, the ME and RMSE of SRTM were 3.126m and 3.6521m, respectively, as presented in Table 3.

Evaluation of the geoid model
In this study, a gravimetric geoid model for northern region of Peninsular Malaysia has been computed using KTH approach in which go_cons_gcf_2_dir_r3 GGM model and AW3D-30 GDEM model were utilized. To verify the accuracy of the geoid obtained from this study, the model (N NGM17 ) was being validated with geometrical geoid (N geo ) derived from ellipsoidal height (H gnss ) using 38 GNSS observation on benchmark and mean sea level height (MSL) at the benchmark as follows: In addition, the accuracy of the resulted geoid (N NGM17 ) was also compared with the local gravimetric geoid available from DSMM (derived previously by R-C-R method with Gravsoft software) and both models are shown in Figure 8(a) and (b) respectively. The geoid accuracy in this case is evaluated based on ME and RMSE using the following equation: Both models show the geoid height at the northern part is higher compared to the southern region . Figure 6 shows a plot of absolute orthometric height differences of the NGM17 and the DSMM geoid model relative to the msl. The statistical analysis of NGM17 and DSMM geoid model is shown in more consistent compared to NGM17 as shown in Table 4 with minimum and maximum difference of 1.0310m and 1.2430m respectively. However, NGM17 geoid model shows better accuracy with the mean error of NGM17 and DSMM geoid model were 0.2568m and 1.1648m, respectively while the RMSE is ±0.2686m and ±1.1656m respectively.  more consistent compared to NGM17 as shown in Table 4 with minimum and maximum difference of 1.0310m and 1.2430m respectively. However, NGM17 geoid model shows better accuracy with the mean error of NGM17 and DSMM geoid model were 0.2568m and 1.1648m, respectively while the RMSE is ±0.2686m and ±1.1656m respectively.

Conclusion
The main aim of this paper is to present the modeling of the geoid model (NGM17) over northwestern region of Peninsular Malaysia (Kedah and Perlis) using KTH method by combining terrestrial gravity data, AW3D30 GDEM, DTU 10 Global Gravity Field and WGM2012 Earth's gravity anomalies. Based on the statistical figure analysis done shows that the geoid model obtained in this study (i.e., using the KTH method) giving slightly a better result compared to the Peninsular Malaysia's geoid model computed using RCR method with RMSE of ±0.2686m. The results from this present study more consistent compared to NGM17 as shown in Table 4 with minimum and maximum difference of 1.0310m and 1.2430m respectively. However, NGM17 geoid model shows better accuracy with the mean error of NGM17 and DSMM geoid model were 0.2568m and 1.1648m, respectively while the RMSE is ±0.2686m and ±1.1656m respectively.

Conclusion
The main aim of this paper is to present the modeling of the geoid model (NGM17) over northwestern region of Peninsular Malaysia (Kedah and Perlis) using KTH method by combining terrestrial gravity data, AW3D30 GDEM, DTU 10 Global Gravity Field and WGM2012 Earth's gravity anomalies. Based on the statistical figure analysis done shows that the geoid model obtained in this study (i.e., using the KTH method) giving slightly a better result compared to the Peninsular Malaysia's geoid model computed using RCR method with RMSE of ±0.2686m. The results from this present study