Method and system for dynamic monitoring of soil erosion based on remote sensing image

Abstract

The application provides a kind of water and soil loss dynamic monitoring method and system based on remote sensing image, first, the multi-temporal remote sensing image dataset of target area is acquired, it contains multiple time period remote sensing image subset and surface cover information, then the surface feature extraction is carried out to the multi-temporal remote sensing image dataset, after obtaining the characteristics of vegetation coverage, terrain slope, soil exposure, construct the water and soil loss prediction model based on space-time correlation, input the above characteristics for prediction, generate water and soil loss grade distribution map, according to this water and soil loss grade distribution map, identify target loss risk area, generate management priority ranking and vegetation recovery strategy, finally feedback to monitoring platform triggers regional management task allocation operation, so as to realize the dynamic monitoring of water and soil loss, effectively plan management task.

Description

Technical Field

[0001] This invention relates to the field of artificial intelligence technology, and more specifically, to a method and system for dynamic monitoring of soil erosion based on remote sensing images. Background Technology

[0002] Soil erosion, a global ecological and environmental problem, poses a serious threat to land productivity, water resource quality, and ecosystem balance. It not only leads to decreased soil fertility and reduced arable land but can also trigger a series of chain reactions such as siltation and flooding, severely impacting regional economic sustainability and the human living environment. Therefore, timely and accurate understanding of the dynamic changes in soil erosion is crucial for formulating scientifically sound soil and water conservation strategies and ecological restoration plans.

[0003] However, traditional methods for monitoring soil erosion mainly rely on ground surveys and field measurements, such as setting up runoff plots and calculating soil erosion. While these methods can provide relatively accurate local data, they have significant limitations. On the one hand, ground surveys require substantial manpower, resources, and time, and their monitoring scope is limited, making it difficult to achieve large-scale, high-frequency dynamic monitoring. On the other hand, field measurements are easily constrained by natural conditions such as topography and climate, making it difficult to conduct effective monitoring in some remote and dangerous areas. Summary of the Invention

[0004] In view of the aforementioned problems, and in conjunction with the first aspect of the present invention, embodiments of the present invention provide a method for dynamic monitoring of soil erosion based on remote sensing imagery, the method comprising:

[0005] Acquire a multi-temporal remote sensing image dataset of the target area, wherein the multi-temporal remote sensing image dataset includes remote sensing image subsets of multiple time periods, and each remote sensing image subset contains land cover information of at least one spectral band.

[0006] The surface feature extraction process is performed on the multi-temporal remote sensing image dataset to obtain the vegetation cover features, topographic slope features, and soil exposure features of each remote sensing image subset.

[0007] A spatiotemporal correlation-based soil erosion prediction model is constructed. The vegetation coverage characteristics, terrain slope characteristics and soil exposure characteristics are input into the soil erosion prediction model to predict soil erosion and generate a soil erosion level distribution map of the target area.

[0008] Based on the soil erosion level distribution map, target erosion risk areas are identified, and a priority ranking of treatment and vegetation restoration strategies are generated for the target erosion risk areas.

[0009] The governance priority ranking and vegetation restoration strategy are fed back to the monitoring platform to trigger the regional governance task allocation operation.

[0010] In another aspect, embodiments of the present invention also provide a dynamic monitoring system for soil erosion based on remote sensing imagery, including a processor and a machine-readable storage medium connected to the processor. The machine-readable storage medium is used to store programs, instructions, or code, and the processor is used to execute the programs, instructions, or code in the machine-readable storage medium to implement the above-described method.

[0011] Based on the above, this application embodiment achieves comprehensive, accurate, and dynamic monitoring and management guidance for soil erosion in the target area. First, a multi-temporal remote sensing image dataset of the target area is acquired, covering subsets of remote sensing images from multiple time periods. Each subset contains land cover information for at least one spectral band, capturing subtle changes in the land surface over different time periods, thus more accurately reflecting the dynamic process of soil erosion. Based on this, surface feature extraction processing is performed on the multi-temporal remote sensing image dataset to obtain vegetation cover features, topographic slope features, and soil exposure features. This characterizes the potential correlation between land surface conditions and soil erosion from different perspectives. Vegetation cover reflects the protective capacity of surface vegetation for the soil, topographic slope affects the erosive effect of water flow on the soil, and soil exposure directly reflects the degree to which the exposed soil is susceptible to erosion. Furthermore, a spatiotemporal correlation-based soil erosion prediction model is constructed. The extracted features are input into the model to predict soil erosion, generating a soil erosion level distribution map of the target area. This model fully considers the impact of temporal and spatial factors on soil erosion, accurately simulating the evolution trend of soil erosion at different temporal and spatial scales. The generated soil erosion level distribution map visually displays the severity and spatial distribution of soil erosion in the target area. Based on this map, target erosion risk areas can be accurately identified, and priority ranking and vegetation restoration strategies for these areas can be generated, thus avoiding blind governance and resource waste, and improving governance efficiency and effectiveness. By feeding back the priority ranking and vegetation restoration strategies to the monitoring platform to trigger regional governance task allocation, seamless integration of monitoring and governance is achieved. In other words, through the acquisition of multi-temporal remote sensing image data, accurate extraction of surface features, construction of a spatiotemporal correlation-based prediction model, and the formulation and feedback of targeted governance strategies, comprehensive and dynamic monitoring and efficient governance of soil erosion are achieved. Attached Figure Description

[0012] Figure 1 This is a schematic diagram of the execution flow of the dynamic monitoring method for soil erosion based on remote sensing images provided in an embodiment of the present invention.

[0013] Figure 2This is a schematic diagram of exemplary hardware and software components of the dynamic monitoring system for soil erosion based on remote sensing images provided in an embodiment of the present invention. Detailed Implementation

[0014] The present invention will now be described in detail with reference to the accompanying drawings. Figure 1 This is a flowchart illustrating a method for dynamic monitoring of soil erosion based on remote sensing images, provided in one embodiment of the present invention. The following is a detailed description of this method for dynamic monitoring of soil erosion based on remote sensing images.

[0015] Step S110: Obtain a multi-temporal remote sensing image dataset of the target area. The multi-temporal remote sensing image dataset includes remote sensing image subsets for multiple time periods, and each remote sensing image subset contains land cover information for at least one spectral band.

[0016] For example, the target area could be a large area of ​​mountains, forests, farmland, or other areas with different land cover types. Taking a mountainous area as an example, this area may experience varying degrees of soil erosion. In order to accurately monitor the dynamic changes in soil erosion, it is necessary to acquire multi-temporal remote sensing images of the area.

[0017] A multi-temporal remote sensing image dataset consists of subsets of remote sensing images from multiple time periods. If we choose to monitor soil erosion in a mountainous area for 10 years, with each year as a time period, then there will be 10 subsets of remote sensing images from different time periods. Each subset of remote sensing images from different time periods is a collection of remote sensing images acquired at a specific time point within that time period.

[0018] Each subset of remote sensing imagery contains land cover information for at least one spectral band. Common spectral bands include blue, green, red, near-infrared, and shortwave infrared bands. Different spectral bands exhibit different reflectance characteristics for different land cover types. For example, vegetation has high reflectance in the near-infrared band, while water bodies have unique reflectance characteristics in the blue band. Information from these spectral bands can help identify and analyze different land cover types, such as vegetation, water bodies, and soil.

[0019] There are various ways to acquire multi-temporal remote sensing image datasets. Currently, many sources provide abundant remote sensing image data resources. In this embodiment, suitable remote sensing images can be selected from these data sources based on the geographical location of the target area and the required time range.

[0020] Taking the Landsat 8 satellite as an example, it is equipped with Operational Land Imager (OLI) and Thermal Infrared Sensor (TIRS), capable of providing image data in multiple spectral bands, including blue, green, red, near-infrared, and short-wave infrared. This embodiment uses the EarthExplorer platform to input the latitude and longitude range of the target area and the desired time range to search and download Landsat 8 image data that meets the criteria. Assuming the latitude and longitude range of the target mountainous area is 110°-112°E and 30°-32°N, this embodiment can select the summer (June-August) of each year from 2010 to 2020 as the data acquisition time range, because vegetation growth is vigorous during this period, which can more accurately reflect the land cover situation. Through searching and filtering, this embodiment can obtain 10 subsets of remote sensing images for different time periods, each subset containing Landsat 8 image data for the summer of the corresponding year.

[0021] After acquiring remote sensing image data, preprocessing is required to ensure data quality and usability. Preprocessing steps include radiometric correction, geometric correction, and atmospheric correction. Radiometric correction aims to eliminate the influence of sensor errors and atmospheric scattering and absorption on image radiometric values, ensuring that the image's radiometric values ​​accurately reflect the true reflectivity of the Earth's surface. Geometric correction matches the image's geographic coordinates with the actual geographic coordinates, eliminating geometric distortions. Atmospheric correction further eliminates the influence of the atmosphere on the image, improving image quality and analytical accuracy.

[0022] Step S120: Perform surface feature extraction processing on the multi-temporal remote sensing image dataset to obtain the vegetation cover features, topographic slope features, and soil exposure features of each remote sensing image subset.

[0023] After obtaining the multi-temporal remote sensing image dataset of the target area, the next step is to perform surface feature extraction processing on these datasets to obtain the vegetation cover features, topographic slope features, and soil exposure features of each remote sensing image subset. These features are of great value for analyzing soil erosion.

[0024] Step S121: Perform radiometric correction processing on the subset of remote sensing images to obtain standardized remote sensing images.

[0025] Before extracting surface features, a subset of remote sensing images needs to be radiometrically corrected. The purpose of radiometric correction is to eliminate the influence of sensor errors, atmospheric scattering, and absorption on the image radiometric values, so that the image radiometric values ​​can accurately reflect the true reflectivity of the surface.

[0026] Radiometric correction can be divided into relative radiometric correction and absolute radiometric correction. Relative radiometric correction mainly eliminates radiometric differences within an image, making the radiometric values ​​of different areas within the same image comparable. Absolute radiometric correction, on the other hand, converts the radiometric values ​​of the image into the true reflectance of the Earth's surface.

[0027] Taking Landsat 8 imagery data as an example, the radiometric values ​​in Landsat 8 imagery data are usually expressed as digital quantization (DN) values. In order to perform absolute radiometric correction, the DN values ​​need to be converted into radiance values, and then the radiance values ​​need to be converted into surface reflectance.

[0028] First, the DN values ​​are converted to radiance values. The metadata file of Landsat 8 imagery data contains the radiometric calibration coefficients for each band, including gain and bias. The radiance value can be calculated as follows: Radiance value = Gain × DN value + Bias. For example, for the blue band of Landsat 8, assuming a gain of 0.01, a bias of -0.1, and a DN value of 100 for a certain pixel, the radiance value of that pixel is 0.01 × 100 + (-0.1) = 0.9.

