Effective Identification of Terrain Positions from Gridded DEM Data Using Multimodal Classification Integration

Abstract

Terrain positions are widely used to describe the Earth’s topographic features and play an important role in the studies of landform evolution, soil erosion and hydrological modeling. This work develops a new multimodal classification system with enhanced classification performance by integrating different approaches for terrain position identification. The adopted classification approaches include local terrain attribute (LA)-based and regional terrain attribute (RA)-based, rule-based and supervised, and pixel-based and object-oriented methods. Firstly, a double-level definition scheme is presented for terrain positions. Then, utilizing a hierarchical framework, a multimodal approach is developed by integrating different classification techniques. Finally, an assessment method is established to evaluate the new classification system from different aspects. The experimental results, obtained at a Loess Plateau region in northern China on a 5 m digital elevation model (DEM), show reasonably positional relationship, and larger inter-class and smaller intra-class variances. This indicates that identified terrain positions are consistent with the actual topography from both overall and local perspectives, and have relatively good integrity and rationality. This study demonstrates that the current multimodal classification system, developed by taking advantage of various classification methods, can reflect the geographic meanings and topographic features of terrain positions from different levels.

1. Introduction

As an important factor describing the Earth’s topographic features, terrain position composes of spatially adjacent units with similar topographic features under a certain scale [1,2,3,4]. On a natural slope, different terrain positions have specific topographic and geomorphic features. Different geographic processes (e.g., soil and hydrology, etc.) on various terrain positions will exhibit characteristics of diversification [4,5]. Therefore, terrain position classification plays an important role in many fields, such as landform evolution [6,7], digital soil mapping [8,9,10], hydrological modeling [11,12] and risk analysis [13].

To identify terrain positions, local terrain attributes (e.g., elevation, slope and curvature) are often utilized individually or compositely in different ways according to definite rules [14,15,16,17]. In order to overcome the defect of applying local terrain attributes with a fixed-size window, Weiss [18] presented the topographic position index (TPI) to classify terrain positions based on the fact that landforms are mostly a scale-dependent phenomenon. Although TPI compares the elevation of each cell in its surrounding area with two different window sizes, it is essentially a local terrain attribute. In addition, Jasiewicz et al. [19] proposed a geomorphons approach that serves many ternary patterns as the archetypes of landform elements via counting the elevation differences in a certain neighborhood. To reflect the terrain information of a location in a region, Skidmore [1] developed a regional terrain attribute, namely the relative position index (RPI), to provide an approximate estimate of how far a location is from a ridge or a valley. This approach is more suitable for reflecting the geographical meanings of terrain positions.

Development of reliable rules is the key to terrain position classification. We can build exact classification rules from expert knowledge or a definition of terrain positions, but sometimes these rules are vague and unavailable. Thus, referencing to the image classification technology, supervised [20,21,22] and unsupervised [23,24] classification methods are introduced to implicitly obtain classification rules. Supervised classification is achieved by establishing a discriminant function using selected samples based on statistical identification. However, unsupervised classification directly groups cells into “clusters” based on their properties and then assigns classes without involving any samples.

To take advantage of the segmentation technique over pixel by pixel classification, Drăguţ et al. [25] demonstrated a classification approach based on object-oriented image analysis. This method delineates relatively homogenous objects through image segmentation, and then classifies these objects into landform elements using pre-defined rules. The object-oriented classification considers both the variability of object values and the spatial structural relationships between geographic objects, which can result in better integrity of the classes. As such, object-oriented approaches have become very popular in geomorphometry and widely applied in terrain position classification [26,27,28].

According to different classification criteria, the existing approaches may be roughly grouped into local terrain attribute (LA)-based and regional terrain attribute (RA)-based, or rule-based and supervised and unsupervised, or pixel-based and object-oriented. The characteristics of each of these classification methods vary significantly [29,30]. The LA-based classification usually selects multiple local terrain attributes to compose an attribute set for classification, but the uncertainty associated with the selection of local terrain attributes may result in complex classification processes. Other the other hand, the RA-based classification is simple and easy to implement, and provides better integrity of results, but it usually adopts only one type of terrain attributes for classification, leading to relatively simple classified terrain positions. Generally, in practical applications, the RA-based method is used to simply classify an overall slope from the top to the bottom, whereas the LA-based method is mostly adopted to express the complicated conformations of local slopes effectively.

