Next Article in Journal
Assessment and Benchmarking of Spatially Enabled RDF Stores for the Next Generation of Spatial Data Infrastructure
Next Article in Special Issue
Spatiotemporal Change Analysis of Earthquake Emergency Information Based on Microblog Data: A Case Study of the “8.8” Jiuzhaigou Earthquake
Previous Article in Journal
Anomalous Urban Mobility Pattern Detection Based on GPS Trajectories and POI Data
Previous Article in Special Issue
A New Agent-Based Methodology for the Seismic Vulnerability Assessment of Urban Areas

VS30 Seismic Microzoning Based on a Geomorphology Map: Experimental Case Study of Chiang Mai, Chiang Rai, and Lamphun, Thailand

Department of Architecture and Building Engineering, School of Environment and Society, Tokyo Institute of Technology, Yokohama 226-8502, Japan
Survey Division, Department of Survey, Electricity Generating Authority of Thailand, 53 Moo 2 Charansanitwong Road, Bang Kruai, Nonthaburi 11130, Thailand
Department of Civil Engineering, Faculty of Engineering, Thammasat University, Pathumthani 12120, Thailand
Geospatial Information Authority of Japan, Geography and Crustal Dynamics Research Center, Ibaraki 305-0811, Japan
Author to whom correspondence should be addressed.
ISPRS Int. J. Geo-Inf. 2019, 8(7), 309;
Received: 13 May 2019 / Revised: 8 July 2019 / Accepted: 13 July 2019 / Published: 18 July 2019
(This article belongs to the Special Issue Geomatics and Geo-Information in Earthquake Studies)


Thailand is not known to be an earthquake-prone country; however, in 2014, an unexpected moderate earthquake caused severe damage to infrastructure and resulted in public panic. This event caught public attention and raised awareness of national seismic disaster management. However, the expertise and primary data required for implementation of seismic disaster management are insufficient, including data on soil character which are used in amplification analyses for further ground motion prediction evaluations. Therefore, in this study, soil characterization was performed to understand the seismic responses of soil rigidity. The final output is presented in a seismic microzoning map. A geomorphology map was selected as the base map for the analysis. The geomorphology units were assigned with a time-averaged shear wave velocity of 30 m (VS30), which was collected by the spatial autocorrelation (SPAC) method of microtremor array measurements. The VS30 values were obtained from the phase velocity of the Rayleigh wave corresponding to a 40 m wavelength (C(40)). From the point feature, the VS30 values were transformed into polygonal features based on the geomorphological characteristics. Additionally, the automated geomorphology classification was explored in this study. Then, the seismic microzones were compared with the locations of major damage from the 2014 records for validation. The results from this study include geomorphological classification and seismic microzoning. The results suggest that the geomorphology units obtained from a pixel-based classification can be recommended for use in seismic microzoning. For seismic microzoning, the results show mainly stiff soil and soft rocks in the study area, and these geomorphological units have relatively high amplifications. The results of this study provide a valuable base map for further disaster management analyses.
Keywords: microzoning map; geomorphology classification; time-averaged shear wave velocity at 30 m (VS30); soil characteristics microzoning map; geomorphology classification; time-averaged shear wave velocity at 30 m (VS30); soil characteristics

1. Introduction

Thailand is not a highly earthquake-prone country compared to other countries. Relatively higher frequency and magnitude seismic activities are constrained to the northern and western parts of the country [1]. This research covered the study area of Chiang Rai, Chiang Mai, and some parts of Lumphun. There are many active faults located in these provinces (see Figure 1). In 2014, a moderate seismic eruption occurred in the Chiang Rai Province, however, the damage was severe due to inefficient seismic disaster management [2]. Past disaster management was not as effective as it should have been due to fewer experience of seismic activity, which has resulted in insufficient data, understanding, and expertise for seismic hazard planning in terms of mitigation and response. However, Nakasu (2017) reported statistical information about disasters in Thailand, including fatalities and economic losses. The results of this study showed that financial losses due to damage caused by earthquakes are higher than for other hazards (e.g., flooding, storms, and heavy rainfall) that occur with greater frequency [3].
Historically, seismic events did not cause noticeable damage, and thus, seismic disaster management awareness was overlooked. The moderate seismic activity that occurred in 2014 provoked strong local recognition of seismic disaster management. However, seismic analysis data and ground-truth data that provide information on the characteristics of each soil layer are not available in the area. Therefore, in this research, the soil characteristics were explored with field work and using time-averaged shear wave velocity data at 30 m (VS30) from past research for seismic amplification analyses.
VS30 is widely applied to understand soil characteristics and there are several methods for performing VS30 field observations. In this study, microtremor measurements were obtained using the microtremor array measurement with spatial autocorrelation (SPAC) method [4,5,6,7,8]. This technique is popular due to its convenience. The VS30 results in the field survey represent the subsoil under the sensor. The values were then transformed to seismic microzoning with the geomorphological base map. The relationship between VS30 and the amplification factor were referenced from the United States National Earthquake Hazards Reduction Program (NEAR) provisions. The geomorphological units were applied to create the polygon feature from the point feature of VS30, representative of the soil characteristics. A geomorphological map of the study area was not available. Therefore, the geomorphology classification proposed by Iwahashi and Pike (2007) [9] was followed in this study, with different samples in the areas of interest (AOI) and resolutions to observe the influences. The classification was created with the automated decision tree method with a digital elevation model (DEM) as input. Additionally, the classification methodology of pixel-based classification (PBC) and object-based classification (OBC) were explored. Both PBC and OBC are methodologies for putting pixels of images into categories. OBC categorizes pixels based on the spatial relationship with the surrounding pixels, while PBC categorizes pixels based on information from a single pixel. The differences between PBC and OBC outputs are also compared with the geological map in the Discussion section; this geological map is currently the map most commonly used for seismic analysis and other applications in Thailand. Additionally, the sites severely damaged during the 2014 earthquake were overlaid on the seismic hazard map that was developed, for visual evaluation. In conducting this study, no specific grants were received from funding agencies in the public, commercial, or nonprofit sectors.

2. Seismic Microzoning in Thailand

