Disturbance unit division method in soil loss amount measurement and calculation process

Abstract

The invention relates to the technical field of water and soil conservation, and discloses a disturbance unit division method in a soil loss amount measuring and calculating process, which comprises the following steps of: acquiring a multi-period remote sensing image and a multispectral image of a project area, and performing spatial registration preprocessing on the multispectral image; multi-dimensional disturbance feature extraction is carried out; performing multi-scale image segmentation, and intelligently identifying a disturbance unit in combination with a segmentation result; eliminating classification noise, and extracting a disturbance unit boundary; smoothing the boundary of the disturbance unit; judging the disturbance type of each disturbance unit; carrying out area calculation and attribute labeling based on a judgment result; performing disturbance unit division on the multi-period remote sensing image, and performing change detection according to a division result; performing space-time analysis on a detection result; standardized result output and precision verification are carried out; through multi-source remote sensing data fusion and an intelligent classification algorithm, technicians are liberated from heavy manual sketching work, and the division efficiency is improved compared with that of a traditional manual method.

Description

Technical Field

[0001] This invention relates to the field of soil and water conservation technology, and more specifically, to a method for dividing disturbance units in the process of soil loss measurement. Background Technology

[0002] Construction projects inevitably disturb the existing land surface during construction, creating different types of disturbed surfaces, including general disturbed surfaces, excavation faces, and engineering accumulations. These disturbed surfaces are major sources of soil erosion, and accurately measuring their soil loss is crucial for soil and water conservation planning, monitoring and assessment, and ecological restoration decisions. According to the guidelines for soil loss measurement in construction projects, the first step in soil loss measurement is to delineate disturbed units, which involves scientifically classifying and delineating the boundaries of disturbed surfaces within the project area caused by construction activities based on characteristics such as disturbance mode, disturbance intensity, soil type, meteorological conditions, and spatial continuity.

[0003] Currently, in practice, the delineation of disturbance units mainly relies on the manual visual interpretation and experience-based judgment of technical personnel. Through on-site surveys, combined with basic data such as topographic maps and engineering design drawings, technical personnel manually delineate the boundaries of different disturbance types based on their personal experience. This traditional method suffers from problems such as low delineation efficiency, difficulty in guaranteeing delineation accuracy and consistency, insufficient ability to identify complex disturbance characteristics, and poor traceability and verifiability, making it difficult to meet the needs of modern precise control and information-based management of soil and water conservation.

[0004] With the rapid development of remote sensing and artificial intelligence technologies, it has become possible to automatically divide disturbed units using multi-source remote sensing data and intelligent classification algorithms. However, existing technologies mainly focus on remote sensing identification and classification of disturbed areas, failing to strictly follow guidelines to classify disturbed units into types such as general disturbed surfaces, engineering excavation faces, and engineering accumulation bodies. They also fail to combine digital elevation model differential data for quantitative determination of disturbance intensity, and do not achieve refined optimization of disturbed unit boundaries and automatic labeling of scale levels, making them difficult to directly apply to soil loss measurement processes. Summary of the Invention

[0005] This invention provides a method for dividing disturbance units in the soil loss measurement process, which solves the technical problems of low efficiency, unstable accuracy, strong subjectivity, and poor traceability in the division of disturbance units in related technologies.

[0006] This invention provides a method for dividing disturbance units in the soil loss measurement process, including:

[0007] Multiple remote sensing images and multispectral images of the project area were collected. Spatial registration preprocessing was performed on the multispectral images to obtain radiometrically normalized fused image data.

[0008] Based on radiometrically normalized fused image data, multidimensional perturbation features are extracted to obtain a standardized multidimensional feature dataset.

[0009] Multi-scale image segmentation is performed based on a standardized multi-dimensional feature dataset, and perturbation units are intelligently identified by combining the segmentation results to obtain preliminary perturbation unit classification results data.

[0010] Based on the preliminary perturbation unit classification results, classification noise is eliminated and perturbation unit boundaries are extracted; the perturbation unit boundaries are then smoothed to obtain perturbation unit vector boundary data.

[0011] Based on the vector boundary data of the perturbation unit and the preliminary classification results of the perturbation unit, the perturbation type of each perturbation unit is determined; based on the determination results, the area is calculated and the attributes are labeled to obtain perturbation unit data with attributes.

[0012] Disturbance units are divided into multiple remote sensing images, and change detection is performed based on the division results; the detection results are then subjected to spatiotemporal analysis to obtain spatiotemporal change data of the disturbance units.

[0013] Based on attributed perturbation unit data and spatiotemporal variation data of perturbation units, standardized results are output and accuracy is verified to obtain standardized perturbation unit result data.

[0014] In a preferred embodiment, the spatial registration preprocessing of the multispectral image includes:

[0015] Radiometric calibration is performed using sensor calibration parameters to convert the digital quantization values ​​of multispectral images into physical radiance values. A fast atmospheric correction algorithm is used to eliminate the influence of atmospheric scattering and absorption on surface reflectivity.

[0016] Subpixel-level spatial registration is performed using an automatic registration algorithm based on feature point matching. Scale-invariant feature transformation feature points are extracted from the two multispectral images, and the nearest neighbor distance ratio method is used for feature point matching. False matching points are removed using a random sampling consensus algorithm.

[0017] A principal component transformation fusion algorithm is used to perform principal component transformation on multispectral images, extract the first principal component, perform histogram matching between high-resolution panchromatic images and the first principal component, replace the first principal component with the matched panchromatic images, and perform inverse principal component transformation.

[0018] The terrain factors were extracted using a third-order inverse distance weighted difference algorithm, and the radiometric consistency of the multi-temporal images was processed using a relative radiometric normalization method to obtain radiometrically normalized fused image data.

[0019] In a preferred embodiment, the multidimensional perturbation feature extraction includes:

[0020] Based on radiometrically normalized fused image data, normalized vegetation index and soil-adjusted vegetation index were calculated to obtain spectral feature layer data.

[0021] Based on radiometrically normalized fused image data, texture features are extracted at multiple scales using the gray-level co-occurrence matrix method to obtain texture feature layer data.

[0022] Based on the digital elevation model data before and after construction, the digital elevation model is preprocessed with median filtering to calculate the terrain change and obtain terrain change feature layer data.

[0023] The Z-score normalization method was used to normalize the spectral feature layer data, texture feature layer data, and terrain change feature layer data. The normalization results were then fused, and principal component analysis was used to reduce the dimensionality of the multidimensional features to obtain a standardized multidimensional feature dataset.

[0024] In a preferred embodiment, the multi-scale image segmentation and perturbation identification unit based on the segmentation results includes:

[0025] Based on a standardized multidimensional feature dataset, a multi-scale segmentation algorithm is used to generate image objects, resulting in image segmentation object data.

[0026] Based on image segmentation object data, a classification training sample set was constructed using a sample selection method. The training sample set includes samples from four categories: general disturbed surface, engineering excavation face, engineering accumulation body, and undisturbed area. An improved synthetic minority class oversampling technique was used to augment the samples, resulting in an augmented training sample dataset.

[0027] The expanded training sample dataset is trained using a random forest classifier to obtain a trained random forest classification model.