The rule-based classification can provide exact and unique classes, but the definite rules are sometimes difficult to define. With artificial participation, the supervised classification may increase classification accuracy by controlling sample selection and training samples repetitively. But this may also cause certain issues in that the classification results may not be unique and may be susceptible to subjective factors. Because of the unresolved uncertainty associated with the unsupervised classification results, the application of the unsupervised approach has gradually decreased.

The pixel-based classification is easy to implement. However, when the pixel value in a layer is not sufficiently continuous, the pixel-based method may give fragmented classification results. Object-oriented classification resolves the problem of fragmented classification results by means of conducting firstly object segmentation and then classification. However, due to the uncertainty in setting segmentation parameters, the object-oriented method may be subject to the problems of insufficient segmentation and missing relatively small segmented objects.

In summary, for the classification of terrain positions, the prevailing classification approaches all have their own strengths and weaknesses. Few attempts have been made to take advantage of different methods to develop a multimodal classification system. In addition, the hierarchical conceptualization is commonly considered in landform classification, which provides the possibility of applying different classification methods adaptive to terrain position features at different levels. The aim of this work is to present a multimodal classification system to integrate different approaches by considering different levels of terrain features/components for more effective terrain position classification. The remainder of the work is organized as follows: Section 2 describes the methodology of terrain positon classification as well as the study area and data; the results and discussion are presented in Section 3; and finally several remarks are summarized in Section 4.

2. Materials and Methods

2.1. Study Area and Data

The Jiuyuangou region of the Loess Plateau in the northern Shaanxi Province, China, is selected as the study area for the current study, as shown in Figure 1. The case study site is located in the middle reach of the Wuding River and its topography consists of loess hilly-gully landforms featured with severe erosion, fully developed gullies and broken terrain. The site covers about 100 km² of area and the relative altitude difference from hill top to gully bottom is 374 m, ranging from 814 to 1188 m. The average ground surface slope is 29°, and the density of gullies is 6.52 km/km². A digital elevation model (DEM) of 5 m resolution is available from the National Fundamental Geographic Information Database to support the current study.

2.2. Definition Scheme of Terrain Positions

Herein, a two-level definition scheme is first designed for clarifying terrain positions (

Table 1. Definition of terrain positions (ND—not defined).

Terrain Position		RPI	Curvature (1/m)		Slope (°)	Description
First Level	Second Level		Plan	Profile
Upper slope	Ridge	[β, 1]	ND	ND	ND	Mostly with large covering area and strip-shaped cross-sections
	Hillock		ND	ND	ND	Mainly with small covering area and short cross-sections
Mid slope	Shoulder	(α, β)	≤0	>0	ND	Convex element
	Nose slope		>0	>0	ND	Convex slope
	Head slope		<0	<0	ND	Concave slope
	Back slope		≥0	<0	ND	Concave element
	Side slope		(−0.2,0.2)	(−0.2,0.2)	ND	Rectilinear slope
Lower slope	Valley	[0, α]	ND	ND	>2	Lower element receiving water
	Lower flat		ND	ND	≤2	Flat element in the lower topographic position

Note: α and β are the break values for the upper, mid and lower slope classifications; RPI is relative position index.

). The first level of the scheme focuses on the geographic meanings via considering the spatial context along a downslope profile and classifies the slope by the regional terrain attribute. At the second level, the classification mainly emphasizes the morphological features of terrain positions and adopts the local terrain attributes to refine the first-level classification results.

