A Multi-Sensor Comparative Analysis on the Suitability of Generated DEM from Sentinel-1 SAR Interferometry Using Statistical and Hydrological Models

Digital elevation model (DEM) plays a vital role in hydrological modelling and environmental studies. Many essential layers can be extracted from this land surface information, including slope, aspect, rivers, and curvature. Therefore, DEM quality and accuracy will affect the extracted features and the whole process of modeling. Despite freely available DEMs from various sources, many researchers generate this information for their areas from various observations. Sentinal-1 synthetic aperture radar (SAR) images are among the best Earth observations for DEM generation thanks to their availabilities, high-resolution, and C-band sensitivity to surface structure. This paper presents a comparative study, from a hydrological point of view, on the quality and reliability of the DEMs generated from Sentinel-1 data and DEMs from other sources such as AIRSAR, ALOS-PALSAR, TanDEM-X, and SRTM. To this end, pair of Sentinel-1 data were acquired and processed using the SAR interferometry technique to produce a DEM for two different study areas of a part of the Cameron Highlands, Pahang, Malaysia, a part of Sanandaj, Iran. Based on the estimated linear regression and standard errors, generating DEM from Sentinel-1 did not yield promising results. The river streams for all DEMs were extracted using geospatial analysis tool in a geographic information system (GIS) environment. The results indicated that because of the higher spatial resolution (compared to SRTM and TanDEM-X), more stream orders were delineated from AIRSAR and Sentinel-1 DEMs. Due to the shorter perpendicular baseline, the phase decorrelation in the created DEM resulted in a lot of noise. At the same time, results from ground control points (GCPs) showed that the created DEM from Sentinel-1 is not promising. Therefore, other DEMs’ performance, such as 90-meters’ TanDEM-X and 30-meters’ SRTM, are better than Sentinel-1 DEM (with a better spatial resolution).


Introduction
Digital elevation model (DEM) is one of the essential geospatial data tools used in geological, geographical, geomorphological, and environmental applications. It allows extracting essential terrain features such as slope, aspect, elevation, and many more [1][2][3][4]. The accuracy of such features is directly associated with the accuracy of DEM. There are various geospatial technologies with different technical complexity, geographical coverage, geometric and geographic accuracy, cost, user-effort, and time requirements for generating DEMs. Ground surveying, photogrammetry, light detection and ranging (LiDAR) systems, satellite optical imagery, and SAR interferometry are the best operational technologies that have been used to generate the DEMs for different applications. Compared with the remote sensing (RS) techniques, ground surveying is a more accurate method to create a DEM [5,6]. However, RS-based technologies are relatively less expensive timely and can cover larger areas shortly [7].
Among the RS sensor technologies, SAR Interferometry has shown its advantages for DEM generation, such as sensitivity of SAR signal to the surface physical and topographic properties, independence of weather and illumination conditions, and robust and accurate analytical approaches [8,9]. With the open data strategy by the national space and Earth observation agencies, both scientific and commercial end-users have had access to high spatial and temporal SAR imagery from various wavelengths, polarization, and orbits. Thus, remote sensing has received lots of attention in the last decade [7,10].
Because of the different product types and sensor modes, Sentinel-1 has many practical applications such as subsidence, landslide and, land cover flood detection [11][12][13]. Simultaneously, many researchers have been applying Sentinel-1 for DEM creation (because of freely availability) [14]. Since this study, for the first time, tries to test and compare the accuracy and validity of a derived DEM from Sentinel-1 data with the other DEMs, it can be significant. Therefore, this study opens the researchers' eyes towards the DEM generation from Sentinel-1 satellite data.
As a very successful example, many researchers worldwide tend to create DEMs from Sentinel-1 SAR data. Braun [45] created a DEM from Sentinel-1 data conducted with the European Space Agency (ESA); since then, it has been a key source for researchers worldwide, but there is no comparative study yet to test its accuracy and validity the other products, especially freely available DEMs. For generating a reliable DEM, the perpendicular baseline (PB) of SAR image data should be between 150 to 300 m [46], while in the best condition, the PB of Sentinel-1 cannot exceed 150 m. Therefore, this can be a significant problem. For this reason, the researchers should be careful about using Sentinel-1 to generate DEM. It is worth mentioning that the other presented DEMs' capability and applicability are not in the scope of the current study and are just used to test the accuracy of the created DEM.
The main objective of the current study was to generate a DEM from Sentinel-1 satellite imagery by using the InSAR technique and also to compare it with DEMs of airborne synthetic aperture Radar (AIRSAR) product (with 10 × 10 cell size), advanced land observing satellite phase array L-band synthetic aperture radar (ALOS-PALSAR) imagery (with 12.5 × 12.5 m spatial resolution), TanDEM-X (with 90 × 90 cell size) and SRTM (with 30 × 30 cell size) for two different study areas.

