An Improved Method for the Evaluation and Local Multi-Scale Optimization of the Automatic Extraction of Slope Units in Complex Terrains

: Slope units (SUs) are sub-watersheds bounded by ridge and valley lines. A slope unit reﬂects the physical relationship between landslides and geomorphological features and is especially useful for landslide sensitivity modeling. There have been signiﬁcant algorithmic advances in the automatic delineation of SUs. But the intrinsic difﬁculties of determining input parameters and correcting for unreasonable SUs have hindered their wide application. An improved method of the evaluation and local multi-scale optimization for the automatic extraction of SUs is proposed. The Sus’ groups more consistent with the topographic features were achieved through a stepwise approach from a global optimum to a local reﬁning. First, the preliminary subdivisions of multiple SUs were obtained based on the r.slopeunit software. The optimal subdivision scale was obtained by a collaborative evaluation approach capable of simultaneously measuring objective minimum discrepancies and seeking a global optimum. Second, under the selected optimal scale, unreasonable SUs such as over-subdivided slope units (OSSUs) and under-subdivided slope units (USSUs) were further distinguished. The local average similarity ( LS ) metric for each SU was designed based on calculating the SU’s area, common boundary and neighborhood variability. The inﬂection points of the cumulative frequency curve of LS were calculated as the distinguishing intervals for those unrealistic SUs by maximum interclass variance threshold. Third, a new effective optimization mechanism containing the re-subdivision of USSUs and merging of OSSUs was put into effect. We thus obtained SUs composed of terrain subdivisions with multiple scales, which is currently one of the few available methods for non-single scales. The statistical distributions of density, size and shapes demonstrate the excellent performance of the reﬁned SUs in capturing the variability of complex terrains. Beneﬁting from the sufﬁcient integrating approach of diverse features for each object, it is a signiﬁcant advantage that the processing object can be transferred from general entirety to each precise individual.


Introduction
Regional landslide risk assessment has been shown to be an effective tool for mitigating casualties and property damage caused by landslide hazards [1]. To obtain reliable assessments of landslide risks, it is crucial to first select an appropriate mapping unit for extracting and distributing geoenvironmental data [2]. Common mapping units mainly include grid units, slope units (SUs), terrain units and unique condition units [3]. SUs are associated with the geomorphological process shaping the natural morphology, and thus reflect the physical relationships between landslides and the environmental information [4]. Therefore, SUs are especially suitable for landslide sensitivity modeling [5,6]. The acquisition and application of SUs have received increased attention [7]. devoted to the continuous optimization of subsequent unrealistic objects after one "optimal" segmentation [29,30]. The identification and combination of multiple scales have been shown to be help resolve undesirable objects [31]. The concepts of over-subdivision and under-subdivision are equally applicable to classifying undesirable SUs. Unfortunately, the problem of continuous optimization for the subsequently unreasonable SUs has always been ignored by researchers, and the corresponding solutions are not tried.
To improve the production of objective and reproducible SUs and make up for the lack of subsequent amending for unrealistic SUs, an improved quality evaluation and local multi-scale optimization method for automatically extracting SUs in complex terrains is proposed. The optimal scale criterion for delineating SUs, and effective integration of diverse spatial features into refinements of unrealistic SUs are demonstrated in this paper. It is the first application for improvement of individual SUs' objects from automatic subdivision procedures and from which SUs consisted of terrain subdivisions with multiple scales are obtained.

