A Proposed Merging Methods of Digital Elevation Model Based on Artificial Neural Network and Interpolation Techniques for Improved Accuracy

Accurate terrain elevation is crucial for numerous geoscientific engineering applications. digital elevation model (DEM) is one of the most important sources of elevation (Hawker, Neal, and Bates 2019). In geostatistical analysis, DEM is used to quantitatively represent the bare Earth’s surface, which is used to describe the spatial patterns of various surfaces, e.g., surface water, ground surface, and canopy (Mukherjee et al. 2013). This can be represented in a raster or vector model in geographical information system (GIS) (Habib et al. 2020, Huang and Yang 2011). Aside the DEM, digital surface model (DSM) and digital terrain model (DTM) are other elevation surfaces that are frequently used for the ground terrain (Kim et al. 2020). However, DTM represents the Earth’s terrain, i.e., bare ground (Xiong, Wang, and Wessel 2017), while DSM includes objects on the ground, such as buildings, vegetation, and human-made structures (Nwilo et al. 2022). Apparently, most data providers and professional users use the term “digital elevation model” as a general term for digital topography (Badura and Przybylski 2005), (Varga and Bašić 2015). It is important to note that there are four different representations of DEM. They are regular grid, outlines, triangulated grid, and triangulated irregular network (TIN) (Khemiri et al. 2013, Leitão, Prodanović, and Maksimović 2016). Generally, DEM can be created using remote sensing (Habib et al. 2020), field surveying (Mesa-Mingorance and Ariza-López 2020), aerial photography (Arun 2013), satellite images (Yap et al. 2019), interferometry, and laser surveying. This can also be achieved through scanning and digitizing topographic maps (Khemiri et al. 2013) and field surveying, e.g., the leveling, total station, and global navigation satellite system (GNSS) survey (Nwilo et al. 2022).

Recently, Light Detection and Ranging (LIDAR), Radio Detection and Ranging (RADAR) altimetry (Yap et al. 2019), and Interferometric Synthetic Aperture Radar (InSAR) (Zhou et al. 2012) have been the most advanced and widely used technologies for acquiring three-dimensional information over large areas with high horizontal resolution and high vertical accuracy (Jain et al. 2007). These techniques are expensive, i.e., high instrument and survey costs (Shaikh, Yadav and Manekar 2021, Capolupo 2021). Thus, they are not well integrated in developing countries (Jakovljevic et al. 2019, Kulp and Strauss 2018). Consequently, most users in the developing countries resort to freely accessible open-source products.

There are several global DEMs (GDEMs), such as the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) and the Shuttle Radar Topography Mission (SRTM), created for free and commercial use (Reuter, Strobl and Mehl 2011). These GDEMs represent the new-generation models, which have been significantly improved in comparison with the first released models like the global 30 arc-second elevation (GTOPO30) (Mukherjee et al. 2013) and one arc-minute global relief model (ETOPO1) (Varga and Bašić 2015).

In this study, the SRTM DEM (1 arc-second) is considered in the investigation. It is a joint international project of the National Aeronautics and Space Administration (NASA), the National Geospatial-Intelligence Agency, and the German and Italian Space Agencies to collect three-dimensional digital mapping (Yang, Meng and Zhang 2011). The SRTM DEM is mostly used as the essential alternative data when there is no DEM availability, or there is budget constraint to purchase costly high-resolution satellite data (Wendi et al. 2016). Moreover, the GDEMs are vital sources of topographic information for many geoscience studies and applications (Mirza, Dawod and Al-Ghamdi 2011). Some examples include, the geodetic survey (Yamazaki et al. 2017), geological survey, water resources and hydrology (Hawker et al. 2019), evaluation of natural hazards and vegetation survey (Wendi et al. 2016), and mineral and petroleum industries (Badura and Przybylski 2005). In addition, they are applied in engineering works, e.g., highway construction and wind turbine location optimization (Akturk and Altunel 2019); land surface modeling, such as management of flood risk (Jakovljevic et al. 2019), volcanology, ecology, and glaciology (Kim et al. 2020, Mukherjee et al. 2013). Despite these achievements, their accuracy, still, does not meet the requirements for precise applications due to DEM errors (Yamazaki et al. 2017). Therefore, the accuracy assessment of DEMs is critical in these applications (Varga and Bašić 2015).

