A geomorphology based approach for digital elevation model fusion – case study in Danang City , Vietnam

Introduction Conclusions References


Introduction
A digital elevation model (DEM) is a digital model representing a surface which is presently used in many applications such as hydrology, geomorphology, geology and disaster risk mitigation.It is one of the essential inputs in modeling or simulating landscapes as well as dynamic natural phenomena such as flooding, soil erosion and landslides.Due to the important role of DEMs in terrain-related research and applications, it is necessary to create high-quality DEMs at various levels of details.DEM can be generated using photogrammetry, interferometry, ground and laser surveying and other techniques (Mukherjee et al., 2013a).Usually, aerial photos, high-resolution satellite data or field-surveyed spot height and light detection and ranging (lidar) data are used as inputs to generate high-resolution/high-quality DEMs.Surveying T. A. Tran et al.: A geomorphology-based approach for digital elevation model fusion have been several attempts to develop methodologies for enhancing quality of these global, free DEMs.
Several authors (e.g., Li et al., 2013;Ravibabu et al., 2010;Zhao et al., 2011;Suwandana et al., 2012;Mukherjee et al., 2013a;Czubski et al., 2013) have evaluated the accuracy of GDEM as well as SRTM and carried out comparative evaluation of two DEMs.Results from these studies indicated that, due to the inherent difficulties in acquiring satellite data both with the optical stereoscopic and the interferometric synthetic aperture radar (InSAR) technologies, global DEMs are not complete in and of themselves (Yang and Moon, 2003).Some authors (e.g., Reuter et al., 2007;Mukherjee et al., 2013a;Czubski et al., 2013;Fuss, 2013) have also evaluated the accuracy of global DEMs based on terrain characteristic.The vertical accuracy of these quasi-global DEMs vary depending on the terrain and land cover (Czubski et al., 2013).The main purpose of these studies was to verify the quality of global DEMs.However the unique characteristics and different factors affecting the vertical accuracy of optical stereoscopy and InSAR provide an opportunity for DEM fusion (Kaab, 2005).
This study proposes a geomorphological approach for DEM fusion based on evaluation of the accuracy of GDEM and SRTM in mountain slopes, valleys and flat areas.This approach was used to combine DEMs from different sources with appropriate weights to generate a fused elevation data.This could be an effective method to enhance the quality of global DEMs that have not been attempted in previous studies on DEM fusion (e.g., Corsetto and Crippa, 1998;Kaab, 2005;Karkee et al., 2008;Papasaika et al., 2011;Lucca, 2011;Fuss, 2013) 2 Study area This study was conducted in Danang city in the middle of central Vietnam (Fig. 1).The test site of 950 km 2 covers the inland area of Danang city and is characterized by elevation ranging from 0 to 1664 m a.m.s.l.Danang city is located on the Eastern Sea coast, extending from 15 • 55 N to 16 • 14 N and 107 • 18 E to 108 • 20 E. The topography of this area has great variation from flat to mountainous regions.Due to varying of topography and geomorphology, the optical stereoscopy technique used to generate GDEM as well as the InSAR technique used in SRTM show different representation on DEM data, and contain inherent anomalies that need to be detected and minimized.
There are few studies in this area using global, free DEMs such as GDEM or SRTM.Ho and Umitsu (2011) and Ho et al. (2013) developed a landform classification method and flood hazard assessment of the Thu Bon alluvial plain, central Vietnam.In their study, the authors used SRTM as an input DEM source and applied bias elimination method to correct surface elevation data to the height of bare-earth surface.However, SRTM with low resolution (90 m) may not give sufficient terrain information.Also, the InSAR technique used in SRTM may fail to provide reliable estimate of elevation if images contain layovers, nonlinear distortion of the images due to slanted geometry of the radar sensing and shadows, or suffer from temporal decorrelation and changes in atmospheric conditions between two acquisitions (Karkee et al., 2008).Although Ho et al. (2013) already masked the high and upland areas and focused only on a low-lying alluvial plain, their research did not discuss methods to enhance accuracy of free DEMs, especially in the areas that have high topographic relief.