In Thailand, interest in seismic hazard studies has been limited due to the low number of such seismic events. However, recently, the number of such studies has increased due to higher awareness of seismic hazards. In this section, seismic microzoning and other seismic-related studies conducted in Thailand are reviewed.
Ashford et al. (2006) studied the amplification of peak ground acceleration and spectral acceleration in Bangkok [10]. In the same year, Warnitchai et al. (2006) also studied seismic effects in Bangkok due to a distant earthquake [11]. Arai and Yamazaki (2007) performed a study of VS30 based on frequency-wavenumber (f-k) spectral analyses [12]. The results from those studies were then used for seismic microzoning in the Bangkok metropolitan area. Poovarodom and Plalinyot (2011, 2012) conducted a similar study using the SPAC method for the analysis process at the same study site [13,14]. Shibuya et al. (2003) studied VS30 using the seismic cone penetration test [15]. Tuladhar (2003) also determined VS30 values from boreholes from the Bangkok area [16]. Chantamaras (2007) performed a comparison between the multichannel analysis of surface wave (MASW) and borehole methods for obtaining VS30 values in the Bangkok metropolitan area using NEHRP standards [17]. Palsri et al. (2009) studied VS30 with N values obtained from standard penetration test observations in the Chiang Mai, Chiang Rai, and Bangkok metropolitan areas [18]. Pitakwong (2010) compared 2sSPAC and horizontal-to-vertical (H/V) spectral ratio observations in the Chiang Mai, Chiang Rai, and Bangkok metropolitan areas [19,20]. The above studies examining soil characteristics produced different types of hazard maps. However, past studies did not extensively explore certain areas, particularly Chiang Rai. In this study, we conducted field surveys and used VS30 values from previous studies in an area with a high frequency of earthquake incidents [21]. Seismic studies in Thailand began later than in other earthquake-prone countries. As such, certain areas remain to be addressed, especially in the northern part of the country where there are many active faults.
There are many past studies that observed soil characteristics to create seismic hazard maps. Thitimakorn et al. (2010, 2012) applied the MASW method to acquire VS30 values in Bangkok in 2010, and downtown Chiang Rai in 2012. The results were based on the NEHRP classes by spatial interpolation [22,23]. In comparison with this research, the current study covers a larger site. Therefore, data collection during the field survey was challenging. The SPAC method was selected for its convenience in terms of time, labor, and other resource consumption. Another past study that represents soil characteristics as a seismic hazard map was that by Tuladhar et al. (2003) [24]. This study used more than 150 stations to determine the H/V spectral ratio observation in the Bangkok metropolitan area. The results of the natural period ranges of the sites, their effects on tall buildings, and the long natural periods of the site effects on tall infrastructure and long-span bridges were discussed. These past studies represented the soil characteristics in the form of polygon-based features from the point-based feature. The representation of soil characteristics based on some proxy is another main consideration point in the current research.
Matsuoka et al. (2006) suggested the use of geomorphology maps for VS30 mapping [25,26]. The results from this study were used throughout Japan for site classification and applied to the web-based quick estimation system of strong ground motion. Later, Yamazaki et al. (2002) and Matsuoka et al. (2012) also demonstrated the relationship between geomorphologic units and site amplification ratios in Japan [27,28]. Furthermore, geomorphology was also applied in a liquefaction potential analysis and other [29,30,31]. In the current research, geomorphology units were chosen for their accuracy in representing the subsoil characteristics.
In contrast to past seismic research in Thailand, this study covered a wider area, where seismic activity is high and where there is a high level of seismic disaster awareness among the local population. The methodology of seismic microzoning was based on geomorphology units, which have rarely been used in Thailand.

3. Study Area

This study concentrated on Chiang Mai, Chiang Rai, and some parts of the Lamphun province in the north of Thailand (see Figure 1), which is surrounded by active faults [1]. The frequency of earthquake events is relatively higher than in other parts of the country. However, the major focus in terms of disaster management is on Chiang Rai, where a major seismic event occurred in 2014 (Phayao fault). In addition, in the same area and located relatively close to the center of the 2014 event, is the Mea Chan fault. This fault has the highest slip rate of all faults in Thailand. The local population is highly interested in what potential damage may occur in a future earthquake.
From a social perspective, this study area has rapidly increased in terms of its economic status. According to the National Economic and Social Development Board of Thailand (NESDB), the gross provincial products (GPPs) of Chiang Mai and Chiang Rai are $55,218.56 and $37,189.54, respectively (2.58% and 1.74% of the national gross domestic product, respectively) [32]. Chiang Mai’s economic activity and development are equivalent to those of Bangkok (the capital of Thailand). In 2017, the Digital Economy Promotion Agency (DEPA) announced investments worth 36.5 million baht for Chiang Mai’s “smart city” project [33]. The project includes smart transportation, smart agriculture, smart healthcare, and smart tourism. In addition, although the current economic value is not as high as that of Chiang Mai, many development activities are planned for Chiang Rai. In the 2015 Prime Minister’s official announcement, Chiang Rai was selected to be in phase II of the development of special economic zones (SEZs), with Asian Development Bank support [33]. Chiang Rai was selected as the connection point for cooperation between six countries: Cambodia, the People’s Republic of China, Lao People’s Democratic Republic, Myanmar, Thailand, and Vietnam (Greater Mekong subregion). The primary purpose of developing SEZs is to activate local trading competition, reduce the uneven development distribution, and promote regional economic stability to prepare for the Association of Southeast Asian Nations summit.
The potential investments will expand the infrastructure and increase population density, which will result in higher exposure and vulnerability in terms of risk management, making it essential to obtain data for seismic risk analysis. This study aimed to support the sustainable and resilient development of the area with seismic microzoning, which will be a useful tool in urban planning and disaster management, i.e., preparation, mitigation, response, and recovery.

4. Geomorphological Classification

