Soil erosion analysis method and system

Abstract

The invention provides a soil erosion analysis method and system, and belongs to the field of water and soil conservation monitoring. The method comprises the steps of inputting environmental data of a target area into a soil erosion calculation model to generate a current soil erosion grade map of the target area, and obtaining a soil erosion calculation result of the target area; inputting the environmental data of the target area into a soil erosion prediction model to predict a soil erosion grade map of the target area in a preset future time length to obtain soil erosion prediction data of the target area; and performing soil erosion analysis based on the soil erosion calculation result of the target area, the soil erosion prediction data and a preset database to obtain a soil erosion analysis result. Therefore, soil erosion dynamic monitoring is carried out on the cross-regional small watershed, automatic calculation and soil erosion change prediction are carried out on soil erosion of different regions, the purpose of analyzing the future change trend of the soil erosion of the different regions is achieved, and support is provided for generating water and soil conservation measure suggestions.

Description

Technical Field

[0001] The present invention relates to the field of soil and water conservation monitoring, and in particular to a soil erosion analysis method, a soil erosion analysis system, a machine-readable storage medium and an electronic device. Background Art

[0002] At present, soil erosion areas are mainly concentrated in the Loess Plateau, the hills in the middle and upper reaches of the Yangtze River, and the plains of Northeast China. Soil erosion is widely distributed, covers a large area, has many types, and has complex causes. In recent years, with the rise of artificial intelligence, existing technologies can realize the dynamic calculation of soil erosion in a single small watershed, but the calculations are mostly for fixed areas, and cross-regional calculations require a lot of time, manpower, and capital investment to rebuild the model. In addition, existing technologies require a large amount of long-term measured data. These on-site measurements are expensive and time-consuming, and cannot meet the needs of enterprises' medium- and short-term soil and water conservation work.

[0003] Enterprises need to invest a lot of money in soil and water conservation every year, and soil erosion, as an important indicator for dynamic monitoring of soil and water loss, can not only verify the effectiveness of the enterprise's soil and water conservation work but also provide support for future soil and water conservation work plans. Therefore, how to achieve dynamic monitoring of soil erosion in multiple cross-regional small watersheds, automatically calculate soil erosion in different regions and predict soil erosion changes, so as to analyze the future trend of soil erosion in different regions and provide support for the generation of soil and water conservation measures is an urgent problem to be solved. Summary of the invention

[0004] The purpose of the embodiments of the present invention is to provide a soil erosion analysis method and system to at least solve the above-mentioned problem of how to realize dynamic monitoring of soil erosion in multiple cross-regional small watersheds, so as to automatically calculate soil erosion in different regions and predict soil erosion changes.

[0005] In order to achieve the above-mentioned purpose, the first aspect of the present invention provides a soil erosion analysis method, comprising: inputting the preprocessed environmental data of the target area into a trained soil erosion calculation model to generate a current soil erosion grade map of the target area, and obtaining a soil erosion calculation result of the target area; inputting the preprocessed environmental data of the target area into a trained soil erosion prediction model to predict the soil erosion grade map of the target area within a preset future time length, and obtaining soil erosion prediction data of the target area; wherein the preset future time length represents a future time of a preset time length with the current moment as the time starting point; performing soil erosion analysis based on the soil erosion calculation result of the target area, the soil erosion prediction data and a preset database, and obtaining a soil erosion analysis result.

[0006] Optionally, the acquisition rules of the environmental data of the above-mentioned target area include: inputting the retrieval range data of the target area into a preset database for matching to obtain the environmental data of the target area; wherein the retrieval range data of the target area includes the geometric boundary parameters of the target area, data collection start and end date parameters, data resolution parameters, cloud cover threshold parameters and file format parameters; the preset database includes multiple data sources; the environmental data of the target area includes meteorological data, image data, hydrological data, soil data and land management data of the target area.

[0007] Optionally, the preprocessing rules for the environmental data of the target area include: preprocessing the meteorological data, image data, hydrological data, soil data and land management data of the target area respectively; wherein, the preprocessing rules for the meteorological data of the target area include: batch reading the folders corresponding to the meteorological data of the target area to construct a meteorological data list, and performing data cleaning on the meteorological data in the meteorological data list; the preprocessing rules for the image data of the target area include: when it is determined that the image data of the target area is complete, performing image processing on the image data of the target area to obtain image processing data of the target area; wherein, the image processing includes one or more of band stacking, geometric correction, atmospheric correction, image stitching and image cropping; performing projection conversion and resampling on the image processing data of the target area in turn to obtain resampled image data. The method comprises the following steps: first, performing a projection transformation and resampling on the soil data of the target area in turn to obtain the resampled soil data; second, performing a projection transformation and resampling on the soil data of the target area in turn to obtain the resampled soil data; third, generating and saving a soil data map based on the resampled soil data; fourth, performing a projection transformation and resampling on the land management data of the target area in turn to obtain the resampled land management data; and fourth, generating and saving a land management data map based on the resampled land management data.

[0008] Optionally, the training rules of the above-mentioned soil erosion calculation model include: obtaining multiple first preset object data and soil erosion data corresponding to each first preset object data in a preset database, and performing projection conversion, resampling, cropping and normalization processing on each first preset object data in turn to superimpose each first preset object data into a multi-channel raster file; performing data enhancement processing on the multi-channel raster file, generating a data sample set based on the soil erosion data corresponding to each first preset object data and the corresponding multi-channel raster file after data enhancement processing, and dividing the data sample set into a first training set and a first test set according to a preset ratio; based on the Bayesian optimization algorithm and a pre-set hyperparameter search space, performing hyperparameter optimization on the first training set and the first test set to obtain the optimal hyperparameter combination; based on the optimal hyperparameter combination, training the soil erosion calculation model to obtain a trained soil erosion calculation model; wherein the soil erosion calculation model is constructed based on a convolutional neural network algorithm.

[0009] Optionally, the training rules of the above-mentioned soil erosion prediction model include: obtaining multiple second preset object data and soil erosion data corresponding to each second preset object data in a preset database, and performing projection conversion, resampling, cropping, normalization and data enhancement processing on each second preset object data in turn; wherein the time length between the collection time of each second preset object data is a preset time length; based on each second preset object data after data enhancement processing and the soil erosion data corresponding to each second preset object data, a training sample set is generated; based on the training sample set, the soil erosion prediction model is trained to obtain a trained soil erosion prediction model; wherein the soil erosion prediction model is constructed based on a spatiotemporal convolutional neural network algorithm.

[0010] Optionally, the above soil erosion analysis includes spatial distribution analysis, change trend analysis, overlay analysis, erosion risk analysis and factor association analysis; wherein, the rules of spatial distribution analysis include: based on the soil erosion calculation results and soil erosion prediction data of the target area, drawing a hydraulic erosion spatial distribution map and a wind erosion spatial distribution map; comparing the hydraulic erosion spatial distribution map and the wind erosion spatial distribution map to determine the average erosion value, maximum erosion value, minimum erosion value, the location of the area reaching severe erosion conditions and the spatial distribution characteristics of the location of the area; the rules of change trend analysis include: based on the soil erosion calculation results, soil erosion prediction data and preset database of the target area, reading the soil erosion data of the target area in different years, and based on the soil erosion data of the target area in each year, counting the corresponding light erosion values ​​greater than the preset values ​​for each year. The erosion area of ​​the target area meets the intensity threshold; based on the erosion area corresponding to each year that is greater than the preset mild erosion intensity threshold, draw a trend map of the erosion area change in the target area; the rules of superposition analysis include: based on the spatial distribution map of hydraulic erosion and the spatial distribution map of wind erosion, mark the areas only affected by hydraulic erosion, the areas only affected by wind erosion, and the areas affected by both hydraulic erosion and wind erosion respectively; the rules of erosion risk analysis include: based on the soil erosion calculation results of the target area and the preset erosion risk level threshold, evaluate the erosion risk level of each area in the target area; the rules of factor association analysis include: based on the soil erosion calculation results of the target area, soil erosion prediction data, preset database and pre-determined influencing factors, use correlation analysis and / or regression analysis to determine the influencing factors that have the greatest impact on soil erosion.

[0011] Optionally, after obtaining the soil erosion analysis results, the method further includes: determining a soil and water conservation plan based on the soil erosion analysis results.

[0012] The second aspect of the present invention provides a soil erosion analysis system, comprising: a soil erosion calculation module, which is used to input the pre-processed environmental data of the target area into a trained soil erosion calculation model to generate a current soil erosion grade map of the target area and obtain a soil erosion calculation result of the target area; a soil erosion prediction module, which is used to input the pre-processed environmental data of the target area into the trained soil erosion prediction model to predict the soil erosion grade map of the target area within a preset future time length and obtain soil erosion prediction data of the target area; wherein the preset future time length represents a future time of a preset time length with the current moment as the time starting point; a soil erosion analysis module, which is used to perform soil erosion analysis based on the soil erosion calculation result of the target area, the soil erosion prediction data and a preset database to obtain a soil erosion analysis result.

[0013] In a third aspect of the present invention, a machine-readable storage medium is provided, on which instructions are stored. When the instructions are executed by a processor, the processor is configured to execute the soil erosion analysis method.

[0014] In a fourth aspect of the present invention, an electronic device is provided. The electronic device includes a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the computer program, the soil erosion analysis method is implemented.

[0015] Through the above technical scheme, a soil erosion analysis method and system are provided to input the pre-processed environmental data of the target area into the trained soil erosion calculation model to generate the current soil erosion grade map of the target area, thereby obtaining the soil erosion calculation result of the target area. The pre-processed environmental data of the target area is input into the trained soil erosion prediction model to predict and generate the soil erosion grade map of the target area within a preset future time length, thereby obtaining the soil erosion prediction data of the target area. Soil erosion analysis is performed in combination with the soil erosion calculation results, soil erosion prediction data and preset database of the target area to obtain the soil erosion analysis results. Thereby, the purpose of dynamically monitoring soil erosion in multiple cross-regional small watersheds is achieved, and soil erosion in different regions is automatically calculated and soil erosion changes are predicted, thereby achieving the purpose of analyzing the future change trend of soil erosion in different regions, so as to provide support for the generation of soil and water conservation measures.

[0016] Other features and advantages of the embodiments of the present invention will be described in detail in the following detailed description. BRIEF DESCRIPTION OF THE DRAWINGS