[0028] Based on the trained random forest classification model, a multi-classifier ensemble system is constructed by combining a support vector machine classifier and a gradient boosting decision tree classifier. The ensemble decision adopts a weighted voting mechanism to obtain the ensemble classification prediction results.

[0029] Morphological post-processing is performed on the integrated classification prediction results. Opening operation is used to remove small noise patches, closing operation is used to fill small holes, and a classification confidence evaluation mechanism is introduced to obtain preliminary perturbation unit classification result data.

[0030] In a preferred embodiment, the elimination of classification noise includes:

[0031] Based on the preliminary perturbation unit classification results, the morphological opening operation method is used to eliminate isolated noise points in the classification results, and the morphological filtering results are obtained.

[0032] Based on the morphological filtering results, the Canny edge detection algorithm is used to extract the boundaries of perturbation cells. The algorithm includes four steps: Gaussian filtering, gradient calculation, non-maximum suppression, and double threshold detection. The edge connection method is used to generate continuous boundary lines of perturbation cells, thus obtaining boundary line data.

[0033] In a preferred embodiment, the smoothing process of the perturbation unit boundary includes:

[0034] Based on boundary line data, a vectorization method is used to convert raster boundaries into vector boundaries, the cell coordinates of the edge chain are converted into geographic coordinates, the Douglas-Peucker algorithm is used to simplify the boundary points, and the Bézier curve smoothing method is used to smooth the simplified boundaries to obtain smooth vector boundary data.

[0035] Spatial overlay analysis is performed based on smooth vector boundary data. Disturbance cells located outside the project boundary and those overlapping with known non-disturbance areas are removed. Topology checks are then performed, and automatic repair methods are used to process any topology errors found, resulting in the obtained vector boundary data of the disturbance cells.

[0036] In a preferred embodiment, determining the disturbance type of each disturbance unit includes:

[0037] Establish the spatial correspondence between the perturbation unit vector polygon and the raster classification results;

[0038] Count the number distribution of each category of pixels within each vector polygon, traverse all pixels within the polygon range, and read the category label of each pixel in the integrated classification prediction result;

[0039] The dominant disturbance type of the disturbance unit is determined by majority voting.

[0040] Extract the slope values ​​of all pixels within the range of each perturbation cell, and calculate the average slope of the perturbation cell using the arithmetic mean method;

[0041] The determined disturbance type and calculated average slope are assigned as attribute information to the corresponding vector polygon to establish a disturbance unit type attribute data table, thereby obtaining the disturbance unit type attribute data.

[0042] In a preferred embodiment, the area calculation and attribute labeling based on the determination result includes:

[0043] For disturbance units with a general disturbance type of surface, the vector polygon area algorithm is used, which calculates the polygon area based on the geographic coordinates of the polygon vertices using the shoelace algorithm.

[0044] For disturbance units whose disturbance type is engineering excavation face and engineering accumulation, based on the topographic change feature layer data, the absolute value of the topographic change of all pixels within the disturbance unit is calculated, multiplied by the pixel area, and summed to obtain the total volume of the disturbance unit;

[0045] The disturbance cells are classified into three levels: large, medium, and small, using a scale-level determination method. An attribute assignment method is then used to generate a vector map of the disturbance cells with attributes.

[0046] In a preferred embodiment, the process involves dividing multi-period remote sensing images into disturbance units, detecting changes based on the division results, and performing spatiotemporal analysis on the detection results.

[0047] Based on multiple periods of remote sensing image data, the processing flow of multi-source remote sensing data acquisition and spatial registration preprocessing, multi-dimensional perturbation feature extraction, object-oriented image segmentation and supervised learning classification, boundary optimization and vectorization processing, area and volume calculation and attribute annotation was repeated to obtain the perturbation unit division results for each period.

[0048] Spatial overlay analysis was used to detect changes by overlaying vector maps of disturbance units from the early and later periods to generate an overlay result layer. The changes in disturbance units were identified by comparing the attribute values ​​of the two periods.

[0049] The change type determination method is used to identify the change type of the disturbance element, and the change area and change volume calculation method is used to quantify the degree of change of the disturbance element;

[0050] The evolution of perturbation units is tracked using time-series trajectory analysis. A unique tracking identifier is assigned to each perturbation unit. Based on the spatial overlap and attribute similarity of perturbation units in adjacent periods, the correspondence between perturbation units in the early and late periods is established.

[0051] Spatiotemporal evolution attributes of disturbance units are generated using a spatiotemporal attribute generation method.

[0052] In a preferred embodiment, the standardization output and accuracy verification include:

[0053] The data is output in a standard format using data format conversion methods.

[0054] The spatial integrity of the disturbance element vector data was verified using a topology check method.

[0055] Based on field survey data of the project area, an accuracy verification method was used to evaluate the classification accuracy of the disturbance unit division results;

[0056] The quality of the perturbation unit partitioning result is judged by the accuracy evaluation method. When the overall accuracy of the perturbation unit partitioning result is higher than the preset accuracy threshold and the Kappa coefficient is higher than the preset Kappa threshold, the partitioning result is considered to meet the application requirements.

[0057] Metadata files for the results data are created using metadata generation methods.

[0058] The beneficial effects of this invention are as follows:

[0059] This invention automates the entire process of disturbance unit division. By fusing multi-source remote sensing data and using intelligent classification algorithms, it frees technicians from tedious manual drawing work. The division efficiency is improved compared to traditional manual methods. It is particularly suitable for soil and water conservation monitoring and management of large-scale, linear production and construction projects, and effectively responds to the timeliness requirements of modern engineering management.

[0060] This invention constructs a multi-dimensional disturbance feature extraction system, which comprehensively utilizes multi-dimensional information such as spectrum, texture, and topographic changes to identify and classify disturbance units. It overcomes the limitations of traditional methods that rely solely on a single feature or human experience, and can objectively and accurately identify indicators that require quantitative judgment, such as disturbance intensity. This ensures the accuracy, scientific nature, and consistency of the classification results among different projects, and provides standardized and normalized data support for subsequent sampling of typical units and calculation of soil loss. Attached Figure Description

[0061] Figure 1 This is the main flowchart of the present invention;

[0062] Figure 2 This is a detailed flowchart of the present invention. Detailed Implementation

[0063] The subject matter described herein will now be discussed with reference to exemplary embodiments. It should be understood that these embodiments are discussed only to enable those skilled in the art to better understand and implement the subject matter described herein, and changes may be made to the function and arrangement of the elements discussed without departing from the scope of this specification. Various processes or components may be omitted, substituted, or added as needed in the examples. Furthermore, some features described in the examples may be combined in other examples.

[0064] At least one embodiment of the present invention discloses a method for dividing disturbance units in the soil loss measurement process, such as... Figures 1 to 2 As shown, it includes the following steps:

[0065] Step 1: Collect multiple remote sensing images and multispectral images of the project area, perform spatial registration preprocessing on the multispectral images, and obtain radiometrically normalized fused image data.

[0066] Step 1.1, Acquisition of multi-source remote sensing data;