The term geomorphology refers to landform features, which result from physical, biological, and chemical processes. There are several methods used to interpret geomorphological maps, e.g., manual interpretation using orthophoto visualization and field surveys by geomorphologists or other experts, and automated classification without human expertise utilizing remote sensing data. In this research, automatic classification was applied using a DEM based on geo-signatures, i.e., slope, convexity, and texture referenced from the past study of Iwahashi and Pike (2007) [9].
Few geomorphological maps are available in Thailand. The current practice in seismic studies in Thailand bases the soil characteristics on geology. The concept of geomorphology in Thailand was first introduced by Japanese researchers. The purpose of early geomorphology studies in Thailand was for flood disaster planning and mitigation. The map production method involves manual interpretation by specialists based on aerial images, which is combined with available maps or other data [34,35]. In past studies, geomorphological maps of Thailand were available only in Bangkok watershed areas. Since no geomorphological map is available in the study area, the geomorphological classification was also covered in this study.
Later, studies of geomorphological mapping pivoted toward automatic classification based on DEMs, owing to less time being required and without the essential need for expertise in interpretations. The geo-signatures that appear in past studies comprise slopes, aspects, elevations, convexities (curvatures), textures (peaks, pits, ridges, and passes) and the normalized difference vegetation index (NDVI) [36,37,38,39,40,41,42,43,44,45,46,47,48,49]. Different studies proposed the use of the methodology with different techniques, which require different DEM resolutions, filter window sizes, and numbers of classifying units. The literature review in this research scaled down the factors of classification to the classification method, i.e., PBC, OBC, AOI (such as the coverage scene size of the study area), and DEM resolution.
The studies of Iwahashi and Pike (2007), and Iwahashi et al. (2018), were selected in this study [9,50], as these researchers utilized the same geo-signature indices with the PBC and OBC methods, respectively. The indices include slope, convexity, and texture (i.e., the density of pits and peaks). In the PBC presented in the Iwahashi and Pike (2007) study, the digital number values (DN) from the three indices were classified by a decision tree with the mean value threshold of each signature [9]. As the classification is based on the mean value, the number of pixels involved, and the resolution affect the result. The methodologies of two AOIs, i.e., provincial and continental coverage scales, were applied in the current study. AOI1 covers Chiang Mai, Chiang Rai, and part of Lumphun, while AOI2 covers continental Asia-Oceania (see Figure 2). Both scenes were available in 1000 and 250 m resolutions. The classification of four samples resulted in 16 classes. The description of classified categories of geomorphologic units (see Figure 3 and Figure 4) is applicable to the global landscape broad definition, which should be referred to with caution. To obtain 16 geomorphology units, four major groups of geomorphology units were first classified, i.e., mountain, high-land slope, weak rock plain, and low-land plain; changes in AOIs and resolutions affect the pixels in these four major classes. Thermal change detection was applied to detect the movement of the class changing for each pixel between each sample, according to the changes in the AOIs and resolutions.
The results in this study that followed the methodology of automating PBC from the past research of Iwahashi and Pike (2007) using four samples, suggested that the best factors were in PBC (see Section 6.1). The best PBC results were then compared with the OBC results from the study of Iwahashi et al. (2018) to select the geomorphology map for seismic microzoning (see Section 6.2). The classified geomorphology units were compared with the VS30 observations in the districts of Chiang Mai, Chiang Rai, and (some parts of) Lamphun, using statistical significance and correlation tests, to obtain a seismic microzoning map. The geomorphological classification practice used in this paper can possibly be used as a guideline for geomorphological classifications at other sites in Thailand for different purposes.

5. Observation of Time-Averaged Shear Wave Velocity

VS30 indicates the rigidity of the soil. This measure assumes a particular relationship with site amplification, which is included in the ground motion prediction equation [26]. The average spectral amplification factor also corresponds to VS30. There are many data collection methodologies. The method with the highest accuracy for observing VS30 is borehole logging, which can directly examine the soil characteristics in each sediment layer. However, drilling a borehole to a depth of 30 m is costly and resource intensive. To cover the entire study area, with limited time and workforce, the microtremor utilizing technique was applied in this study. The term “microtremor” refers to the small waves that propagate through the ground. These originate from several sources near ground zero such as human activities, wind, and ocean tides.
In total, observations from 73 stations in the districts of Chiang Rai, Chiang Mai, and some parts of Lamphun (see Figure 1) were considered in this study. PAC2016 (red dots in Figure 1) represent data obtained from 23 stations during the current research. The observations from previous studies include 41 stations from NKN2014 and 6 stations from NKN2015 (blue and yellow dots, respectively, in Figure 1) [51]. For the PAC2016 dataset, the SPAC method was applied to analyze microtremors [4,5,6,7,8]. In this method, four sensors were placed in an array to detect the phase velocity from tremors that arrive from every direction. The array of the sensors consisted of one sensor in the middle and three sensors at the 1 m radians (see Figure 5). This research took about one week to conduct 27 stations (PAC2016 dataset). There are some factors that can cause disturbance to the measured signal. However, the site that were located in remote areas had higher systematic noise due to insufficient tremors. Four stations were eliminated. From sensor’s configuration, the duration of the observation period was 20 min. For data treatment, the first 5 min and last 3 min were discarded. The waveform was split for every 20 s. The periods that had high amplitudes due to noise during the observation were removed. Noise refers to abnormal activity that creates higher than average amplitudes, for example, due to heavy vehicles passing. The VS30 values collected in this study (PAC2016) were obtained using the concept of the phase velocity of the Rayleigh wave corresponding to a 40-m wavelength (C(40)). VS30 can be obtained without creating a soil profile from the phase dispersion curve inversion.
Ekrem et al. (2010) also showed a high correlation between the VS30 values and phase velocity C(40) [52]. This approach was proposed by Konno et al. (2000), who examined 85 sites around Tokyo [53]. A method of VS30 estimation was proposed in this study using the phase velocity of the Rayleigh wave fundamental mode obtained from microtremor array observations. The results from that study showed that VS30 was nearly equal to the phase velocity at wavelengths ranging from 35 to 40 m in the small array. Validation in the current study is was not possible at this site owing to the lack of ground-truth data, which was the main obstacle and the reason we used the microtremor observation technique along with the existing methodology. However, the validation of VS30 values from the past study showed that the regression coefficients of VS30-C(40) were 0.9204, 17.86, (±13.83) and the correlation coefficient was 0.982. It was concluded that VS30 ≈ C(35, 40) [53].
The obtained VS30 values classify the sites based on NEHRP provisions [54,55,56]. The VS30 values obtained as the point features from stations were later compared with the classified geomorphological unit (see Section 4) in which they landed: the correlations between VS30 values and geomorphology units are presented in Section 6.2. In other words, this research aimed to identify the VS30 of geomorphological units. However, not every geomorphological unit location could be accessed for VS30 measurements as site accessibility and time constraints restricted access. For the geomorphological units that did not have a VS30 station available, NEHRP classes that have similar geo-signatures were used from the unit classes.