Description of the Study Area
When a comparative study is done, every aspect, even the study area's selection, must be considered. For example, regarding SAR data and DEM generation, collecting data on hot days Sensors 2020, 20, 7214 3 of 24 and dense vegetation coverage (especially Band-C) result in low coherence extraction. In this study, two study areas of a part of the Cameron Highlands, Pahang, Malaysia, and Sanandaj county, Kurdistan, Iran, were selected to carry out this comparative study. The Cameron Highlands is mostly covered by dense and sparse vegetation, while vegetation coverage for Sanandaj county is negligible. The Cameron Highlands is located between Longitudes 101 • 20 00" to 101 • 27 10" East and Latitudes 4 • 23 30" to 4 • 31 10" North ( Figure 1). This scope covers an area of 81.249 km 2 and is placed in Peninsular Malaysia's central part and the main range of geological units of Malaysia [47,48]. There are two prominent outcrops in the Cameron Highlands, the first one is acid intrusive (undifferentiated), which covers most of the study area of 60.99 km 2 , and the second formation (Silurian-Ordovician) consists of slate, schist, limestone with minor intercalation of sandstone, and phyllite [49]. Based on the tropical rainfall measuring mission (TRMM) data, the Cameron Highlands' annual rainfall fluctuates between 3800 mm to 4200 mm annually [49,50].
considered. For example, regarding SAR data and DEM generation, collecting data on hot days and dense vegetation coverage (especially Band-C) result in low coherence extraction. In this study, two study areas of a part of the Cameron Highlands, Pahang, Malaysia, and Sanandaj county, Kurdistan, Iran, were selected to carry out this comparative study. The Cameron Highlands is mostly covered by dense and sparse vegetation, while vegetation coverage for Sanandaj county is negligible. The Cameron Highlands is located between Longitudes 101°20′00″ to 101°27′10″ East and Latitudes 4°23′30″ to 4°31′10″ North ( Figure 1). This scope covers an area of 81.249 km 2 and is placed in Peninsular Malaysia's central part and the main range of geological units of Malaysia [47,48]. There are two prominent outcrops in the Cameron Highlands, the first one is acid intrusive (undifferentiated), which covers most of the study area of 60.99 km 2 , and the second formation (Silurian-Ordovician) consists of slate, schist, limestone with minor intercalation of sandstone, and phyllite [49]. Based on the tropical rainfall measuring mission (TRMM) data, the Cameron Highlands' annual rainfall fluctuates between 3800 mm to 4200 mm annually [49,50].
With a population of almost 432 thousand, the city of Sanandaj is the biggest in Kurdistan Province, Iran. Sanandaj is located in the Longitudes 46°24′00″ to 47°19′00″ East and Latitudes 35°3′00″ to 35°39′00″ North ( Figure 1). The city has a cold and semi-arid climate but moderate weather during spring. The maximum and the minimum temperature fluctuate between 44 °C and −13.5 °C. In this study area, vegetation coverages are mostly grass, while trees and high vegetation coverage are less marked. Therefore, the two study areas have different geomorphologies and are appreciated for this comparative study.  With a population of almost 432 thousand, the city of Sanandaj is the biggest in Kurdistan Province, Iran. Sanandaj is located in the Longitudes 46 • 24 00" to 47 • 19 00" East and Latitudes 35 • 3 00" to 35 • 39 00" North ( Figure 1). The city has a cold and semi-arid climate but moderate weather during spring. The maximum and the minimum temperature fluctuate between 44 • C and −13.5 • C. In this study area, vegetation coverages are mostly grass, while trees and high vegetation coverage are less marked. Therefore, the two study areas have different geomorphologies and are appreciated for this comparative study.

Image Selection for InSAR DEM Generation
Sentinel-1 satellite data are the first program by Copernicus, conducted by the ESA. Sentinel-1A and Sentinel-1B are two platforms for it [51]. This product has a 5.7 cm wavelength in C-band [12,52]. Crisis management, natural hazards and interferometry studies are among the most applications of the Sentinel-1 mission [53]. Interferometric Wide swath (IW), strip map (SM), extra-wide swath (EW) and wave (WV) are four sensor modes for the Sentinel-1 mission [51]. Single look complex (SLC), ground range detected (GRD), and ocean (OCN) are three types of data products for Sentinel-1 [12,52]. SAR interferometry needs the combination of two images from the same scene but at different times and sensor position [54]. It is worth mentioning that for Sentinel-1 if one wants to collect data with more PB, the TB will dramatically go higher. Table 1 shows the details of the data used in this study.

Pre-Processing and Processing of Data Used
InSAR technique is the comparison of the phase values for a pair of SAR images [49,55]. The phase information of two or more SAR imageries, acquired over the same region (same geographical position), either simultaneously or at different times, due to the scattered echoes' monochromatic nature, can be combined [3]. Figure 2 details the DEM creation's comprehensive methodology from Sentinel-1 satellite imagery and its comparison with the other products. Sentinel-1 satellite data are the first program by Copernicus, conducted by the ESA. Sentinel-1A and Sentinel-1B are two platforms for it [51]. This product has a 5.7 cm wavelength in C-band [12,52]. Crisis management, natural hazards and interferometry studies are among the most applications of the Sentinel-1 mission [53]. Interferometric Wide swath (IW), strip map (SM), extra-wide swath (EW) and wave (WV) are four sensor modes for the Sentinel-1 mission [51]. Single look complex (SLC), ground range detected (GRD), and ocean (OCN) are three types of data products for Sentinel-1 [12,52]. SAR interferometry needs the combination of two images from the same scene but at different times and sensor position [54]. It is worth mentioning that for Sentinel-1 if one wants to collect data with more PB, the TB will dramatically go higher. Table 1 shows the details of the data used in this study.