Study Area
The Yuqu River Basin is at the southeastern margin of the Tibetan Plateau and is a tributary on the left bank of the Salween River (latitude/longitude: 28 • (Figure 1a). The region forms a typical high mountain and canyon landform. The total study area is 9190 km 2 , with a river drop of 2122 m (Figure 1b). It can be divided into a high hilly zone, a middle high-mountain zone, a gorge transition zone and an alpine gorge according to altitude and topography [32,33] (Figure 1c). The high hilly zone (I) alternates the distribution of an open intermountain basin and rounded mountains, of which the valley bottom width is 800~5000 m. The middle high-mountain zone (II) is a widely developed glacier erosion and accumulation landform, located in the northeast highlands at 4300~5800 m elevation. The gorge transition zone (III) is an erosional mountain with a "U" type river valley, with an average gradient of 0.69%. River floodplains and multi-level terraces are generally developed and have a valley width of 50~800 m. The alpine gorge zone (IV) is located downstream of the river and has an elevation of 1800~3200 m. The width of the deep "V" valley is only 40~80 m and the flow is turbulent, with an average gradient of 1.12%.

Data Preparation
The Advanced Land Observing Satellite-1 from the Japan Aerospace Exploration Agency (JAXA) provides ALOS-12m DEM data with a spatial resolution of 12.5 m ( Figure  2a). Remote sensing images were obtained from the Gaofen-2 satellite data acquired on 20 January 2020, and a true color image covering the Yuqu River Basin at 1-m spatial resolution was obtained (Figure 2b). The hillshade of the DEM and images unified into the projection coordinate system WGS 1984/UTM Zone 47N, were imported in ArcScene 10.6.1

Data Preparation
The Advanced Land Observing Satellite-1 from the Japan Aerospace Exploration Agency (JAXA) provides ALOS-12m DEM data with a spatial resolution of 12.5 m (Figure 2a). Remote sensing images were obtained from the Gaofen-2 satellite data acquired on 20 January 2020, and a true color image covering the Yuqu River Basin at 1-m spatial resolution was obtained (Figure 2b). The hillshade of the DEM and images unified into the projection coordinate system WGS 1984/UTM Zone 47N, were imported in Arc-Scene 10.6.1 to obtain the three-dimensional terrain representation. Actual comparison photos of typical locations are obtained by using an unmanned air vehicle (Phantom 4pro V2.0 by DJI Innovation Technology Co., Ltd., Shenzhen, China, flight altitude 500 m) in field investigations (Figure 2c).

Reference SUs
The reference SUs play a crucial role in the evaluation of optimal scale for the automatically extracted SUs. To obtain normative reference samples, the following procedures are performed: (1) The forward and reverse DEM hydrological analysis method is selected to obtain the sub-basins divided by drainage lines and dividing lines. The method conforms to researchers' original definition of SUs, and detailed and common rules can be referenced to obtain undisputed results [8][9][10]. Based on the ArcGIS Hydrological analysis tool, the extraction of ridgelines mainly includes depression filling, flow direction extracting, flow accumulation setting, river network extracting, obtaining catchment area [10]. The extraction of valley lines is obtained through reversing DEM, performing the same steps as for ridgelines [10,34]. According to previous studies on the adjacent regions [35,36], the appropriate catchment area thresholds are at the axis between river network density and catchment area [18,36]. As is shown in Figure 3b, the catchment area threshold of positive topography is 3000 pixels, and the negative topography is 3500 pixels. The filling thresholds are set as the recommended defaults (Table 1). (2) In order to pick out the single slopes and reject multiple slopes or even an entire watershed from the subbasins, the circular variance of aspect is taken as measure of homogeneity for each subregion. Referring to experience about the subdivision of SUs from Alvioli et al. [14,15] and Jacobs et al. [5], only the sub-basins of aspect circular variance ( c ) between 0.20~0.30 are selected as candidate reference SUs. (3) To improve the validity and typicality of those samples, regular hexagonal grid zones are created to implement an equally spaced sampling strategy. The abnormal parallel pseudo valleys (Figure 3a(B)) are abandoned and manually labelled as non-sampling zones, which is the limitation inherent in the hydro-

Reference SUs
The reference SUs play a crucial role in the evaluation of optimal scale for the automatically extracted SUs. To obtain normative reference samples, the following procedures are performed: (1) The forward and reverse DEM hydrological analysis method is selected to obtain the sub-basins divided by drainage lines and dividing lines. The method conforms to researchers' original definition of SUs, and detailed and common rules can be referenced to obtain undisputed results [8][9][10]. Based on the ArcGIS Hydrological analysis tool, the extraction of ridgelines mainly includes depression filling, flow direction extracting, flow accumulation setting, river network extracting, obtaining catchment area [10]. The extraction of valley lines is obtained through reversing DEM, performing the same steps as for ridgelines [10,34]. According to previous studies on the adjacent regions [35,36], the appropriate catchment area thresholds are at the axis between river network density and catchment area [18,36]. As is shown in Figure 3b, the catchment area threshold of positive topography is 3000 pixels, and the negative topography is 3500 pixels. The filling thresholds are set as the recommended defaults (Table 1). (2) In order to pick out the single slopes and reject multiple slopes or even an entire watershed from the sub-basins, the circular variance of aspect is taken as measure of homogeneity for each subregion. Referring to experience about the subdivision of SUs from Alvioli et al. [14,15] and Jacobs et al. [5], only the subbasins of aspect circular variance (c) between 0.20~0. 30  zones are created to implement an equally spaced sampling strategy. The abnormal parallel pseudo valleys (Figure 3a(B)) are abandoned and manually labelled as non-sampling zones, which is the limitation inherent in the hydrological analysis method. Taking each hexagon of sampling zones as the center (Figure 3a(A)), the candidate SUs with the largest size in the statistical area are selected as reference SUs. As is shown in Figure 3c,d, the sub-basins C with uniform aspect and proper size are more suitable as reference objects in that hexagon grid zone. (4) Based on the commonly used criteria for estimating the sampling quantity (the confidence level, 95%; the confidence interval, 3%) [37] and 8875 sub-basins as references for estimating the total number of SUs, an appropriate sampling size of reference SUs should be 952. Accordingly, the research area was divided into 1011 square hexagonal grids, 59 of which were marked as non-sampling zones, leaving 952 sampling zones to generate reference SUs (Figure 3a).    [14,15,39] Area The maximum area of candidate reference SU in each sampling zone [14]

Methodology
Given the challenges in the determination of the optimal scale and refining of unre alistic SUs, this study proposes an iterative process to achieve SUs more consistent with the terrain (

Methodology
Given the challenges in the determination of the optimal scale and refining of unrealistic SUs, this study proposes an iterative process to achieve SUs more consistent with the terrain (

Initial Subdivision of SUs
Due to its extensive use and demonstrated skill in capturing the aspect variability of the landscape [5,15,39], the r.slopeunits software v1.0 (Ivan Marchesini and Massimiliano Alvioli, Perugian, Italy) was adopted to obtain the multiple groups of the preliminary subdivision of SUs. Given a DEM and some input parameters, the algorithm first partitions the digital topography into several larger half basins (HBs); then, the sub-basins are continuously divided under the reduced factor (r) of the initial flow accumulation area threshold (t) in each iteration [39]. When the child of an HB matches the user-defined parameters of minimum area (a) and circular variance (c), it is selected as a candidate SU and the iterative procedure ends. The software of r.slopeunits v1.0 was obtained from the portal of the Geomorphology Research Group (http://geomorphology.irpi.cnr.it/tools/slopeunits (accessed on 10 June 2022)). Users need to define the following parameters: the initial flow accumulation area threshold (t), minimum surface area (a), minimum circular variance (c), reduction factor (r), threshold value for the cleaning procedures (Cleansize).

Determination of the Optimal Subdivision Scale
The ability to explicitly calculate the error between the automatically extracted SUs

Initial Subdivision of SUs
Due to its extensive use and demonstrated skill in capturing the aspect variability of the landscape [5,15,39], the r.slopeunits software v1.0 (Ivan Marchesini and Massimiliano Alvioli, Perugian, Italy) was adopted to obtain the multiple groups of the preliminary subdivision of SUs. Given a DEM and some input parameters, the algorithm first partitions the digital topography into several larger half basins (HBs); then, the sub-basins are continuously divided under the reduced factor (r) of the initial flow accumulation area threshold (t) in each iteration [39]. When the child of an HB matches the user-defined parameters of minimum area (a) and circular variance (c), it is selected as a candidate SU and the iterative procedure ends. The software of r.slopeunits v1.0 was obtained from the portal of the Geomorphology Research Group (http://geomorphology.irpi.cnr.it/tools/ slope-units (accessed on 10 June 2022)). Users need to define the following parameters: the initial flow accumulation area threshold (t), minimum surface area (a), minimum circular variance (c), reduction factor (r), threshold value for the cleaning procedures (Cleansize).

Object-Level Consistency Error (OCE)
The ability to explicitly calculate the error between the automatically extracted SUs and the reference SUs from the ground truth is of prior importance in our assessment framework. It is the first step for eliminating the unreasonable SUs result. The common error measures are constructed by regarding the assessment as a process of intersection (correct) or difference (wrong) pixel labeling [40,41]. Consequently, it cannot accurately distinguish the error originated from under-subdivision or over-subdivision. When the combination of multiple extracted SUs is exactly the same as reference SUs, the error measure would consider the two to be consistent. Therefore, the object level consistency error (OCE) was introduced to evaluate the performance of extracted SUs objectively [25]. The error measures are based on object-by-object comparisons of extracted SUs and reference SUs. Compared to existing error measures, it can take into account the size, shape, and position of each recorded SUs at the object level. Moreover, it is sensitive to both under-subdivision and over-subdivision, which contributes to make reasonable inferences from stepwise errors evolution analyses. The OCE is calculated as: where E g,s represents the part error of automatically extracted SUs to reference SUs, and E s,g represents the part error of reference SUs compared to automatically extracted SUs. The part error E g,s is measured as: where I g = {A 1 , A 2 , . . . , A M } is the reference SUs, and A j is the j th SU in I g ; I s = {B 1 , B 2 , . . . , B N } is the automatically extracted SUs, and B i is the i th SU in I s ; A j represents the number of grid cells in A; A j ∩ B i and A j ∪ B i denotes the intersection and combination of A j and B i , respectively; W ji weights each B i that intersects with A j according to the size of B i relative to all grid cells in I s that intersect with A j ; W j weights the importance of A j relative to I g ; δ(x) is the delta function. In addition, E g,s is calculated by replacing A j and B i in Equation (2). OCE is normalized between [0, 1], where 0 means completely consistent without error, 1 means completely mismatched. It is generally considered that OCE < 0.35 has good consistent matching [25,27]. The standardized acquisition of ground truth (reference) objects is fundamental to the proper functioning of OCE. The lack of explicit correspondence between historical landslides and extracted SUs is not recommended as reference objects. Experts should widely accept high quality reference SUs following detailed processing guidelines to eliminate subjectivity and enhance reproducibility.

Global Optimal Heterogeneity
A moderate scale of subdivision of the terrain is also required to pursue a global optimal optimum. All the SUs should maximize the internal and external heterogeneity. The SUs closest to this assumption should be selected as the global optimal scale. The Global Moran's Index (MI) and global variance (V) are adopted to make a straight- forward evaluation of external heterogeneity and internal homogeneity, respectively [42]. The two quantities are calculated by Equations (3) and (4): where m is the total number of SUs; w ih is an indicator for spatial proximity, whose value equals 1 when SU h and SU i share a common border, 0 otherwise; y i is the average aspect of SU i , y is the average aspect of the whole terrain aspect; s i and c i are the area and circular variance of SU i . Note that the angle is needed to convert to radians. The average values and the difference should be intended vectorially following the equation in Alvioli et al. [15].
The global heterogeneity score (GS) is calculated by Equation (5): The group of highest GS is considered to achieve the maximum balance of the internal homogeneity and external heterogeneity and identified as the best subdivision. However, as the candidate SU scale sets gradually increase, the maximum or minimum of MI and V may change accordingly and so does the highest GS value. To compensate for the uncertainty in confirming the optimal scale in existing methods, a combinative use of the OCE strategy is implemented. The measures of OCE are first performed to filter out unrealistic terrain subdivisions, and then qualified SU groups are left to participate in the calculation of GS values.

Identification of USSUs and OSSUs
Accurately distinguishing these unreasonable SUs is critical for subsequent providing refining objects. In this section, a new metric is designed to automatically detect those unrealistic SUs and categorize them as USSUs and OSSUs. Their specificities are summarized from field investigation and three-dimensional topographic analysis. From a general impression, the USSUs form a larger area than their adjacencies, and a loose boundary along the landform. They have stronger aspect differences between their adjacencies and low intra-unit similarity. Correspondingly, the OSSUs have smaller and more fragmentized sizes than normal units. They are usually mixed with surrounding SUs and share a relatively long common boundary. Thus, the OSSUs have higher intra-unit similarity and few differences in aspect with their adjacencies.
To compare the similarity and difference of each unit in a more reasonable way, we took the indirect value of the neighborhood variability of aspect difference or standard deviation rather than their direct values. The average value of all the adjacent SUs was taken as the reference standard of the normal level. The ratio of the SU value to the average value was the neighborhood variability. When the ratio was far from 1, the magnitude of the difference from the normal value was highlighted. Firstly, the interior similarity and adjacent differences for each SU were determined by the neighborhood variability of aspect difference (D a ) and the neighborhood variability of aspect standard deviation (D sd ), as shown in Equations (6) and (8), respectively.
where Y i and Y i is the aspect difference of SU i and the average value difference with its adjacent SUs, respectively; y i , w im are defined as in Equation (3); SD i and SD i is the standard deviation of SU i and the average value of its adjacent SUs, respectively. These SUs conflicted with the original terrain discontinuity, and had spatial features that were too large, too small, or oddly shaped. Secondly, to improve the accuracy of the identification, the diversity of the SUs's size and boundary were considered by introducing the area ratio (r a ) and the length ratio (r l ), respectively. The area ratio r a for each SU was calculated as in Equation (9), where a i and a m was the s area of SU i and its adjacent SUs, respectively. The length ratio r l was the maximum length of the common boundary to its perimeter, calculated using the Equation (10), indicating the degree of coincidence of the common boundary. Where p i and l i n were the perimeter of SU i and the length of common boundary of its adjacent SUs, respectively: Subsequently, the area, boundary and aspect difference were integrated into local heterogeneity he L , which was calculated with Equation (11). Similarly, the area and aspect uniformity were integrated into local homogeneity ho L and calculated with Equation (12).
The design structure of the formulas was ingenious. Taken he L for example, the form of the square power function could enhance specificity with D a larger than 1. The part of (r a i + 1 1+r l i ) was a combinatorial representation of the SU's local area feature and boundary morphology, acting as the coefficient of D a . The single excessive D a was neutralized with smaller coefficients when the area and common boundary of the SUs were not in accord with USSUs, thereby reducing the probability of being misjudged.
Finally, after performing a normalization procedure for he L and ho L , a general algebraic formula was employed to characterize the quantitative comparison of two parts [27]. The metric of local average similarity (LS) for each SU was calculated with Equation (13), where the LS was in the range [−1, 1].
It is easily inferred that a higher LS represented the stronger nature of ho L (local homogeneity) and was preferred to be identified as an OSSU. Conversely, a smaller or even close to −1 LS, represented the SUs with weak ho L but high he L (local heterogeneity), which is in accord with an USSU.
The metric of local similarity was able to transform the evaluation object from a general entirety to each precise individual, which was an obvious improvement in evaluation scales. For demonstrating the role of spatial boundary features and neighborhood variability in the identification of OSSUs and USSUs, three local similarities (LS) from different variables were prepared: (1) No A condition, where the factors related to difference D a or standard deviation D sd were replaced with the actual value; (2) No B condition, where the factor related to spatial size and boundary as r a and r l were discarded; (3) A and B condition, where both of the spatial features and neighborhood variability were taken into account.
The histogram frequency diagram of LS was established to demonstrate the distribution of the SUs with different properties at the optimal scale. Through the dual inspection of OCE and GS value, the selected optimal scale can make the majority of SUs consistent with the subdivision of the terrain and they were defined as moderate-subdivided SUs (MSSUs). The remaining small part of specific local heterogeneity or homogeneity was concentrated near the 1 or −1 of LS. The X-axis is the LS value, and 40 intervals of 0.05 were assigned in the range [−1, 1]. The discontinuity points of curve shapes on [−1, 0] and [0, 1] were used as threshold intervals to distinguish the USSUs, MSSUs and OSSUs. The discontinuities were calculated by the maximum between-cluster variance (Otsu) method [43].

The Optimization of USSUs and OSSUs
Setting separately variable (smaller or larger) scales for certain parts during the delineation of SUs remained a conceptual problem with operational difficulty [44]. After being marked as unrealistic SUs, a new effective optimization mechanism containing the re-subdivision of USSUs and merging of OSSUs was put into effect.
Firstly, the regions related to USSUs were substituted by corresponding SUs in multiple finer subdivisions, following which the substitution of the highest GS was considered as the best. Nevertheless, the optimal substitution could not guarantee to be the proper terrain resubdivision for all coexisting USSUs. Some substitutions might generate new subdivisions. Therefore, the re-subdivision of SUs will be marked as OSSUs for examination in the next optimization stage.
Then, although those aggregated and adjacent SUs could be amended conveniently and effectively through merging, the technical difficulties were merging sequences and preventing excessive merging. Referring to Section 3.3, the sub-region merging metrics of local similarity change (LSC) were defined. The principles of LSC merger guidelines were as follows: on the one hand, local heterogeneity (he L ) had a fundamental effect on determining the merging order. Among two or more neighboring OSSUs, the one of lower local heterogeneity had merging priority. On the other hand, the internal homogeneity changes between the new SUs and the original SUs needed to be estimated. The merging between quite different SUs could generate new SUs with complex internal structures. It is the situation of over-merging and should be avoided. The formula mode is shown in Equation (14), and each OSSU can obtain the value by merging with its adjacent and largest SUs.
where SU C is the new SU generated by the SU A merges with its adjacent SU B . y is the average aspect of SU and acts as the weight of the standard deviation, which are combined to represent the uniformity of SU. A given threshold of LSC T is used for quantifying those excessive merging situations. The LSC T is consulted from a large number of LSC produced by an artificial training merging process and is presented detailed in optimization results. In our study, we used a Debian GNU/Linux10.0 system (https://www.debian.org/releases/ stable/ (accessed on 10 June 2022)) and GRASS GIS 7.6.0 software to invoke the vector commands of v.category, v.edit and v.dissolve (https://grass.osgeo.org/grass80/manuals/ vector/ (accessed on 10 June 2022)) and execute the multi-scale optimization of USSUs and OSSUs.

Global Optimal Subdivision Scale
Considering the diversity of the Yuqu River Basin, the t parameter attempted to set additional multiple values, from 50 × 10 4 m 2 to 150 × 10 4 m 2 with an interval of 25 × 10 4 m 2 . The other parameters were set as the recommended defaults [5,14,15,39] ( Table 2). There were 30 valid terrain subdivisions obtained from different combinations of the (t, c) parameters. Figure 5 shows 8 of the 30 results, from the finest to the coarsest SU partitioning, from the upper left corner of (50, 0.1) to the lower right corner of (150, 0.6). The level of detail of the SUs' subdivisions depends heavily on the changes in parameter combinations of (t, c), especially the parameter t. However, the parameter t was considered as having no explicit geomorphological meaning and it was suggested to have a large value (500 × 10 4 m 2 ) in previous studies [14,15,39]. Compared with the reference SUs, extremely large or extremely small t (such as 50 or 150 × 10 4 m 2 ) cannot produce a proper subdivision of the landscape.  The quantitative indicators of OCE and GS were thoroughly exploited to determine the optimal scale. Firstly, the errors between 952 reference SUs and the corresponding extracted SUs of different scales were calculated. As is shown in Figure 6a, the errors were at a consistently high level at the finer scale. With the gradual increase of t from 50 to 100 × 10 4 m 2 , the errors decreased accordingly until they stabilized. However, the errors increased again as t continuously increased to 150 × 10 4 m 2 . The evolution features of OCE The quantitative indicators of OCE and GS were thoroughly exploited to determine the optimal scale. Firstly, the errors between 952 reference SUs and the corresponding extracted SUs of different scales were calculated. As is shown in Figure 6a, the errors were at a consistently high level at the finer scale. With the gradual increase of t from 50 to 100 × 10 4 m 2 , the errors decreased accordingly until they stabilized. However, the errors increased again as t continuously increased to 150 × 10 4 m 2 . The evolution features of OCE are significant because they are consistent with the three stages of over-subdivision, optimal-subdivision and under-subdivision as scale increases scale. The evolution trend analysis indicates that the t of 100 × 10 4 m 2 represents an appropriate scale.

4, x FOR PEER REVIEW
14 of 24 groups increasing from 15 to 30 or even 100, the maximum or minimum of MI and V may show new changes and eventually an uncertain optimal scale is generated. It is originated from the absence of an effective screening mechanism and an inherent drawback of unsupervised evaluation. Therefore, the candidate groups can be confirmed by retaining only meaningful subdivisions when jointing the supervised evaluation of OCE , thus producing a stable and consistent optimal scale. This complementary evaluation strategy is a significant improvement over the existing method.

Identification Result of USSUs and OSSUs
In this section, the OSSUs, MSSUs and USSUs were distinguished at the selected optimal scale. The frequency distribution with three kinds of LS are shown in Figure 7. In the No A condition (Figure 7a), the frequency of LS appeared to be relatively concentrated in the range of [−1.0, 0], and the frequency of the cumulative curve did not show discontinuity characteristics, as expected. In Figure 7b,c, instead, the frequency followed Furthermore, there were 15 groups of SUs selected as meaningful subdivisions to perform the normalization step, when OCE < 0.35 was set as the screening criterion. They are SUs' groups with the parameter combinations of t (75, 100, 125) and c (0.2, 0.3, 0.4, 0.5, 0.6). The GS value of those candidate groups is shown in Figure 6b. The maximum GS is 1.88 with the parameter combination of (100, 0.3), indicating that both internal homogeneity and external heterogeneity are the greatest. This is selected global optimal SUs' subdivision in this paper. As a comparison, under the method of Alvioli et al. [14,15] and Jacobs et al. [5] without the inspection of the OCE, the GS was calculated in Figure 6c. The highest GS value of all 30 groups was 1.83 produced by the parameter combinations of (100, 0.1), which would be judged as optimal under the existing method. However, the SUs of (100, 0.1) are distinctly over-subdivided for the landscape ( Figure 5) and the OCE is accordingly as high as 0.38. The scale of (100, 0.3) from the new proposed method is more appropriate for the selection of the optimal scale.
To understand why different optimal scales are generated by the two methods, each MI and V in the calculation of GS were collected in Figure 6d. The MI and V indices of different scales SUs groups varied widely. Following the gradual increase of the scale parameters, the holistic trend of MI index went through reduction, stabilization and increase successively, and the V index increased monotonically. With the candidate SUs' groups increasing from 15 to 30 or even 100, the maximum or minimum of MI and V may show new changes and eventually an uncertain optimal scale is generated. It is originated from the absence of an effective screening mechanism and an inherent drawback of unsupervised evaluation. Therefore, the candidate groups can be confirmed by retaining only meaningful subdivisions when jointing the supervised evaluation of OCE, thus producing a stable and consistent optimal scale. This complementary evaluation strategy is a significant improvement over the existing method.

Identification Result of USSUs and OSSUs
In this section, the OSSUs, MSSUs and USSUs were distinguished at the selected optimal scale. The frequency distribution with three kinds of LS are shown in Figure 7. In the No A condition (Figure 7a), the frequency of LS appeared to be relatively concentrated in the range of [−1.0, 0], and the frequency of the cumulative curve did not show discontinuity characteristics, as expected. In Figure 7b,c, instead, the frequency followed a normal distribution, and emerged with good discontinuity at the ends. Further, compared to A and B condition, the frequency of No B condition was more concentrated in the middle interval of LS, and less distributed at high values, implying that more SUs may be identified as MSSUs and insufficient SUs are classified into OSSUs.
Consequently, the thresholds of discontinuity points were calculated from the frequency cumulative curve via the Otsu method [43]. The results were summarized in Table 3 and typical examples of classification were given in Figure 7d-f. For the No A condition, due to the absence of neighborhood variability, poor discriminative power was exhibited for all three types of SUs, both in terms of quantity and accuracy (Figure 7d). Subsequently, the recognition capability of USSUs was improved in No B condition, as only the SU with distinct differences from its surrounding SUs were identified as USSUs, and the proportion descended from 61.14% to 3.14% (Table 3). However, many SUs with narrow and long boundaries were still wrongly selected as MSSUs and the recognition ability for OSSUs was insufficient (Figure 7e). Once spatial weight and neighborhood variability both functioned during the identification process, more SUs was marked as OSSUs, with an increase from 8.31% to 14.6%. The USSUs are also more reasonable, with large SUs correctly distinguished and small SUs excluded (Figure 7f   Consequently, the thresholds of discontinuity points were calculated from the frequency cumulative curve via the Otsu method [43]. The results were summarized in Table  3 and typical examples of classification were given in Figure 7d-f. For the No A condition, due to the absence of neighborhood variability, poor discriminative power was exhibited for all three types of SUs, both in terms of quantity and accuracy (Figure 7d). Subsequently, the recognition capability of USSUs was improved in No B condition, as only the SU with distinct differences from its surrounding SUs were identified as USSUs, and the proportion descended from 61.14% to 3.14% (Table 3). However, many SUs with narrow and long boundaries were still wrongly selected as MSSUs and the recognition ability for OSSUs was insufficient (Figure 7e). Once spatial weight and neighborhood variability both functioned during the identification process, more SUs was marked as OSSUs, with an increase from 8.31% to 14.6%. The USSUs are also more reasonable, with large SUs correctly distinguished and small SUs excluded (Figure 7f Figure 8 shows the spatial distribution of each type of SU in the Yuqu River Basin. A total of 6514 SUs, 106 cases of USSUs with an average area of 3.47 × 10 6 m 2 , and 952 cases of OSSUs with an average area of 2.41 × 10 5 m 2 . As details are shown in local regions, the large SUs that divided the whole gentle valley into one SU are correctly classified as USSUs (Figure 8a), and the narrow polygons formed by the interference of broken steep terrain are labeled as OSSUs (Figure 8b,c). Therefore, the method proposed for identifying USSUs and OSSUs is practicable.

Multi-Scale Optimization of Undesirable SUs
According to the multi-scale optimization method proposed in Section 3.4, ascertaining the merging threshold of the OSSUs is crucial. Some of the samples are shown in Figure 9a Table 4). The frequency distribution of LSC in each iteration is shown in Figure 9c; the main interval gradually increases with upper limit. Therefore, the frequency cumulative curve is constructed by summarizing the LSC of five merges. A value below 95% of all LSC is selected as the final merging criteria, that is, LSC T < 0.50.  Figure 8 shows the spatial distribution of each type of SU in the Yuqu River Basin. A total of 6514 SUs, 106 cases of USSUs with an average area of 3.47 × 10 6 m 2 , and 952 cases of OSSUs with an average area of 2.41 × 10 5 m 2 . As details are shown in local regions, the large SUs that divided the whole gentle valley into one SU are correctly classified as US-SUs (Figure 8a), and the narrow polygons formed by the interference of broken steep terrain are labeled as OSSUs (Figure 8b,c). Therefore, the method proposed for identifying USSUs and OSSUs is practicable.

Multi-Scale Optimization of Undesirable SUs
According to the multi-scale optimization method proposed in Section 3.4, ascertaining the merging threshold of the OSSUs is crucial. Some of the samples are shown in Figure 9a,b. There are 392 MSSUs obtained at the selected optimal scale (100, 0.3) and 1009 SUs assumed as OSSUs in the corresponding region at the finer scale (75, 0.3). A total of 1637 reference LSC values gained pass through five simulation merging processes ( Table  4). The frequency distribution of LSC in each iteration is shown in Figure 9c; the main interval gradually increases with upper limit. Therefore, the frequency cumulative curve is constructed by summarizing the LSC of five merges. A value below 95% of all LSC is selected as the final merging criteria, that is, T LSC < 0.50.    Subsequently, the gradual optimization process is shown in Figure 10. The original SUs at the selected optimal scale of (100, 0.3) (Figure 10a) are categorized as OSSUs, MSSUs, USSUs ( Figure 10b). As the second operation, the USSUs were replaced with the corresponding three finer SUs of (100, 0.2), (75, 0.2) and (75, 0.3), respectively. The substitution with the largest GS value of (100, 0.2) was considered the appropriate re-subdividing SUs. As is shown in Figure 10c, those USSUs of super-large area were effectively partitioned into two to four more refined SUs. There are 756 finer SUs generated by the re-subdividing process of the 106 USSUs and re-classified as OSSUs (Table 5). In compliance with the merging criteria of LSC T < 0.5, 1708 SUs were merged for the first time and 867 SUs required a next merge iteration. Finally, the same merging iterations were executed four times until the new LSC did not meet the condition. 1058 undesirable SUs were refined into 1117 appropriate SUs with the scale combinations of (100, 0.2), (100, 0.3) and other coarser scales. The OCE was improved to as low as 0.22 with a decrease of 27%.    With the aid of the field photography, a detailed comparison between the SUs of a single optimal scale and refined SUs after multi-scale recombination was carried out in the red box marked in Figure 10. As is shown in Figure 11c, almost each SU has achieved accurate matching with the terrain discontinuous line. Those tiny OSSUs without any terrain meaning have been merged with adjacent larger SUs. The USSUs of watersheds or large valleys have been re-subdivided along the inherent aspect turnings. Through comparing the SUs of yellow dashed lines, the refined SUs have better performance than the original SUs ( Figure 11b) in terms of consistency with the actual scene. Therefore, the proposed multi-scale optimization method has well improved the limitation of matchless between the SUs and the geomorphologic background under a single scale.

Rationality of the Evaluation and Post-Processing for the Optimal Scale
The existing measures cannot address issues in determining the input scale parameter for the automatic delineation of SUs. Additionally, there is a lack of subsequent solution for unrealistic SUs generated from a single optimal subdivision, especially in complex terrains. This paper applied an effective subdivision quality evaluation and multi-scale refinement approach to improve consistency in a complex natural landscape. The parameters and criteria involved in the improved method are objective and can be established independently from the geographical extent. It is the first application for the improvement of individual SUs objects from automatic subdivision procedures in the Tibetan Plateau.
The advantages of unsupervised evaluation cooperated with supervised means for determining the optimal scale were significant. The method enriched the scientific justification for matching between the global optimal scale and the actual land surface discontinuities. The uncertainty in calculating the GS value was removed through defining qualified candidate SUs' groups. Compared to previous research [14,15,39], the OCE

Rationality of the Evaluation and Post-Processing for the Optimal Scale
The existing measures cannot address issues in determining the input scale parameter for the automatic delineation of SUs. Additionally, there is a lack of subsequent solution for unrealistic SUs generated from a single optimal subdivision, especially in complex terrains. This paper applied an effective subdivision quality evaluation and multi-scale refinement approach to improve consistency in a complex natural landscape. The parameters and criteria involved in the improved method are objective and can be established independently from the geographical extent. It is the first application for the improvement of individual SUs objects from automatic subdivision procedures in the Tibetan Plateau.
The advantages of unsupervised evaluation cooperated with supervised means for determining the optimal scale were significant. The method enriched the scientific justification for matching between the global optimal scale and the actual land surface discontinuities. The uncertainty in calculating the GS value was removed through defining qualified candidate SUs' groups. Compared to previous research [14,15,39], the OCE could serve as a new supervised evaluation tool for comparing terrain subdivision algorithms and adjusting the subdivision parameters. The 952 reference samples used for this study are compliant with the standard SU definition, sampling from the entire basin at equal intervals. They are capable of representing the local topography and actual land surface discontinuities. The OCE can be used to reliably quantify the error with the subdivision SUs with different levels of details. The discrepancy value during the process of the overgrowth parameters can be sensitively confirmed, as the OCE accordingly went through the decreasing, plateauing and increasing phases. Moreover, the cooperative evaluation application was sequential. The OCE, acting as supervised evaluation, was first performed to filter out unrealistic terrain subdivisions, and qualified SUs groups were left to participate in the calculation of GS values. Consequently, the very fine scale (100, 0.1) improperly determined as the optimal subdivision by the previous method was eliminated. In addition, the parameter t contributed the fundamental upslope area for the calculation of the orientation average value and circular variance and should be dealt with as prudently as parameter c.
Further, the spatial characteristics of each SU's area, boundary and aspect were sufficiently emphasized and effectively integrated during the processes of identification and optimization. The measuring of local homogeneity and heterogeneity for each SU was ingenious and practicable. In particular, the variability of aspect beyond the normal level of the surrounding regions could more appropriately indicate the differences and similarities between each SU and its adjacent SUs. Due to the scarcity of available statistical properties, the diversity of the SUs size and boundary were taken into account and acted as a neutralization coefficient. Only when the difference or similarity of aspect were in accord with the features of boundaries and regions can the specificities of local heterogeneity or homogeneity be outstanding. Thus, this linkage mechanism could reduce the probability of misclassification. Although the assumption that the discontinuity of the curve was consistent with the threshold interval of the SUs' classification lacks sufficient theoretical basis, it was an efficient and reproducible objective method. Compared with the single optimum parameter optimal, those regions with special landforms demanded other coarser or finer parameters. Both the re-subdivided with the finer scale for USSUs and the re-merged with adjacent SUs for OSSUs were defined as detailed operation guidelines. An effective examination criterion was implemented for the two refinement stages, and the reference knowledge was obtained through many manual merging training processes. After a series of identification, re-subdivision and merging processes, the SUs' groups consisting of multiple scales were achieved. This is currently one of the few available methods for non-single scales.

Terrain-Adaptive Performance of the SUs Subdivision
Owing to the complex terrains of the Yuqu River Basin and the diversity of thousands of SUs, evaluating the performance of the SUs' subdivision against the corresponding topographic feature was not an easy target. The density, spatial scale, and spatial morphology of the SUs in four topographic patterns were statistically analyzed. In general, the slopes tended to be long and large at the bigger relative relief with severe geomorphic cutting, so the density, size and shapes of SUs should be distinguished from those of small and low elevation differences. Firstly, the I, II, III and IV geomorphic regions were divided into 2~4 intervals according to the 300 m interval of elevation difference. As shown in Figure 12a,b, the average density of SUs (the ratio of the number of SU to the area of the statistical interval) decreases linearly with the increase of the elevations' relative relief, and the average area of SUs increases. In the high hilly zone (I), the height difference was between 0.5~1.1 km, the SUs presented the most intensive distribution with average density values (1.47 SU/km 2 ), which corresponded to a large number of homogeneous low-lying hills and valleys. In the alpine gorge zone (IV), the height difference was between 1.4~2.3 km, the average area of SUs was shown as a maximum value of 0.99 km 2 , which is in accordance with a tall valley formed by the rapid river downcutting. landform, corresponding to a narrow strip shape change of SUs. Therefore, the multi-scale optimization of SUs was able to capture the morphological variability of the landscape and divided the study area into SUs with different shapes and sizes.

Conclusions
A set of appropriate SUs that properly partition a complex landscape into reasonable terrain subdivisions is crucial for landslide sensitivity modeling. Difficulties in determining the input scale parameters and the absence of subsequent amending procedures have restrained the wide application of existing automatic extraction methods. An improved automatic subdivision quality evaluation and multi-scale refinement method were proposed in this paper. An evaluation approach coordinated with supervised and unsupervised means were exploited. The OCE value is an effective discrepancy metric that can inhibit the uncertainty existing in optimal scale from the combination of global variance and global MI. Multiple space geometry features and aspects were effectively integrated into the process of the identification and refinements of undesirable SUs. In the example region, 6514 SUs were automatically obtained at the optimal subdivision scale, of which 106 cases were distinguished as USSUs and 956 were OSSUs. To better deal with unreasonable SUs at the selected single optimal scale, an effective optimization mechanism was established. The final SUs' groups composed of terrain subdivisions with multiple scales were achieved, which is an infrequently available method for non-single scales. Field investigations and statistical distribution characteristics have demonstrated the excellent performance of the SUs for the corresponding geomorphological reality. With improved organization and distribution of geoenvironmental data, the refining SUs have great application potential in landslide sensitivity modeling and other situations that require the identification of homogeneous terrain domains. In addition, the shapes of SUs were also sensitive to various topographic features. The shape index (R index) was calculated using the equation from Hang et al. [18,45], which was the ratio of the square of the perimeter to the corresponding area. For more narrow-flat or more strip-shaped polygons, the R index increased, such as the triangle is 20.78 and the rectangle that the length to width of 6:1 is 32.67. As shown in Figure 12c, the shape index of SUs from each geomorphic region were distinguishable. The average minimum value of the R index in the high hilly zone (I) was about 27.72, showing that the geometric shape more resembled a triangle. Even in the same height interval of 0.8~1.1 km, the R index of the middle-high mountain zone (II) was as high as 31.5, indicating that the shape of SUs were narrower and longer. It conformed to the characteristics of mountain landscape shaped by glacier erosion. The sensitive change of the R index also existed in the interior of geomorphic regions, especially the gorge transition zone (III). The R index in the interval [0.9-1.2 km] connecting the broad and gentle landform in the upper reach was 27.01, implying a more uniform shape of aspect ratio of the SUs. However, the R index was 34.05, up 25% in the interval of [1.2~1.5 km] transition to deep canyon landform, corresponding to a narrow strip shape change of SUs. Therefore, the multi-scale optimization of SUs was able to capture the morphological variability of the landscape and divided the study area into SUs with different shapes and sizes.

Conclusions
A set of appropriate SUs that properly partition a complex landscape into reasonable terrain subdivisions is crucial for landslide sensitivity modeling. Difficulties in determining the input scale parameters and the absence of subsequent amending procedures have restrained the wide application of existing automatic extraction methods. An improved au-tomatic subdivision quality evaluation and multi-scale refinement method were proposed in this paper. An evaluation approach coordinated with supervised and unsupervised means were exploited. The OCE value is an effective discrepancy metric that can inhibit the uncertainty existing in optimal scale from the combination of global variance and global MI. Multiple space geometry features and aspects were effectively integrated into the process of the identification and refinements of undesirable SUs. In the example region, 6514 SUs were automatically obtained at the optimal subdivision scale, of which 106 cases were distinguished as USSUs and 956 were OSSUs. To better deal with unreasonable SUs at the selected single optimal scale, an effective optimization mechanism was established. The final SUs' groups composed of terrain subdivisions with multiple scales were achieved, which is an infrequently available method for non-single scales. Field investigations and statistical distribution characteristics have demonstrated the excellent performance of the SUs for the corresponding geomorphological reality. With improved organization and distribution of geoenvironmental data, the refining SUs have great application potential in landslide sensitivity modeling and other situations that require the identification of homogeneous terrain domains.

Data Availability Statement:
The GNU/Linux system was obtained from (https://www.debian. org/releases/stable/, accessed on 10 April 2022) and GRASS GIS 7.8.1 was provided by project of the Open Source Geospatial Foundation (OSGeo) (https://grass.osgeo.org/download/windows/, accessed on 10 April 2022). The Software of r.slopeunits v1.0 was obtained from the Portal of the Geomorphology Research Group (http://geomorphology.irpi.cnr.it/tools/slope-units, accessed on 10 April 2022). The Gaofen-2 satellite data and ALOS-12m DEM data are provided by a third party, not a public dataset. The source data used in the analysis of this paper have been uploaded on as Supplementary Table S1.

Conflicts of Interest:
The authors declare no conflict of interest.