6. Analysis and Results

The approach in this study was broadly divided into two parts: geomorphology classification and VS30 investigation. The geomorphology classification involved studying the classification factors for automating geomorphological classifications from a DEM. The research issues in this part included choosing between the PBC and OBC methodologies, determining the appropriate AOI for the provincial scale of classifications, selection of the data resolution, and the AOI and resolution factor effects on the classification results. Notably, the findings of this study serve only the research in this particular study area (see Section 4). For the VS30 research parts, the findings of this study are applicable to the available methodologies to explore the soil characteristics of the study area. The aims were to cover most of the geomorphology units with fewer resources. The data examined were matched with the geomorphological classifications from the first part of the approach to result in seismic microzoning, which references the NEHRP provision (see Section 5).
The results are presented in this section. The first part of the research results is presented in Section 6.1. The relevant factors for geomorphological classifications in this particular study area are determined. The classified geomorphology map is paired with the VS30 data from the fieldwork and past studies to identify the VS30 values for each geomorphological unit, which is presented in Section 6.2 with the statistical checks. In Section 6.3, the final seismic microzoning based on the geomorphological units of the study area is shown.

6.1. AOIs and Resolution Factors in Geomorphology Classification

Because there is no geomorphological map available for this study area, this research included exploration of geomorphological classification. In this section, the factors for classification were tested with four samples: AOIs 1 and 2, which differed in terms of scene coverage, and the resolutions for each scene, i.e., 250 and 1000 m (see Figure 2). Each DEM sample was classified using the methodology of Iwahashi and Pike (2007); that is, the automated DEM pixel-base geomorphological classification based on slope, curvature, and texture. For the results, the units from each sample covered the difference coverage pixel. Thermal change detection was applied to detect the movement of the class change in each pixel between each sample, according to the changes in the AOIs and resolutions. The number of class changes in every direction case was counted; for example, the class of the pixel changed from 1 to 2 or 2 to 1, with the number of class changes and the change direction for each case.
Changing the resolution mainly affects the changes within a major geomorphological group. For example, for the major geomorphological group of the weak rock plain (see Figure 3), which included unit numbers 11, 15, 9, and 13 (see Figure 4), changing the resolution affects the changes between those unit numbers. From VS30 in Section 6.2, classes 11 and 15 have significantly different VS30 values, reflecting a substantial difference in soil properties. Changes in the AOIs mainly affect the changes between the major geomorphology groups of the geomorphological units (see Figure 3), e.g., it shifts from mountain to weak rock plain, mountain to the high-land slope, or from weak rock plain to low-land plain. The AOI suitability is understood by comparing the results with an available geological map or other maps to confirm its boundaries. Therefore, a 250 m resolution AOI2 was selected for PBC.

6.2. Geomorphology Maps for VS30 Representation

In this section, the suitability of the PBC and OBC base maps for representing VS30 is discussed. The geomorphology map from PBC in this study was generated based on the methodology from the study of Iwahashi and Pike (2007) with different classification factors (see Section 6.1), while the geomorphology map from OBC was from Iwahashi et al.’s (2018) results [9,50]. Both the PBC and OBC geomorphological maps utilized the same geo-signature indices, i.e., slope, convexity, and texture. However, the procedural differences resulted in a different number of geomorphological units in the output. The comparisons and descriptions of the PBC and OBC units are shown in Figure 4. The number of classes in high sedimentation areas is higher in the PBC map.
Each unit from PBC and OBC was plotted against the VS30 values (see Figure 6 and Figure 7). Notably, not all of the geomorphological units can be determined because the observations were insufficient due to site accessibility. Some of the geomorphological units that were statistically unusable, which were mostly located on steep slopes or in forested areas, were excluded. The VS30 plot focused on the NEHRP provision class due to the same amplification behavior property. The standard deviation of the unit classes is considered to be a better proxy for identifying the property based on VS30 values. Hence, the statistical analysis was also applied.
From the observations, the PBC geomorphological units (see Figure 7) have lower standard deviations than the OBC (see Figure 8). Therefore, the PBC geomorphological map was selected for further analysis. A significance test was performed to test whether the selected representative unit was meaningful. Based on the data characteristics, a nonparametric significance test and the Kruskal-Wallis test were applied for the geomorphology groups representing abilities [57]. At a confidence limit of 95% (p-value < 0.005), the p-values of the significance tests for the units from PBC and OBC were 0.004319 and 0.03096, respectively. However, some units were eliminated due to an insufficient number of observations, and some units were also eliminated if their values were higher than the absolute standard deviation. The significance test results showed that the PBC units represented the VS30 value at a significant level. Hence, the correlation between the PBC and OBC geomorphological units and the VS30 values was obtained using the statistical computing package Cocor in R programming language [58]. The correlation between the VS30 and PBC geomorphological units was 0.01539, while that of the VS30 and OBC geomorphological units was −0.354. A correlation test was also conducted at a confidence level of 95%. The null hypothesis was “if the correlation numbers between the VS30 and PBC geomorphological unit and the VS30 and OBC geomorphological unit are equal”. From Cocor’s set of hypothesis testing based on the results of the Cocor method, the null hypothesis was rejected. Therefore, the correlation between the VS30 and the PBC geomorphological unit was higher than that of the VS30 and the OBC geomorphological unit.

6.3. Seismic Microzoning Based on Geomorphology with NEHRP Standards

In Section 6.1 and Section 6.2, the study of geomorphological classifications and matching of the geomorphological units to the VS30 value was reported. The point feature data were transformed into polygonal features to create a seismic microzoning map according to the NEHRP provisions. The NEHRP defines five classes of soil, i.e., A, B, C, D, and E, based on the soil rigidity, which represents the VS30 value.
The VS30 values gathered are represented with the classified PBC geomorphological units. The available NEHRP classes in the study area were C and D. The C class has a VS30 value of between 360 and 760 m/s (referring to dense soil or soft rock), while the D class is between 180 and 360 m/s (referring to stiff soil). When referencing NEHRP class D geomorphological units, weak rock, incised terraces next to eroded rock, delta plains, and low-lying fluvial plain areas (unit numbers 11, 12, 13, 14, 15, and 16, respectively; see Figure 7) are included, which are mostly areas along the river due to high sedimentation processes. Class C covers the incised terrace next to the foot of the mountain, moderately eroded mountains, and volcanic regions (units 7, 8, and 15, respectively; see Figure 7).
The final output, Figure 9, shows NEHRP class D in red, and the orange area is class C. The gray regions present mountain areas which could not be covered by the microtremor stations. However, when comparing the geo-signature with the geomorphological unit that included microtremor data, it was similar to the units in class C.

