DEM-Based Vs30 Map and Terrain Surface Classification in Nationwide Scale—A Case Study in Iran

Abstract

Different methods have been proposed to create seismic site condition maps. Ground-based methods are time-consuming in many places and require a lot of manual work. One method suggests topographic data as a proxy for seismic site condition of large areas. In this study, we mainly focused on the use of an ASTER 1c digital elevation model (DEM) to produce Vs30 maps throughout Iran using a GIS-based regression analysis of Vs30 measurements at 514 seismic stations. These maps were found to be comparable with those that were previously created from SRTM 30c data. The Vs30 results from ASTER 1c estimated the higher velocities better than those from SRTM 30c. In addition, a combination of ASTER 1c and SRTM 30c amplification maps can be useful for the detection of geological and geomorphological units. We also classified the terrain surface of six seismotectonic regions in Iran into 16 classes, considering three important criteria (slope, convexity and texture) to extract more information about the location and morphological characteristics of the stations. The results show that 98% of the stations are situated in six classes, 30% of which are in class 12, 27% in class 6, 17% in class 9, 16% in class 3, 4% in class 3and the rest of the stations are located in other classes.

1. Introduction

Recent developments in urban environments and settlement areas, as well as their complexity has led to increased demand for site condition maps [1,2,3]. In order to connect the effective parameters related with the ground motion, coordinates and geospatial analysis play an important role [4,5,6,7]. In fact, coordinates establish a unique opportunity for us to understand our surrounding areas better than we have been able to in the past. Generally, Geographical Information System (GIS) is a technology which offers potential for spatial analysis. Using either local or global coordinate systems, it connects effective but somehow discrete factors, such as Vs30 (average shear wave velocity between a 0 and 30 meter depth) and elevation data. GIS calculations based on Vs30 have been used frequently for geotechnical investigations and geomorphologic classifications on national scales [8,9,10,11].

