Quantification of Loess Landforms from Three-Dimensional Landscape Pattern Perspective by Using DEMs

Quantitative analysis of the differences and the exploration of the evolution models of different loess landform types are greatly important to the in-depth understanding of the evolution process and mechanism of the loess landforms. In this research, several typical loess landform areas in the Chinese Loess Plateau were selected, and the object-oriented image analysis (OBIA) method was employed to identify the basic loess landform types. Three-dimensional (3D) landscape pattern indices were introduced on this foundation to measure the morphological and structural features of individual loess landform objects in more detail. Compared with the traditional two-dimensional (2D) landscape pattern indices, the indices consider the topographic features, thereby providing more vertical topographic information. Furthermore, the evolution modes between different loess landform types were discussed. Results show that the OBIA method achieved satisfying classification results with an overall accuracy of 88.12%. There are evident differences in quantitative morphological indicators among loess landform types, especially in indicators such as total length of edge, mean patch size, landscape shape index, and edge dimension index. Meanwhile, significant differences are also found in the combination of loess landform types corresponding to different landform development stages. The degree of surface erosion became increasingly significant as loess landforms developed, loess tableland area rapidly reduced or even vanished, and the dominant loess landform types changed to loess ridge and loess hill. Hence, in the reconstruction and management of the Loess Plateau, the loess tableland should be the key protected loess landform type. These preliminary results are helpful to further understand the development process of loess landforms and provide a certain reference for regional soil and water conservation.


Introduction
Landforms are comprehensive systems that show the shape, origin, distribution, and evolution of the earth's surface [1,2]. The morphology of landforms is the outward manifestation of the interaction of the earth's internal and external forces with the surface material under certain time conditions. Expressing the differences between various landforms through morphological description is essential to understanding the geographical environment and evolution [3,4]. Therefore, studying the differences of landforms has always been a crucial research content of geography [5,6]. Among various landforms worldwide, loess landforms have attracted much attention due to their particular forms. Loess landforms are complex and diverse combinations formed by following unique development modes based on the underlying geological strata and under the joint action of soil erosion [7]. The loess landforms are composed of loess tablelands, loess ridges, loess hills, and crisscross loess gullies. Among them, loess tablelands, loess ridges, and loess hills are the reflections at different stages in the development and evolution process of loess landforms [8]. They are direct evidence of the development process of loess landforms while hiding the loess gully erosion laws. Studying the differences in individual morphology and population distribution of these landform types helps to understand the spatial differentiation pattern of the loess landforms. More importantly, it is conducive to speculate the development process of loess landforms and provides the basis for the engineering control of soil erosion in the Chinese Loess Plateau. For a long time, people have been committed to improving the methods of describing and expressing landforms to satisfy the requirements of practical applications, such as textual description, pictographic mapping, and contour maps [9]. Nevertheless, the previous methods are mainly for finer visualization of the landforms, and most of the descriptions of the landforms are qualitative [10,11]. Subsequently, with the development of computer technology and GIS, accurate and quantitative descriptions of landforms gradually became possible. Currently, the research on loess landforms mainly focuses on four aspects, namely (1) classification methods of loess landforms [12][13][14], (2) extraction and analysis of terrain derivatives of loess landforms [15,16], (3) analysis on development and evolution process of loess landforms [5,7,17], and (4) construction of quantitative indices of loess landforms [18][19][20]. More specifically, regarding the classification methods of loess landforms, the current classification methods can be divided into two categories: pixel-based and object-based landform classification [13,21]. Despite the long-term dominance of pixel-based image analysis methods in remote sensing image processing, object-based image analysis is becoming increasingly popular [22]. The object-oriented image analysis (OBIA) method and deep learning (DL) classification method are typical representatives of these two categories and have been widely used in landform classification researches. Nevertheless, these studies have mainly focused on sample areas or watersheds at small scales [12,23,24] and therefore have limited applicability and guidance for large-scale regions. In addition, some studies have compared the relative performance of different classification algorithms using pixel-based and/or object-based image analysis methods without reaching a uniform conclusion [25][26][27]. In the terrain derivative extraction of loess landforms, researchers have made great efforts in the algorithm and proposal of terrain derivatives. More than 100 terrain derivatives have been proposed to date, and these derivatives are sufficient to describe the topographic features and differences of the loess landforms in detail. In the analysis of the evolution and development process of the loess landforms, different methods have been applied to the prediction and simulation of the evolution process of the loess landforms development, including geological dating [28], feature expression method [29], space for time theory [30], subbasin monitoring simulation method [31]. Scholars strive to realize the scientific cognition of "past-modern-future" in the development and evolution of loess landforms based on these methods and approaches [32]. Indeed, the research on quantitative indicators of loess landforms, as a fundamental approach to evaluate the landform differences and explore the landform development process, has attracted significant attention from scholars. In addition, many quantitative indicators were proposed and applied to the analysis of loess landforms.
In the quantitative research of loess landforms, mathematical morphological indices and terrain factors extracted from DEMs are widely used indicators. Scholars have widely used these types of indicators in loess landforms classification [33,34], loess soil erosion research [35,36], and loess landslide susceptibility assessment [37,38]. However, mathematical morphological indicators can reflect the morphological characteristics of individuals but are unable to express the aspects of spatial distribution and topological relationships among different landform types. Additionally, few morphological indices are available to support the detailed description of landform elements. Topographical factors can reflect the topographical characteristics of the loess landform. Nevertheless, topographic factors pay more attention to local features due to their calculation principles while ignoring the overall spatial relationship [14]. Terrain factors are primarily obtained through grid-by-grid calculation or sliding window algorithm, which causes the problem of "short-sightedness" in the analysis process. Thus, the landform is challenging to describe from a macro perspective [39]. Moreover, this research introduces the 3D landscape pattern indices, prevalent indicators in landscape ecology, to solve this bottleneck problem.
Three-dimensional landscape pattern indices are highly concentrated information of landscape patterns, reflecting the structural composition and spatial configuration characteristics of landscape patterns [40,41]. The 3D landscape pattern indices are calculated by improving the patch area and perimeter algorithm based on the 2D landscape pattern indices [42]. The traditional landscape pattern indices only consider the two-dimensional structure of the landscape, and its research objects are derived from the bird's eye view of the landscape, thus ignoring the topographic characteristics and the differences in the vertical direction of the research objects [43]. Especially in mountainous areas with complex terrain, the patch area and circumference obtained using 2D plane information are much smaller than the actual patch surface area and surface circumference, leading to inaccurate calculation results and weakening the differences between the research objects [44]. As a vital attribute of landform patterns, topographic features must be considered when quantifying landform differences. Therefore, the 3D landscape pattern indices are derived by substituting the surface area and surface edge length for the area and edge length in the 2D landscape indices [45]. Compared with the 2D landscape pattern indices, the 3D landscape pattern indices consider topographic information and express richer facade landscape information. As a unique landscape pattern, the introduction of 3D landscape pattern indices can provide more quantitative indicators for analyzing spatiotemporal heterogeneity of landform patterns. On this basis, the 3D landscape pattern indices are adopted in this paper to quantify the differences of the loess landforms.
In summary, despite the fact that previous studies have made many explorations in the classification and quantification of loess landforms and the analysis of the evolution model, there are still some limitations. Firstly, although various methods have been successfully implemented for the information extraction of loess landforms, the small size of the sample area and the identification of specific loess landform types have led to a lack of clarity in the identification of loess landforms on a large regional scale. On the other hand, studies on the quantification of yellow landforms have focused on topographic and individual morphological features but neglected the description of spatial distribution features, and the quantification indicators are few and simple. These deficiencies limit the clear understanding of the development and evolution of the loess landforms. To fill this gap, more in-depth researches are still required. Therefore, in this contribution, four typical sample areas of different zones in the Chinese Loess Plateau were selected to verify the applicability of the OBIA method for loess landforms classification at large scales. Meanwhile, the landscape pattern indices have been innovatively implemented into the loess landforms analysis, hoping to provide a new index system for the quantitative analysis of the geomorphology. Additionally, an attempt has been made to discuss the developmental processes of loess landforms based on the variability of the different loess landform types. This research aims to provide an essential reference for further understanding the development and evolution of loess landforms and regional soil and water conservation.
The objectives of this study are to (1) classify the loess landforms in the selected sample areas and evaluate the performance of the OBIA and DL methods, (2) introduce the three-dimensional landscape pattern indices to construct a quantitative index system of loess landforms, and aim to provide new index options for the quantitative expression of landforms, and (3) exploring the evolution process of the loess landforms based on the morphological changes of the loess landform types, providing certain support for related theories.