Pre-Processing and Processing of Data Used
InSAR technique is the comparison of the phase values for a pair of SAR images [49,55]. The phase information of two or more SAR imageries, acquired over the same region (same geographical position), either simultaneously or at different times, due to the scattered echoes' monochromatic nature, can be combined [3]. Figure 2 details the DEM creation's comprehensive methodology from Sentinel-1 satellite imagery and its comparison with the other products.  3.2.1. TOPS Split and Applying Orbit File SLC product has 3 Sub-Swaths of IW1, IW2, and IW3, which each sub swath is for an adjacent acquisition by the TOPS mode ( Figure 3). Sentinel-1 TOPS Split is used to select only those bursts which are needed for the analysis. Currently, each Sub-Swath should be processed separately, unless they Sensors 2020, 20, 7214 5 of 24 must be merged. Orbit auxiliary data consists of information related to the satellite's position during data acquisition [55]. This information should be added to the images' metadata using the application orbit file operator.
Sensors 2020, 18, x FOR PEER REVIEW 5 of 24 3.2.1. TOPS Split and Applying Orbit File SLC product has 3 Sub-Swaths of IW1, IW2, and IW3, which each sub swath is for an adjacent acquisition by the TOPS mode ( Figure 3). Sentinel-1 TOPS Split is used to select only those bursts which are needed for the analysis. Currently, each Sub-Swath should be processed separately, unless they must be merged. Orbit auxiliary data consists of information related to the satellite's position during data acquisition [55]. This information should be added to the images' metadata using the application orbit file operator.

Co-Registration and Enhanced Spectral Diversity
The back geocoding (co-registration) operator in SNAP software was used to co-register the two split satellite data based on the orbital information. In order to extract the phase difference, two acquisitions must be stacked first. This is an essential step in interferogram formation, as it ensures that each ground object has the same range and azimuth pixel in either image [8,46]. Enhanced spectral diversity (ESD) can be employed to increase the quality of the co-registration process. ESD corrects range and azimuth's shift and is only applied for more than one burst [55], which results in a better co-registered image.

Interferogram Formation and Coherence Estimation
An interferogram is calculated by cross multiplication of the master and slave images [57]. Here, the phase band shows the difference of a pixel in two images [39,40,46]. An interferogram ranging from -π to +π and is displayed in a rainbow color scale, which has variations from topography and surface deformation [57]. On the other hand, these patterns are called fringes and a full 2π cycle, where each cycle is demonstrating half of the sensor's wavelength [45,57].

