MAETEC: AN OPTIMIZED MATLAB-BASED APPROACH FOR ACCURATE AND EFFICIENT TERRAIN CORRECTION IN COMPLETE BOUGUER ANOMALY

2.1. Data and Tools

This study uses the global Digital Elevation Model (DEM) data from SRTM (Shuttle Radar Topography Mission) USGS Earth Explorer with a resolution of 1 arc second or about 30 meters (OpenTopography, 2013). The study case area is around Telomoyo Mountain, part of the geothermal exploration area in the middle of Java Island, Indonesia (Abdillah et al., 2024), where the DEM covers about 120 by 120 km, and the synthetic gravity station is located in the middle of the study area (see Figure 2). According to the geological map of the Magelang and Semarang (Thanden et al., 1996), the Telomoyo area's lithological unit mainly consists of Quaternary volcanic breccia and lava, surrounded by Tertiary sedimentary rocks overlaid by Quaternary volcaniclastic products. According to Saibi et al. (2012), the average rock density value in the Ungaran area is 2470 kg/m³ or 2.47 gram/cm³, and this value is used as the terrain density (ρ) input for terrain correction in the Telomoyo area due to their similar geological settings and volcanic history.

Figure 2. Available digital elevation model (DEM) data covers an area of 18,200 km² acquired from SRTM USGS Earth Explorer with a square gridded resolution of 50 meters on the Central of Java Island, Indonesia, which results in 7,280,000 pixels (red rectangle). The inset map shows the location of the 861 synthetic gravity station (black circle) for terrain correction around Mt. Telomoyo geothermal exploration area that covers an area of 200 km² (black rectangle).

The code is written in MATLAB software, and all the computational processes performed using 11^th Gen Intel® Core (TM) i7-1160G7@1.20GHz (8 CPUs) running on a 64-bit Windows 11 operating system with a memory of 16 GB of RAM. Using these computer specs, the code can run up to about 7.3 million points of DEM data, and it takes about 1-5 minutes to calculate the terrain correction of 861 gravity stations.

This research uses GTeC (Gravity Terrain Correction), the same MATLAB-based software, for benchmarking because GTeC provides a reliable and flexible solution for accurately calculating terrain corrections. GTeC's ability to integrate multiple DEM with progressively finer resolutions improves accuracy, especially in areas with steep slopes or near coastlines. Additionally, its use of tessellation-based algorithms for inner zones allows precise modelling of terrain effects, even in complex and rugged regions. GTeC software has been validated using synthetic topographic data and applied to a microgravity survey real-case (Cella, 2015). The validation results demonstrated excellent agreement of GTeC results with those obtained from commercial software, showcasing its reliability and accuracy.

2.2. Input Data and Setting Parameter

The input data is DEM data with a horizontal resolution of 50 meters for detailed geology prospecting depending on computer specs where the finer resolution will result in better accuracy. Still, it demands more memory and computational time. The DEM input data file should be in the *.xyz format consisting of coordinates X and Y (in meters), UTM Zone (e.g. WGS 1984) format, and the Z-elevation unit (in meters). This code also accepts bathymetric DEM data (if available) in the negative elevation values. We only use positive values in the Mt. Telomoyo terrain correction test case because TCDEMF does not accept bathymetric data. However, we include negative DEM values (bathymetric) for the regional terrain correction test in the Atlas Mountain regional case. In general, DEM for terrain correction should cover a minimum of 50 km in all directions from the terrain gravity study area (see Figure 2). The data input format is in *.txt format, with the information of gravity station name, longitude (can be in UTM or decimal degrees), latitude (can be in UTM or decimal degrees), and elevation of the gravity stations (in meters).

The parameter setting in MAETEC needs to be set as the input, including the average terrain density values of the study area (in gram/cm³). The maximum terrain correction radius search (maxTCrange) should be defined in meters. An input of 50,000 meters for a maximum terrain correction radius should be sufficient; beyond this range, the effect will not be significant to the measured point. For the zonation, MAETEC divides zonation into two main zones: the near and far zones. Each zone consists of two sub-zones, Zones 1 and 2 in the near zone and Zones 3 and 4 in the far zone (see Figure 3).