DEM data sets
The global, free DEMs used in this study include GDEM Version 2 (http://earthexplorer.usgs.gov)and SRTM Version 4.1 (http://www.cgiar-csi.org).GDEM Version 2 released in October 2011 has the resolution of 30 m. GDEM data were compiled from over 1.2 million scene-based DEMs covering land surface between 83 • N and 83 • S latitudes.GDEM was generated from ASTER optical satellite images using stereoscopy technique with differing sensor look angles.The Terra spacecraft used in ASTER GDEM is capable of collecting in-track stereo using nadir-and aft-looking near-infrared cameras (ASTER GDEM Validation Team, 2011).DEM from such optical satellite images as GDEM usually contains some height errors because of cloud coverage.ASTER GDEM Version 2 was improved with respect to Version 1 (released on June 2009) due to a better data-processing algorithm and additional data used during the processing.However, the revised version still contains anomalies and artifacts which need to be corrected before being used in any application, especially on a local scale (ASTER GDEM Validation Team, 2011).
SRTM Version 4.1 has been obtained from the Consortium for Spatial Information (CGIAR-CSI; http://www.cgiar-csi.org).The DEM data were derived from the 11-day Shuttle Radar Topographic Mission that flew in February 2000, and have provided publicly available elevation surface data for approximately 80 % (from 60 • N to 56 • S) of the world's land surface area (Reuter et al., 2007).The SRTM elevation data are derived from X-band and C-band InSAR sensor.The first release of SRTM was provided in 1 • DEM tiles in 2003.When the data were processed by NASA and the USGS, they were made available at 1 arcsec resolution (approximately 30 m) for the United States, and 3 arcsec resolution (approximately 90 m) for the rest of the world.The Consortium for Spatial Information of the CGIAR (CGIAR-CSI) is offering post-processed 3 arcsec DEM data for the globe.The original SRTM has been subjected to a number of processing steps to provide seamless and complete elevation surface for the globe.In its original release, SRTM data contained regions of no data, specifically over water bodies (lakes and rivers) and in areas where insufficient textural detail was available  projected in a Vietnamese projection named VN2000.In this study, the DEM generated from contour and spot height elevation is referred to as the "reference" DEM.Firstly a DEM was generated from the contour map using the regularized spline with tension (RST) algorithm.The RST interpolation is considered as one of the effective interpolation methods available for elevation data (Hofierka et al., 2002).The RST method is based on the assumption that the approximation function should pass as closely as possible to the given data and should be as smooth as possible (Mitasova et al., 1995).RST interpolation was carried out in GRASS GIS open-source software (http://grass.osgeo.org).However, the contour lines do not cover the whole area of Danang city.In flat areas with elevation less than 10 m, there are no contour lines.A large number of spot height data are available for flat areas (more than 190 000 elevation points), and inverse distance weighting (IDW) interpolation was applied to generate the DEM where contour data are not available and merged with DEM generated using RST with contour data for hilly areas.This reference DEM was also generated at a resolution of 30 m.The RMSE of the reference DEM comparing to spot height data is 1.66 m.Some statistical data on the global DEMs and reference DEM are shown in Table 1.The mean elevation and standard deviation (SD) in GDEM and SRTM are analogous to the reference DEM.Due to some artifacts located in GDEM, the maximum elevation value of GDEM (8016 m) shows significant dissimilarity.Compared to GDEM, the SD of SRTM (304.6 m) is almost similar to the reference DEM (302.6 m).

Methodology
SRTM was interpolated from 90 to 30 m resolution in order to compare with other DEM sources.The artifacts in GDEM were eliminated using the fill and feather method (Dowding et al., 2004).DEM alignment was also carried out in order to co-register GDEM and interpolated SRTM with respect to the reference DEM.Next, both GDEM and SRTM were evaluated in terms of vertical and horizontal accuracy.The quality of each DEM was also assessed according to different topographic conditions.The result of evaluation has been used to devise an appropriate DEM fusion method considering various factors responsible for degradation of data quality.Basically, there is a difference between the digital surface model (DSM) like GDEM, SRTM and the digital terrain model (DTM) that refers to the bare-earth surface.The overestimations as well as underestimated elevation values in GDEM and SRTM need to be detected and corrected by comparing these elevation data to the reference DEM on the basis of the geomorphology and land cover map.In the case of land cover category, the offsets were calculated by taking mean values of the difference in elevation between the global DEMs and reference DEM.The corrected GDEM and SRTM were used as input data for the DEM fusion process.The landform classification map was generated from SRTM to determine the area suitable for different fusion methods.The algorithm used in the DEM fusion process is weighted averaging based on geomorphologic characteristics.In relatively flat areas, the higher weight was used for SRTM and lower weights for GDEM.In the mountainous areas, SRTM and GDEM were weighted equally.The higher weight was applied for GDEM in the valley areas, because of the limitation of SRTM in those areas.The output-fused DEM was filtered using a denoising algorithm according to Sun et al. (2007).
Finally, the fused DEM was compared to the reference DEM to assess the efficiency of the DEM fusion method.
The data processing described above is shown in Fig. 2. The data fusion workflow includes four main steps, namely, pre-processing, DEM quality assessment, bias elimination and DEM fusion.

Pre-processing
It is observed that SRTM has anomalies in the coastal area and some small areas inland with negative values; 377 pixels show negative values and cover about 0.34 km 2 .These pixels were filled by averaging elevation of 3 × 3 neighboring pixels.SRTM and GDEM have been converted from geographic coordinates to UTM_WGS84_zone 49N projection.The reference DEM was also converted from VN2000 to UTM_WGS84_zone 49N projection.The vertical datums used in the global DEMs and reference DEM are different.The global DEMs use EGM96 vertical datum, while the reference DEM uses the Vietnamese vertical datum named Hon Dau-Hai Phong, which is related to m.s.l. in Hon Dau Island, Hai Phong province, Vietnam.An offset 1.5 m downwards was applied to convert the global DEMs from EGM96 to Hon Dau-Hai Phong vertical datum.
SRTM was interpolated from 90 to 30 m using the RST algorithm, which is available in GRASS GIS as the r.resamp.rstfunction.RST interpolation not only re-samples the DEM to higher resolution but also reduces the staircase effect in the original SRTM and smoothens the DEM surface.Figure 5a    the reference DEM before and after interpolation.The interpolated SRTM also has better RMSE and correlation to the reference DEM than the original 90 m data (Table 3).
GDEM has some artifacts in the western mountain part of Danang city, with extreme high-elevation values.These artifacts may be caused due to cloud coverage that is very common in optical satellite data.These artifacts are the main reason for high RMSE (75.6 m) observed in raw GDEM (Table 2).The artifacts in GDEM need to be eliminated before further processing.Several algorithms for void filling have been proposed, such as kriging, spline, IDW (Reuter et al., 2007), moving window (Karkee et al., 2008), fill and feather (Dowding et al., 2004) and delta surface fill (Grohman et al., 2006).All the void-filling algorithms can be categorized into three groups, namely, interpolation, moving window and fill and feather (F & F).The F & F method proposed by Dowding et al. (2004) was applied in this study to fill artifacts in GDEM.In the F & F approach, an artifact is replaced with the most accurate digital elevation source available with the void-specific perimeter bias removed (Grohman et al., 2006).The artifacts were detected by overlaying the slope map of GDEM and the difference elevation map between GDEM and the reference DEM, and digitizing from the anomalies that can be visualized from the overlaying display.SRTM was chosen as auxiliary data to fill the artifacts for GDEM.
After filling these artifacts, the surface is feathered to mitigate any abrupt change (Grohman et al., 2006).In this case study, DEM surface will be feathered in the final step of data processing using a filtering algorithm.As the result, GDEM after filling artifacts has a RMSE of only 14.9 m.The scatter plot of GDEM after applying F & F also shows a good correlation to the reference DEM, while the original one has several outliers (Fig. 3).Comparing to original GDEM, it can also be seen that most of the artifacts were eliminated.

DEM quality assessment
The horizontal accuracy of the global DEMs was evaluated by comparing the extracted stream networks (Fig. 4).Stream networks extracted from the reference DEM, GDEM and SRTM indicate that SRTM has a horizontal difference of about 15 m, and GDEM has a difference of around 30 m with respect to the reference DEM.Therefore, GDEM was shifted one pixel to the east, and SRTM was shifted half a pixel to the west, in order to align all input DEMs before fusion process.Figure 5 compares the profiles of GDEM, SRTM and the reference DEM before and after shifting.The ridge lines as well as canyon bottoms in GDEM and SRTM become more similar to the reference DEM.In Table 2, GDEM after shifting shows better RMSE and correlation with the reference DEM as compared to before shifting.
In this study area, the RMSE of GDEM and SRTM with respect to the reference DEM was observed as 14.9 and 14.8 m, respectively (Tables 2 and 3).The correlation coefficient (R 2 ) of GDEM in the whole area is 0.9976, while this value in the original SRTM is 0.9979.The accuracy of the individual DEM should be considered based on the different topographic condition.Figure 6 shows the correlation coefficients of each global DEM in flat and mountain areas.In mountain areas, GDEM and SRTM have a similar correlation with the reference DEM (0.9966 and 0.9969, Fig. 6b).However, in some specific areas, especially in the steep valleys, GDEM provides better accuracy than SRTM.The circled areas in Fig. 5 show that GDEM preserves the considerable details of topography in the valley areas, while SRTM is ineffective in those areas.In such valley areas, SRTM seems to suffer from layover and shadow effects.In the case of a very steep slope, targets in the valley have a larger slant range than related mountain tops; consequently the fore-slopes are "reversed" in the slant range image.This is referred to as layover effect when the ordering of surface elements on the radar image is the opposite of the ordering on the ground (European Space Agency, https://earth.esa.int/applications/data_util/SARDOCS/spaceborne/RadarCourses/Radar_CourseIII/layover.htm).Radar shadow is caused when a slope is away from the radar illumination with an angle that is steeper than the sensor depression angle (European Space Agency, https://earth.esa.int/applications/data_util/SARDOCS/spaceborne/RadarCourses/Radar_CourseIII/shadow.htm).In such areas, SRTM may not provide sufficient information, compared to GDEM or other DEM sources.In relatively flat areas, the correlation coefficient between SRTM and the reference DEM (R 2 = 0.8504) is better than GDEM (R 2 = 0.5578) (Fig. 6a).This is because degradation of the elevation estimate of GDEM in the area has low topographic relief.In the profile of Fig. 7, it can be seen that GDEM has many spikes and unstable elevation values in this flat area, while SRTM shows similar trends to the reference DEM.
The difference elevation maps of the global DEMs were also generated by subtracting GDEM and SRTM values from the reference DEM.Both GDEM and SRTM show high vertical error in mountain areas, and lower vertical error in flat areas (Fig. 8).These errors occur because of the forest cover in mountain areas and due to some limitations of the sensing techniques used to generate DEM in high-relief areas.The profile of SRTM from the difference elevation map in flat areas is closer to 0 m line (Fig. 8), while GDEM contains higher difference and spikes that affect the quality of GDEM significantly.

Minimizing DEM bias effect
The topographic height variation between the global DEMs and reference DEM is caused due to the differences in vertical datum used and in primary data collection methods.Vertical datum is one of the reasons for difference in elevation between the global DEMs and reference DEM.In addition, both GDEM and SRTM, which were generated from satellite data, are DSMs, while the reference DEM is considered a bare-earth DTM; this difference also introduces the bias offsets depending on the land cover.
Firstly, the global DEMs were converted to the Hon Dau-Hai Phong vertical datum.According to the Vietnam Land Administration, the global EGM96 model is almost similar to the Vietnamese vertical datum; 97 % of data shows the height difference around 1.5 m, while only 3 % of data shows higher than 1.5 m (Nguyen and Le, 2002).Therefore, an offset of 1.5 m was subtracted from the global DEMs, considering height difference between EGM96 and Vietnamese vertical datum.Secondly, the height offsets of the global DEMs were determined based on the land cover map.Because the SRTM data were derived in 2000 and GDEM data were collected from millions of ASTER images from 1999 to 2009, a land cover map of Danang city in 2001 was used to calculate the height offsets for the global DEMs.These offsets were calculated based on the difference elevation maps of GDEM and SRTM with respect to the reference DEM considering land cover.This was done using the r.statistics function in GRASS GIS.The mean elevation differences on each land cover type were calculated, and used as offsets to verify elevation for GDEM and SRTM (Table 4).As the result, GDEM has the highest difference in the water body (4 m).This error is  is too small to create a radar image.From global assessment of the SRTM data, voids were found to be very common in mountainous areas, as well as in very flat areas especially in deserts (Zandbergen, 2008).SRTM V4 used in this study already dealt with the water body problem using a number of interpolation techniques and void-filling algorithms (Reuter et al., 2007).Therefore, the error of SRTM in water bodies currently is only 0.4 m (Table 4).Based on the above investigations, the elevations for GDEM and SRTM with respect to the reference DEM were corrected by subtracting GDEM and SRTM from the elevation offsets for each land cover type (Table 4).The calculation was executed by the r.mapcalc function in GRASS GIS software using the land cover map as the base.The corrected GDEM and SRTM were used as input data for DEM fusion processing.
After removing the offsets, GDEM and SRTM were compared to the reference DEM again to make better input for DEM fusion processing.The mean value of GDEM and SRTM with respect to each elevation value in the reference DEM was calculated.Figure 9a shows the behavior of the global DEMs with respect to the reference DEM, from flat to mountainous areas.In the A and C area (Fig. 9b and d), the mean elevation of SRTM is closer to the reference DEM, while the profile of GDEM shows higher error.In the case of the B area (Fig. 9c), both SRTM and GDEM show the good correlation to the reference DEM.In Fig. 9e, the profile of GDEM is comparable to the reference DEM in this mountainous area.From this analysis, it is evident that using a global data fusion for the whole area is not a good solution.Appropriate weights for the DEM fusion process need to be considered depending upon the topographic context, and is used as the basis for DEM fusion in this study.

DEM fusion algorithm
Both GDEM and SRTM contain intrinsic errors due to primary data acquisition technology and processing methodology in relation with a particular terrain and land cover type (Mukherjee et al., 2013a).The optical stereoscopy technique used in GDEM is limited by the cloud coverage, radiometric variation and low levels of texture (Karkee et al., 2008), while the InSAR technique used in SRTM may not work well in the case of shadowing, layovers or complex dielectric constant (Reuter et al., 2007).Combination of two data can take into account the advantages of each DEM source and provide complementary inputs to enhance the quality for the global DEMs.DEM fusion workflow combines the weighted averaging and denoising algorithm (Sun et al., 2007).

Weighted averaging
Several authors have proposed fusion methods for digital elevation data.Karkee et al. (2008) carried out a fusion between GDEM and SRTM using fast Fourier transformation (FFT) combining with frequency domain filtering.Papasaika et al. (2011) have proposed an approach that performs DEM fusion using sparse representations.Lucca (2011) examined different DEM fusion methods, such as weighted averaging and collocation prediction, and compared the result to lidar DSM to assess the improvement of DEM fusion.Fuss (2013) has developed a DEM fusion algorithm from multiple, overlapping DEMs, using slope thresholding, K-means clustering and filtering of elevations.Tran et al. (2013a, b) have given a fusion method by selecting appropriate DEM-sourcebased geomorphological conditions.The most frequent DEM fusion method that has been suggested is weighted averaging.The weighted mean (x) of a nonempty set of data {x 1 , x 2 , . . ., x n } with nonnegative weights {ω 1 , ω 2 , . . ., ω n } (Papasaika, 2012) is shown: where x 1 , x 2 , . . ., x n are the input DEMs.ω 1 , ω 2 , . . ., ω n are the weights for DEM fusion.However, weighted averaging applied in previous studies referred to in the earlier section considers weights based on the accuracy of the whole raster DEM source.Each raster DEM x 1 , x 2 , . . ., x n is used as one input data for weighted averaging.Actually, the DEM accuracy also changes depending upon the topographic context.Therefore, in this research, a new method for DEM fusion using weighted averaging based on geomorphologic characteristics was proposed.Firstly, a landform map was extracted from SRTM.This landform classification method was done according to Dickson and Beier (2006).The algorithm is based on the topographic position index (TPI) and slope map.In general, TPI allows classifying landscape into discrete landform categories by comparison of individual cell heights with an average height of neighboring cells (Czubski et al., 2013).The TPI-based landform classification method according to Dickson and Beier (2006)  In this study, three categories demarcated from the landforms classification result, namely, mountain slopes (include ridge lines and steep slopes), valleys and flat areas (Fig. 10).
In order to determine the weight for the global DEMs on each landform class, the following equation (Hengl and Reuter, 2009) was applied: where w i is the weight for each DEM source for a given landform unit and a is the given accuracy parameter for the DEM for a given landform unit.
Terrain-related parameters were used to determine the weighting scheme for DEM fusion.Firstly, slope error (difference in slope between the global DEMs and reference DEM) was used to compare the accuracy of GDEM and SRTM in flat, valley and mountain slope areas.On each landform unit, the mean of absolute error (MAE) from the slope error map was calculated.The result is shown in Table 5.
In flat areas, GDEM has many overestimates and unstable elevation Therefore slope error of GDEM is larger than SRTM in this area.The weight used for GDEM can be determined according to Eq. ( 2): w 1 = 1/(2.1) 2 = 0.22; and the weight for SRTM can be shown as w 2 = 1/(1.6) 2 = 0.39.It can  be seen that w 2 ≈ 2 • w 1 ; therefore the following formula was applied for DEM fusion in flat areas: In mountain slope areas, the similar way was applied to calculate weight for DEM fusion, using MAE of slope error.
In this case, GDEM SRTM have almost same MAE (6.08 and 6.1 • ).Therefore, the same weights were applied for GDEM and SRTM in mountain slope areas (w 1 = w 2 ).
The following equation was used in mountain slopes: In valley, GDEM and SRTM also have the similar MAE of slope error (5.8 and 5.7 • ).However, considering the topographic characteristic in some steep valleys, it can be seen that SRTM is ineffective in representing the valley bottom, while GDEM is still more correlative to the reference DEM (Fig. 5).In the case of valley landforms, Slope Variability (SV) (Popit and Verbovsek, 2013) was used to determine weight for DEM fusion.SV was calculated by the distance between maximum and minimum slope in the neighborhood of 3 × 3 pixels.SV errors of GDEM and SRTM with respect to the reference DEM was calculated.GDEM has a MAE of SV error of about 5.6, and SRTM has an error of about 7.3 • .The weight for GDEM was calculated according to Eq. ( 2): w 1 = 1/(5.6) 2 = 0.032; and the weight for SRTM is calculated as w 2 = 1/(7.3) 2 = 0.018.It can be observed that w 1 ≈ 2 • w 2 ; therefore the following formula was used for DEM fusion in valley: Fused DEM = (GDEM • 2 + SRTM)/3. (5) The weighted averaging method based on the landform classification map is shown in Fig. 11.

Filtering the noises for the fused DEM
The fusion of different DEMs is problematic, since the DEMs are obtained from different sources and have different resolutions as well as accuracies (Lucca, 2011).The bias eliminations for GDEM and SRTM also use different offsets depending upon the land cover.Different weights have been used for DEM fusion in each landform type.Therefore, it is essential to filter the fusion DEM to reduce the mismatched and noisy data.In this study, the denoising algorithm (Sun et al., 2007) was used to minimize the noise effect.The level of denoising is controlled by two parameters, namely, the threshold (T ) that controls the sharpness of the features to be preserved, and the number of iterations (n) that controls how much the data are smoothed.The optimum settings depend upon the nature of the topography and of the noise to be removed (Stevenson et al., 2009).Sun's algorithm (Sun et al., 2007) has been implemented in GRASS GIS as an add-on (r.denoise).In this denoising process, the topographic feature need to be preserved as far as possible in the fused DEM, so the parameters that were used are T = 0.95 and n = 5.As the result, the fused DEM becomes more smooth and the mismatched surfaces are minimized.The profile of the fused  DEM is also very much comparable to the reference DEM (Fig. 12).

Results and discussions
Weighted averaging based on the landform classification map has been verified as an effective method for DEM fusion.
The accuracy of the fused DEM can be evaluated by statistical analysis such as RMSE, MAE and linear regression.The MAE and RMSE of the fused DEM were much improved compared to the available global DEMs.The RMSE was reduced from 75.6 m in original GDEM, 14.9 m in GDEM after removing artifacts and 13 m in GDEM after bias elimination to 11 m in the fused DEM.In SRTM, the RMSE was reduced from 14.8 m in the original SRTM and 11.4 m in the processed SRTM to 11 m in the fused DEM (Table 6).
The linear regression between the fused DEM and reference DEM also shows the significant correlation between two DEMs, with R 2 = 0.9986 (Fig. 13).Comparing to original data with correlation coefficient for GDEM and SRTM of 0.9976 and 0.9979, respectively, it can be, therefore, be concluded that the fused DEM shows better correlation with the reference DEM.
Statistical comparison of vertical accuracy of GDEM, SRTM and the fused DEM is shown in Table 6.The minimum error, maximum error, MAE and RMSE of the fused DEM show improvement when compared with GDEM and  there is significant improvement in quality of global DEMs using the proposed DEM fusion algorithm.
The slope, profile curvature and tangential curvature maps were extracted from GDEM, SRTM and the fused DEM.Then the error maps with respect to the reference DEM were created in each terrain parameter (Table 7).Comparing to GDEM and SRTM, the fused DEM has smaller MAE and SD, and better correlation with the reference DEM. Figure 16 shows the slope, profile curvature and tangential curvature maps from the fused DEM.In these DEM derivative parameters, no major anomaly or terrace artifacts can be seen in the transition zones between landform classes.
Aspect is calculated as circular degrees clockwise from 0 to 360 • , and it is therefore difficult to compare quantitatively (Deng et al., 2007).In order to assess the accuracy in aspect as well as slope, unit normal vector (NV) of topographic surface was considered.The NVs of the global DEMs and fused DEM were computed from slope and aspect values of respective DEMs.The NVs from these DEMs then were compared with the reference DEM to determine the angular difference between two NVs (Fig. 17).The NV of the terrain surface (T ) can be calculated as below as suggested by Hodgson and Gaile (1999): where x = sin(aspect) • sin(slope), y = cos(aspect) • sin(slope) and z = cos(slope).
To derive the three-dimensional angular difference between two unit NVs (T and S) pointing away from the same origin, the following formula (Hodgson and Gaile, 1999) was applied: The result of angular differences of NV is shown in Table 8.As a result, the fused DEM has a smaller mean error than GDEM and SRTM, and the SD of the fused DEM is also comparable with the global DEMs.
The topographic roughness index (TRI) was also considered to assess the quality of the fused DEM.In this study, the TRI was used as the amount of elevation difference among the adjacent cells of a DEM (Mukherjee et al., 2013b).The residuals in elevation between a grid cell and its eight neighbors were derived, and the RMSE of the elevation differences was calculated as the TRI.The TRI of the reference DEM and GDEM, SRTM as well as the fused DEM show a correlation coefficient of 0.71, 0.75 and 0.76, respectively (Table 7).The TRI derived from the fused DEM compares well with the reference DEM as compared with GDEM and SRTM.(Hodgson and Gaile, 1999).

Conclusions
Global, free DEMs generated from remote-sensing data always have some vertical and horizontal errors.Assessing the quality of global DEMs and validating their accuracy before use in any application is very important.In this study, the accuracies of GDEM and SRTM were determined based on height differences with the reference DEM.The artifacts with extreme high-elevation values in GDEM were eliminated by using SRTM as auxiliary data.River networks extracted from both DEMs that were used to detect and correct the horizontal errors for global DEMs can make better co-registration.The bias effect caused by tree-top canopy and building on global DEMs was also calculated by comparing these DSMs with the elevation from the reference DEM.A land cover map of Danang city in 2001 was used to calculate the height difference of GDEM and SRTM on each land cover type.
Once the bias offsets were determined, effort was made to correct the elevation of these DEMs with respect to the bareland surface.
Based on global DEM assessment in Danang city, it is observed that the accuracy of GDEM and SRTM varies depending upon the geomorphological characteristics of the target area.Fusion between two global DEMs using a geomorphological approach is an appropriate solution to enhance the quality of free DEMs for Danang city, Vietnam.The data fusion technique was applied by weighted averaging of GDEM and SRTM based on the topographic context.The weighting scheme was determined according to accuracy parameters, including MAE of slope and slope variability.The weights used for each DEM were changed locally according to the landform types.The results were compared with the reference DEM to discuss accuracy and impact of landform in variation on DEM quality.Terrain-related parameters such as slope, curvature, TRI and NV of topographic surface were considered to assess seriously the quality of the fused DEM.Results indicate that the fused DEM has improved accuracy compared to individual global DEM and most artifacts are successfully eliminated.The proposed method supports the effective utilization for the areas where the better-quality DEM is not available.
In future work, the more robust weighting scheme needs to be considered by defining a greater number of landform types.In this regard, the landform classification method may also need to be improved further.In the future, we plan to investigate landform classification using r.geomorphon, a new add-on that is available in GRASS 7. A "geomorphon" is a relief-invariant, orientation-invariant and size-flexible abstracted elementary unit of terrain (Stepinski et al., 2011).This landform classification map will not only be a good way to compare the height errors in micro-geomorphological classes, but will also help to compare terrain parameters extracted from fused global DEMs and the reference DEM.
The difference in elevation between DEM and DSM is useful for estimating the canopy height, especially in areas used for silviculture.Further investigation on bias effect introduced by land cover and silviculture needs to be carried out.The relationship between land cover and geomorphology also need be studied in the future, to understand the impact of topographic condition on land cover change.Several new satellite data, including ALOS-2 PRISM and PALSAR-2 (http://www.eorc.jaxa.jp/ALOS/en/index.htm),need to be incorporated to enhance the methods for multiresolution DEM fusion based on a better understanding of characteristics of DEMs derived from multiple sources.

Figure 1 .
Figure 1.Location of study area and topographic overview.

Figure 3 .
Figure 3. Correlation between GDEM and the reference DEM before (left panel) and after (right panel) filling voids.

Figure 4 .
Figure 4. Comparing stream networks of the global DEMs and reference DEM before (top panels) and after (bottom panels) shifting DEM: (a) GDEM, (b) SRTM.

Figure 5 .
Figure 5. Comparing GDEM and SRTM to the reference DEM: (a) before re-interpolation SRTM and shifting data, (b) after reinterpolation SRTM and shifting data.

Figure 8 .
Figure 8. Difference elevation of GDEM and SRTM with respect to the reference DEM from mountain to flat areas.

Figure 12 .
Figure 12.Result of a denoising algorithm (Sun et al., 2007) on the fused DEM.

Figure 13 .
Figure 13.Correlation between the fused DEM and reference DEM.

Figure 14 .Figure 15 .
Figure 14.Difference in elevation between the fused DEM and reference DEM.

Figure 17 .
Figure 17.Normal vector of a topographic surface (a) and the angular difference between two normal vectors (b)(Hodgson and Gaile, 1999).

Table 1 .
General information on the global DEMs and reference DEM (all the negative values were filled by neighboring pixels).

Table 2 .
Results of GDEM after filling artifacts and shifting.

Table 3 .
SRTM before and after interpolation to 30 m.

Table 4 .
The mean errors of GDEM and SRTM according to the land cover map.Unit: m.

Table 5 .
Mean of absolute error (MAE) from slope error maps of GDEM and SRTM on each landform unit.

Table 6 .
General statistics for the error of GDEM, SRTM and the fused DEM.Unit: m.

Table 7 .
Comparison of differences in some terrain parameters of GDEM, SRTM and the fused DEM with respect to the reference DEM.

Table 8 .
Result of angular difference of unit NVs between the global DEMs, fused DEM and reference DEM.