7. Discussion

This section includes discussion of the geomorphological map and geology map comparisons for suitability in representing the VS30 value (see Section 7.1). The geological map was used for the primary manual visualization check with the classified geomorphological map. This check was based on the simple hypothesis that their boundaries must not be different. Additionally, the majority of sites damaged in the 2014 earthquake were overlaid on the resulting seismic hazard map for visual evaluation (see Section 7.2).

7.1. Geomorphological Map from PBC, OBC, Geology Map Comparison

For seismic microzoning, this study aimed to find the best representative for VS30. VS30 varies based on the soil’s properties. According to current practice, geological units are the base map, which is a limitation for the PBC and OBC geomorphological maps. There are two important points of discussion in this subsection: consideration of the PBC and OBC geomorphological map suitability and the ability of this map to become the proxy for VS30 between geological and geomorphological units.
In general, the PBC and OBC methods of classification in photogrammetry have been debated. Past studies have made many comparisons. Adam et al. (2016) and Weih et al. (2010) reported that OBC is suitable for high-resolution image classifications, while PBC is better for landcover classifications [59,60]. As stated in previous sections, this study followed the geomorphological classification method proposed by Iwahashi and Pike (2007) with an adjusted AOI and input resolution. The PBC geomorphological map was selected to obtain the VS30 representing the seismic microzoning because it can define the VS30 values more precisely (see Figure 6 and Figure 7). Figure 10 presents the zoomed-in geological map, the OBC geomorphological map, and the PBC geomorphological map. This figure shows that the boundaries are somewhat agreed (Figure 10) based on visual verification, and that the representational ability results are convincing (Section 6.2). The unit boundary from the OBC geomorphological map agrees better with the geological unit. However, in lowlands, the PBC geomorphological map class units of 11 and 15 (see Figure 10(iii)), which are shown in navy and dark-navy colors, could not be separated from the geological map and OBC geomorphological map (see Figure 10(i,ii)). From VS30, units 11 and 15 are highly different (see Figure 7). The output suggested that the PBC map results in higher details than the OBC map, especially in lowland areas, where seismic microzoning has a greater focus.
Hence, the locations of the VS30 stations are available in six geological units, i.e., fluvial deposit (flood plain [Qff]), terrace and colluvium [Qt], alluvium (Qa), rhyolite [PTrv], semiconsolidate [Tmm], and granite [Trgr] [61]. From Figure 10, the eroded plain of weak rock (major class eroded rock plain) and fluvial plain (in yellow tone) of the geomorphological map is under the fluvial deposit (i.e., Qff, Qt, and Qa) unit of the geology map, which covers NEHRP classes C and D. Therefore, the geology units have limitations in distinguishing the subsoil characteristics based on NEHRP standards.

7.2. Seismic Microzoning and Major Site Damage Comparisons According to the 2014 Incident

In this subsection, the output of the seismic microzoning map is overlaid with the damaged site locations. Figure 11 shows the location of the majorly damaged sites in the 2014 earthquake incident. “Major damage” in this study refers to structures that collapsed. However, the infrastructure in the study area was locally constructed without consideration of standards. There are many active faults in the area, but awareness was lower in the past because no significant strong ground motions or intensive seismic damage had previously occurred. The local houses are mainly two stories with a combination of a concrete frame with a fill-in brick wall for the first floor and a wooden structure on the second floor. Typically, a Thai house is raised off the ground (in preparation for flooding or intensive rainfall), and the first floor is usually built later. Therefore, the damage level may be influenced by the construction standard and the vulnerability of the location. From Figure 11, the damage location mostly corresponds to NEHRP class D. When considering the geomorphological units, all damage sites were found to be located in incised terraces and eroded plains of weak rock (classes 11 and 13). However, the damage that was located in NEHRP class C (one location from Figure 11) may have been influenced more by the construction quality. From an elementary visual evaluation, seismic microzoning based on geomorphological units correlates well with the construction damage.

8. Conclusions

To support the continuing economic growth in the study area, considerations of disaster management are necessary. In this research, the aim was to create a seismic microzoning map in the study area by using the existing methodologies in base map classifications and VS30 observations. The study area comprised Chiang Mai, Chiang Rai, and Lumphun (some parts), which are provinces in northern Thailand where seismic activity is relatively higher than in other areas of the country.
In this research, the methodology from a past study was applied to create a geomorphological map with different factors. This map was later assigned VS30 values from fieldwork observations and data gathered from previous studies. Since there are no ground-truth data available in the study area, this methodology can be used to obtain data within both time and resource constraints. The approach in this study can be divided into two main parts, i.e., microtremor measurements and geomorphological classification. The seismic microzoning in this research refers to subsoil characterization according to NEHRP standards. VS30 values were obtained from 73 stations using microtremor array measurements and the SPAC method. For post-processing, the phase velocity methodology of the Rayleigh wave corresponding to a 40-m wavelength by Konno et al. (2000) [57] was applied without PS logging, which was not available in the study area. Then, the VS30 values were transformed from point features into area maps using geomorphological units. The geomorphological map used in this study was interpreted using a DEM with the methodology from Iwahashi and Pike (2007) [5], which is the automated PBC from geo-signatures, including slope, convexity, and texture. The factors that affect the classification results in PBC were explored. Hence, the representative abilities of the PBC and OBC geomorphological maps and geology maps were also compared.
The results of this study suggest that the PBC method is a suitable classification method for landcover classifications on a provincial scale. The range of VS30 values compared to the PBC geomorphology unit can represent the unit with better precision than OBC. The VS30 values were paired with the geomorphological units for the statistical correlation check. From 16 classes of PBC geomorphological units, the results were re-classified into two NEHRP classes. Only NEHRP class C and D existed in the study area. The output of this study can be used as a reference for future urban development in the study area. For example, it provides information to understand soil conditions for the appropriate construction codes, and to prepare for seismic disaster management. The geomorphological classifications can be used as a guideline for classification in nearby areas to serve other uses.