As a regional terrain attribute, RPI, defined as the ratio between the distance from the point to the nearest valley and the total distance from the point to the nearest valley and ridge, has been widely applied to obtain a relatively ideal quantification of spatial transitions or spatial gradations between slope positions from the top to the bottom [1,31,32,33]. Hence, RPI is selected to classify a slope into three terrain positions, i.e., upper-, mid- and lower-slope, at the first level of terrain position classification. This taxonomy is similar to the highest level of the hierarchical system of landform classification as proposed by MacMillan et al. [16] or Drăguţ et al. [25]. It is worth pointing out that ridges and valleys identified previously for computing RPI can also be employed as the prototypes for ridges and valleys at the second level, and the first-level rules are defined by RPI to be consistent with the generation of other classes. At the second level, we conduct fine-grained classifications on the upper, mid and lower slopes, respectively. For the upper slopes, we may find that there are two types of terrain positions with significantly different morphologies, which include large parts and small parts. The cross-sections of the relatively large parts are mostly strip-shaped and features with ridge-like distributions. In contrast, the covering area of the relatively small parts is small, and their cross-sections are much shorter and exhibit discrete distributions. These morphological variances are detailed later on in Section 3.2. Adopting the definition of Speight et al. [8], we refer to the aforementioned large and small terrain positions as ridges and hillocks, respectively. As there are significant morphological differences between ridges and hillocks, we select two indexes based on the area and aspect ratio to reclassify the upper slopes. The mid slopes have slopes of more complex landforms. Accordingly, we introduce two local terrain attributes—plan curvature and profile curvature—to process the fine-grained classification. Utilizing the classification systems presented by Drăguţ et al. [25], the mid slopes are reclassified into shoulders, side slopes, nose slopes, head slopes and back slopes in term of different plan and profile curvatures. Herein, we replace the term of “negative contact” defined by Drăguţ et al. [25] by “back slope”. In fact, it is hard to find a terrain position with the plan and profile curvature both equal to 0. Hence, the fuzzy interval (−0.2, 0.2) is applied instead of the absolute value 0 for side slopes. The morphology of lower slopes is relatively simple and so the lower slope classes, namely valleys and lower flats, are simply differentiated according to whether the corresponding slope gradient is larger than 2° or not. The quantitative indicators for each terrain position are shown in Table 1.

2.3. Methodology for Classifying Terrain Positions

In this work, we integrate various terrain-position classification methods and design a multimodal classification approach as illustrated in Figure 2. Firstly, the regional and local terrain attributes are calculated from a DEM to provide the input for the following two-level classification. Secondly, combining the RA-based, rule-based and pixel-based classification methods, the slopes are classified into upper-, mid- and lower-slopes from the top to the bottom. Thirdly, by integrating the LA-based, supervised and object-oriented classification methods, ridges and hillocks are reclassified from the upper slopes. Meanwhile, the mid slopes are reclassified into shoulders, side slopes, nose slopes, head slopes and back slopes via comprehensive application of the LA-based, rule-based and object-oriented classification methods. Similar to classifying mid slopes, the lower slopes are reclassified into valleys and lower flats. Finally, the rationality assessment of the classification is performed by considering the distribution characteristics of terrain attributes at different terrain positions.

To resolve the fine terrain positions, different methods are adopted according to the characteristics of terrain positions and the data to be classified. As a regional terrain attribute, the RPI distribution is more continuous than other local attributes (e.g., slope and curvature) calculated via a fixed-size window. Therefore, the rule-based and pixel-based classification methods are comprehensively applied at the first level. Although the differences between ridges and hillocks may be denoted by the area and aspect ratio of each patch, it is difficult to define the classification thresholds of the indexes explicitly. Therefore, it is better to combine the supervised and object-oriented classification methods to reclassify the upper slopes at the second level. As the definite rules of local terrain attributes are used, the rule-based and object-oriented classification methods are employed to avoid fragment classification results when reclassifying the mid and lower slopes at the second level.

2.3.1. Calculation of Terrain Attributes