[0017] The accompanying drawings are used to provide a further understanding of the embodiments of the present invention and constitute a part of the specification. Together with the following specific embodiments, they are used to explain the embodiments of the present invention, but do not constitute a limitation on the embodiments of the present invention. In the accompanying drawings: Figure 1 is a flow chart of a soil erosion analysis method provided by one embodiment of the present invention; Figure 2 is a flow chart of obtaining environmental data of a target area provided by an embodiment of the present invention; Figure 3 is a flow chart of preprocessing environmental data of a target area provided by one embodiment of the present invention; Figure 4 It is a structural schematic diagram of a CNN model provided by one embodiment of the present invention; Figure 5 It is a training flow chart of a soil erosion calculation model provided by one embodiment of the present invention; Figure 6 is a schematic diagram of data acquisition and preprocessing provided by one embodiment of the present invention; Figure 7 is a schematic diagram of data enhancement provided by an embodiment of the present invention; Figure 8 is an accuracy curve diagram provided by an embodiment of the present invention; Fig. 9 is a loss value curve diagram provided by an embodiment of the present invention; Fig.10 It is a schematic diagram of hydraulic erosion calculated by a CNN model provided by one embodiment of the present invention; Fig.11 It is a residual graph of a model calculation result provided by an embodiment of the present invention; Fig.12 It is a schematic diagram of the structure of a STGCN model provided by an embodiment of the present invention; Fig.13 is another loss value curve graph provided by one embodiment of the present invention; Fig.14 It is a schematic diagram of hydraulic erosion prediction obtained by a STGCN model provided by one embodiment of the present invention; Fig.15 is an actual value and prediction residual graph provided by an embodiment of the present invention; Fig.16 It is a schematic diagram of a change in hydraulic erosion area provided by an embodiment of the present invention; Fig.17 It is a spatial distribution map of hydraulic and wind erosion provided by one embodiment of the present invention; Fig.18 is another spatial distribution map of hydraulic and wind erosion provided by one embodiment of the present invention; Fig.19 It is another spatial distribution map of hydraulic and wind erosion provided by one embodiment of the present invention; Fig. 20 is another spatial distribution map of hydraulic and wind erosion provided by an embodiment of the present invention; Fig.21 Another spatial distribution map of hydraulic and wind erosion provided by an embodiment of the present invention; Fig. 22 is another schematic diagram of an erosion area change curve provided by an embodiment of the present invention; Fig.23 It is a schematic diagram of erosion superposition provided by one embodiment of the present invention; Fig.24is a schematic diagram of erosion risk levels provided by an embodiment of the present invention; Fig.25 It is a block diagram of a soil erosion analysis system provided by one embodiment of the present invention. DETAILED DESCRIPTION

[0018] The specific implementation of the embodiment of the present invention is described in detail below in conjunction with the accompanying drawings. It should be understood that the specific implementation described here is only used to illustrate and explain the embodiment of the present invention, and is not used to limit the embodiment of the present invention.

[0019] It should be noted that, in the absence of conflict, the embodiments of the present invention and the features in the embodiments may be combined with each other.

[0020] The present invention will be described in detail below with reference to the accompanying drawings and in conjunction with embodiments.

[0021] Figure 1 FIG. 1 is a flow chart of a soil erosion analysis method provided by one embodiment of the present invention. Figure 1 As shown, an embodiment of the present invention provides a soil erosion analysis method, comprising: S110: inputting the preprocessed environmental data of the target area into the trained soil erosion calculation model to generate a current soil erosion grade map of the target area and obtain a soil erosion calculation result of the target area; It should be noted that the soil erosion grade map in this application belongs to the soil erosion type map. The target area can be multiple small watersheds across the region.

[0022] In some implementations of the present embodiment, the acquisition rules of the environmental data of the above-mentioned target area include: inputting the retrieval range data of the target area into a preset database for matching to obtain the environmental data of the target area; wherein the retrieval range data of the target area includes the geometric boundary parameters of the target area, the data collection start and end date parameters, the data resolution parameters, the cloud cover threshold parameters and the file format parameters; the preset database includes multiple data sources; the environmental data of the target area includes the meteorological data, image data, hydrological data, soil data and land management data of the target area.

[0023] Please refer to Figure 2 , Figure 2It is a flowchart of obtaining environmental data of a target area provided by an embodiment of the present invention. Specifically, the steps of collecting environmental data of the target area in the method include: first, obtaining the geometric boundary of the task area (it should be noted that the task area in this application is the target area) through a user interactive map tool or uploading a geographic file (such as Shapefile, GeoJSON, KML, etc.). Then, the data collection start and end date parameters, data resolution parameters, cloud threshold parameters and file format parameters are obtained through user input. According to the information collected above, it is matched with a preset database (the preset database includes a local database, a public data platform API, and an external data source) to obtain relevant data. The relevant data finally obtained include: daily meteorological data and time-sharing meteorological data of a meteorological observation station, multispectral images (including at least four bands of blue, green, red and near infrared), digital ground model (DEM), soil moisture content image, surface runoff vector, soil distribution vector, land use status vector, human disturbance plot vector, soil and water conservation measures vector, soil erodibility factor K value vector, etc.

[0024] In some implementations of the present embodiment, the preprocessing rules for the environmental data of the target area include: preprocessing the meteorological data, image data, hydrological data, soil data and land management data of the target area respectively; wherein, the preprocessing rules for the meteorological data of the target area include: batch reading the folders corresponding to the meteorological data of the target area to construct a meteorological data list, and performing data cleaning on the meteorological data in the meteorological data list; the preprocessing rules for the image data of the target area include: when it is determined that the image data of the target area is complete, performing image processing on the image data of the target area to obtain image processing data of the target area; wherein, the image processing includes one or more of band stacking, geometric correction, atmospheric correction, image stitching and image cropping; performing projection conversion and resampling on the image processing data of the target area in sequence , and obtain resampled image data; based on the data type of the image data of the target area, match the corresponding image processing scheme, and based on the resampled image data, execute the corresponding image processing scheme to obtain the preprocessed image data of the target area; the preprocessing rules of the soil data of the target area include: when it is determined that the soil data of the target area is complete, the soil data of the target area is projected and resampled in sequence to obtain the resampled soil data; based on the resampled soil data, a soil data map is generated and saved; the preprocessing rules of the land management data of the target area include: when it is determined that the land management data of the target area is complete, the land management data of the target area is projected and resampled in sequence to obtain the resampled land management data; based on the resampled land management data, a land management data map is generated and saved.

[0025] Please refer to Figure 3, Figure 3 This is a flow chart of preprocessing environmental data of a target area provided by an embodiment of the present invention. The method preprocesses the environmental data of the target area, including: preprocessing meteorological data, multispectral image data, hydrological data, soil data and land management data accordingly. Specifically, the meteorological data, multispectral image data, hydrological data, soil data and land management data are processed as follows: ① Preprocessing of meteorological data: Meteorological data is in TXT, CSV, EXCEL, JSON and other formats. Meteorological data types are mainly divided into daily meteorological data and time-sharing meteorological data. First, use the Pandas library in Python to batch read the meteorological data folder to form a list, then clean the data, and finally save the historical rainfall table and historical wind speed table according to the meteorological observation station.

[0026] Daily meteorological data processing takes the historical daily rainfall data SURF_CLI_CHN_MUL_DAY_PRE of China Meteorological Data Network as an example: SURF_CLI_CHN_MUL_DAY_PRE is a collection of monthly rainfall records of national stations, the data format is TXT, and the data structure is district station number, latitude (unit: degree, minute), longitude (unit: degree, minute), observation field altitude (unit: m), year, month, day, 20-8:00 precipitation (night) (unit: mm), 8-20:00 precipitation (daytime) (unit: mm), 20-20:00 cumulative precipitation (unit: mm), 20-8:00 precipitation quality control code, 8-20:00 cumulative precipitation quality control code, 20-20:00 precipitation quality control code. 1. Filter files according to the start and end dates entered by the user, and select the record files in the specified time range. 2. Batch read SURF_CLI_CHN_MUL_DAY_PRE files to build a list. When reading, filter by the station number of the mission area observation station. First check the data integrity and throw an exception prompt for missing data. In addition, combine the three columns of year, month, and day into a new column "date" in the standard timestamp format: YYYY-MM-DD:HH-MM-SS. 3. Select the station number, latitude (unit: degree, minute), longitude (unit: degree, minute), date (YYYY-MM-DD:HH-MM-SS), and cumulative precipitation from 20:00 to 20:00 (unit: mm) in the list, and create a new table. In the new table, replace the column name "cumulative precipitation from 20:00 to 20:00" with "cumulative rainfall". 4. Group the list according to the station number and split the list. 5. De-duplicate the list data according to the date column and sort it in ascending order. Check whether the date is continuous and not missing. If missing, throw an exception prompt. 6. Convert the characteristic value codes (32766, 32700, 32XXX, 31XXX, 30XXX) in the cumulative rainfall column, convert missing values, null values ​​and abnormal values ​​to zero, and correct the cumulative rainfall value according to the scale factor. 7. Save the processed list data as a CSV file according to the district station number.

[0027] The processing of time-sharing meteorological data takes the INTEGRATED SURFACEDATA data of NOAA (National Oceanic and Atmospheric Administration) as an example: INTEGRATED SURFACE DATA is a collection of meteorological data of meteorological observation stations throughout the year. The data format is CSV, and the data structure is STATION meteorological observation station area station number, DATE date (format: YYYY-MM-DD:HH-MM-SS), LATITUDE latitude (unit: degree), LONGITUDE longitude (unit: degree), ELEVATION elevation (unit: meter), NAME observation station name, REPORT_TYPE report format, CALL_SIGN call sign code, QUALITY_CONTROL quality control code, WND wind speed data group (content: wind observation direction angle, wind observation direction quality code, wind force observation type code, wind speed, wind speed observation quality code), etc. 1. Filter files according to the start and end dates entered by the user, and select the record files in the specified time range. 2. Batch read files to build a list. When reading, filter by the station number of the mission area observation station. Check the data integrity first, and throw an exception prompt for missing data. 3. Split the new column "WIND SPEED" wind speed from the WND column, select the STATION meteorological observation station number, DATE date (format: YYYY-MM-DD:HH-MM-SS), LATITUDE latitude (unit: degree), LONGITUDE longitude (unit: degree), WIND SPEED wind speed (unit: meter / second) to create a new table. 4. The data interval in the original data is 3 hours. Use the DATE date column to deduplicate the list data and arrange it in ascending order. Check whether the date is continuous and not missing. If missing, throw an exception prompt. 5. Remove the "099999" suffix in the STATION column, replace the "9999" data in the WIND SPEED column with zero, and correct the wind speed value according to the scaling factor. 6. Insert new rows at 1-hour intervals, temporarily define the wind speed values ​​in the new rows as null values, and sort them in ascending order using the DATE column. Use linear interpolation to fill in all null wind speed values, and save the processed list data as a CSV file by district station number.