[0029] Then, the radiance value is converted into surface reflectance. This requires considering factors such as the solar zenith angle and the Earth-Sun distance. The solar zenith angle can be obtained from the image's metadata file, and the Earth-Sun distance can be calculated based on the date. The surface reflectance can be calculated using the following formula: Surface reflectance = (π × radiance value × Earth-Sun distance squared) / (Solar spectral irradiance × cos(solar zenith angle)). Assuming the solar spectral irradiance is 1000, the solar zenith angle is 30°, the Earth-Sun distance is 1.01 astronomical units, and the radiance value of the above pixel is 0.9, then the surface reflectance of this pixel is (π × 0.9 × 1.01²) / (1000 × cos(30°)) ≈ 0.0033.

[0030] By performing the above calculations on all pixels in each band, the DN values ​​of the entire remote sensing image subset can be converted into surface reflectance, resulting in a standardized remote sensing image. The radiometric values ​​of the standardized remote sensing image can more accurately reflect the true reflectance characteristics of the surface.

[0031] Step S122: Call the pre-trained vegetation index calculation model, and calculate the vegetation coverage feature based on the reflectance difference between the near-infrared band and the red band in the standardized remote sensing image, wherein the vegetation coverage feature is positively correlated with the reflectance difference.

[0032] In this embodiment, after obtaining the standardized remote sensing image, a pre-trained vegetation index calculation model can be used to calculate vegetation cover characteristics. A vegetation index is a numerical value calculated by combining reflectance values ​​from different spectral bands; it reflects the growth status and coverage of vegetation.

[0033] A commonly used vegetation index is the Normalized Difference Vegetation Index (NDVI), which is calculated based on the difference in reflectance between the near-infrared band and the red band. The formula for calculating NDVI is: NDVI = (Near-infrared band reflectance - Red band reflectance) / (Near-infrared band reflectance + Red band reflectance).

[0034] Taking standardized remote sensing imagery as an example, this embodiment can extract reflectance data from the near-infrared and red bands of the imagery. Assuming that the reflectance of a certain pixel in the near-infrared band is 0.8 and the reflectance of the same pixel in the red band is 0.2, then the NDVI value of that pixel is (0.8-0.2) / (0.8+0.2)=0.6.

[0035] The pre-trained vegetation index calculation model can be a simple program that takes near-infrared and red band reflectance data from standardized remote sensing imagery as input and calculates the NDVI value for each pixel according to the formula described above. By calculating for all pixels in the entire image, an NDVI image can be obtained, where the value of each pixel represents the vegetation cover characteristics at that location.

[0036] NDVI values ​​typically range from -1 to 1. A value closer to 1 indicates higher vegetation cover, while a value closer to -1 indicates lower vegetation cover. For example, an NDVI value of 0.8 indicates vigorous vegetation growth and high coverage in the area; an NDVI value of 0.2 indicates low vegetation cover and potential soil exposure.

[0037] Step S123: Perform digital elevation model inversion processing on the standardized remote sensing image, and calculate the terrain slope feature based on the elevation change rate between adjacent pixels, wherein the elevation change rate threshold of the terrain slope feature is dynamically adjusted according to a preset terrain classification rule.

[0038] To obtain terrain slope features, it is necessary to perform digital elevation model (DEM) inversion processing on standardized remote sensing imagery. A digital elevation model is a model representing surface elevation information and can be obtained through various methods, such as radar interferometry and lidar measurement. In this embodiment, DEM inversion can be performed using stereo image pairs from remote sensing imagery.

[0039] A stereo image pair refers to two remote sensing images of the same area taken from different angles. By matching and calculating stereo image pairs, the elevation information of each pixel can be obtained, thereby generating a digital elevation model.

[0040] Once the digital elevation model is obtained, the terrain slope characteristics can be calculated based on the rate of elevation change between adjacent pixels. Terrain slope refers to the degree of inclination at a point on the Earth's surface, reflecting the undulation of the terrain.

[0041] A common method for calculating terrain slope is a difference algorithm based on a digital elevation model (DEM). In this embodiment, for each cell in the DEM, the elevation difference between it and its neighboring cells can be calculated, and then the slope of that cell can be calculated based on these differences.

[0042] Suppose a pixel A in a digital elevation model has eight neighboring pixels B1-B8. In this embodiment, the elevation difference between pixel A and each of its neighboring pixels can be calculated separately. For example, if the elevation of pixel A is 100 meters and the elevation of pixel B1 is 102 meters, then their elevation difference is 102-100=2 meters.

[0043] Then, the slope of pixel A is calculated based on these elevation differences. The slope calculation can be performed through the following steps: First, calculate the elevation change rate of pixel A in both the horizontal and vertical directions. The horizontal elevation change rate can be obtained by calculating the elevation difference between pixel A and its horizontally adjacent pixels, and the vertical elevation change rate can be obtained by calculating the elevation difference between pixel A and its vertically adjacent pixels. Then, based on the horizontal and vertical elevation change rates, the slope of pixel A is calculated using trigonometric functions.

[0044] For example, assuming that the elevation change rate of pixel A in the horizontal direction is 0.1 and the elevation change rate in the vertical direction is 0.2, the slope of pixel A can be calculated by the arctangent function: slope = arctan(√(0.1² + 0.2²)) ≈ 11.3°.

[0045] The elevation change rate threshold for terrain slope characteristics can be dynamically adjusted according to preset terrain classification rules. These preset rules can be set based on different terrain types and application requirements. For example, for mountainous terrain, the elevation change rate threshold can be set to a higher value to distinguish between steep slopes and relatively gentle valleys; for plains terrain, the threshold can be set to a lower value to more accurately reflect subtle terrain undulations.

[0046] Step S124: Extract the reflectance ratio of the shortwave infrared band to the visible light band in the standardized remote sensing image, and combine it with the regional soil reflectance baseline in the soil type database to determine the soil exposure characteristics, wherein the degree of soil exposure is negatively correlated with the reflectance ratio.

[0047] To determine the characteristics of soil exposure, the reflectance ratio of the shortwave infrared band to the visible light band in standardized remote sensing images can be extracted and analyzed in conjunction with the regional soil reflectance baseline in the soil type database.

[0048] Soil exhibits different reflectance characteristics across different spectral bands. The shortwave infrared band is more sensitive to soil moisture content and mineral composition, while the visible light band reflects information such as soil color and texture. By calculating the reflectance ratio of the shortwave infrared band to the visible light band, a numerical value reflecting soil characteristics can be obtained.

[0049] Taking standardized remote sensing imagery as an example, the reflectance ratio can be calculated using the red band within the shortwave infrared band and the visible light band. Assuming the reflectance of a pixel in the shortwave infrared band is 0.3 and the reflectance of the same pixel in the red band is 0.2, then the reflectance ratio of that pixel is 0.3 / 0.2 = 1.5.

[0050] The soil type database contains baseline information on reflectance for different regions and soil types. This baseline information was obtained through extensive field measurements and data analysis, and reflects the typical reflectance characteristics of different soil types.

[0051] By combining the regional soil reflectance baseline in the soil type database, the degree of soil exposure for each pixel can be determined. If the reflectance ratio of a pixel is close to or higher than the regional soil reflectance baseline, it indicates that the area corresponding to that pixel may have a high degree of soil exposure; if the reflectance ratio is much lower than the regional soil reflectance baseline, it indicates that the area may be covered by vegetation and the degree of soil exposure is low.

[0052] For example, assuming the baseline soil reflectance in the soil type database for this area is 1.2, and the reflectance ratio of a certain pixel is 1.5, since 1.5 is greater than 1.2, it indicates that the area corresponding to this pixel may have a certain degree of soil exposure. By performing the above analysis on all pixels of the entire image, the soil exposure characteristics of each pixel can be obtained, thereby determining the soil exposure status of the entire target area.

[0053] Step S130: Construct a soil and water loss prediction model based on spatiotemporal correlation, input the vegetation coverage characteristics, terrain slope characteristics and soil exposure characteristics into the soil and water loss prediction model to predict soil and water loss, and generate a distribution map of soil and water loss levels in the target area.

[0054] After obtaining the vegetation cover features, topographic slope features, and soil exposure features of each remote sensing image subset, the next step is to construct a soil erosion prediction model based on spatiotemporal correlation, and input these features into the model to predict soil erosion in order to generate a distribution map of soil erosion levels in the target area.

[0055] Step S210: Construct a soil erosion prediction model based on spatiotemporal correlation, including:

[0056] Step S211: Obtain a historical soil erosion labeling dataset, which includes multiple historical remote sensing image samples and their corresponding erosion level labels. The erosion level labels are determined based on field survey data and soil and water conservation indicators.

[0057] To construct a soil erosion prediction model based on spatiotemporal correlation, it is first necessary to obtain a historical soil erosion labeled dataset. This historical soil erosion labeled dataset is the foundation for model training, and it contains multiple historical remote sensing image samples and their corresponding erosion level labels.

[0058] Historical remote sensing image samples can be obtained from historical remote sensing image data resources, similar to the method for obtaining multi-temporal remote sensing image datasets in step S110. For example, in this embodiment, remote sensing images from the past 10-20 years can be selected as historical remote sensing image samples from Landsat series satellite data.

[0059] The soil erosion level label is determined based on a combination of field survey data and soil and water conservation indicators. Field survey data can be obtained through methods such as manual field investigation and drone mapping, and it can accurately reflect the actual situation of soil erosion in the target area. Soil and water conservation indicators include multiple factors such as soil erosion modulus, vegetation cover, and topographic slope, which can be obtained through relevant monitoring equipment and calculation methods.

[0060] For example, by conducting on-site surveys of the target area, data such as soil erosion modulus, vegetation coverage, and terrain slope of different areas can be obtained, and the soil and water loss situation can then be divided into four levels: mild, moderate, severe, and extremely severe. For example, when the soil erosion modulus is less than 500 tons / (km²·year), the vegetation coverage is greater than 80%, and the terrain slope is less than 10°, the soil erosion level of the area can be marked as mild; when the soil erosion modulus is between 500 and 2500 tons / (km²·year), the vegetation coverage is between 60% and 80%, and the terrain slope is between 10° and 25°, the soil erosion level of the area can be marked as moderate; when the soil erosion modulus is between 2500 and 5000 tons / (km²·year), the vegetation coverage is between 30% and 60%, and the terrain slope is between 25° and 35°, the soil erosion level of the area can be marked as severe; when the soil erosion modulus is greater than 5000 tons / (km²·year), the vegetation coverage is less than 30%, and the terrain slope is greater than 35°, the soil erosion level of the area can be marked as extremely severe.

[0061] By conducting on-site surveys and analyzing soil and water conservation indicators on multiple historical remote sensing image samples, each sample can be labeled with a corresponding level of soil erosion, thus obtaining a labeled dataset of soil erosion from historical periods.

[0062] Step S212: Perform the surface feature extraction process on each historical remote sensing image sample to obtain historical vegetation cover features, historical topographic slope features, and historical soil exposure features.

[0063] In this embodiment, after obtaining the historical soil erosion labeled dataset, it is necessary to perform surface feature extraction processing on each historical remote sensing image sample to obtain historical vegetation cover features, historical topographic slope features, and historical soil exposure features.