Study Area
The Jinghe Basin is one of the most important basins in the Chinese Loess Plateau; it is in the hinterland of the Loess Plateau and has typical loess landform features. This basin is located at the junction of the Shaanxi, Gansu, and Ningxia provinces (106 • 14 ~108 • 42 E, 34 • 46 ~37 • 19 N), covering a total area of 45,421 km 2 ( Figure 1). With the gradual change of the topography from the northwest to the southeast, the loess landform types also show apparent regional differences. The basin has a typical temperate continental climate, with annual precipitation of 350-600 mm [46]. The uneven rainfall distribution in time and space is reflected in the gradual decrease in space from south to north. The precipitation in May-September accounts for more than 70% of the annual rainfall. The yearly average temperature of the basin is 8 • C, and the spatial distribution of temperature also presents regional differences between high in the south and low in the north [47]. This basin is an important sediment export source of the Yellow River due to its severe surface erosion and soil loss [48]. Typical loess landform types, such as loess tableland, loess ridge, and loess hill, are distributed in this area, and evident regional differences are found in the north-south direction. It is an ideal sample area for studying the erosion process and spatial differences of loess landforms due to the differences and hierarchical distribution of landform types in this watershed.

Study Area
The Jinghe Basin is one of the most important basins in the Chinese Loess Plateau; it is in the hinterland of the Loess Plateau and has typical loess landform features. This basin is located at the junction of the Shaanxi, Gansu, and Ningxia provinces (106°14 '～ 108 ° 42' E, 34° 46 '～ 37° 19' N), covering a total area of 45421 km 2 ( Figure 1). With the gradual change of the topography from the northwest to the southeast, the loess landform types also show apparent regional differences. The basin has a typical temperate continental climate, with annual precipitation of 350-600 mm [46]. The uneven rainfall distribution in time and space is reflected in the gradual decrease in space from south to north. The precipitation in May-September accounts for more than 70% of the annual rainfall. The yearly average temperature of the basin is 8 ℃, and the spatial distribution of temperature also presents regional differences between high in the south and low in the north [47]. This basin is an important sediment export source of the Yellow River due to its severe surface erosion and soil loss [48]. Typical loess landform types, such as loess tableland, loess ridge, and loess hill, are distributed in this area, and evident regional differences are found in the north-south direction. It is an ideal sample area for studying the erosion process and spatial differences of loess landforms due to the differences and hierarchical distribution of landform types in this watershed.

Data Sources
The data used in the study mainly include DEMs data, remote sensing image data, hydrological data, and topographic factors derived from DEMs data ( Table 1). DEMs data include the 12.5 m resolution of the Advanced Land Observing Satellite (ALOS) downloaded from the National Aeronautics and Space Administration (NASA) official website. ALOS is currently one of the largest earth observation satellites, and the collected elevation data can achieve good accuracy in both horizontal and vertical directions. Compared with the 30 m resolution ASTER Global Digital Elevation Model (GDEMV2), ALOS can reflect more refined terrain texture information. The terrain features and contour lines at this resolution can be observed more clearly. The image data are obtained from Google Earth's 17-level satellite remote sensing images with a resolution of 2.38 m. As the primary source of landform classification, the data can provide rich image spectral information and texture information. Meanwhile, some terrain factors are calculated and added to the landform classification process to provide more details about the terrain texture [49]. River vector data are an essential reference for auxiliary image segmentation.