[0028] ② Preprocessing of image data: Image data is in the formats of GeoTIFF, HDF, NetCDF, JPEG2000, ENVI, etc. The main types are multispectral images, DEM data, and soil moisture content images.

[0029] Multispectral image data preprocessing includes: checking the integrity of each image data and band stacking; using ground control points to perform geometric correction on the image; performing atmospheric correction on the image to reduce or eliminate atmospheric interference on the image to obtain a surface reflectance image; checking and removing cloud cover; performing projection conversion according to the coordinate system input by the user; resampling according to the resolution input by the user; calculating the Normalized Difference Vegetation Index (NDVI) for the processed image; calculating the vegetation cover (FVC, Fractional Vegetation Cover) from NDVI through an empirical formula, and saving the FVC map (the horizontal and vertical coordinates are the geographic spatial position, and the pixel value is the vegetation cover).

[0030] The processing of multispectral image data takes the HLS Landsat Operational Land Imager Surface Reflectance and TOA Brightness DailyGlobal land remote sensing satellite Sentinel-2 data of NASA Earthdata (NASA Earth Database Project) as an example (because the HLS.30L data has been geometrically corrected and atmospherically corrected, these two steps are skipped): 1. Stack the 13 band rasters such as B01, B02, B03, B04, and B05. Check the integrity of each image data, record the images with missing data, and output the missing data list after processing. 2. The image after band stacking uses the threshold segmentation method to detect the cloud area, and then uses the mean interpolation method to fill the cloud area. 3. Perform projection conversion according to the coordinate system input by the user, resample according to the resolution input by the user, and save the image processing results. 4. Calculate the Normalized Difference Vegetation Index (NDVI): = In the formula, NDVI is the Normalized Difference Vegetation Index; NIR is the reflectance of the near infrared band; R is the reflectance of the visible red band. 5. Calculate the Fractional Vegetation Cover (FVC): FVC = , where NDVIsoil is the NDVI of areas with bare soil or no vegetation, and NDVIveg represents the NDVI value of pixels completely covered by vegetation, that is, the value of pure vegetation pixels. Save FVC data (horizontal and vertical coordinates are geographic spatial positions, units: ° or m, pixel values ​​are vegetation coverage, units: %).

[0031] DEM data preprocessing includes: checking image data integrity and image stitching; performing projection conversion according to the coordinate system input by the user; resampling according to the resolution input by the user; filling depressions and leveling, and saving the processed DEM data; using the elevation change rate to calculate the slope and save the slope map (the horizontal and vertical coordinates are the geographic spatial position, unit: degree or meter, and the pixel value is the slope value, unit: degree); calculating the flow direction by the D8 flow method (Eight-Direction Flow Model, by determining that each grid cell flows to the cell with the lowest elevation among its eight adjacent cells, and defining codes for the eight adjacent cells to indicate the flow direction.), D-Infinity flow method (D∞, by calculating the slope between each grid cell and its adjacent cells, and determining the angle of the water flow along the maximum descent direction of the slope.) or MFD flow method (a flow algorithm that distributes water flow from a grid cell to multiple adjacent cells based on slope and weight.) Direction); Through the reverse flow model, the upslope length is calculated by recursively accumulating the distance from each grid cell along the direction of maximum slope rise, and the upslope length map is saved (the horizontal and vertical coordinates are the geographic space position, unit: ° or m, the pixel value is the upslope length value, unit: m).

[0032] DEM data processing takes NASADEMMerged DEM Global 1 arc second V001 data processing of NASA Earthdata (National Aeronautics and Space Administration Earth Database Project) as an example: 1. Use the Gdal library to batch read HGT data, calculate the range of the raster file after splicing, check the integrity of the image data, and throw an exception prompt when the mission area data is missing. 2. Perform projection conversion according to the coordinate system input by the user, resample according to the resolution input by the user, and save the image processing results. 3. Use the Pysheds library to fill and flatten the image, and save the processed DEM data. 4. Use the difference method to calculate the elevation change rate of each grid cell in the x and y directions. The slope calculation formula is as follows: S = in, and are the rates of change of elevation in the x-direction and the y-direction respectively. The slope is expressed in degrees as follows: Slope(degrees) =arctan(S) , save the slope map. 5. Use the slope map to calculate the flow direction using the D8 flow method. The direction of the water flow is determined by comparing the elevation difference between each grid cell and its eight adjacent cells. This method assumes that water flows from cells with higher elevations to adjacent cells with the lowest elevations. The flow direction is represented by the codes of eight possible directions: east (1), southeast (2), south (4), southwest (8), west (16), northwest (32), north (64), and northeast (128). Save the flow map. 6. Calculate the upslope length using the flow map, and calculate the reverse flow direction based on the flow direction data. That is, reverse the direction of the water flow and determine that the water flow in each cell flows from downstream to upstream. Initialize the upslope length value of each grid cell to 0. Starting from the downstream cell, gradually accumulate the upslope length value of each cell along the reverse flow direction until you reach the ridge or watershed divide. Save the upslope length map.

[0033] Soil moisture content image preprocessing includes: checking image data integrity and image cropping; performing projection conversion according to the coordinate system input by the user; resampling according to the resolution input by the user; data cleaning, and saving soil moisture content maps (horizontal and vertical coordinates are geographic spatial positions, units: ° or m, pixel values ​​are soil moisture content, units: cm³ / cm³).

[0034] The soil moisture content image processing takes the SMAP Enhanced L3 Radiometer Global and Polar Grid Daily 9 km EASE-Grid SoilMoisture V006 data processing of NASA Earthdata (National Aeronautics and Space Administration Earth Database Project) as an example: 1. Use the Gdal library to batch read HGT data, calculate the range of the raster file after splicing, check the integrity of the image data, and throw an exception prompt when the mission area data is missing. 2. Perform projection conversion according to the coordinate system input by the user, resample according to the resolution input by the user, and save the image processing results. 3. Clean the image data, process outliers and null values, and perform replacement, interpolation or elimination. 4. Save the processed soil moisture content map. The saved image uses the geographic spatial position (° or m) as the horizontal and vertical coordinates, and the pixel value represents the soil moisture content (unit: cm³ / cm³) ③ Preprocessing of surface runoff data: Surface runoff data is in CSV, Excel, Shapefile, GeoTIFF, NetCDF and other formats. The main types are time series flow data, watershed raster data, and point source runoff data. The preprocessing process of surface runoff data includes: checking the integrity of each period of runoff data, processing missing values ​​and outliers; using interpolation methods to convert point source runoff data into continuous raster data; performing projection conversion according to the user input coordinate system; resampling according to the user input resolution; saving the processed runoff raster data (horizontal and vertical coordinates are geographic spatial locations, units: degrees or meters, pixel values ​​are runoff in the area, units: m 3 / s). The surface runoff data processing takes the basin point source runoff data provided by USGS as an example (the data format is a CSV file): 1. Check the integrity of each period of runoff data, handle abnormal values ​​and null values, perform replacement, interpolation or elimination, and record the list of missing values ​​after processing. 2. Perform projection conversion according to the user input coordinate system, resample according to the user input resolution, and save the image processing results. 3. Save the processed runoff raster data as a GeoTIFF file, and each pixel value represents the runoff volume of the area.

[0035] ④ Preprocessing of soil data: Soil data mainly includes soil texture, soil thickness, soil permeability, soil density, soil erodibility factor, soil erosion observation data, etc. The data file formats are CSV, Excel, Shapefile, GeoTIFF, etc. Soil data preprocessing includes: checking the integrity of each period of soil data, handling abnormal values ​​and null values, replacing, interpolating or removing them; performing projection conversion according to the coordinate system input by the user; resampling according to the resolution input by the user; generating and saving soil data maps (horizontal and vertical coordinates are geographic spatial positions, and pixel values ​​are soil attribute values). Soil texture data preprocessing takes the data processing of FAO (Food and Agriculture Organization of the United Nations) World Soil Database 2.0 (HWSD v2.0) as an example: 1. Check the integrity of the three fields 'SAND', 'SILT', and 'CLAY' in the database, handle abnormal values ​​and null values, replace, interpolate or remove them, and record the missing value rows after processing. 2. Extract the three fields 'SAND', 'SILT', and 'CLAY' from the database to form a multi-band raster file, perform projection conversion according to the coordinate system input by the user, resample according to the resolution input by the user, and save the raster processing results. 3. Generate a soil texture distribution map: a three-band raster map of sand, silt, and clay content, and save the soil texture distribution map data (the horizontal and vertical coordinates are the geographic spatial position, unit: ° or m, the pixel value is the sand, silt, and clay content of the area, unit: %).

[0036] ⑤ Preprocessing of land management data: Land management data mainly includes land use status, human disturbance, soil and water conservation measures, farming measures, etc. The data file formats are GeoJSON, Shapefile, GeoTIFF, etc. Land management data preprocessing includes: checking the integrity of land management data, processing outliers and null values, replacing or removing them, etc.; performing projection conversion according to the coordinate system input by the user; resampling according to the resolution input by the user; generating and saving land management data maps (horizontal and vertical coordinates are geographic spatial positions, and pixel values ​​represent the corresponding land management types.). The preprocessing of land use status data takes the processing of the third national land survey data as an example: 1. Check the integrity of the two fields 'DLBM' and 'DLMC' in the database, process outliers and null values, replace and remove them, and record the missing value rows after processing. 2. Extract the two fields 'DLBM' and 'DLMC' in the database, assign values ​​to different land use types to form a single-band raster file, perform projection conversion according to the coordinate system input by the user, resample according to the resolution input by the user, and save the raster processing results. Save the land use status map data (the horizontal and vertical coordinates are the geographic space position, unit: ° or m, each pixel value represents the corresponding land use type).