To control the radius of each zone in the setting, we need to determine the DEM resolution (dx) as the input (unit in meters). Zone 1 is the closest area when the gravity station is not located exactly on the DEM point, and Zone 1 is fixed as four nodes closest to the gravity stations (see Figure 3). Zone 2 (zone2rad) and Zone 3 (zone3rad) can be set based on resolution and the number of desired grids as the input where the far zone (Zone 3 and 4) DEM resolution is resampled to a lower resolution as input resampling to save the computer memory. For example, if the dx input is 50 meters and the resampling input is 4 then the new DEM resolution for Zone 3 and 4 will be 200 meters (see Figure 3). The maximum radius of Zone 4 is taken from the parameter of maxTCrange. The variation of grid size or resolution of the DEM and the zonation radius significantly impacts the accuracy of the terrain correction and computational time.

Figure 3. The left figure shows DEM zonation determination for calculating the terrain effect in each gravity station. DEM is divided into two main zones for each observation point: the near zone (red grid inside the red circle) and the far zone (yellow grid outside of the red circle). The near zone consists of Zones 1 (DEM point in black square) and 2 (DEM point in red square). The far zone consists of Zone 3, located between the red and blue circle (DEM point in blue square), and Zone 4, located outside of the blue circle, up to maxTCrange distance (DEM point in yellow square). The right figure illustrates how to calculate terrain correction in a gravity station for each DEM point (red rectangular prism with eight corner) in Zone 2 and Zone 3, which applies the flat top squared prism.

2.3. Algorithm Implementation

2.3.2. Far-zone (zon2rad < R_obs < maxTCrange)

For the far zone, the DEM needs to be resampled (by resampling parameter) with a larger grid size (see Figure 3) to reduce the computational time. However, it is important to note that decreasing the DEM resolution (increasing grid size or cross-sectional area) may lead to lower terrain correction accuracy. Many terrain correction programs, such as FA2BOUG, GTeC, and Oasis Montaj, use this resampling strategy for different zone radius (Fullea et al., 2008; Cella, 2015). In Zone 3, MAETEC implements the same flat-top square prism as Zone 2 (see Equation 2) but with a larger grid size or resolution to calculate the terrain effect ([TeX:] g_{zone3}$g_{zone3}$).

Since the flat-top squared prism can be simplified to a simple equation for a large distance (Blais and Ferland, 1984), we employ Ferland’s equation to calculate terrain correction ([TeX:] g_{zone4}$g_{zone4}$) for Zone 4:

[TeX:] g_{zone4} = G\rho A(\frac{1}{R} - \frac{1}{\sqrt{R^{2} + H^{2}}})$g_{zone4}=G\rho A(\frac{1}{R}-\frac{1}{\sqrt{R^2+H^2}})$ (3)

A represents the cross-sectional area, grid size, or resolution of the resampled DEM. R and H are the same as the distance and DEM elevation in Equation 1. Setting the minimum distance for Zone 4 is important to optimize the algorithm's computational time since the computational cost for the flat-top square prism is more expensive than Equation 3. The minimum distance for Zone 4 (zone3rad) is suggested at 10000 meters regardless of the height of the prism for the DEM resolution less than 2000 meters or less than 4[TeX:] \times$\times$10⁶ square meters cross-sectional area (A). However, increasing the DEM resolution or cross-sectional area (A) requires further radius or distance (R) to reach a high accuracy (lower error), where the best distance for Zone 4 (zone3rad) is selected when less than 10 microgal accuracies are achieved by considering the elevation difference (H) and DEM resolution or cross-sectional area (A) (see Figure 4).

Figure 4. A plot of the difference gravity effect of a single DEM point calculated using the simplified squared prism (Equation 3) and Flat-top squared prism (Equation 2) considering the distance of the gravity station with the DEM point (R), the elevation of the DEM point (H), and the DEM resolution or the cross-sectional area (A) (R, H, and A referring to Equation 3). The red line is distance (R) = 10 km.

3.1. MAETEC optimum setting