[0067] Based on the project area coordinates and construction timeline, multiple phases of remote sensing and multispectral imagery covering the project area were acquired from a remote sensing data platform. These phases included imagery data from key periods such as before, during, and after construction. The spatial resolution of the multispectral imagery was set to be better than one meter, and it included visible light, near-infrared, and shortwave infrared bands. Simultaneously, digital elevation model data for both before and after construction, as well as a vector map of the project construction layout, were acquired. During data acquisition, it was ensured that the cloud cover was below 10%, and the imaging time was matched to the construction schedule.

[0068] Step 1.2, remote sensing image preprocessing;

[0069] Based on the multispectral imagery from step 1.1, radiometric calibration is performed using sensor calibration parameters to convert the digital quantization values ​​of the imagery into physical radiance values, eliminating the influence of sensor response differences. A fast atmospheric correction algorithm is employed to eliminate the impact of atmospheric scattering and absorption on surface reflectance. This algorithm, based on an atmospheric radiative transfer model, estimates atmospheric parameters using dark target pixels in the imagery to retrieve the true surface reflectance. Geometric correction is then performed based on a digital elevation model to eliminate geometric distortions caused by topographic relief, ensuring that the position of each pixel in the image accurately corresponds to its actual position on the ground, resulting in radiometrically corrected and geometrically corrected image data.

[0070] Step 1.3, spatial registration of multispectral images;

[0071] Based on the geometrically corrected image data from step 1.2, a sub-pixel-level spatial registration is performed using an automatic registration algorithm based on feature point matching. Scale-invariant feature transformation points are extracted from both images, and the nearest neighbor distance ratio method is used for feature point matching. A random sampling consensus algorithm is used to remove mismatched points, and the registration accuracy is set to be better than three pixels from zero. The registration accuracy setting rule is determined according to the accuracy requirements of subsequent analysis, generally taking one-third to one-half of the pixel size, to obtain the spatially registered image data.

[0072] Step 1.4, Multispectral image fusion processing;

[0073] Based on the spatially registered image data from step 1.3, a principal component transformation (PCT) fusion algorithm is used to perform PCT on the multispectral image to extract the first principal component. Histogram matching is then performed between the high-resolution panchromatic image and the first principal component. Finally, the matched panchromatic image replaces the first principal component, and an inverse PCT transformation is performed. The fused image simultaneously possesses spatial detail information from the high-resolution image and spectral information from the multispectral image, resulting in high-resolution fused image data.

[0074] Step 1.5, Topographic factor extraction and radiation normalization;

[0075] Based on the digital elevation model data from step 1.1 and the high-resolution fused image data from step 1.4, a third-order inverse distance weighted difference algorithm is used to extract topographic factors such as slope and aspect. This algorithm uses the elevation values ​​of eight neighboring pixels around the central pixel to calculate the elevation change rates in the east-west and north-south directions, and then calculates the slope value. The aspect is calculated based on the arctangent values ​​of the elevation change rates in the east-west and north-south directions. A relative radiometric normalization method is used to perform radiometric consistency processing on the multi-temporal images. Pseudo-invariant feature regions are selected on both the reference image and the image to be normalized. A linear transformation is used to adjust the radiometric values ​​of the image to be normalized to the same level as the reference image, ensuring the reliability of subsequent time-series analysis, resulting in radiometrically normalized fused image data and topographic factor raster data.

[0076] Step 2: Based on the radiometrically normalized fused image data, perform multidimensional perturbation feature extraction to obtain a standardized multidimensional feature dataset;

[0077] Step 2.1, extraction of spectral features of vegetation cover;

[0078] Based on the radiometrically normalized fused image data from step 1, the Normalized Difference Vegetation Index (NDVI) and the Soil-Adjusted Vegetation Index (SDI) are calculated. Specifically, the NDVI is calculated based on reflectance data from the near-infrared and red bands, effectively characterizing surface vegetation cover. The value ranges from -1 to +1; a higher value indicates better vegetation cover, while a lower value indicates poorer vegetation cover or non-vegetated features such as bare soil or buildings. The calculation process involves extracting the near-infrared and red band reflectance values ​​for each pixel in the image, calculating the difference between the near-infrared and red band reflectance, and then calculating the sum of the near-infrared and red band reflectance values. Dividing the difference by the sum yields the NDVI.

[0079] The Soil-Adjusted Vegetation Index (SVEI) effectively reduces the impact of soil background on vegetation indices by introducing a soil adjustment factor. The soil adjustment factor is determined based on the vegetation cover of the study area, with a setting rule of 0.5 for medium vegetation cover, 0.25 for high vegetation cover, and 0.75 for low vegetation cover. The calculation process is as follows: The soil adjustment factor value is determined based on the vegetation cover of the study area. Then, the near-infrared and red band reflectance values ​​of each pixel in the image are extracted. The difference between the near-infrared reflectance and the product of the soil adjustment factor and the red band reflectance is calculated. The sum of these products is then calculated, and the soil adjustment factor is added. Finally, the difference is divided by the sum to obtain the SVEI. This index has higher recognition accuracy in sparse vegetation areas and can more accurately identify slightly disturbed areas, obtaining spectral feature layer data.

[0080] Step 2.2, Multi-scale texture feature extraction;

[0081] Based on the radiometrically normalized fused image data from step 1, texture features were extracted at multiple scales using the gray-level co-occurrence matrix method. Distance parameters were set to one, two, and four pixels, and direction parameters were set to 0 degrees, 45 degrees, 90 degrees, and 135 degrees. Texture parameters such as contrast, homogeneity, entropy, and second moment of angle were calculated. The distance parameter was set according to the typical size of the target feature, generally taking multiple values ​​within the range of one to five pixels. Contrast reflects the roughness of the texture and is calculated by statistically analyzing the sum of squares of gray-level differences between pixel pairs; homogeneity reflects the uniformity of the texture and is calculated by statistically analyzing the distribution of similar gray-level pixel pairs; entropy reflects the complexity of the texture and is calculated by statistically analyzing the randomness of the gray-level distribution. Excavated surfaces are smooth with low contrast and high homogeneity; excavated surfaces are rough with high contrast and high entropy; generally, disturbed surface texture features fall between these two. The fusion of multi-scale texture features employed a feature vector concatenation method, combining texture features from different scales into a multi-dimensional feature vector to enhance the ability to distinguish different types of disturbance, resulting in texture feature layer data.

[0082] Step 2.3, Calculation of terrain change characteristics;

[0083] Based on the digital elevation model (DEM) data before and after construction obtained in step 1.1, median filtering preprocessing was performed on the DEM. The filtering window was set to 3x3 pixels to effectively eliminate the influence of data noise and outliers. The filtering window size was determined according to the spatial resolution and noise level of the DEM, generally ranging from 3x3 to 7x7 pixels. The topographic change was calculated, i.e., the elevation after construction minus the elevation before construction. Positive values ​​represent fill areas, negative values ​​represent cut areas, and the absolute value represents the thickness of the fill and cut areas. Topographic change is the core basis for judging the intensity of disturbance, enabling quantitative identification of the excavation face and the accumulated material. Simultaneously, the topographic change intensity index was calculated, i.e., the ratio of the absolute value of the topographic change to the standard deviation of the neighboring elevation. This index better reflects the intensity of disturbance relative to the surrounding environment, effectively distinguishing between real disturbances and natural terrain undulations, thus obtaining topographic change feature layer data.

[0084] Step 2.4, Feature normalization and fusion optimization;

