*2.5. Canopy Height Modeling*

The nonvegetation area in the land cover classification realized in the above steps was extracted to extract the nonvegetation points from the whole point cloud data. Nonvegetation points were merged with points of GCP on the field survey as modeling samples. As the geostatistical method, the empirical Bayesian kriging (EBK) was adopted after cross-validation with other four common models (S3) to construct a geostatistical model (i.e., digital terrain model (DTM)) for predicting the entire terrain elevation. Before modeling with these nonvegetation points, some points within other classes that could have been misclassified into nonvegetation points needed to be filtered using a slope-based filter method [46]. Furthermore, before modeling by EBK, declustering was performed and the global first-order trend was removed from the data (S3). The point clouds were also used to derive the digital surface model (DSM), a raster-based description of the surface that includes objects on the terrain, such as trees and shrubs. In this work, multilevel B-spline interpolation was used to generate the DSM, which is a high-speed interpolation algorithm for calculating the continuous surface from points based on irregular regional samples [47]. Finally, a superimposed subtraction of the two sets of data produced a canopy height model (CHM) by the grid difference, which is a raster computing tool used to subtract the DTM from the DSM to produce the CHM. The geostatistical analysis was accomplished by the ArcGIS Pro software. The processing of DSM and CHM data was achieved by the SAGAGIS software.

### *2.6. Extraction and Classification of Forests–Ecotone–Abandoned Landscape*

In this study, the ecotones between forest and abandoned agricultural land refer to mixed vegetation above the grass layer at abandoned land but below the forest layer, and they formed height gradient transitional zones. This vegetation height difference of transitions between abandoned agricultural land and forest ecosystems provided a meaningful result in this case and can be used as the characteristic for a three-dimensional landscape [39,48].

For extracting the ecotones between forest and abandoned agricultural land, two prerequisites are used: vegetation height and proportion of ecotones. In this study, the ecotone

detection method developed by Hou and Walz [48] was adopted. First, the land cover map obtained using OBIC described above is used as the thematic layer. The forest and agricultural land classes were overlaid with the CHM (DSM-DTM). The height difference is used as a criterion of ecotones. Thus, it was necessary to define the threshold values of height for the distinguish of ecotones. At the Baijixun plot, the average heights of the vegetation in the forest and abandoned agricultural land were computed by the zonal statistic. Second, the proportion of nonecotone pixels (i.e., forest pixels or abandoned land pixels) within a fixed size "window" was assumed to be between 40% and 60% based on the model from Riitters et al. [49]. Unlike the Hou and Walz method, CHM extracted from the photographic point cloud could be completely matched with the land cover map based on the high-resolution classification in our study. Therefore, combined with Hou and Walz's method and the perspective of Riitters, the classification results obtained from the above classification system needed to be shrunken only. In our case, that meant inside a moving window (61 by 61 pixels) surrounding the ecotone pixel, the proportion of both forest and abandoned land pixels should be less than 60%. Combined with the result of classification and CHM data, the average height of forest (Hf) and the average height of abandoned land (Ha) in the landscape were calculated. The detailed step was to form transitions by shrinking grass along the border pixels where CHM ≥ Ha m and the relative area of grass pixels in a moving window of (61 by 61) pixels < 60% and shrinkage of the forest where CHM ≤ Hf m and the relative area of forest pixels in a moving window of (61 by 61) pixels < 60% (Figure 3). Finally, we assigned ecotone class to the unclassified class and involved eliminating impurities inside ecotones. It should be emphasized that considering for the shrinking from the nongrass land class (i.e., shadow class and bare land class), the shadow class and bare land class were reassigned into the interior of forest class or the interior of grass class according to the proportion of classes in the nearest pixel before do the processing of shrinking. The shrunken work was completed with a tool that is 'pixel-based object resizing' in eCognition Developer software.

**Figure 3.** Shrinking process diagram. The moving window along the border between the abandoned land and the forest was used as the judgment basis for the shrinking into the interior of the forest and the interior of the abandoned land where *Pa* denotes the proportion of pixels in the whole window of the abandoned land. Moreover, the *Pf* denotes the proportion of pixels in the whole window of the forest. Ha and Hf denote the average height of abandoned land and forest, respectively.

### *2.7. Quantifying the Landscape Containing Ecotones*

Compared with the traditional ecotone identification method, this is an obvious advantage of the spatial continuous detection results, that is, the spatial quantization of the ecotone can be realized. Mapping the landscape pattern, including the ecotones, not only shows where the ecotones occur but also makes it possible to quantify the scale of the ecotones. At the CLASS level, a series of landscape metrics are computed in an R package (landscapemetrics [50]) to assess the landscape characteristics of ecotones. The quantitative indicators selected are shown in Table S3. Furthermore, the total area (TA) reflects the size of the corresponding class in the landscape and measures the composition of the landscape. The percentage of landscape (PLAND) is a relative measure that quantifies the landscape weight of each class in the landscape, which is also an important indicator for measuring the landscape composition. The total edge (TE) reflects the total edge length of a class, and longer edge lengths may be associated with larger edge effects. The shape index (SI) is the most straightforward measure of complexity, and the changes in a landscape may reflect the extent of disturbance. In fact, these common indicators are also rich in ecological significance. TA is a direct measure of the scale of landscape elements, and the area of vegetation has an impact on the species in it, that is, the area effect [51]. While PLAND reflects the weight of corresponding classes in the landscape, it also indirectly reflects the contribution of this class to the spatial heterogeneity of landscape elements within the landscape. It is assumed that the larger the PLAND of a class is, the more outstanding the contribution of the corresponding class to the homogeneity of the habitat types within the landscape is, which may influence the biodiversity within the landscape, because it may lead to concentrations of rare species in marginal habitats [52]. TE and SI measure the edge characteristics of land cover types, which may be related to the edge effect of the corresponding land cover class [53].