. Loess Landform Types Extraction
The OBIA and DL methods mentioned above are two approaches widely used in landform classification [12,50]. The OBIA method was adopted in this research for landform classification since an individual landform object has more geomorphological meaning than a pixel and more conforms to the actual situation of the loess surface. Furthermore, the classification experiment based on the DL classification method was also implemented as a reference to analyze the difference between the object-oriented classification method and pixel-based classification method. The overall loess landforms classification flowchart is shown in Figure 2. The OBIA classification process can be divided into three parts: image segmentation, the establishment of interpretation rule sets, and the post-processing of classification [51]. In this research, the commonly used multiresolution segmentation (MRS) method based on the eCognition platform (manufacturer: Trimble; version: eCognition Developer 9.0) was adopted for image segmentation.
The principle of MRS is to satisfy the conditions of minimum average heterogeneity of pixels and maximum homogeneity between pixels in the object; this objective is realized by combining adjacent pixels or small segmentation objects [52]. Multiresolution segmentation is a bottom-up region merging technology starting from a pixel object [53,54]. Small image objects can be merged into larger objects. Adjacent image objects are merged if they satisfy the defined minimum growth criteria for heterogeneity, otherwise the merging process stops. The image object generated after the segmentation process is the smallest unit for further image interpretation. The specific algorithm of MRS is as follows [55]: represents the spectral standard deviation of object 2 before merging, w cmpt is the weight of compactness, h cmpt is compactness, and h smo is smooth.
where l merge represents the perimeter of the merged object, b merge represents the perimeter of the circumscribed rectangle of the merged object, l obj1 and l obj2 represent the perimeter of object 1 and object 2 before merging, and b obj1 and b obj2 represent the perimeter of the circumscribed rectangle of the object 1 and object 2 before merging.  The principle of MRS is to satisfy the conditions of minimum average heterogeneity of pixels and maximum homogeneity between pixels in the object; this objective is realized by combining adjacent pixels or small segmentation objects [52]. Multiresolution segmentation is a bottom-up region merging technology starting from a pixel object [53,54]. Small image objects can be merged into larger objects. Adjacent image objects are merged if they satisfy the defined minimum growth criteria for heterogeneity, otherwise the merging process stops. The image object generated after the segmentation process is the smallest unit for further image interpretation. The specific algorithm of MRS is as follows [55]: According to the formula, the factors influencing the multiresolution segmentation effect in eCognition software include the weight of different bands, the segmentation scale parameter, shape weight, and compactness weight. When applying MRS, considering the strong control effect of the terrain factor on the segmentation process, the weight of the terrain feature is set to 2, whereas the weight of the spectral feature is set to 1. The color and shape factors are set to 0.4 and 0.6, respectively, and compactness and smoothness parameters are set at 0.5. The optimal segmentation scale is determined by the estimation scale parameter (ESP) tool. ESP is mainly based on the concept of local variance. The local variance of a specific scale is automatically calculated as the average standard deviation of all image objects under its scale. The optimal segmentation scale is selected based on the curve of local variance and its scale parameters [56]. When the local variance reaches the extremum and the change rate curve turns, the spatial heterogeneity of all image objects in this scale is the largest. The optimal scale of image segmentation is finally determined by analyzing and comparing different extremum points [57].
The DL classification model used in this paper was trained by Li et al. and has been successfully applied to the loess landforms. In this model, a widely used structure for a deep learning network, U-Net (Figure 3), is the fundamental network structure to execute the classification process due to its effectiveness. We optimized the U-Net network structure to make it more suitable for geomorphological research. A more effective multi-channel feature fusion structure was proposed to replace the single-channel input of the traditional U-Net network. Moreover, the concatenation between the contracting and expanding paths has also been optimized. In the optimized format, the channel of remote images was selected to connect with corresponding layers of the expanding path. The optimized network structure is proven to obtain better classification accuracy in the landform classification of the Chinese Loess Plateau. A more detailed and supported introduction to the adopted DL classification model and its training process could be found in reference [34]. object 1 and object 2 before merging, and bobj1 and bobj2 represent the perimeter of the circumscribed rectangle of the object 1 and object 2 before merging.
According to the formula, the factors influencing the multiresolution segmentation effect in eCognition software include the weight of different bands, the segmentation scale parameter, shape weight, and compactness weight. When applying MRS, considering the strong control effect of the terrain factor on the segmentation process, the weight of the terrain feature is set to 2, whereas the weight of the spectral feature is set to 1. The color and shape factors are set to 0.4 and 0.6, respectively, and compactness and smoothness parameters are set at 0.5. The optimal segmentation scale is determined by the estimation scale parameter (ESP) tool. ESP is mainly based on the concept of local variance. The local variance of a specific scale is automatically calculated as the average standard deviation of all image objects under its scale. The optimal segmentation scale is selected based on the curve of local variance and its scale parameters [56]. When the local variance reaches the extremum and the change rate curve turns, the spatial heterogeneity of all image objects in this scale is the largest. The optimal scale of image segmentation is finally determined by analyzing and comparing different extremum points [57].
The DL classification model used in this paper was trained by Li et al. and has been successfully applied to the loess landforms. In this model, a widely used structure for a deep learning network, U-Net (Figure 3), is the fundamental network structure to execute the classification process due to its effectiveness. We optimized the U-Net network structure to make it more suitable for geomorphological research. A more effective multi-channel feature fusion structure was proposed to replace the single-channel input of the traditional U-Net network. Moreover, the concatenation between the contracting and expanding paths has also been optimized. In the optimized format, the channel of remote images was selected to connect with corresponding layers of the expanding path. The optimized network structure is proven to obtain better classification accuracy in the landform classification of the Chinese Loess Plateau. A more detailed and supported introduction to the adopted DL classification model and its training process could be found in reference [34].

3D Landscape Pattern Indices.
Accurately obtaining the surface area and surface edge length of landscape patches is particularly important in calculating 3D landscape indices. This study adopts the calculation method based on DEMs data developed by Jenness [58]. This method calculates the surface area of each grid cell through the sliding frame algorithm and the triangle method. Each triangle is connected in 3D space to the center of the grid to be computed as well as the center points of the two neighboring grids (Figure 4). As a result, the Pythagorean theorem may be used to determine the triangle's surface edge length and the area of each