The terrain attributes, including RPI, slope, plan curvature and profile curvature, used for terrain position classification are derived from the DEM. The last three attributes are directly calculated using the ArcGIS software. To calculate the RPI, the ridges and valleys are needed. In order to enhance the accuracy and automation of the RPI calculation, several algorithms are chosen to calculate the corresponding elements. Firstly, the negative terrains relatively lower than the adjacent area are delineated using the algorithm presented by Yan et al. [34], which is implemented in a C++ and Matlab programming environment. Then, the ridges are extracted through the following three steps: (1) calculate initial ridge cells using flow accumulations estimated by the D8 algorithm of the ArcGIS software, which, after setting zero as the basic condition, returns 1 for ridge cells and 0 otherwise; (2) screen those initial ridge cells with values greater than 0.5 by using a 3 by 3 average window to filter discrete cells; and (3) determine final ridge cells by deleting the cells having negative terrains and the small ridge patches where their eight neighbors all are less than 4. Finally, when extracting the valleys, the flooding algorithm proposed by Rueda et al. [35] is adopted to identify the valley cells using the GDAL’s C++ library. The flooding algorithm is able to detect watercourses with a width greater than one cell. To implement the algorithm, the initial water depth of each cell is set to 100 mm and the simulation finished when the water transferred in an iteration fell below 1% of the total amount of water initially dropped on the DEM. After setting a threshold of drainage accumulations, the valley cells are reached following a clipping operation by negative terrains. Likewise, visual inspection by overlaying the terrain hillshade from the DEM is applied to further correct some valley parts such as wide river beds and dams. After extracting the ridges and valleys, we can calculate the distances of each grid cell from its location to the nearest valley and ridge and then acquire the corresponding RPI. In summary, the information related to the calculation of terrain attributes is listed in

Table 2. Information of the calculation of terrain attributes.

Terrain Attribute	Algorithm	Implementation
Slope	Third-order finite difference weighted by reciprocal of squared distance	ArcGIS software
Plan and profile curvature	Fourth-order polynomial	ArcGIS software
Negative terrain	Algorithm from Yan et al. [34]	Programming in C++ with Matlab
Ridge	D8 algorithm	ArcGIS software
Valley	Algorithm from Rueda et al. [35]	Programming in C++ with GDAL
Terrain hillshade	Single-directional hillshade algorithm	ArcGIS software

.

2.3.2. Determining the First-Level Parameters

The crucial task of the first-level classification is to choose suitable values for α and β. In fact, the critical values of RPI for the upper, mid and lower slope classifications are related to the topographic characteristics of the research area. Hence, the critical values are determined for the first-level classification by analyzing the relationship between RPI and the slope gradient. For the parameter α, the starting value of RPI is fixed from 0 with a step of 0.01, and then the average slope of each interval is counted in turn to generate the RPI-slope curve. For the parameter β, the RPI step is the same, but the ending value of RPI is fixed to 1. Generally, the RPI-slope curve for the parameter α moves from gentle to upward trend and then tends to be gentle, while the curve for the parameter β goes from gentle to downward and then to gentle again. As such, the parameter α and β can be determined at the change point from the gentle trend or decline to upward or gentle trend.

2.3.3. Determining the Second-Level Parameters

The second-level classification is to reclassify the upper, mid and lower slopes, in which the objects of the upper, mid and lower slopes are the primary elements. The upper slope objects (or patches) are grouped by the adjacent upper slope cells, and the mid and lower slope objects are generated by the multi-scale segmentation approach available in the eCognition software. Based on the upper slope patches, the feature space is constructed through selecting training samples and then the maximum likelihood method, a well-known and highly efficient supervised classification approach, is employed to reclassify the upper slope patches into ridges and hillocks. According to the classification rules composed of the indexes of plan and profile curvature and slope, the mid and lower slopes can be classified. It is worth pointing out that, in the process of mid slope classification, all objects are firstly classified into shoulders, nose slopes, head slopes and back slopes, and then the side slopes are further reclassified from the result of the four terrain positions.

• Upper slope patch acquisition

To distinguish the ridges and hillocks, the patches are obtained by merging the adjacent upper slope grid cells in the eight directions, followed by the calculation of the area and aspect ratio of each patch. The area is counted by multiplying the cell number of the patch by the pixel area, and the aspect ratio is calculated via the skeleton-line measurement. The skeleton-line measurement is achieved through three steps: (1) obtain the skeleton line of each patch using a thinning tool in the ArcGIS software; (2) calculate the shortest distance from each patch boundary to its skeleton line and take two times of the distance as the average width of the patch; and (3) further acquire the aspect ratio of each patch from the ratio of the average length and average width, where the average length is equal to its skeleton-line length.

• Multi-scale segmentation