[0085] Based on the spectral feature layer data from step 2.1, the texture feature layer data from step 2.2, and the terrain change feature layer data from step 2.3, the Z-score standardization method was used to normalize various features, eliminating dimensional differences between different features. The standardization process involves subtracting the mean from the original feature values ​​and then dividing by the standard deviation, ensuring all features have the same numerical range and distribution characteristics. The normalized spectral features, texture features, and terrain change features were then fused to construct a twelve-dimensional feature vector, including four dimensions for spectral features, six dimensions for texture features, and two dimensions for terrain features. Principal component analysis was used to reduce the dimensionality of the multidimensional features, retaining principal components with a cumulative contribution rate of 95%. The cumulative contribution rate threshold was set according to the feature information retention requirements, generally between 90% and 98%, thus maintaining the integrity of feature information while reducing computational complexity, resulting in a standardized multidimensional feature dataset.

[0086] Step 3: Perform multi-scale image segmentation based on standardized multi-dimensional feature dataset, and intelligently identify perturbation units by combining the segmentation results to obtain preliminary perturbation unit classification results data;

[0087] Step 3.1, multi-scale image segmentation;

[0088] Based on the standardized multidimensional feature dataset from step 2, a multi-scale segmentation algorithm is used to generate image objects. The multi-scale segmentation algorithm is a bottom-up region merging and segmentation method. Starting from a single pixel, it iteratively merges adjacent regions based on spectral and shape similarity until a preset heterogeneity threshold is reached. The segmentation scale parameter controls the size of the image objects, and the setting rule is determined based on the typical size of the perturbation unit, generally ranging from 50 to 200. The shape weight parameter controls the relative importance of spectral and shape information in the merging process, and the setting rule is determined based on the regularity of the ground feature boundaries, generally ranging from 0.1 to 0.3. The compactness weight parameter controls the compactness of the image objects, and the setting rule is determined based on the shape complexity of the ground features, generally ranging from 0.3 to 0.7, resulting in image segmentation object data.

[0089] Step 3.2, Training Sample Construction and Expansion;

[0090] Based on the image segmentation object data from step 3.1 and existing survey data of the project area, a classification training sample set was constructed using a sample selection method. The training sample set includes samples from four categories: generally disturbed surfaces, engineering excavation faces, engineering accumulation bodies, and undisturbed areas. The number of samples in each category was set according to the rule of meeting statistical significance requirements, generally no less than fifty samples per category, and the spatial distribution of the samples should cover different areas and different disturbance intensity levels within the project area. An improved synthetic minority class oversampling technique was used to augment the samples, effectively solving the problems of insufficient samples and class imbalance. For minority class samples, synthetic samples were generated within their five nearest neighbors. The new samples were generated by random interpolation between the original samples and their nearest neighbors. After quality assessment, the generated synthetic samples were added to the training sample set. The sample quality assessment was based on two indicators: feature similarity and class consistency. The quality assessment threshold was set according to the rule of ensuring the reliability of the augmented samples, generally between 0.8 and 0.9, resulting in the augmented training sample dataset.

[0091] Step 3.3, Random Forest Classifier Construction and Training;

[0092] Based on the expanded training sample dataset obtained in step 3.2, a random forest classifier is used for training. The random forest classifier is an ensemble learning method that constructs multiple decision trees and votes on their predictions. The main parameters of the random forest classifier include the number of decision trees and the maximum depth of each tree. The number of decision trees is determined based on the size of the training sample set and the feature dimension, typically ranging from one hundred to five hundred trees; the maximum depth is set to prevent overfitting, typically between ten and twenty layers. A grid search method is used to optimize the classifier parameters, and cross-validation is used to evaluate the classification performance of different parameter combinations. The parameter combination with the highest validation accuracy is selected as the final model parameters, resulting in the trained random forest classification model.

[0093] Step 3.4, multi-classifier ensemble optimization;

[0094] Based on the random forest classification model in step 3.3, a multi-classifier ensemble system is constructed by combining a support vector machine classifier and a gradient boosting decision tree classifier to improve classification accuracy and robustness. Each classifier is trained using the same training samples but employs different feature weighting strategies: the random forest classifier treats all features with equal weight; the support vector machine classifier emphasizes nonlinear feature relationships through kernel functions; and the gradient boosting decision tree classifier emphasizes difficult-to-classify samples through iterative optimization. The ensemble decision uses a weighted voting mechanism, with weights dynamically adjusted based on the classification accuracy of each classifier on the validation set. The weight calculation rule is the accuracy of a single classifier divided by the sum of the accuracies of all classifiers. For each image object, the three classifiers output predicted probabilities for four categories: generally disturbed surface, engineering excavation face, engineering accumulation, and undisturbed area. The final decision selects the category with the highest weighted probability as the classification result for the image object, yielding the ensemble classification prediction result.

[0095] Step 3.5: Post-processing and quality control of classification results;

[0096] Based on the integrated classification prediction results from step 3.4, morphological post-processing is performed on the classification results to eliminate salt-and-pepper noise and unreasonable classification patches. Opening operations are used to remove small-area noise patches, with the structuring element set to a 3x3 circle; closing operations are used to fill small-area holes, with the structuring element set to a 5x5 circle. The size of the structuring element is determined based on the minimum mapping unit, generally ranging from 3x3 to 7x7 pixels. Simultaneously, connectivity analysis is performed, merging objects smaller than the minimum mapping unit into adjacent, larger objects. The minimum mapping unit is determined based on the application accuracy requirements, generally ranging from 0.1 to 0.5 hectares. A classification confidence assessment mechanism is introduced, calculating the confidence score for each classification result. The confidence score is based on the consistency of the multi-classifier prediction probabilities. Regions with confidence scores below a set threshold are marked as uncertain regions requiring manual review. The confidence score threshold is determined based on quality control requirements, generally ranging from 0.7 to 0.9, resulting in preliminary perturbation unit classification results data.

[0097] Step 4: Based on the preliminary perturbation unit classification results, eliminate classification noise and extract the perturbation unit boundaries; smooth the perturbation unit boundaries to obtain the perturbation unit vector boundary data;

[0098] Step 4.1, morphological filtering processing;

[0099] Based on the preliminary perturbation unit classification results from step 3, morphological opening is used to eliminate isolated noise points in the classification results. Morphological opening is a combination of erosion and dilation, which can eliminate isolated regions smaller than the structuring element while maintaining the shape and area of ​​larger regions. The structuring element is circular, and its size is determined based on the typical size of the noise points, generally ranging from three to five pixels. Morphological closing is then used to fill small holes in the classification results. Closing is a combination of dilation and erosion, which can fill holes smaller than the structuring element while maintaining the shape and area of ​​larger regions, resulting in the morphologically filtered data.

[0100] Step 4.2, edge detection and boundary extraction;