Mount Telomoyo geothermal exploration area is selected to test the accuracy and optimum setting of the MAETEC software. 861 synthetic gravity stations were calculated using 50-meter grid size resolution DEM from SRTM (OpenTopography, 2013) (see Figure 2). Utilizing the same DEM resolution input, GTeC is considered the reference for other terrain correction software because it enables triangular tessellation to mesh the DEM. The triangular tessellation to mesh calculates the terrain effect using sloped triangular prisms, ensuring greater precision (Götze and Lahmeyer, 1988). Theoretically, triangular tessellation is closer to the real topographic condition, where it reduces the “staircase effect” due to the flat top squared prism (Cella, 2015). The GTeC settings to get good accuracy for radius zones 1, 2, 3, 4, and 5 are 100 meters, 1000 meters, 1600 meters, 14000 meters, and 50000 meters with the step or grid resolution of 10 meters, 50 meters, 100 meters, 200 meters, and 400 meters, respectively (see Table 1).

To optimize the settings for MAETEC and ensure computational accuracy, we tested four zonation radius and grid size scenarios (see Table 1). Scenarios 1, 2, and 3 progressively decrease terrain correction accuracy by altering zonation parameters and DEM resolution, whereas Scenario 4 aims to balance accuracy with computational time. We assessed accuracy by comparing the terrain corrections from each scenario with the reference GTeC at each gravity station location and computing its root mean square error (RMSE).

In Scenario 1, Zone 2 (zon2rad) is set to 1,600 meters. The DEM resampling for the far zones (Zones 3 and 4) is performed at four times the original grid spacing, resulting in a 200-meter grid resolution. The radius of Zone 3 (zon3rad) is set at 14,000 meters to enhance accuracy, and the radius of Zone 4 is set at 50,000 meters (see Table 1).

In Scenarios 2 and 3, Zone 2 (zon2rad) is reduced to 800 meters to evaluate the improvement in computational time and its effect on near-zone (Zones 1 and 2) accuracy. Additionally, the radius of Zone 3 (zon3rad) is decreased to 10,000 meters in both scenarios to further improve computational efficiency (see Table 1). The primary difference between Scenarios 2 and 3 lies in the DEM resampling of the far zones: Scenario 2 uses a resampling factor of four times the original resolution (4*dx), while Scenario 3 uses a factor of eight times (8*dx) to assess the impact on accuracy.

Finally, Scenario 4 maintains the near zone (Zones 1 and 2) radius at 1,600 meters (zon2rad), as in Scenario 1. However, the DEM resampling for the far zones (Zones 3 and 4) is increased to eight times the DEM resolution (dx), resulting in a 400-meter grid resolution (see Table 1). In Scenario 4, the radius of Zone 3 (zon3rad) is set to 10,000 meters to reduce computational time while maintaining accuracy since the far zones (particularly Zone 4) usually have less influence on terrain gravity.

Table 1. The comparison of the computational time average and root mean square error (RMSE) of different available terrain correction software with different scenarios, including zonation ranges, terrain correction methods, and DEM resolution. GTeC, TCDEMF, and Oasis Montaj zone radius and DEM grid size have been adjusted to the MAETEC zonation setting for simplification.

3.2. Analysis of MAETEC Scenarios

The result of terrain correction from scenario 1, where the larger Zone 2 radius of 1600 meters (zon2rad) results in better accuracy in the near zone (Zone 1 and 2) while minimizing the resampling factor at 4 times, leading to better accuracy in the far zone (Zone 3 and 4). Determining the radius of Zone-3 (zon3rad) at 14000 meters will also increase the accuracy in the far zone (Zone 3 and 4), where the accuracy test of Eq. 3 in Zone 4 suggests that the distance is more than 14000 meters at a cross-sectional area of less than 10⁶ square meters (or grid resolution of 1000 meters) at any height will have the error less than 1 microgal per DEM point (see Figure 4). Therefore, this scenario results in the highest where the root means square error (RMSE) of terrain correction compared to the GTeC as the reference is only 65 µgal (see Table 1). However, scenario 1 spent around 330 seconds, the longest time compared to other MAETEC scenarios.