3D Landscape Pattern Indices
Accurately obtaining the surface area and surface edge length of landscape patches is particularly important in calculating 3D landscape indices. This study adopts the calculation method based on DEMs data developed by Jenness [58]. This method calculates the surface area of each grid cell through the sliding frame algorithm and the triangle method. Each triangle is connected in 3D space to the center of the grid to be computed as well as the center points of the two neighboring grids (Figure 4). As a result, the Pythagorean theorem may be used to determine the triangle's surface edge length and the area of each triangle. Finally, the surface area of the central grid is calculated using eight triangles. The calculation formula is as follows: where a s is the surface area, d s is the surface edge length, c is the cell size, S i is the slope of pixel i, a i is the relative altitude difference between the center point of the pixel i and the center point of the adjacent pixel, and 2c 2 is the square of the distance between the centers of two neighboring pixels.
triangle. Finally, the surface area of the central grid is calculated using eight triangles. The calculation formula is as follows: where as is the surface area, ds is the surface edge length, c is the cell size, Si is the slope of pixel i, ai is the relative altitude difference between the center point of the pixel i and the center point of the adjacent pixel, and 2c 2 is the square of the distance between the centers of two neighboring pixels. According to the research requirements, different representative indices are selected from the morphological and landscape levels to analyze the difference of loess landform types. The selected indices are shown in Table 2.

Total length of edge (TE) =
Ei is the surface edge length of patch i, n is the total number of patches of landform type This indicator can reflect the total edge length of a particular landform type (unit: km).
Mean patch size (MPS)

= × 10
Ai is the total surface area of loess landform i, and Ni is the total number of patches of landform type i.
This indicator can reflect a specific landform type's average patch area size (unit: pi is the surface circumference of patch i, ai is the surface area of patch i, n is the total number of patches i.
The value range of MPFD is between 1-2, and the larger the value, the more complex the shape of the patch. According to the research requirements, different representative indices are selected from the morphological and landscape levels to analyze the difference of loess landform types. The selected indices are shown in Table 2.
the total surface area of loess landform i, and N i is the total number of patches of landform type i.
This indicator can reflect a specific landform type's average patch area size (unit: km 2 ).

Mean patch fractal dimension (MPFD)
ln(a i ) n p i is the surface circumference of patch i, a i is the surface area of patch i, n is the total number of patches i.
The value range of MPFD is between 1-2, and the larger the value, the more complex the shape of the patch.
Landscape shape index (LSI) E is the total length of the surface edge of a certain loess landform type, and A i is the total surface area.
This indicator reflects the elongation of the patch. The larger the value of LSI, the longer the shape of the patch.
Circularity index (CI) This index is to quantify the shape of the patch. The closer the value is to 1, the closer the patch is to circle.

Slope Calculated in ArcGIS platform
The index represents the topographic features of the patch ( • ).
Edge dimension index (EDI) EDI = 2ln p i ln a i p i is the surface circumference of patch i, a i is the surface area of patch i.
This index reflects the complexity of the shape of the patch.

Extracted Loess Landform Types
Remote sensing image segmentation is carried out through the multiresolution segmentation tool in eCognition. The ESP tool is employed to test the optimal segmentation scale. The test results of the four sample areas under different segmentation scales are shown in Figure 5. The most significant extreme point in the local variance change rate curve is chosen as an alternative. The ideal segmentation scales of the four plots are eventually determined to be 36, 40, 32, and 36, respectively, based on the visual assessment of the segmentation results. The segmentation results show that although the four sample areas are in the same subbasin and have the same DEMs and remote sensing image resolution, the optimal segmentation scales of the four sample areas are still different, but the difference is not significant. The local segmentation effect of the four sample areas ( Figure 6) shows that under the optimal segmentation scale, the regions with greater homogeneity in the images were segmented into one or several complete objects. The boundary of the loess landform types was well detected. Meanwhile, the vector boundary of the image object formed by segmentation is more consistent with the actual objects.
topographic features of the patch (°).
Edge dimension index (EDI) = 2 pi is the surface circumference of patch i, ai is the surface area of patch i.
This index reflects the complexity of the shape of the patch.

Extracted Loess Landform Types
Remote sensing image segmentation is carried out through the multiresolution segmentation tool in eCognition. The ESP tool is employed to test the optimal segmentation scale. The test results of the four sample areas under different segmentation scales are shown in Figure 5. The most significant extreme point in the local variance change rate curve is chosen as an alternative. The ideal segmentation scales of the four plots are eventually determined to be 36, 40, 32, and 36, respectively, based on the visual assessment of the segmentation results. The segmentation results show that although the four sample areas are in the same subbasin and have the same DEMs and remote sensing image resolution, the optimal segmentation scales of the four sample areas are still different, but the difference is not significant. The local segmentation effect of the four sample areas ( Figure  6) shows that under the optimal segmentation scale, the regions with greater homogeneity in the images were segmented into one or several complete objects. The boundary of the loess landform types was well detected. Meanwhile, the vector boundary of the image object formed by segmentation is more consistent with the actual objects.  The classification tree method based on the feature threshold is applied to classify loess landform types. The first level classification is to distinguish the loess gullies from other loess landform types. Due to the significantly low brightness value of the loess gullies on the remote images and much lower elevations than surrounding areas, the main body of the loess gullies can be determined based on the mean brightness and DEMs. The remaining part is a mixture of three types of landforms. Following the definition and description of these three loess landform types, different thresholds can be set on the area, slope, and shape characteristics to distinguish them. It should be emphasized that auxil- The classification tree method based on the feature threshold is applied to classify loess landform types. The first level classification is to distinguish the loess gullies from other loess landform types. Due to the significantly low brightness value of the loess gullies on the remote images and much lower elevations than surrounding areas, the main body of the loess gullies can be determined based on the mean brightness and DEMs. The remaining part is a mixture of three types of landforms. Following the definition and description of these three loess landform types, different thresholds can be set on the area, slope, and shape characteristics to distinguish them. It should be emphasized that auxiliary visual judgment is also necessary. The detailed classification parameter settings are shown in Table 3. Table 3. Loess landform classification rule set.