[0101] Based on the morphological filtering results from step 4.1, the Canny edge detection algorithm is used to extract the boundaries of perturbation cells. The Canny edge detection algorithm includes four steps: Gaussian filtering, gradient calculation, non-maximum suppression, and dual-threshold detection. Gaussian filtering smooths the image and reduces the impact of noise on edge detection. The standard deviation of the Gaussian kernel is set according to the image noise level, typically one to two. Gradient calculation uses the Sobel operator to calculate the gradients in the horizontal and vertical directions. Dual-threshold detection uses two parameters: a high threshold and a low threshold. The high threshold is used to identify strong edge points, and the low threshold is used to identify weak edge points. The threshold setting rule is based on the boundary intensity distribution; the high threshold is generally taken as the 80ths percentile of the gradient magnitude distribution, and the low threshold is taken as one-third to one-half of the high threshold. A continuous boundary line of the perturbation cells is generated using an edge-connection method, resulting in boundary line data.

[0102] Step 4.3, Vectorization and Boundary Smoothing;

[0103] Based on the boundary line data from step 4.2, a vectorization method is used to convert the raster boundaries into vector boundaries. During vectorization, the pixel coordinates of the edge chains are converted to geographic coordinates. Then, the Douglas-Peucker algorithm is used to simplify the boundary points, reducing the number of boundary points while maintaining the boundary shape characteristics. The simplification tolerance setting rule is determined according to the boundary accuracy requirements, generally taking one to two times the image spatial resolution. Finally, vector polygons are generated, with each polygon representing a perturbation unit. The simplified boundaries are then smoothed using a Bézier curve smoothing method to generate more natural and smoother perturbation unit boundaries, resulting in smoothed vector boundary data.

[0104] Step 4.4, Spatial overlay analysis and quality check;

[0105] Based on the smoothed vector boundary data from step 4.3 and the project construction layout vector map from step 1.1, spatial overlay analysis is performed to eliminate obviously misjudged areas. The purpose of spatial overlay analysis is to eliminate disturbance units located outside the project boundary line and disturbance units overlapping with known non-disturbed areas. Spatial overlay analysis uses vector intersection operations to retain disturbance units located within the project boundary line that do not overlap with known non-disturbed areas. Simultaneously, topology checks are performed to verify the spatial integrity of the disturbance unit vector data, including polygon closure checks, polygon overlap checks, and polygon gap checks. For topology errors discovered during the checks, an automatic repair method is used to process them, resulting in disturbance unit vector boundary data.

[0106] Specifically, the polygon closure check verifies whether the coordinates of the starting node and the ending node are completely consistent by traversing the boundary nodes of each disturbance unit polygon, and checks whether the boundary line segments form a complete closed loop. For any unclosed polygons found, the coordinates of the ending node are automatically adjusted to the coordinates of the starting node to ensure the complete closure of the polygon boundary.

[0107] The polygon overlap detection method uses spatial intersection analysis to compare the geometric relationships of adjacent disturbed polygons one by one. By calculating the intersection points of polygon boundaries and the area of ​​overlapping regions, it identifies pairs of overlapping polygons. For the detected overlapping regions, the overlapping parts are merged into the polygons with higher priority according to the size of the overlapping area and the priority of polygon attributes. At the same time, the boundary of the other polygon is adjusted to eliminate the overlap.

[0108] Polygon gap inspection constructs an adjacency graph of the perturbation unit polygons, analyzes the distance relationship between the boundaries of adjacent polygons, identifies areas where the boundary spacing exceeds a preset gap threshold, and fills the gaps by boundary extension or interpolation polygon generation based on the geometry of the gaps and the attribute characteristics of the surrounding polygons, ensuring the spatial continuity and integrity of the perturbation unit partitioning results.

[0109] The automatic repair method includes three levels: node coordinate correction, boundary segment reconstruction, and polygon geometry reconstruction. Node coordinate correction adjusts the position of offset nodes using the least squares method. Boundary segment reconstruction generates smooth connecting segments using spline interpolation. Polygon geometry reconstruction restores the reasonable geometric shape of the polygon using convex hull algorithm and morphological operations.

[0110] Step 5: Based on the perturbation unit vector boundary data and the preliminary perturbation unit classification results, determine the perturbation type of each perturbation unit; based on the determination results, perform area calculation and attribute labeling to obtain perturbation unit data with attributes;

[0111] Step 5.1, Determination of perturbation unit category and feature statistics;

[0112] Based on the perturbation unit vector boundary data from step 4 and the preliminary perturbation unit classification results from step 3, the perturbation type of each perturbation unit is determined using a category determination method.

[0113] A spatial correspondence between the perturbation unit vector polygon and the raster classification results is established. Through vector raster overlay analysis, all raster pixels contained within the boundary of each vector polygon are determined, and a mapping table between vector polygon identifiers and pixel coordinates is established. The number distribution of each category of pixels within each vector polygon is counted one by one. All pixels within the polygon are traversed, and the category label of each pixel in the integrated classification prediction results in step 3 is read. The number of pixels in four categories—generally disturbed surface, engineering excavation face, engineering accumulation body, and undisturbed area—is counted separately, and the proportion of each category of pixels to the total number of pixels within the polygon is calculated. The dominant perturbation type of the perturbation unit is determined by majority voting, and the category with the most pixels is taken as the perturbation unit. The final disturbance type of the element is determined by classifying the elements according to the priority order of engineering excavation face, engineering accumulation body, general disturbed surface, and undisturbed area when there are two or more categories with the same number of elements. Based on the topographic factor raster data in step 1.5, the slope values ​​of all elements within each disturbance unit are extracted, and the average slope of the disturbance unit is calculated using the arithmetic mean method, which is used as the topographic feature attribute of the disturbance unit. The determined disturbance type and the calculated average slope are assigned as attribute information to the corresponding vector polygon to establish a type attribute data table for the disturbance unit. The average slope, as an important factor affecting soil loss, provides key input parameters for subsequent soil loss measurement, thus obtaining the type attribute data of the disturbance unit.

[0114] Step 5.2, Calculation of general disturbed surface area;

[0115] Based on the disturbance units identified in step 5.1 as having a general disturbed surface, their areas are calculated using an area calculation method. The area calculation employs a vector polygon area algorithm, which uses the geographic coordinates of the polygon vertices and a shoelace algorithm to calculate the polygon area. The shoelace algorithm decomposes the polygon into several triangles and sums the areas of each triangle to obtain the total polygon area. The unit of the area calculation result is square meters or hectares, and the unit conversion is performed as needed. For disturbance units with complex shapes, an adaptive mesh subdivision method is used to improve the accuracy of the area calculation. The mesh subdivision density is set according to the polygon complexity, generally taking one-quarter to one-half of the original pixel size, ensuring that the calculation results meet the accuracy requirements for soil erosion measurement, thus obtaining the area calculation result data.

[0116] Step 5.3, Calculation of the volume of the excavated face accumulation;

[0117] Based on the disturbance units identified in step 5.1 as either engineering excavation faces or engineering accumulation bodies, a volumetric calculation method is used to calculate their volumes. The volume calculation is based on the terrain change feature layer data from step 2.3. The absolute values ​​of terrain changes for all pixels within the disturbance unit's range are statistically analyzed, multiplied by the pixel area, and summed to obtain the total volume of the disturbance unit. During the volume calculation process, for each pixel within the disturbance unit's range, the absolute value of its terrain change is multiplied by the pixel area, and all these products are summed to obtain the total volume of the disturbance unit. The unit of the volume calculation results is cubic meters, and unit conversion is performed as needed. To improve the accuracy of the volume calculation, a three-dimensional interpolation method is used to refine the digital elevation model. The interpolation method selection rule is determined based on terrain complexity: bilinear interpolation is used for flat areas, and cubic spline interpolation is used for complex terrain, ensuring the reliability of the volume calculation results and obtaining the volume calculation result data.