The MAETEC algorithm shows that the near zone is the most important zone for acquiring good accuracy, which can be seen by comparing scenarios from the radius of Zone 2 of scenarios 1 and 2 with scenarios 3 and 4. Thus, the computational time is spent mostly in Zone 2, where the flat-top squared prism is implemented. As the flat-top squared prism is also applied in Zone 3, the computational time will increase if the radius of Zone 3 (zon3rad) gets larger. However, accuracy will also increase as the Zone 3 radius increases, as shown by the RMSE and computational time in scenarios 3 and 4. Thus, Scenario 4 is the optimum setting for MAETEC, which gives accuracy similar to scenarios 1 and 2 but faster (2-5 times faster) due to the DEM resampling in the far zone (Zone 3 and 4) (see Table 1).

3.3. Comparison with Available Software

The accuracy of the MAETEC needs to be verified by other terrain correction software. This study compares three software with different methods and strategies to evaluate MAETEC accuracy and computational time. The first software is TCDEMF (Nishijima, 2009), which separates zonation into two zones and implements a single slope octant for the closest zone (Kane, 1962) and a flat-top squared prism for the rest of the DEM data (see Table 1). The second software is GTeC (Cella, 2015), which separates the zone into five zones that implement flat-top squared prism in Zones 3, 4, and 5 and triangular tessellation in Zones 1 and 2 (see Table 1). To optimize the GTeC accuracy, the zonation distance setting, different input DEM resolutions (steps), selection of terrain correction method in zones 4 and 5 (near zone), and parallel computing are available in the input setting. The last software is the commercial Geosoft Oasis Montaj, which implements algorithm methods described by Kane (1962) and Nagy (Nagy, 1966; Plouff, 1976; Nagy et al., 2000), which apply advanced grid-mesh interpolation, zonation, and resampling techniques (see Table 1).

Figure 5. Maps of the terrain correction calculated at the points stations are shown in Figure 2 using different software: GTeC (a), MAETEC (b), TCDEMF (c), and Oasis Montaj (d). The error of each software is compared with GTeC software as a reference. The MAETEC error is relatively low, with a maximum error of less than 0.5 mgal (e). The TCDEMF error is up to ± 2.5 mgal at high terrain areas (f). The Oasis Montaj error appears highest compared to other software, up to + 3 mgal (g). The plot showing the accuracy comparison between MAETEC and other terrain correction software, using GTeC as the reference. All software demonstrates high correlation coefficients above 0.980, indicating comparable accuracy, with MAETEC achieving the highest correlation of up to 0.999(h).

In the GTeC, the terrain correction value ranges from 1.049 mgal to 24.032 mgal where the maximum value terrain effect is at the peak of Telomoyo Mountain (see Figure 5a). MAETEC uses scenario 4 to compare with other software by plotting and visualizing the result, which results in a range of terrain correction from 1.056 mgal to 24.052 mgal (see Figure 5b). As shown in Table 1, the RMSE of MAETEC with the reference software (GTeC) is less than 0.1 mgal, but some areas still have a difference of about 0.2 mgal (see Figure 5e). TCDEMF ranges from 1.086 mgal to 24.721 mgal (see Figure 5c), which differs with reference software at some locations from -2.592 mgal to 0.693 mgal (see Figure 5f). Oasis Montaj terrain correction result ranges from 0.247 mgal to 22.883 mgal (see Figure 5d); this results in about 0.483 mgal to 2.773 mgal lower value of terrain correction than the reference software (see Figure 5g).

Overall, MAETEC can produce an accurate terrain calculation with an average coefficient correlation of 0.999 with GTeC, 0.987 with TCDEMF, and about 0.985 with Oasis Montaj software (see Figure 5h), despite using a different method and algorithm than other software (see Table 1). The significant difference in the calculation of MAETEC with TCDEMF might be due to the resampling or gridding system of the DEM data during the conversion of the coordinate from decimal degree to UTM, whereas in the MAETEC and GTeC input can be directly in meters coordinate and the TCDEMF input is in decimal degrees. The significant difference in the Oasis Montaj is due to the resampling and zonation technique where, according to the explanation of the software to improve processing efficiency, the far zones are sequentially resampled by a factor of two then set to the average of the grid values they cover (see Table 1). Zone 2 averages every four initial neighboring DEM cells, while Zone 3 averages every sixteen initial cells and applies to the next zones.