As one of the current main-stream approaches for processing image data, the object-oriented classification partitions an image into multiple homogeneous segments and takes the segmented object as the classification unit. For the segmentation, the single-band or multi-band image provides the basic material. In this study, the data layers of RPI, plan curvature, profile curvature and slope are normalized to 0–255 for segmentation. Here, the normalization is to ensure the weight balance of the four layers and to develop a transferable classification system applicable to different datasets. As shown in Figure 2, the multi-band image composed of the four normalized layers is used for mid slope segmentation, while the single-band image from the normalized slope layer is applied for lower slope segmentation.

In this process, the multi-scale segmentation parameters determine the size and shape of the segmented objects. The optimal segmentation scale is a critical parameter. A relatively large segmentation scale will lead to the loss of relatively small objects, whilst a relatively small scale may cause the segmentation results to be over fragmented. In addition to segmentation scale, other parameters, such as weights of shape/color (with a sum of 1) and compactness/smoothness (with a sum of 1), will also have dominant effects on the segmentation results. The optimal segmentation parameters should exhibit an essential consistency between the segmented object and corresponding target object without any excessive or missing segmentation. To obtain the optimal segmentation parameters for mid and lower slope classification, we set the step sizes for both the weights of shape and compactness to 0.1, which begins from 0 to detect each optimal segmentation scale via the rate of change of local variance of object heterogeneity [36]. For every combination of these three parameters, the results of segmentation series are compared by visual inspection with overlaying the terrain hillshade to finally obtain the optimal segmentation parameters.

2.4. Rationality Assessment of Terrain Position Classification

It is crucial to assess the accuracy of terrain position classification. However, due to the lack of ground true data, quantitative evaluation is still challenging and no standard procedures exist [25,29,37]. In addition to the visual analysis with the terrain hillshade and 10 m contour, the rational evaluation is carried out from the following two aspects although we believe that more work needs to be done to improve it. Firstly, the average elevation and RPI of each terrain position are counted to evaluate the consistency of their geographic meanings. As illustrated in Figure 3, the elevation or RPI relationship among each type of terrain positions could be approximately described as follows: ridges > head slopes > “side” slopes > nose slopes > valley. It is worth noting that the “side” slopes, as shown in Figure 3, include shoulders, back slopes and side slope as defined in Table 1.

On the other hand, based on its definition, terrain position should have small intra-class variances and large inter-class differences. As a topographic factor that can accurately depict topographic changes and their influence on soil moisture distribution, topographic wetness index has been widely applied in different fields related to soil, hydrology and landform [39,40,41]. Thus, the classification rationality is also evaluated by measuring the intra- and inter-class differences of topographic wetness index for terrain positions. The topographic wetness index is calculated using the method reported by Wilson et al. [42]. The intra-class variance for each terrain position class is measured through the following steps: (1) divide the values of topographic wetness index over the whole area into 10 grades using the natural breaks method available in the ArcGIS software; (2) sort the first 10 patch samples from large to small according to the area of the terrain positon patch; (3) derive the frequency curve of topographic wetness index values for each sample; and (4) adopt a metric to quantify the variances between every two frequency curves and finally take the average to denote the intra-class variance. Through testing the most popular metrics of Manhattan distance, Euclidean distance and Pearson correlation coefficient, the Manhattan distance (Equation (1)) is employed due to its simplicity and efficiency. A very similar process can be applied to measure the inter-class difference between one terrain position class and other classes. The difference is that the frequency curve of each class is counted together in the 10 samples and the Manhattan distances are calculated from one class to the other 8 classes. The equation used to calculate the Manhattan distance is given as follows:

$$d(x,y)={\displaystyle \sum _{i=1}^N\left|x_i-y_i\right|}$$

(1)

where d(x,y) is the Manhattan distance, x and y are the one-dimensional frequency matrices corresponding to the different frequency curves of topographic wetness index values, and N is the number of elements in matrices x and y.

3. Results and Discussion

3.1. First-Level Classification

