Erosion risk assessment method and system based on erosion degradation of black soil

Abstract

The application discloses an erosion degradation-based black soil erosion risk assessment method and system, relates to the technical field of black soil erosion assessment, and comprises the following steps: collecting high-resolution soil structure tomographic images of a target black soil area, long-term climatic and hydrological time series data, a digital model of surface microtopography, and a historical vector diagram of land use change; analyzing the soil structure tomographic images for pore connectivity, extracting a soil aggregate breakdown index and a vertical crack development density, calculating an effective rainfall erosion force factor and a seasonal freeze-thaw cycle intensity factor, and fusing the microtopography and land use change data to generate a surface disturbance intensity layer; inputting multiple parameters into a black soil erosion process coupling model, and generating an erosion degradation risk distribution map under different future scenarios through spatial iteration operation; and the method finely represents the microstructure of black soil and the surface composite disturbance characteristics, and improves the fine level and scene adaptation capability of black soil erosion risk assessment.

Description

Technical Field

[0001] This invention belongs to the field of black soil erosion assessment technology, specifically a method and system for assessing the risk of black soil erosion based on erosion and degradation. Background Technology

[0002] Conventional black soil erosion risk assessments often rely on macroeconomic factors such as rainfall erosion capacity, topographic slope, and land use type as the core assessment basis. Data sources mainly consist of conventional ground monitoring data and ordinary remote sensing data. The assessment process often adopts a single-factor overlay analysis model, and the treatment of seasonal freeze-thaw effects is often simplified. A coupled analysis framework of climate factors and soil characteristics has not been constructed.

[0003] Current assessment methods lack detailed analysis of the internal pore structure of black soil, failing to quantify the extent of soil aggregate damage caused by tillage and erosion, and also failing to characterize the water infiltration channels formed by the development of vertical fissures in black soil. The coupled calculation methods for freeze-thaw cycles and rainfall erosion are rather crude and do not accurately reflect the intrinsic mechanisms of black soil erosion and degradation. Surface disturbance analysis relies solely on topographic data or static land use data, without incorporating land use change history and tillage activity trajectories, thus failing to reflect the complex impact of multiple surface elements on soil stability. Furthermore, the assessment models lack spatial iterative computation capabilities and cannot extrapolate the distribution of black soil erosion risks under different future scenarios.

[0004] It is necessary to rely on soil structure tomography images to complete pore connectivity analysis, obtain soil aggregate destruction index and vertical fracture development density, and at the same time integrate surface micro-topography and land use change historical vector data to construct a surface disturbance intensity layer. Through spatial iterative calculation of coupled models, the risk distribution of black soil erosion and degradation under multiple scenarios can be obtained. Summary of the Invention

[0005] This invention aims to solve at least one of the technical problems existing in the prior art;

[0006] Therefore, this invention proposes a method for assessing the erosion risk of black soil based on erosion and degradation, including:

[0007] Collect multi-source monitoring data of the target black soil area. The multi-source monitoring data includes high-resolution soil structure tomography images, long-term climate and hydrological time series data, digital model of surface micro-topography, and historical vector map of land use change.

[0008] The high-resolution soil structure tomography image is subjected to pore connectivity analysis to extract the soil aggregate destruction index and vertical fissure development density. The soil aggregate destruction index represents the degree of structural degradation caused by tillage or erosion, and the vertical fissure development density represents the complexity of water infiltration channels.

[0009] Based on the aforementioned long-term climate and hydrological time-series data, the effective rainfall erosivity factor and the seasonal freeze-thaw cycle intensity factor are calculated.

[0010] By integrating the digital model of surface micro-topography with the historical vector map of land use change, a surface disturbance intensity layer is generated. The surface disturbance intensity layer comprehensively reflects the combined effects of slope and slope length, surface cover change frequency and cultivation activity trajectory on soil stability.

[0011] The soil aggregate destruction index, the vertical fissure development density, the effective rainfall erosivity factor, the seasonal freeze-thaw cycle intensity factor, and the surface disturbance intensity layer are input into the black soil erosion process coupling model. Through spatial iterative calculation of the model, an erosion and degradation risk distribution map of the target black soil area under different future scenarios is generated.

[0012] Furthermore, the high-resolution soil structure tomography image is subjected to pore connectivity analysis to extract the soil aggregate disruption index and vertical fracture development density, including:

[0013] A digital image segmentation algorithm is used to separate soil particles, pores, and organic matter regions from the high-resolution soil structure tomography image;

[0014] The separated pore regions are skeletonized and topologically analyzed to calculate the connectivity path length and bottleneck width of the pore network and identify the preferred flow paths.

[0015] Based on the geometric morphology and distribution uniformity of soil particle regions, the average equivalent diameter and shape factor dispersion of aggregates are calculated, and combined with historical benchmark images, the soil aggregate destruction index is quantified.

[0016] On a vertical slice, the number of pores with connectivity reaching a preset threshold within a unit area is counted, and the angle between their orientation and the vertical direction is measured. Pores with an angle less than a preset angle are selected, and their number is the vertical fracture development density.

[0017] The output is a thematic map of soil structure parameters, which includes pixel coordinates, the soil aggregate destruction index, and the corresponding relationship between the vertical fissure development density.

[0018] Furthermore, based on the aforementioned long-term climate and hydrological time-series data, the effective rainfall erosivity factor and the seasonal freeze-thaw cycle intensity factor are calculated, including:

[0019] Extract minute-by-minute rainfall intensity sequences and daily average temperature sequences from the long-term climate and hydrological time-series data;

[0020] For the minute-by-minute rainfall intensity sequence, the rainfall events are divided into events, the total kinetic energy and the maximum half-hour rainfall intensity of each rainfall event are calculated, and the total kinetic energy and the maximum half-hour rainfall intensity are accumulated to obtain the effective rainfall erosivity factor on an annual scale;

[0021] From the daily average temperature sequence, the days in which the daily maximum temperature is below freezing and the daily minimum temperature is above freezing are identified and defined as one freeze-thaw cycle. The total number of freeze-thaw cycles in a year is counted as the number of freeze-thaw days.

[0022] For each freeze-thaw cycle, the daily temperature difference amplitude is calculated, and the daily temperature difference amplitudes of all freeze-thaw cycles are weighted and accumulated, where the weight is the change in soil moisture content before and after each freeze-thaw cycle, and finally the seasonal freeze-thaw cycle intensity factor is calculated.

[0023] Establish a time-series database of the effective rainfall erosivity factor and the seasonal freeze-thaw cycle intensity factor, which are bound to timestamps and geographic coordinates.

[0024] Furthermore, by fusing the digital model of surface micro-topography with the land use change history vector map, a surface disturbance intensity layer is generated, including:

[0025] Hydrological analysis was performed on the digital model of the surface micro-topography to extract slope layer, slope length layer and topographic humidity index layer;

[0026] Spatiotemporal overlay analysis was performed on the land use change history vector map to calculate the number of land use type changes for each spatial unit within the assessment period, and a land cover change frequency layer was generated.

[0027] From the land use change history vector map, extract the patches related to cultivation activities, calculate the overlap of their boundary interannual movement trajectories and the consistency with the cultivation direction, and generate a cultivation disturbance accumulation layer.

[0028] Principal component analysis was used to reduce the dimensionality and fuse the slope layer, slope length layer, topographic humidity index layer, land cover change frequency layer, and cultivated disturbance accumulation layer.

[0029] The fused principal component values ​​are normalized and reclassified to generate a comprehensive surface disturbance intensity layer, whose pixel values ​​reflect the comprehensive intensity of the disturbance affecting surface stability.