[0118] Step 5.4, Classification of disturbance unit scale levels;

[0119] Based on the area calculation results from step 5.2 and the volume calculation results from step 5.3, a scale classification method is used to divide the disturbed units into three levels: large, medium, and small, according to the guidelines for soil loss measurement in production and construction projects. For general disturbed surfaces, the scale classification is based on area thresholds: areas greater than five hectares are large, areas between five hectares and 0.5 hectares are medium, and areas less than 0.5 hectares are small. For engineering excavation faces and engineering accumulation bodies, the scale classification is based on volume thresholds: volumes greater than 100,000 cubic meters are large, volumes between 100,000 and 10,000 cubic meters are medium, and volumes less than 10,000 cubic meters are small. The threshold setting rules are determined according to the project type and disturbance intensity level, and can be adjusted according to actual conditions. The scale classification results provide an important basis for subsequent sampling of typical disturbed units and soil loss measurement, obtaining scale attribute data.

[0120] Step 5.5, Spatial location information extraction and attribute annotation;

[0121] Based on the perturbation unit vector boundary data from step 4, a spatial location extraction method is used to obtain the spatial location information of each perturbation unit. The spatial location information includes the centroid coordinates, the minimum bounding rectangle range, and the administrative division to which the perturbation unit belongs. The centroid coordinates are obtained using a polygon centroid calculation method; the minimum bounding rectangle range is obtained using a minimum area bounding rectangle algorithm; and the administrative division to which the perturbation unit belongs is obtained by spatially connecting it to the administrative division vector map. Based on the perturbation unit type attribute data from step 5.1, the scale level attribute data from step 5.4, the area calculation results from step 5.2, the volume calculation results from step 5.3, and the spatial location information, an attribute-assigned method is used to generate an attribute-bearing perturbation unit delineation vector map. During the attribute assignment process, attribute fields are added to each perturbation unit vector polygon, including perturbation unit number, perturbation type, scale level, area, volume, average slope, centroid longitude, centroid latitude, and the administrative division to which the perturbation unit belongs. The disturbance unit number is automatically generated using a combination of project code and serial number. The numbering rule is project code plus year plus four-digit serial number to ensure the uniqueness and traceability of the number, thus obtaining disturbance unit data with attributes.

[0122] Step 6: Divide the multi-period remote sensing images into disturbance units, and perform change detection based on the division results; perform spatiotemporal analysis on the detection results to obtain spatiotemporal change data of the disturbance units;

[0123] Step 6.1, multi-stage disturbance unit division processing;

[0124] Based on the multi-period remote sensing image data obtained in step 1, the processing steps 1 through 5 are repeated to obtain the perturbation unit division results for each period. The time interval between the multi-period images is determined according to the project construction progress and monitoring needs, generally ranging from one to three months. The perturbation unit division results for each period have the same attribute field structure, facilitating time-series comparative analysis. To ensure the consistency of the multi-period data, a unified spatial reference system, resolution specification, and processing parameters are used to obtain the multi-period perturbation unit division dataset.

[0125] Step 6.2, Spatial Overlay Change Detection;

[0126] Based on the multi-phase perturbation unit division dataset from step 6.1, spatial overlay analysis is used for change detection. Spatial overlay analysis employs vector intersection operations to overlay the vector maps of perturbation units from the earlier and later phases, generating an overlay result layer. Each polygon in the overlay result layer contains perturbation type attributes from both the earlier and later phases; by comparing the attribute values ​​of the two phases, the changes in perturbation units can be identified. The spatial matching threshold for change detection is set according to the positioning accuracy, generally one to two times the pixel size, to ensure the accuracy of change detection and obtain the spatial overlay change detection results.

[0127] Step 6.3, Determining and quantifying the type of change;

[0128] Based on the spatial overlay change detection results in step 6.2, a change type determination method is used to identify the change type of the perturbation unit. The change types include four categories: addition, disappearance, type change, and no change. Addition refers to a perturbation unit that exists in a later period but not in an earlier period; disappearance refers to a perturbation unit that exists in an earlier period but not in a later period; type change refers to a perturbation unit that exists in both periods but its perturbation type changes; no change refers to a perturbation unit that exists in both periods and has the same perturbation type. The degree of change of the perturbation unit is quantified using the change area and change volume calculation methods. For added perturbation units, the change area equals their area, and the change volume equals their volume; for disappeared perturbation units, the change area equals the negative of their area, and the change volume equals the negative of their volume; for perturbation units undergoing type change, the area difference and volume difference before and after the change are calculated to obtain the change type and quantified change data.

[0129] Step 6.4, Time-series trajectory analysis;

[0130] Based on the change type and quantification data from step 6.3, a time-series trajectory analysis method is used to track the evolution of the disturbance units. A unique tracking identifier is assigned to each disturbance unit. The correspondence between preceding and following disturbance units is established based on the spatial overlap and attribute similarity between adjacent disturbance units. Spatial overlap is calculated using the intersection-union ratio (IUGR), which is equal to the ratio of the area of ​​the intersection to the area of ​​the union of two polygons. Attribute similarity is calculated based on the consistency of attributes such as disturbance type and scale level. When both spatial overlap and attribute similarity are higher than preset thresholds, the preceding and following disturbance units are considered to be different periods of the same disturbance unit, and the tracking identifier remains unchanged; otherwise, it is considered that a disturbance unit has been added or disappeared. The rules for setting the overlap threshold and similarity threshold are determined according to the change detection accuracy requirements, generally set to 0.5 and 0.8 respectively, to obtain the time-series trajectory data of the disturbance units.

[0131] Step 6.5, Spatiotemporal attribute generation and trend analysis;

[0132] Based on the time-series trajectory data of the disturbance units obtained in step 6.4, spatiotemporal evolution attributes of the disturbance units are generated using a spatiotemporal attribute generation method. These attributes include the first appearance time, last existence time, duration, and change history of the disturbance units. The first appearance time is the time when the disturbance unit is first identified in the time-series image; the last existence time is the time when the disturbance unit is last identified in the time-series image; the duration is the difference between the last existence time and the first appearance time; and the change history records the disturbance type and scale level of the disturbance unit in each period. A simple trend analysis method is used to predict the trend of area and volume changes of the disturbance units. A trend model is established using linear regression or multinomial fitting to predict the disturbance state for the next one to two monitoring periods. The selection rule for the trend analysis model is determined based on the change patterns of historical data: linear regression is used when the changes are stable, and quadratic multinomial fitting is used when the changes are complex, thus obtaining the spatiotemporal change data of the disturbance units.

[0133] Step 7: Based on the attributed perturbation unit data and the spatiotemporal variation data of the perturbation unit, perform standardized result output and accuracy verification to obtain standardized perturbation unit result data;

[0134] Step 7.1, standard format data output and topology check;