[0064] The method for surface feature extraction is the same as that used in steps S120-S124 above for surface feature extraction from the multi-temporal remote sensing image dataset. First, radiometric correction is performed on the historical remote sensing image samples to obtain standardized historical remote sensing images. Then, a pre-trained vegetation index calculation model is used to calculate historical vegetation cover characteristics based on the reflectance difference between the near-infrared and red bands in the standardized historical remote sensing images. Next, digital elevation model inversion processing is performed on the standardized historical remote sensing images, and historical topographic slope characteristics are calculated based on the elevation change rate between adjacent pixels. Finally, the reflectance ratio of the shortwave infrared band to the visible light band in the standardized historical remote sensing images is extracted, and combined with the regional soil reflectance baseline in the soil type database, the historical soil exposure characteristics are determined.

[0065] For example, for a historical remote sensing image sample, through the above processing steps, this embodiment can obtain the historical vegetation cover features (such as NDVI value), historical topographic slope features (such as slope angle), and historical soil bareness features (such as reflectance ratio) for each pixel in the sample. These features will be used as input data for model training.

[0066] Step S213: Arrange the historical vegetation cover characteristics, historical topographic slope characteristics, and historical soil exposure characteristics in chronological order to generate a spatiotemporal feature sequence.

[0067] After obtaining the historical vegetation cover characteristics, historical topographic slope characteristics, and historical soil exposure characteristics of each historical remote sensing image sample, these characteristics need to be arranged in chronological order to generate a spatiotemporal feature sequence.

[0068] Suppose we have 10 historical remote sensing image samples, each corresponding to a different time point in the past 10 years. For each historical remote sensing image sample, we have obtained its corresponding historical vegetation cover characteristics, historical topographic slope characteristics, and historical soil exposure characteristics. These characteristics can then be arranged in chronological order to form a three-dimensional spatiotemporal feature sequence.

[0069] For example, for the first historical remote sensing image sample (corresponding to year 1), this embodiment obtains a feature vector composed of vegetation cover features, terrain slope features, and soil bareness features. For the second historical remote sensing image sample (corresponding to year 2), this embodiment also obtains a corresponding feature vector. By arranging these feature vectors in chronological order, a spatiotemporal feature sequence containing 10 feature vectors can be obtained.

[0070] Spatiotemporal feature sequences can reflect the changes in surface features of a target area over time, which is of great significance for soil erosion prediction. For example, trends such as decreased vegetation cover, increased topographic slope, and increased soil exposure may all indicate an increased risk of soil erosion. By analyzing spatiotemporal feature sequences, these changes can be captured.

[0071] Step S214: Construct a bidirectional temporal convolutional network model, input the spatiotemporal feature sequence into the bidirectional temporal convolutional network model for forward and backward feature propagation, and capture the dependency relationship of features changing over time.

[0072] To better process spatiotemporal feature sequences and uncover the dependencies between features over time, it is necessary to construct a bidirectional temporal convolutional network model. This model combines forward and backward information processing capabilities, enabling it to capture information from time series more comprehensively.

[0073] The specific steps for constructing a bidirectional temporal convolutional network model are as follows:

[0074] Step S2141: Set the forward convolutional layer to perform sliding window convolution operation on the spatiotemporal feature sequence in forward temporal order, and output a forward feature map containing forward temporal dependencies; set the backward convolutional layer to perform sliding window convolution operation on the spatiotemporal feature sequence in reverse temporal order, and output a backward feature map containing inverse temporal dependencies.

[0075] In this embodiment, the forward convolutional layer processes the spatiotemporal feature sequence in a forward temporal order. Taking a spatiotemporal feature sequence containing 10 time steps as an example, the forward convolutional layer starts from the first time step and performs convolution operations on the feature vectors of each time step sequentially. The convolution operation uses a sliding window, which slides across the spatiotemporal feature sequence, covering a set number of time steps and feature dimensions each time.

[0076] Assuming the sliding window size is 3 time steps, during the forward convolution process, when the window is at time steps 1-3, convolution is performed on the feature vectors of these 3 time steps to obtain a new feature vector containing the positive temporal dependencies of time steps 1-3. Then the window slides forward one time step, covering time steps 2-4, and convolution is performed again to obtain a new feature vector. This process continues until the window has traversed the entire spatiotemporal feature sequence, ultimately outputting a positive feature map containing the positive temporal dependencies.

[0077] The backward convolutional layer processes the spatiotemporal feature sequence in reverse temporal order. Taking a spatiotemporal feature sequence with 10 time steps as an example, the backward convolutional layer starts from the 10th time step and performs convolution operations on the feature vectors of each time step sequentially. The sliding window size and convolution calculation method are the same as the forward convolutional layer, only the direction is reversed. Through the backward convolution operation, the output is a backward feature map containing inverse temporal dependencies.

[0078] For example, an eigenvector in a forward feature map may reflect a trend of gradually decreasing vegetation cover over a past period, while an eigenvector in a reverse feature map may reflect a trend of potentially increasing soil bareness over a future period.

[0079] Step S2142: Set a depth convolution module at the output of the forward convolutional layer and the backward convolutional layer, and perform single-channel spatial convolution operation on the forward feature map and the backward feature map respectively to generate a depth convolution feature map.

[0080] In this embodiment, the depthwise convolution module further extracts and enhances features from the forward and inverse feature maps. Single-channel spatial convolution refers to performing convolution calculations independently on each channel, without merging channels.

[0081] For the forward feature map, the depthwise convolution module performs convolution operations on the feature map of each channel. Assuming the forward feature map has 10 channels, and each channel's feature map size is 100×100 pixels, the depthwise convolution module uses a convolution kernel (e.g., a 3×3 kernel) to perform sliding convolution on the feature map of each channel, calculating new feature values ​​through convolution, thus generating the depthwise convolutional feature map. Similarly, the same single-channel spatial convolution operation is performed on the inverse feature map to obtain the corresponding depthwise convolutional feature map.

[0082] Deep convolution operations can extract more local and subtle feature information, enhancing the expressive power of feature maps. For example, there may be some local features about changes in vegetation cover in the forward feature map, and deep convolution operations can further highlight these features.

[0083] Step S2143: Set a pointwise convolution module at the output of the depth convolution module, and perform a channel fusion operation with a 1×1 convolution kernel on each depth convolution feature map to generate a dimension-reduced fused feature map.

[0084] In this embodiment, the pointwise convolution module uses a 1×1 convolution kernel to perform channel fusion on the depthwise convolution feature map. The function of the 1×1 convolution kernel is to linearly combine the channels without changing the spatial size of the feature map, thereby achieving dimensionality reduction of the number of channels.

[0085] Assuming a depthwise convolutional feature map has 50 channels, a pointwise convolution module using a 1×1 kernel can be used to perform convolution calculations, fusing the feature information from these 50 channels to obtain a fused feature map with fewer channels (e.g., 20 channels). This reduces the dimensionality of the feature map, lowers computational complexity, and retains important feature information.

[0086] For example, multiple channels in the original deep convolutional feature map may contain some redundant feature information. By performing channel fusion operation through pointwise convolution, this redundant information can be integrated to obtain a more compact and efficient fused feature map.

[0087] Step S2144: Set a channel attention module at the output of the pointwise convolution module to concatenate the forward fusion feature map and the inverse fusion feature map along the channel dimension to generate a concatenated feature map.

[0088] The function of the channel attention module is to allocate attention to the channels of the feature map, highlighting important channel information. First, the forward fused feature map and the inverse fused feature map are concatenated along the channel dimension.

[0089] Assuming the forward fusion feature map has 20 channels and the backward fusion feature map also has 20 channels, concatenating them along the channel dimension results in a concatenated feature map with 40 channels. This concatenation operation integrates the forward and backward feature information, allowing the model to simultaneously consider both forward and backward temporal dependencies.

[0090] Step S2145: Calculate the weight coefficients of each channel in the stitched feature map using the channel attention module, and perform a weighted summation operation on the stitched feature map according to the weight coefficients to generate an attention-weighted interactive feature map.

[0091] In this embodiment, the channel attention module can calculate the weight coefficient of each channel in the concatenated feature map. The specific calculation process is as follows: First, a global average pooling operation is performed on the concatenated feature map to compress the feature map of each channel into a scalar value, resulting in a vector of one channel dimension. Then, a fully connected layer is used to perform a non-linear transformation on this vector to obtain a new vector. Next, the new vector is processed by the Sigmoid activation function, mapping each element in the vector to the range of 0-1, thus obtaining the weight coefficient of each channel.

[0092] For example, a concatenated feature map has 40 channels. After processing with global average pooling and a fully connected layer, a vector containing 40 elements is obtained. Then, after processing with the Sigmoid activation function, 40 weight coefficients are obtained, with each weight coefficient corresponding to one channel.

[0093] Finally, a weighted summation operation is performed on the stitched feature maps based on these weight coefficients. For each channel in the stitched feature map, its feature map is multiplied by the corresponding weight coefficient, and then the weighted feature maps of all channels are summed to obtain an attention-weighted interaction feature map. This highlights the feature information of important channels and suppresses the information of unimportant channels.

[0094] Step S2146: Set a multi-scale pooling module at the output of the channel attention module, and perform max pooling operations of different scales on the interactive feature map to generate a multi-scale pooled feature map.

[0095] In this embodiment, the purpose of the multi-scale pooling module is to extract feature information at different scales, and to perform max pooling operations at different scales on the attention-weighted interactive feature map respectively.

[0096] For example, three pooling windows of different sizes are set: 2×2, 4×4, and 8×8. For the 2×2 pooling window, it slides across the interactive feature map, taking the maximum value within the window each time as the new feature value, resulting in a feature map after 2×2 max pooling. Similarly, max pooling operations are performed using 4×4 and 8×8 pooling windows to obtain the corresponding feature maps.

[0097] Pooling operations at different scales can capture feature information of varying sizes. A 2×2 pooling window can extract more local features, while an 8×8 pooling window can extract more global features. Through multi-scale pooling operations, multi-scale pooled feature maps containing feature information at different scales can be obtained.

[0098] Step S2147: Concatenate the multi-scale pooling feature maps according to the channel dimension to generate a feature vector containing multi-scale spatiotemporal features.

[0099] The multi-scale pooling feature maps obtained by max pooling at different scales are concatenated along the channel dimension. Assuming that the feature maps obtained by 2×2, 4×4 and 8×8 max pooling have 10, 15 and 20 channels respectively, after concatenating them along the channel dimension, the resulting feature vector has 10+15+20=45 channels.

[0100] The aforementioned feature vectors, which incorporate multi-scale spatiotemporal features, integrate feature information from different scales and time dependencies, enabling a more comprehensive description of soil erosion in the target area.

[0101] In this embodiment, step S215: supervise the training of the bidirectional temporal convolutional network model based on the erosion level label until the error between the predicted erosion level output by the bidirectional temporal convolutional network model and the erosion level label is lower than a preset threshold, thereby obtaining the trained soil and water loss prediction model.

[0102] After constructing the bidirectional temporal convolutional network model, it is necessary to use a historical soil erosion labeled dataset to conduct supervised training on the model in order to adjust the model's parameters so that it can accurately predict the level of soil erosion.