DEM error sources include data collection error (gross error), lacking orientation of stereo-images (systematic error), data processing methodology (Yue et al. 2007), georeferencing of stereo-images, grid spacing, interpolation techniques for DEM generation, characteristics of the terrain surface (Chen and Yue 2010), and random error, which cannot be avoided (Varga and Bašić 2015, Mukherjee et al. 2013). As a result, several investigations on DEM have been conducted by researchers leading to many methodical propositions to reduce errors and increase accuracy. For example, Yue et al. (2007), Cook et al. (2012), Wilson (2012), Sari et al. (2019), Hugenholtz et al. (2013), Chen and Yue (2010), Warriner (2005), and Zhou et al. (2012) suggested different methods for generating improved DEM. Further, DEM accuracy was evaluated by Hawker et al. (2019) over Croatia. Varga and Bašić (2015) investigated floodplain sites, while Akturk and Altunel (2019) considered a study area that was highly broken and vegetated terrain. Furthermore, other studies evaluated the vertical accuracy of open-source DEM, i.e., SRTM and ASTER (Yap et al. (2019), Miliaresis and Paraschou (2005), Suwandana et al. (2012), Mirza et al. (2011), Mukherjee et al. (2013), Falorni et al. (2005), and Shaikh et al. (2021)). To improve the accuracy of the DEM, Yamazaki et al. (2017) produced a high-accurate map of global terrain elevations, Capolupo (2021) enhanced the accuracy of GDEM of differences in the Google Earth Engine for 3D change detection analysis, Muhadi, Mohd Kassim, and Abdullah (2019) Improved DEM using data fusion technique for oil palm replanting phase. In addition, Xu and Zhou (2016) suggested recovering distorted DEMs of regular terrain while Deilami et al. (2012) recommended reducing pit and flat error for improving DEM accuracy.

In constructing DEMs, interpolation methods make use of an effective factor. Consequently, many investigations have been achieved for selecting the suitable method for enhancing DEM accuracy. Such as, Habib et al. (2020) studied the impact of interpolation techniques on the accuracy of the large-scale DEM. Rishikeshan, Katiyar and Mahesh (2014) also evaluated DEM interpolation methods in GIS. Moreover, Jana (2011), Arun (2013), and Schwendel, Fuller, and Death (2012) established the appropriate interpolation methods for creating DEMs.

In the past few years, artificial neural networks (ANNs) have received increased attention for achieving unprecedented results in many complicated geomatics problems. ANN has been applied in the creation of DEMs, leading to some contributions like the SRTM DEM improvement (Kim et al. 2020, Wendi et al. 2016, Kulp and Strauss 2018), landslide detection (Gorsevski et al. 2016, Kawabata and Bandibas 2009), calibration and verification of DEM parameters (Ye et al. 2019, Benvenuti, Kloss and Pirker 2016), DTM creation and flood risk mapping (Jakovljevic et al. 2019), and the analysis of remotely sensed data for producing DEM (Mas and Flores 2008).

Recently, an effective merging technique was applied to the elevation data sources. This has been suggested due to the availability of elevation data sets for more than one area. In some cases, these datasets have low accuracy in some areas and high accuracy levels in other areas. Therefore, using a low-accuracy dataset for the whole area when parts of the area have high-accuracy data is not recommended. Consequently, merging the most accurate data sources (e.g., DEMs) to produce a single DEM that covers the whole area of interest with the highest possible resolution and accuracy was proposed. For example, Gruber et al. (2013), Choussiafis, Karathanassi and Nikolakopoulos (2012), and Jain et al. (2007) applied mosaic methods to improve DEM accuracy.

In this study, we hypothesize that, if there is available data with high accuracy (i.e., GNSS observations or high-accuracy DEM) for a small area S located in a large area B, a highly accurate DEM for the large area B can be created. This can be achieved by proposing DEM merging methods based on ANN and interpolation techniques. These techniques are employed for generating a correction surface (Vs) for the large area. The SRTM DEM is adjusted by adding the correction surface and merging with the small area S. Note that the small area S is the overlap area between the SRTM DEM and the large area B. Therefore, the main objectives of this study are to (i) generate a highly accurate DEM for a large area based on GNSS observations for a small portion of that area, (ii) mutually evaluate the proposed DEM merging methods (PDMMs) in comparison with the conventional DEM merging methods (CDMMs) available in ArcGIS software (in terms of elevations, slopes, and aspects), and (iii) compare and analyze the accuracy of the merged DEMs with SRTM DEM.