[0135] Based on the attribute-bearing perturbation unit data from step 5, a data format conversion method is used to output result data in multiple standard formats. The data format conversion supports outputting common GIS formats such as Shapefile, GeoJSON, and KML, as well as formats conforming to the National Geospatial Data Exchange Standard, facilitating data exchange with different GIS platforms and soil and water conservation management platforms. The integrity of attribute fields and data accuracy are maintained during the data format conversion process. Based on the output result data, a topology check method is used to verify the spatial integrity of the perturbation unit vector data. Topology checks include polygon closure checks, polygon overlap checks, and polygon gap checks. Polygon closure checks verify whether the boundaries of each perturbation unit polygon form closed loops; polygon overlap checks verify whether there are overlapping areas between adjacent perturbation units; polygon gap checks verify whether there are uncovered gap areas between adjacent perturbation units. For topology errors discovered during the checks, an automatic repair method is used to process them, resulting in standard format result data.

[0136] Step 7.2, verify the accuracy of the partitioning results;

[0137] Based on field survey data or high-precision reference data of the project area, an accuracy verification method is used to evaluate the classification accuracy of the perturbation unit segmentation results. Accuracy verification employs the confusion matrix method, comparing the segmentation results with the reference data pixel-by-pixel or object-by-object, counting the number of correctly classified and incorrectly classified items for each category, and calculating accuracy indices such as overall accuracy, producer accuracy, user accuracy, and Kappa coefficient. Overall accuracy reflects the average classification accuracy across all categories; producer accuracy reflects the probability of a particular category being correctly classified; user accuracy reflects the probability that a pixel or object classified as a particular category actually belongs to that category; the Kappa coefficient comprehensively considers classification accuracy and random consistency, and is an important indicator for evaluating the reliability of the classification results. Simultaneously, a cross-validation method is used to evaluate the generalization ability of the segmentation method. The training sample set is randomly divided into several subsets, and one subset is selected as the validation set each time, while the remaining subsets are used as the training set. Model training and accuracy verification are repeated. The validation set ratio is determined based on the total sample size, generally between 20% and 30%, to obtain the accuracy verification results.

[0138] Step 7.3: Application requirements for accuracy evaluation and result verification;

[0139] Based on the accuracy verification results from step 7.2, an accuracy evaluation method is used to determine whether the perturbation unit division results meet the application requirements. The accuracy evaluation is based on the accuracy requirements of the guidelines for soil loss measurement in production and construction projects. When the overall accuracy is higher than 85% and the Kappa coefficient is higher than 0.8, the division results are considered to meet the application requirements and can be used for subsequent sampling of typical perturbation units and soil loss measurement; otherwise, the reasons for insufficient accuracy need to be analyzed, classification parameters adjusted, or training samples supplemented, and the perturbation unit division re-performed. The accuracy threshold setting rules are determined according to the application accuracy requirements and industry standards, and can be adjusted according to the specific needs of the project. Based on the verified perturbation unit division results, a metadata file for the result data is created using a metadata generation method. The metadata file records the basic information, data quality information, spatial reference information, and processing information of the result data. Basic information includes data name, data format, data range, and data time; data quality information includes classification accuracy, spatial accuracy, and attribute completeness; spatial reference information includes coordinate system, projection method, and spatial resolution; and processing information includes data source, processing method, processing parameters, and processing time. Metadata files provide necessary descriptive information for the management, sharing, and reuse of output data, ensuring data traceability and reproducibility, and obtaining standardized perturbation unit output data.

[0140] This invention, through the organic combination of the above seven steps, achieves intelligent, high-precision, and predictable division of disturbance units. Compared with traditional manual division methods, this invention has advantages such as high efficiency, stable accuracy, strong objectivity, and good traceability. It can meet the actual needs of soil and water conservation monitoring and management in large-scale production and construction projects, provide reliable basic data for soil loss measurement, and provide strong technical support for precise soil and water conservation control and forward-looking decision-making.

[0141] In one embodiment of the present invention, a specific example is provided:

[0142] The method described in this invention is applied and verified using dynamic monitoring of soil and water conservation during the construction period of a large-scale highway construction project as an example. The project is located in a city in Province A, with a total length of approximately 120 kilometers, traversing various terrains including mountains, hills, and plains. The construction process generates significant surface disturbance due to excavation, filling, and temporary soil stockpiling. The project owner commissioned a soil and water conservation monitoring unit to conduct dynamic monitoring of soil and water conservation during the construction period, requiring multi-phase division and change detection of disturbance units within the project area.

[0143] Using the method described in this invention, high-resolution remote sensing images and multispectral images of the project area for March, June, and September 2025, as well as digital elevation model data before and after construction, were acquired. Following the steps described in this invention, the three images were divided into disturbance units, and change detection and dynamic updates were performed.

[0144] Table 1 shows an example of the results of dividing a typical section of disturbance units:

[0145] Table 1: Example data of disturbance unit division results for a typical section;

[0146]

[0147] Table 2 shows an example of the detection results of disturbance unit changes in a typical section:

[0148] Table 2: Example data of disturbance unit change detection results in a typical section;

[0149]

[0150] Application results show that the present invention can automatically, quickly and accurately divide disturbance units and effectively identify the spatiotemporal evolution characteristics of disturbance units, providing reliable data support for dynamic monitoring and optimization of soil and water conservation measures during the construction period.

[0151] The embodiments of the present invention have been described above. However, the embodiments are not limited to the specific implementation methods described above. The specific implementation methods described above are merely illustrative and not restrictive. Those skilled in the art can make more equivalent embodiments under the guidance of the present embodiments, and all of them are within the protection scope of the present embodiments.

Claims

1. A method for dividing disturbance units in a soil loss measurement process, characterized in that, include: Multiple remote sensing images and multispectral images of the project area were collected. Spatial registration preprocessing was performed on the multispectral images to obtain radiometrically normalized fused image data. Based on radiometrically normalized fused image data, multidimensional perturbation features are extracted to obtain a standardized multidimensional feature dataset. Multi-scale image segmentation is performed based on a standardized multi-dimensional feature dataset, and perturbation units are intelligently identified by combining the segmentation results to obtain preliminary perturbation unit classification results data. Based on the preliminary perturbation unit classification results, classification noise is eliminated and perturbation unit boundaries are extracted; the perturbation unit boundaries are then smoothed to obtain perturbation unit vector boundary data. Based on the vector boundary data of the disturbance unit and the preliminary classification results of the disturbance unit, the disturbance type of each disturbance unit is determined; Based on the judgment results, area calculation and attribute labeling are performed to obtain attribute-perturbed cell data; Disturbance units are divided into multiple remote sensing images, and change detection is performed based on the division results; the detection results are then subjected to spatiotemporal analysis to obtain spatiotemporal change data of the disturbance units. Based on attributed perturbation unit data and spatiotemporal variation data of perturbation units, standardized results are output and accuracy is verified to obtain standardized perturbation unit result data.