[0030] Furthermore, the soil aggregate failure index, the vertical fracture development density, the effective rainfall erosivity factor, the seasonal freeze-thaw cycle intensity factor, and the surface disturbance intensity layer are input into the black soil erosion process coupling model. Through spatial iterative calculations of the model, erosion and degradation risk distribution maps of the target black soil area under different future scenarios are generated, including:

[0031] The soil aggregate destruction index and the vertical fracture development density are spatially superimposed, and the basic soil erodibility map is calculated using an empirical function.

[0032] The effective rainfall erosivity factor and the seasonal freeze-thaw cycle intensity factor are weighted and fused over time to generate a dynamic map of climate erosion potential.

[0033] The surface disturbance intensity layer is multiplied pixel by pixel with the soil erodibility baseline map to obtain the comprehensive erosion sensitivity map under the current state.

[0034] In the coupled model of the black soil erosion process, the comprehensive erosion sensitivity map is used as the initial field, and the dynamic map of climate erosion potential is used as the driving field that changes over time, to carry out iterative simulation based on physical processes.

[0035] Within each iteration step, the amount of surface material migration and the increment of soil structure weakening are calculated, and the comprehensive erosion sensitivity map is dynamically updated.

[0036] The simulation runs until a preset future time point, and the final total amount of soil structure weakening is compared with a preset threshold range to divide erosion and degradation risk zones of different levels, forming the erosion and degradation risk distribution map. The erosion and degradation risk distribution map expresses spatial differences in the form of risk level grids.

[0037] Furthermore, the soil aggregate failure index and the vertical fracture development density are spatially superimposed, and a basic soil erodibility map is calculated using an empirical function, including:

[0038] Establish a lookup table relationship between the soil aggregate destruction index and the soil cohesion attenuation coefficient;

[0039] Establish a lookup table relationship between the vertical fracture development density and the soil permeability rate correction coefficient;

[0040] The soil aggregate destruction index and the vertical fracture development density at each spatial location are input into the soil erodibility empirical function, which is expressed as a product of the cohesion attenuation coefficient, the permeability correction coefficient, and the standard soil erodibility factor.

[0041] Calculate the soil erodibility values ​​at all spatial locations within the target black soil region;

[0042] Spatially interpolate the soil erodibility values ​​at all spatial locations to generate a continuous baseline map of soil erodibility.

[0043] Furthermore, the simulation runs to a preset future time point, compares the final total amount of soil structure weakening with a preset threshold range, and divides different levels of erosion and degradation risk zones, including:

[0044] After the simulation, the total amount of soil structure weakening in each grid cell output by the coupled model of the black soil erosion process was extracted;

[0045] Multiple threshold ranges are set to characterize the severity of erosion and degradation, with each threshold range corresponding to a risk level;

[0046] The total amount of soil structure weakening in each grid cell is matched with the threshold range to determine the risk level of the grid cell.

[0047] Assign a unique identifier and color code to each risk level;

[0048] All grid cells are rendered according to their risk level and color code to generate the erosion and degradation risk distribution map with spatial distribution information.

[0049] Furthermore, the method for constructing the coupled model of the black soil erosion process includes:

[0050] Based on the principles of hydrodynamics and soil mechanics, a multi-process erosion dynamics equation set is established that couples surface runoff erosion, gully erosion and freeze-thaw stripping. The multi-process erosion dynamics equation set includes a runoff shear force sub-equation describing the sediment yield of surface runoff, a gully headward erosion sub-equation describing the gully development rate, and a freeze-thaw stripping sub-equation describing the stripping of surface soil caused by freeze-thaw action.

[0051] The multi-process erosion dynamics equations are spatiotemporally discretized using a numerical discretization method, transforming the continuous equations into difference equations suitable for grid computing, and determining the iteration step size and boundary conditions of the model.

[0052] A soil shear strength decay function is preset in the model. The soil shear strength decay function uses the soil aggregate failure index, the vertical crack development density and the seasonal freeze-thaw cycle intensity factor as dynamic input parameters to update the soil's erosion resistance in each iteration step.

[0053] A material migration and deposition function is preset in the model. The material migration and deposition function uses the micro-topography and surface cover reflected by the surface disturbance intensity layer as dynamic input parameters to simulate the transport and deposition process of eroded materials on the slope.

[0054] The key parameters in the model are calibrated, including the Manning roughness coefficient in the runoff shear force sub-equation, the soil erodibility coefficient in the gully headward erosion sub-equation, and the latent heat coefficient of water phase change in the freeze-thaw stripping sub-equation.

[0055] The calibrated model is coupled and integrated with the data structures of the soil aggregate destruction index, the vertical fissure development density, the effective rainfall erosivity factor, the seasonal freeze-thaw cycle intensity factor, and the surface disturbance intensity layer to complete the construction of the coupled model of the black soil erosion process.

[0056] Furthermore, the method also includes the step of generating a risk mitigation strategy plan based on the erosion and degradation risk distribution map:

[0057] Spatial areas with high and extremely high risk levels identified in the erosion and degradation risk distribution map are defined as priority remediation areas;

[0058] For the priority treatment area, the dominant erosion driving factors for the corresponding area are extracted from the multi-source monitoring data. The dominant erosion driving factors include dominant climate factors, dominant topographic factors and dominant anthropogenic disturbance factors.

[0059] Based on the combination type of the dominant erosion driving factors, the corresponding sets of engineering measures, biological measures, and management measures are matched from the strategy knowledge base, which stores the mapping relationship between different combinations of driving factors and control measures.

[0060] Based on the principles of spatial proximity and measure compatibility, the matched sets of engineering measures, biological measures, and management measures are optimized and combined to form a customized governance solution package for each priority governance area.

[0061] The customized governance packages for all priority governance areas are spatially correlated with the aforementioned erosion and degradation risk distribution map to generate a risk mitigation strategy plan map that includes spatial location, risk level, and specific measures.

[0062] The step of matching corresponding sets of engineering measures, biological measures, and management measures from the strategy knowledge base based on the combination type of the dominant erosion driving factors includes:

[0063] The dominant climate factor, dominant topographic factor, and dominant anthropogenic disturbance factor are encoded to form a three-dimensional driving factor code;

[0064] In the strategy knowledge base, a record is queried that exactly matches the three-dimensional driving factor code. The record stores a list of recommended engineering measures, a list of biological measures, and a list of management measures.

[0065] If no perfect match is found, a fuzzy matching algorithm is used to search for records in the knowledge base whose driving factor coding similarity exceeds a preset threshold, and to extract the list of measures associated with them.

[0066] The extracted lists of measures are deduplicated and merged to form preliminary sets of engineering measures, biological measures, and management measures.

[0067] Based on the on-site constraints of the priority treatment area, the preliminary set of measures is screened for feasibility, resulting in the final matching set of engineering measures, biological measures, and management measures.

[0068] Furthermore, the present invention also includes a black soil erosion risk assessment system based on erosion and degradation, the system including a memory, a processor, and a computer program stored in the memory and running on the processor, wherein the processor, when executing the computer program, implements the steps of the black soil erosion risk assessment method based on erosion and degradation described above.

[0069] Compared with the prior art, the beneficial effects of the present invention are:

[0070] By performing pore connectivity analysis on high-resolution soil structure tomography images, the soil aggregate destruction index and vertical fissure development density can be accurately extracted. The soil aggregate destruction index can directly characterize the structural degradation state of black soil caused by tillage and erosion, transforming the degree of degradation of the internal structure of black soil into a quantifiable parameter. The vertical fissure development density can fully depict the development and complexity of water infiltration channels within black soil, achieving a quantitative expression of the microstructural characteristics of black soil. This processing method can directly link changes in the internal structure of soil with the erosion and degradation process, fully presenting the influence of soil microstructure on water transport and soil stability, refining the characterization dimensions of soil structural attributes in erosion risk assessment, and making the characterization of erosion driving factors more closely aligned with the physical structural characteristics of black soil itself.