[0103] Step S2151: Divide the spatiotemporal feature sequence into continuous time segments according to time order, randomly select some time segments as training set, and use the remaining time segments as validation set. Each time segment contains a preset number of historical remote sensing image samples.

[0104] In this embodiment, assuming the spatiotemporal feature sequence contains 100 time steps, it is divided into time segments of 10 time steps each, resulting in 10 time segments. Each time segment contains feature information corresponding to 10 historical remote sensing image samples.

[0105] Seven time segments are randomly selected as the training set for model training, while the remaining three time segments are used as the validation set to evaluate model performance. This partitioning ensures that the model learns feature information from different time stages during training, while the validation set is used to test the model's generalization ability.

[0106] Step S2152: Input the training set into the forward and backward convolutional layers of the bidirectional temporal convolutional network model to generate the forward and backward feature maps of the training set.

[0107] In this embodiment, the spatiotemporal feature sequence of the training set can be input into the forward and backward convolutional layers of a bidirectional temporal convolutional network model. The forward convolutional layer performs a sliding window convolution operation on the spatiotemporal feature sequence of the training set in forward temporal order, generating a forward feature map containing forward temporal dependencies. The backward convolutional layer performs a sliding window convolution operation on the spatiotemporal feature sequence of the training set in reverse temporal order, generating a backward feature map containing inverse temporal dependencies.

[0108] For example, if the training set contains 7 time segments, each with 10 time steps, the forward convolutional layer will sequentially convolve the feature vectors at each time step, ultimately generating 7 forward feature maps, each corresponding to one time segment. Similarly, the backward convolutional layer will generate 7 backward feature maps.

[0109] Step S2153: Process the forward and inverse feature maps through the depthwise convolution module and the pointwise convolution module to generate the fused feature map of the training set. Process the fused feature map through the channel attention module and the multi-scale pooling module to generate the multi-scale feature vector of the training set. Input the multi-scale feature vector into the fully connected layer and output the predicted churn level probability distribution of the training set.

[0110] In this embodiment, the forward and inverse feature maps undergo single-channel spatial convolution operations via a depthwise convolution module to obtain a depthwise convolution feature map. Then, a pointwise convolution module performs a channel fusion operation with a 1×1 convolution kernel on the depthwise convolution feature map to generate a fused feature map for the training set.

[0111] The channel attention module concatenates the fused feature maps along the channel dimension, calculates the weight coefficients of each channel, performs a weighted summation operation, and generates an attention-weighted interactive feature map. The multi-scale pooling module performs max pooling operations at different scales on the interactive feature maps to generate multi-scale pooled feature maps, and concatenates them along the channel dimension to obtain the multi-scale feature vector of the training set.

[0112] Finally, the multi-scale feature vectors are input into the fully connected layer. The fully connected layer performs linear transformations and softmax activation on the multi-scale feature vectors, outputting the probability distribution of predicted soil erosion levels for the training set. Assuming that soil erosion levels are divided into four levels—slight, moderate, severe, and extremely severe—the fully connected layer will output a probability distribution vector containing four elements, where each element represents the probability of predicting the corresponding level.

[0113] Step S2154: Based on the true category index of the churn level label, extract the probability value of the corresponding category in the predicted churn level probability distribution, calculate the difference between the probability value and 1, and generate the prediction error coefficient for each training sample.

[0114] In this embodiment, for each sample in the training set, its churn level label corresponds to a true category index. For example, a mild churn level corresponds to index 0, a moderate churn level corresponds to index 1, a severe churn level corresponds to index 2, and a very severe churn level corresponds to index 3.

[0115] Based on this true category index, the probability value of the corresponding category is extracted from the predicted churn level probability distribution. Assuming a sample's true churn level is moderate, its true category index is 1, and the predicted churn level probability distribution is [0.1, 0.3, 0.4, 0.2], then the extracted corresponding category probability value is 0.3.

[0116] Then, the difference between this probability value and 1 is calculated to obtain the prediction error coefficient for the training sample. In the example above, the prediction error coefficient is 1 - 0.3 = 0.7.

[0117] Step S2155: Perform an exponential operation on the prediction error coefficient and multiply it by a preset focusing coefficient to generate a sample-level weight adjustment factor. Then, multiply the weight adjustment factor by the cross-entropy loss to generate the focus loss value for each training sample.

[0118] In this embodiment, the prediction error coefficient is exponentially calculated. Assuming the exponent is 2, the prediction error coefficient of 0.7 in the above example becomes 0.49 after exponential calculation.

[0119] The preset focus coefficient is a hyperparameter used to adjust the degree of attention given to samples of different difficulty. Assuming the focus coefficient is 2, then the sample-level weight adjustment factor is 0.49 × 2 = 0.98.

[0120] Cross-entropy loss is a loss function that measures the difference between the predicted probability distribution and the true label. For each training sample, its cross-entropy loss is calculated. Assuming the cross-entropy loss of the above sample is 0.5, multiplying the weight adjustment factor by the cross-entropy loss yields a focal loss value of 0.98 × 0.5 = 0.49 for this training sample.

[0121] The focus loss value can effectively solve the problem of sample imbalance by giving higher weights to samples with large prediction errors, making the model pay more attention to these difficult-to-classify samples.

[0122] Step S2156: Set the category weight coefficient according to the severity of the churn level label, multiply the focus loss value by the category weight coefficient to generate a weighted focus loss value, and perform batch averaging on the weighted focus loss value to generate the total loss value of the current training batch.

[0123] In this embodiment, a category weight coefficient is set according to the severity of the churn level label. For example, the category weight coefficient is set to 1 for the mild churn level; 1.5 for the moderate churn level; 2 for the severe churn level; and 2.5 for the extremely severe churn level.

[0124] For each training sample, its focus loss value is multiplied by the corresponding class weight coefficient to obtain the weighted focus loss value. Assuming the focus loss value of the above sample is 0.49, its true churn level is moderate, and the class weight coefficient is 1.5, then the weighted focus loss value is 0.49 × 1.5 = 0.735.

[0125] In a training batch, there are multiple training samples. The weighted focus loss values ​​of these samples are averaged to obtain the total loss value of the current training batch. Assuming a training batch has 10 samples, the weighted focus loss values ​​of these 10 samples are added together and then divided by 10 to obtain the total loss value.

[0126] Step S2157: Calculate the gradient of the total loss value with respect to the model parameters during backpropagation, detect the norm of the gradient vector, and when the norm of the gradient vector exceeds a preset clipping threshold, scale the gradient vector proportionally until the norm is equal to the clipping threshold, and use the scaled gradient vector to update the model parameters to complete one training iteration.

[0127] In this embodiment, during backpropagation, the gradient of the total loss value with respect to each parameter in the model is calculated. The gradient represents the direction and rate of change of the total loss value with respect to the parameters.

[0128] The norm of the gradient vector is detected; the norm reflects the magnitude of the gradient. A preset clipping threshold is a pre-defined value used to control the magnitude of the gradient and prevent gradient explosion.

[0129] When the norm of the gradient vector exceeds a preset clipping threshold, the gradient vector is scaled proportionally until the norm equals the clipping threshold. For example, if the preset clipping threshold is 5 and the norm of the gradient vector is 8, then each element of the gradient vector is multiplied by 5 / 8 to obtain the scaled gradient vector.

[0130] The model parameters are updated using the scaled gradient vector. For example, for a parameter W in the model, the update formula is W = W - learning rate × scaled gradient. The learning rate is a hyperparameter that controls the step size of parameter updates. One training iteration is completed in this way, continuously adjusting the model parameters to gradually reduce the total loss.

[0131] Step S2158: Repeat the forward propagation process on the validation set, calculate the prediction accuracy of the validation set, and stop training and save the current model parameters when the prediction accuracy of the validation set does not improve after N consecutive training iterations.

[0132] In this embodiment, the forward propagation process is repeated on the validation set after each training iteration. The spatiotemporal feature sequence of the validation set is input into the trained bidirectional temporal convolutional network model, which passes through forward convolutional layers, backward convolutional layers, depthwise convolutional modules, pointwise convolutional modules, channel attention modules, multi-scale pooling modules, and fully connected layers, and outputs the predicted churn level probability distribution of the validation set.

[0133] Based on the predicted churn level probability distribution and the actual churn level labels in the validation set, calculate the prediction accuracy of the validation set. Prediction accuracy refers to the proportion of correctly predicted samples out of the total number of samples.

[0134] A threshold N is set. When the prediction accuracy on the validation set does not improve after N consecutive training iterations, it indicates that the model has reached a relatively stable state, and continued training may lead to overfitting. At this point, training is stopped, and the parameters of the current model are saved, resulting in the trained soil erosion prediction model.

[0135] Step S131: Perform differential calculations on the vegetation cover characteristics, topographic slope characteristics, and soil bareness characteristics of the current time period with the corresponding characteristics of historical time periods to obtain the changes in vegetation cover, topographic slope, and soil bareness.

[0136] In this embodiment, after obtaining the trained soil and water loss prediction model, the vegetation coverage characteristics, terrain slope characteristics, and soil exposure characteristics of the current time period are compared with the corresponding characteristics of historical time periods using differential calculation.

[0137] Assume the current time period is year 11 and the historical time period is year 10. Taking NDVI value as an example, for vegetation cover characteristics, if the NDVI value of a certain pixel is 0.6 in year 10 and 0.5 in year 11, then the change in vegetation cover for that pixel is 0.5-0.6=-0.1, indicating that the vegetation cover has decreased.

[0138] Regarding topographic slope characteristics, if the slope of a pixel is 15° in year 10 and 16° in year 11, then the change in topographic slope for that pixel is 16° - 15° = 1°, indicating an increase in topographic slope. Regarding soil exposure characteristics, taking the reflectance ratio of the shortwave infrared band to the visible light band as an example, if the reflectance ratio of a pixel is 1.2 in year 10 and 1.3 in year 11, then the change in soil exposure for that pixel is 1.3 - 1.2 = 0.1, meaning an increase in soil exposure. By performing this difference calculation on all pixels of the entire image, we can obtain the changes in vegetation cover, topographic slope, and soil exposure for each pixel, thus obtaining feature maps of vegetation cover change, topographic slope change, and soil exposure change for the target area. The aforementioned change feature map can intuitively reflect the changes in surface features of the target area in the current time period relative to historical time periods, providing key input information for soil and water loss prediction.

[0139] Step S132: Input the vegetation cover change, topographic slope change, and soil exposure change into the bidirectional temporal convolutional network model of the soil and water loss prediction model to generate a loss probability score for each pixel, wherein the loss probability score reflects the increase in loss risk in the current time period relative to historical periods.