[0037] In some implementations of the present embodiment, the training rules of the above-mentioned soil erosion calculation model include: obtaining multiple first preset object data and soil erosion data corresponding to each first preset object data in a preset database, and performing projection conversion, resampling, cropping and normalization processing on each first preset object data in turn to superimpose each first preset object data into a multi-channel raster file; performing data enhancement processing on the multi-channel raster file, generating a data sample set based on the soil erosion data corresponding to each first preset object data and the corresponding multi-channel raster file after data enhancement processing, and dividing the data sample set into a first training set and a first test set according to a preset ratio; based on the Bayesian optimization algorithm and a pre-set hyperparameter search space, performing hyperparameter optimization on the first training set and the first test set to obtain the optimal hyperparameter combination; based on the optimal hyperparameter combination, training the soil erosion calculation model to obtain a trained soil erosion calculation model; wherein the soil erosion calculation model is constructed based on a convolutional neural network algorithm.

[0038] In some implementations of this embodiment, the soil erosion calculation model includes a hydraulic erosion model and a wind erosion model (wherein the hydraulic erosion model and the wind erosion model can be respectively the hydraulic erosion model and the wind erosion model in the "Technical Guide for Dynamic Monitoring of Soil and Water Loss in 2023" compiled by the Soil and Water Conservation Monitoring Center of the Ministry of Water Resources). Wherein: The hydraulic erosion model is: ; Where A represents the soil erosion modulus, the unit is ; R represents the rainfall erosivity factor, the unit is ; K represents soil erodibility factor, unit is ; L represents the slope length factor, which is a dimensionless parameter; S represents the slope factor, which is a dimensionless parameter; B represents the vegetation cover and biological measures factor, which is a dimensionless parameter; E represents the engineering measures factor, which is a dimensionless parameter; T represents the tillage measures factor, which is a dimensionless parameter.

[0039] ① Rainfall erosivity factor R: The calculation formula of rainfall erosivity is as follows: ; ; ;in, Indicates the average annual rainfall erosivity over many years, in units of MJ·mm·hm -2 ·h -1 ·a -1 ; k represents 1, 2, ..., 24, which means that a year is divided into 24 and a half months; Half month k represents the rainfall erosivity in the kth half month, in MJ·mm·hm -2 ·h -1 ; i represents 1, 2, ... N; N represents the time series input by the user, such as 1980-2024; j represents 0, 1, ... m; m represents the number of erosive rainfall days in the kth half month of the i-th year (erosive rainfall days refer to daily rainfall greater than or equal to 10 mm); Pi, j, k represents the jth erosive rainfall in the kth half month of the i-th year, in mm; if there is no erosive rainfall in a certain half month of a certain year, that is, j=0, then Pi, 0, k=0; α represents a parameter, α in the warm season (5 represents September) is 0.3937, and α in the cold season (10 represents December, 1 represents April) is 0.3101; Half-month k represents the average rainfall erosivity of the kth half-month ( Half-month k) accounts for the average annual rainfall erosivity over many years ( ) ratio.