[0071] By integrating digital models of surface micro-topography with historical land use change vector maps to generate a surface disturbance intensity layer, this method can comprehensively integrate elements such as slope and slope length, frequency of surface cover change, and cultivation activity trajectories to form a quantitative layer reflecting the combined effects of multiple surface elements. This fully demonstrates the cumulative effect of various surface disturbance factors on soil stability, eliminating the limitations of single topographic data or static land use data in characterizing surface disturbance states. By inputting soil aggregate failure index, vertical fracture development density, effective rainfall erosiveness factor, seasonal freeze-thaw cycle intensity factor, and the surface disturbance intensity layer into a coupled model of black soil erosion processes, and through spatial iterative calculations of the model, the collaborative calculation of multi-source parameters can be completed. This outputs erosion and degradation risk distribution maps of the target black soil area under different future scenarios, enabling spatial extrapolation and multi-scenario presentation of erosion risk. This ensures that the assessment results of erosion risk are consistent with the actual evolution process of black soil erosion and degradation, improving the completeness and relevance of the erosion risk distribution characterization. Attached Figure Description

[0072] Figure 1 This is a flowchart illustrating the steps of the black soil erosion risk assessment method based on erosion and degradation described in this invention.

[0073] Figure 2 A flowchart for the analytical processing of pore connectivity;

[0074] Figure 3 A flowchart for generating a surface disturbance intensity layer. Detailed Implementation

[0075] The technical solution of the present invention will be clearly and completely described below with reference to the embodiments. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0076] See Figure 1 This invention provides a method for assessing the erosion risk of black soil based on erosion and degradation. The overall implementation scheme includes the following steps:

[0077] Multi-source monitoring data was collected from the target black soil region, including high-resolution soil structure tomography images, long-term climate and hydrological time-series data, digital models of surface micro-topography, and historical land use change vector maps. Pore connectivity analysis was performed on the high-resolution soil structure tomography images to extract the soil aggregate destruction index and vertical fissure development density. The soil aggregate destruction index characterizes the degree of structural degradation caused by tillage or erosion, while the vertical fissure development density characterizes the complexity of water infiltration channels. Based on long-term climate and hydrological time-series data, the effective rainfall erosivity factor and seasonal freeze-thaw cycle intensity factor were calculated. The surface micro-topography digital model and historical land use change vector map were fused to generate a surface disturbance intensity layer, which comprehensively reflects the combined effects of slope and slope length, surface cover change frequency, and tillage activity trajectory on soil stability. The soil aggregate destruction index, vertical fissure development density, effective rainfall erosivity factor, seasonal freeze-thaw cycle intensity factor, and surface disturbance intensity layer are input into the black soil erosion process coupling model. Through spatial iterative calculation of the model, the erosion and degradation risk distribution map of the target black soil area under different future scenarios is generated.

[0078] In one embodiment of the present invention, see [reference] Figure 2 This study employs a digital image segmentation algorithm to separate soil particles, pores, and organic matter regions from high-resolution soil structure tomography images. The separated pore regions are then skeletonized and subjected to topological analysis to calculate the connectivity path length and bottleneck width of the pore network, identifying preferred flow paths. Based on the geometric morphology and distribution uniformity of soil particle regions, the average equivalent diameter and shape factor dispersion of aggregates are calculated. Combined with historical baseline images, a soil aggregate destruction index is quantified. On vertical slices, the number of pores with connectivity reaching a preset threshold per unit area is counted, and the angle between their orientation and the vertical direction is measured. Pores with an angle smaller than a preset angle are selected, and their number represents the vertical fracture development density. The output is a thematic map of soil structure parameters containing pixel coordinates, the soil aggregate destruction index, and the corresponding relationship between the vertical fracture development density and the output.

[0079] Minute-by-minute rainfall intensity sequences and daily average temperature sequences were extracted from long-term climate and hydrological time-series data. For the minute-by-minute rainfall intensity sequences, rainfall events were segmented by event, and the total kinetic energy and maximum half-hour rainfall intensity of each event were calculated. The total kinetic energy and maximum half-hour rainfall intensity were accumulated to obtain the effective rainfall erosivity factor on an annual scale. From the daily average temperature sequences, dates in which the daily maximum temperature is below freezing and the daily minimum temperature is above freezing were identified and defined as one freeze-thaw cycle. The total number of freeze-thaw cycles per year was counted as the number of freeze-thaw days. For each freeze-thaw cycle, its daily temperature range amplitude was calculated. The daily temperature range amplitudes of all freeze-thaw cycles were weighted and accumulated, with the weight being the change in soil moisture content before and after each freeze-thaw cycle. Finally, the seasonal freeze-thaw cycle intensity factor was calculated. A time-series database of effective rainfall erosivity factors and seasonal freeze-thaw cycle intensity factors, bound to timestamps and geographic coordinates, was established.

[0080] In practical implementation, the processing of high-resolution soil structure tomography images and long-term climate and hydrological time-series data is detailed through specific example scenarios and data comparisons. In one example scenario, the target black soil region is located in the Northeast Plain of China. The acquired high-resolution soil structure tomography images have a spatial resolution of 50 micrometers and a coverage depth of 1 meter. The long-term climate and hydrological time-series data spans the past 30 years, including minute-by-minute rainfall records and daily average temperature records. By comparing images and data from different years or different plots, the degradation of soil structure and changes in climate-driven erosion potential can be quantified.

[0081] In the specific implementation, a digital image segmentation algorithm is used to separate soil particle regions, pore regions, and organic matter regions from high-resolution soil structure tomography images. The digital image segmentation algorithm is based on threshold segmentation and region growing methods. The separated pore regions undergo skeletonization and topological analysis. Skeletonization uses a median transformation algorithm, and topological analysis involves calculating the length of connected paths and bottleneck widths in the pore network to identify preferred flow paths. Based on the geometric morphology and distribution uniformity of the soil particle regions, the average equivalent diameter and shape factor dispersion of aggregates are calculated. Combined with historical baseline images, the soil aggregate destruction index is quantified. On vertical slices, the number of pores with connectivity reaching a preset threshold per unit area is counted, and the angle between the pore orientation and the vertical direction is measured. Pores with an angle smaller than a preset angle are selected, and their number represents the vertical fracture development density. The output is a thematic map of soil structure parameters containing pixel coordinates, the soil aggregate destruction index, and the corresponding relationship between the vertical fracture development density and the data. It is understandable that the parameter settings of image segmentation algorithms need to be adjusted according to the soil type. For example, in clayey black soil, the gray threshold of threshold segmentation may be set to 0.3 to accurately separate pores.

[0082] In some embodiments, minute-by-minute rainfall intensity sequences and daily average temperature sequences are extracted from long-term climate and hydrological time-series data, with data sources including meteorological station observations and remote sensing inversion products. For the minute-by-minute rainfall intensity sequences, rainfall events are segmented by event, with a segmentation criterion of rainfall intervals exceeding 6 hours considered independent events. The total kinetic energy and maximum half-hour rainfall intensity of each rainfall event are calculated, and the total kinetic energy and maximum half-hour rainfall intensity are summed to obtain the effective rainfall erosivity factor on an annual scale. From the daily average temperature sequence, dates where the daily maximum temperature is below freezing and the daily minimum temperature is above freezing are identified and defined as one freeze-thaw cycle. The total number of freeze-thaw cycles in a year is counted as the number of freeze-thaw days. For each freeze-thaw cycle, the daily temperature difference amplitude of the freeze-thaw cycle is calculated. The daily temperature difference amplitudes of all freeze-thaw cycles are weighted and accumulated, with the weight being the change in soil moisture content before and after each freeze-thaw cycle, ultimately calculating the seasonal freeze-thaw cycle intensity factor. A time-series database of effective rainfall erosivity factors and seasonal freeze-thaw cycle intensity factors, bound to timestamps and geographic coordinates, is established. In practical implementation, the effective rainfall erosivity factor can be calculated using the following formula:

[0083]

[0084] in: Indicates the effective rainfall erosivity factor. This indicates the number of rainfall events during the year. Indicates the first The total kinetic energy of the rainfall, Indicates the first The maximum half-hour rainfall intensity of the event. It can be understood that the total kinetic energy in the formula... The calculation is based on the integration of rainfall intensity and duration, with the maximum half-hour rainfall intensity... Extract from a sliding window in a rainfall event.

[0085] In practice, the calculation of the seasonal freeze-thaw cycle intensity factor involves a weighted accumulation process. The weights for this accumulation are based on changes in soil moisture content, which are obtained from soil moisture observations or model estimations in long-term climate and hydrological time-series data. Comparing data from different years, for example, the number of freeze-thaw cycles may increase in colder years, leading to an increase in the seasonal freeze-thaw cycle intensity factor. Optionally, in the identification of freeze-thaw cycles, the freezing point temperature is set to 0 degrees Celsius, and the daily temperature range is the difference between the daily maximum and minimum temperatures. In some embodiments, the time-series database is constructed using a relational database, with each record associated with a timestamp, geographic coordinates, effective rainfall erosivity factor value, and seasonal freeze-thaw cycle intensity factor value, facilitating subsequent spatial analysis and model input.

[0086] Optionally, in the image segmentation algorithm, when separating soil particle regions, pore regions, and organic matter regions, a machine learning classifier can be combined to improve accuracy, such as using a support vector machine to classify image pixels. In specific implementations, the preset threshold for pore connectivity analysis is set according to the pore size distribution; for example, the connectivity threshold is set to a pore diameter greater than 10 micrometers. The preset angle for selecting vertical slices is set to 15 degrees to accurately reflect the density of vertical fracture development. Comparing image data at different soil depths, the density of vertical fracture development in surface soil is generally higher than that in deeper soils. It can be understood that the output format of the soil structure parameter thematic map is a raster image, with each pixel containing the soil aggregate destruction index and the vertical fracture development density value, and the spatial resolution is consistent with the original high-resolution soil structure tomography image.

[0087] In one embodiment of the present invention, see [reference] Figure 3 Hydrological analysis was performed on the digital model of surface micro-topography to extract slope, slope length, and topographic humidity index layers. Spatiotemporal overlay analysis was conducted on the land use change history vector map to calculate the number of land use type changes for each spatial unit within the assessment period, generating a land cover change frequency layer. From the land use change history vector map, patches related to cultivation activities were extracted, and the overlap of their interannual boundary movement trajectories and the consistency with cultivation direction were calculated to generate a cultivation disturbance accumulation layer. Principal component analysis was used to reduce the dimensionality and fuse the slope, slope length, topographic humidity index, land cover change frequency, and cultivation disturbance accumulation layers. The fused principal component values ​​were normalized and reclassified to generate a comprehensive surface disturbance intensity layer, whose pixel values ​​reflect the comprehensive intensity of the disturbance impact on surface stability.

[0088] In the specific implementation, a surface micro-topography digital model and a land use change history vector map are integrated to generate a surface disturbance intensity layer. This process is explained in detail through specific example scenarios and data comparisons. In one example scenario, the target area is a typical black soil hilly farmland. The surface micro-topography digital model is obtained by UAV LiDAR scanning with a spatial resolution of 0.5 meters. The land use change history vector map is derived from the annual remote sensing interpretation results of the past ten years, including polygonal patches of cultivated land, forest land, and grassland, along with their year attributes. By comparing data from different periods, the spatiotemporal changes in land cover and the trajectory of human activities can be quantified. Hydrological analysis is performed on the surface micro-topography digital model, using the D8 algorithm to calculate the direction of water flow, and then extracting slope, slope length, and topographic humidity index layers. Spatiotemporal overlay analysis is performed on the land use change history vector map. In the geographic information system software, the vector map layers from previous years are overlaid to calculate the number of land use type changes for each fixed-size spatial grid unit within the ten-year assessment period. The number of changes refers to the number of years in which the type has changed, generating a land cover change frequency layer. From the land use change history vector map, all cultivated land types are extracted as patches related to cultivation activities. The overlap of the interannual movement trajectories of cultivated land patch boundaries and the consistency with the cultivation direction are calculated. The cultivation direction is represented by the azimuth of the major axis of the patch, generating a cumulative cultivation disturbance layer. In some embodiments, the slope layer and slope length layer are calculated based on the elevation matrix of the surface micro-topography digital model, and the topographic moisture index represents the spatial distribution trend of soil moisture. It is understood that the size setting of the spatial grid cells needs to match the resolution of the surface micro-topography digital model and the subsequent model calculation scale, for example, set to 10 meters × 10 meters.

[0089] In practice, principal component analysis (PCA) is used to reduce the dimensionality and fuse the slope layer, slope length layer, topographic humidity index layer, land cover change frequency layer, and cultivated disturbance accumulation layer. PCA treats the five input layers as multiple variables and calculates their covariance matrix and eigenvectors. The fused principal component values ​​are then normalized and reclassified. Normalization linearly transforms the values ​​to between 0 and 1, while reclassification uses the natural discontinuity method to divide continuous values ​​into five discrete levels, generating a comprehensive surface disturbance intensity layer. The pixel values ​​of this layer reflect the comprehensive intensity of the disturbance impact on surface stability. Optionally, the first principal component is retained as the fusion result in PCA, as it typically contains the largest variance information from the original variable set. In practice, the calculation process of PCA can be described as finding a new orthogonal basis that maximizes the projection variance of the original multivariate data onto this basis. The calculation of the first principal component can be expressed as:

[0090]

[0091] in: This represents the score of the first principal component. This represents the standardized values ​​of the slope layer. The normalized values ​​represent the slope length layer. This represents the standardized value of the terrain humidity index layer. The standardized values ​​of the layer representing the frequency of land cover change. This represents the normalized value of the cumulative tillage disturbance layer. These are the eigenvector coefficients corresponding to the first principal component. It can be understood that the eigenvector coefficients are obtained by performing eigenvalue decomposition on the covariance matrix of the standardized data from the five input layers.

[0092] In some embodiments, the generation of the cumulative tillage disturbance layer involves calculating the overlap of interannual movement trajectories. The overlap is defined as the ratio of the area of ​​the spatial intersection to the area of ​​the union of tillage patches from adjacent years. The consistency of tillage direction is measured by calculating the standard deviation of the major axis azimuth of tillage patches over consecutive years; a smaller standard deviation indicates a more consistent tillage direction and a more stable disturbance pattern. In specific implementations, the values ​​of the five layers are standardized to eliminate the influence of dimensions. The standardization process uses the Z-score method. By comparing different sub-regions, such as the upslope region and the downslope region, the slope layer value of the upslope region is larger. If its land cover change frequency is also high, it may obtain a higher pixel value in the surface disturbance intensity layer. Optionally, in the spatiotemporal overlay analysis, the spatial grid cells can also adopt a size similar to the minimum patch area of ​​the land use change history vector map to reduce boundary effects.