Following the multiple comparative analysis, we adopt 20,000 mm as the flow accumulation threshold to extract the channel network. After the extraction of ridges and valleys, the RPI distribution is finally obtained, as shown in Figure 4. As the regional terrain attribute, the continuity of RPI is found to be relatively good and exhibits gradual change from the upper to lower slope parts, which complies with the requirements of the first-level classification. The relationship between RPI and slope gradient is illustrated in Figure 5. Figure 5a shows that the average slope remains almost unchanged when RPI is less than 0.03; after that, it starts to increase rapidly. As depicted in Figure 5b, the variation curve of the average slope becomes gentle when the corresponding RPI exceeds 0.97. Therefore, 0.03 and 0.97 are selected as the break/critical values for the upper, mid and lower slope classifications. In another word, the values of α and β in Table 1 are 0.03 and 0.97, respectively. Accordingly, the terrain positions at the first level are acquired using the definite rules via the pixel-by-pixel calculation. As presented in Figure 6, the slope is divided into three terrain positions from the upper to lower parts, followed by the distinct boundaries of the terrain positions with almost no fragmentation. Conforming to the basic geo-cognition, the proportion of mid slope area is the largest and the proportion of upper slope area is the smallest. The above results indicate that it is reasonable to comprehensively apply the RA-based, rule-based and pixel-based classification methods to carry out terrain position classification at the first level.

3.2. Second-Level Classification

Figure 7 shows that the skeleton lines are located at the center of the patches, reflecting the patch shapes fairly well. This indicates that it is reasonable to adopt the skeleton-line method to calculate the geometric parameters of each patch (i.e., the average length and average width). It is worth noting that it will be difficult to obtain the skeleton line when the number of grid cells in a patch is relatively small, which leads to the relatively small aspect ratio. By analyzing the multiple experimental results, it is found that most of the patches with no skeleton lines belong to the hillocks. This is because that the area and average length of hillocks are both relatively small, which complies with the characteristics of patches for which the skeleton lines cannot be acquired. Therefore, before performing supervised classification to the upper slopes, we first reclassify those patches with no skeleton lines as hillocks. As necessary for the supervised classification, 30 ridge samples and 30 hillock samples, uniformly distributed in the study region, are chosen to train the rules based on the feature space constructed by the area and aspect ratio indexes. By using the maximum likelihood method, the classification of ridges and hillocks is completed (Figure 8). It is evident that the long ridge-like distributed patches with relatively large aspect ratio are reclassified as ridges, and those scattered small patches around the ridges are reclassified as hillocks, which is consistent with the characteristics described in Table 1. The experimental results demonstrate that, by combining the supervised and object-oriented classification methods, upper slopes can be reclassified fairly well into ridges and hillocks at the second-level classification.

When the object-oriented classification method is applied to reclassify the mid and lower slopes at the second level, the object segmentation should go first. Through the multiple comparative analysis, the optimal segmentation parameters (i.e., segmentation scale, weights of shape and compactness) are acquired for the mid and lower slope segmentations, which are 10, 0.3, and 0.5, and 6, 0.2 and 0.5, respectively. Based on these segmented objects, the mid and lower slopes are successively reclassified into the fine-grained items as defined in Table 1 via the specific rules. At last, the final classification results at the second level are shown in Figure 9. It is found that every type of terrain positions is not fragmented and has good integrity of their boundary, which complies well with the geographical cognition.

3.3. Rationality Analysis of Terrain Positions

Based on the visual analysis with the terrain hillshade and 10 m contour, it may be observed from Figure 9 that: (1) overall, the upper, mid, and lower slopes exhibit a distribution from the upper to lower parts of the slopes, and the veined lower slopes are cross-distributed with the upper slopes to constitute the skeleton of the terrain; (2) locally, in the upper slopes, the isolated bulges with the short cross-section are identified as hillocks and the relatively broad bulges with the strip-shaped cross-sections are reclassified as ridges. For the mid slopes, each fine-grained terrain position is distributed between upper and lower slopes. According to the definition of terrain position, the approximately inclined flat slopes are divided into side slopes; the lateral raised slopes are mostly reclassified as nose slopes; and the side concave slopes along upper slopes are mostly labelled as head slopes. Meanwhile, the relatively flat parts of lower slopes are reclassified as flats. Figure 10 shows the terrain positions classified by the pixel-based classification method using the same rules defined at the second level. The result of classification is too fragmented, which does not accord with the geographical cognition of terrain positions. Compared Figure 9 with Figure 10, the advantage of the method proposed in this paper is more obvious. Moreover,