[0140] In this embodiment, after obtaining the feature maps of vegetation cover change, topographic slope change, and soil exposure change, these feature maps are used as input data and fed into the bidirectional temporal convolutional network model of the trained soil erosion prediction model. The bidirectional temporal convolutional network model processes the input data according to its internal structure and parameters. First, the forward convolutional layer performs a sliding window convolution operation on the input feature maps in forward temporal order to capture positive temporal dependencies and generate positive feature maps. For example, for the vegetation cover change feature map, the forward convolutional layer analyzes the trend information of vegetation cover change from the past to the present. The backward convolutional layer then performs a convolution operation on the input feature maps in reverse temporal order to generate inverse feature maps, capturing inverse temporal dependencies, such as predicting possible future trends.

[0141] Next, the deep convolution module performs single-channel spatial convolution operations on both the forward and inverse feature maps to further extract subtle local features and generate deep convolution feature maps. The pointwise convolution module performs channel fusion operations with 1×1 kernels on the deep convolution feature maps to reduce their dimensionality and generate a dimensionality-reduced fused feature map. The channel attention module concatenates the forward and inverse fused feature maps along their channel dimensions to generate a concatenated feature map. It then calculates the weight coefficients of each channel in the concatenated feature map, performs a weighted summation operation, and generates an attention-weighted interactive feature map that highlights the features of important channels.

[0142] The multi-scale pooling module performs max pooling operations at different scales on the interactive feature maps, generating multi-scale pooled feature maps. It extracts feature information of different sizes and then concatenates these multi-scale pooled feature maps along the channel dimension to generate a feature vector containing multi-scale spatiotemporal features. Finally, this feature vector is input into a fully connected layer, which performs a linear transformation and a Softmax activation function on it, outputting a churn probability score for each pixel.

[0143] Assuming the erosion probability score ranges from 0 to 1, the closer the score is to 1, the greater the increase in erosion risk for that pixel in the current time period compared to historical periods. For example, a erosion probability score of 0.8 indicates a significant increase in the risk of soil erosion for that pixel in the current time period compared to historical periods; while a score of 0.2 indicates a relatively smaller increase in erosion risk. By processing all pixels in the entire image in this way, a erosion probability score can be obtained for each pixel, forming a erosion probability score map. This map clearly shows the increase in soil erosion risk at various locations within the target area.

[0144] Step S133: Map the churn probability score to discrete churn levels according to the preset churn level classification rules.

[0145] In this embodiment, the preset erosion level classification rule is based on a large amount of actual data and experience, and is used to convert continuous erosion probability scores into discrete erosion levels. Assuming that soil erosion levels are divided into four levels—mild, moderate, severe, and extremely severe—the corresponding erosion probability score ranges are set as follows: 0-0.2 for mild erosion, 0.2-0.5 for moderate erosion, 0.5-0.8 for severe erosion, and 0.8-1 for extremely severe erosion.

[0146] For each pixel in the erosion probability scoring map, it is mapped to the corresponding erosion level based on the interval in which its erosion probability score falls. For example, a pixel with a erosion probability score of 0.15 falls within the range of 0-0.2, so its corresponding erosion level is mild. Another pixel with a erosion probability score of 0.6 falls within the range of 0.5-0.8, so its corresponding erosion level is severe. By performing this mapping operation on all pixels in the entire erosion probability scoring map, the erosion level of each pixel can be obtained, forming a erosion level map. This map visually displays the soil and water erosion situation at various locations in the target area in a discrete level format.

[0147] Step S134: Associate the erosion level with the corresponding geographic coordinates to generate a water and soil erosion level distribution map with spatial resolution, wherein the spatial resolution is consistent with the resolution of the remote sensing image subset.

[0148] In this embodiment, after obtaining the erosion level map, it is necessary to associate the erosion level of each pixel with its corresponding geographic coordinates to generate a spatially resolved soil erosion level distribution map. A subset of remote sensing images has a certain spatial resolution, which represents the actual ground area represented by each pixel in the image. For example, a subset of remote sensing images with a spatial resolution of 30 meters × 30 meters means that each pixel in the image corresponds to a 30-meter × 30-meter area on the ground.

[0149] Each pixel in the image has its corresponding row and column number. Using the image's georeferenced information (such as projection method, coordinate origin, pixel size, etc.), the pixel's row and column number can be converted into actual geographic coordinates. For each pixel in the loss level map, its loss level is associated with its corresponding geographic coordinates. For example, a pixel located in the 10th row and 20th column of the image has a moderate loss level. Using georeferenced information, the pixel's row and column number can be converted into geographic coordinates (110.123°E, 30.456°N), and then the moderate loss level is associated with these geographic coordinates.

[0150] By integrating the erosion levels of all pixels with their corresponding geographic coordinates, a spatially resolved distribution map of soil erosion levels can be generated. This map has the same spatial resolution as a subset of remote sensing imagery, accurately displaying the soil erosion levels at various locations within the target area, providing a clear basis for subsequent soil erosion control and management.

[0151] Step S140: Identify target erosion risk areas based on the soil erosion level distribution map, and generate a priority ranking of governance and vegetation restoration strategies for the target erosion risk areas.

[0152] Step S141: Identify the target erosion risk area based on the soil erosion level distribution map, including:

[0153] Step S1411: Perform spatial cluster analysis on the soil erosion level distribution map and extract continuously distributed areas with severe erosion levels as candidate erosion areas.

[0154] In this embodiment, spatial clustering analysis is an analytical method that groups spatially adjacent objects with similar attributes together. For a soil erosion level distribution map, a spatial clustering algorithm (such as the DBSCAN algorithm) is used to perform clustering analysis on pixels of severe erosion level. The DBSCAN algorithm performs clustering based on the spatial distance and density between pixels.

[0155] First, set a neighborhood radius and a minimum point threshold. For each severely eroded cell in the soil erosion grade distribution map, check the number of other severely eroded cells within its neighborhood radius. If this number is greater than or equal to the minimum point threshold, these cells are grouped into a cluster. For example, set the neighborhood radius to a distance of 3 pixels and the minimum point threshold to 5. If a severely eroded cell has 6 other severely eroded cells within a 3-pixel radius, then these cells will be grouped into a cluster.

[0156] In this way, all pixels at the severe erosion level are clustered to extract continuously distributed areas at the severe erosion level as candidate erosion areas. These candidate erosion areas are spatially continuous regions with relatively severe soil erosion, and are areas that require key attention and remediation.

[0157] Step S1412: Obtain land use type data of the candidate lost areas, and filter out areas belonging to cultivated land, forest land or grassland as areas to be verified.

[0158] In this embodiment, to further determine the areas requiring remediation, land use type data of candidate loss areas is obtained. Land use type data can be obtained in various ways, such as high-resolution remote sensing image interpretation and geographic information system (GIS) database queries.

[0159] High-resolution remote sensing image interpretation involves analyzing high-resolution remote sensing images to identify different land use types. For example, deep learning algorithms can be used to classify high-resolution remote sensing images into different types such as cultivated land, forest land, grassland, and construction land. GIS database queries retrieve land use type information from existing geographic information system databases.

[0160] For each candidate soil erosion area, it is overlaid with land use type data for analysis, and areas belonging to cultivated land, forest land, or grassland are selected as areas to be verified. This is because cultivated land, forest land, and grassland are land use types closely related to the ecological environment and soil and water conservation. Soil and water erosion in these areas can have a significant impact on the ecosystem and agricultural production, so it is necessary to further verify their soil and water erosion situation.

[0161] Step S1413: Call the pre-trained erosion gully detection model to calculate the gully density of the remote sensing image of the area to be verified, and mark the area to be verified where the gully density exceeds the first density threshold as the confirmed loss area.

[0162] In this embodiment, step S1413 calls a pre-trained erosion gully detection model to calculate the gully density of the remote sensing image of the area to be verified. The specific detailed steps are as follows:

[0163] Step S14131: Obtain multispectral remote sensing images of the area to be verified, and extract reflectance data in the visible light band and near-infrared band.

[0164] To accurately calculate the gully density of the area to be verified, multispectral remote sensing imagery of that area is first required. Multispectral remote sensing imagery contains information from multiple different spectral bands, among which the visible light and near-infrared bands are crucial for identifying erosion gullies. Multispectral remote sensing imagery of the area to be verified can be acquired through specialized remote sensing satellite data platforms, such as the Landsat series and Sentinel series satellites. Assuming Landsat 8 satellite imagery is used, it contains visible light bands (blue, green, red, etc.) and near-infrared bands. After acquiring the imagery, the image data is analyzed to extract reflectance data for the visible light band (e.g., the red band) and the near-infrared band. Taking a remote sensing imagery of the area to be verified as an example, with an image size of 1000×1000 pixels, the reflectance values ​​in the visible light and near-infrared bands can be extracted for each pixel, forming two 1000×1000 reflectance data matrices, representing the reflectance data for the visible light and near-infrared bands respectively.

[0165] Step S14132: Perform edge enhancement filtering on the reflectivity data of the visible light band to generate an enhanced surface texture feature map.

[0166] The purpose of edge enhancement filtering is to highlight edge information in visible light reflectance data, thereby better identifying surface texture features and providing clearer features for subsequent trench detection. For example, the Sobel operator can be used for edge enhancement filtering. The Sobel operator is a commonly used edge detection operator that detects edges by calculating the gradients of the image in the horizontal and vertical directions. For the extracted visible light reflectance data matrix, the Sobel operator is applied in convolution operations in both the horizontal and vertical directions. Assuming the visible light reflectance data matrix is ​​A, the horizontal Sobel operator is Gx, and the vertical Sobel operator is Gy, the convolution operation yields the horizontal gradient matrix Gx_result and the vertical gradient matrix Gy_result. Then, the gradient magnitude of each pixel is calculated using the Pythagorean theorem; that is, the gradient magnitude of each pixel is equal to the square root of the sum of the squares of the corresponding values ​​in the horizontal and vertical gradient matrices. In this way, a new matrix is ​​obtained, where each pixel value represents the edge intensity of the corresponding pixel in the original visible light reflectance data. Normalizing this matrix to ensure pixel values ​​range from 0 to 255 generates an enhanced surface texture feature map. This feature map more clearly displays surface texture and edge information, aiding in subsequent identification of erosion gullies.

[0167] Step S14133: Perform multi-scale morphological opening operation based on the reflectivity data of the near-infrared band and the enhanced surface texture feature map to eliminate isolated noise points and retain continuous gully structures, thereby obtaining the processed surface texture feature map.

[0168] Multi-scale morphological opening combines near-infrared reflectance data with enhanced surface texture features to effectively eliminate isolated noise points in images while preserving continuous gully structures. Morphological opening includes two operations: erosion and dilation. First, structuring elements of different sizes are selected for multi-scale processing. For example, 3×3, 5×5, and 7×7 structuring elements are chosen. For each structuring element, an erosion operation is performed on the enhanced surface texture feature map. The erosion operation involves sliding the structuring element across the image; if the structuring element is completely contained within the foreground region (i.e., a non-zero pixel region), the pixel corresponding to the center of the structuring element is retained as a foreground pixel; otherwise, it is set as a background pixel (zero pixel). Erosion removes small noise points and fine connections in the image. Then, dilation is performed on the eroded image. Dilation involves sliding the structuring element across the image; if the structuring element intersects with the foreground region, the pixel corresponding to the center of the structuring element is set as a foreground pixel. Dilation can recover the eroded gully structure while preserving continuous gullies. During multi-scale morphological opening operations, for each scale's structuring element, the enhanced surface texture feature map undergoes erosion and dilation operations. The results from different scales are then fused. A simple weighted averaging method can be used for fusion; for example, the results processed at 3×3, 5×5, and 7×7 scales are assigned weights of 0.3, 0.3, and 0.4 respectively, and these weights are summed to obtain the final result. This yields a processed surface texture feature map that eliminates isolated noise points while preserving continuous gully structures, making it more suitable for subsequent gully identification.