[0093] In one embodiment of the present invention, the soil aggregate destruction index and vertical fracture development density are spatially superimposed, and a basic soil erodibility map is calculated using an empirical function. The effective rainfall erosivity factor and the seasonal freeze-thaw cycle intensity factor are weighted and fused over time to generate a dynamic climate erosion potential map. The surface disturbance intensity layer is multiplied pixel-by-pixel with the basic soil erodibility map to obtain a comprehensive erosion sensitivity map under the current state. In the coupled model of black soil erosion process, the comprehensive erosion sensitivity map is used as the initial field, and the dynamic climate erosion potential map is used as the driving field that changes over time, for iterative simulation based on physical processes. Within each iteration step, the amount of surface material migration and the increment of soil structure weakening are calculated, and the comprehensive erosion sensitivity map is dynamically updated. The simulation runs to a preset future time node, and the final total amount of soil structure weakening is compared with a preset threshold range to divide erosion and degradation risk zones of different levels, forming an erosion and degradation risk distribution map. The erosion and degradation risk distribution map expresses spatial differences in the form of a risk level raster.

[0094] A lookup table was established to correlate soil aggregate failure index with soil cohesion attenuation coefficient. A lookup table was also established to correlate vertical fracture development density with soil permeability correction coefficient. The soil aggregate failure index and vertical fracture development density for each spatial location were input into an empirical function for soil erodibility, which was expressed as a product of the cohesion attenuation coefficient, the permeability correction coefficient, and a standard soil erodibility factor. Soil erodibility values ​​were calculated for all spatial locations within the target black soil region. Spatial interpolation of the soil erodibility values ​​for all spatial locations was then performed to generate a continuous baseline map of soil erodibility.

[0095] After the simulation, the total soil structure weakening of each raster cell in the coupled model of black soil erosion process is extracted. Multiple threshold intervals characterizing the severity of erosion and degradation are set, each corresponding to a risk level. The total soil structure weakening of each raster cell is matched with the threshold intervals to determine the risk level of the raster cell. A unique identifier and color code are assigned to each risk level. All raster cells are rendered according to their risk level and color code to generate an erosion and degradation risk distribution map with spatial distribution information.

[0096] In the specific implementation, the calculation of the basic soil erosibility map, the generation of the dynamic map of climate erosion potential, the iterative simulation of the model, and the creation of the risk distribution map are explained in detail through specific example scenarios and data comparisons. In one example scenario, the target black soil region contains multiple geomorphic units, and the time series databases of effective rainfall erosivity factors and seasonal freeze-thaw cycle intensity factors have been prepared. The time span of the model simulation is set to the next 20 years, with an iteration step of 1 month.

[0097] In the specific implementation, the soil aggregate failure index and vertical fissure development density are spatially superimposed. This spatial superposition is performed using raster algebra operations under uniform geographic coordinates and spatial resolution. A basic soil erodibility map is obtained through empirical function calculations. The calculation process includes establishing a lookup table relationship between the soil aggregate failure index and the soil cohesion attenuation coefficient, based on regression analysis of indoor soil mechanics experiments and image parameters. A lookup table relationship is also established between the vertical fissure development density and the soil infiltration rate correction coefficient, based on regression analysis of infiltration test data and image parameters. The soil aggregate failure index and vertical fissure development density at each spatial location are input into the soil erodibility empirical function, which is expressed as a product of the cohesion attenuation coefficient, the infiltration rate correction coefficient, and the standard soil erodibility factor. Soil erodibility values ​​are calculated for all spatial locations within the target black soil region. Finally, spatial interpolation of all spatial location soil erodibility values ​​is performed using the Kriging method to generate a continuous basic soil erodibility map. It is understandable that the standard soil erodibility factor is a regional background constant derived from historical soil survey data.

[0098] In practical implementation, the empirical function of soil erodibility takes the following form:

[0099]

[0100] in: Indicates the first Soil erodibility values ​​at each spatial location This represents the first value obtained from the table of soil aggregate destruction index. Soil cohesion attenuation coefficient at each spatial location This represents the number obtained from looking up the table of vertical fracture development density. Soil infiltration rate correction factor for each spatial location This represents the standard soil erodibility factor. Comparing different plots, the soil cohesion attenuation coefficient is shown in areas with high soil aggregate destruction index. The value is relatively large, and the calculated value is... The value is also correspondingly higher.

[0101] In practical implementation, the effective rainfall erosivity factor and the seasonal freeze-thaw cycle intensity factor are fused using a time-series weighted fusion. This time-series weighted fusion involves weighted summation of the two factor values ​​for each pixel and each future simulated month, generating a dynamic climate erosion potential map, which is a raster sequence that changes over time. The surface disturbance intensity layer is multiplied pixel-by-pixel with the soil erodibility baseline map. This pixel-by-pixel multiplication is performed in a raster calculator to obtain the comprehensive erosion sensitivity map under the current state. In the coupled model of black soil erosion processes, the comprehensive erosion sensitivity map is used as the initial field, and the dynamic climate erosion potential map is used as the driving field that changes over time, for iterative simulation based on physical processes. Within each iteration step, the amount of surface material migration and the increment of soil structure weakening are calculated, and the comprehensive erosion sensitivity map is dynamically updated. When the simulation reaches a preset future time node, the final total amount of soil structure weakening is compared with a preset threshold range to delineate erosion and degradation risk zones of different levels, forming an erosion and degradation risk distribution map. This map expresses spatial differences in risk level raster format. Optionally, in time-series weighted fusion, the weights of the effective rainfall erosivity factor and the seasonal freeze-thaw cycle intensity factor can be determined based on statistical analysis of historical erosion events, for example, giving the seasonal freeze-thaw cycle intensity factor a higher weight during the spring snowmelt period.

[0102] In some embodiments, after the simulation, the total soil structure weakening of each raster cell output by the black soil erosion process coupling model is extracted. The total soil structure weakening is a dimensionless cumulative simulated value. Multiple threshold intervals characterizing the severity of erosion and degradation are set, each corresponding to a risk level. The total soil structure weakening of each raster cell is matched with the threshold intervals to determine the risk level to which the raster cell belongs. A unique identifier and color code are assigned to each risk level. All raster cells are rendered according to their risk level and color code to generate an erosion and degradation risk distribution map with spatial distribution information. In specific implementations, the threshold intervals can be set with reference to historical gully development rates or surface soil thickness loss data. See Table 1, which shows an exemplary threshold interval division relationship:

[0103] Table 1: Threshold Interval Division Relationship Table

[0104] Range of total soil structure weakening (S) Risk level Color codes (RGB) S<1.0 Low risk (0,255,0) 1.0≤S<2.5 Medium risk (255,255,0) 2.5≤S<5.0 High risk (255,165,0) S≥5.0 Extremely high risk (255,0,0)

[0105] It is understood that the total range of soil structure weakening values ​​in the table above are for illustrative purposes only, and actual applications require adjustments based on regional characteristics and model calibration results. In some embodiments, the logic for dynamically updating the integrated erosion sensitivity map is to use the increment of soil structure weakening calculated in each iteration as a decay coefficient, multiplied by the integrated erosion sensitivity map value of the previous iteration. Optionally, the risk level classification can be determined by analyzing the statistical distribution of the total soil structure weakening using the natural breakpoint method. By comparing simulation results under different future climate scenarios, such as high-emission and low-emission scenarios, the spatial patterns of the generated erosion and degradation risk distribution maps will show significant differences.