A review of the research literature indicates that the effects of site amplifications (stratigraphy + topography) on seismic ground motion and consequently, their potential for earthquake damages, have been studied in previous work [12,13,14,15,16,17,18]. Vs30 is now one of the standard indicators for mapping seismic site conditions in most building codes of earthquake-prone countries. For instance, similarly to NEHRP (National Earthquake Hazards Program) classifications in the U.S., Iranian standard (Standard No. 2800) under the authority of the BHRC (Building & Housing Research Center) considers soil classifications as a basis for estimating Vs30 [19,20]. However, the quality and density of Vs30 measurements vary from one region to another. Additionally, sometimes, geological information of the ground surface is lacking, proves to be unavailable, or is unreliable. In this situation, continuous Vs30 measurements (cell-by-cell) from the satellite-based elevation data gain importance and become an ongoing topic in order to produce the site condition maps related with prompt seismic damage assessments. A glimpse of site classification maps for vast areas (e.g., on a national scale), shows that each of the models use specific formulation approaches and can be categorized into two main classes: (1) hybrid models which use geomorphologic, geologic and topographic data, and (2) single models which use only one item. A principal example for the first category is a nationwide site classification map of Japan developed by Matsuoka et al. [5], which made use of a multiple regression analysis to calculate the regression coefficients of 21 geomorphologic units. The Japanese QuiQuake system (a Quick estimation system for earthquake maps triggered by observation records) uses the method proposed by Matsuoka et al. [5] to provide amplification capability maps throughout Japan [21]. The site condition map of the continental United States by Allen and Wald [16] (updated by Allen and Wald [17]) is an example of the second category which utilized only topographic slope data from the Shuttle Radar Topography Mission (SRTM) 30 arcsec (~ 1km resolution near the equator) as a proxy for the site amplification. The Prompt Assessment of the Global Earthquake for Response (PAGER) used the method described by Wald and Allen [15] as a core idea for real-time damage assessments (http://earthquake.usgs.gov/eqcenter/pager/S). In recent studies for nationwide purposes, Forte et al. [22] provided seismic soil classification maps for Italy based on site-specific measurements, including Vs30 measurements and geological information. Stewart et al. [23] produced geologic and terrain-based site condition maps for California using geological information and 1, 3, 9 and 30 arcsec DEMs. The results produced by three arcsec datapoints were more accurately relative to the real measurements than other DEMs. Later, six proxy-based models for California (using geological information, topography, and geomorphological classification) were compared with the 503 Vs30 measurements and NEHRP site classes [24].

In this study, we used digital elevation data (ASTER 1 arcsec), gathered by a ASTER sensor (Advanced Spaceborne Thermal Emission and Reflection Radiometer), to understand how DEM resolutions from different sources can affect DEM-based Vs30 maps. Based on Vs30 measurements of 514 locations with strong motion seismometers in Iran, linear regression formulas for different slope classes were developed within six seismotectonic regions of Iran. Recognition of the small features can be improved in medium-resolution topographic gradient, but the derived Vs30 maps do not show significant improvements in comparison with those from low-resolution DEM.

2. Study Area

2.1. Seismotectonic Setting

Iran is trapped between the Arabian and Eurasian plates, with approximately 22 mm/yr of northward shortening [25]. The convergence of these plates shows the emergence of faults and different seismotectonic regions in Iran (Figure 1). Different seismotectonic regions have been identified for Iran, typically from the geographical and historical distributions of active faults and structural trends [26,27,28,29]. The main seismotectonic regions in Iran are as follows: (1) Alborz-Azarbaijan, (2) Central Iran, (3) Kopet Dagh, (4) Lut, (5) Makran, and (6) Zagros. Alborz-Azarbaijan and Zagros are the most populated blocks which absorb the main portion of continental shortening and consequently, the number of earthquakes in these blocks is considerably higher than elsewhere [25]. From the Bayesian estimation, the seismic hazard is lowest in the Lut and Central Iran blocks, meaning that these blocks are craton [30]. Slope classes in active tectonic regions are different from craton parts. Ideally, in the regression analyses, active and stable areas should be separated during Vs30 calculations, but due to large investments on geotechnical and seismological investigations in highly populated regions, a significant portion of our Vs30 measurements are from the Alborz-Azarbaijan and Zagros regions. On the other hand, the large deserts of the Middle East (Dasht-e Lut and Dasht-e Kavir) are situated in the Central Iran and Lut blocks. The density of the population in these blocks is sparse and in some places, close to zero. Since the main goal of nationwide Vs30 maps is prompt damage estimation, we considered all Iranian territory as an active region in Section 3.

2.2. Iranian Strong Motion Network (ISMN)

In 1981, the BHRC started to run the network with 261 accelerometers and today, there are over 1000 digital accelerometers recording earthquakes [31]. Although the ISMN has more than 1000 stations, some of these stations are no longer active. Accordingly, measured Vs30 data have been obtained only from 514 ISMN stations (Figure 1). As mentioned, the majority of stations are constructed in the heavily populated active blocks. The number of stations for Alborz-Azarbaijan, Central Iran, and Zagros blocks is 106, 128 and 190, respectively and the remaining 90 stations belong to the Lut (54), Makran (2), and Kopet Dagh (34) blocks. Except for marginal parts of the Central Iran and Lut blocks, there is no station in these blocks due to lack of population (Figure 1).

3. Data

3.1. Vs30 Measurements

Most of the Vs30 data (347 stations out of 514 stations) represent relatively low gradients of less than about 6% ($grade = \left(\frac{rise}{run}\right)\times 100$), or a slope of about 3 degrees (angle = arctangent (0.06)). However, in most samples, Vs30 data within ISMN were measured at hard-rock sites rather than the soil sites, implying that slope and Vs30 do not always have a straightforward relationship and must be formulated by applying some constraints, and assigning starting-and-ending values where the data is fewer. Even quite different performance can be expected from rock site measurements [32]. Holzer et al. [14] produced a Vs30 map of San Francisco Bay (140 km²) using Vs30 data of 210 stations, which was derived from (1) the average shear wave velocity values for surficial geologic units, (2) the thickness of each layer, and (3) NEHRP classifications (

Table 1. Topographic slope ranges for NEHRP classes adopted by Allen and Wald [16].

NEHRP Site Classification	Vs30 Range (m/s)	Slope Range (m/m)		General Description
		Active Tectonic Regions	Craton Regions
$E$	$<180$	$<3\times {10}^{-4}$	$<1\times {10}^{-4}$	Soil profile with soft clay
$D$	$180-240$	$3\times {10}^{-4}-3.5\times {10}^{-3}$	$1\times {10}^{-4}-8.5\times {10}^{-3}$	Stiff soil
$D$	$240-300$	$3.5\times {10}^{-3}-0.010$	$4.5\times {10}^{-3}-8.5\times {10}^{-3}$	Stiff soil
$D$	$300-360$	$0.010-0.024$	$8.5\times {10}^{-3}-0.013$	Stiff soil
$C$	$360-490$	$0.024-0.08$	$0.013-0.022$	Very dense soil and soft rock
$C$	$490-620$	$0.08-0.14$	$0.022-0.03$	Very dense soil and soft rock
$C$	$620-760$	$0.14-0.20$	$0.03-0.40$	Very dense soil and soft rock
$B$	$>760$	$>0.20$	$>0.40$	Rock & Hard rock

). However, it must be considered that in large countries such as Iran (~1.6 million km²), it is difficult to use this approach because it might require a huge pool of Vs30 measurements [33]. The obtained 514 samples do not seem enough to show the significance of the DEM-based Vs30 maps—there is still significant scatter. The regression trend is considered sufficient to be used between gradient of topography and Vs30. Figure 2 shows the statistical results of Vs30 measurements in different NEHRP classes. The number of stations in category B with higher velocities (>760 m/s) is considerably higher than the other classes. Moreover, Tthe standard deviation of this category is high, but in the other classes, it becomes smaller by reducing Vs30.

3.2. SRTM 30c and ASTER 1c

For topography, we considered two types of either low-resolution (SRTM 30c) or medium-resolution (ASTER 1c) elevation data. The SRTM contains the following two radars: (1) the C-band ($\lambda = 5.6 cm$) Imaging Radar, and (2) the X-Band ($\lambda = 3.0 cm$) Synthetic Aperture Radar sensors, indicating the elevation values between DSM (Digital Surface Model) and DTM (Digital Terrain Model). They were used on board the space shuttle for gathering data about the Earth’s environment with different spatial accuracies (30, 90, and 1000 m). Wald and Allen [15] used SRTM 30c for rather global-purpose Vs30 maps. We used SRTM 30c for the country, freely available by USGS (http://dds.cr.usgs.gov/srtm/version2_1/SRTM30//srtm30). The ASTER sensor gathers data in 14 spectral bands. ASTER 1arcsec is composed of 22,600 tiles, and in this research, over 300 ASTER GeoTIFF tiles (1°×1°) were obtained and merged (http://gdem.ersdac.jspacesystems.or.jp). The key difference between the two datasets is that although ASTER 1c benefits from a stereo imagery technique which gives it the unique ability to create digital elevation models with better accuracy, but the elevations from ASTER sensor constitute a mostly surficial model, whereas SRTM 30c is an interferometry-based elevation model (Allen and Wald, 2007 [15]). Neither SRTM 30c nor ASTER 1c will completely satisfy the topographic gradient of the near surface geology, rather, a bare-ground elevation model will better suit our application. However, owing to the semiarid and arid climate and short vegetation of the study area, the SRTM 30c and ASTER 1c are ideal for slope classification [34]. Note that different resolution topographic data result in different slope values; therefore, some pixel-based resizing can be applied before Vs30 comparison.

4. Vs30 Proxy Map

4.1. Regression Analysis and Vs30 Estimation

As the first step of research methodology, topographic attributes (i.e., elevation and slope of topography) were derived and compared with Vs30 data before applying the relevant regression analysis. Every pixel in the output raster (slope map) has a grade slope value ($grade = \left(\frac{rise}{run}\right)\times 100$) as a unitless ratio (meter per meter). The lowest slope value represents the flattest terrain and conversely, the highest one represents the steepest ground. The existence of a linear correlation between real Vs30 measurements and topographic attributes at the location of each station is the main concept before the regression analysis, irrespective of whether topographic attributes are useful in the creation of the site condition map or not. Here, the correlation of the variables was examined in both ordinary and logarithmic manners (Figure 3). As shown in Figure 3, the logarithmic correlation between Vs30 and the topographic data is better reflected than the simple correlation between them. The correlation between the height values and Vs30 measurements is too small in comparison with the slope values and the Vs30 measurements (Figure 3a). The slope of the topography describes the relief of the land using a percentage of change over a specified distance. After examining the ASTER 1c data, we noticed that some flat areas have height values near zero. Since the log (0) is not defined, we manually assigned a 1 m pixel value for flat areas. Figure 3b,d manifests that higher Vs30 measurements are in good agreement with the higher slope positions, which is likely related with the competent materials in such higher positions.

Three assumptions are needed before the regression analysis: (1) A logarithmic correlation between the height and the slope with Vs30 measurements should exist. This is essential to perform an independent or joint regression analysis to yield an acceptable outcome. Here, we are only considering slope values since elevation values near zero could be problematic in several ways. (2) We assign a unique Vs30 value (199.06 m/s) for the slope values near to zero in Equation (1). (3) the Vs30 of 600 m/s is assigned to water-covered areas, but inland water bodies can be masked as an alternative option. The regression model of Vs30 and slope is as follow:

$$\log \left(vs30\right) = a\times (\log \left(Slope\right))+b$$

(1)

where $a$ and $b$ are the coefficients of regression analysis for NEHRP classification. For each NEHRP class, the corresponding regression coefficient obtained from a different slope range is given in

Table 2. Corresponding coefficients for active regions compatible with NEHRP slope classes [37].

Slope Range (m/m)	Vs30 Range (m/s)	$\mathbf{log}\left(\mathit{v}\mathit{s}30\right) = \mathit{a}\times (\mathbf{log}\left(\mathit{S}\mathit{l}\mathit{o}\mathit{p}\mathit{e}\right))+\mathit{b}$		${\mathit{R}}^2$
		$\mathit{a}$	$\mathit{b}$
$<3\times {10}^{-4}$	$<180$	$0.46$	$2.229$	$0.136$
$3\times {10}^{-4}-3.5\times {10}^{-3}$	$180-240$	$-0.18$	$2.293$	$0.008$
$3.5\times {10}^{-3}-0.010$	$240-300$	$0.12$	$2.443$	$0.025$
$0.010-0.024$	$300-360$	$0.12$	$2.531$	$0.035$
$0.024-0.08$	$360-490$	$0.001$	$2.638$	$0.000$
$0.08-0.14$	$490-620$	$0.003$	$2.747$	$0.002$
$0.14-0.20$	$620-760$	$0.002$	$2.841$	$0.001$
$>0.20$	$>760$	$0.006$	$3.001$	$0.000$

, which shows that the coefficient of determination for the areas with a gentle slope is better estimated. This indicates that the database of Wald and Allen [15] was formed from mainly flat sites with relatively lower Vs30 values. Differently from Wald and Allen [15] and Allen and Wald [16], most of the Vs30 measurements in this study are more than 800 m/s, demonstrating that the seismic sites are probably chosen in rock sites (Figure 2).

4.2. Vs30 and Amplifiation Map

Raster calculation flow allows us to execute mathematical calculations, set up selection queries, and build a model. Once the relationship between regression coefficients and slope ranges is established, the Vs30 map of the country can be estimated [35]. The unitless slope values were assigned to each pixel according to the NEHRP classification (Table 2). As soon as we choose a specific range of slope values in the Raster Calculator tool, all pixels will be rated 0 (unacceptable) or 1 (acceptable), which means that the original values will be lost in the new binary map. One way to recover original pixel value for each class is to multiply the binary map and original map again. Using this method, we were able to recover the real values estimated from ASTER 1c and SRTM 30c slope classes while 0 values remained for void pixels.

Wald et al. [36], following Wald and Allen [15] more fully, described a methodology of a global-based Vs30 grid map from SRTM 30c with larger cells with dimensions of about $0.008\times 0.008$. But the size of each grid of Vs30 values from ASTER 1c was $0.001\times 0.001$. In order to compare new and previous results, their grid sizes must be equal. Thus, we applied the raster generalization to either clean up small erroneous cells or smooth out unnecessary details. We produced 8 discrete classes of the Vs30 map, which should be mosaicked. This procedure tends to produce a very large raster dataset; however, the generalized map along with mosaicking, allow us to take all the raster datasets and combine them into a single, seamless raster dataset between 8 adjacent raster datasets [35]. Nevertheless, there are some overlaps of the raster dataset edges that are being mosaicked together. These overlapping areas have been solved using a weight-based algorithm. Figure 4 and Figure 5 show Vs30 maps deduced from the log slope and the regression coefficients (Table 2).

Next, the amplification map of the study area was calculated using the method suggested by Borcherdt [38]:

$$F_v={\left(\frac{\overline{v}}{v}\right)}^{m_v}$$

(2)

where $F_v$ is the amplification factor, $\overline{v}$ represents the average value of the shear wave velocity for different sites, $v$ is an individual shear wave velocity in each cell or pixel, and $m_v$ is the constant for the corresponding ground motion. Here, Equation (2) has the flexibility to adopt different methods to yield both input ground motion and amplification factors. We calculated the amplification map for midperiods of ground motion by assuming a uniform peak ground acceleration. We took the ratio of the two maps (i.e., Vs30_30c and Vs30_1c) to examine the relative difference in amplification of specific parts of NW Iran (Figure 6). The calculation flow for the Vs30 and amplification ratio maps derived from slope classes is summarized in Figure 7.

4.3. Vs30 Results

The results from the regression analysis and raster calculation indicate that the estimated Vs30 values from SRTM 30c and ASTER 1c can be compared with the ground measurements. At first glimpse, the Vs30 measurements in Iran correlate better with 30c than 1c gradients (Figure 8). The reason for the better correlation between Vs30 from the coarser topography and real Vs30 values is complex, but it seems that it is because of nongeomorphic features. The scatter plots of Vs30 versus Vs30_1c and Vs30 versus log (slope_1c) do not show significant trends (Figure 8a,c). Figure 8 summarizes the trends that we could find from SRTM 30c and ASTER 1c. The overall trend indicates that higher velocities from either ASTER 1c or SRTM 30c are underestimated, mainly because most of the Vs30 data used by Wald and Allen [15] represent relatively low gradients of less than about 7%. Overall, the higher-resolution datasets had better estimation relative to SRTM 30c data. Figure 9a represents the obtained histograms of ASTER 1c for comparison with those produced from SRTM 30c (Figure 9b). The computed logarithmic standard deviation from ASTER 1c (0.23) was larger than that calculated from SRTM 30c (0.20), but the logarithmic median of the first one was smaller than the second, which was a better correlation between the measured and estimated values. The ratio of measured Vs30 and estimated Vs30 is compared for both ASTER 1c and SRTM 30c maps in Figure 8a,b.

5. Terrain Classification and Site Characterization

In this study, the terrain classification is referred to as geomorphological classification using elevation data. We followed an automated scale-dependent technique developed by Iwahashi and Pike that considers “geometric signature” as a tool of numerical terrain classification [39]. The three components used in the above-mentioned classification method are as follows: slope gradient, surface texture, and local convexity. The slope gradient described in Section 3 is a comprehensive element that describes the surficial processes. Iwahashi and Pike believed that the presence of slope gradient component is enough to describe the elevation. This meant that in the terrain classification process, mere elevation data were not included [39]. Surface texture as a second element was extracted from the ASTER 1c DEM too. The term “texture” indicates the coarseness of the surface or grain sizes. For example, fine textures in local areas have more “peaks” and “troughs” while coarse textures are flatter. Here, the statistical noises which are inherent in digital elevation models were reduced using a $3\times 3$ median filter. One important contribution of the texture information is that the texture can enrich the slope information extracted from the first element, in which, some lithologies can be differentiated from the others. For instance, local areas located in coarse texture and a planar slope may explain a Quaternary plain. Local convexity as a third element measures convexity of the features in which the value of 1 represents no holes, and the values less than 1 represent concave features. Convexity can be calculated as follow:

$$convextiy = \left(\frac{l}{L}\right)$$

(3)

where $l$ and $L$ are length of convex hull and Euclidian length, respectively. Positive values are upward convexity and negative values are concave areas. The classification algorithm is able to classify the terrain into 8, 12, and 16 classes. We classified the terrain into 16 classes to have more diverse classes. The description of different geomorphological information is represented in a colored matrix (Figure 10a). Geomorphological descriptions of classes in pixel-based and object-based manners are given in Iwahashi and pike [39] and Iwahashi et al. [40]. However, the pixel-based approach, in some cases (e.g., geomorphological map of Thailand), has shown reliable results [41]. Therefore, we also followed the pixel-based method. Each family of classes had rather the same color and characteristics. For instance, the colors near to pink (units 2, 4, 6, and 8 ((Family 2)) describe volcanic landforms in a pixel-based manner (Figure 10b).

6. Terrain Classification Results

To determine where exactly the seismic stations are constructed, we carried out a geomorphological classification for the whole study area with a 10-pixel window size. It must be noted that we only classified the study area using ASTER 1c because terrain classification in global scale had been done using 90 m and 280 m elevation data in the previous studies [39,40]. The aim was to add more attributes for the existing seismic sites with an affordable DEM-based approach. Unlike the Vs30 regression analysis, we did not assign pre-defined values for water bodies (e.g., lakes, wetlands); these areas were removed from ASTER 1c. The majority of the stations in different seismotectonic zones were categorized as follows:

In the Kopet Dagh zone, eight stations out of 34 were categorized in class 12, in the Alborz-Azarbaijan zone, 24 stations out of 103 were categorized in class 9, in the Lut block, 23 stations out of 55 were categorized in class 6, in the Makran zone, 3 stations out of 7 were categorized in class 12, in the Zagros zone, 72 stations out of 190 were categorized in class 12, and in the Central Iran block, 28 stations out of 124 were categorized in class 12. In addition, the results show that 98% of the stations are situated in six classes 30% of which are in class 12, 27% in class 6, 17% in class 9, 16% in class 3, 4% in class 4, and the rest of the stations are located in other classes (Figure 11). This implies that the stations are not only located in volcanic sites, but also in desert plains of weak rocks. We eliminated other classes with insufficient Vs30 observations and tried to plot only the main six units (3, 6, 7, 9, 10, and 12) versus the measured Vs30 values in Figure 12. Although most of the stations in Figure 12 belong to class C in NEHRP classification, the standard deviation bars for class 3 and 12 have a tendency to the upward NEHRP class (i.e., class D) and the downward NEHRP class (class B). The median value for mountain classes like 3 and 7 is relatively high, whilst for smoother classes like 12, it is relatively low (Figure 12). Intermediate classes like basaltic lava plains also showed lower Vs30 values.

7. Discussion

The regression analysis of Vs30 shows that the correlation between slope and Vs30 values in the Iranian territory is low. The lack of a meaningful trend is not unexpected because sedimentary slopes have a better correlation with Vs30 measurements (Wald and Allen [15]) and most of the stations used in this study are probably constructed in rock slopes. Thus, the lack of correlation in rock units may obscure the trend interpretations. The terrain surface classification also showed that most of the stations were constructed in mountainous or intermediate positions (Figure 12). Generally, this trend is not important because steep hill slopes and ridge positions are not suitable lands for building structures, whereas flat and gentle slopes tend to consist of more fertile soils suitable for inhabitation. Thus, the calculated Vs30 values in lower slopes can be used profitably in the development of damage/loss models. The histogram resulting from the ASTER 1c map shows a wider distribution in the comparison of the observed and estimated Vs30 (Figure 9a). However, unlike most other regions studied by Wald and Allen [15] (e.g., California, Italy, and Taiwan), many of the Iranian sites were constructed on rock, so the distribution is biased toward higher-velocity Vs30 data. Accordingly, the correlation of the measured Vs30 and estimated Vs30 were modified for both ASTER 1c and SRTM 30c maps by ignoring Vs30 data larger than 1000 m/s. The logarithmic median values in the new histograms are now narrower and closer to zero, which means that the misfit between the measured and estimated values for both ASTER 1c and SRTM 30c was reduced (Figure 9c and d), but the standard deviation of the residuals is not significantly improved in SRTM 30c. This procedure shows that the ISMN stations are designed on rock sites to avoid noisy records of accelerometers. Thus, the 514 Vs30 measurements are expected to reflect high Vs30 values. Probably, rock slope is not indicative of Vs30, but the correlation between Vs30 and slope within sediments (e.g., young alluviums) reflects stream power.

8. Conclusions

In this study, we focused on the use of two different digital elevation models with two different resolutions (medium and low resolutions). Although the results of the regression analysis cannot explain the Vs30 measurements perfectly, the combination of ASTER 1c and SRTM 30c amplification maps can be useful in the detection of geological and geomorphological units. For the amplification ratio (SRTM 30c/ASTER 1c) of NW Iran (Figure 6), the blue pixels indicate where the SRTM-based amplification is higher and the red pixels indicate where the ASTER-based amplification is higher. The white-to-yellow pixels represent a similar amplification and are thus commonly related with low gradients. For example, the North Tabriz fault (NW-SE trend) with red pixel values and the Tabriz basin with blue pixel values can be explicitly identified using a combined amplification map (Figure 6). Thus, while medium-resolution DEM is able to delineate the position of small slope changes, low-resolution DEM results in more realistic maps in flat areas.

In addition, the geomorphological analysis showed that under a pre-defined scale or radius, the stations are not only located in volcanic sites, but also in desert plains of weak rocks. However, significant changes in the scale could affect the number of stations in each geomorphological unit. Therefore, the scale should be selected carefully according to the purposes.