Table 3. Average RPI and elevation for terrain positions.

	Shoulder	Nose Slope	Head Slope	Back Slope	Side Slope	Lower Flat	Valley	Ridge	Hillock
RPI	0.45	0.32	0.61	0.55	0.56	0.00	0.00	1.00	1.00
Elevation (m)	1003.71	981.77	1013.72	1002.63	1008.39	891.93	938.77	1034.75	1016.87

clearly lists the positional relationship of terrain positions as follows: ridges > hillocks > head slopes > side slopes ≈ shoulders ≈ back slopes > nose slopes > valleys > lower flats, which has a good consistency with their geographic meanings.

To further confirm the rationality of the classification results, according to the measurement of the intra- and inter-class differences by topographic wetness index for terrain positions, we compare the classification results with those obtained by the well-known classification method as proposed by Drăguţ et al. [25] (i.e., the Drăguţ method).

Table 4. Distances of intra- and inter-class for terrain positions.

Method in this Paper			Drăguţ Method
Terrain Position	Inter-Class Distance	Intra-Class Distance	Terrain Position	Inter-Class Distance	Intra-Class Distance
Ridge	0.90	0.38	Peak	0.94	0.58
Hillock	1.06	0.76	Steep slope	1.05	0.36
Shoulder	0.77	0.36	Shoulder	0.88	0.17
Side slope	1.06	0.53	Side slope	0.97	1.22
Nose slope	0.80	0.20	Nose slope	0.85	0.55
Head slope	0.75	0.24	Head slope	0.93	0.30
Back slope	0.73	0.11	Negative contact	0.97	0.17
Valley	1.01	0.26	Toe slope	1.51	1.15
Lower flat	1.88	0.98	Flat or gentle slope	1.65	1.34
Mean	1.00	0.42	Mean	1.08	0.65

lists the results of the comparative study. The current method gives the intra-class distance of terrain positions as 0.42 and the inter-class distance as 1.00, while the Drăguţ method produced the intra-class distance as 0.65 and the inter-class distance as 1.08. Comparatively, the terrain positions classified by the current method comply better with the classification criteria of small intra-class and large inter-class differences, and better reflect the topographic changes and their influence on the spatial distribution of soil moisture. For each type of terrain position, the intra-class distances of different terrain positions obtained using the current method are smaller and more stable, except for the flats and hillocks where the intra-class distances are all smaller than the intra-class mean of the Drăguţ method. The relatively largest intra-class distances of flats from both methods may be caused by the computational error of topographic wetness index in the flat or gentle slope areas. In addition, the stability of the two methods on quantifying the inter-class distances of different terrain positions is similar. However, the current method returns 77.78% of the differences between the inter- and intra-class distances exceeding 0.5, whilst the Drăguţ method only returns 44.44% and its intra-class distance unexpectedly appears to be greater than the inter-class distance. In general, the terrain position classification results obtained using the current method is satisfactory and the proposed method out-performs the Drăguţ method in several aspects.

4. Conclusions

According to the experimental results obtained from a 5 m gridded DEM in a loess hilly-gully region in north-west China, the following conclusions can be drawn:

(1)

The positional relationship of terrain positions has a good consistency with their geographic meanings. The inter-class difference of terrain positions is relatively large, and the intra-class variances are relatively small. Both above viewpoints indicate that terrain positions are consistent with the actual topography from both overall and local perspectives.

(2)

The two-level definition scheme effectively complies with the geographical cognition and topographic features of terrain positions from different levels, and sufficiently matches with the multimodal classification.

In summary, the new multimodal classification system can integrate the advantages of various classification methods effectively in order to identify terrain positions with relatively better integrity and rationality. However, it should be also noted that the classification method proposed in this work involves more parameters (e.g., the threshold of flow accumulations, RPI break values and segmentation parameters), which may affect the accuracy of the classification results. In future, a further research effort will be made to estimate these parameters for regions of different landforms to give guidance of selecting the parameter values. In addition, the automation and computation efficiency should also be further optimized.