[0169] Step S14134: Input the processed surface texture feature map into the pre-trained gully semantic segmentation model and output a binary mask image of the gully pixels.

[0170] The pre-trained gully semantic segmentation model is trained on a large amount of remote sensing image data with gully annotations, and can accurately identify gully pixels in the images. The processed surface texture feature map is input into this gully semantic segmentation model to classify each pixel in the surface texture feature map, determining whether it belongs to a gully pixel. The model's output is a binary mask image, where pixels with a value of 1 represent gully pixels, and pixels with a value of 0 represent non-gully pixels. For example, if the processed surface texture feature map is 1000×1000 pixels in size, after processing by the pre-trained gully semantic segmentation model, the output binary mask image is also a 1000×1000 pixel matrix, where each element is either 0 or 1, clearly identifying the location of gullies in the image.

[0171] Step S14135: Mark the connected components of the binary mask image and count the groove length and width data of each connected component.

[0172] In this embodiment, connected component labeling involves dividing the interconnected groove pixels in the binary mask image into different connected components and assigning a unique label to each connected component. A classic seed-fill algorithm can be used for connected component labeling. Starting from the top left corner of the binary mask image, pixels are traversed row by row and column by column. When an unlabeled groove pixel (pixel value 1) is encountered, it is used as a seed point. From this seed point, the algorithm expands to its neighboring pixels, labeling the adjacent groove pixels with the same label as the seed point. This process is repeated until all interconnected groove pixels are labeled with the same label. In this way, the grooves in the binary mask image are divided into different connected components. For each connected component, its groove length and width are calculated. The groove length can be obtained by calculating the longest path among the pixels in the connected component, while the groove width can be obtained by calculating the average width of the connected component perpendicular to its length. For example, for a connected component, the longest path is found by tracing its boundary pixels, and the number of pixels on the path is taken as the groove length; the width of the connected component is measured at regular intervals in the direction perpendicular to the longest path, and the average value is taken as the groove width.

[0173] Step S14136: Calculate the gully extension density per unit area based on the gully length and width data, and map it to the gully density level using a preset density grading table.

[0174] In this embodiment, the furrow extension density per unit area is an important indicator for measuring the density of furrow distribution. First, the total area of ​​the region to be verified is calculated. Assuming the remote sensing image size of the region to be verified is 1000×1000 pixels, and each pixel represents an actual ground area of ​​1 square meter, then the total area of ​​the region to be verified is 1000×1000=1,000,000 square meters. Then, the furrow lengths of all connected domains are added together to obtain the total furrow extension length. Dividing the total furrow extension length by the total area of ​​the region to be verified yields the furrow extension density per unit area, expressed in meters per square meter (m / m²). A preset density grading table is established based on extensive actual data and experience, classifying furrow extension densities into different levels, such as light, moderate, and heavy. Based on the calculated furrow extension density per unit area, the corresponding level is looked up in the preset density grading table and mapped to a furrow density level. For example, if the calculated furrow extension density is 0.01 m / m², according to the density grading table, this value corresponds to a light furrow density level.

[0175] Step S14137: Overlay the gully density level with the geographical boundary of the candidate loss area to generate a gully density distribution layer with spatial attributes.

[0176] In this embodiment, the geographic boundary information of the candidate lost areas can be obtained from a Geographic Information System (GIS), typically in the form of vector data, such as a polygon vector layer. The calculated gully density level is overlaid with the geographic boundary of the candidate lost areas. In GIS software, the polygon vector layer representing the candidate lost areas can be overlaid with a raster layer containing gully density level information. After overlay, each candidate lost area is assigned corresponding gully density level information, forming a gully density distribution layer with spatial attributes, which can intuitively display the gully density distribution within the candidate lost areas.

[0177] Based on the above, a first density threshold is set, and the gully density of the area to be verified is compared with the first density threshold. If the gully density exceeds the first density threshold, the area to be verified is marked as a confirmed loss area. A confirmed loss area can be understood as an area with more serious soil erosion and which has been further verified. It is an area that needs to be focused on for treatment.

[0178] Step S1414: Perform spatial overlay analysis on the confirmed loss area and the candidate loss area, and remove abnormal areas with terrain abrupt changes or artificial building interference to obtain the final set of target loss risk areas.

[0179] In this embodiment, spatial overlay analysis is performed on the confirmed loss areas and candidate loss areas, that is, the ranges of the confirmed loss areas and the candidate loss areas are analyzed to overlap. During this process, some areas may exhibit abrupt terrain changes or interference from man-made structures.

[0180] Abrupt topographic changes can be caused by drastic alterations in terrain due to natural disasters such as landslides and earthquakes. In such cases, soil erosion may be caused by specific geological events and differs from normal soil erosion. Human-caused disturbances refer to the presence of man-made structures such as buildings and roads within a region, which may affect the natural process of soil erosion.

[0181] For anomalous areas with abrupt topographic changes or artificial disturbances, they are removed from the confirmed loss area. For example, in a confirmed loss area, if a portion of the area has experienced a landslide, a dramatic change in topography, or a highway runs through it, then that portion is considered an anomalous area and removed from the confirmed loss area.

[0182] After such screening and elimination, the final set of target water loss risk areas is obtained. The areas in this set are the areas that truly need water and soil erosion control. The influence of some interfering factors has been eliminated, making the control work more targeted.

[0183] Step S142: Generate a governance priority ranking and vegetation restoration strategy for the target loss risk area, including:

[0184] Step S1421: Obtain historical rainfall data, soil erosion rate data, and vegetation restoration cost data for the target loss risk area.

[0185] In this embodiment, in order to formulate a reasonable governance priority ranking and vegetation restoration strategy, it is necessary to obtain historical rainfall data, soil erosion rate data, and vegetation restoration cost data for the target loss risk area.

[0186] Historical rainfall data can be obtained from meteorological monitoring stations. Meteorological departments have set up multiple monitoring stations in or around the target area, which regularly record rainfall information. By compiling and analyzing the historical rainfall data from these stations, information such as the multi-year average rainfall in the target loss risk area and the distribution of rainfall in different seasons can be obtained. For example, the average annual rainfall in a target loss risk area over the past 10 years may be 800 mm, with rainfall mainly concentrated in the summer.

[0187] Soil erosion rate data can be obtained through two methods: field monitoring and model estimation. Field monitoring involves setting up soil erosion monitoring points in the target loss-risk area and directly measuring the soil erosion rate using methods such as sediment deposition measurement and slope runoff observation. Model estimation, on the other hand, uses soil erosion models (such as the Universal Soil Loss Equation (USLE)) to estimate the soil erosion rate based on factors such as the topography, soil type, and vegetation cover of the target area. For example, through field monitoring and model estimation, the soil erosion rate of a target loss-risk area is estimated to be 3000 tons per square kilometer per year.

[0188] Cost data for vegetation restoration can be obtained through market research and project budget estimation. Market research involves understanding information such as seedling prices, planting and maintenance costs in the vicinity of the target area. Project budget estimation, on the other hand, estimates the total cost required for vegetation restoration based on factors such as the area of ​​the target region and the vegetation restoration plan. For example, through market research and project budget estimation, it is found that the cost of vegetation restoration in a target area at risk of vegetation loss is approximately 50,000 yuan per hectare.

[0189] Step S1422: Calculate the natural recovery potential score based on the soil erosion rate data and rainfall data. The natural recovery potential score is negatively correlated with the erosion rate and positively correlated with the rainfall.

[0190] In this embodiment, the natural restoration potential score is used to assess the ability of a target loss-risk area to restore its ecology under natural conditions. Since a higher soil erosion rate indicates more severe soil erosion and greater difficulty in natural restoration, the natural restoration potential score is negatively correlated with the erosion rate; while more rainfall is more conducive to vegetation growth and restoration, so the natural restoration potential score is positively correlated with rainfall.

[0191] One method for calculating natural recovery potential scores is to first standardize the soil erosion rate and rainfall data, transforming them to the same numerical range (e.g., 0-1). For example, using the min-maximum standardization method, for the soil erosion rate data, find its minimum and maximum values, subtract the minimum value from each data point, and then divide by the difference between the maximum and minimum values ​​to obtain the standardized soil erosion rate data. The same standardization process is applied to the rainfall data.

[0192] Then, the weight of soil erosion rate was set to 0.6, and the weight of rainfall was set to 0.4. For each target loss risk area, the standardized soil erosion rate data was multiplied by 0.6, and the standardized rainfall data was multiplied by 0.4. These two values ​​were then added together to obtain the natural recovery potential score. For example, if the standardized soil erosion rate data for a target loss risk area is 0.8 and the standardized rainfall data is 0.6, then the natural recovery potential score for that area is 0.8 × 0.6 + 0.6 × 0.4 = 0.72. The higher the score, the greater the natural recovery potential of the area.

[0193] Step S1423: Construct a multi-objective optimization function based on the vegetation restoration cost data and natural restoration potential score, and solve for the comprehensive priority index of each target loss risk area.

[0194] In this embodiment, the purpose of constructing a multi-objective optimization function is to comprehensively consider both vegetation restoration costs and natural restoration potential to determine the governance priority of each target area at risk of vegetation loss. The multi-objective optimization function can be expressed as: Comprehensive Priority Index = Natural Restoration Potential Score / Vegetation Restoration Cost.

[0195] For each target area at risk of vegetation loss, its calculated natural restoration potential score is divided by the vegetation restoration cost to obtain the area's comprehensive priority index. For example, if a target area at risk of vegetation loss has a natural restoration potential score of 0.7 and a vegetation restoration cost of 50,000 yuan per hectare, then the comprehensive priority index for that area is 0.7 / 5 = 0.14. The higher the comprehensive priority index, the more priority the area should be given to restoration, considering both its natural restoration potential and vegetation restoration costs.

[0196] Step S1424: Sort the target loss risk areas from high to low according to the comprehensive priority index to generate a governance priority ranking.

[0197] For example, if there are three target areas at risk of soil erosion, with area A having a comprehensive priority index of 0.2, area B having a comprehensive priority index of 0.15, and area C having a comprehensive priority index of 0.1, then the priority ranking for soil erosion control is: area A > area B > area C. This ranking provides a clear order for soil erosion control work, prioritizing the control of areas with higher comprehensive priority indices, which can more effectively utilize resources and improve the effectiveness of control efforts.

[0198] Step S1425: Based on the terrain slope characteristics and soil type of the target loss risk area, match the preset vegetation restoration scheme library to generate a vegetation restoration strategy that includes planting density and maintenance cycle.