Landform Types Classification Rules Description
Loess tableland (a) Ses = 36, Bri > 134, DEMs > 1335, slo < 7 The loess tableland is a large and flat plain that remains after the original loess plain is cut by gullies. (d) Ses = 36, Bri > 319, DEMs > 1115, slo < 6 Loess ridge (a) Area < 10 km 2 , LWR > 3 The loess ridge is a long strip of high elevation with a long shape and a narrow width. The area is significantly smaller than the loess tableland and surrounded by loess gullies. Loess hill Area < 1 km 2 , LWR < 2.5, slo < 10 • The loess hill is a round, nearly circular loess mound, an independent patch with an area usually less than 1 km 2 .
Note: ses, segmentation scale; bri, brightness; slo, slope; LWR, length-width ratio. After the initial classification workflow is completed, the manually assisted postprocessing is also implemented as a supplement to improve the classification accuracy. The main post-processing includes the discrimination and classification of tiny patches, internal holes filling, and manual correction of image noise areas (e.g., cloud occlusion areas and spectral color uneven areas). The classification accuracy is obtained by checking the consistency between the OBIA classification results and the manually sketched classification results. Twelve plots of 4 km × 4 km from the four sample areas are selected to sketch the boundaries of the loess landform types manually. The principle of sample selection is that there must be at least two loess landform types in the plot area. The overall classification accuracy is calculated by counting the ratio of the correctly classified area to the total manually sketched area. The classification accuracy evaluation results of each sample area are shown in Table 4. The results of the classification accuracy evaluation show that the classification accuracy of each sample area ranges from 81.23% to 93.47%. The overall classification accuracy is 88.12%, which shows that the employed method has achieved satisfying classification results. The results of the accuracy evaluation show that sample d has the best classification accuracy (93.47%), followed by sample a (87.47%), and samples b and c have the relatively worst classification accuracy. The widely distributed loess ridges and hills in samples b and c caused the regional surface to be fragmented, making it challenging to identify the loess landform types. In contrast, the dominant landform type in sample area d is relatively complete loess tablelands. The boundary between landform types in this area is more obvious, making it easier to obtain higher classification accuracy. The final classification results of loess landform types in four sample areas are shown in Figure 7. A total of 20 loess tablelands, 109 loess ridges, and 126 loess hills were interpreted in the four sample areas, providing sufficient samples to analyze the difference of loess landform types. For better visualization, the loess gullies are set to be no color. The area proportions of different landform types in the four sample areas were calculated, and the results are shown in Figure 8. Evidently, the loess gully is the dominant landform type in the four sample areas. Loess tableland still exists in samples a and d, and the proportion of loess tableland in the two sample areas is 11.28% and 38.98%. Among the four sample areas, sample d has the largest loess tableland area, and the loess tableland is still the primary landform type, which is second only to loess gully. This finding indicates that this sample is the area with the weakest surface erosion. The erosion degree of sample a is slightly higher than that of sample d. The loess tablelands are divided into loess ridges, making the area ratio of loess ridge (31.17%) second only to loess gully, but uncut loess tablelands are still found. Sample areas b and c have only two landform types, namely, loess ridge and loess hill, indicating that the surface erosion in these sample areas is more serious. Among them, the proportion of the loess gully in sample c has exceeded 70%, indicating that this sample area is the most severely eroded. The area proportions of different landform types in the four sample areas were calculated, and the results are shown in Figure 8. Evidently, the loess gully is the dominant landform type in the four sample areas. Loess tableland still exists in samples a and d, and the proportion of loess tableland in the two sample areas is 11.28% and 38.98%. Among the four sample areas, sample d has the largest loess tableland area, and the loess tableland is still the primary landform type, which is second only to loess gully. This finding indicates that this sample is the area with the weakest surface erosion. The erosion degree of sample a is slightly higher than that of sample d. The loess tablelands are divided into loess ridges, making the area ratio of loess ridge (31.17%) second only to loess gully, but uncut loess tablelands are still found. Sample areas b and c have only two landform types, namely, loess ridge and loess hill, indicating that the surface erosion in these sample areas is more serious. Among them, the proportion of the loess gully in sample c has exceeded 70%, indicating that this sample area is the most severely eroded. ISPRS Int. J. Geo-Inf. 2021, 10, x FOR PEER REVIEW 13 of 22

Calculated Quantitative Indices
The 3D landscape pattern indices calculation method mentioned above was used to calculate the indices value of each patch separately, and the average values of the 3D landscape pattern indices of different landform types were counted. The results are shown in Table 5. The calculation results of the 3D landscape pattern indices indicate the evident morphological differences in patches of different loess landform types, leading to noticeable differences in the landscape pattern indices. The MPS of loess tableland is the largest, with an average area size of 63.31 km 2 , followed by the loess ridge (4.47 km 2 ). The MPS of the loess hill is much smaller than that of the loess tableland and loess ridge, with an average area of 0.16 km 2 . This phenomenon is caused by the continuous erosion of surface materials with the development of loess landforms. From the perspective of the TE, the loess tableland has the most significant value due to its larger patch area and longer edges,