[0106] In one embodiment of the present invention, the method for constructing a coupled model of black soil erosion includes the following steps. Based on the principles of hydrodynamics and soil mechanics, a multi-process erosion dynamics equation set is established, coupling surface runoff erosion, gully erosion, and freeze-thaw stripping. This multi-process erosion dynamics equation set includes a runoff shear force sub-equation describing the sediment yield capacity of surface runoff, a gully headward erosion sub-equation describing the gully development rate, and a freeze-thaw stripping sub-equation describing the stripping of surface soil caused by freeze-thaw action. A numerical discretization method is used to perform spatiotemporal discretization of the multi-process erosion dynamics equation set, transforming the continuous equations into difference equations suitable for raster computation, and determining the iteration step size and boundary conditions of the model. A soil shear strength attenuation function is preset in the model, using the soil aggregate failure index, vertical fissure development density, and seasonal freeze-thaw cycle intensity factor as dynamic input parameters to update the soil's erosion resistance within each iteration step. A material migration and deposition function is preset in the model, using the micro-topography and surface cover reflected by the surface disturbance intensity layer as dynamic input parameters to simulate the transport and deposition process of eroded materials on the slope. Key parameters in the model were calibrated, including the Manning roughness coefficient in the runoff shear force sub-equation, the soil erodibility coefficient in the gully headward erosion sub-equation, and the latent heat coefficient of water phase change in the freeze-thaw stripping sub-equation. The calibrated model was then coupled and integrated with data structures for soil aggregate failure index, vertical fracture development density, effective rainfall erosivity factor, seasonal freeze-thaw cycle intensity factor, and surface disturbance intensity layer to complete the construction of a coupled model for black soil erosion.

[0107] In practical implementation, the construction of the coupled model of black soil erosion process is explained in detail through specific example scenarios and data comparisons. In one example scenario, the target black soil region is located in a cold region. The purpose of constructing the coupled model of black soil erosion process is to simulate the combined process of surface runoff erosion, gully erosion, and freeze-thaw stripping in this region. The input parameters required for model construction include soil mechanical parameters, hydrodynamic parameters, and thermodynamic parameters, which are obtained through literature review, indoor experiments, and field monitoring.

[0108] In practical implementation, based on the principles of hydrodynamics and soil mechanics, a multi-process erosion dynamics equation set coupling surface runoff erosion, gully erosion, and freeze-thaw stripping is established. This multi-process erosion dynamics equation set includes a runoff shear force sub-equation describing the sediment yield capacity of surface runoff, a gully headward erosion sub-equation describing the gully development rate, and a freeze-thaw stripping sub-equation describing the topsoil stripping caused by freeze-thaw action. A numerical discretization method is used to spatiotemporally discretize the multi-process erosion dynamics equation set. The numerical discretization method employs the finite difference method to transform the continuous equations into difference equations suitable for raster computation, and the iteration step size and boundary conditions of the model are determined. A soil shear strength attenuation function is preset in the model, using the soil aggregate failure index, vertical fracture development density, and seasonal freeze-thaw cycle intensity factor as dynamic input parameters to update the soil's erosion resistance within each iteration step. A material migration and deposition function is also preset in the model, using the micro-topography and surface cover reflected by the surface disturbance intensity layer as dynamic input parameters to simulate the transport and deposition process of eroded material on the slope. Key parameters in the model were calibrated, including the Manning roughness coefficient in the runoff shear force sub-equation, the soil erodibility coefficient in the gully headward erosion sub-equation, and the latent heat coefficient of water phase change in the freeze-thaw stripping sub-equation. The calibrated model was then coupled with data structures for soil aggregate failure index, vertical fracture development density, effective rainfall erosivity factor, seasonal freeze-thaw cycle intensity factor, and surface disturbance intensity layer to complete the construction of a coupled model for black soil erosion. It can be understood that the establishment of the multi-process erosion kinetic equation set referenced the modified general soil loss equation, the gully development model, and the freeze-thaw stripping physical model.

[0109] In some embodiments, the freeze-thaw stripping equation is used to quantify the loss of interparticle cohesion and surface stripping due to water phase change. An example of the freeze-thaw stripping equation is as follows:

[0110]

[0111] in: This represents the rate of soil stripping per unit area due to freeze-thaw cycles. This represents the latent heat of phase change of water. Indicates soil volumetric moisture content. Indicates time, Indicates soil temperature, This indicates the freezing point temperature of water. It is a unit step function, when the soil temperature Below freezing point The time function value is 1 if the time is constant, and 0 otherwise. The formula describes the relationship between soil moisture content changes, temperature conditions, and stripping rate. In practice, the parameter calibration process uses historical erosion data and adjusts key parameters to minimize the error between the model simulation results and the observed data. Table 2 shows the range of Manning roughness coefficient values ​​attempted in a single parameter calibration process and their corresponding land cover types.

[0112] Table 2: Range of Manning's Roughness Coefficient and Land Cover Types

[0113] Land cover type The range of possible values ​​for Manning's roughness coefficient (n) Final rate constant dense bare soil 0.010-0.020 0.015 Sparse vegetation 0.030-0.060 0.045 Dense grassland 0.150-0.240 0.180 Residual straw mulch 0.040-0.080 0.065

[0114] It is understandable that the final calibration values ​​in the table above are the results of optimization based on specific observation data of the target area, and need to be recalibrated when applied to other areas. Comparing different soil types, the soil erodibility coefficient of clayey black soil is usually lower than that of sandy soil, and different optimal parameter values ​​will be exhibited in the calibration.

[0115] In practical implementation, spatiotemporal discretization divides the computational domain into regular raster cells. The time iteration step size is determined based on rainfall data and the stability conditions of the freeze-thaw process; for example, a minute-level step size is used in the heavy rainfall-driven simulation phase, while an hour-level step size is used during the non-rainfall period. The soil shear strength decay function is expressed as a linear or nonlinear combination of the soil aggregate failure index, vertical fracture development density, and seasonal freeze-thaw cycle intensity factor, used to calculate the decay ratio of soil shear strength at each time step. The material migration and deposition function is calculated based on the water flow sediment transport capacity and topographic slope; deposition occurs when the sediment supply at a certain point on the slope exceeds the water flow sediment transport capacity. Optionally, the model boundary conditions can be set as follows: upstream is a flow boundary, downstream is a free outflow boundary, and both sides are stagnant boundaries. The coupling and docking process ensures that the spatial extent, projected coordinate system, and raster resolution of all input layers are completely consistent, and the model directly reads parameter values ​​from these spatial data layers according to the raster position for calculation. By comparing a simplified model that only considers runoff erosion with a coupled model of black soil erosion process, the latter can simulate a significant contribution of freeze-thaw stripping during the spring snowmelt period, thus more comprehensively reflecting the actual erosion process in the black soil region.

[0116] In one embodiment of the present invention, a risk mitigation strategy plan is generated based on an erosion and degradation risk distribution map. This process is described in detail through specific example scenarios and data comparisons. In one example scenario, an erosion and degradation risk distribution map of the target black soil area has been generated through model simulation. The map divides the risk into four levels: low, medium, high, and extremely high. The high-risk and extremely high-risk areas exhibit discontinuous patchy distributions. Differentiated governance strategy plans need to be formulated for these areas to compare the matching results of measures under different combinations of driving factors.

[0117] In practical implementation, spatial areas with high and extremely high risk levels on the erosion and degradation risk distribution map are identified and defined as priority remediation areas. This identification process is completed through raster reclassification and vectorization procedures using Geographic Information System (GIS) software. For priority remediation areas, the dominant erosion driving factors are extracted from multi-source monitoring data. These dominant erosion driving factors include dominant climate factors, dominant topographic factors, and dominant anthropogenic disturbance factors. The extraction method involves statistically analyzing the mode or mean of the values ​​of each driving factor layer within the polygon of each priority remediation area. Based on the combination type of the dominant erosion driving factors, corresponding sets of engineering measures, biological measures, and management measures are matched from the strategy knowledge base, which stores the mapping relationship between different combinations of driving factors and remediation measures. Based on the principles of spatial proximity and measure compatibility, the matched sets of engineering measures, biological measures, and management measures are optimized and combined to form customized remediation plan packages for each priority remediation area. All customized remediation plan packages for priority remediation areas are spatially correlated with the erosion and degradation risk distribution map to generate a risk mitigation strategy plan map containing spatial location, risk level, and specific measures. It is understandable that the principle of spatial proximity refers to recommending coordinated and consistent governance measures for priority governance areas that are geographically adjacent and have similar driving factors, while the principle of measure compatibility means that the various engineering, biological and management measures combined should not conflict with each other when implemented.