2. The method for dividing disturbance units in the soil loss measurement process according to claim 1, characterized in that, The spatial registration preprocessing of the multispectral image includes: Radiometric calibration is performed using sensor calibration parameters to convert the digital quantization values ​​of multispectral images into physical radiance values. A fast atmospheric correction algorithm is used to eliminate the influence of atmospheric scattering and absorption on surface reflectivity. Subpixel-level spatial registration is performed using an automatic registration algorithm based on feature point matching. Scale-invariant feature transformation feature points are extracted from the two multispectral images, and the nearest neighbor distance ratio method is used for feature point matching. False matching points are removed using a random sampling consensus algorithm. A principal component transformation fusion algorithm is used to perform principal component transformation on multispectral images, extract the first principal component, perform histogram matching between high-resolution panchromatic images and the first principal component, replace the first principal component with the matched panchromatic images, and perform inverse principal component transformation. The terrain factors were extracted using a third-order inverse distance weighted difference algorithm, and the radiometric consistency of the multi-temporal images was processed using a relative radiometric normalization method to obtain radiometrically normalized fused image data.

3. The method for dividing disturbance units in the soil loss measurement process according to claim 1, characterized in that, The multidimensional perturbation feature extraction includes: Based on radiometrically normalized fused image data, normalized vegetation index and soil-adjusted vegetation index were calculated to obtain spectral feature layer data. Based on radiometrically normalized fused image data, texture features are extracted at multiple scales using the gray-level co-occurrence matrix method to obtain texture feature layer data. Based on the digital elevation model data before and after construction, the digital elevation model is preprocessed with median filtering to calculate the terrain change and obtain terrain change feature layer data. The Z-score normalization method was used to normalize the spectral feature layer data, texture feature layer data, and terrain change feature layer data. The normalization results were then fused, and principal component analysis was used to reduce the dimensionality of the multidimensional features to obtain a standardized multidimensional feature dataset.

4. The method for dividing disturbance units in the soil loss measurement process according to claim 1, characterized in that, The multi-scale image segmentation and perturbation identification unit based on the segmentation results includes: Based on a standardized multidimensional feature dataset, a multi-scale segmentation algorithm is used to generate image objects, resulting in image segmentation object data. Based on image segmentation object data, a classification training sample set was constructed using a sample selection method. The training sample set includes samples from four categories: general disturbed surface, engineering excavation face, engineering accumulation body, and undisturbed area. An improved synthetic minority class oversampling technique was used to augment the samples, resulting in an augmented training sample dataset. The expanded training sample dataset is trained using a random forest classifier to obtain a trained random forest classification model. Based on the trained random forest classification model, a multi-classifier ensemble system is constructed by combining a support vector machine classifier and a gradient boosting decision tree classifier. The ensemble decision adopts a weighted voting mechanism to obtain the ensemble classification prediction results. Morphological post-processing is performed on the integrated classification prediction results. Opening operation is used to remove small noise patches, closing operation is used to fill small holes, and a classification confidence evaluation mechanism is introduced to obtain preliminary perturbation unit classification result data.

5. The method for dividing disturbance units in the soil loss measurement process according to claim 1, characterized in that, The elimination of classification noise includes: Based on the preliminary perturbation unit classification results, the morphological opening operation method is used to eliminate isolated noise points in the classification results, and the morphological filtering results are obtained. Based on the morphological filtering results, the Canny edge detection algorithm is used to extract the boundaries of perturbation cells. The algorithm includes four steps: Gaussian filtering, gradient calculation, non-maximum suppression, and double threshold detection. The edge connection method is used to generate continuous boundary lines of perturbation cells, thus obtaining boundary line data.

6. The method for dividing disturbance units in a soil loss measurement process according to claim 1, characterized in that, The smoothing process for the boundary of the disturbance unit includes: Based on boundary line data, a vectorization method is used to convert raster boundaries into vector boundaries, the cell coordinates of the edge chain are converted into geographic coordinates, the Douglas-Peucker algorithm is used to simplify the boundary points, and the Bézier curve smoothing method is used to smooth the simplified boundaries to obtain smooth vector boundary data. Spatial overlay analysis is performed based on smooth vector boundary data. Disturbance cells located outside the project boundary and those overlapping with known non-disturbance areas are removed. Topology checks are then performed, and automatic repair methods are used to process any topology errors found, resulting in the obtained vector boundary data of the disturbance cells.

7. The method for dividing disturbance units in a soil loss measurement process according to claim 1, characterized in that, The determination of the disturbance type for each disturbance unit includes: Establish the spatial correspondence between the perturbation unit vector polygon and the raster classification results; Count the number distribution of each category of pixels within each vector polygon, traverse all pixels within the polygon range, and read the category label of each pixel in the integrated classification prediction result; The dominant disturbance type of the disturbance unit is determined by majority voting. Extract the slope values ​​of all pixels within the range of each perturbation cell, and calculate the average slope of the perturbation cell using the arithmetic mean method; The determined disturbance type and calculated average slope are assigned as attribute information to the corresponding vector polygon to establish a disturbance unit type attribute data table, thereby obtaining the disturbance unit type attribute data.

8. The method for dividing disturbance units in the soil loss measurement process according to claim 1, characterized in that, The area calculation and attribute labeling based on the judgment result include: For disturbance units with a general disturbance type of surface, the vector polygon area algorithm is used, which calculates the polygon area based on the geographic coordinates of the polygon vertices using the shoelace algorithm. For disturbance units whose disturbance type is engineering excavation face and engineering accumulation, based on the topographic change feature layer data, the absolute value of the topographic change of all pixels within the disturbance unit is calculated, multiplied by the pixel area, and summed to obtain the total volume of the disturbance unit; The disturbance cells are classified into three levels: large, medium, and small, using a scale-level determination method. An attribute assignment method is then used to generate a vector map of the disturbance cells with attributes.

9. The method for dividing disturbance units in a soil loss measurement process according to claim 1, characterized in that, The process involves dividing multi-period remote sensing images into disturbance units, detecting changes based on the division results, and performing spatiotemporal analysis on the detection results. Based on multiple periods of remote sensing image data, the processing flow of multi-source remote sensing data acquisition and spatial registration preprocessing, multi-dimensional perturbation feature extraction, object-oriented image segmentation and supervised learning classification, boundary optimization and vectorization processing, area and volume calculation and attribute annotation was repeated to obtain the perturbation unit division results for each period. Spatial overlay analysis was used to detect changes by overlaying vector maps of disturbance units from the early and later periods to generate an overlay result layer. The changes in disturbance units were identified by comparing the attribute values ​​of the two periods. The change type determination method is used to identify the change type of the disturbance element, and the change area and change volume calculation method is used to quantify the degree of change of the disturbance element; The evolution of perturbation units is tracked using time-series trajectory analysis. A unique tracking identifier is assigned to each perturbation unit. Based on the spatial overlap and attribute similarity of perturbation units in adjacent periods, the correspondence between perturbation units in the early and late periods is established. Spatiotemporal evolution attributes of disturbance units are generated using a spatiotemporal attribute generation method.

10. The method for dividing disturbance units in a soil loss measurement process according to claim 1, characterized in that, The standardization of output results and accuracy verification includes: The data is output in a standard format using data format conversion methods. The spatial integrity of the disturbance element vector data was verified using a topology check method. Based on field survey data of the project area, an accuracy verification method was used to evaluate the classification accuracy of the disturbance unit division results; The quality of the perturbation unit partitioning result is judged by the accuracy evaluation method. When the overall accuracy of the perturbation unit partitioning result is higher than the preset accuracy threshold and the Kappa coefficient is higher than the preset Kappa threshold, the partitioning result is considered to meet the application requirements. Metadata files for the output data are created using metadata generation methods.

Images