3.4. Regional Terrain Correction

MAETEC is also tested for regional data with topography or DEM data resolution of the 1’ x 1’ ETOPO Global database, gridded at about 2000 meters resolution, taken from supplementary data in Cella (2015). The study case uses the same area as the one investigated by Fullea et al. (2008) and Cella (2015) from the Atlas Mountain, Morocco (see Figure 6). This terrain correction can be used for regional studies such as crustal thickness or tectonic studies. As mentioned by Cella (2015), the error of gravity data by FA2BOUG software for the regional data is up to ± 20 mgal in the higher slope gradient, which is very big considering the gravity anomaly in the subsurface. Therefore, it is important to create the best scenarios to ensure the accuracy of the computational cost of the regional terrain correction by MAETEC (see Table 2).

Table 2. The comparison of input parameters of MAETEC with GTeC for the regional terrain correction test, the computational time, and the accuracy (RMSE). GTeC is the reference for calculating RMSE.

To evaluate the regional gravity correction in Atlas Mountain, GTeC parameter was set for a density of 2.67 g/cm³ (Cella, 2015) with search radius zones 1, 2, 3, 4, and 5, which are 2000 m, 6000 m, 38000 m, 86000 m, and 167000 m with the grid size of 250 m, 500 m, 1000 m, 2000 m, and 4000 m, respectively (see Table 2). With this parameter, the GTeC produces gravity values ranging from 1.137 mgal to 164.194 mgal (see Figure 7a), where most of the contributors to gravity terrain correction come from Zone 1 and 2, which apply a sloped triangle prism method.

For the MAETEC setting, the grid size of the near zone (Zones 1 and 2) is resampled to 250 meters (dx), with the radius of Zone 2 being 4000 meters (zon2rad), to which the MAETEC applies a linear gridding algorithm. The far zone (Zone 3 and 4) grid size is resampled 8 (resampling) times of the near zone, corresponding to a 2000-meter grid size with the zon3rad and maxTCrange set to 20000 meters and 167000 meters as the radius for the Zone 3 and 4, respectively (see Table 2). Compared to GTeC, the MAETEC terrain calculation produces an RMSE of 1.009 mgal and needs only 80 seconds to compute all 4716 survey points (see Table 2). MAETEC code produces gravity terrain correction values ranging from 0.787 mgal to 165.039 mgal (see Figure 7b). The difference of this terrain correction with GTeC code is considered low with a maximum error of not more than ± 3.15 milligal (see Figure 7c).

Figure 6. Regional case terrain correction in North Western Africa with DEM resolution of 2000 meters, data taken from ETOPO Global Database (NOAA National Geophysical Data Center, 2009). A total of 4717 synthetic gravity stations are located onshore (Atlas Mountain) and offshore (at sea level), covering an area of 18200 km² (inside red rectangle).

MAETEC code can calculate the terrain correction accurately and efficiently for regional gravity terrain correction purposes. The plot of GTeC versus MAETEC correlation coefficient is 0.999, which is considered a very good result (see Figure 7d). The high error from the plot is between 0 and 20 mgal value, corresponding to the shallow marine area (see Figure 7c). This high error is affected by the topography of the DEM in the near zone due to the continental slope of the African continent.

In the regional test, DEM resampling must be performed to produce a smoother topography that leads to more accurate terrain correction. Resampling becomes crucial for regional terrain gravity correction when the input DEM resolution exceeds 100 meters. According to Cella (2015), down-sampling DEM resolution will increase the accuracy of the terrain correction. Thus, it is recommended to perform a down-sampling of the original DEM to 1/8 of its original size for large DEM cases (from 2000 meters to 250 meters) to balance accuracy and computational time. Meanwhile, the near zone radius (zon2rad) suggested a maximum of 32 times the resolution down-sampled DEM, which is 4000 meters (see Table 2). For the far zone, zon3rad was set at 20000 meters to get an accuracy of less than 0.1 microgal for the 2000-meter grid size DEM input (4[TeX:] \times$\times$10⁶ m² cross-sectional area) (see Figure 4).