Calculated Quantitative Indices
The 3D landscape pattern indices calculation method mentioned above was used to calculate the indices value of each patch separately, and the average values of the 3D landscape pattern indices of different landform types were counted. The results are shown in Table 5. The calculation results of the 3D landscape pattern indices indicate the evident morphological differences in patches of different loess landform types, leading to noticeable differences in the landscape pattern indices. The MPS of loess tableland is the largest, with an average area size of 63.31 km 2 , followed by the loess ridge (4.47 km 2 ). The MPS of the loess hill is much smaller than that of the loess tableland and loess ridge, with an average area of 0.16 km 2 . This phenomenon is caused by the continuous erosion of surface materials with the development of loess landforms. From the perspective of the TE, the loess tableland has the most significant value due to its larger patch area and longer edges, followed by loess ridge and loess hill. The MPFD and EDI indices of the three loess landform types are ranked as loess tableland, loess ridge, and loess hill, confirming that the patch shape complexity of loess tableland is the most complicated, followed by loess ridge and loess hill.
Theoretically, the loess tableland has the lowest degree of erosion; the shape of the loess tableland should be the most complete, and the boundary complexity should be less than that of the loess ridge. However, the calculation results show the opposite due to the long border of the loess tableland, which is one of the most vulnerable parts to erosion. Therefore, the boundary erosion effect of the loess tableland is more serious. The loess tableland is constantly eroded to form the residual loss tableland. The residual loess tableland is a transitional form from loess tableland to loess ridge, with a more complex shape. The loess tableland in the selected sample areas is dominated by the residual loess tableland. Thus, the shape complexity of the loess tableland is higher than that of the loess ridge. The shape of the loess hill is the simplest, mostly near-circular or elliptical; thus, the MPFD and EDI values of loess hill are the lowest. The LSI reflects the length-width ratio of the loess landform type patch, that is, the elongation of the patch. The loess ridge has the most significant LSI value due to its typical long strip shape, followed by the loess tableland and the loess hill. The CI reflects the similarity between the patch shape and the circle. Evidently, the shape of the loess hill is the closest to the circle, resulting in the lowest CI value. The continuous advancement of the gullies leads to the evolution mode of the loess ridge mainly extending longitudinally; thus, the CI value is greater than that of the loess tableland. The slope calculation results show that the slope of the loess tableland is the smallest (4.08 • ), followed by the loess ridge (5.25 • ) and the loess hill (6.95 • ), thereby indirectly indicating that the loess tableland is most suitable for human habitation and agricultural production. The loss of surface material leads to steeper terrain as the surface erosion intensifies, resulting in a gradual increase in the slope value.

Comparison of Loess Landform Classification Methods
To evaluate the performance of the commonly used OBIA method and DL method in extracting loess landform types, sample area d was selected as the target region for implementing classification experiments based on these two approaches. The parameters of the OBIA method have been described in detail above, and the DL model utilized was previously trained and successfully implemented in the relevant study [34]. The accuracy assessment method introduced above was employed to calculate the classification accuracy of the DL method. The classification results of the two approaches and their comparison are shown in Figure 9. The results show that the classification accuracy based on the OBIA method is 93.47%, while the value of the DL method is 89.33%. Under the DL approach, the classification results of the identical loess landform patches are more fractured (Figure 9a), and a complete landform object was divided into several parts, resulting in the separated portions being categorized as other loess landform types by mistake. This phenomenon was due to the processing flow of the DL method. The DL method initially cut the remote sensing image into a specific image size and then distinguished these images based on the previously trained model [34]. After the model determined all the segmented remote sensing images, they were spliced to obtain the final classification results. Therefore, the results of landform types classified by DL methods are often relatively fragmented and prone to misclassification and omission. In comparison, the OBIA method divides the landforms into a reasonably complete patch, with a high degree of overlap between the classification boundary and the actual boundary (Figure 9b), and the overall classification result was better than the DL method. According to the viewpoint proposed by Blaschke and Strobl based on currently available applications, the OBIA method is in many instances superior to the traditional pixel-based method because the analysis unit is transferred from pixel to meaningful object [59]. The OBIA method considers the spectral difference of remote sensing images and incorporates terrain features into the segmentation and classification process [60]. A geomorphic unit with minimal heterogeneity for different geomorphic types is formed to maintain the integrity of the objects by considering the image's mean square error value and continuously testing the segmentation scale [61]. Meanwhile, the introduction of terrain factors can overcome the disparity of regional topography and enhance the stability of the classification system [62].
ical landform types in China [12]. These disputes indicate that there are uncertainties in the classification of complex landforms. Factors such as classification method, scale selection, data resolution, and land cover will affect the classification accuracy [62]. We believe that our contribution is to validate the applicability of the OBIA method in large-scale loess landform classification. This finding provides a reference for the extraction of landform information and the classification of landform types. We also expect to apply this method to other landform categories in further research to gain a wider application prospect. The entrance of human activities is another important contributor to the high classification accuracy of the OBIA method. The changes in loess landforms caused by human activities have made the boundaries of different landform types more evident, such as the construction of terraces and cultivation. The emergence of artificial landform improves the low accuracy of landform edge classification to a certain extent and obtains a more precise image segmentation. When the results of the OBIA and the DL methods are superimposed, evident differences in the boundary of landform types were found (Figure 9c). The reason for this phenomenon is that the OBIA method judges the segmented independent landform units when classifying, and the DL method judges a single pixel. The OBIA method can identify edges based on image color differences, but the DL method fails to show good adjustable recognition capabilities when encountering areas with significant differences in color or elevation at the edge of the landform.
This research reflects that the OBIA approach showed superior performance than the DL approach in classifying large-scale loess landforms. This conclusion is consistent with the findings of Zhao et al. After comparing the classification performance of the OBIA approach with two machine learning methods, he found that the OBIA approach obtained the best results in feature selection and terrain recognition [26]. However, controversy also exists. The experimental results of Li et al. at the small-scale level concluded that the classification effect of the DL method is superior to the OBIA method [34]. A multi-modal DL model presented by Du et al. also achieved convincing performance in classifying six typical landform types in China [12]. These disputes indicate that there are uncertainties in the classification of complex landforms. Factors such as classification method, scale selection, data resolution, and land cover will affect the classification accuracy [62]. We believe that our contribution is to validate the applicability of the OBIA method in large-scale loess landform classification. This finding provides a reference for the extraction of landform information and the classification of landform types. We also expect to apply this method to other landform categories in further research to gain a wider application prospect.