Interpolation is the process of predicting the values of unsampled locations based on collected samples, i.e., neighboring measurements (Sacchi, Ulrych and Walker 1998). The GIS software presents powerful analytical tools for various interpolation techniques (Habib et al. 2020). Some examples include kriging, inverse distance weighted (IDW), natural neighbor, triangulation with linear interpolation (TLI), spline, spline with barriers, topo to raster, and trend. However, the most widely used interpolation methods are Kriging, IDW, and spline (Yang et al. 2004). The IDW and spline are classified as deterministic interpolation methods. They are directly based on the surrounding measurements or specified mathematical formulas that determine the smoothness of the resulting surface. Kriging is classified as a geostatistical method (i.e., based on statistical models). This includes the statistical relationships among the measured points. Geostatistical techniques create a prediction surface and evaluate it at appreciable accuracy (Heine 1986). These methods are illustrated in the following sections.

Artificial neural network (ANN) is one type of machine learning algorithm that has, lately, received increased attention for achieving promising results in different geoscience applications (Marmanis et al. 2015, Mas and Flores 2008). This algorithm is implemented as a system of interconnected processing elements that are functionally analogous to biological neurons (Kim et al. 2020). ANN consists of three layers: the input, hidden, and output layers (Figure 1).

Figure 1.

As shown in Figure 1, the input layer includes input P, which has a single vector of R elements that relate to the hidden layer by synaptic weights, WNh,R′h$$W_{{N^h},R'}^h$$. The synaptic weights WNo,Nho$$W_{{N^o},{N^h}}^o$$ make connections between the hidden and output layers (Gurney 2018). The superscripts h and o represent the hidden and output layer, respectively. The first index indicates the specific neuron destination for that weight; the second index denotes the source of the signal fed to the neuron; while N^h, N⁰ represent the number of neurons (NON) in the hidden and output layers, respectively. These neurons represent the processing units in ANN, where each neuron has a bias b, summation S, activation function f, and output a. The hidden layer’s output a^h is the input for the output layer, which produces an output a⁰. The bias b is like a weight, but it has a constant input of 1. The summation S can be calculated in the hidden and output layers using Equations (4) and (5) (Hagan, Demuth and Beale 1996): 4Sjh=∑i=1RWjihPi+bjh$$S_j^h = \sum\nolimits_{i = 1}^R {W_{ji}^h{P_i} + b_j^h} $$ 5Sko=∑j=1NWkjoajh+bko$$S_k^o = \sum\nolimits_{j = 1}^N {W_{kj}^oa_j^h + b_k^o} $$ where Sjh$$S_j^h$$ and Sko$$S_k^o$$ represent the weighted inputs summed with the bias at the hidden neuron j and output neuron k, respectively.

This summation is fed to the activation function (i.e., transfer function), which is transformed to specific values depending on the type of activation function considered. There are three activation functions that are most commonly used (Han and Moraga 1995, Haykin 2009). The first function (Equation 6) is a linear transfer function (Purelin), which produces output equal to its input, where: 6f(S)=S$$f(S) = S$$

The second function (Equation 7) is the log sigmoid (Logsig); this transfer function takes the input, which lies between plus and minus infinity, and transforms the output into the range 0–1, using the expression: 7f(S)=11+e−S$$f(S) = {1 \over {1 + {e^{ - S}}}}$$

Equation (8) presents the third function, which is the hyperbolic tangent sigmoid (Tansig); this function is like the log sigmoid function, but the output here is transformed into the range −1 to 1 using the expression: 8f(S)=eS−e−SeS+e−Sf(S)=11+e−S$$f(S) = {{{e^S} - {e^{ - S}}} \over {{e^S} + {e^{ - S}}}}f(S) = {1 \over {1 + {e^{ - S}}}}$$

Although there are different types of ANNs, the feed-forward neural network, the most straightforward type of ANN, is frequently considered. It generates a relationship between the inputs and outputs using a backpropagation (BP) algorithm (Seiffert 2001, Schmidhuber 2015). The Levenberg–Marquardt (LM) algorithm, often the fastest BP algorithm, is applied for training datasets in the network to minimize errors (cost function) using the mean squared error (Levenberg 1944, Marquardt 1963). The training is continued until the stopping criteria meet the requirements. The dataset, in most cases, is divided into 70% for training, 15% for validation test, and 15% for independent testing of the network’s performance.

GIS software provides functions for merging two or more grid-based (raster) datasets. These methods include the first, the last, the minimum, the maximum, the mean, and the blend. They are based on the assumption that grid-based DEMs have the same spatial resolutions (cell size) as well as a coordinate system.

The First method determines the pixel value from the first raster dataset, i.e., high accuracy DEM, in the mosaic list. It does not operate any elevation adjustment on the DEMs; the DEMs are just superimposed. Then, the existing raster dataset is considered first. Figure 2 shows the results when there are four overlapping pixels, and the First option is selected. The values of the first raster dataset (on the left) are superior to another mosaicked dataset. Therefore, the result is the same as the first set of overlapping pixels.

Figure 2.

The Last method determines the pixel values from the overlapping last raster dataset. Figure 3 shows the result of a mosaic where the Last option is selected. When the two rasters are mosaicked, the overlapping values from the second raster dataset become the output mosaic.

Figure 3.

Figure 4 shows the Minimum method; it determines the lower pixel value from the two overlapping raster datasets. The illustration below shows the result of a mosaic where the Minimum option is selected. When the two rasters are mosaicked, the output mosaic results from the minimum values of the two raster datasets.

Figure 4.

The Maximum method determines the higher pixel values from the two overlapping raster datasets. Figure 5 shows the result of a mosaic where the Maximum option is selected. When the two rasters are mosaicked, the output mosaic results from the maximum values of the two raster datasets.

Figure 5.

The Mean method determines the average pixel values from the two overlapping raster datasets. As shown in Figure 6, the result of a mosaic resulting from the Mean option is shown; the output pixel type is float. When the two rasters are mosaicked, the output results from averaging the two overlapping values. If many raster datasets are overlapping, only two raster datasets are processed at one time. Also, if the output pixel type is integer, then the values are rounded.

Figure 6.

Figure 7 displays the Blend method; it uses a distance-weighted algorithm to determine the value of overlapping pixels. The output cell is a blend of the overlapping areas. This method is the most computationally intensive option for mosaicking.

In the illustration in Figure 7, the diagram shows two overlapping raster datasets. The cell where the x is located has two values: the value of the pixel in dataset R1 (outlined red) and the value of the pixel in dataset R2 (outlined blue). Since the x is closer to dataset R2, the value of the R2 pixel is more heavily weighted in the output. For more details, see ESRI (2020).

Figure 7.

As was previously stated, the resulting DEM by min., max., first, and last methods may generate terrain surface discontinuities (cliffs) along the boundaries. These discontinuities are created due to the elevation differences in (variation in accuracies) DEMs, leading to errors in the slope and aspect values (Leitão et al. 2016). On the contrary, Mean and Blend methods determine the average elevation value within the overlapping area of the two DEMs. Thus, a smooth transition between the merged DEMs is achieved. This results in new elevations.

The drawbacks identified in the conventional DEM merging methods reveal the need for new methods. These new methods can retain high-accuracy DEM data while creating smooth transitions between the two merged DEMs.

The study area is a region with an area equal to 3276400 m2 located in the Assiut governorate desert, Egypt. The elevations of the area’s topography vary between 88 m and 134 m. A reference station (R), fixed by the Egyptian survey authority, was selected for the GNSS survey. Figure 8 and Table 1 show the coordinates of the corner points of the study area and the reference point R.

Figure 8.

Two types of data were collected.

This study is based much on statistical comparisons. Thus, relevant statistical and engineering software and programs were employed for highly accurate results. These programs are: -

Trimble Business Centre (TBC),

The investigation was performed in five stages, which are shown in Figure 9 and discussed in the following sections.

Figure 9.

2.2.1.

Data processing and preparation

The GNSS observations were processed using Trimble Business Center (TBC) software based on the known station (R). The processed data was in the geodetic coordinate system (Latitude, Longitude, Ellipsoidal height); thus, the coordinates were converted to projected coordinates (Easting, Northing, Orthometric Height) according to the Universal Transverse Mercator (UTM) ZONE 36N.

The TBC outputs were used (in Autodesk Civil 3d and Excel programs) to produce grid data (5 m x 5 m). Also, the border coordinates (E, N) of the whole area with zeros and H elevations were assigned. Then, the overlap area coordinates (E, N, H) were extracted from the grid. These overlap coordinates were added to the border coordinates in two cases (Figure 10): the border with zeros elevations (zeros border) and the border with H elevations (H border). The grid files of the overlap and the whole area were converted to DEM files in the ArcGIS environment. For achieving a high level of accuracy in the statistical stage, especially in comparisons, the horizontal coordinates of the whole area grid (E, N) were used for extracting the elevations of the DEMs under test, which means the same horizontal coordinates were applied for all DEMs. These coordinates were extracted in an Excel file and converted to a shape file named Extraction File (E.F.) in ArcGIS. Accordingly, the E.F. was used to extract the elevations from the DEMSRTM, thereby converting DEMSRTM (raster data) to vector data. After that, the following files were retrieved: -

DEM file for the whole area (DEMGNSS), which was used as a reference surface in the statistics stage (Figure 11),

-

DEM file for the overlap area (DEM Overlap),

-

DEM file for the whole area from the USGS website, i.e., DEMSRTM (Figure 12),

-

Grid file for the whole area (GNSS processed data),

-

Grid file for the whole area (extracted from DEMSRTM data),

-

Grid file for the overlap area,

-

Grid files for the overlap area with the two cases: zeros and H borders.

Figure 10.

Figure 11.

Figure 12.

The DEM _Overlap and its location on the DEM_SRTM are shown in Figure 13 and Figure 14, respectively.

Figure 13.

Figure 14.

2.2.2.

Correction surface generation

This is the basic stage of the study. As mentioned in Section 1.1, a small area of the overlapping DEM_GNSS and DEM_SRTM was selected from the study area. The differences (corrections) between the DEM_SRTM and DEM _Overlap were calculated based on two techniques: (1) interpolation technique and (2) ANN technique. This was used as the correction surface, which was extended to the entire area.

2.2.2.1.

Interpolation techniques

The most popular interpolation methods, considered, were kriging, IDW, and spline. It is important to note that the external borders do not have data to interpolate. Thus, additional data is required at the borders of the whole area to possibly execute the interpolation operation. Consequently, the two scenarios mentioned in Section 2.2.1 were applied. The first scenario assumed that the differences between the elevations of DEM_SRTM and the border elevations are zeros (zeros borders) (Warriner 2005). The second scenario was proposed to test the effect of the same boundary with H elevations (H border). In this case, the differences between the elevations of DEM_SRTM and that of the borders were estimated. After that, the two scenario outputs were added to the overlap differences.

2.2.2.2.

ANN technique

Also, the ANN algorithm was formulated in MATLAB software for generating the correction surface (Figure 15 9). According to the ANN components, the horizontal coordinates of the overlap area (E_o, N_o) were used as inputs, while the elevation differences between the overlap area and DEMSRTM were the target of outputs. The horizontal coordinates of DEM_SRTM (E_D, N_D) were the sample data.

Figure 15.

It is important to mention that ANN initialization produces different results. Thus, the ANN initializations are often repeated several times to guarantee more accurate results. Consequently, the number of initializations was set at 100; this iterated automatically. One hidden layer with ten neurons was selected, according to the recommendation of Alemam, Yong, and Mohammed (2022). During the development of this ANN algorithm, the results indicated that the linear activation function (pure-line), in the hidden and output layers, is suitable for this network and data compared to the other activation functions. The Train L.M. was selected as the training function. The involved training parameters according to the Train L.M. function and learning rate are shown in Table 3. Unlike the interpolation method, the ANN does not need a zeros border. Apparently, two scenarios were applied; ANN only and ANN with H elevations (H border) boundary. This network produced the DEM_SRTM corrections (V_D) as the output. The average value of these outputs was estimated by applying the least squares method (LSM) to minimize the weighted residuals (v) by applying the normal equation: 9L+v=AX$$L + v = AX$$ and the solution model: 10X=(ATA)−1ATL$$X = {\left( {{A^T}A} \right)^{ - 1}}{A^T}L$$ in the output data, L is the output matrix, A is the design matrix, and X is the unknown matrix.

Table 3.

The parameters involved in the ANN algorithm

Epoch	Goal	Max_fail	Min_fail	Mu	Learning rate
1000	0	6	1e-7	0.001	0.01

Epoch – maximum number of epochs to train; Goal – performance goal; Max_fail – maximum validation failrues; Min_fail – minimum validation failures; Mu – adaptive value

This resulted in eight types of corrections surfaces (V_S). Consequently, DEM_SRTM elevations were modified according to the equation: 11HSRTMmodifled i=HSRTMUSGS +VSi$${H_{SRT{M_{{\rm{modifled}}{{\rm{}}_i}}}}} = {H_{SRT{M_{{\rm{USGS}}}}}} + {V_{{S_i}}}$$ where V_si is the ith corrections surface.

The eight modified HSRTM were exported to ArcGIS software to generate DEM files.

2.2.5.

Statistical analysis

The processed data of the whole area (study area), whether the vector (i.e., E, N, and H coordinates) or raster (DEM_GNSS), represent the reference (comparison) surface for the merged DEMs. The results were evaluated according to the differences in elevations, slopes, and aspects between the merged DEMs and the comparison surface. Likewise, the differences between the DEM_SRTM and DEM_GNSS were estimated. These differences were analyzed in the case of the absolute value and considering the signs. The absolute values statistics were achieved in terms of mean absolute error (MAE), Maximum value (Max.), Minimum value (Min.), and root mean square error (RMSE) as expressed in Equations 12 and 13: 12MAE=∑i=1n| vi |n$$MAE = {{\sum\nolimits_{i = 1}^n {\left| {{v_i}} \right|} } \over n}$$ 13RMSE=∑i=1n(vi)2n$$RMSE = \sqrt {{{\sum\nolimits_{i = 1}^n {{{\left( {{v_i}} \right)}^2}} } \over n}} $$ where v_i is the ith value of v, which represents the elevation differences between the GNSS processed data and the different DEMs produced; n is the number of observations.

In the case of considering the sign, the data were analyzed based on the box and whisker plot (Figure 16). It is very significant when a large number of datasets are analyzed. This is a standardized way of displaying the distribution of the dataset based on the five-number summary: the sample median, the lower and upper quartiles, any outliers (computed using the interquartile range), and the minimum and maximum values that are not outliers. Each box’s top and bottom edges are the upper and lower quartiles. The distance between the top and bottom edges is the interquartile range (IQR). The upper quartile corresponds to the 0.75 quantile, and the lower quartile corresponds to the 0.25 quantile. Outliers are the values more than 1.5 IQR away from the top or bottom of the box. The whiskers are lines that extend above and below each box. One whisker connects the upper quartile to the nonoutlier maximum value, and the other connects the lower quartile to the nonoutlier minimum value. For more information on this, see Wickham and Stryjewski (2011), Thomas (2010), and Dekking et al. (2005). To clarify, assume the following elevations samples: 54, 66, 69, 70, 71, 73, 73, 80, 84, 85, 86, 88, 90, 91, and 95 m. Accordingly: -

Min. = 54;

-

Max. = 91;

-

Lower Quartile = 70;

-

Upper Quartile = 88;

-

Median = 80.

Figure 16.

The elevation differences between DEMGNSS and DEM_SRTM are shown in Figure 17 and illustrated in Table 4. It can be observed that, for the data distribution, the median is closer to the top of the box (upper quartile); thus, the distribution is negatively skewed (skewed left). In contrast, Figure 18, which represents the elevation differences in the case of CDMMs, shows a positive skewed in all methods. Further, comparing the median value in this case (Table 5) with that obtained in DEM_SRTM, disregarding the sign, it shows that all the median values obtained by the CDMMs are less than the DEM_SRTM median. Furthermore, the upper and lower quartile values indicate that DEM_SRTM results are more dispersed than those from the CDMMs.

Figure 17.

Figure 18.

Considering the absolute value, the statistics demonstrate that the mean absolute error (MAE) for all CDMMs is better than the MAE in the case of DEM_SRTM. These results have been confirmed through the RMSEs, which are shown in Tables 4, 5, and Figure 19. Excluding the last method’s results, whose RMSE is equal to that produced by DEM_SRTM, all the CDMMs achieved higher accuracy than DEM_SRTM. The first and blend methods have provided the best results in terms of the median, MAE, and RMSE. This is because the first method takes advantage of the highly accurate DEM through the merging operation and the distance-weighted algorithm that is applied in the blend method.

Figure 19.

On the other hand, Figure 20 a–c shows the results of the PDMMs in the case of the zeros border. The results reveal a negative skew (skewed left), but the results of ANN in the same case (Figure 20 d) show that the median is closer to the bottom of the box, which means the distribution is positively skewed (skewed right). Besides, comparing these results (Table 6 a) and the best results of the CDMMs demonstrate that the data dispersion from the lower and upper quartile is significantly less than those of the CDMMs. The median, MAE, and RMSE show improvement in the merged DEM than the CDMMs. The MAE of the spline method is less than that produced by the IDW and kriging methods. However, the RMSE of the kriging method is better than both spline and IDW. Moreover, the ANN achieved the best results in the median, MAE, and RMSE.

Figure 20.

Table 6 b and Figure 21 display the results of the PDMMs in the case of the H border. The results show that the lower and upper quartile as well as the median values indicate fewer scatter than in all previous cases. Also, while the PDMMs related to the interpolation methods show a negative skew, the ANN shows a positive skew. According to the absolute values (Table 6 b), it can be concluded that the PDMMs, by the interpolation methods, achieved higher accuracy for the merged DEMs than that produced using the zeros border. Additionally, the kriging method achieved better MAE and RMSE than the spline method, which improved the IDW results except for MAE. Generally, the ANN proved its superiority as it outperformed all in terms of MAE and RMSE. The RMSEs of the PDMMs are illustrated in Figure 22.

Figure 21.

Figure 22.

To confirm the previous results, the elevation differences for all applied methods were classified into four categories. These categories were between less than one meter (highly accurate) and more than five meters (least accurate).

Table 7 and Figure 23 display the classification results of the DEMSRTM. In addition, the results of the CDMMs are shown in Figure 24 and illustrated in Table 8. The CDMMs achieved higher accuracy (diff. ≤ 1 m) than the DEMSRTM, which achieved a high percentage that is equal to the mean method in the medium accuracy category (i.e., 1 m < diff. ≤ 3 m). The results also show that the first and blend methods achieved the highest percentage in the first category, while the last and minimum methods achieved a high percentage in the fourth category (very low accuracy).

Figure 23.

Figure 24.

The PDMMs classifications are shown in Table 9 and presented in Figs. 25 and 26. The results for the zeros border (Table 9 a) show a higher percentage in the first category than the CDMMs, which produced a high percentage in the fourth category. Further, the spline method provided a better percentage in the first category than the kriging and IDW methods. The IDW produced the highest percentage in the second category (medium accuracy). The ANN outperformed all the other methods; it achieved the highest and lowest percentages in the first and fourth categories, respectively.

Figure 25.

Figure 26.

Referring to Figure 26 and Table 9 b, it can be noted that the interpolation methods, i.e., H border, achieved superiority over the case of the zeros border. Again, the spline method yielded a high percentage in the second group, but the kriging method was the best in the high accuracy and fourth categories. ANN and kriging achieved similar percentages in the fourth category, but ANN outperformed kriging in the high-accuracy category.

Generally, it can be concluded that, in interpolation methods, kriging performed best in the two applied cases (i.e., zeros border and H border). This could be due to kriging’s ability to take into account the spatial structure of the data, as proven by Rishikeshan et al. (2014) and Arun (2013). Notwithstanding, ANN outperformed kriging in the merged DEMs.

The elevation accuracy directly impacts the slope and aspects accuracies. This was investigated to validate the claim. The differences in slope and aspect degrees according to DEM_GNSS were calculated based on RMSEs for comparison among the different DEMs merging methods.

Figs. 27 and 28 show the slope and aspects distributions of the DEM_GNSS and DEM_SRTM, respectively. Also, the slopes and aspects of the DEMs obtained by the CDMMs are displayed in Figs. 29 and 30, respectively. Tables 10 and 11 provide the slopes and aspects statistics of the DEM_SRTM and CDMMs. The results show that the slope accuracy resulting from DEM_SRTM is better than the last, minimum, and maximum methods. Also, the accuracy of the aspects for DEM_SRTM is better than the mean method. While the blend method achieved the best slope in terms of RMSE, the minimum method was the worst.

Figure 27.

Figure 28.

Figure 29.

Figure 30.

Similarly, the results of the aspect show that the last method achieved the best accuracy, while the largest error was seen in the mean method. Furthermore, the first and blend methods reached the same results of aspects. Consequently, it can be noted that the applied technique that achieved the best slopes may not necessarily achieve the best aspects. The RMSEs of the slopes and aspects are displayed in Figs. 31 and 32, respectively.

Figure 31.

Figure 32.

The slopes and aspects distribution of the PDMMs in the two cases (zeros and H border) are shown in Figs. 33, 34, 35 and 36, respectively. Table 12, Figure 37, and Figure 38 display the RMSEs of the slope and aspects for the PDMMs. As shown in Table 12 a, zeros border case, although the slope results establish that kriging and IDW are better than the CDMMs, the ANN achieved the best RMSE of the slope differences. The aspect results of the ANN were the best among the PDMMs and CDMMs, except for the last method, which achieved the minimum RMSE. In addition, kriging, IDW, and spline methods resulted in aspects like the first and blend methods, but better than the mean method. However, they were the worst compared to the remaining CDMMs.

Figure 33.

Figure 34.

Figure 35.

Figure 36.

Figure 37.

Figure 38.

Table 12 b, which represents the PDMMs (H border case), shows that the slope and aspects accuracies are the best. Although kriging achieved the highest accuracy among the selected interpolation techniques, the ANN achieved the best RMSE for the slope and aspect differences.

The improvements in the elevations, slopes, and aspects were analyzed according to the DEM_SRTM for both CDMMs and PDMMs.

The improvements in the CDMMs are shown in Figure 39 and Table 13 a6. The results show no improvement in the elevation and slope for the last method. Also, some CDMMs show no improvements in accuracy, i.e., the minimum and maximum methods in the slopes and the mean method in the aspects. Nonetheless, some improvements are seen in other indicators. For example, the first and blend method in the elevations achieved the best improvement, while the blend improved in the slopes, and the last and minimum methods in the aspects.

Figure 39.

Table 13 b and Figure 40 illustrate the percentages of improvement in the case of the PDMMs. The results show that all the proposed methods significantly improved in the three terms, i.e., the elevations, slopes, and aspects. The elevation and slope achieved more improvements than those of the CDMMs. kriging, IDW, and spline methods (H border case) improved more in the elevations and aspects than for the ANN (zeros border case), which rather improved in the interpolation methods. Interestingly, kriging (in the two cases) improved than the IDW and spline methods.

Figure 40.

It is also important to note that the ANN (zeros border case) achieved better results than the applied interpolation methods in the two cases (i.e., the zeros border and H border) for the PDMMs, for slope accuracy improvement. Noticeably, the ANN (H border case) outperformed all other methods in the elevations, slopes, and aspects, with improvements reaching 76%, 55%, and 52%, respectively.

In view of the results, it can be noted that the PDMMs were better than the CDMMs in all statistical analyses. This could be due to the correction surface (Vs) applied for the PDMMs, which had a highly accurate overlapping zone before merging the DEMs. Thus, the CDMMs applied the merging technique only at the overlap area and retained the remaining zones without changes.

The ANN technique, in the case of the zeros border, performed the correction surface without additional points at the borders, but the interpolation methods cannot perform the interpolation properly without additional points at the border. This is because of the power of the backpropagation algorithm and further the ANN mechanism that depends on establishing a relation (i.e., mapping) between the inputs and outputs data based on the activation function. This mapping presents an ANN simulation that can produce any data out of the input range.

In the case of the H border, the interpolation methods realized the superiority of the ANN without additional points at the border. This is because the interpolation technique depends, basically, on the surrounding points (Oliver and Webster 1990). So, the GNSS points at the borders increase the accuracy and the confidence level of the interpolated points more than in the case of zero points.

Overall, the ANN reached the best results with additional GNSS points at the border. This is because the ANN simulation performs best with a high volume of data (Alemam et al. 2022). Therefore, the ANN is preferred as its cost-effective. However, if there exist additional observations at the borders (GNSS points), the interpolation methods can equally provide reasonable results.