[0040] The calculation steps of rainfall erosivity are as follows: 1. Call the pre-processed daily rainfall data in the local database. 1. The multi-year average semi-monthly erosivity is calculated using the semi-monthly k formula. 2. The multi-year average 1 to 24 semi-monthly rainfall erosivity is interpolated using the ordinary kriging spatial interpolation method to generate 24 semi-monthly rainfall erosivity raster data. The coordinate system and resolution are set according to user input. Save the 24 semi-monthly rainfall erosivity raster data as GeoTIFF (the horizontal and vertical coordinates are the geographic spatial position in degrees or meters, and the value of each pixel is the multi-year average semi-monthly erosivity in units of MJ·mm·hm -2 ·h -13. Accumulate the 24 and a half months of rainfall erosivity grid data into annual rainfall erosivity grid data (horizontal and vertical coordinates are geographic spatial positions, units are ° or m, and each pixel value is the annual average rainfall, units are MJ·mm·hm -2 ·h -1 ·a -1 ). 4. Divide the 24 and a half months of rainfall erosivity grid data by the annual rainfall erosivity grid data to obtain the grid data of the proportion of the 24 and a half months of rainfall erosivity to the annual rainfall erosivity (the horizontal and vertical coordinates are the geographic space position, the unit is ° or m, and the value of each pixel is the proportion of the half-month average rainfall erosivity, unit: %).

[0041] ②Soil erodibility factor K: Call the pre-processed soil erodibility factor grid in the local database. Perform projection conversion according to the user input coordinate system, resample according to the user input resolution, and save the grid processing results (horizontal and vertical coordinates are geographic spatial positions, units are ° or m, and each pixel value is the soil erodibility factor value of the corresponding area, units are t·h·MJ -1 mm -1 ).

[0042] ③ Slope length factor L and slope factor S: Slope length factor calculation formula: ; where λi and λi-1 represent the uphill length of the i-th and i-1-th slope sections respectively, in meters; m represents the slope length index, which varies with the slope and is a dimensionless parameter. The slope value is θ, in degrees, The slope factor calculation formula is: ; Where S represents the slope factor, which is a dimensionless parameter; θ represents the slope, in units of º. Call the preprocessed upslope length map and slope map in the local database. Perform matrix calculations according to the slope length factor formula and slope factor to obtain the slope length factor grid (horizontal and vertical coordinates are geographic spatial positions, in units of ° or m, and each pixel value is the slope length factor value of the corresponding area, which is a dimensionless parameter) and the slope factor grid (horizontal and vertical coordinates are geographic spatial positions, in units of ° or m, and each pixel value is the slope factor value of the corresponding area, which is a dimensionless parameter).

[0043] ④ Vegetation coverage and biological measure factor B: Vegetation coverage and biological measure factor B is calculated based on land use type and vegetation coverage, and the calculation results are combined to generate the B factor raster layer. Gardens, woodlands and grasslands are calculated using formulas, and other land use types are directly assigned by looking up Table 1. Calculation formula for garden, woodland and grassland B factors: ; Among them, WR i It indicates the ratio of the calculated rainfall erosivity of the i-th half month to the annual erosivity, with a value range of 0-1; SLR iIt represents the soil loss ratio of garden, woodland and grassland in the i-th half month. It is a dimensionless parameter with a value range of 0-1. i The calculation formula is: ; Orchards, other gardens, wooded land and other wooded land SLR i The calculation formula is: ; Grassland SLR i The calculation formula is: ; Among them, FVC represents the vegetation coverage calculated based on NDVI, with a value range of 0-1; GD represents the understory coverage of arbor forests, with a value range of 0-1, including the understory coverage composed of all vegetation (shrubs, herbs and litter) except the canopy of arbor forests. The understory coverage is assigned based on the positioning observation results or field surveys. As shown in Table 1, Table 1 is the B factor assignment table for non-garden land, woodland, and grassland: Table 1 B factor assignment table for non-garden land, forest land and grassland

[0044] The calculation steps of vegetation coverage and biological measure factor B are as follows: 1. Call the long-term vegetation coverage grid pre-processed in the local database, first group all vegetation coverage data by year according to the date, and then divide it into 24 half-months, and then calculate the half-month vegetation coverage mean, and save the half-month vegetation coverage mean grid of each year. 2. Call the land use grid data and half-month rainfall erosivity ratio grid data pre-processed in the local database, use the land use grid as a mask to extract the half-month vegetation coverage mean matrix, extract the vegetation coverage matrix of tea gardens and shrub forests, and calculate the vegetation coverage matrix according to the year. Calculate the B factor value, and process the orchard, other gardens, wooded land and other wooded land matrices, and grasslands in the same way. Finally, merge the raster matrices and save the B factor raster data of gardens, woodlands, and grasslands. 3. Replace the values ​​of non-gardens, woodlands, and grasslands in the land use raster matrix with the values ​​in Table 1, and replace the values ​​of other land classes with 0, and save the B factor raster data of non-gardens, woodlands, and grasslands. 4. Merge the B factor raster of gardens, woodlands, and grasslands with the B factor raster of non-gardens, woodlands, and grasslands to obtain the B factor raster of semi-monthly vegetation coverage and biological measure factors (the horizontal and vertical coordinates are the geographic spatial position, units: ° or m, and the value of each pixel is the B factor value of the semi-monthly vegetation coverage and biological measure factors, which is a dimensionless parameter). 5. Accumulate the semi-monthly B factor raster to obtain the B factor raster of a certain year, and save the annual B factor raster (the horizontal and vertical coordinates are the geographic spatial position, units: ° or m, and the value of each pixel is the B factor value of the annual vegetation coverage and biological measure factors, which is a dimensionless parameter).

[0045] ⑤ Soil and water conservation engineering measure factor E: Soil and water conservation engineering measure factor E is obtained by looking up the soil and water conservation engineering measure factor assignment table according to the engineering measure type of the soil and water conservation measure raster data, as shown in Table 2, to obtain the soil and water conservation engineering measure factor value.

[0046] Table 2 Soil and water conservation engineering measures factor assignment table

[0047] Call the preprocessed soil and water conservation measures map in the local database. Extract the soil and water conservation measures grid matrix, replace the engineering measures types in Table 1 with the E factor value, and replace the measures types not in Table 1 with 1, and save the soil and water conservation engineering measures factor E factor grid (horizontal and vertical coordinates are geographic spatial positions, units: ° or m, each pixel value is the soil and water conservation engineering measures factor E, which is a dimensionless parameter).

[0048] ⑥ Tillage measure factor T: Tillage measure factor T obtains the tillage measure factor value according to the tillage measure map and the tillage rotation measure assignment table, replaces the pixel value with the T factor value, and saves the T factor raster data (the horizontal and vertical coordinates are the geographic space position, unit: ° or m, each pixel value is the tillage measure factor T, which is a dimensionless parameter).

[0049] ⑦ Hydraulic erosion calculation: According to the hydraulic erosion model A=RKLSBET, the grids of the factors R, K, L, S, B, E, and T are multiplied to obtain the grid data of the hydraulic erosion modulus A. The hydraulic erosion modulus grid is saved (the horizontal and vertical coordinates are the geographic space position, unit: ° or m, and the value of each pixel is the hydraulic erosion modulus, unit: t·hm -2 ·a -1 ). Wind erosion can be calculated based on the land use type by selecting the corresponding sandy land (desert), grass (shrub) land, and cultivated land wind erosion models. The applicable scope of the wind erosion model is shown in Table 3: Table 3 Applicable scope of wind erosion model

[0050] The wind erosion model of cultivated land is: ; Among them, Q fa is the wind erosion modulus of cultivated land in every half month, t·hm -2 ·a -1 ; Z0 is the surface roughness, cm; the meanings of other parameters are the same as above. The wind erosion model for grass (shrub) land is: ; Among them, Q fg is the wind erosion modulus of sandy land (desert) every half month, t·hm -2 ·a -1 ; The meanings of other parameters are the same as above. The sandy land (desert) wind erosion model is: ; Among them, Qfs The wind erosion modulus of sandy land (desert) every half month, in t·hm -2 ·a -1 ; W is the surface soil moisture factor within every half month, between 0 and 1; T j is the cumulative time of each wind speed level in every half month, in min; V is the vegetation coverage, in %; U j is the average wind speed of the jth level, in m / s. For example, if the wind speed level is 5~6m / s, U j =5.5m / s.

[0051] ⑧ Wind factor T: The cumulative time corresponding to each wind speed level within every half month (T j ) is calculated as follows: ; Among them, T j It represents the cumulative time corresponding to the jth wind speed level in each meteorological station every half month, in min. For example, if the wind speed at the hour belongs to a certain wind speed level, the time t corresponding to the wind speed level is accumulated. jmik , otherwise it will not be accumulated; t jmik Indicates whether the wind speed at the kth moment on the i-th day in a half month of the m-th year at each meteorological station belongs to the j-th wind speed level. If yes, t jmik =1, otherwise t jmik =0; j represents the wind speed level number, which is divided into 35 levels at intervals of 1m / s within 5~40m / s, with values ​​of 1, 2, ..., 35; N represents the number of years for wind speed data collection; m takes values ​​of 1, 2, ..., N; L represents the number of days corresponding to each half month, with 15 days in the first half of each month and the rest in the second half of the month (13, 14, 15 or 16 days); i takes values ​​of 1, 2, ..., L; k represents a value in 24 hours of a day, with values ​​of 1, 2, ..., 24. The calculation time range of wind force factor is from January to December of the whole year. Calculate the cumulative time of each wind speed level. If the wind speed collected is at the hour of 24 o'clock every day, it is directly calculated based on 60 minutes per hour. If the wind speed is collected four times a day, the Tj value is specified as follows: for sandy land (desert), when the vegetation coverage is greater than 80%, Tj is 0; when the vegetation coverage is ≤80%, only the critical erosion wind speed (U25) corresponding to the vegetation coverage in Table 4 is exceeded. jt ) is accumulated over time for each level of wind speed. When the land use is grassland and shrub land, and the vegetation coverage is greater than 70%, Tj is taken as 0; when the vegetation coverage is ≤70%, only the critical erosion wind speed (U jt ) are accumulated over time for each level of wind speed. Table 4 is the accumulation of critical erosion wind speed Ujt under different vegetation coverage of sandy land (desert) and grass (shrub) land, as shown in Table 4: Table 4 Critical erosion wind speed U for sandy land (desert) and grass (shrub) land at different vegetation coverage jt

[0052] The calculation process of wind factor T is as follows: 1. Call the pre-processed time-sharing wind speed data in the local database, according to The wind speed is divided into 35 levels and the corresponding cumulative time in each half month is counted according to the level, and the cumulative time of each level in the half month is saved as a CSV file. 2. The half-month wind speed cumulative time uses the ordinary Kriging spatial interpolation method to generate 35 wind speed half-month cumulative time t jmik Grid data. The coordinate system and resolution are set according to user input. 3. Call the pre-processed vegetation coverage data and land use data in the local database, use the land use grid as a mask to extract the half-month vegetation coverage mean matrix of grass and shrub land, calculate the mean vegetation coverage of each block, and obtain the critical erosion wind speed U of each block according to Table 4. jt , and then Ujt corrects the cumulative time t of each block jmik , perform the same calculation for cultivated land and sandy land, and save the revised half-month cumulative time grids of wind speed at all levels. 5. Accumulate the half-month cumulative time grids of wind speed at all levels according to the wind speed level to obtain the annual cumulative time grid T of each wind speed level. j , save the annual cumulative time grid T of each level of wind speed j (The horizontal and vertical coordinates are the geographic spatial position, unit: ° or m, each pixel value is the annual cumulative time grid of wind speed T j , unit: min).

[0053] ⑨Surface soil moisture factor W: The calculation formula of surface soil moisture factor W is as follows: ; Where W represents the surface soil moisture factor, which is a dimensionless parameter. Call the pre-processed soil moisture content image in the local database and follow Perform matrix calculation to obtain the topsoil moisture factor grid and save the W factor grid (horizontal and vertical coordinates are geographic space positions, units: ° or m, each pixel value is the topsoil moisture factor, dimensionless).

[0054] ⑩Surface roughness factor Z0: Call the pre-processed land use data in the local database, and assign the surface roughness (Z0) according to the type of cultivated land in the cultivated land category according to Tables 5 to 7. Among them, Table 5 is the surface roughness (Z0) table of plowed and harrowed cultivated land, Table 6 is the surface roughness (Z0) table of plowed but not harrowed cultivated land, and Table 7 is the surface roughness (Z0) table of stubble cultivated land. Convert the assigned cultivated land vector into a raster and save it (the horizontal and vertical coordinates are the geographic space position, unit: ° or m, and the value of each pixel is the surface soil moisture factor, unit: cm).

[0055]

[0056] ⑪ Wind erosion calculation Q: Call the annual cumulative time grid Tj of each level of wind speed, the half-month vegetation coverage mean grid V, the surface soil humidity factor grid W, and the surface roughness factor Z0, according to the formula , and Calculate and obtain wind erosion grid data of cultivated land, grass and shrubs, sandy land, etc., and obtain wind erosion data of all land types by grid matrix superposition. Save the wind erosion modulus grid (horizontal and vertical coordinates are geographic spatial positions, units: ° or m, and each pixel value is the wind erosion modulus, units: t·hm -2 ·a -1 ).

[0057] Specifically, the method uses the CNN model to construct a soil erosion calculation model, and finally selects the soil erosion calculation model with the best performance through the evaluation of the trained soil erosion calculation model to calculate the soil erosion in the task area. Among them, the target data of the soil erosion calculation model can adopt the calculation results of the hydraulic erosion model and the wind erosion model in the "Technical Guidelines for Dynamic Monitoring of Soil and Water Loss in 2023" compiled by the Soil and Water Conservation Monitoring Center of the Ministry of Water Resources (that is, the output object of the soil erosion calculation model (training target or target data) is consistent with the output object of the hydraulic erosion model and the wind erosion model (that is, the calculation results of the hydraulic erosion model and the wind erosion model)). The input data (feature data) of the soil erosion calculation model uses the annual average rainfall raster, the annual cumulative time raster of various levels of wind speed, the vegetation coverage raster, DEM, soil moisture content raster, surface runoff raster, soil texture raster, soil thickness raster, soil permeability raster, soil distribution raster and / or land use status raster. The process of constructing a soil erosion calculation model using the CNN model includes: ① Data preparation: call the local database of hydraulic erosion grid, wind erosion grid, annual average rainfall grid, annual cumulative time grid of wind speed at all levels, vegetation coverage grid, DEM after filling, soil moisture content grid, surface runoff grid, soil density grid, soil texture grid, soil thickness grid, soil porosity grid, and land use status grid. The specific process is as follows: 1. According to the coordinate system, resolution and range of the target data hydraulic erosion grid and wind erosion grid, all input data are projected, resampled and cropped. 2. The input data is normalized and all input data are superimposed into a multi-channel raster file. 3. Use the ImageDataGenerator in the Keras library for data enhancement, that is, generate new data samples through rotation, translation, scaling and flipping. 4. Group the input data with the target data, including training set (80%) and test set (20%).

[0058] ②Model training and evaluation: Construct a convolutional neural network (CNN) model, such as Figure 4 As shown, Figure 4 It is a structural schematic diagram of a CNN model provided by an embodiment of the present invention. The CNN model structure includes an input layer (Input Layer), multiple convolution layers (Conv Layer), a pooling layer (Pooling Layer) and a random dropout layer (Dropout Layer). Finally, the multi-dimensional data is flattened into one dimension through a flattening layer (Flatten Layer), and the data is processed through a fully connected layer (Dense Layer), and the output layer finally generates a prediction result. The convolution layer uses a 3x3 convolution kernel, the activation function is Sigmoid, and regularization is introduced to prevent overfitting. The pooling layer uses a 2x2 maximum pooling layer to reduce the data dimension, and the Dropout layer randomly discards a certain proportion of neurons during the training process to further prevent overfitting. The CNN model optimizer uses RMSprop to dynamically adjust the learning rate of each parameter, and sets the mean square error (MSE) as the loss function. Before training, first use the Keras Tuner library to define Bayesian optimization, set the hyperparameter search space, and then pass the training and validation data sets through the tuner.search method to obtain the best hyperparameter combination, and finally rebuild and train the model based on these hyperparameters. In addition, in order to control the training process and prevent overfitting, early stopping and learning rate scheduler (ReduceLROnPlateau) are used. After training, the performance of the model is measured by evaluating the mean square error (MSE) and accuracy on the training set and test set. In addition, the mean absolute error (MAE) on the test set is calculated to evaluate the prediction accuracy of the model in practical applications. The model training process is as follows Figure 5 As shown, Figure 5 The present invention provides a training flow chart of a soil erosion calculation model according to an embodiment of the present invention.

[0059] ③Soil erosion calculation model deployment: After the soil erosion calculation model is trained, save it for subsequent loading and use. Save format: HDF5 file (.h5) contains complete model information, or TensorFlow Serving format. 1. Save the model as an HDF5 file, including saving the model structure, weights, and training configuration; save the model in TensorFlow SavedModel format for use by the service system TensorFlow Serving. When used in a production environment, you need to restore the saved model, load the model to read the saved file, and restore the model's structure and weights. 1. Load the model from the HDF5 file, read the saved .h5 file, and restore the model structure and weights; load the model from the TensorFlow SavedModel format, read the saved SavedModel directory, and restore the model structure and weights.

[0060] Local deployment uses Flask (Flask is a lightweight Python web framework that can be used to create REST APIs, making the model available through HTTP requests.) 1. Create a Flask application and load the saved model. 2. Create a prediction endpoint that receives input data and returns prediction results. 3. Start the Flask application and listen to the specified port.

[0061] Cloud deployment uses the TensorFlow Serving service system and makes prediction requests through REST API or gRPC. 1. Upload the model to the cloud service platform. 2. Configure and deploy the model service, set endpoints and resources. 3. Use the API provided by the platform to make prediction requests.

[0062] ④ Example of soil erosion calculation model: Taking the construction of a hydraulic erosion calculation model (i.e. soil erosion calculation model) in a certain area in Northwest China as an example, the soil erosion calculation model is trained. The specific process is as follows: 1. Use the mission area range vector to obtain the original data from the local database. The original data includes: hydraulic erosion raster, annual average rainfall raster, vegetation coverage raster, DEM after filling, surface runoff raster, soil density raster, soil texture raster, soil thickness raster, soil porosity raster, and land use status raster. The target data is set to a coordinate system of CGCS2000 and a resolution of 10m. First, convert the original data projection to CGCS2000, and then resample it to 10m. Finally, cut the original data according to the mission area and align the raster matrix, such as Figure 6 As shown, Figure 6 It is a schematic diagram of data acquisition and preprocessing provided by one embodiment of the present invention.

[0063] 2. First, normalize the 12 rasters, and then overlay the data into a multi-channel raster file.

[0064] 3. Use the ImageDataGenerator in the Keras library to perform a series of random transformations on the input data to generate new data samples. These transformations include but are not limited to: Rotation: Rotate the image at a random angle. Translation (width_shift, height_shift): Shift the image horizontally or vertically. Zoom: Randomly scale the image. Flip (horizontal_flip): Randomly flip the image horizontally. This improves the generalization ability of the model and prevents overfitting of the model, such as Figure 7 As shown, Figure 7 It is a schematic diagram of data enhancement provided by an embodiment of the present invention.

[0065] 4. Group the input data and target data into training set (80%) and test set (20%).

[0066] 5. Build a CNN model. Define Hyperparameter Search Space, and use Keras Tuner to set the search range of hyperparameters: the number of filters of the convolution layer (16 to 512), the size of the convolution kernel (1×1 to 8×8), the activation function (ReLU, Sigmoid, Tanh, etc.), the maximum pooling window (2×2 to 4×4), the Dropout random loss rate (0.1 to 0.9), the regularization parameter (0.0001 to 0.1), etc. Use the tuner.search method to search for hyperparameters and finally obtain the best hyperparameter combination: the number of filters of the convolution layer is 32, 64, 128, 256, the convolution kernel size is 3×3, the activation function is Sigmoid, the maximum pooling window is 2×2, the Dropout random loss rate is 0.6, 0.6, 0.6, 0.6, 0.5, and the regularization parameter is 0.01.

[0067] 6. Load the best hyperparameter combination into the model and start training. The EarlyStopping training round number patience is preset to 10. The ReduceLROnPlateau learning rate scheduler factor is preset to 0.1, and the training round number patience is preset to 5. The upper limit of the model training epochs is 100000, and the gradient update sample size batch_size is 32. During the convergence process, the accuracy rate gradually increases and the loss rate gradually decreases. Early Stopping is triggered after the number of training times reaches 936. Output the final accuracy curve and loss value curve, such as Figure 8and Fig. 9 As shown, Figure 8 is an accuracy curve diagram provided by an embodiment of the present invention, Fig. 9 It is a loss value curve diagram provided by an embodiment of the present invention.

[0068] 7. After training is complete, call the model.save() method to save the model as a .h5 file. Create Flask to load and use the saved model for prediction. Load the original data of a mission area in the northwest in 2023 into the Flask application for calculation, and compare the calculation results of the original model with those of the CNN model. Fig.10 As shown, Fig.10 is a schematic diagram of hydraulic erosion calculated by a CNN model provided by one embodiment of the present invention. In addition, the residual values ​​of the two data are compared, such as Fig.11 As shown, Fig.11 It is a residual graph of model calculation results provided by an embodiment of the present invention.

[0069] S120: inputting the preprocessed environmental data of the target area into the trained soil erosion prediction model to predict the soil erosion grade map of the target area within a preset future time length, and obtaining soil erosion prediction data of the target area; wherein the preset future time length represents a future time of a preset time length with the current moment as the time starting point; In some implementations of the present embodiment, the training rules of the above-mentioned soil erosion prediction model include: obtaining multiple second preset object data and soil erosion data corresponding to each second preset object data in a preset database, and performing projection transformation, resampling, cropping, normalization and data enhancement processing on each second preset object data in turn; wherein the time length between the collection time of each second preset object data is a preset time length; based on each second preset object data after data enhancement processing and the soil erosion data corresponding to each second preset object data, a training sample set is generated; based on the training sample set, the soil erosion prediction model is trained to obtain a trained soil erosion prediction model; wherein the soil erosion prediction model is constructed based on the spatiotemporal convolutional neural network algorithm.

[0070] Specifically, this method uses a spatiotemporal convolutional network (STGCN) to build a soil erosion prediction model, trains the STGCN model (i.e., soil erosion prediction model) through certain time series data, evaluates the training effect based on the average loss and mean square error, and selects the model with the best performance as the soil erosion prediction model. Among them, the input data and target data use the calculation results of the hydraulic erosion model and wind erosion model in the "Technical Guidelines for Dynamic Monitoring of Soil and Water Loss in 2023" compiled by the Soil and Water Conservation Monitoring Center of the Ministry of Water Resources. The process of building a prediction model using a spatiotemporal convolutional network (STGCN) is as follows: ① Data preparation: call the water erosion grid and wind erosion grid in the local database. 1. Perform projection conversion, resampling, and cropping on all input data according to the coordinate system, resolution, and range input by the user. 2. Normalize the input data. 3. Enhance the data by random noise, flipping, and rotation. 4. Group the input data and the target data into a training set (80%) and a test set (20%).

[0071] ② Training and evaluation of soil erosion prediction model: Construct STGCN model, such as Fig.12 As shown, Fig.12 It is a schematic diagram of the structure of an STGCN model provided by an embodiment of the present invention. The STGCN model includes an input layer (Input), a graph convolutional network (GCN), and a temporal convolutional network (TCN). Finally, the data is processed by a fully connected layer (FC Layer) to generate a prediction result. The input layer sets the input data shape to [batch_size, num_time_steps, num_channels, height, width], where batch_size is the batch size, num_time_steps is the time step, num_channels is the number of channels, height is the image height, and width is the image width. The graph convolutional network is divided into multiple layers of graph convolution (GCNConv), and each layer uses batch normalization (Batch Normalization) and random dropout (Dropout) to ensure stability. The temporal convolutional network (TCN) is a one-dimensional convolution (Conv1d), batch normalization (Batch Normalization) and random dropout (Dropout). The fully connected layer outputs the final result of the shape of [batch_size, height, width, out_channels] through Linear. The model optimizer uses ReduceLROnPlateau to dynamically adjust the learning rate of each parameter. First, create the edge index matrix (EdgeConnections) and the weight matrix (Weight Martix), and then start training the model. In addition, the residual connection (Residual Connection) is used to retain feature information and improve performance. After the training is completed, the performance of the model is measured by evaluating the mean square error (MSE) and loss value (Loss) on the training set and the test set. In addition, the mean absolute error (MAE) on the test set is calculated to evaluate the prediction accuracy of the soil erosion prediction model in practical applications.

[0072] ③Soil erosion prediction model deployment: After the soil erosion prediction model is trained, save it for subsequent loading and use. Save format: HDF5 file (.h5) contains complete model information, or TensorFlowServing format. 1. Save the model as an HDF5 file, including the model structure, weights, and training configuration; save the model in TensorFlow SavedModel format for use by the service system TensorFlow Serving.

[0073] When used in a production environment, you need to restore the saved model, load the model, read the saved file, and restore the model's structure and weights. 1. Load the model from the HDF5 file, read the saved .h5 file, and restore the model structure and weights; load the model from the TensorFlow SavedModel format, read the saved SavedModel directory, and restore the model structure and weights.

[0074] Local deployment uses Flask (Flask is a lightweight Python web framework that can be used to create REST APIs, making the model available through HTTP requests.) 1. Create a Flask application and load the saved model. 2. Create a prediction endpoint that receives input data and returns prediction results. 4. Start the Flask application and listen to the specified port.

[0075] Cloud deployment uses the TensorFlow Serving service system and makes prediction requests through REST API or gRPC. 1. Upload the model to the cloud service platform. 2. Configure and deploy the model service, set endpoints and resources. 3. Use the API provided by the platform to make prediction requests.

[0076] ④ Example of soil erosion prediction model: Taking the construction of a hydraulic erosion prediction model (i.e. soil erosion prediction model) in a certain area in Northwest China as an example, the soil erosion prediction model is trained. The specific process is as follows: 1. Use the mission area range vector to obtain the original data from the local database. The original data includes: 2019-2023 hydraulic erosion grid. The target data is set to CGCS2000 coordinate system and 10m resolution. First, convert the original data projection to CGCS2000, and then resample to 10m. Finally, cut the original data according to the mission area and align the grid matrix.

[0077] 2. Normalize the raster.

[0078] 3. Use torchvision.transforms in the PyTorch library to perform a series of random transformations on the input data to generate new data samples. These transformations include but are not limited to: Noise: Add noise to the image. Rotation: Rotate the image at a random angle. Translation (width_shift, height_shift): Shift the image horizontally or vertically. Zoom: Randomly scale the image. Flip (horizontal_flip): Randomly flip the image horizontally. This improves the generalization ability of the model and prevents overfitting of the model.

[0079] 4. The 2019-2022 raster is used as input data, and the 2023 data is used as target data. The data are grouped into training set (80%) and test set (20%).

[0080] 5. Build the STGCN model. The GCN network consists of three layers of graph convolution. Each layer includes GCN convolution (GraphConvolution), batch normalization (Batch Normalization), and random dropout (Dropout). The input is 32, 64, and 128 dimensions, respectively, and the random dropout rate is 0.6. The TCN network includes one-dimensional convolution (Conv1d), batch normalization (BatchNormalization), random dropout (Dropout), the convolution kernel size kernel_size is 3, and the padding is 1. The residual connection (Residual Connection) performs a two-dimensional convolution (Conv2d) on the input data and adds it to the TCN network output. The fully connected layer (FC) uses a linear method.

[0081] 6. First create the graph structure and weight matrix, and then start training. During the training process, the preset ReduceLROnPlateau learning rate scheduler factor is 0.5, and the number of training rounds patience is 3. The upper limit of the model training epochs is 1000, and the gradient update sample size batch_size is 32. During the convergence process, the accuracy rate gradually increases and the loss rate gradually decreases. Output the loss value curve, such as Fig.13 As shown, Fig.13 Another loss value curve diagram provided by an embodiment of the present invention.

[0082] 7. After training is complete, call the model.save() method to save the model as a .h5 file. Create a Flask application to load and use the saved model for prediction. Load the original data of a mission area in the northwest from 2019 to 2023 into the Flask application for calculation, and compare the 2023 hydraulic erosion data with the calculation results of the prediction model. Fig.14 As shown, Fig.14 is a schematic diagram of hydraulic erosion prediction obtained by a STGCN model provided by an embodiment of the present invention. In addition, the residual values ​​of the two data are compared, such as Fig.15 As shown, Fig.15 is an actual value and prediction residual diagram provided by an embodiment of the present invention, predicting the change of the area of ​​mild or above hydraulic erosion from 2024 to 2028, such as Fig.16 As shown, Fig.16 It is a schematic diagram of the change of hydraulic erosion area provided by an embodiment of the present invention.

[0083] S130: Perform soil erosion analysis based on the soil erosion calculation results, soil erosion prediction data and a preset database of the target area to obtain soil erosion analysis results.

[0084] Specifically, the method inputs the preprocessed environmental data of the target area into the trained soil erosion calculation model to generate the current soil erosion grade map of the target area, thereby obtaining the soil erosion calculation result of the target area. The preprocessed environmental data of the target area is input into the trained soil erosion prediction model to predict and generate the soil erosion grade map of the target area within a preset future time length, thereby obtaining the soil erosion prediction data of the target area. Soil erosion analysis is performed in combination with the soil erosion calculation results, soil erosion prediction data and preset database of the target area to obtain the soil erosion analysis results. Thereby, the purpose of dynamic monitoring of soil erosion in multiple cross-regional small watersheds is achieved, and soil erosion in different regions is automatically calculated and soil erosion changes are predicted, thereby achieving the purpose of analyzing the future change trend of soil erosion in different regions, so as to provide support for the generation of soil and water conservation measures.

[0085] In some implementations of the present embodiment, the soil erosion analysis includes spatial distribution analysis, change trend analysis, overlay analysis, erosion risk analysis and factor association analysis; wherein the rules of spatial distribution analysis include: based on the soil erosion calculation results and soil erosion prediction data of the target area, drawing a hydraulic erosion spatial distribution map and a wind erosion spatial distribution map; comparing the hydraulic erosion spatial distribution map and the wind erosion spatial distribution map to determine the average erosion value, maximum erosion value, minimum erosion value, the location of the area reaching severe erosion conditions and the spatial distribution characteristics of the location of the area; the rules of change trend analysis include: based on the soil erosion calculation results, soil erosion prediction data and preset database of the target area, reading the soil erosion data of the target area in different years, and based on the soil erosion data of the target area in each year, counting the corresponding soil erosion values ​​greater than the predicted soil erosion values ​​in each year. Set the erosion area of ​​the mild erosion intensity threshold; draw the erosion area change trend map of the target area based on the erosion area corresponding to each year that is greater than the preset mild erosion intensity threshold; the rules of superposition analysis include: based on the spatial distribution map of hydraulic erosion and the spatial distribution map of wind erosion, mark the areas only affected by hydraulic erosion, the areas only affected by wind erosion, and the areas affected by both hydraulic erosion and wind erosion respectively; the rules of erosion risk analysis include: based on the soil erosion calculation results of the target area and the preset erosion risk level threshold, evaluate the erosion risk level of each area in the target area; the rules of factor association analysis include: based on the soil erosion calculation results of the target area, soil erosion prediction data, preset database and pre-determined influencing factors, use correlation analysis and / or regression analysis to determine the influencing factors that have the greatest impact on soil erosion.

[0086] Specifically, after completing the calculation and prediction of hydraulic and wind erosion, this method can use the results (i.e. soil erosion calculation results, soil erosion prediction data) and the original data of the database (i.e. soil erosion data in the preset database) for data analysis: ① Spatial distribution analysis: 1. Draw the spatial distribution maps of hydraulic erosion and wind erosion for intuitive comparison. Identify the areas with severe erosion and analyze their spatial distribution characteristics. 2. Count and compare the erosion intensity of different areas, and calculate the average erosion value, maximum erosion value, and minimum erosion value of each area. Through these statistics, evaluate the severity and distribution characteristics of hydraulic erosion and wind erosion.

[0087] ② Trend analysis: Read the erosion data of different years and count the erosion area that is greater than a certain threshold (such as mild erosion intensity) each year. Draw these data into a trend chart to analyze the changing pattern of erosion area between different years.

[0088] ③ Overlay analysis: Overlay the hydraulic erosion map and the wind erosion map to identify areas affected by both types of erosion. Different colors or symbols can be set to identify areas affected only by hydraulic erosion, areas affected only by wind erosion, and areas affected by both types of erosion.

[0089] ④ Erosion risk analysis: Evaluate the erosion risk level of each area based on the erosion value. Different risk level thresholds can be set to divide the area into low-risk, medium-risk and high-risk areas. Through risk assessment, identify high-risk areas that need to be focused on governance and protection.

[0090] ⑤ Factor correlation analysis: Combine other influencing factors (such as precipitation, wind speed, land use type, topography, etc.) to analyze the relationship between these factors and erosion. The key factors that have the greatest impact on erosion can be identified through correlation analysis, regression analysis, and other methods.

[0091] ⑥Data analysis example: Taking a mission area in the northwest as an example, the hydraulic and wind erosion data from 2019 to 2023 and the forecast data from 2024 to 2028 are selected for analysis. The specific process is as follows: 1. Spatial distribution analysis: Identify areas with severe erosion through the spatial distribution map of hydraulic and wind erosion from 2019 to 2023. Figures 17 to 21 As shown, Fig.17 is a spatial distribution map of hydraulic and wind erosion provided by an embodiment of the present invention, Fig.18 is another spatial distribution map of hydraulic and wind erosion provided by one embodiment of the present invention, Fig.19 is another spatial distribution map of hydraulic and wind erosion provided by an embodiment of the present invention. Fig. 20 is another spatial distribution map of hydraulic and wind erosion provided by an embodiment of the present invention, Fig.21 Another spatial distribution map of water and wind erosion provided by an embodiment of the present invention. When the erosion value is greater than 1000t / km 2 When the erosion reaches a level above mild erosion, corresponding engineering measures and plant measures must be taken.

[0092] 2. Trend analysis: Analyze the changing pattern of erosion area between different years, and count the annual erosion area greater than 1000t / km 2 ·a erosion area, generate erosion area change curve, such as Fig. 22 As shown, Fig. 22 It is another schematic diagram of the erosion area change curve provided by an embodiment of the present invention, wherein 2019-2023 is the actual erosion area change curve, and 2024-2028 is the predicted erosion area change curve.

[0093] 3. Overlay analysis: Identify the areas affected by both hydraulic erosion and wind erosion. Overlay the hydraulic erosion map and the wind erosion map, marking the area mainly affected by wind erosion as 1, the area mainly affected by hydraulic erosion as 2, and the area affected by both erosion as 3. Fig.23 As shown, Fig.23 It is a schematic diagram of erosion superposition provided by one embodiment of the present invention. Comprehensive management: In areas affected by both hydraulic and wind erosion, comprehensive management measures are adopted, combining vegetation restoration, engineering measures and farmland improvement methods. Key protection: In key protection areas, multi-level and diversified protection measures are adopted to improve the soil and water conservation effect of the region.

[0094] 4. Erosion risk analysis: Based on the erosion value, assess the erosion risk level of each area and identify high-risk areas, such as Fig.24 As shown, Fig.24 It is a schematic diagram of erosion risk levels provided by an embodiment of the present invention, wherein 2 or above is a high-risk area, less than 2 and greater than 1 is a medium-risk area, and less than 1 is a low-risk area. Management of high-risk areas: Priority is given to implementing soil and water conservation projects in high-risk areas, such as building retaining walls, drainage ditches and sedimentation ponds. Prevention in medium-risk areas: Take preventive measures in medium-risk areas, such as vegetation restoration and surface cover, to reduce erosion risks. Maintenance of low-risk areas: Maintain existing soil and water conservation measures in low-risk areas to prevent increased erosion risks.

[0095] 5. Factor correlation analysis: The correlation between hydraulic erosion and related factors from 2019 to 2023 is calculated, as shown in Table 8. Table 8 is a correlation table of each factor: Table 8 Correlation table of each factor

[0096] Among them, the correlation results show that the correlation coefficient between rainfall and hydraulic erosion is between 0.2516 and 0.3206, which is positively correlated. This indicates that when rainfall increases, the intensity of hydraulic erosion also increases. This is in line with general expectations, because rainfall directly affects surface runoff, which leads to increased soil erosion. The correlation coefficient between topography and hydraulic erosion is between 0.3134 and 0.3206, which is positively correlated. This shows that topography (mainly slope and elevation difference) has a significant effect on hydraulic erosion. The steeper the terrain, the faster the water flow, and the greater the possibility of soil erosion. The correlation coefficient between land use and hydraulic erosion is between 0.4248 and 0.4299, showing a strong positive correlation. Different types of land use have different effects on soil erosion. For example, cultivated land and bare land are more prone to soil erosion than forests and grasslands. Changes in land use types have a significant impact on hydraulic erosion, and reasonable land use can effectively slow down soil erosion. The correlation coefficient between vegetation coverage and hydraulic erosion ranges from -0.0530 to -0.0422, showing a weak negative correlation. This indicates that when the vegetation coverage is higher, the intensity of hydraulic erosion is lower. Vegetation reduces soil erosion by slowing down runoff velocity and increasing soil stability, so increasing vegetation coverage is an effective soil and water conservation measure.

[0097] In some implementations of this embodiment, after obtaining the soil erosion analysis results, the method further includes: determining a soil and water conservation plan based on the soil erosion analysis results.

[0098] Exemplary soil and water conservation programs may include: 1. Rationally plan land use to avoid farming and construction activities in areas with large slopes. These areas can be planned as forests or grasslands to reduce soil erosion. 2. Implement agricultural measures such as crop rotation and intercropping to improve the sustainability of land use and reduce the exposure time of bare land. 3. Increase vegetation coverage and promote measures such as afforestation and grass soil consolidation, especially in areas with severe soil and water loss, to slow down soil erosion by increasing vegetation coverage. 4. Adopt biological measures such as slope protection vegetation and bank protection plants to strengthen the protection of slopes and river banks. 5. Strengthen the construction of soil and water conservation projects to build terraces on slopes and mountains to reduce the slope and reduce the water flow rate, thereby reducing soil erosion. Construct soil and water conservation facilities such as interception ditches and drainage ditches to reasonably guide and discharge runoff and prevent water from eroding the soil in a concentrated manner. 6. Establish a long-term monitoring system to regularly monitor factors such as rainfall, terrain changes, land use and vegetation coverage to promptly discover and solve soil and water loss problems.

[0099] Fig.25 FIG. 1 is a block diagram of a soil erosion analysis system provided by an embodiment of the present invention. Fig.25 As shown, an embodiment of the present invention provides a soil erosion analysis system, comprising: A soil erosion calculation module is used to input the pre-processed environmental data of the target area into the trained soil erosion calculation model to generate a current soil erosion grade map of the target area and obtain the soil erosion calculation result of the target area; A soil erosion prediction module is used to input the pre-processed environmental data of the target area into the trained soil erosion prediction model to predict the soil erosion grade map of the target area within a preset future time length, and obtain the soil erosion prediction data of the target area; wherein the preset future time length represents the future time of the preset time length with the current moment as the time starting point; The soil erosion analysis module is used to perform soil erosion analysis based on the soil erosion calculation results, soil erosion prediction data and preset database of the target area to obtain soil erosion analysis results.

[0100] Specifically, the system inputs the preprocessed environmental data of the target area into the trained soil erosion calculation model to generate the current soil erosion grade map of the target area, thereby obtaining the soil erosion calculation results of the target area. The preprocessed environmental data of the target area is input into the trained soil erosion prediction model to predict and generate the soil erosion grade map of the target area within a preset future time length, thereby obtaining the soil erosion prediction data of the target area. Soil erosion analysis is performed in combination with the soil erosion calculation results, soil erosion prediction data and preset database of the target area to obtain the soil erosion analysis results. Thereby, the purpose of dynamic monitoring of soil erosion in multiple cross-regional small watersheds is achieved, and soil erosion in different regions is automatically calculated and soil erosion changes are predicted, thereby achieving the purpose of analyzing the future change trend of soil erosion in different regions, so as to provide support for the generation of soil and water conservation measures.

[0101] The optional implementation modes of the embodiments of the present invention are described in detail above in conjunction with the accompanying drawings. However, the embodiments of the present invention are not limited to the specific details in the above implementation modes. Within the technical concept of the embodiments of the present invention, various simple modifications can be made to the technical scheme of the embodiments of the present invention, and these simple modifications all belong to the protection scope of the embodiments of the present invention.

[0102] It should also be noted that the various specific technical features described in the above specific embodiments can be combined in any suitable manner without contradiction. To avoid unnecessary repetition, the embodiments of the present invention will not further describe various possible combinations.

[0103] Those skilled in the art can understand that all or part of the steps in the above-mentioned embodiment method can be completed by instructing the relevant hardware through a program, and the program is stored in a storage medium, including a number of instructions to enable a single-chip microcomputer, a chip or a processor to execute all or part of the steps of the method described in each embodiment of the present application. The aforementioned storage medium includes: U disk, mobile hard disk, read-only memory (ROM, Read-Only Memory), random access memory (RAM, Random Access Memory), disk or optical disk and other media that can store program codes.

[0104] In addition, various implementations of the embodiments of the present invention may be arbitrarily combined, and as long as they do not violate the concept of the embodiments of the present invention, they should also be regarded as the contents disclosed in the embodiments of the present invention.

Claims

1. A soil erosion analysis method, characterized in that: include: The pre-processed environmental data of the target area is input into the trained soil erosion calculation model to generate a current soil erosion grade map of the target area and obtain the soil erosion calculation result of the target area; The pre-processed environmental data of the target area is input into the trained soil erosion prediction model to predict the soil erosion grade map of the target area within a preset future time length, so as to obtain the soil erosion prediction data of the target area; wherein the preset future time length represents the future time of the preset time length with the current moment as the time starting point; Soil erosion analysis is performed based on the soil erosion calculation results, soil erosion prediction data and preset database of the target area to obtain soil erosion analysis results.

2. The soil erosion analysis method according to claim 1, characterized in that: The acquisition rules of the environmental data of the target area include: Input the search range data of the target area into the preset database for matching to obtain the environmental data of the target area; wherein, The search range data of the target area includes geometric boundary parameters of the target area, data collection start and end date parameters, data resolution parameters, cloud threshold parameters and file format parameters; The preset database includes multiple data sources; The environmental data of the target area includes meteorological data, image data, hydrological data, soil data and land management data of the target area.

3. The soil erosion analysis method according to claim 2, characterized in that: The preprocessing rules of the environmental data of the target area include: The meteorological data, image data, hydrological data, soil data and land management data of the target area are preprocessed respectively; The preprocessing rules for the meteorological data of the target area include: Batch read the folders corresponding to the meteorological data of the target area to build a meteorological data list, and clean the meteorological data in the meteorological data list; The preprocessing rules of the image data of the target area include: When it is determined that the image data of the target area is complete, image processing is performed on the image data of the target area to obtain image processing data of the target area; wherein the image processing includes one or more of band stacking, geometric correction, atmospheric correction, image stitching and image cropping; Performing projection transformation and resampling on the image processing data of the target area in sequence to obtain resampled image data; Based on the data type of the image data of the target area, a corresponding image processing scheme is matched, and based on the resampled image data, the corresponding image processing scheme is executed to obtain the pre-processed image data of the target area; The preprocessing rules for soil data in the target area include: When it is determined that the soil data of the target area is complete, projection transformation and resampling are performed on the soil data of the target area in sequence to obtain resampled soil data; Generate and save soil data maps based on the resampled soil data; The preprocessing rules for land management data in the target area include: When it is determined that the land management data of the target area is complete, projection transformation and resampling are performed on the land management data of the target area in sequence to obtain resampled land management data; Generate and save land management data maps based on the resampled land management data.

4. The soil erosion analysis method according to claim 1, characterized in that: The training rules of the soil erosion calculation model include: Acquire a plurality of first preset object data and soil erosion data corresponding to each first preset object data in a preset database, and sequentially perform projection conversion, resampling, cropping and normalization processing on each first preset object data, so as to superimpose each first preset object data into a multi-channel raster file; Performing data enhancement processing on the multi-channel raster file, generating a data sample set based on the soil erosion data corresponding to each first preset object data and the corresponding multi-channel raster file after the data enhancement processing, and dividing the data sample set into a first training set and a first test set according to a preset ratio; Based on the Bayesian optimization algorithm and the pre-set hyperparameter search space, the hyperparameters of the first training set and the first test set are optimized to obtain the best hyperparameter combination; Based on the optimal hyperparameter combination, the soil erosion calculation model is trained to obtain a trained soil erosion calculation model; wherein the soil erosion calculation model is constructed based on a convolutional neural network algorithm.

5. The soil erosion analysis method according to claim 1, characterized in that: The training rules of the soil erosion prediction model include: Acquire multiple second preset object data and soil erosion data corresponding to each second preset object data in a preset database, and sequentially perform projection conversion, resampling, cropping, normalization and data enhancement processing on each second preset object data; wherein the time length between the collection time of each second preset object data is a preset time length; Generate a training sample set based on each second preset object data after data enhancement processing and the soil erosion data corresponding to each second preset object data; The soil erosion prediction model is trained based on the training sample set to obtain a trained soil erosion prediction model; wherein the soil erosion prediction model is constructed based on a spatiotemporal convolutional neural network algorithm.

6. The soil erosion analysis method according to claim 1, characterized in that: The soil erosion analysis includes spatial distribution analysis, change trend analysis, overlay analysis, erosion risk analysis and factor association analysis; The rules of the spatial distribution analysis include: Based on the soil erosion calculation results and soil erosion prediction data of the target area, draw the spatial distribution map of hydraulic erosion and the spatial distribution map of wind erosion; Compare the spatial distribution map of hydraulic erosion with the spatial distribution map of wind erosion to determine the average erosion value, maximum erosion value, minimum erosion value, location of the area reaching severe erosion conditions, and the spatial distribution characteristics of the location of the area; The rules for the change trend analysis include: Based on the soil erosion calculation results, soil erosion prediction data and preset database of the target area, the soil erosion data of the target area in different years are read, and based on the soil erosion data of the target area in each year, the erosion area corresponding to each year that is greater than the preset light erosion intensity threshold is counted; Based on the erosion area greater than the preset mild erosion intensity threshold corresponding to each year, draw a trend chart of the erosion area change in the target area; The rules of the overlay analysis include: Based on the spatial distribution map of hydraulic erosion and the spatial distribution map of wind erosion, the areas affected only by hydraulic erosion, the areas affected only by wind erosion, and the areas affected by both hydraulic erosion and wind erosion are marked respectively; The rules for erosion risk analysis include: Based on the soil erosion calculation results of the target area and the preset erosion risk level threshold, the erosion risk level of each location in the target area is evaluated; The rules of the factor association analysis include: Based on the soil erosion calculation results, soil erosion prediction data, preset database and pre-determined influencing factors in the target area, the influencing factors with the greatest impact on soil erosion are determined using correlation analysis and / or regression analysis.

7. The soil erosion analysis method according to claim 1, characterized in that: After obtaining the soil erosion analysis results, the method further includes: Determine soil and water conservation plans based on soil erosion analysis results.

8. A soil erosion analysis system, characterized in that: include: A soil erosion calculation module is used to input the pre-processed environmental data of the target area into the trained soil erosion calculation model to generate a current soil erosion grade map of the target area and obtain the soil erosion calculation result of the target area; A soil erosion prediction module is used to input the pre-processed environmental data of the target area into the trained soil erosion prediction model to predict the soil erosion grade map of the target area within a preset future time length, and obtain soil erosion prediction data of the target area; wherein the preset future time length represents the future time of the preset time length with the current moment as the time starting point; The soil erosion analysis module is used to perform soil erosion analysis based on the soil erosion calculation results, soil erosion prediction data and preset database of the target area to obtain soil erosion analysis results.

9. A machine-readable storage medium having instructions stored thereon, characterized in that: When the instruction is executed by a processor, the processor is configured to execute the soil erosion analysis method according to any one of claims 1 to 7.

10. An electronic device, comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein: When the processor executes the computer program, the soil erosion analysis method according to any one of claims 1 to 7 is implemented.

Images