Author Contributions

M.M. advised the project, the main conceptual ideas, and proof outline. P.T. worked out the technical details and wrote the manuscript with support from M.M. and J.I. P.T. and M.M. performed the microtremor measurement of the 2016 dataset. N.P. provided the 2014 and 2015 datasets for the analysis. P.T. analyzed all the datasets under the suggestions from J.I. and M.M. J.I. also provided the object-based geomorphological map for comparison.


This research was supported in a part by the grant-in-aid for scientific research (KAKENHI No.15H05134).

Conflicts of Interest

This research has no conflict of interest.


  1. Department of Mineral Resources. Active Fault Map in Thailand and Report; Department of Mineral Resources, 2006. Available online: (accessed on 15 June 2018). (In Thai)
  2. Vervaeck, A. Very Strong Deadly Earthquake Close to Chiang Rai, Thailand—At least 1 Dead and 32 Injuries + Huge Number of Aftershocks. Available online: (accessed on 15 June 2018).
  3. Nakasu, T. Natural disasters and disaster management in Thailand: Status, risks, and trends. In Proceedings of the 13th International Conference on Thai Studies: Globalized Thailand? Connectivity, Conflict and Conundrums of Thai Studies, Chiang Mai, Thailand, 15–18 July 2017. [Google Scholar]
  4. Department of Public Works and Town & Country Planning, Ministry of Interior. Earthquake Resistance Building Control Standards; Department of Public Works and Town & Country Planning, Ministry of Interior: Bangkok, Thailand, 2009. (In Thai)
  5. Kanai, K.; Tanaka, T.; Osada, K. Measurement of the microtremor. Bull. Earthq. Res. Inst. Tokyo 1954, 32, 199–209. (In Japanese) [Google Scholar]
  6. Aki, K. Space and time spectra of stationary stochastic waves, with special reference to microtremors. Bull. Earthq. Res. Inst. 1957, 35, 415–457. [Google Scholar]
  7. Horike, M. Inversion of the phase velocity of long-period microtremors to the S-wave-velocity structure down to the basement in urbanized areas. J. Phys. Earth 1985, 33, 59–96. [Google Scholar] [CrossRef]
  8. Okada, H.; Matsushima, T.; Moriya, T.; Sasatani, T. An exploration technique using long-period microtremors for determination of deep geological structures under urbanized areas. Butsuri-Tansa (Geophys. Explor.) 1990, 43, 402–417. (In Japanese) [Google Scholar]
  9. Iwahashi, J.; Pike, R.J. Automated classification of topography from DEMs by an unsupervised nested-means algorithm and a three-part geometric signature. Geomorphology 2007, 86, 409–440. [Google Scholar] [CrossRef]
  10. Ashford, S.A.; Jakrapiyanun, W.; Lukkunaprasit, P. Amplification of earthquake ground motions in Bangkok. In Proceedings of the 12th World Conference on Earthquake Engineering, Auckland, New Zeland, 30 January–4 February 2000. [Google Scholar]
  11. Warnitchai, P.; Sangarayakul, C.; Ashford, S.A. Seismic hazard in Bangkok due to long-distance earthquakes. In Proceedings of the 12th World Conference on Earthquake Engineering, Auckland, New Zeland, 30 January–4 February 2000. [Google Scholar]
  12. Arai, H.; Yamazaki, F.; (Earthquake Disaster Mitigation Research Center). Exploration of S-Wave Velocity Structure Using Microtremor Arrays in the Greater Bangkok, Thailand; EDM Technical Report No.15; Earthquake Disaster Mitigation Research Center: Bangkok, Thailand, 2002. [Google Scholar]
  13. Poovarodom, N.; Plalinyot, N. Site Characterization in the Greater Bangkok Area by Microtremor Observations. J. Earthq. Eng. 2013, 17, 209–226. [Google Scholar] [CrossRef]
  14. Jirasakjamroonsri, A.; Poovarodom, N.; Warnitchai, P. Seismic site characteristics of shallow sediments in Bangkok metropolitan region and their inherent relations. Bull. Eng. Geol. Environ. 2019, 78, 1327–1343. [Google Scholar] [CrossRef]
  15. Shibuya, S.; Tamrakar, S.B.; Manakul, W. Geotechnical Hazard in Bangkok—Present and Future. Lowl. Technol. Int. 2003, 5, 1–13. [Google Scholar]
  16. Tuladhar, R. Classification of Soil Profile and Seismic Response Analysis (Elastic) in the Greater Bangkok Area; Internal Report, Structural Engineering Field of Study; Asian Institute of Technology: Khlong Luang District, Thailand, 2003. [Google Scholar]
  17. Chantamas, P. The Comparison the Shear Wave Velocity of Bangkok Soil Deposit from Multi Channel Signal Receiver and Borehole. Doctoral Thesis, Department of Geology, Faculty of Science, Chulalongkorn Unitversity, Bangkok, Thailand, 2007. (In Thai). [Google Scholar]
  18. Palsri, J.; Reungrasimi, A. The relationship between the shear wave velocity and the N value with the unsaturated soil for shear force analysis in the north of Thailand, GTE50407. In Proceedings of the 14th National Civil Engineer Conference, Suranaree University of Technology, Muang District, Thailand, 2009. (In Thai). [Google Scholar]
  19. Pitikwong, K. Estimation of Shear Wave Velocity Using Microtremor Simultaneously. Master’s Thesis, Thammasat University, Phra Nakhon, Thailand, 2010. [Google Scholar]
  20. Poovarodom, N.; Pitakwong, K. Exploration of site characteristics of subsoils by microtremor observations. Res. Dev. J. 2010, 21, 3. [Google Scholar]
  21. The Thailand Research Fund. Full Report of Reduction Hazard from Earthquake in Thailand; Phase 2; The Thailand Research Fund: Bangkok, Thailand, 2010. [Google Scholar]
  22. Thitimakorn, T. Comparison of shear-wave velocity profiles of Bangkok subsoils from multi-channel analysis of surface wave and downhole seismic methods. J. Appl. Sci. Res. 2010, 6, 1953–1959. [Google Scholar]
  23. Thitimakorn, T.; Channoo, S. Shear wave velocity of Soils and NEHRP site classification map of Chiangrai City, northern Thailand. EJGE 2012, 17, 2891–2904. [Google Scholar]
  24. Tuladhar, R.; Yamasaki, F.; Warntichai, P.; Saita, J. Seismic microzonation of the greater Bangkok area using microtremor observation. Earthq. Eng. Struct. Dyn. 2003, 33, 211–225. [Google Scholar] [CrossRef]
  25. Matsuoka, M.; Wakamatsu, K.; Fujimoto, A.; Midorikawa, S. Average Shear-wave velocity mapping using Japan Engineering Geomorphologic classification map. Struct. Eng./Earthq. Eng. JSCE 2006, 23, 57s–68s. [Google Scholar] [CrossRef]
  26. Matsuoka, M.; Yamamoto, N. Web-based quick estimation system of strong ground motion maps using engineering geomorphologic classification map and observed seismic records. In Proceedings of the 15th World Conference on Earthquake Engineering, Lisbon, Portugal, 24–28 September 2012. [Google Scholar]
  27. Yamazaki, F.; Wakamatsu, K.; Onishi, J.; Shabestari, K.T. Relationship between geomorphological land classification and site amplification ratio based on JMA strong motion records. Soil Dyn. Earthq. Eng. 2000, 19, 41–53. [Google Scholar] [CrossRef]
  28. Matsuoka, M.; Wakamatsu, K. Nationwide site amplification zonation study using Japan engineering geomorphologic classification map. In Proceedings of the 13th World Conference on Earthquake Engineering, Vancouver, BC, Canada, 1–6 August 2004. [Google Scholar]
  29. Kotoda, K.; Wakamatsu, K.; Oya, M. Mapping liquefaction potential based on geomorphological land classification. In Proceedings of the 9th World Conference on Earthquake Engineering, Kyoto, Japan, 8–9 August 1988. [Google Scholar]
  30. Wakamatsu, K.; Matsuoka, M.; Hasegawa, K. GIS-based nationwide hazard zoning using the Japan engineering geomorphologic classification map. In Proceedings of the 8th U.S. National Conference on Earthquake Engineering, San Francisco, CA, USA, 18–22 April 2006. [Google Scholar]
  31. Lyell, C. Principles of Geology; Penguin Classics: London, UK, 1854. [Google Scholar]
  32. Boonnoon, J. (2017-01-27); DE Ministry Pushing for Nationwide; The Nation: Bangkok, Thailand, (Retrieved March 2018).
  33. Prime Minister’s Office of Thailand. Special Economic Zone Board Policy; Designated Special Economic development zone Phase II; Prime Minister’s Office of Thailand: Bangkok, Thailand, 2014. (In Thai)
  34. Oya, M. Applied Geomorphology for Mitigation of Natural Hazards; Advances in Natural Technological Hazards Research; Springer: Berlin/Heidelberg, Germany, 2001. [Google Scholar]
  35. Yamada, Y. Flood analysis using satellite data and geomorphological survey map showing classification of flood-inundated areas. Int. Arch. Photogramm. Remote Sens. 2000, 33 (Suppl. B4), 96–100. [Google Scholar]
  36. Blaszczynski, J.S. Landform characterization with geographic information systems. Photogramm. Eng. Remote Sens. 1997, 63, 183–191. [Google Scholar]
  37. Crevenna, A.C.; Rodriguez, V.T.; Sorani, V.; Frame, D.; Ortiz, M.A. Geomorphometric analysis for characterizing landforms in Morelos State, Mexico. Cienc. UNAM 2005, 67, 407–422. [Google Scholar]
  38. Dikau, R.; Brabb, E.E.; Mark, R.K. Landform classification of New Mexico by Computer. Open-File Rep. (U. S. Geol. Surv.) 1991, 15, 91–634. [Google Scholar]
  39. Dymond, J.R.; Derose, R.C.; Harmsworth, G.R. Automated mapping of land components from digital elevation data. Earth Surf. Process Landf. 1995, 20, 131–137. [Google Scholar] [CrossRef]
  40. Giles, P.T. Geomorphological signatures: Classification of aggregated slope unit objects from digital elevation and remote sensing data. Earth Surf. Process Landf. 1998, 23, 581–594. [Google Scholar] [CrossRef]
  41. Graff, L.H.; Usery, E.L. Automated classification of generic terrain features in digital elevation models. Photogramm. Eng. Remote Sens. 1993, 59, 1409–1417. [Google Scholar]
  42. Hengl, T.; Rossiter, D.G. Supervised landform classification to enhance and replace photo-interpretation in semi-detailed soil survey. Soil Sci. Soc. Am. J. 2003, 67, 1810–1822. [Google Scholar] [CrossRef]
  43. Irvin, B.J.; Ventura, S.J.; Slater, B.K. Fuzzy and isodata classification of landform elements from digital terrain data in Pleasant Valley, Wisconsin. Geoderma 1997, 77, 137–154. [Google Scholar] [CrossRef]
  44. Miliaresis, G.C.; Argialas, D.P. Segmentation of physiographic features from the global digital elevation model/GTOPO30. Comput. Geotech. 1999, 25, 715–728. [Google Scholar] [CrossRef]
  45. Prima, O.D.A.; Echigo, A.; Yokoyama, R.; Yoshida, T. Supervised landform classification of northeast Honshu from DEM-derived thematic maps. Geomorphology 2006, 78, 373–386. [Google Scholar] [CrossRef]
  46. MacMillan, R.A.; Martin, T.C.; Earle, T.J.; McNabb, D.H. Automated analysis and classification of landforms using high-resolution digital elevation data: Applications and issues. Can. J. Remote Sens. 2003, 29, 592–606. [Google Scholar] [CrossRef]
  47. Drăguţ, L.; Eisanka, C. Automated object-based classification of topography from SRTM data. Geomorphology 2012, 141–142, 21–33. [Google Scholar] [CrossRef] [PubMed]
  48. Piloyan, A.; Konečný, M. Semi-automated classification of landform elements in Armenia based on SRTM DEM using K-means unsupervised classification. Quaest. Geogr. 2017, 36, 93–103. [Google Scholar] [CrossRef]
  49. Wayne, W.D. Applied Nonparametric Statistics; Cengage Learning: Boston, MA, USA, 1990. [Google Scholar]
  50. Iwahashi, J.; Kamiya, I.; Matsuoka, M.; Yamazaki, D. Global terrain classification using 280 m DEMs: Segmentation, clustering, and reclassification. Prog. Earth Planet. Sci. 2018, 5, 1. [Google Scholar] [CrossRef]
  51. Poovarodom, N. Full Report of Reduction Hazard from Earthquake in Thailand; Phase I and III, internal report. (In Thai)
  52. Ekrem, Z.; Serdar, O.Z.; Aylin, K.; Tapırdamaz, M.C.; Suna, C.O.Z.; Adil, T.; Bora, E. Shear wave velocity structure of the İzmit bay area (Turkey) estimated from active–passive array surface wave and single-station microtremor methods. Geophys. J. Int. 2010, 182, 1603–1618. [Google Scholar]
  53. Konno, K.; Kataoka, S. An estimating method for the average S-wave velocity of ground from the phase velocity of Rayleigh wave. J. Jpn. Soc. Civ. Eng. 2000, 647/I–51, 415–423. (In Japanese) [Google Scholar]
  54. FEMA. NEHRP Recommended Provisions for Seismic Regulations for New Buildings and Other Structures; FEMA: Washington, DC, USA, 2003; pp. 450–451.
  55. FEMA. NEHRP Recommended Seismic Provisions for New Buildings and other Structures; FEMA: Washington, DC, USA, 2009; p. 750.
  56. Dobry, R.; Borcherdt, R.D.; Crouse, C.B.; Driss, I.M.; Joyner, W.B.; Martin, G.R. New site coefficients and site classification system used in recent building seismic code provisions. Earthq. Spectra 2000, 16, 41–67. [Google Scholar] [CrossRef]
  57. Daniel, W. Kruskal–Wallis one-way analysis of variance by ranks. In Applied Nonparametric Statistics, 2nd ed.; PWS-Kent: Boston, MA, USA, 1990; pp. 226–234. [Google Scholar]
  58. Diedenhofen, B.; Musch, J. Cocor: A Comprehensive Solution for the Statistical Comparison of Correlations. PLoS ONE 2015, 10, e0121945. [Google Scholar] [CrossRef] [PubMed]
  59. Adam, H.E.; Csaplovics, E.; Elhaja, M.E. A comparison of pixel-based and object-based approaches for land use land cover classification in semi-arid areas, Sudan. IOP Conf. Ser. Earth Environ. Sci. 2016, 37, 12061. [Google Scholar] [CrossRef]
  60. Weih, R.C., Jr.; Riggan, N.D., Jr. Object-based classification vs. pixel-based classification: Comparative importance of multi-resolution imagery. Int. Arch. Photogramm. Remote Sens. Spat. Inf. Sci. 2010, 38, C7. [Google Scholar]
  61. Department of Mineral Resources. Thailand Geology Map Scale of 1:250,000; Department of Mineral Resources: Bangkok, Thailand, 1999.