[0118] In some embodiments, the specific steps for matching the corresponding engineering measure set, biological measure set, and management measure set from the strategy knowledge base include: encoding the dominant climate factor, dominant topographic factor, and dominant anthropogenic disturbance factor to form a three-dimensional driving factor code, for example, the code "A1-B2-C3" represents a specific dominant climate, topography, and anthropogenic disturbance type, respectively. In the strategy knowledge base, records that perfectly match the three-dimensional driving factor code are queried. These records store recommended lists of engineering measures, biological measures, and management measures. If no perfectly matching record is found, a fuzzy matching algorithm is used to search for records in the knowledge base whose driving factor code similarity exceeds a preset threshold, and extract their associated measure lists. All extracted measure lists are deduplicated and merged to form a preliminary engineering measure set, biological measure set, and management measure set. Combined with the on-site constraints of the priority governance area, the preliminary measure set is subjected to feasibility screening to obtain the final matched engineering measure set, biological measure set, and management measure set. In specific implementations, the fuzzy matching algorithm is implemented by calculating the similarity between driving factor codes. The formula for calculating the similarity S can be expressed as:

[0119]

[0120] in: This represents the similarity score between the query code and the code of a record in the knowledge base. Indicates the first of the query codes Bit character, This indicates the first digit of the encoded record in the knowledge base. Bit character, The function returns 1 if the two are the same, otherwise it returns 0. It is to give the first The weights of the bit-driven factors. Only when the similarity score... Only when the threshold value exceeds a preset threshold (e.g., 0.8) is the knowledge base record considered a similar match. By comparing different priority governance areas, an area coded as "heavy rainfall-steep slope-over-cultivation" and an area coded as "freeze-thaw-long slope-vegetation destruction" will be associated with different sets of measure records in the knowledge base through a fuzzy matching algorithm.

[0121] In practical implementation, the set of engineering measures may include specific engineering methods such as terracing, constructing millstones, and building drainage ditches; the set of biological measures may include specific biological methods such as planting shrub hedges, promoting crop rotation, and constructing protective forest belts; and the set of management measures may include specific management methods such as implementing fallow rotation, controlling grazing intensity, and adjusting fertilization systems. Optional on-site constraints include slope, soil thickness, property ownership, and project cost budget. In feasibility screening, for example, engineering measures involving large-scale mechanical construction will be excluded in areas with slopes greater than 25 degrees. It can be understood that the strategy knowledge base is constructed in the form of a relational database, with each record containing a driving factor encoding field and associated engineering, biological, and management measure text description fields. The optimization and combination process of customized governance solution packages may merge identical measures matching in two adjacent priority governance areas to form a larger contiguous governance engineering unit, thereby improving governance efficiency and economies of scale. By comparing the final generated risk mitigation strategy plan map with the original erosion and degradation risk distribution map, the differentiated governance measure layout planned for different spatial risk areas can be intuitively seen.

[0122] The above embodiments are only used to illustrate the technical methods of the present invention and are not intended to limit it. Although the present invention has been described in detail with reference to preferred embodiments, those skilled in the art should understand that modifications or equivalent substitutions can be made to the technical methods of the present invention without departing from the spirit and scope of the technical methods of the present invention.

Claims

1. A method for assessing the erosion risk of black soil based on erosion and degradation, characterized in that, The method includes: Collect multi-source monitoring data of the target black soil area. The multi-source monitoring data includes high-resolution soil structure tomography images, long-term climate and hydrological time series data, digital model of surface micro-topography, and historical vector map of land use change. The high-resolution soil structure tomography image is subjected to pore connectivity analysis to extract the soil aggregate destruction index and vertical fissure development density. The soil aggregate destruction index represents the degree of structural degradation caused by tillage or erosion, and the vertical fissure development density represents the complexity of water infiltration channels. Based on the aforementioned long-term climate and hydrological time-series data, the effective rainfall erosivity factor and the seasonal freeze-thaw cycle intensity factor are calculated. By integrating the digital model of surface micro-topography with the historical vector map of land use change, a surface disturbance intensity layer is generated. The surface disturbance intensity layer comprehensively reflects the combined effects of slope and slope length, surface cover change frequency and cultivation activity trajectory on soil stability. The soil aggregate failure index, the vertical fracture development density, the effective rainfall erosivity factor, the seasonal freeze-thaw cycle intensity factor, and the surface disturbance intensity layer are input into a coupled model of black soil erosion process. Through spatial iterative calculations of the model, erosion and degradation risk distribution maps of the target black soil area under different future scenarios are generated, including: The soil aggregate destruction index and the vertical fracture development density are spatially superimposed, and the basic soil erodibility map is calculated using an empirical function. The effective rainfall erosivity factor and the seasonal freeze-thaw cycle intensity factor are weighted and fused over time to generate a dynamic map of climate erosion potential. The surface disturbance intensity layer is multiplied pixel by pixel with the soil erodibility baseline map to obtain the comprehensive erosion sensitivity map under the current state. In the coupled model of the black soil erosion process, the comprehensive erosion sensitivity map is used as the initial field, and the dynamic map of climate erosion potential is used as the driving field that changes over time, to carry out iterative simulation based on physical processes. Within each iteration step, the amount of surface material migration and the increment of soil structure weakening are calculated, and the comprehensive erosion sensitivity map is dynamically updated. The simulation runs until a preset future time point, and the final total amount of soil structure weakening is compared with a preset threshold range to divide erosion and degradation risk zones of different levels, forming the erosion and degradation risk distribution map. The erosion and degradation risk distribution map expresses spatial differences in the form of risk level grids.

2. The method for assessing the risk of black soil erosion based on erosion and degradation according to claim 1, characterized in that, The high-resolution soil structure tomography image was subjected to pore connectivity analysis to extract the soil aggregate disruption index and vertical fracture development density, including: A digital image segmentation algorithm is used to separate soil particles, pores, and organic matter regions from the high-resolution soil structure tomography image; The separated pore regions are skeletonized and topologically analyzed to calculate the connectivity path length and bottleneck width of the pore network and identify the preferred flow paths. Based on the geometric morphology and distribution uniformity of soil particle regions, the average equivalent diameter and shape factor dispersion of aggregates are calculated, and combined with historical benchmark images, the soil aggregate destruction index is quantified. On a vertical slice, the number of pores with connectivity reaching a preset threshold within a unit area is counted, and the angle between their orientation and the vertical direction is measured. Pores with an angle less than a preset angle are selected, and their number is the vertical fracture development density. The output is a thematic map of soil structure parameters, which includes pixel coordinates, the soil aggregate destruction index, and the corresponding relationship between the vertical fissure development density.