### *2.8. Transect-Based Analysis*

The ecotones occur between adjacent ecosystems frequently with high heterogeneity. In general, the variability of biodiversity is sensitive to environmental change. The species importance value (IV) is computed in every species per small sample quadrats. In addition, the mean values of the importance values of plants included in each quadrat are calculated to represent the importance values at the level of the quadrat. This work (i.e., the measurement of IV, pH, and OM) has been conducted in previous experiments (Table S1). For analysis of the variation along transects from abandoned land to the interior of the forest, the MSW alone transect was adopted in this study. The MSW is a common approach to reflect the habitat change or indicate how habitats are separated [54]. The discrepancy coefficient was computed between two subwindows (1/2 window) using the squared Euclidean distance (SED) (Formula (1))

$$SED\_n = \left(\overline{X}\_{inv} - \overline{X}\_{ibw}\right)^2\tag{1}$$

where *n* denotes the location of the medium point of two subwindows and *a* and *b* respectively represent the two subwindows when the window is *n*. The *w* is the size of a window. The *m* is the number of parameters in a quadrat [55]. The times of removing are determined by subtracting the size of the window from the number of total quadrats and adding 1. In addition, the experiment was carried out several times with a variety of window sizes so that any potential bias caused by the impact of window scale may be eliminated. According to the findings of our earlier study, the position and breadth of the ecotones may be identified more precisely when the size of a window is set to 6. [56]. In this study, eight medium points were determined (The distance along the one-dimensional transect from the starting point is 17.5 m, 25 m, 32.5 m, 40 m, 47.5 m, 55 m, 65 m, and 75 m, respectively) by the window size and the number of all quadrats.

Meanwhile, the obtained results were normalized. The calculated distance coefficient (i.e., SED) was plotted along the transect. The dependent variable was the value of the normalized result, and the independent variable was the position of the starting distance from the transect. Similarly, the width of ecotones occurrence was expressed at the same location of the transects. An interpolation line was set up through the midpoint of each small quadrat, the starting point and ending point of the entire transect. The landscape map containing ecotones was binarized before computing by interpolation line, with one as the ecotones and zero as the nonecotone. Based on the graph, the occurrence of the location of ecotones was detected according to crest and estimate the scale of ecotone according to the width of wave.

### *2.9. Accuracy Assessment for Ecotone Detection*

Two aspects of accuracy need to be verified in this study. One aspect was the classification accuracy of the landcover of biotopes, which will directly affect the ecotones' location of occurrence. Another aspect was the delineation accuracy of the ecotones, which tests the accuracy of the ecotones' width based on our method, reflecting the scale (width) of the transition.

In this study, field surveys are carried out to establish the characteristics and variability of land cover and to acquire reliable field data for training the classifiers of the UAS image data and evaluating the accuracy of the classification results. The ground truth data were combined with the UAS images to assess the accuracy of the classification. The F1-measure was adapted for assessing the accuracy of each class [48]. This indicates the harmonic mean between precision (*pi*) and recall (*ri*) for each class *i*, because recall and precision are evenly weighted [57].

$$F\_1 = \frac{2p\_ir\_i}{p\_i + r\_i} \tag{2}$$

*pi* denotes objects correctly classified as class *i*/(ground reference objects in class *i*) and *ri* denotes objects correctly classified as class *i*/(number of objects correctly classified as class *i* + objects falsely classified as class *i*).

Meanwhile, the F1-measure was also used to test the accuracy of the scale of ecotones occurrence extracted based on proposed method.

$$F\_1 = \frac{2p\_j r\_j}{p\_j + r\_j} \tag{3}$$

$$p\_{\dot{j}} = \frac{TP\_{\dot{j}}}{TP\_{\dot{j}} + FP\_{\dot{j}}} \tag{4}$$

$$r\_{\dot{j}} = \frac{TP\_{\dot{j}}}{TP\_{\dot{j}} + FN\_{\dot{j}}} \tag{5}$$

*TPj* denotes the overlap width of the ecotones derived from the UAS data and the ecotones detected from the analysis of the transects based on situ investigation at transect *j*. *FPj* denotes the width of the part not detected by situ data but detected in the UAS data. *FNj* denotes the width of the part not detected by UAS data but detected in the situ data. The *pj* represents Precision, and the *rj* represents Recall.

### **3. Results**

### *3.1. Canopy Height Estimate*

In the Baijixun sample site, the unbiasedness and uncertainty of multiple geostatistical models were compared by cross-validation method, and the EBK was finally selected to model the DTM. A CHM was successfully used for modeling, and it was used to interpolate the height value across the study area (Figure 4). From the results, the predicted slope direction was consistent with the actual slope direction, and the vegetation height changed obviously at the boundary between abandoned land and forest (in the north) (Figures 1 and 4). According to the results of OBIC, we calculated that the average height of the forest class was 3.84 m, and that of the abandoned land class was 0.27 m.

**Figure 4.** The digital surface model (DSM) (**a**), digital terrain model (DTM) (**b**), and canopy height model (CHM) (**c**) of the Baijixun plot.