Sensitivity of the 3D Landscape Pattern Indices to Topography
Existing studies show that the 3D landscape pattern indices perform significantly better than the 2D landscape pattern indices in expressing complex landscape patterns [63,64]. The most noticeable advantage of the 3D landscape pattern indices is to improve the measurement algorithm of the area and perimeter of the object in the traditional 2D landscape pattern indices, making the calculated value more reliable [58]. Despite the widespread application of 2D and 3D landscape pattern indices [65,66], studies that addressed the differences are still limited. Thus, the four sample areas' 3D and 2D landscape pattern indices were calculated and compared to explore their difference and sensitivity in expressing landscape pattern information. As shown in Table 6, a certain difference was found between 3D and 2D landscape pattern indices. Among the selected indicators, the index with evident differences includes MPS, TE, EDI, and LSI, indicating that these indices are more sensitive to topographical changes. In contrast, slope, CI, and MPFD are less susceptible to terrain information than these indices. The difference in the sensitivity of different indices to terrain information is attributed to the strengthening or weakening of the weight of terrain information in the calculation process. In addition, the calculation results of the 3D landscape pattern indices are slightly larger than the 2D landscape pattern indices except for the LSI index. The explanation for this result can be summarized as the selected sample areas have a significant topography, and the use of 2D area and perimeter underestimates the actual area and perimeter of the terrain type. The result is that the calculated 2D landscape index does not accurately reflect the actual value. After overlaying the terrain information, the area and perimeter values are corrected, resulting in the calculated 3D landscape pattern indices larger than the 2D landscape pattern indices. Nevertheless, when the area and circumference increase simultaneously, their ratio may decrease, resulting in a greater LSI index in 2D than in 3D. Exceptional cases also exist in flat areas with relatively low topographic relief (e.g., plain river network region) where the 3D area and perimeter are approximately equal to the 2D area and perimeter. Therefore, no remarkable difference exists between the calculated 2D and 3D landscape pattern indices. Moreover, the numerical differences between 3D and 2D landscape pattern indices are insignificant because the loess tableland and the loess ridge have a certain flat surface area, and only the terrain at the edge of the landform type fluctuates wildly. The terrain information has little influence on its area and perimeter because the loess hill is relatively small. The selection of scale is an essential factor influencing the differences between the 2D and 3D landscape pattern indices [67]. In a previous study, Hoechstetter selected an extent of 2000 × 2000 m in eastern Germany to test the sensitivity and explanatory power of the 3D landscape pattern indices and found that area and landscape shape index had minor differences between 2D and 3D space. This result was considered to be caused by the tiny fluctuations in the DEMs data [68]. Some scholars have pointed out that this difference will increase as DEMs fluctuations increase, especially in mountainous landscapes with large topography undulations [44]. This research provides supporting evidence for this viewpoint. Furthermore, it needs to be emphasized that the calculation method of the 3D landscape pattern indices employed in this study is not unique. Some other methods of expressing 3D pattern information have also been proposed, such as the 3D indicators based on the lidar point clouds data [69], the roughness index [70], and geometric approaches [71]. These studies help to reveal the influence of topographical fluctuations on landscape patterns or landscape processes, provide more accurate indicators selection for geomorphological research and landscape quantification, and have potential application prospects. We aim to demonstrate the advantages of the 3D landscape pattern indices in quantifying the landform differences through this case study and to inspire relevant scholars for their future research.

Analysis of Evolution Process between Loess Landform Types
Numerous studies have focused on the development and evolution of loess landforms [7,17,72]. There are also mature and acceptable theories on the evolution of the loess landforms, for example, space for time theory [30], loess gully profile combination [73], and geo-information Tupu theory [74]. Previous studies have provided some insights into the evolution of the loess landforms; for example, Xiong pointed out that the development of the loess landforms is closely related to the underlying bedrock, which directly affects the evolution of the loess landforms. The surface processes, such as water erosion and wind erosion, can only act on the surface loess, which is difficult to affect the underlying bedrock because the loess is directly accumulated on the underlying bedrock [7]. Tong extracted the morphological indicators of the loess tableland and established the morphological parameter indices system through the topographic analysis method, which provided a specific foundation for the study of the refinement of the landform [18]. These studies provide a valuable reference for the time calibration of the morphological sequence of loess landform types. However, these theories and explorations mostly speculate on the evolution process of the entire loess landform from a macroscopic holistic perspective. More detailed evidence is needed to supplement the transformation process between loess landform types, which is also the concern of this research.
To date, it is generally accepted that the evolution mode of the loess landforms is from the loess tableland to the loess ridge and finally forms the loess hill ( Figure 10). However, the characteristic changes of loess landform types during the development process have not been fully explained. This study attempts to deeply discuss the development process of different loess landform types based on the calculated 3D indicators. As can be seen, the main landform type is broad and complete loess tableland in the initial stage of loess landform development. The loess tableland has a large area and a simple shape, which can be cultivated and inhabited by human beings, and the edge of the tableland only has a small amount of water erosion (Figure 10a). The gullies began to invade the loess tableland due to large-scale tectonic action, and the residential areas and farmland began to gather inside the loess tableland. The original complete and large tableland was gradually divided to form residual loess tableland with the continuous extension of the gullies. After the end of the tectonic movement, the surface process is dominated by erosion, the gully gradually widens and deepens, and the slope of the edge of the residual loess tableland becomes increasingly steep. When the gully completely cuts through the residual loess tableland, the strip-shaped loess ridge is formed (Figure 10b). This conclusion can be confirmed by the calculation results of landscape pattern indices. Compared with the loess tableland, the CI index of the loess ridge is more prominent, and the MPFD and ED indices are decreased. This finding indicates that the shape of landform units becomes more uncomplicated. Meanwhile, the average slope of the loess ridge has increased compared with the loess tableland, indicating that the shape of the landform units is slenderer, and the edge erosion effect is more serious. Loess hill is the typical representative of the old stage of loess landforms' evolution. With the continuous action of water erosion, the edge erosion effect of the long loess ridge (loess residual tableland) is more remarkable. The shape also begins to shrink to the inside, and the area and perimeter decrease simultaneously. When the connection between the loess ridge and the loess tableland is completely broken, a relatively independent loess hill landform unit is formed (Figure 10c). Compared with the loess ridge, the CI value of the loess hill gradually approaches 1, and the shape slowly evolves into a nearly circular shape.
representative of the old stage of loess landforms' evolution. With the continuous action of water erosion, the edge erosion effect of the long loess ridge (loess residual tableland) is more remarkable. The shape also begins to shrink to the inside, and the area and perimeter decrease simultaneously. When the connection between the loess ridge and the loess tableland is completely broken, a relatively independent loess hill landform unit is formed (Figure 10c). Compared with the loess ridge, the CI value of the loess hill gradually approaches 1, and the shape slowly evolves into a nearly circular shape. Water flow erosion is the main factor driving the development of the loess landforms. Meanwhile, the surface erosion of the Loess Plateau is the main source of surface fragmentation and sedimentation of the Yellow River [75]. Therefore, several policies have been implemented to reduce the surface erosion of the Loess Plateau, such as gully consolidation and highland projection project and the conversion of cropland to forest policy [76]. As the main gathering area of human production and living activities, the loess tableland should be protected, especially the loess residual tableland [19]. This condition is attributed to the fact that the loess tableland is in the young stage of loess landforms development, with a large area, gentle terrain, slight soil erosion, and abundant nutrients in the surface soil; it is more conducive to large-scale human habitation and farming [77]. On the contrary, loess ridge and loess hill are not dominant in quantitative indicators such as shape, topography, and area. Therefore, the loess tableland is a geomorphic unit that humans rely on for survival. It is also a key area where the evolution from the loess tableland to the loess ridge occurred and needs to be specially protected.