DEM flat disp atm noise
Interferogram formation represents a raw DEM resolution Ω via an independent number of looks used in the processing (Equation (2)

Co-Registration and Enhanced Spectral Diversity
The back geocoding (co-registration) operator in SNAP software was used to co-register the two split satellite data based on the orbital information. In order to extract the phase difference, two acquisitions must be stacked first. This is an essential step in interferogram formation, as it ensures that each ground object has the same range and azimuth pixel in either image [8,46]. Enhanced spectral diversity (ESD) can be employed to increase the quality of the co-registration process. ESD corrects range and azimuth's shift and is only applied for more than one burst [55], which results in a better co-registered image.

Interferogram Formation and Coherence Estimation
An interferogram is calculated by cross multiplication of the master and slave images [57]. Here, the phase band shows the difference of a pixel in two images [39,40,46]. An interferogram ranging from -π to +π and is displayed in a rainbow color scale, which has variations from topography and surface deformation [57]. On the other hand, these patterns are called fringes and a full 2π cycle, where each cycle is demonstrating half of the sensor's wavelength [45,57].
The interferogram consists of phase variation φ from contributing factors, the Earth's curvature (phase of flat-earth) φ f lat, the topographic surface (phase of topographic) φDEM, conditions of the atmosphere (e.g., pressure change, temperature, and humidity between the two acquisitions) φatm, and noise φnoise (e.g., change of the scatters, volume scattering and different look angle), and finally the surface deformation φdisp (Equation (1)).
Interferogram formation represents a raw DEM resolution Ω r via an independent number of looks used in the processing (Equation (2)) [31]: where n az and n rg are azimuth and range independent number of looks, respectively, while δ gr az and δ gr rg are single SAR pixel azimuth and range ground resolution, respectively. Coherence is a useful source to know about the quality of an interferogram [57]. Coherence is formed as a separate raster, which shows how similar a pixel can be between the two images ranging from 0 to 1 [57]. Areas of high and low coherence will appear in a bright and dark color, respectively.

TOPS Debursting
TOPS data has 9 bursts, which is typically separated by some specific zones [51]. However, any values within these zones are defined as invalid data and must be debursted. To remove the lines between the bursts, the Sentinel-1 TOPS Debursting operator using SNAP software should be applied to the interferogram.

Phase Filtering and Multilooking
For a proper unwrapping process, phase filtering and multi-looking are two essential steps. Noise from geometric and temporal decorrelation and volume scattering can corrupt the interferometric phase [46,57]. Restoring phase information in decorrelated areas is quite hard, but applying proved phase filters (e.g., the Goldstein filter), could enhance the fringes' quality in the interferogram [57]. It is a preprocessing technique that considerably enhances the phase unwrapping accuracy and decreases the residues in the phase unwrapping step. The method was first proposed by Goldstein and Werner [58]. For Goldstein filter, adaptive filter exponent was 1 (its value lies in range of 0 and 1, which the larger the value the stronger the filtering process will be), for fast Fourier transformation (window length size of 64) and windows size (size of 3) were remained as default. It is worth mentioning that the coherence mask was not applied.

Phase Unwrapping
Unwrapping was applied in this study by using SNAPHU software [59]. Phase unwrapping (PU) is associated with the convergence of the phase difference from point to point by adding the integer number of cycles, which minimizes the phase difference [8]. PU recovers the integer number of cycles n to be applied to the wrapped phase ϕ. Therefore, the unambiguous phase value ψ can be eventually obtained for each pixel (Equation (3)): The interferometric phase is somewhat ambiguous and within the scale of 2π can be known only [57]. In order to relate the phase to the topographic height, it must be unwrapped first. PU solves the ambiguity through integrating phase differences among neighboring pixels [57].

Phase to Elevation Conversion and Terrain Correction (TC)
The unwrapped phase is not a metric measure yet. To convert the radian into the meter unit, the phase to elevation operator should be applied. The acquisition geometry for both images (master and slave), phase noise and possible unwrapping errors significantly affect this step [57].
TC will geocode the final bands by correcting geometric distortions of SAR images using a ready to download DEM and, lastly, produce a projected map. On the other hand, TC (geocoding) converts an image from either slant range (SR) or ground range (GR) into a coordinate system. Finally, this step using a DEM corrects geometric distortions of foreshortening, layover, and shadow.

Standard Errors of the Estimate
Generally, the accuracy assessment is vital in approving the quality of extracted information from remotely sensed data [60]. Regression analysis is associated with a collection of statistical techniques, where the behavior of random variables (dependent variables) using quantitative variables (independent variables) are described [61,62]. The correlation coefficient (R) has the values between +1 and −1, where 1 determines a total positive linear correlation, 0 represents no linear correlation, and −1 shows the total negative linear correlation [63,64]. The accuracy of predictions is measured by the standard error of the estimate and defined by the following equation [65,66]: where, σ est is the standard error of the estimate,Ŷ is estimate variables, Y is an actual variable, and N is the number of GCPs. The numerator is squared differences between the actual variables and the estimated one. As a whole, the lower the numerator, the higher the accuracy will be [65,66].

Hydrological Delineation
By delineating the stream networks for the created DEM from Sentinel-1 and the other products of AIRSAR, TanDEM-X, ALOS-PALSAR and SRTM in the GIS environment using the hydrology toolbox, the study was well-validated and justified as well. Based on the quality of the decorrelation of phase and spatial resolution for each DEM, many stream orders can be extracted in this method. The higher the spatial resolution, the more the stream orders can be delineated. Based on the system's structure, rivers of the order number "1" are the highest tributaries. Thus, if two streams of the outermost tributaries merge, it gives a higher order. Therefore, whenever two streams with different orders merge, it gives a higher-order number than them.
The DEMs were used to delineate streams one by one. For example, the created DEM from Sentinel-1 was imported into GIS, and a hydrology toolbox was employed for extracting stream orders for it. In the first step, the filling operator was utilized for filling any possible gap in the DEM. In the second step, the flow direction operator was applied for obtaining all directions of water flow in the DEM. At the third stage, using the flow accumulation raster threshold, those in the same direction or path were connected. Next, by using a stream link, the streams in the same directions were linked together. Then, a stream order operator was applied to order the streams from 1 to whatever. Finally, the extracted rivers were converted to feature using stream to feature operator. It is worth mentioning that these processes were applied for all the DEMs.

Validation Using Ground Control Points (GCPs)
The vertical accuracy of the DEMs was validated via comparing the interpolated elevations from each DEM and global navigation satellite system (GNSS) ground control points (GCPs). All DEMs of the both study areas were validated through available GCPs. The static GCPs (validation points) were measured using a GPS (geodetic receivers). Trimble device in Malaysia and Leica device in Iran were placed for minimum three hours until the signals were received from the satellites (based on the nearest continuously operating reference stations). The horizontal and vertical accuracy of each gcpare not exceeded ±2 cm and ±5 cm, respectively.
Available eight GCPs for Malaysia and six for Iran were used for validation. Normally, available GCPs are rare for example, Halim, et al. [67] used 45 GCPs (in 198,000 km 2 meaning 1 GCP per 4400 km 2 ) for a study performed in East of Malaysia. And in another study done in Perlis region in Malaysia, Pa'suya, et al. [68], used 38 GCPs (in 810 km 2 meaning 1 GCP per 21 km 2 ). The obtained elevations from DEMs compared with the corresponding elevation values of the GCPs. The outlier test using interquartile range (IQR) was done for the discrepancy between each DEM and the GCPs based on Equation 5 [69].
where, IQR is Interquartile range, Q1 represents the lower quartile, and Q3 is the upper quartile. Accordingly, the Root Mean Square Error (RMSE) for the discrepancy between each DEM and the GCPs were also computed using Equation 6 [70].
Where, X(measured) indicates the elevation of the GCPs, X(predicted), defines the elevation extracted from the DEMs, and n is the total number of GCPs.

DEM Creation from Sentinel-1 Using InSAR Technique
Using InSAR technique, for each study area, a DEM was extracted from Sentinel-1 data. Because of low coherence, vegetated areas typically will produce poor elevation value in C-band [44,71].
Our results also showed that DEM from L-band (ALOSPALSAR) is more suitable than DEM from C-band (Sentinel-1). Because of the inherent characteristics of ALOS PALSAR L-band data, the accuracy of its resulted DEM is normally better than the Sentinel-1 and TanDEM-X DEMs. However, in Iran's selected study area, mostly covered by rangelands, the DEM accuracy was also slightly better than Malaysia's study site, where forests are dominant. Figure 4 shows the created DEM from Sentinel-1 for both areas. Accordingly, the Root Mean Square Error (RMSE) for the discrepancy between each DEM and the GCPs were also computed using Equation 6 [70].
Where, indicates the elevation of the GCPs, , defines the elevation extracted from the DEMs, and n is the total number of GCPs.

DEM Creation from Sentinel-1 Using InSAR Technique
Using InSAR technique, for each study area, a DEM was extracted from Sentinel-1 data. Because of low coherence, vegetated areas typically will produce poor elevation value in C-band [44,71]. Our results also showed that DEM from L-band (ALOSPALSAR) is more suitable than DEM from C-band (Sentinel-1). Because of the inherent characteristics of ALOS PALSAR L-band data, the accuracy of its resulted DEM is normally better than the Sentinel-1 and TanDEM-X DEMs. However, in Iran's selected study area, mostly covered by rangelands, the DEM accuracy was also slightly better than Malaysia's study site, where forests are dominant. Figure 4 shows the created DEM from Sentinel-1 for both areas.

The Linear Regression and Standard Errors of the Estimate
DEMs vary in quality and pixel sizes depending on the creation methods and data used [72]. The linear regression as one of the reliable ways to validate raster layers [73] was used to evaluate the relative accuracy of the created DEM. The extracted DEM from Sentinel-1 imagery was 14-m cell

The Linear Regression and Standard Errors of the Estimate
DEMs vary in quality and pixel sizes depending on the creation methods and data used [72]. The linear regression as one of the reliable ways to validate raster layers [73] was used to evaluate the relative accuracy of the created DEM. The extracted DEM from Sentinel-1 imagery was 14-m cell size.
Sensors 2020, 20, 7214 9 of 24 Therefore, a 12.5 × 12.5 spatial resolution DEM of the ALOS-PALSAR products was selected as the reference for validating using linear regression.
The linear regression and standard errors of the estimate were employed to validate the extracted DEM's relative accuracy. For this reason, 8120 and 5334 systematic points for Malaysia and Iran were selected, respectively ( Figure 5) were created and the elevation values from DEMs were exported into them. These points were employed to validate (relative accuracy) the DEMs using the linear regression (Table 2) and scatter plots ( Figure 6) for the both areas. The correlation coefficient was 99% and 100% for 2 study areas of Malaysia and Iran, respectively. This fact confirms that the variables were fully correlated (positive linear correlation). The estimated standard error was quite high for either study area, but for the Cameron Highlands, this was more noticeable (almost 7 m); because on one hand, vegetation coverages mostly covers the study area and C-band SAR data cannot penetrate dense vegetation coverages, on the other hand, the short PB can be another reason. As it can be seen from the scatter plot, for the study area of Iran, where is mostly covered by spare vegetation coverages, the values are more correlated than the denser vegetated area of Malaysia.
Sensors 2020, 18, x FOR PEER REVIEW 9 of 24 size. Therefore, a 12.5 × 12.5 spatial resolution DEM of the ALOS-PALSAR products was selected as the reference for validating using linear regression. The linear regression and standard errors of the estimate were employed to validate the extracted DEM's relative accuracy. For this reason, 8120 and 5334 systematic points for Malaysia and Iran were selected, respectively ( Figure 5) were created and the elevation values from DEMs were exported into them. These points were employed to validate (relative accuracy) the DEMs using the linear regression (Table 2) and scatter plots ( Figure 6) for the both areas. The correlation coefficient was 99% and 100% for 2 study areas of Malaysia and Iran, respectively. This fact confirms that the variables were fully correlated (positive linear correlation). The estimated standard error was quite high for either study area, but for the Cameron Highlands, this was more noticeable (almost 7 m); because on one hand, vegetation coverages mostly covers the study area and C-band SAR data cannot penetrate dense vegetation coverages, on the other hand, the short PB can be another reason. As it can be seen from the scatter plot, for the study area of Iran, where is mostly covered by spare vegetation coverages, the values are more correlated than the denser vegetated area of Malaysia.

Validation Using GCPs
To assess the results' quality and validity, several available GCPs points were collected to be compared with all DEMs including the created DEM from Sentinel-1 SAR satellite imagery. Figure 7 shows eight GCPs for Malaysia and six GCPs for Iran over land cover type map. The results showed differences between them, meaning that based on the GCP points, the DEM created from Sentinel-1 could not be reliable. Therefore, the researchers should be cautious about using it as a base DEM for their studies. The outlier test results show that there are no outlier concerning the obtaining elevation values from the DEMs for both study areas (Table 3).

Validation Using GCPs
To assess the results' quality and validity, several available GCPs points were collected to be compared with all DEMs including the created DEM from Sentinel-1 SAR satellite imagery. Figure 7 shows eight GCPs for Malaysia and six GCPs for Iran over land cover type map.

Validation Using GCPs
To assess the results' quality and validity, several available GCPs points were collected to be compared with all DEMs including the created DEM from Sentinel-1 SAR satellite imagery. Figure 7 shows eight GCPs for Malaysia and six GCPs for Iran over land cover type map. The results showed differences between them, meaning that based on the GCP points, the DEM created from Sentinel-1 could not be reliable. Therefore, the researchers should be cautious about using it as a base DEM for their studies. The outlier test results show that there are no outlier concerning the obtaining elevation values from the DEMs for both study areas (Table 3). The results showed differences between them, meaning that based on the GCP points, the DEM created from Sentinel-1 could not be reliable. Therefore, the researchers should be cautious about using it as a base DEM for their studies. The outlier test results show that there are no outlier concerning the obtaining elevation values from the DEMs for both study areas (Table 3). In fact, establishing GCPs is a time-consuming, hard job, costly, and even impossible sometimes [74]. The GCPs are essential for validating a DEM based on their actual elevations using Root Mean Square Error (RMSE) [72,75]. The validation technique using GCPs was earlier employed by the United States and Japanese researchers for the accuracy assessment of the Advanced Space-borne Thermal Emission and Reflection Radiometer (ASTER) [76]. Results indicate that the ALOSPALSAR for Iran and AIRSAR for Malaysia were the best DEMs with total RMSE of ±5.2 m and ±6.4 m, respectively (Tables 4 and 5).  Accuracy assessment using the available GCPs was performed for the created DEM from Sentinel-1 based on the land cover types ( Table 6). The interferometric measurements is normally affected by land covers [77]. To assess the impacts of different land covers on the accuracy of the created DEM using Sentinel-1, the GCP points were grouped based on the corresponding land cover types. In the study area of Malaysia, GCPs number 2, 4, and 8 were located in vegetation and florification, GCPs number 1, 3, 5, and 7 were lies on forests and GCP number 6 was located on township. At the same time, in the study area of Iran, GCPs number 1, 5, and 6 were placed in rangeland, and GCPs number 2, 3, and 4 were located in township, agriculture and dryfarming, respectively. Overall, RMSE values are quite high and ranges from five to thirteen meters in the study area of Malaysia, were mostly covered by dense forest and vegetation. In Malaysia, the highest RMSE of 13 m can be seen in township (this maybe because there are high rising hotels in the Cameron Highlands with sparse trees coverages). Then, the highest value was calculated for forests of 11.7 m, followed by vegetation and florification with 5.6 m. Because of low vegetation coverages these values were quite lower in Iran than Malaysia, but because of the technical issues the RMSE for Iran also showing that Sentinel-1 failed to create an accurate DEM. In Iran, sparse and bare lands are contributed to lower RMSE meaning that agriculture, dryfarming, and rangeland showed lower RMSE of 7, 5, and 3.5 m, respectively. The township has RMSE of 9 m.

Hydrological Analysis
Hydrological units are used for the planning and management of natural resources. To better understand the low accuracy of the extracted DEM, the hydrological analysis was performed as well. Therefore, the created DEMs were hydrologically compared with AIRSAR, ALOS-PALSAR. TanDEM-X, and SRTM products. The hydrological networks were delineated for all four DEMs using the hydrology toolbox in the GIS environment. The concern with the extracting stream networks, given the same accumulation raster layer, the higher threshold value contributes to less dense stream orders and vice versa [78]. Figures 8 and 9 clearly show the differences between the extracted DEM from Sentinel-1 data with the other DEMs for both study areas. Because of the close spatial resolution and the same accumulation threshold in GIS, the equal number of stream orders was delineated for extracted DEM from Sentinel-1 data and the ALOS-PALSAR in either study area. Moreover, because of the higher spatial resolution, AIRSAR DEM contributes to the highest stream orders, While, regarding the lower spatial resolution than the other, DEM of SRTM has only five stream orders in the Cameron Highlands. Regarding Iran's study area, because of the lower spatial resolution, the lower stream orders were extracted from TanDEM-X and SRTM. The wavelength, PB and TB have an essential role in the accuracy of a DEM. More importantly, because of poor coherence and interferogram resulted from the low PB and short wavelength, the extracted DEMs from Sentinel-1 (C-band) showed too many noises in-stream orders, exceptionally orders number one (Figures 8a and 9a). Simultaneously, such noises cannot be seen in the other DEMs, even DEM of SRTM and TanDEM-X with the lower cell size.

Limitations of InSAR Technique
The nature of the InSAR technique is noisy and this approach is limited to coarse scene descriptions. At the same time, in built-up areas it result in layover and shadow [79]. InSAR limits depends on the spatial extent and magnitude of features, which depends upon geometric attributes of data used (L-band or C-band) [79,80]. Although it is technically possible to create an interferogram, InSAR is not an ideal technique for dense vegetated areas [80], especially when the images are also geometrically unusable for such areas. To obtain better results from InSAR technique, it is better to use the longer wavelength images such as L-band and more importantly using it in less vegetated areas.

Multi-Temporal InSAR and PSInSAR Approaches in DEM Generation over Vegetated Areas
For such areas like the Cameron Highlands (Malaysia), where mostly covered by dense vegetation coverages, multi-temporal InSAR technique, which integrates different estimators, resulting in a self-adaptive and reliable data processing [81,82]. Another advantage of this model is that, it is self-adaptive and the algorithm works well under various situations (more accurate coherence estimation); therefore, it offers more promising results than InSAR for surface deformation studies and DEM creation [81,83]. Another technique that can be used for DEM generation in vegetated areas is the Permanent scatterer (PSInSAR), which is different completely from the traditional InSAR that firstly developed by Ferretti, et al. [84]. The results of this model is by far more accurate than InSAR [85,86]. For reliable results, unlike InSAR, this technique requires many interferometric data pairs (more than 20, the more the data the more the result will be accurate) [86]. The PSInSAR calculates the motion of scatterers and the propagation delay due to changes [86,87]. This technique has been proved effective in many surface deformation-monitoring studies [87].

Atmospheric Delays
Short temporal baseline, suitable perpendicular baseline, and suitable atmospheric conditions for data acquisition are very essential and should be considered [46,88,89]. In this study, all above conditions were considered carefully. For example, water vapor in the atmosphere result in phase delays and decreases the quality of the measurement [46,88]. Therefore, it is advisable to collect images obtained during dry days and also during no rainfall days. Hanssen [89], suggested an alternative for this and mentioned that try to collect nighttime's data, which are normally less affected by atmospheric conditions. Due to the good temporal characteristics of the data used in this study, we did not apply atmospheric delays correction. Date of the data used for the study area (Malaysia) are in months of late February and early March. In addition, for the study area (Iran) the data are in months of late July and mid-August. Refer to the weather condition in the study areas, we can find out that the images were not much affected by hot days (water vapor) and rainfall (wet days). Because there are two main monsoon seasons in Malaysia from late May to September and from October to February [90]. In Iran also the data were collected during the times that neither there is rainfall nor very hot days [91].
Although, based on the above-mentioned statements, this paper did not used any method for removing atmospheric path delays, but suggests a few models to be used in the future studies (especially, where researchers cannot collect suitable data), including direct observations of atmospheric properties (e.g., continuous GPS observations), empirical corrections [92], weather models [93], global atmospheric models [94], ERA5 global atmospheric model [95], and MERIS water vapor data and elevation-dependent interpolation model [96]. Atmospheric delays can be classified as neutral and ionospheric components, which ionospheric component is associated with the free electrons in the atmosphere [92,97]. The neutral atmospheric component is usually caused by a combination of non-dipole (hydrostatic) and dipole (water vapor) components [92,97]. The non-dipole delay dominates with delays of almost 2 m, while the dipole is associated with approximately a few tens of centimeters [97].

Effects of Different Wavelengths on DEM Products
The X-, C-and L-band SAR signals have a different interaction with surface objects. Typically, X-band signals backscatter from the top of the trees canopy, while C-band signals penetrate deeper into the canopy [98]. In contrast, L-band signals provide a greater penetration rate than X-and C-bands [71,[99][100][101]. The estimated interferometric phase and the coherence between two images (master and slave) is generally an indicator for assessing the phase information quality, which have direct influences on DEM generation. The coherence value ranges from 0-1. The more the value the better the quality of the interferogram will be [44,102]. The different wavelengths have a different interaction with different land covers [44,102,103]. For example, the coherence coefficient calculated from L-band ALOSPALSAR is normally higher than C-band Sentinel-1 and X-band TanDEM-X in vegetated areas [44]. When SAR signals could not penetrate vegetation, then coherence will be poor (dark) [104]. Therefore, coherence loss causes poor interferometric results [45,71], which result in a low quality DEM. Short wavelength SAR sensors of C-and X-bands (because of low penetration rate) have poor coherence (dark) in vegetation areas (this differs based on the density of vegetation); while having a slightly higher coherence (bright) in township areas. Compared to the C-and X-bands that become uncorrelated quickly, because of a better coherence, L-band interferometry is more suitable for denser vegetated areas [98,103]. The interferogram is generally shown in a color scale that ranges from -π to +π (also, called "fringes" and it is displayed as a set of arbitrary colors' cycles) [45]. The dense interferometric fringe is generated by X-band, while the relatively sparse fringes are derived from Sentinel-1 and ALOSPALSAR in order [44]. To create a high quality DEM, these fringes should be present throughout the image [45]. Like the wavelength, the land cover also affects the DEM quality [41,44]. Accordingly, phase information (interferogram) over vegetated areas, which have low coherence will not produce accurate elevation measures. For the denser vegetated coverage areas, it is better to use L-band satellite data to create a DEM rather than C and X-bands [67,103].

Co-Registration of Different DEMs Using the Least Trimmed Squares (LTS) Estimator
Different DEMs only with common coordinate systems can be compared [105][106][107]. In this study, the least trimmed squares (LTS) estimator was used to co-register the DEMs, which was also applied by Zhang and Cen [105]. The model is a simple but high-performance technique to co-register the different DEMs, which was initially introduced by [108]. The algorithm does not need any prior information about the terrain changes [105]. In addition, identifying the terrain changes with a self-adaptive threshold and very high matching and change detection accuracy are among advantages of the model [105].
The optimum perpendicular baseline (PB) for DEM creation should be between 150 to 300 m [45,46]. The higher and the lower PB can cause loss of coherence and leading to decorrelation (due to high sensitivity to phase noise and atmospheric effects). Simultaneously, the temporal baseline (TB) should be short and usually 6 to 12 days at best [46]. The Sentinel-1 mission was basically designed for the deformation studies not for DEM generation [100]. Therefore, PB above 100 m are hard to find [45]. Most of the image pairs have a baseline below 30 m and pairs with short TB and large PB can be hard to find [45,100]. In fact, this is the main reason why it cannot be applied for DEM generation. However, because of the abovementioned statements and knowing that Sentinel-1 is a C-band SAR data, it cannot be employed for DEM generation in vegetated areas. Accordingly, because of the low PB, creating DEM from Sentinel-1 over bare lands will also produce bad results.
In terms of data acquisition and generation processes, creating a DEM using SAR data is a time-consuming issue [109]. It is worth mentioning that if a researcher can collect data with the optimum PB and TB; the best spatial resolution for Sentinel-1 will be around 13 m [110]. The quality of the DEM generation is strongly affected by atmospheric disturbances and coherence. From Sentinel-1, a PB higher than 100 m is hard to find. Therefore, creating the DEM is a challenging issue. Despite the higher spatial resolution, the short PB and wavelength result in errors, so the created DEM has not better quality than of freely available DEMs, including SRTM.
Based on findings using evaluation techniques of the linear regression, ground control points and hydrological models, the created DEM has low accuracy. The linear regression showed more error for dense vegetated areas. Moreover, the hydrological networks for the created DEM, showed bad results, especially in stream orders number one. Additionally, results by GCPs were also confirmed that creating DEM from Sentinel-1 data is not giving promising result, especially over the vegetated areas.

Conclusions
It is noted that DEM generation due to numerous usages for extracting different attributes, such as slope, aspect, elevation, curvature, etc., is one of the essential spatial information tools used in the earth sciences. Therefore, optical and radar satellite images are the two main sources to produce DEM. DEM generation with the SAR interferometry approach is one of the most common methods and It is very interesting for geosciences researchers. Furthermore, one of the most recent Sentinel-1 missions is SAR imagery. Sentinal-1 images are among the best earth observations for DEM generation due to acquire images from the whole earth surface once every six-day, availabilities, high-resolution, and C-band sensitivity to the surface structure.
The objective of this study was to produce a DEM from Sentinel-1 satellite imagery using the InSAR technique along with to compare it with DEMs of AIRSAR, ALOS-PALSAR, TanDEM-X, and SRTM. For this reason, Sentinel-1 products for two different study areas were acquired, from which DEM with 14-m pixel size was obtained for a part of the Cameron Highlands, Pahang, Malaysia, and a part of Sanandaj, Iran. Based on the analysis and comparisons with the other proven DEMs using Linear regression and hydrological models, the quality of derived DEM for the two study areas was less than the other DEMs, even 90-and 30-m cell size DEM of TanDEM-X and SRTM, respectively. At the same time, results of RMSE using a few available GCPs showed that the created DEMs were not giving promising results. The validation results showed that Sentinel-1 could not produce accurate DEMs for two reasons. Firstly, the C-band signals cannot penetrate the denser vegetation coverage, especially canopies and trees. Secondly, the relatively low accuracy can be due to the small perpendicular baseline of Sentinel-1 data, which is rarely longer than 100 m could be found. For creating a DEM, the baseline should generally be between 150 to 300 m. This study recommends researchers to use freely available DEMs (e.g., SRTM) instead of creating DEM from Sentinel-1. Finally, if a DEM with a higher spatial resolution is required, they can use other products such as ALOS-PALSAR and LiDAR (Light Detection and Ranging) data, which is a technique to measure distances by target illuminating using a laser and calculating the backscatters with a sensor.