[0199] For example, the specific steps of this process are as follows:

[0200] Step S14251: Extract the digital elevation model data of the target loss risk area and calculate the slope aspect distribution map based on the elevation gradient.

[0201] Digital Elevation Model (DEM) data can be acquired through various methods, such as lidar measurement and Interferometric Synthetic Aperture Radar (InSAR). After acquiring DEM data for areas at risk of target loss, it is processed to calculate the aspect distribution map. Slope aspect refers to the orientation of a slope, which has a significant impact on vegetation growth and distribution. The method for calculating slope aspect is based on the elevation gradient in the DEM data. For each pixel in the DEM data, its elevation change rate in the horizontal and vertical directions is calculated, i.e., the elevation gradient. Assuming the DEM data is a two-dimensional matrix, each element in the matrix represents the elevation value of a pixel. For a pixel, the elevation change rate in the horizontal and vertical directions is obtained by calculating the elevation difference between it and its neighboring pixels. Then, based on the elevation change rate in the horizontal and vertical directions, the slope aspect angle of that pixel is calculated. The slope aspect angle ranges from 0 to 360 degrees, where 0 degrees represents true north, 90 degrees represents true east, 180 degrees represents true south, and 270 degrees represents true west. By performing this calculation on each pixel in the DEM data, a slope aspect angle matrix is ​​obtained. Visualizing this matrix generates a slope aspect distribution map.

[0202] Step S14252: Spatially overlay the terrain slope characteristics with the slope aspect distribution map to divide the sunny slope area and the shady slope area.

[0203] The terrain slope features, extracted in previous steps, reflect the slope magnitude at each location within the target loss risk area. The terrain slope feature map is then spatially overlaid with the aspect distribution map. In GIS software, these two layers can be overlaid to ensure spatial alignment. Based on the relationship between slope angle and sunlight exposure, areas with slope angles within a certain range (e.g., 0-90 degrees and 270-360 degrees) are classified as sunny slope areas, while areas with slope angles outside this range (e.g., 90-270 degrees) are classified as shady slope areas. Furthermore, combining the terrain slope features, the classification of sunny and shady areas under different slope grades is further refined. For example, sunny areas with slopes less than 10 degrees are marked as low-slope sunny areas, and shady areas with slopes greater than 30 degrees are marked as high-slope shady areas. This spatial overlay and classification allows for a more accurate understanding of the sunlight and slope conditions in different areas within the target loss risk area, providing more detailed information for subsequent vegetation restoration strategy development.

[0204] Step S14253: Obtain soil water-holding capacity parameters from the regional soil type database, and generate a soil-slope water-holding capacity matrix by combining the slope regional type.

[0205] The regional soil type database contains soil type information and corresponding soil water-holding capacity parameters for the target erosion risk area and its surrounding areas. Soil water-holding capacity refers to the soil's ability to retain moisture, which is related to factors such as soil texture and structure. Water-holding capacity parameters for different soil types within the target erosion risk area are obtained from the regional soil type database. Then, combined with previously defined slope area types (such as low-slope sunny areas, high-slope shady areas, etc.), a soil-slope water-holding capacity matrix is ​​generated. The rows of this matrix represent different soil types, the columns represent different slope area types, and each element in the matrix represents the soil water-holding capacity under the corresponding combination of soil type and slope area type. For example, for a certain soil type in a low-slope sunny area, its corresponding soil water-holding capacity is a specific numerical value, which is filled into the corresponding position in the matrix. By generating the soil-slope water-holding capacity matrix, the soil water-holding capacity under different combinations of soil types and slope areas can be clearly displayed, providing an important basis for selecting suitable vegetation.

[0206] Step S14254: Query the list of suitable plant species in the vegetation restoration scheme library according to the soil-slope water holding capacity matrix, and select plant species.

[0207] The vegetation restoration scheme library contains a list of suitable plant species for different soil conditions, slope conditions, and light conditions. Based on the generated soil-slope water-holding capacity matrix, a search is performed in the vegetation restoration scheme library. For each element in the matrix—that is, each combination of soil type and slope region type—a list of plant species suitable for that condition is retrieved from the vegetation restoration scheme library. For example, for low-slope, sunny areas with strong soil water-holding capacity, suitable plant species, such as some drought-tolerant and sun-loving herbaceous plants and shrubs, are found from the vegetation restoration scheme library. All plant species that meet the criteria are summarized, and plant species suitable for different areas of the target loss-risk zone are selected. During the selection process, factors such as the plant's ecological adaptability, growth rate, and stress resistance can also be considered to ensure that the selected plant species can grow well in the target area.

[0208] Step S14255: Based on the growth cycle data of the selected plant species and the accumulated temperature data of the target loss risk area, calculate the optimal planting time window for each plant species in different slope areas.

[0209] The selected plant species all have their own growth cycle data, including germination, growth, flowering, and fruiting stages. Accumulated temperature data for the target loss-risk area can be obtained from meteorological monitoring data. Accumulated temperature refers to the sum of the daily average temperatures over a certain period, reflecting the thermal conditions of the area. Based on the plant species' growth cycle data and the accumulated temperature data of the target loss-risk area, the optimal planting time window for each plant species in different slope areas is calculated. Thermal conditions may vary in different slope areas; the accumulated temperature in sunny slopes is usually higher than in shady slopes. For each plant species, the minimum accumulated temperature requirement for its growth is determined. Then, based on the changes in accumulated temperature in different slope areas of the target loss-risk area, the time period that meets the minimum accumulated temperature requirement for plant growth is identified. This time period is the optimal planting time window for that plant species in that slope area. For example, for a plant species that requires 1000 degrees of accumulated temperature to germinate and grow, in a low-slope sunny area, according to the accumulated temperature data, the accumulated temperature can reach 1000 degrees in April and May each year. Therefore, April and May are the optimal planting time window for this plant species in a low-slope sunny area.

[0210] Step S14256: Generate a planting plan table divided by slope based on the optimal planting time window and maintenance cycle data, and associate the spatial coordinates of each zone with the plant species configuration scheme.

[0211] Maintenance cycle data refers to the length of time and maintenance measures required for each plant species during its growth process. Based on the calculated optimal planting time window and maintenance cycle data for each plant species in different slope areas, a planting plan table by slope is generated. The planting plan table includes information such as the plant species, quantity, and maintenance measures that should be planted in each slope area at different times. For example, in a low-slope, sunny area, a certain herbaceous plant is planted in April and May at a density of 10 plants per square meter, with a maintenance cycle of 3 months, including regular watering and fertilization. Simultaneously, the spatial coordinates of each slope area are associated with the corresponding plant species configuration scheme. In GIS software, the planting plan table can be linked to the slope zoning layer of the target loss-risk area, with each slope zone corresponding to a specific plant species configuration scheme and planting time schedule. This forms a complete, slope-zoning vegetation restoration planting plan, providing detailed guidance for actual vegetation restoration work.

[0212] Step S150: Feed back the governance priority ranking and vegetation restoration strategy to the monitoring platform to trigger the regional governance task allocation operation.

[0213] In this embodiment, after generating the priority ranking of the target loss-risk areas and the vegetation restoration strategy, the above information is fed back to the monitoring platform. The monitoring platform is a system that integrates data management, analysis, and task allocation functions.

[0214] The priorities for remediation and vegetation restoration strategies are uploaded to the monitoring platform as data files. Upon receiving this information, the monitoring platform processes and analyzes it. For example, firstly, based on the priorities, the remediation order for each target area at risk of loss is determined. Then, combined with the vegetation restoration strategies, a detailed remediation task plan is developed.

[0215] For example, for the highest-priority target loss risk areas, the monitoring platform will allocate corresponding human and material resources for vegetation planting and maintenance work based on the planting density and maintenance cycle in the vegetation restoration strategy. At the same time, the monitoring platform will assign the governance tasks to the relevant governance teams or departments, clearly defining the tasks and responsibilities of each team or department.

[0216] After receiving a governance task, the governance team or department can carry out specific governance work according to the vegetation restoration strategy and task requirements. During the governance process, the monitoring platform will monitor the progress and effects in real time, and make adjustments and optimizations based on the actual situation to ensure that the soil and water conservation work can achieve the expected results.

[0217] Figure 2The illustration shows exemplary hardware and software components of a remote sensing image-based dynamic soil erosion monitoring system 100 that can implement the ideas of this application, according to some embodiments of this application. For example, a processor 120 can be used in the remote sensing image-based dynamic soil erosion monitoring system 100 and to perform the functions in this application.

[0218] The soil erosion dynamic monitoring system 100 based on remote sensing imagery can be a general-purpose server or a special-purpose server; both can be used to implement the soil erosion dynamic monitoring method based on remote sensing imagery of this application. Although only one server is shown in this application, for convenience, the functions described in this application can be implemented in a distributed manner on multiple similar platforms to balance the processing load.

[0219] For example, a dynamic soil erosion monitoring system 100 based on remote sensing imagery may include a network port 110 connected to a network, one or more processors 120 for executing program instructions, a communication bus 130, and various forms of storage media 140, such as a disk, ROM, or RAM, or any combination thereof. Exemplarily, the dynamic soil erosion monitoring system 100 based on remote sensing imagery may also include program instructions stored in ROM, RAM, or other types of non-transitory storage media, or any combination thereof. The methods of this application can be implemented according to these program instructions. The dynamic soil erosion monitoring system 100 based on remote sensing imagery also includes an input / output (I / O) interface 150 between the computer and other input / output devices.

[0220] For ease of explanation, only one processor is described in the image-based dynamic soil erosion monitoring system 100. However, it should be noted that the image-based dynamic soil erosion monitoring system 100 of this application may also include multiple processors. Therefore, the steps performed by one processor as described in this application may also be performed jointly by multiple processors or individually. For example, if the processor of the image-based dynamic soil erosion monitoring system 100 performs steps A and B, it should be understood that steps A and B may also be performed jointly by two different processors or individually by one processor. For example, the first processor performs step A, the second processor performs step B, or the first processor and the second processor jointly perform steps A and B.

[0221] Furthermore, this embodiment of the invention also provides a readable storage medium, which has computer-executable instructions pre-set in it. When the processor executes the computer-executable instructions, the above-mentioned method for dynamic monitoring of soil erosion based on remote sensing images is realized.

[0222] It should be noted that, in order to simplify the description of the present invention and thus help to understand one or more embodiments of the invention, multiple features may sometimes be grouped into one embodiment, drawing or description thereof in the foregoing description of the embodiments of the present invention.

Claims