Figure 7. Maps of the regional terrain correction calculated at the points stations that are shown in Figure 6 of GTeC (a), MAETEC (b), and the difference (error) of the MAETEC compared with GTeC as the reference software (c). The plot for each gravity station regional terrain correction of GTeC and MAETEC shows a very good correlation (d).

4.1. Computational Efficiency

One of the most notable achievements of MAETEC is its ability to significantly reduce computation times without compromising accuracy. This is primarily due to the utilization of MATLAB's parallel computing capabilities and an optimized zonation strategy. The parallel computing toolbox allows the distribution of computational tasks across multiple CPU cores, effectively speeding up the processing of large DEM datasets.

The comparison of different zonation scenarios (see Table 1) indicates that MAETEC's Scenario 4, which combines a larger near-zone radius with effective DEM resampling in the far zone, offers the best balance between accuracy and computational time. This scenario achieves an RMSE of 65 µgal compared to GTeC, the reference terrain correction software, while completing the calculations five times faster than Scenario 1. This setting strategy can also be applied to the down-sampling case of the low-resolution DEM input for a regional case study, which results in an RMSE of only 1 mgal for the study case but five times faster than the GTeC code.

MAETEC also shows better efficiency than Oasis, TCDEMF, and GTeC by choosing the right terrain correction equation for each zonation. In this case, the most important part is that implementing the simplified flat-top prism in Equation 3 (Blais and Ferland, 1984) on the right distance radius of the far zone will save much computational time yet maintain high accuracy compared to resampling the DEM grid size and calculating all data by the flat-top prism which is applied in many terrain correction software.

4.2. Accuracy and Validation

The accuracy of MAETEC was validated against other terrain correction software, such as GTeC, TCDEMF, and Oasis Montaj. The results demonstrate a high correlation between MAETEC and GTeC, with a coefficient of 0.999. Minor discrepancies up to 0.2 mgal are within acceptable limits for most geological applications. For regional practice, MAETEC also shows great accuracy with a maximum error of around 3 mgal compared to software such as FA2BOUG, which has errors up to 20 mgal in the steep slope area, referring to Cella (2015).

The differences between TCDEMF and Oasis Montaj are attributed to their distinct resampling and zonation techniques. TCDEMF shows a variation due to the conversion of coordinates and potential interpolation errors during DEM gridding. Oasis Montaj, on the other hand, employs more aggressive resampling strategies to enhance processing efficiency, which can lead to underestimation of value in some areas. Meanwhile, the higher FA2BOUG error (Cella, 2015) might be caused by using only two zonation techniques without applying a down-sampling technique, which produces big errors in the higher slope gradient.

4.3. Practical Applications

MAETEC's ability to handle high-resolution DEM data efficiently makes it particularly useful for detailed geological and geophysical surveys, such as geothermal exploration with higher terrain complexity, as demonstrated in the Telomoyo Mountain case study. The algorithm's flexibility in setting zonation parameters and DEM resolutions allows users to adapt the processing according to their specific requirements and available computational resources. For regional studies, such as the Moroccan region case, MAETEC proved to be capable of processing lower-resolution DEM data while maintaining acceptable accuracy levels.

The input process for MAETEC is straightforward and similar to the settings used in Oasis Montaj, FA2BOUG, and TCDEMF. It requires a DEM input that uses the same UTM reference as the coordinate location of the gravity survey. Additionally, MAETEC can accept gravity survey data in latitude and longitude (decimal degrees) and provide outputs in UTM coordinates (Palacios, 2006). MAETEC can also generate UTM references for the gravity stations, which can then be used to determine the appropriate UTM reference for exporting the input DEM topography data. The zoning strategy is easy to use, with simple zonation in both near and far zones. As the input setting, it only needs the DEM size resolution (dx) for the near zone and the DEM resampling for the far zone and simple radius determination techniques, that are zon2rad, zon3rad, and maxTCrange for Zone 2, 3, and 4, respectively. This simplicity suggests that MAETEC is effective for both local and regional scale gravity studies and also easy to use, making it a versatile geophysical tool.