3. The method for assessing the risk of black soil erosion based on erosion and degradation according to claim 2, characterized in that, Based on the aforementioned long-term climate and hydrological time-series data, the effective rainfall erosivity factor and the seasonal freeze-thaw cycle intensity factor are calculated, including: Extract minute-by-minute rainfall intensity sequences and daily average temperature sequences from the long-term climate and hydrological time-series data; For the minute-by-minute rainfall intensity sequence, the rainfall events are divided into events, the total kinetic energy and the maximum half-hour rainfall intensity of each rainfall event are calculated, and the total kinetic energy and the maximum half-hour rainfall intensity are accumulated to obtain the effective rainfall erosivity factor on an annual scale; From the daily average temperature sequence, the days in which the daily maximum temperature is below freezing and the daily minimum temperature is above freezing are identified and defined as one freeze-thaw cycle. The total number of freeze-thaw cycles in a year is counted as the number of freeze-thaw days. For each freeze-thaw cycle, the daily temperature difference amplitude is calculated, and the daily temperature difference amplitudes of all freeze-thaw cycles are weighted and accumulated, where the weight is the change in soil moisture content before and after each freeze-thaw cycle, and finally the seasonal freeze-thaw cycle intensity factor is calculated. Establish a time-series database of the effective rainfall erosivity factor and the seasonal freeze-thaw cycle intensity factor, which are bound to timestamps and geographic coordinates.

4. The method for assessing the risk of black soil erosion based on erosion and degradation according to claim 3, characterized in that, By integrating the aforementioned digital model of surface micro-topography with the aforementioned land use change history vector map, a surface disturbance intensity layer is generated, including: Hydrological analysis was performed on the digital model of the surface micro-topography to extract slope layer, slope length layer and topographic humidity index layer; Spatiotemporal overlay analysis was performed on the land use change history vector map to calculate the number of land use type changes for each spatial unit within the assessment period, and a land cover change frequency layer was generated. From the land use change history vector map, extract the patches related to cultivation activities, calculate the overlap of their boundary interannual movement trajectories and the consistency with the cultivation direction, and generate a cultivation disturbance accumulation layer. Principal component analysis was used to reduce the dimensionality and fuse the slope layer, slope length layer, topographic humidity index layer, land cover change frequency layer, and cultivated disturbance accumulation layer. The fused principal component values ​​are normalized and reclassified to generate a comprehensive surface disturbance intensity layer, whose pixel values ​​reflect the comprehensive intensity of the disturbance affecting surface stability.

5. The method for assessing the risk of black soil erosion based on erosion and degradation according to claim 4, characterized in that, The soil aggregate failure index and the vertical fracture development density are spatially superimposed, and a basic soil erodibility map is calculated using an empirical function, including: Establish a lookup table relationship between the soil aggregate destruction index and the soil cohesion attenuation coefficient; Establish a lookup table relationship between the vertical fracture development density and the soil permeability rate correction coefficient; The soil aggregate destruction index and the vertical fracture development density at each spatial location are input into the soil erodibility empirical function, which is expressed as a product of the cohesion attenuation coefficient, the permeability correction coefficient, and the standard soil erodibility factor. Calculate the soil erodibility values ​​at all spatial locations within the target black soil region; Spatially interpolate the soil erodibility values ​​at all spatial locations to generate a continuous baseline map of soil erodibility.

6. The method for assessing the risk of black soil erosion based on erosion and degradation according to claim 5, characterized in that, The simulation runs until a preset future time point, then compares the final total amount of soil structure weakening with a preset threshold range to divide erosion and degradation risk zones of different levels, including: After the simulation, the total amount of soil structure weakening in each grid cell output by the coupled model of the black soil erosion process was extracted; Multiple threshold ranges are set to characterize the severity of erosion and degradation, with each threshold range corresponding to a risk level; The total amount of soil structure weakening in each grid cell is matched with the threshold range to determine the risk level of the grid cell. Assign a unique identifier and color code to each risk level; All grid cells are rendered according to their risk level and color code to generate the erosion and degradation risk distribution map with spatial distribution information.

7. The method for assessing the risk of black soil erosion based on erosion and degradation according to claim 6, characterized in that, The method for constructing the coupled model of the black soil erosion process includes: Based on the principles of hydrodynamics and soil mechanics, a multi-process erosion dynamics equation set is established that couples surface runoff erosion, gully erosion and freeze-thaw stripping. The multi-process erosion dynamics equation set includes a runoff shear force sub-equation describing the sediment yield of surface runoff, a gully headward erosion sub-equation describing the gully development rate, and a freeze-thaw stripping sub-equation describing the stripping of surface soil caused by freeze-thaw action. The multi-process erosion dynamics equations are spatiotemporally discretized using a numerical discretization method, transforming the continuous equations into difference equations suitable for grid computing, and determining the iteration step size and boundary conditions of the model. A soil shear strength decay function is preset in the model. The soil shear strength decay function uses the soil aggregate failure index, the vertical crack development density and the seasonal freeze-thaw cycle intensity factor as dynamic input parameters to update the soil's erosion resistance in each iteration step. A material migration and deposition function is preset in the model. The material migration and deposition function uses the micro-topography and surface cover reflected by the surface disturbance intensity layer as dynamic input parameters to simulate the transport and deposition process of eroded materials on the slope. The key parameters in the model are calibrated, including the Manning roughness coefficient in the runoff shear force sub-equation, the soil erodibility coefficient in the gully headward erosion sub-equation, and the latent heat coefficient of water phase change in the freeze-thaw stripping sub-equation. The calibrated model is coupled and integrated with the data structures of the soil aggregate destruction index, the vertical fissure development density, the effective rainfall erosivity factor, the seasonal freeze-thaw cycle intensity factor, and the surface disturbance intensity layer to complete the construction of the coupled model of the black soil erosion process.

8. The method for assessing the risk of black soil erosion based on erosion and degradation according to claim 7, characterized in that, The method further includes the step of generating a risk mitigation strategy plan based on the erosion and degradation risk distribution map: Spatial areas with high and extremely high risk levels identified in the erosion and degradation risk distribution map are defined as priority remediation areas; For the priority treatment area, the dominant erosion driving factors for the corresponding area are extracted from the multi-source monitoring data. The dominant erosion driving factors include dominant climate factors, dominant topographic factors and dominant anthropogenic disturbance factors. Based on the combination type of the dominant erosion driving factors, the corresponding sets of engineering measures, biological measures, and management measures are matched from the strategy knowledge base, which stores the mapping relationship between different combinations of driving factors and control measures. Based on the principles of spatial proximity and measure compatibility, the matched sets of engineering measures, biological measures, and management measures are optimized and combined to form a customized governance solution package for each priority governance area. The customized governance packages for all priority governance areas are spatially correlated with the aforementioned erosion and degradation risk distribution map to generate a risk mitigation strategy plan map that includes spatial location, risk level, and specific measures. The step of matching corresponding sets of engineering measures, biological measures, and management measures from the strategy knowledge base based on the combination type of the dominant erosion driving factors includes: The dominant climate factor, dominant topographic factor, and dominant anthropogenic disturbance factor are encoded to form a three-dimensional driving factor code; In the strategy knowledge base, a record is queried that exactly matches the three-dimensional driving factor code. The record stores a list of recommended engineering measures, a list of biological measures, and a list of management measures. If no perfect match is found, a fuzzy matching algorithm is used to search for records in the knowledge base whose driving factor coding similarity exceeds a preset threshold, and to extract the list of measures associated with them. The extracted lists of measures are deduplicated and merged to form preliminary sets of engineering measures, biological measures, and management measures. Based on the on-site constraints of the priority treatment area, the preliminary set of measures is screened for feasibility, resulting in the final matching set of engineering measures, biological measures, and management measures.

9. A risk assessment system for black soil erosion based on erosion and degradation, comprising a memory, a processor, and a computer program stored in the memory and running on the processor, characterized in that, When the processor executes the computer program, it implements the steps of the black soil erosion risk assessment method based on erosion and degradation as described in any one of claims 1 to 8.

Images