1. A method for dynamic monitoring of soil erosion based on remote sensing imagery, characterized in that, The method includes: Acquire a multi-temporal remote sensing image dataset of the target area, wherein the multi-temporal remote sensing image dataset includes remote sensing image subsets of multiple time periods, and each remote sensing image subset contains land cover information of at least one spectral band. The surface feature extraction process is performed on the multi-temporal remote sensing image dataset to obtain the vegetation cover features, topographic slope features, and soil exposure features of each remote sensing image subset. A spatiotemporal correlation-based soil erosion prediction model is constructed. The vegetation coverage characteristics, terrain slope characteristics and soil exposure characteristics are input into the soil erosion prediction model to predict soil erosion and generate a soil erosion level distribution map of the target area. Based on the soil erosion level distribution map, target erosion risk areas are identified, and a priority ranking of treatment and vegetation restoration strategies are generated for the target erosion risk areas. The governance priority ranking and vegetation restoration strategy are fed back to the monitoring platform to trigger the regional governance task allocation operation; The step of identifying target erosion risk areas based on the soil erosion level distribution map includes: Spatial cluster analysis was performed on the soil erosion level distribution map to extract continuously distributed areas of severe erosion level as candidate erosion areas; Obtain land use type data of the candidate lost areas, and filter out areas belonging to cultivated land, forest land or grassland as areas to be verified. The pre-trained erosion gully detection model is invoked to calculate the gully density of the remote sensing image of the area to be verified, and the area to be verified where the gully density exceeds the first density threshold is marked as the confirmed loss area. The confirmed loss areas and the candidate loss areas are spatially overlaid and analyzed to remove abnormal areas with abrupt terrain changes or interference from man-made buildings, thus obtaining the final set of target loss risk areas. The step of calling the pre-trained erosion gully detection model to calculate the gully density on the remote sensing image of the area to be verified includes: Acquire multispectral remote sensing images of the area to be verified, and extract reflectance data in the visible and near-infrared bands; Edge enhancement filtering is applied to the reflectivity data in the visible light band to generate an enhanced surface texture feature map; Based on the reflectivity data of the near-infrared band and the enhanced surface texture feature map, a multi-scale morphological opening operation is performed to eliminate isolated noise points and retain continuous gully structures, resulting in a processed surface texture feature map. The processed surface texture feature map is input into a pre-trained gully semantic segmentation model, which outputs a binary mask image of the gully pixels. Connected component labeling is performed on the binary mask image, and the groove length and width data of each connected component are statistically analyzed; The gully extension density per unit area is calculated based on the gully length and width data, and mapped to a gully density level using a preset density grading table. The gully density level is overlaid with the geographical boundary of the candidate loss area to generate a gully density distribution layer with spatial attributes; The generation of governance priority ranking and vegetation restoration strategies for the target loss-risk areas includes: Acquire historical rainfall data, soil erosion rate data, and vegetation restoration cost data for the target erosion risk area; A natural recovery potential score is calculated based on the soil erosion rate data and rainfall data. The natural recovery potential score is negatively correlated with the erosion rate and positively correlated with the rainfall. Based on the vegetation restoration cost data and natural restoration potential score, a multi-objective optimization function is constructed, and the comprehensive priority index of the loss risk area for each objective is obtained by solving the function. The target areas at risk of loss are sorted from high to low according to the comprehensive priority index to generate a governance priority ranking. Based on the topographic slope characteristics and soil type of the target loss-risk area, a pre-set vegetation restoration scheme library is matched to generate a vegetation restoration strategy that includes planting density and maintenance cycle, specifically including: Extract digital elevation model data of the target loss risk area and calculate the slope aspect distribution map based on the elevation gradient; The terrain slope features are spatially overlaid with the slope aspect distribution map to divide the sunny slope area into a shady slope area. Obtain soil water-holding capacity parameters from the regional soil type database, and generate a soil-slope water-holding capacity matrix by combining the slope region type; Based on the soil-slope water-holding capacity matrix, query the list of suitable plant species in the vegetation restoration scheme database and filter out plant species; Based on the growth cycle data of the selected plant species and the accumulated temperature data of the target loss risk area, the optimal planting time window for each plant species in different slope areas is calculated; a planting plan table by slope is generated according to the optimal planting time window and maintenance cycle data, and the spatial coordinates of each zone are associated with the plant species configuration scheme.

2. The method for dynamic monitoring of soil erosion based on remote sensing imagery according to claim 1, characterized in that, The surface feature extraction process of the multi-temporal remote sensing image dataset yields vegetation cover features, topographic slope features, and soil bareness features for each remote sensing image subset, including: Radiometric correction is performed on the subset of remote sensing images to obtain standardized remote sensing images; The pre-trained vegetation index calculation model is invoked to calculate the vegetation coverage feature based on the reflectance difference between the near-infrared band and the red band in the standardized remote sensing image. The vegetation coverage feature is positively correlated with the reflectance difference. The standardized remote sensing image is subjected to digital elevation model inversion processing, and the terrain slope feature is calculated based on the elevation change rate between adjacent pixels, wherein the elevation change rate threshold of the terrain slope feature is dynamically adjusted according to a preset terrain classification rule. The reflectance ratio of the shortwave infrared band to the visible light band in the standardized remote sensing image is extracted, and combined with the regional soil reflectance baseline in the soil type database, the soil exposure characteristics are determined, wherein the degree of soil exposure is negatively correlated with the reflectance ratio.

3. The method for dynamic monitoring of soil erosion based on remote sensing imagery according to claim 1, characterized in that, The construction of the soil erosion prediction model based on spatiotemporal correlation includes: A historical soil erosion labeling dataset is obtained, which includes multiple historical remote sensing image samples and their corresponding erosion level labels. The erosion level labels are determined based on field survey data and soil and water conservation indicators. For each historical remote sensing image sample, the aforementioned surface feature extraction process is performed to obtain historical vegetation cover features, historical topographic slope features, and historical soil exposure features. The historical vegetation cover characteristics, historical topographic slope characteristics, and historical soil exposure characteristics are arranged in chronological order to generate a spatiotemporal feature sequence. A bidirectional temporal convolutional network model is constructed, and the spatiotemporal feature sequence is input into the bidirectional temporal convolutional network model for forward and backward feature propagation to capture the dependencies of features changing over time. The bidirectional temporal convolutional network model is trained under supervision based on the erosion level label until the error between the predicted erosion level output by the bidirectional temporal convolutional network model and the erosion level label is lower than a preset threshold, thus obtaining the trained soil and water loss prediction model.

4. The method for dynamic monitoring of soil erosion based on remote sensing imagery according to claim 1, characterized in that, The process of inputting the vegetation cover characteristics, terrain slope characteristics, and soil exposure characteristics into the soil and water loss prediction model to predict soil and water loss and generate a distribution map of soil and water loss levels in the target area includes: The vegetation cover characteristics, topographic slope characteristics, and soil bareness characteristics of the current time period are compared with the corresponding characteristics of historical time periods to obtain the changes in vegetation cover, topographic slope, and soil bareness. The changes in vegetation cover, topographic slope, and soil exposure are input into the bidirectional temporal convolutional network model of the soil erosion prediction model to generate a erosion probability score for each pixel, wherein the erosion probability score reflects the increase in erosion risk in the current time period relative to historical periods. According to the preset churn level classification rules, the churn probability score is mapped to discrete churn levels; The erosion levels are associated with the corresponding geographic coordinates to generate a spatially resolved distribution map of soil erosion levels, wherein the spatial resolution is consistent with the resolution of the subset of remote sensing images.

5. The method for dynamic monitoring of soil erosion based on remote sensing imagery according to claim 3, characterized in that, The construction of the bidirectional temporal convolutional network model includes: The forward convolutional layer performs a sliding window convolution operation on the spatiotemporal feature sequence in forward temporal order, and outputs a forward feature map containing forward temporal dependencies. The backward convolutional layer performs a sliding window convolution operation on the spatiotemporal feature sequence in reverse temporal order, and outputs a backward feature map containing inverse temporal dependencies. A depth convolution module is set at the output of the forward convolution layer and the backward convolution layer to perform single-channel spatial convolution operation on the forward feature map and the backward feature map respectively to generate a depth convolution feature map; A pointwise convolution module is set at the output of the depth convolution module to perform a channel fusion operation with a 1×1 convolution kernel on each depth convolution feature map to generate a dimension-reduced fused feature map. A channel attention module is set at the output of the pointwise convolution module to concatenate the forward fusion feature map and the inverse fusion feature map along the channel dimension to generate a concatenated feature map. The weight coefficients of each channel in the stitched feature map are calculated by the channel attention module, and a weighted summation operation is performed on the stitched feature map according to the weight coefficients to generate an attention-weighted interactive feature map. A multi-scale pooling module is set at the output of the channel attention module to perform max pooling operations of different scales on the interactive feature map to generate a multi-scale pooling feature map. Multi-scale pooling feature maps are concatenated along the channel dimension to generate feature vectors containing multi-scale spatiotemporal features.

6. The method for dynamic monitoring of soil erosion based on remote sensing imagery according to claim 5, characterized in that, The process of supervising the training of the bidirectional temporal convolutional network model based on the erosion level label until the error between the predicted erosion level output by the bidirectional temporal convolutional network model and the erosion level label is lower than a preset threshold, thereby obtaining the trained soil erosion prediction model, includes: The spatiotemporal feature sequence is divided into continuous time segments in chronological order. Some time segments are randomly selected as the training set, and the remaining time segments are used as the validation set. Each time segment contains a preset number of historical remote sensing image samples. The training set is input into the forward and backward convolutional layers of the bidirectional temporal convolutional network model to generate the forward and backward feature maps of the training set. The forward and inverse feature maps are processed by the deep convolution module and the pointwise convolution module to generate the fused feature map of the training set. The fused feature map is then processed by the channel attention module and the multi-scale pooling module to generate the multi-scale feature vector of the training set. The multi-scale feature vector is then input into the fully connected layer to output the probability distribution of the predicted churn level of the training set. Based on the true category index of the churn level label, extract the probability value of the corresponding category in the predicted churn level probability distribution, calculate the difference between the probability value and 1, and generate the prediction error coefficient for each training sample. The prediction error coefficient is exponentially operated on and multiplied by a preset focusing coefficient to generate a sample-level weight adjustment factor. The weight adjustment factor is then multiplied by the cross-entropy loss to generate the focus loss value for each training sample. The category weight coefficient is set according to the severity of the churn level label. The focus loss value is multiplied by the category weight coefficient to generate a weighted focus loss value. The weighted focus loss value is then averaged in batches to generate the total loss value for the current training batch. During backpropagation, the gradient of the total loss value with respect to the model parameters is calculated, and the norm of the gradient vector is detected. When the norm of the gradient vector exceeds a preset clipping threshold, the gradient vector is scaled proportionally until the norm equals the clipping threshold, and the scaled gradient vector is used to update the model parameters to complete one training iteration. Repeat the forward propagation process on the validation set and calculate the prediction accuracy of the validation set. If the prediction accuracy of the validation set does not improve after N consecutive training iterations, stop training and save the current model parameters.

7. A dynamic monitoring system for soil erosion based on remote sensing imagery, characterized in that, The dynamic monitoring system for soil erosion based on remote sensing imagery includes a processor and a memory, the memory being connected to the processor. The memory is used to store programs, instructions, or code, and the processor is used to execute the programs, instructions, or code in the memory to implement the dynamic monitoring method for soil erosion based on remote sensing imagery as described in any one of claims 1-6.

Images