Figure 1. Study area: the districts of Chiang Mai, Chiang Rai, and Lamphun with the VS30 observation stations.
Figure 1. Study area: the districts of Chiang Mai, Chiang Rai, and Lamphun with the VS30 observation stations.
Ijgi 08 00309 g001
Figure 2. Digital elevation model (DEM) scenes for geomorphology classification: (i) AOI1 covers Asia–Oceania (resolution: 1000 m) and (ii) AOI2 covers Chiang Mai, Chiang Rai, and Lamphun. (resolution: 250 m).
Figure 2. Digital elevation model (DEM) scenes for geomorphology classification: (i) AOI1 covers Asia–Oceania (resolution: 1000 m) and (ii) AOI2 covers Chiang Mai, Chiang Rai, and Lamphun. (resolution: 250 m).
Ijgi 08 00309 g002
Figure 3. Major geomorphology groups.
Figure 3. Major geomorphology groups.
Ijgi 08 00309 g003
Figure 4. Descriptions of classified geomorphology units from pixel-based classification (PBC) and object-based classification (OBC) [9,50].
Figure 4. Descriptions of classified geomorphology units from pixel-based classification (PBC) and object-based classification (OBC) [9,50].
Ijgi 08 00309 g004
Figure 5. Sample of array observation station.
Figure 5. Sample of array observation station.
Ijgi 08 00309 g005
Figure 6. Geomorphology map for AOI2 obtained using a 250-m resolution DEM (Figure 2(ii)).
Figure 6. Geomorphology map for AOI2 obtained using a 250-m resolution DEM (Figure 2(ii)).
Ijgi 08 00309 g006
Figure 7. VS30 based on PBC geomorphological units; the red lines indicate the National Earthquake Hazards Reduction Program (NEHRP) classes.
Figure 7. VS30 based on PBC geomorphological units; the red lines indicate the National Earthquake Hazards Reduction Program (NEHRP) classes.
Ijgi 08 00309 g007
Figure 8. VS30 based on OBC geomorphological units; the red lines indicate the NEHRP classes.
Figure 8. VS30 based on OBC geomorphological units; the red lines indicate the NEHRP classes.
Ijgi 08 00309 g008
Figure 9. VS30 map based on the reclassified geomorphological map with NEHRP soil characteristic standards.
Figure 9. VS30 map based on the reclassified geomorphological map with NEHRP soil characteristic standards.
Ijgi 08 00309 g009
Figure 10. (i) Geology map of Thailand [61]: the geological units are fluvial deposit (flood plain [Qff]), terrace and colluvium [Qt], alluvium (Qa), rhyolite [PTrv], semiconsolidate [Tmm], and granite [Trgr]; (ii) OBC geomorphology map from Iwahashi et al. [50]; and (iii) classified PBC geomorphology map. The geomorphological units for(ii,iii) are based on the descriptions in Figure 4.
Figure 10. (i) Geology map of Thailand [61]: the geological units are fluvial deposit (flood plain [Qff]), terrace and colluvium [Qt], alluvium (Qa), rhyolite [PTrv], semiconsolidate [Tmm], and granite [Trgr]; (ii) OBC geomorphology map from Iwahashi et al. [50]; and (iii) classified PBC geomorphology map. The geomorphological units for(ii,iii) are based on the descriptions in Figure 4.
Ijgi 08 00309 g010
Figure 11. Zoomed-in hazard microzoning map with damage locations in 2014.
Figure 11. Zoomed-in hazard microzoning map with damage locations in 2014.
Ijgi 08 00309 g011
Back to TopTop