Limitations and Future Research
Although the employed OBIA method achieves good results in the classification of loess landforms, further improvements are possible if other factors are considered. The primary consideration is the improvement of the data quality. The topographic Water flow erosion is the main factor driving the development of the loess landforms. Meanwhile, the surface erosion of the Loess Plateau is the main source of surface fragmentation and sedimentation of the Yellow River [75]. Therefore, several policies have been implemented to reduce the surface erosion of the Loess Plateau, such as gully consolidation and highland projection project and the conversion of cropland to forest policy [76]. As the main gathering area of human production and living activities, the loess tableland should be protected, especially the loess residual tableland [19]. This condition is attributed to the fact that the loess tableland is in the young stage of loess landforms development, with a large area, gentle terrain, slight soil erosion, and abundant nutrients in the surface soil; it is more conducive to large-scale human habitation and farming [77]. On the contrary, loess ridge and loess hill are not dominant in quantitative indicators such as shape, topography, and area. Therefore, the loess tableland is a geomorphic unit that humans rely on for survival. It is also a key area where the evolution from the loess tableland to the loess ridge occurred and needs to be specially protected.

Limitations and Future Research
Although the employed OBIA method achieves good results in the classification of loess landforms, further improvements are possible if other factors are considered. The primary consideration is the improvement of the data quality. The topographic information provided by the DEMs is an important foundation of the OBIA method but is limited by the spatial resolution of the DEMs data. Moreover, the limited vertical landscape information provided by low-resolution DEMs data in calculating 3D landscape pattern indices also causes errors in the calculation results. With the reduced acquisition costs of high-resolution DEMs, the widespread application of high-precision topographic data based on low-altitude aerial photogrammetry or airborne LiDAR technology has become possible. High-resolution DEMs data will improve the classification accuracy of the landforms, but the resolution selection will depend on the size of the target area and the field survey. According to the viewpoint of Tarolli, there is no perfect resolution but a suitable resolution [78]. The choice of DEMs data resolution is an issue to be considered in further research.
In addition, four randomly distributed sample areas were selected to obtain sufficient samples of different loess landform types, which is convincing for the study of individual characteristics of loess landform types. However, to better analyze the spatial distribution characteristics of the loess landforms, it is better to take the complete catchment as the study area because the catchment is the basic unit of landform evolution [26]. Furthermore, the description of the evolution process of the loess landforms is based on the speculation of the morphological calculation results, which is reverse reasoning from the result to the process. More details about the evolution process of the loess landforms still need to be explored. An effective method is to analyze the soil erosion process and evolution characteristics of the loess landform by monitoring the terrain changes in multiple periods, which is inseparable from the support of high-precision topographic data. Thus, exploring the loess landform evolution model based on catchment scale still needs further in-depth research.

Conclusions
The quantification and differentiation analysis of complicated loess landforms is an essential prerequisite for in-depth exploration of the pattern and evolution of loess landforms. In this contribution, we verified and compared the performance of the OBIA method in classifying loess landforms at a larger spatial scale. Meanwhile, 3D landscape pattern indices were introduced to analyze the differences between loess landform types quantitatively. Moreover, the conversion mode between loess landform types was also discussed. The results proved the excellent performance of the employed OBIA method in identifying loess landform types, with an overall accuracy of 88.12%. Compared with the DL method, the OBIA method can recognize more complete landform types, and the extracted landform boundaries are more consistent with the actual landform elements. Certainly, this is affected by the quality of remote sensing images, and high-resolution remote sensing images are required to obtain better classification accuracy. Although the classification accuracy of the OBIA method, in this case, is slightly higher than that of the DL method, this is not absolute because the classification accuracy is affected by many factors such as image quality, sample selection, and prior knowledge.
The 3D landscape pattern indices introduced by this research improve the shortcomings of traditional 2D landscape pattern indices that lack topographic information. The results prove that reasonably designed 3D landscape pattern indices are feasible to quantify the differences in loess landforms accurately. This finding provides more accurate indicators for the quantification of loess landform features and helps to better understand the spatial pattern and evolution of loess landforms, which has a great potential for broad application. The results of the quantitative analysis of loess landform types indicated that, in addition to the continuous loss of surface material, a series of indicators such as shape, area, complexity index, and fragmentation index also change regularly during the evolution of loess landforms and eventually lead to differentiated loess landform types and complicated spatial patterns. Compared with the existing theories of loess landform evolution, this contribution focuses more on the differences in loess landform objects and patterns and speculates on the developmental processes of loess landforms through the changes of loess landform entities at a more microscopic level, providing evidence for systematic theories of loess landform evolution.

Data Availability Statement:
The data presented in this study are available on reasonable request from the corresponding author.