Soil and water conservation intelligent monitoring method and system based on multi-source data fusion

Abstract

This invention provides a method and system for intelligent monitoring of soil and water conservation based on multi-source data fusion, belonging to the field of ecological and environmental monitoring technology. First, it collects rainfall remote sensing inversion data and soil property spatial distribution data of the target watershed. Then, it performs spatiotemporal evolution analysis of rainfall erosivity on the rainfall remote sensing inversion data to generate a spatiotemporal dynamic field of rainfall erosivity. It then performs spatial differentiation analysis of soil erosability on the soil property spatial distribution data to generate a spatial heterogeneous field of soil erosability. Next, it inputs the spatiotemporal dynamic field of rainfall erosivity and the spatial heterogeneous field of soil erosability into a soil erosion coupling response model for spatiotemporal overlay and stress coupling processing to generate a spatiotemporal coupled field of soil erosion risk. Finally, it extracts key monitoring areas for soil and water conservation based on this spatiotemporal coupled field of soil erosion risk and generates intelligent monitoring results for soil and water conservation. This invention can comprehensively and accurately monitor soil erosion, helping to improve the efficiency and effectiveness of soil and water conservation work.

Description

Technical Field

[0001] This invention relates to the field of ecological environment monitoring technology, and more specifically, to an intelligent monitoring method and system for soil and water conservation based on multi-source data fusion. Background Technology

[0002] Soil and water conservation plays an irreplaceable role in maintaining ecological balance, ensuring water resource security, and promoting sustainable agricultural development. Traditional soil and water conservation monitoring methods mainly rely on field surveys and manual measurements. While these methods can obtain data with a certain degree of accuracy, they also have many limitations.

[0003] On the one hand, field surveys and manual measurements require significant manpower, resources, and time, and the monitoring scope is limited, making it difficult to conduct comprehensive and timely monitoring of large watersheds. This is especially true in areas with complex terrain and poor transportation, where monitoring work faces enormous challenges. On the other hand, data obtained through traditional methods are often discrete and localized, making it difficult to accurately reflect the spatiotemporal dynamics of soil erosion. For example, rainfall and soil properties vary significantly across different times and spaces, and traditional methods struggle to comprehensively consider the impact of these factors on soil erosion, resulting in inaccurate and incomplete monitoring results.

[0004] With the development of remote sensing technology and geographic information systems, although the efficiency and scope of soil and water conservation monitoring have been improved to some extent, most of them are still based on a single data source. They have failed to fully integrate the advantages of multi-source data, and cannot comprehensively and accurately analyze the causes and risks of soil erosion, making it difficult to meet the high requirements of modern soil and water conservation monitoring. Summary of the Invention

[0005] 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 intelligent monitoring of soil and water conservation based on multi-source data fusion, the method comprising: The remote sensing inversion data of rainfall and the spatial distribution data of soil properties of the target watershed are collected. The remote sensing inversion data of rainfall includes the time series of rainfall intensity and the spatial distribution map of rainfall accumulation in multiple grid cells. The spatial distribution data of soil properties includes soil texture type parameters and soil organic matter content values ​​of multiple soil sampling points. The rainfall intensity time series and rainfall accumulation spatial distribution map in the rainfall remote sensing inversion data are subjected to spatiotemporal evolution analysis of rainfall erosivity to generate a spatiotemporal dynamic field of rainfall erosivity in the target watershed. The spatiotemporal dynamic field of rainfall erosivity includes the rainfall erosivity status identifier and rainfall erosivity peak migration trajectory for each spatiotemporal unit. The soil texture type parameter and soil organic matter content value in the spatial distribution data of soil properties are subjected to spatial differentiation analysis of soil erodibility to generate a spatial heterogeneous field of soil erodibility for the target watershed. The spatial heterogeneous field of soil erodibility includes the soil erodibility status identifier and soil structure stability index of each spatial unit. The spatiotemporal dynamic field of rainfall erosion force and the spatial heterogeneous field of soil erodibility are input into a pre-constructed water and soil loss coupled response model for spatiotemporal superposition and stress coupling processing to generate a spatiotemporal coupled field of water and soil loss risk for the target watershed. The spatiotemporal coupled field of water and soil loss risk includes water and soil loss risk level parameters and water and soil loss sensitive area boundary sequence for each spatiotemporal unit. Based on the soil and water loss risk level parameters and the boundary sequence of soil and water loss sensitive areas in the spatiotemporal coupling field of the soil and water loss risk, the key monitoring areas for soil and water conservation in the target watershed are extracted, and intelligent monitoring results for soil and water conservation, including the boundary polygons of the monitoring areas and the dynamic monitoring frequency of the monitoring areas, are generated.

[0006] Furthermore, embodiments of the present invention also provide a water and soil conservation intelligent monitoring system based on multi-source data fusion, comprising: A processor; a machine-readable storage medium for storing machine-executable instructions of the processor; wherein the processor is configured to execute the above-described intelligent water and soil conservation monitoring method based on multi-source data fusion by executing the machine-executable instructions.

[0007] In another aspect, embodiments of the present invention also provide a computer program product, the computer program product including machine-executable instructions, the machine-executable instructions being stored in a computer-readable storage medium, the processor of the water and soil conservation intelligent monitoring system based on multi-source data fusion reading the machine-executable instructions from the computer-readable storage medium, the processor executing the machine-executable instructions, causing the water and soil conservation intelligent monitoring system based on multi-source data fusion to execute the above-mentioned water and soil conservation intelligent monitoring method based on multi-source data fusion.

[0008] Based on the above, firstly, rainfall remote sensing inversion data and soil property spatial distribution data of the target watershed are collected. Spatiotemporal evolution analysis of rainfall erosivity is performed on the rainfall remote sensing inversion data to generate a spatiotemporal dynamic field of rainfall erosivity, which can present the changes in rainfall erosivity and peak migration trajectory at different times and spaces. Spatial differentiation analysis of soil erosability is performed on the spatial distribution data of soil properties to generate a spatial heterogeneous field of soil erosability, which can accurately reflect the spatial differences in soil erosability and the stability of soil structure. The spatiotemporal dynamic field of rainfall erosivity and the spatial heterogeneous field of soil erosability are input into a pre-constructed water and soil erosion coupled response model for spatiotemporal overlay and stress coupling processing to generate a spatiotemporal coupled field of water and soil erosion risk. This comprehensively considers the interaction between rainfall and soil factors, enabling a more accurate assessment of water and soil erosion risk levels and determination of sensitive area boundaries. Finally, based on the spatiotemporal coupled field of water and soil erosion risk, key water and soil conservation monitoring areas are extracted and intelligent monitoring results are generated, achieving precise monitoring and effective supervision of water and soil erosion, and helping to improve the efficiency and effectiveness of water and soil conservation work. Attached Figure Description

[0009] Figure 1 This is a schematic diagram of the execution flow of the intelligent monitoring method for soil and water conservation based on multi-source data fusion provided in an embodiment of the present invention.

[0010] Figure 2 This is a schematic diagram of exemplary hardware and software components of the intelligent water and soil conservation monitoring system based on multi-source data fusion provided in an embodiment of the present invention. Detailed Implementation

[0011] Figure 1 This is a flowchart illustrating an intelligent monitoring method for soil and water conservation based on multi-source data fusion, provided in one embodiment of the present invention. A detailed description follows.

[0012] Step S110: Collect rainfall remote sensing inversion data and soil property spatial distribution data of the target watershed. The rainfall remote sensing inversion data includes rainfall intensity time series and rainfall accumulation spatial distribution map of multiple grid cells. The soil property spatial distribution data includes soil texture type parameters and soil organic matter content values ​​of multiple soil sampling points.

[0013] In this embodiment, a typical hilly agricultural target watershed located in a humid monsoon climate zone is used as the application scenario. The total area of ​​the target watershed is S, where S = 3200 square kilometers. Remote sensing raw data of the target watershed is acquired using a microwave radiometer and a visible-infrared scanning radiometer from a meteorological satellite. After atmospheric correction, geometric correction, and cloud masking, a rasterized rainfall product covering the entire target watershed is generated. The target watershed is divided into regular grid cells with a spatial resolution of D × D, where D = 1000 meters, resulting in N grid cells, N = 3200. For each grid cell, hourly rainfall intensity values ​​for the past Y consecutive complete hydrological years are extracted from the remote sensing rainfall product, where Y = 10, forming the rainfall intensity time series R(t) for that grid cell. The length of R(t) is L, where L = Y × 365 × 24. Meanwhile, the annual cumulative rainfall values ​​for each grid cell within the same Y hydrological years are extracted from the same remote sensing rainfall product. The annual cumulative rainfall values ​​are arranged according to spatial coordinates to generate a spatial distribution map P_y of the cumulative rainfall for each year. The size of P_y is R_num in the spatial row direction × C_num in the spatial column direction, where R_num=40 and C_num=80. Each grid cell stores a floating-point cumulative rainfall value P_ij_y.

[0014] On the other hand, K soil sampling points were set up within the target watershed according to the principle of uniform grid distribution, K=400, and the spatial coordinates of each soil sampling point were Gaussian projection plane coordinates (X_k, Y_k). At each soil sampling point, a soil sample with a surface depth of H (H=20 cm) was collected using the ring cutter method. After the soil samples were brought back to the laboratory, the soil particle composition was determined using a laser particle size analyzer. According to the soil texture classification standard, the soil texture was divided into four types: sandy loam, loam, clay loam, and clay, generating the soil texture type parameter T_k for each soil sampling point. At the same time, the organic carbon content in the soil sample was determined using the potassium dichromate oxidation external heating method, and the organic matter content value O_k was obtained by multiplying it by the conversion factor 1.724. All the above data collection processes were anonymized for sensitive information involving geographical location, the specific sampling point coordinates were replaced with relative coordinate offsets, and the data storage and transmission process was encrypted.

[0015] Step S120: Perform spatiotemporal evolution analysis on the rainfall intensity time series and rainfall accumulation spatial distribution map in the rainfall remote sensing inversion data to generate a spatiotemporal dynamic field of rainfall erosivity for the target watershed. The spatiotemporal dynamic field of rainfall erosivity includes the rainfall erosivity status identifier and the peak migration trajectory of rainfall erosivity for each spatiotemporal unit.

[0016] Step S121: Extract the spatial coordinates of each grid cell, the rainfall intensity time series corresponding to the grid cell, and the rainfall accumulation value in the spatial distribution map of the rainfall accumulation corresponding to the grid cell from the rainfall remote sensing inversion data. Then, associate and store the spatial coordinates of each grid cell with the rainfall intensity time series and rainfall accumulation value of the grid cell to generate a gridded set of rainfall parameters containing the spatial coordinates of the grid cells and the corresponding rainfall parameter series.

[0017] For a grid cell with spatial row index i and spatial column index j, where i = 1 to R_num and j = 1 to C_num, extract the spatial coordinates (X_ij, Y_ij) of the center point of this grid cell. Extract the rainfall intensity time series R_ij(t) corresponding to this grid cell, where t = 1 to L. Extract the rainfall accumulation value P_ij(y) at the location of this grid cell in the spatial distribution map of rainfall accumulation over Y years, where y = 1 to Y. Associate the spatial coordinates (X_ij, Y_ij), the rainfall intensity time series R_ij(t), and the rainfall accumulation value P_ij(y) according to the spatial index of this grid cell to generate a data record containing the spatial coordinates and the corresponding rainfall parameter series. All data records of N grid cells are organized into a two-dimensional array structure according to the spatial row index and spatial column index. The number of rows of the two-dimensional array is R_num and the number of columns is C_num. Each array element stores the spatial coordinates of a grid cell, a rainfall intensity time series of length L, and a rainfall cumulative value series of length Y. The above two-dimensional array structure is the gridded set of rainfall parameters.

[0018] Step S122: Perform rainfall event segmentation processing on the rainfall intensity time series of each grid unit in the rainfall parameter gridded set. Divide the time period in the rainfall intensity time series where the continuous rainfall intensity value exceeds the preset rainfall intensity threshold into an independent rainfall event. Record the start time point, end time point, and peak rainfall intensity value within the time period of each independent rainfall event to generate a rainfall event feature set for each grid unit.

[0019] For each grid cell, extract its rainfall intensity time series \(R_{ij}(t)\). Set a preset rainfall intensity threshold \(R_{th}\), where \(R_{th} = 0.5\). Starting from \(t = 1\), traverse \(R_{ij}(t)\). When it is detected that \(R_{ij}(t)\geq R_{th}\) and \(R_{ij}(t - 1)<R_{th}\), mark the current time point \(t\) as the starting time point \(T_s\) of a rainfall event. Continue traversing backward. When it is detected that \(R_{ij}(t)<R_{th}\) and \(R_{ij}(t - 1)\geq R_{th}\), mark the current time point \(t - 1\) as the ending time point \(T_e\) of this rainfall event. Record the rainfall intensity values at all time points from \(T_s\) to \(T_e\), and find the maximum value among them as the peak rainfall intensity value \(R_p\) of this independent rainfall event. Repeat the above process until the entire \(R_{ij}(t)\) is traversed to obtain all independent rainfall events of this grid cell. For each independent rainfall event, record its starting time point \(T_s(q)\), ending time point \(T_e(q)\), and peak rainfall intensity value \(R_p(q)\), where \(q\) is the serial number of the independent rainfall event within this grid cell, \(q = 1\) to \(Q_{ij}\), and \(Q_{ij}\) is the total number of independent rainfall events within this grid cell. Organize the records of all the above independent rainfall events into a list structure in the chronological order of event occurrence. Each list element contains three fields: starting time point, ending time point, and peak rainfall intensity value. The above list structure is the rainfall event feature set of this grid cell.

[0020] Step S123: Extract the rainfall intensity time series of each independent rainfall event from the rainfall event feature set of each grid cell, calculate the rainfall erosivity index of each independent rainfall event based on this rainfall intensity time series, and sum up the rainfall erosivity indices of all independent rainfall events within the same grid cell to generate the time period cumulative rainfall erosivity index of this grid cell.

[0021] For each independent rainfall event \(q\), extract the rainfall intensity time sub - series \(R_{ij}(t)\) from \(T_s(q)\) to \(T_e(q)\). Calculate the rainfall erosivity index \(EI(q)\) of this independent rainfall event. \(EI(q)=[\sum_{t = T_s(q)}^{T_e(q)}(R_{ij}(t))^2]\times F_t\), where \(F_t = 1\). Repeat the above calculation for all independent rainfall events within each grid cell to obtain \(EI(q)\) for each independent rainfall event within this grid cell. Then sum up all \(EI(q)\) within this grid cell, \(E_{ij\_total}=\sum_{q = 1}^{Q_{ij}}EI(q)\). Perform the above calculation for each grid cell in the gridded set of rainfall parameters respectively to generate \(E_{ij\_total}\) for each grid cell.

[0022] Step S124: Extract the spatial distribution features of the time-period cumulative rainfall erosivity index of all grid units in the rainfall parameter gridded set, calculate the absolute value of the difference between the time-period cumulative rainfall erosivity index of adjacent grid units, mark the boundary between adjacent grid units where the absolute value of the difference exceeds a preset difference threshold as the spatial abrupt boundary of rainfall erosivity, and record the spatial location coordinate sequence of all spatial abrupt boundary of rainfall erosivity to generate a spatial differentiation boundary map of rainfall erosivity.

[0023] Set a preset difference threshold D_th, D_th=500. For each grid cell (i, j), calculate ΔE_east=|E_ij_total-E_i(j+1)_total|, ΔE_south=|E_ij_total-E_(i+1)j_total|, ΔE_southeast=|E_ij_total-E_(i+1)(j+1)_total|. When ΔE_east≥D_th, mark the common boundary between grid cells (i, j) and (i, j+1) as the spatial abrupt boundary of rainfall erosion force, and record the spatial coordinate sequence of the two endpoints of the common boundary, with the endpoint coordinates being (X_ij, Y_ij) and (X_i(j+1), Y_i(j+1)). When ΔE_south≥D_th, the common boundary between grid cell (i,j) and grid cell (i+1,j) is marked as a spatial abrupt change boundary of rainfall erosivity. The spatial coordinate sequences of the two endpoints of this common boundary are recorded as (X_ij, Y_ij) and (X_(i+1)j, Y_(i+1)j). When ΔE_southeast≥D_th, the diagonal boundary between grid cell (i,j) and grid cell (i+1,j+1) is marked as a spatial abrupt change boundary of rainfall erosivity. The spatial coordinate sequences of the two endpoints of this diagonal boundary are recorded as (X_ij, Y_ij) and (X_(i+1)(j+1), Y_(i+1)(j+1)). The boundary endpoint coordinate sequences of all marked spatial abrupt change boundaries of rainfall erosivity are organized according to their spatial location to generate a set of line features. Each line feature in this set of line features corresponds to a spatial abrupt change boundary of rainfall erosivity. This set of line features is the spatial differentiation boundary map of rainfall erosivity.

[0024] Step S125: Taking each grid cell as the center, extract the time - period cumulative rainfall erosivity indices of all adjacent grid cells within the preset neighborhood range of this grid cell, calculate the statistical average value of these time - period cumulative rainfall erosivity indices as the neighborhood average rainfall erosivity index of this grid cell, compare the time - period cumulative rainfall erosivity index of each grid cell with the preset set of erosivity classification thresholds, and assign a classification code to this grid cell according to the threshold interval it belongs to, generating the rainfall erosivity status identifier of this grid cell.

[0025] Set a neighborhood window size parameter W, where W = 3. That is, taking the current grid cell as the center, extract all grid cells within its surrounding 3×3 neighborhood window. For grid cells located on the spatial boundary, only extract the actually existing adjacent grid cells. Extract the time - period cumulative rainfall erosivity indices E_mn_total of all grid cells within the neighborhood window of the current grid cell (i, j), where m ranges from max(1, i - 1) to min(R_num, i + 1), and n ranges from max(1, j - 1) to min(C_num, j + 1). Calculate the statistical average value E_avg_ij of these E_mn_total, E_avg_ij=(ΣE_mn_total) / N_neighbor. Set the preset set of erosivity classification thresholds, including the first classification threshold G1, the second classification threshold G2, and the third classification threshold G3, where G1 = 1000, G2 = 3000, and G3 = 6000. Compare the E_ij_total of the current grid cell with G1, G2, and G3. When E_ij_total < G1, assign the rainfall erosivity status identifier S_rain_ij of this grid cell as 1. When G1≤E_ij_total < G2, assign S_rain_ij as 2. When G2≤E_ij_total < G3, assign S_rain_ij as 3. When E_ij_total≥G3, assign S_rain_ij as 4. Repeat the above neighborhood average calculation and classification comparison process for all grid cells in the rainfall parameter grid set, generating S_rain_ij for each grid cell.

[0026] Step S126: Repeatedly execute the steps of calculating the time - period cumulative rainfall erosivity index and assigning the rainfall erosivity status identifier for the rainfall parameter grid set obtained at different time sampling points, generating the spatial distribution map of the rainfall erosivity status identifier corresponding to each time sampling point. Overlay and compare the spatial distribution maps of the rainfall erosivity status identifiers of adjacent time sampling points, extract the grid cells where the category conversion of the rainfall erosivity status identifier occurs as the rainfall erosivity status conversion units, and record the spatial coordinate sequences of these rainfall erosivity status conversion units.

[0027] The Y hydrological years are divided into Y time sampling points, with each time sampling point y corresponding to the data of the y-th hydrological year, y=1 to Y. For each time sampling point y, steps S121 to S125 are repeated, but the rainfall intensity time series R_ij(t) and the rainfall accumulation value P_ij(y) in step S121 only use the data of the y-th hydrological year, that is, the value range of t is the 8760 time points of the y-th hydrological year. A spatial distribution map S_rain_ij(y) of rainfall erosivity state identifier corresponding to each time sampling point y is generated. This spatial distribution map is a two-dimensional matrix of size R_num×C_num, and each matrix element stores a classification code value. For two adjacent time sampling points y and y+1, S_rain_ij(y) and S_rain_ij(y+1) are superimposed and compared. All grid cells (i, j) are traversed to determine whether S_rain_ij(y) and S_rain_ij(y+1) are equal. When the two are not equal, the grid cell is marked as a rainfall erosivity state transition cell, and the spatial coordinates (X_ij, Y_ij) of the grid cell are recorded. The spatial coordinates of all marked rainfall erosivity state transition cells are organized into a coordinate sequence set in chronological order, with each coordinate sequence corresponding to the position of the state transition cell within a time sampling interval.

[0028] Step S127: Connect the spatial coordinate sequence of the rainfall erosivity state transformation unit corresponding to each time sampling point in chronological order to generate the rainfall erosivity peak migration trajectory reflecting the spatial transformation trajectory of the rainfall erosivity state identifier. Organize the rainfall erosivity state identifier and the rainfall erosivity peak migration trajectory of each spatiotemporal unit according to spatiotemporal coordinates to generate a rainfall erosivity spatiotemporal dynamic field containing the rainfall erosivity state identifier and the rainfall erosivity peak migration trajectory of each spatiotemporal unit.

[0029] Step S126 yields the spatial coordinate sequence of the rainfall erosivity state transition unit corresponding to each time sampling point interval (from y to y+1). For each rainfall erosivity state transition unit, its spatial coordinates are (X_ij, Y_ij). The rainfall erosivity state of this transition unit at the y-th time sampling point is identified as S_rain_ij(y), and the rainfall erosivity state at the (y+1)-th time sampling point is identified as S_rain_ij(y+1). The spatial coordinates of this transition unit at different time sampling points are connected in chronological order. For the same spatial location (X_ij, Y_ij), the position coordinates from the 1st time sampling point to the 2nd time sampling point are connected to form the first trajectory segment, the position coordinates from the 2nd time sampling point to the 3rd time sampling point are connected to form the second trajectory segment, and so on, until the position coordinates from the (Y-1)th time sampling point to the Y-th time sampling point are connected to form the (Y-1)th trajectory segment. Connecting all the trajectory segments end-to-end forms a complete broken line, which represents the peak migration trajectory T_rain_ij of rainfall erosivity at that spatial location. For spatial locations where the rainfall erosivity state identifier has not undergone a category transformation, its T_rain_ij is an empty sequence. The S_rain_ij(y) and T_rain_ij of each spatiotemporal unit are organized according to spatiotemporal coordinates, which include a time index y and a spatial index (i, j). A four-dimensional data structure F_rain is generated, with dimensions Y×R_num×C_num×2, where the first dimension is the time dimension, the second is the spatial row dimension, the third is the spatial column dimension, and the fourth is the feature dimension. The first position in the feature dimension stores S_rain_ij(y), and the second position stores a reference pointer to T_rain_ij. This four-dimensional data structure F_rain represents the spatiotemporal dynamic field of rainfall erosivity.

[0030] Step S130: Perform spatial differentiation analysis on the soil texture type parameter and soil organic matter content value in the spatial distribution data of soil attributes to generate a spatial heterogeneous field of soil erodibility for the target watershed. The spatial heterogeneous field of soil erodibility includes the soil erodibility status identifier and soil structure stability index of each spatial unit.

[0031] Step S131: Extract the spatial coordinates of each soil sampling point, the soil texture type parameter corresponding to the soil sampling point, and the soil organic matter content value corresponding to the soil sampling point from the spatial distribution data of soil attributes. Then, perform association storage processing on the spatial coordinates of each soil sampling point, the soil texture type parameter, and the soil organic matter content value of the soil sampling point to generate a set of discrete soil attribute points containing the spatial coordinates of the soil sampling points and the corresponding soil attribute parameters.

[0032] For the k-th soil sampling point (k=1 to K, K=400), extract its spatial coordinates (X_k, Y_k). Extract the soil texture type parameter T_k corresponding to this soil sampling point, where T_k can be one of four types: sandy loam, loam, clay loam, or clay. Extract the soil organic matter content value O_k corresponding to this soil sampling point. Associate the spatial coordinates (X_k, Y_k), T_k, and O_k to generate a data record containing the spatial coordinates and corresponding soil attribute parameters. Organize the data records of all K soil sampling points into a list structure according to the sampling point number. The list has a length of K, and each list element stores the spatial coordinates, soil texture type parameter, and soil organic matter content value of a sampling point. This list structure is the set of discrete points for soil attributes.

[0033] Step S132: Perform type encoding conversion processing on the soil texture type parameters in the set of discrete soil attribute points, encode sandy loam type as first texture type code, loam type as second texture type code, clay loam type as third texture type code, and clay type as fourth texture type code, and generate soil texture type code value for each soil sampling point.

[0034] A mapping relationship between soil texture type and its encoded value is established, mapping sandy loam to C_sand=1, loam to C_loam=2, clay loam to C_clay_loam=3, and clay to C_clay=4. Based on this mapping relationship, the T_k of each soil sampling point is converted into the corresponding soil texture type encoded value C_type_k, where C_type_k ranges from {1, 2, 3, 4}. This encoding conversion process is then performed on all K soil sampling points in the discrete set of soil attribute points to generate C_type_k for each soil sampling point.

[0035] Step S133: Extract the soil organic matter content value of each soil sampling point from the set of discrete soil attribute points, compare the soil organic matter content value with a preset organic matter content grading threshold sequence, and use the grading code corresponding to the grading interval into which the soil organic matter content value falls as the soil organic matter grading code value of the soil sampling point. The preset organic matter content grading threshold sequence includes a first grading threshold, a second grading threshold, and a third grading threshold. Sampling points with soil organic matter content values ​​less than the first grading threshold are encoded as the first organic matter grading code, sampling points with soil organic matter content values ​​between the first and second grading thresholds are encoded as the second organic matter grading code, sampling points with soil organic matter content values ​​between the second and third grading thresholds are encoded as the third organic matter grading code, and sampling points with soil organic matter content values ​​greater than the third grading threshold are encoded as the fourth organic matter grading code.

[0036] Set a sequence of preset organic matter content classification thresholds, including H1 = 10, H2 = 20, and H3 = 30. Compare O_k of each soil sampling point with H1, H2, and H3. When O_k < H1, C_om_k = 1. When H1 ≤ O_k < H2, C_om_k = 2. When H2 ≤ O_k < H3, C_om_k = 3. When O_k ≥ H3, C_om_k = 4. Perform the above classification comparison process on all K soil sampling points in the soil property discrete point set to generate C_om_k for each soil sampling point.

[0037] In step S134, perform combined coding on the soil texture type coding value and the soil organic matter classification coding value of each soil sampling point. Concatenate the soil texture type coding value and the soil organic matter classification coding value as a string to generate the comprehensive soil erodibility coding value for each soil sampling point, and establish a mapping relationship table between the comprehensive soil erodibility coding value and the preset soil erodibility status identifier. Query the soil erodibility status identifier corresponding to the comprehensive soil erodibility coding value of each soil sampling point according to the mapping relationship table.

[0038] Perform combined coding on C_type_k and C_om_k. The specific method of combined coding is as follows: Convert the value of C_type_k into a two-digit decimal digit string, convert the value of C_om_k into a one-digit decimal digit string, and concatenate the above two-digit digit string and one-digit digit string in sequence to generate a three-digit digit string. This three-digit digit string is the comprehensive soil erodibility coding value E_code_k. The value of E_code_k ranges from 111 to 444, with a total of 64 possible combinations. Establish a mapping relationship table between the comprehensive soil erodibility coding value and the preset soil erodibility status identifier. This mapping relationship table contains 64 entries, each entry corresponding to a value of E_code_k and a corresponding soil erodibility status identifier S_erode_k. S_erode_k is divided into four categories. The first category represents low erodibility, the second category represents medium erodibility, the third category represents high erodibility, and the fourth category represents extremely high erodibility. Query the corresponding S_erode_k in the above mapping relationship table according to the calculated E_code_k for each soil sampling point. Perform the above combined coding and mapping query process on all K soil sampling points in the soil property discrete point set to generate S_erode_k for each soil sampling point.

[0039] Step S135: Spatial interpolation processing is performed on the soil erodibility status markers of all soil sampling points. The inverse distance weighted interpolation method is used to calculate the interpolation result of the soil erodibility status markers at each grid node in the target watershed, generating a spatial distribution map of soil erodibility status markers covering the entire spatial range of the target watershed. The spatial distribution map of soil erodibility status markers includes the spatial coordinates of each grid node and the interpolation result of the soil erodibility status markers corresponding to that grid node.

[0040] The spatial extent of the target watershed is divided into the same grid system as in step S110. Each grid node is the center point of a grid cell, and the spatial coordinates of each grid node are (X_ij, Y_ij). There are N grid nodes in total, N=3200. From step S134, the S_erode_k and spatial coordinates (X_k, Y_k) of K soil sampling points are obtained. For the target grid node (X_ij, Y_ij), the Euclidean distance D_ijk=[(X_ij-X_k)^2+(Y_ij-Y_k)^2]^{1 / 2} between this grid node and each soil sampling point (X_k, Y_k) is first calculated. A search radius R_search=5000 is set, and only soil sampling points with D_ijk≤R_search are considered for interpolation calculation. For soil sampling points that meet the conditions, the weight W_ijk=1 / (D_ijk)^p of each sampling point is calculated, where p=2. The summation of W_ijk for all sampling points meeting the conditions yields W_total_ij = ΣW_ijk. Then, the weighted soil erodibility status indicator value (W_ijk × S_erode_k) for each meeting the conditions is calculated, and all weighted values ​​are summed to obtain W_sum_ij = Σ(W_ijk × S_erode_k). Finally, the interpolation result of the soil erodibility status indicator at the grid node is calculated as S_interp_ij = W_sum_ij / W_total_ij. S_interp_ij is a continuous floating-point value ranging from 1 to 4. This inverse distance weighted interpolation calculation process is repeated for N grid nodes to generate S_interp_ij for each grid node. All S_interp_ij are organized into a two-dimensional matrix structure according to spatial coordinates (X_ij, Y_ij). The number of rows in this two-dimensional matrix is ​​R_num, the number of columns is C_num, and each matrix element stores a floating-point interpolation result. The above two-dimensional matrix structure is the spatial distribution map of soil erodibility status.

[0041] Step S136: For each grid node in the spatial distribution map of soil erodibility status identifiers, extract the interpolation results of soil erodibility status identifiers of all grid nodes within the preset neighborhood window of the grid node, calculate the statistical variance of these interpolation results of soil erodibility status identifiers as the spatial variation parameter of soil erodibility of the grid node, and organize the spatial variation parameters of soil erodibility of all grid nodes according to spatial coordinates to generate a spatial variation distribution map of soil erodibility.

[0042] For each grid node (i, j), extract the S_interp_mn of all grid nodes within a preset neighborhood window surrounding that grid node. The neighborhood window size parameter W_v=5, meaning that all grid nodes within a 5×5 neighborhood window centered on the current grid node are extracted. For grid nodes located on spatial boundaries, only actual adjacent grid nodes are extracted. Let the number of actual grid nodes within the neighborhood window be N_v_ij. Calculate the arithmetic mean μ_v_ij=(ΣS_interp_mn) / N_v_ij of all S_interp_mn within the neighborhood window. Then calculate the square of the difference between each S_interp_mn and μ_v_ij, obtaining (S_interp_mn-μ_v_ij)^2. Summate all the squares of the differences to obtain SS_v_ij=Σ(S_interp_mn-μ_v_ij)^2. Finally, the statistical variance value V_ij = SS_v_ij / N_v_ij is calculated, where V_ij is the spatial variability parameter of soil erodibility for that grid node. The neighborhood window variance calculation process is repeated for N grid nodes to generate V_ij for each node. All V_ij values ​​are organized into a two-dimensional matrix structure according to spatial coordinates (X_ij, Y_ij). This two-dimensional matrix has R_num rows and C_num columns, with each matrix element storing a floating-point variance value. This two-dimensional matrix structure represents the spatial variability distribution map of soil erodibility.

[0043] Step S137: The interpolation result of the soil erodibility status identifier of each grid node is associated with the soil erodibility spatial variability parameter of that grid node and stored to generate the soil erodibility status identifier and soil erodibility spatial variability parameter of each grid node. The soil erodibility status identifier of each grid node is used as the soil erodibility status identifier of that grid node. The soil erodibility spatial variability parameter of each grid node is normalized after taking the reciprocal to generate the soil structure stability index of that grid node. The soil erodibility status identifier and soil structure stability index of each spatial unit are organized according to spatial coordinates to generate a soil erodibility spatial heterogeneous field containing the soil erodibility status identifier and soil structure stability index of each spatial unit.

[0044] From step S135, obtain S_interp_ij for each grid node (i, j), and from step S136, obtain V_ij for each grid node (i, j). Convert S_interp_ij into a discrete soil erodibility state identifier S_erode_ij, with the following conversion method: when S_interp_ij < 1.5, S_erode_ij = 1; when 1.5 ≤ S_interp_ij < 2.5, S_erode_ij = 2; when 2.5 ≤ S_interp_ij < 3.5, S_erode_ij = 3; when S_interp_ij ≥ 3.5, S_erode_ij = 4. For each grid node, calculate its soil structure stability index F_stable_ij. First, calculate the reciprocal of V_ij, V_inv_ij = 1 / V_ij. Then, V_inv_ij is normalized, and all N grid nodes are traversed to find the minimum value V_inv_min and the maximum value V_inv_max of V_inv_ij. For each grid node, F_stable_ij = (V_inv_ij - V_inv_min) / (V_inv_max - V_inv_min). The value of F_stable_ij ranges from 0 to 1. The S_erode_ij and F_stable_ij of each spatial cell are organized according to the spatial coordinates (i, j) to generate a three-dimensional data structure G_erode. The dimensions of G_erode are R_num × C_num × 2, where the first dimension is the spatial row dimension, the second dimension is the spatial column dimension, and the third dimension is the feature dimension. The first position in the feature dimension stores S_erode_ij, and the second position stores F_stable_ij. The above three-dimensional data structure G_erode is the spatial heterogeneous field of soil erodibility.

[0045] Step S140: The spatiotemporal dynamic field of rainfall erosion force and the spatial heterogeneous field of soil erodibility are input into the pre-constructed water and soil loss coupling response model for spatiotemporal superposition and stress coupling processing to generate the spatiotemporal coupling field of water and soil loss risk of the target watershed. The spatiotemporal coupling field of water and soil loss risk includes the water and soil loss risk level parameters of each spatiotemporal unit and the boundary sequence of water and soil loss sensitive areas.

[0046] Step S141: Extract the rainfall erosion force status identifier and rainfall erosion force peak migration trajectory of each spatiotemporal unit from the spatiotemporal dynamic field of rainfall erosion force; extract the soil erodibility status identifier and soil structure stability index of the same spatial unit from the spatial heterogeneous field of soil erodibility; and perform horizontal splicing processing on the rainfall erosion force status identifier, rainfall erosion force peak migration trajectory, soil erodibility status identifier and soil structure stability index of the same spatiotemporal unit to generate a multidimensional coupled input vector for each spatiotemporal unit. The multidimensional coupled input vector includes spatiotemporal coordinate identifier and corresponding feature components.

[0047] Extract S_rain_ij(y) and T_rain_ij for each spatiotemporal unit from the F_rain generated in step S127, where the spatiotemporal unit is determined by the time index y (y=1 to Y) and the spatial index (i, j). Extract S_erode_ij and F_stable_ij with the same spatial index (i, j) from the G_erode generated in step S137. For each spatiotemporal unit, concatenate the above four features laterally to create a one-dimensional vector V_input. The first component of the vector stores S_rain_ij(y), the second component stores the total length T_len_ij of the T_rain_ij polyline, the third component stores S_erode_ij, and the fourth component stores F_stable_ij. Simultaneously, add spatiotemporal coordinate identifiers to this vector, including the time index y, the spatial row index i, and the spatial column index j. The generated one-dimensional vector V_input is the multidimensional coupled input vector for this spatiotemporal unit. Repeat the above extraction and splicing process for all Y×R_num×C_num spatiotemporal units to generate V_input for each spatiotemporal unit.

[0048] Step S142: Input the multidimensional coupled input vector into the input layer of the soil erosion coupled response model. The soil erosion coupled response model includes an input layer, a spatiotemporal attention mapping layer, a stress coupled response layer, and a risk classification output layer. The input layer performs numerical range compression processing on each feature component in the multidimensional coupled input vector of each spatiotemporal unit to generate a compressed value for each feature component. The compressed values ​​of all spatiotemporal units are arranged in spatiotemporal coordinate order to generate a multi-channel spatiotemporal feature tensor. The multi-channel spatiotemporal feature tensor includes a rainfall erosion force state identification channel, a rainfall erosion force peak migration trajectory channel, a soil erodibility state identification channel, and a soil structure stability index channel.

[0049] The soil erosion coupled response model is a pre-built deep learning model. This model comprises four main components: an input layer, a spatiotemporal attention mapping layer, a stress coupled response layer, and a risk classification output layer. All V_input values ​​generated in step S141 are sequentially input into the model's input layer. The input layer performs numerical range compression on each feature component. For S_rain_ij(y), S_rain_norm = S_rain_ij(y) / 4 is calculated. For T_len_ij, the maximum value T_max and minimum value T_min of T_len_ij in all spatiotemporal units are calculated, and T_norm = (T_len_ij - T_min) / (T_max - T_min). For S_erode_ij, S_erode_norm = S_erode_ij / 4 is calculated. For F_stable_ij, F_norm is directly set to F_stable_ij. After the above numerical range compression processing, the feature components of each spatiotemporal unit are converted into four compressed values. The compressed values ​​of all spatiotemporal units are arranged in spatiotemporal coordinate order to generate a four-dimensional tensor X_input. The dimensions of X_input are Y×R_num×C_num×4, where the four channels in the fourth dimension correspond to the rainfall erosivity state indicator channel, the rainfall erosivity peak migration trajectory channel, the soil erosibility state indicator channel, and the soil structure stability index channel, respectively. This four-dimensional tensor X_input is the multi-channel spatiotemporal feature tensor.

[0050] Step S143: The multi-channel spatiotemporal feature tensor is passed to the spatiotemporal attention mapping layer. The spatiotemporal attention mapping layer calculates the similarity of the changing trend of the rainfall erosivity state identifier between each spatiotemporal unit and its temporally neighboring units, and at the same time calculates the spatial distribution similarity of the soil erodibility state identifier between each spatiotemporal unit and its spatially neighboring units. The similarity of the changing trend and the similarity of the spatial distribution are weighted and summed to generate the spatiotemporal neighborhood influence coefficient of each spatiotemporal unit. The spatiotemporal neighborhood influence coefficient of each spatiotemporal unit is multiplied with the multidimensional coupling input vector of the spatiotemporal unit to generate the enhanced coupling vector of each spatiotemporal unit. The enhanced coupling vectors of all spatiotemporal units are rearranged according to the spatiotemporal coordinates to generate the spatiotemporal attention enhanced feature tensor.

[0051] Step S1431: Extract the rainfall erosion force status identifiers of all spatiotemporal units of the rainfall erosion force status identifier channel from the multi-channel spatiotemporal feature tensor. Arrange the rainfall erosion force status identifiers of different time points at the same spatial location into a rainfall erosion force status time sequence for that spatial location according to the time order. Use the difference between the rainfall erosion force status identifiers of adjacent time points in the rainfall erosion force status time sequence as the state change direction parameter of that spatial location at the corresponding time point.

[0052] Data A_rain, representing the rainfall erosivity state identifier channel, is extracted from the multi-channel spatiotemporal feature tensor X_input generated in step S142. The dimensions of A_rain are Y×R_num×C_num, where Y=10, R_num=40, and C_num=80. For each spatial location (i, j), i=1 to R_num and j=1 to C_num, the data for all time points y in A_rain are extracted and arranged in the order of y=1 to Y to obtain the rainfall erosivity state time series A_rain_ij(y) for that spatial location. For each time point y in this time series (y from 2 to Y), the difference in rainfall erosivity state identifier between adjacent time points is calculated to obtain the state change direction parameter D_rain_ij(y)=A_rain_ij(y)-A_rain_ij(y-1). For y=1, there is no previous time, so D_rain_ij(1) is not calculated. The value range of D_rain_ij(y) is -3 to 3, where negative values ​​indicate a state decrease and positive values ​​indicate a state increase.

[0053] Step S1432: For each spatiotemporal unit, extract the state change direction parameter of the spatial location of the spatiotemporal unit, and extract the state change direction parameter of the previous time point adjacent to the time of the spatiotemporal unit. Calculate the cosine similarity value of the two state change direction parameters as the similarity of the change trend of the rainfall erosion force state identifier between the spatiotemporal unit and its time neighboring units.

[0054] For each spatiotemporal unit (i, j, y), where y ranges from 2 to Y-1, extract the state change direction parameter D_rain_ij(y) at time point y for the spatial location of that spatiotemporal unit. Extract the state change direction parameter D_rain_ij(y-1) at time point y-1 for the spatial location of that spatiotemporal unit. Calculate the cosine similarity value C_sim_ij(y) between these two parameters. The value of C_sim_ij(y) ranges from -1 to 1. The closer the value is to 1, the more consistent the state change trends at the two time points; the closer the value is to -1, the opposite the trends. For the boundary spatiotemporal unit between y=1 and y=Y, there is no complete temporal neighborhood, so C_sim_ij(y) is assigned a value of 0.

[0055] Step S1433: Extract the soil erodibility status identifiers of all spatiotemporal units of the soil erodibility status identifier channel from the multi-channel spatiotemporal feature tensor. Arrange the soil erodibility status identifiers at different spatial locations at the same time point into a spatial distribution matrix of soil erodibility status at that time point according to spatial coordinates. Use the difference between the soil erodibility status identifiers of adjacent spatial locations in the spatial distribution matrix of soil erodibility status as the state space gradient vector at the corresponding spatial location at that time point.

[0056] Data B_erode, representing the soil erodibility state, is extracted from the multi-channel spatiotemporal feature tensor X_input generated in step S142. The dimensions of B_erode are Y×R_num×C_num. For each time point y, the data of all spatial locations (i, j) in B_erode are extracted and arranged according to the spatial row index i and spatial column index j to obtain the spatial distribution matrix B_erode_y(i, j) of the soil erodibility state at that time point. The size of this matrix is ​​R_num×C_num. For each internal spatial location (i, j) in the matrix, where i ranges from 2 to R_num-1 and j ranges from 2 to C_num-1, calculate the state difference between this location and its eastern neighbor: G_east = B_erode_y(i, j+1) - B_erode_y(i, j), the state difference between this location and its northern neighbor: G_north = B_erode_y(i-1, j) - B_erode_y(i, j), the state difference between this location and its western neighbor: G_west = B_erode_y(i, j-1) - B_erode_y(i, j), and the state difference between this location and its southern neighbor: G_south = B_erode_y(i+1, j) - B_erode_y(i, j). The four differences are combined into a four-dimensional vector G_ij(y) = [G_east, G_north, G_west, G_south]. This vector is the state-space gradient vector at the corresponding spatial location at that time point. For boundary locations, only the actual adjacent directions are calculated, and the components of missing directions are assigned a value of 0.

[0057] Step S1434: For each spatiotemporal unit, extract the state space gradient vector at the time point where the spatiotemporal unit is located, and extract the state space gradient vectors of the eight adjacent spatiotemporal units in the eight directions that are spatially adjacent to the spatiotemporal unit. Calculate the dot product of the state space gradient vector of the spatiotemporal unit with the state space gradient vector of each adjacent spatiotemporal unit. Sum all the dot product values ​​and divide by the total number of adjacent spatiotemporal units to generate the spatial distribution similarity of the soil erodibility status markers between the spatiotemporal unit and its spatially neighboring units.

[0058] For each spatiotemporal unit (i, j, y), where i ranges from 2 to R_num-1 and j ranges from 2 to C_num-1, extract the state-space gradient vector G_ij(y) at time point y of that spatiotemporal unit. Extract the state-space gradient vectors of the eight adjacent spatiotemporal units in the eight spatially adjacent directions. The spatial indices corresponding to the eight adjacent directions are (i, j+1) for the east, (i-1, j+1) for the northeast, (i-1, j) for the north, (i-1, j-1) for the northwest, (i, j-1) for the west, (i+1, j-1) for the southwest, (i+1, j) for the south, and (i+1, j+1) for the southeast. For each adjacent direction, obtain the state-space gradient vector G_neighbor at the same time point y of that adjacent location. Calculate the dot product between G_ij(y) and each G_neighbor, and sum the eight dot products to obtain the dot product sum P_sum_ij(y). Divide the dot product sum P_sum_ij(y) by the total number of adjacent spatiotemporal units, 8, to obtain the average dot product value P_avg_ij(y) = P_sum_ij(y) / 8. P_avg_ij(y) represents the spatial distribution similarity between the spatiotemporal unit and its spatial neighbors. For spatiotemporal units near the boundary, where the number of neighboring units is less than 8, only the dot product values ​​of the actual neighboring units are accumulated and divided by the actual number of neighboring units.

[0059] Step S1435: Multiply the similarity of the change trend of each spatiotemporal unit by a preset time weight coefficient to generate a time attention component; multiply the spatial distribution similarity of the spatiotemporal unit by a preset spatial weight coefficient to generate a spatial attention component; and sum the time attention component and the spatial attention component to generate the original attention score of the spatiotemporal unit.

[0060] Set the preset time weight coefficient W_time=0.5 and the preset spatial weight coefficient W_space=0.5. For each spatiotemporal unit (i, j, y), obtain its trend similarity C_sim_ij(y) from step S1432 and its spatial distribution similarity P_avg_ij(y) from step S1434. Calculate the temporal attention component A_time_ij(y)=W_time×C_sim_ij(y). Calculate the spatial attention component A_space_ij(y)=W_space×P_avg_ij(y). Summate the temporal attention component and the spatial attention component to generate the original attention score U_raw_ij(y)=A_time_ij(y)+A_space_ij(y). The value of U_raw_ij(y) ranges from -1 to 1.

[0061] Step S1436: Normalize the original attention scores of all spatiotemporal units, divide the original attention score of each spatiotemporal unit by the sum of the original attention scores of all spatiotemporal units to generate the spatiotemporal neighborhood influence coefficient of each spatiotemporal unit, and perform element-wise multiplication of the spatiotemporal neighborhood influence coefficient of each spatiotemporal unit with the multidimensional coupling input vector of that spatiotemporal unit to generate the weighted coupling vector of each spatiotemporal unit.

[0062] Calculate the cumulative sum of the raw attention scores U_raw_ij(y) of all spatiotemporal units: U_sum = Σ_{y=1}^{Y}Σ_{i=1}^{R_num}Σ_{j=1}^{C_num}U_raw_ij(y). For each spatiotemporal unit (i, j, y), calculate its spatiotemporal neighborhood influence coefficient α_ij(y) = U_raw_ij(y) / U_sum. The value of α_ij(y) ranges from negative to positive, and the sum of α_ij(y) of all spatiotemporal units equals 1. Obtain the multidimensional coupled input vector V_input for each spatiotemporal unit from step S141. V_input contains four components [S_rain_ij(y), T_len_ij, S_erode_ij, F_stable_ij]. The spatiotemporal neighborhood influence coefficient α_ij(y) is multiplied element-wise with each component in V_input to generate a weighted coupling vector V_weighted_ij(y)=[α_ij(y)×S_rain_ij(y), α_ij(y)×T_len_ij, α_ij(y)×S_erode_ij, α_ij(y)×F_stable_ij].

[0063] Step S1437: Arrange the weighted coupling vectors of all spatiotemporal units in the original order of time coordinates and spatial coordinates to generate a spatiotemporal attention enhancement feature tensor, and pass the spatiotemporal attention enhancement feature tensor as the output of the spatiotemporal attention mapping layer to the stress coupling response layer.

[0064] The V_weighted_ij(y) of all spatiotemporal units generated in step S1436 are arranged in the original order of time index y (y=1 to Y), spatial row index i (i=1 to R_num), and spatial column index j (j=1 to C_num) to generate a four-dimensional tensor X_attn. The dimensions of X_attn are Y×R_num×C_num×4, where the four channels of the fourth dimension store the weighted rainfall erosivity state identifier, the weighted rainfall erosivity peak migration trajectory length, the weighted soil erosibility state identifier, and the weighted soil structure stability index, respectively. This X_attn is the spatiotemporal attention enhancement feature tensor. This tensor is used as the output of the spatiotemporal attention mapping layer and passed to the stress coupling response layer.

[0065] Step S144: The spatiotemporal attention-enhanced feature tensor is passed to the stress coupling response layer. The stress coupling response layer extracts the stress coupling relationship between different channels in the spatiotemporal attention-enhanced feature tensor. The state matching coefficient between the rainfall erosivity state identifier channel and the soil erodibility state identifier channel is calculated. The directional response coefficient between the rainfall erosivity peak migration trajectory channel and the soil structure stability index channel is calculated. The state matching coefficient and the directional response coefficient are multiplied to generate the initial value of the soil and water loss risk level parameter for each spatiotemporal unit.

[0066] Step S1441: Extract the rainfall erosion force status identifiers of all spatiotemporal units of the rainfall erosion force status identifier channel from the spatiotemporal attention enhancement feature tensor, and convert the rainfall erosion force status identifier of each spatiotemporal unit into a corresponding state value. Convert the first type of erosion force status identifier into a first state value, the second type of erosion force status identifier into a second state value, and the third type of erosion force status identifier into a third state value to generate the rainfall erosion force state value of each spatiotemporal unit.

[0067] Data from the rainfall erosivity state identifier channel is extracted from the spatiotemporal attention-enhanced feature tensor X_attn generated in step S1437, specifically the data A_attn from the first channel of the fourth dimension of X_attn. A_attn has dimensions Y×R_num×C_num, and each element A_attn_ij(y) stores the weighted compressed value of the rainfall erosivity state identifier. First, A_attn_ij(y) is multiplied by 4 to restore the original rainfall erosivity state identifier value range of 1 to 4, resulting in A_orig_ij(y) = A_attn_ij(y) × 4. The values ​​of A_orig_ij(y) are 1, 2, 3, and 4, corresponding to the first, second, third, and fourth types of erosivity state identifiers, respectively. The first type of erosivity state identifier (value 1) is converted to the first state value A_val_ij(y)=1; the second type of erosivity state identifier (value 2) is converted to the second state value A_val_ij(y)=2; the third type of erosivity state identifier (value 3) is converted to the third state value A_val_ij(y)=3; and the fourth type of erosivity state identifier (value 4) is converted to the fourth state value A_val_ij(y)=4. The rainfall erosivity state value A_val_ij(y) for each spatiotemporal unit is then generated.

[0068] Step S1442: Extract the soil erodibility status identifiers of all spatiotemporal units of the soil erodibility status identifier channel from the spatiotemporal attention enhancement feature tensor, and convert the soil erodibility status identifier of each spatiotemporal unit into the corresponding state value. In the order of soil erodibility status identifiers from the first category to the fourth category, convert them into an increasing sequence of state values ​​to generate the soil erodibility status value of each spatiotemporal unit.

[0069] Data from the soil erodibility status identifier channel is extracted from X_attn generated in step S1437, specifically the data B_attn from the third channel of the fourth dimension of X_attn. The dimension of B_attn is Y×R_num×C_num, and each element B_attn_ij(y) stores the weighted compressed value of the soil erodibility status identifier. First, B_attn_ij(y) is multiplied by 4 to restore the original soil erodibility status identifier value range of 1 to 4, resulting in B_orig_ij(y) = B_attn_ij(y) × 4. The values ​​of B_orig_ij(y) are 1, 2, 3, and 4, corresponding to the first type of soil erodibility status identifier (low erodibility), the second type (medium erodibility), the third type (high erodibility), and the fourth type (extremely high erodibility), respectively. Following the order from category one to category four, the soil erodibility status identifier for category one is converted to the first state value B_val_ij(y)=1, category two to the second state value B_val_ij(y)=2, category three to the third state value B_val_ij(y)=3, and category four to the fourth state value B_val_ij(y)=4. The soil erodibility status value B_val_ij(y) for each spatiotemporal unit is then generated.

[0070] Step S1443: For each spatiotemporal unit, multiply the rainfall erosivity state value and the soil erodibility state value of the spatiotemporal unit to generate a first product value, add the sum of the squares of the rainfall erosivity state value and the soil erodibility state value of the spatiotemporal unit to generate a first sum value, and divide the first product value by the first sum value to generate the original state matching coefficient of the spatiotemporal unit.

[0071] For each spatiotemporal unit (i, j, y), obtain A_val_ij(y) from step S1441 and B_val_ij(y) from step S1442. Calculate the first product value P_mul_ij(y) = A_val_ij(y) × B_val_ij(y). Calculate the square of A_val_ij(y) A_sq_ij(y) = A_val_ij(y)^2, calculate the square of B_val_ij(y) B_sq_ij(y) = B_val_ij(y)^2, and add them together to obtain the first sum value S_sum_ij(y) = A_sq_ij(y) + B_sq_ij(y). Divide the first product value by the first sum value to generate the original state matching coefficient M_raw_ij(y) = P_mul_ij(y) / S_sum_ij(y). The value of M_raw_ij(y) ranges from 0 to 1. When A_val_ij(y) and B_val_ij(y) are equal and neither is zero, M_raw_ij(y) is close to 0.5. When the difference between the two is large, M_raw_ij(y) is smaller.

[0072] Step S1444: Perform nonlinear mapping processing on the original state matching coefficients of all spatiotemporal units, substitute the original state matching coefficients as input parameters into the S-shaped growth curve function for calculation, generate state matching coefficients with values ​​within a preset range, and use the state matching coefficients as the first output parameter of the stress coupling relationship between the rainfall erosivity state identification channel and the soil erodibility state identification channel.

[0073] For each spatiotemporal unit (i, j, y), the M_raw_ij(y) generated in step S1443 is substituted into the sigmoid growth curve function as an input parameter. The sigmoid growth curve function has the form f(x) = 1 / (1 + e^{-k×(x-0.5)}), where k is the curve steepness parameter, k = 10, and e is the natural constant. The state matching coefficient M_final_ij(y) = 1 / (1 + e^{-10×(M_raw_ij(y)-0.5)}) is calculated. This nonlinear mapping maps the range of M_raw_ij(y) from 0 to 1 to a range between 0 and 1, but enhances the distinguishability of intermediate values. M_final_ij(y) is the first output parameter of the stress coupling relationship between the rainfall erosivity state identification channel and the soil erodibility state identification channel.

[0074] Step S1445: Extract the direction vector of the peak migration trajectory of rainfall erosivity for each spatiotemporal unit of the peak migration trajectory of rainfall erosivity from the spatiotemporal attention enhancement feature tensor. At the same time, extract the spatial gradient direction vector of the soil structure stability index for each spatiotemporal unit of the soil structure stability index channel. Calculate the cosine of the angle between the direction vector of the peak migration trajectory of rainfall erosivity and the spatial gradient direction vector of the soil structure stability index for each spatiotemporal unit as the directional consistency coefficient of that spatiotemporal unit.

[0075] The peak migration trajectory data C_attn of rainfall erosivity is extracted from X_attn generated in step S1437. The dimension of C_attn is Y×R_num×C_num. However, C_attn_ij(y) stores the weighted trajectory length and does not contain direction information. Therefore, the direction information is directly extracted from the original T_rain_ij stored in step S127. For each spatiotemporal unit (i, j, y), the overall direction angle θ_ij(y) of the trajectory polyline is extracted from T_rain_ij. The direction angle θ_ij(y) ranges from 0 to 360 degrees and is defined as the angle between the line connecting the trajectory start point and the end point and the due east direction. The direction vector V_rain_ij(y) = [cos(θ_ij(y)), sin(θ_ij(y))] is calculated. Spatial gradient information is extracted from F_stable_ij stored in step S137. For each spatiotemporal unit (i, j, y), calculate the partial derivative of F_stable_ij at that location along the spatial row direction. Partial derivatives in the direction of the column in space Combine these two partial derivatives to form a two-dimensional vector. Calculate the magnitude of the vector. Then calculate the gradient direction vector per unit space. When L_grad is zero, V_soil_ij(y) is assigned the value [0, 0]. The directional consistency coefficient D_cos_ij(y) is calculated as the dot product of V_rain_ij(y) and V_soil_ij(y), i.e., The value of D_cos_ij(y) ranges from -1 to 1, where 1 indicates that the directions are completely the same and -1 indicates that the directions are completely opposite.

[0076] Step S1446: For each spatiotemporal unit, multiply the directional consistency coefficient of the spatiotemporal unit by the value of the soil structure stability index of the spatiotemporal unit to generate a second product value. Add the value of the soil structure stability index of the spatiotemporal unit to a preset constant parameter and take the reciprocal to generate a reciprocal coefficient. Multiply the second product value and the reciprocal coefficient to generate the original directional response coefficient of the spatiotemporal unit.

[0077] For each spatiotemporal unit (i, j, y), obtain D_cos_ij(y) from step S1445. Extract the soil structure stability index channel data D_attn from X_attn generated in step S1437, where D_attn_ij(y) is the weighted soil structure stability index. First, restore D_attn_ij(y) to the original soil structure stability index F_stable_val_ij(y). Since F_stable_ij was compressed but not scaled in step S142, F_stable_val_ij(y) = D_attn_ij(y). Calculate the second product value P_sec_ij(y) = D_cos_ij(y) × F_stable_val_ij(y). Set the preset constant parameter ε = 0.001 and calculate the reciprocal coefficient R_inv_ij(y) = 1 / (F_stable_val_ij(y) + ε). The second product value is multiplied by the reciprocal coefficient to generate the original directional response coefficient R_dir_raw_ij(y) = P_sec_ij(y) × R_inv_ij(y). The range of R_dir_raw_ij(y) depends on D_cos_ij(y) and F_stable_val_ij(y). When F_stable_val_ij(y) is small and D_cos_ij(y) is large, R_dir_raw_ij(y) may be large.

[0078] Step S1447: Normalize the original directional response coefficients of all spatiotemporal units, divide the original directional response coefficient of each spatiotemporal unit by the maximum value of the original directional response coefficients of all spatiotemporal units to generate a directional response coefficient whose value range is within a preset interval, and use this directional response coefficient as the second output parameter of the stress coupling relationship between the peak migration trajectory channel of rainfall erosion force and the soil structure stability index channel.

[0079] Calculate the maximum value of the raw directional responsivity coefficient R_dir_raw_ij(y) for all spatiotemporal units: R_dir_max = max_{y,i,j}R_dir_raw_ij(y). For each spatiotemporal unit (i,j,y), calculate the directional responsivity coefficient R_dir_final_ij(y) = R_dir_raw_ij(y) / R_dir_max. When R_dir_max is zero, assign all R_dir_final_ij(y) values ​​to 0. The value of R_dir_final_ij(y) ranges from 0 to 1. This R_dir_final_ij(y) is the second output parameter representing the stress coupling relationship between the peak rainfall erosivity migration trajectory channel and the soil structure stability index channel.

[0080] Step S145: The initial value of the soil erosion risk level parameter of each spatiotemporal unit is transmitted to the risk classification output layer. The risk classification output layer performs spatiotemporal smoothing filtering on the initial value of the soil erosion risk level parameter of each spatiotemporal unit. The spatiotemporal Gaussian filter kernel is used to perform spatiotemporal neighborhood weighted average calculation on the initial value of the soil erosion risk level parameter to generate the final value of the soil erosion risk level parameter of each spatiotemporal unit. The final value of the soil erosion risk level parameter is divided into three categories: first risk level, second risk level and third risk level according to the preset risk level classification threshold.

[0081] Step S1451: Receive the initial value of the soil erosion risk level parameter for each spatiotemporal unit from the stress coupling response layer, and organize the initial values ​​of the soil erosion risk level parameter for all spatiotemporal units into a three-dimensional risk parameter tensor according to the time dimension and the spatial dimension. The first dimension of the three-dimensional risk parameter tensor is the spatial row coordinate, the second dimension is the spatial column coordinate, and the third dimension is the time coordinate.

[0082] From step S144, generate the initial value of the soil erosion risk level parameter R_risk_initial_ij(y) for each spatiotemporal unit. Organize the R_risk_initial_ij(y) for all spatiotemporal units according to the spatial row index i (i=1 to R_num, R_num=40), spatial column index j (j=1 to C_num, C_num=80), and time index y (y=1 to Y, Y=10). Create a three-dimensional array R_init, with dimensions R_num×C_num×Y, where the first dimension is the spatial row coordinate, the second dimension is the spatial column coordinate, and the third dimension is the time coordinate. The array element R_init(i, j, y) = R_risk_initial_ij(y). This three-dimensional array is the three-dimensional risk parameter tensor.

[0083] Step S1452: Construct a spatiotemporal Gaussian filter kernel function. The spatiotemporal Gaussian filter kernel function includes a spatial Gaussian component and a temporal Gaussian component. The standard deviation parameter of the spatial Gaussian component is used to control the smoothing range of the spatial neighborhood, and the standard deviation parameter of the temporal Gaussian component is used to control the smoothing range of the temporal neighborhood. Multiply the spatial Gaussian component and the temporal Gaussian component to generate the weight coefficient matrix of the spatiotemporal Gaussian filter kernel.

[0084] Construct a spatiotemporal Gaussian filter kernel function. The expression for the spatial Gaussian component is G_s(Δi, Δj)=e^{-(Δi^2+Δj^2) / (2×σ_s^2)}, where Δi is the spatial row offset, Δj is the spatial column offset, and σ_s is the standard deviation parameter of the spatial Gaussian component, set to σ_s=1.5. The expression for the temporal Gaussian component is G_t(Δt)=e^{-Δt^2 / (2×σ_t^2)}, where Δt is the time offset, and σ_t is the standard deviation parameter of the temporal Gaussian component, set to σ_t=1.0. Set the spatial neighborhood window radius W_s=2, meaning the values ​​of Δi and Δj are both in the range of -2 to 2, and set the temporal neighborhood window radius W_t=1, meaning the value of Δt is in the range of -1 to 1. For each combination (Δi, Δj, Δt), calculate G_s(Δi, Δj) and G_t(Δt), and then multiply them to obtain the weight coefficients K(Δi, Δj, Δt) = G_s(Δi, Δj) × G_t(Δt) of the spatiotemporal Gaussian filter kernel. All weight coefficients constitute a three-dimensional weight coefficient matrix with dimensions (2W_s+1) × (2W_s+1) × (2W_t+1) = 5 × 5 × 3.

[0085] Step S1453: Take each spatiotemporal unit in the three-dimensional risk parameter tensor as the current processing center point, and extract the initial values ​​of the soil and water loss risk level parameters corresponding to all spatial units in the spatial neighborhood window and all time points in the temporal neighborhood window with the current processing center point as the center, to form the spatiotemporal neighborhood parameter set of the current processing center point.

[0086] Each cell (i, j, y) in the three-dimensional risk parameter tensor R_init is used as the current processing center point. For this center point, all spatial cells within the spatial neighborhood window are extracted, i.e., the spatial row index range is i-W_s to i+W_s, and the spatial column index range is j-W_s to j+W_s, where W_s=2. All time points within the temporal neighborhood window are extracted, i.e., the time index range is y-W_t to y+W_t, where W_t=1. For center points located on the boundary, when i-W_s is less than 1, the spatial row index range starts from 1; when i+W_s is greater than R_num, the spatial row index range ends at R_num; the spatial column index and time index are processed in the same way. All R_init(i+Δi, j+Δj, y+Δt) extracted from the above range are used to form the spatiotemporal neighborhood parameter set of the current processing center point, where Δi = -W_s to W_s, Δj = -W_s to W_s, Δt = -W_t to W_t, and satisfying that i+Δi is in the range of 1 to R_num, j+Δj is in the range of 1 to C_num, and y+Δt is in the range of 1 to Y.

[0087] Step S1454: Multiply the weight coefficient matrix of the spatiotemporal Gaussian filter kernel with the initial value of each soil and water loss risk level parameter in the spatiotemporal neighborhood parameter set of the current processing center point at corresponding positions to generate the weighted risk parameter value of each spatiotemporal neighborhood position. Sum the weighted risk parameter values ​​of all spatiotemporal neighborhood positions to generate a weighted sum value. Sum all weight coefficients in the weight coefficient matrix of the spatiotemporal Gaussian filter kernel to generate a total weighted sum value.

[0088] For the current processing center point (i, j, y), obtain each R_init(i+Δi, j+Δj, y+Δt) in its spatiotemporal neighborhood parameter set. Obtain the spatiotemporal Gaussian filter kernel weight coefficients K(Δi, Δj, Δt) constructed in step S1452. Multiply each K(Δi, Δj, Δt) by the corresponding R_init(i+Δi, j+Δj, y+Δt) to obtain the weighted risk parameter value W_val(Δi, Δj, Δt) = K(Δi, Δj, Δt) × R_init(i+Δi, j+Δj, y+Δt). Sum the W_val(Δi, Δj, Δt) values ​​at all valid (Δi, Δj, Δt) positions to obtain the weighted sum value W_sum = ΣW_val(Δi, Δj, Δt). Simultaneously, the weights K(Δi, Δj, Δt) at all valid (Δi, Δj, Δt) positions are summed to obtain the total weight value K_sum = ΣK(Δi, Δj, Δt). Note that for the boundary center point, only the actual neighboring positions are included in the calculation, so K_sum may be less than the total weight of the entire window.

[0089] Step S1455: Divide the weighted sum by the total weight to generate the spatiotemporally smoothed risk parameter value of the current processing center point. Use the spatiotemporally smoothed risk parameter value as the intermediate value of the soil and water loss risk level parameter of the current processing center point. Then, repeat the steps of extracting the spatiotemporal neighborhood parameter set, weighted multiplication and accumulation, and weighted average operation for each spatiotemporal unit in the three-dimensional risk parameter tensor to generate the intermediate value of the soil and water loss risk level parameter of all spatiotemporal units.

[0090] For the current processing center point (i, j, y), calculate the spatiotemporally smoothed risk parameter value R_smooth_ij(y) = W_sum / K_sum. Use this R_smooth_ij(y) as the intermediate value of the soil erosion risk level parameter for the current processing center point. Repeat steps S1453 and S1454 sequentially for each spatiotemporal unit (i, j, y) in the three-dimensional risk parameter tensor R_init to generate R_smooth_ij(y) for all spatiotemporal units.

[0091] Step S1456: Perform boundary effect correction processing on the intermediate values ​​of soil erosion risk level parameters of all spatiotemporal units. Replace the intermediate values ​​of soil erosion risk level parameters of spatiotemporal units located on the spatial boundary of the target watershed with the result of one-sided neighborhood smoothing. Replace the intermediate values ​​of soil erosion risk level parameters of spatiotemporal units located at the start and end points of the time series with the result of forward or backward smoothing to generate the corrected intermediate values ​​of soil erosion risk level parameters.

[0092] For spatiotemporal units located on spatial boundaries, i.e., i=1 or i=R_num or j=1 or j=C_num, the smoothing result in step S1455 uses fewer neighborhood units. The above boundary units are re-smoothed, but a one-sided neighborhood approach is adopted, that is, only the neighborhood units located inside the boundary are used, and the units that do not exist outside the boundary are not used. Specifically, for the spatial row boundary i=1, the spatial neighborhood window only takes Δi=0 to W_s; for the spatial row boundary i=R_num, the spatial neighborhood window only takes Δi=-W_s to 0; the same applies to the spatial column boundaries j=1 and j=C_num. The R_smooth_boundary_ij(y) of these boundary units is recalculated and replaced with the original R_smooth_ij(y). For spatiotemporal units located at the time series starting point y=1, the forward smoothing result is used for replacement, that is, only the neighborhood time points from Δt=0 to W_t are used to calculate R_smooth_ij(1). For the spatiotemporal unit located at the time series termination point y=Y, backward smoothing is used for replacement, that is, only the neighborhood time points from Δt=-W_t to 0 are used to calculate R_smooth_ij(Y). After the above boundary effect correction, the intermediate value of the corrected soil erosion risk level parameter R_smooth_corrected_ij(y) is generated.

[0093] Step S1457: The intermediate values ​​of the corrected soil erosion risk level parameters are reorganized according to the original spatiotemporal coordinates to generate the final value of the soil erosion risk level parameters for each spatiotemporal unit, and the final value of the soil erosion risk level parameters is used as the output result of the risk classification output layer.

[0094] The R_smooth_corrected_ij(y) generated in step S1456 is reorganized according to the original spatiotemporal coordinates (i, j, y) to generate the final value of the soil and water loss risk level parameter for each spatiotemporal unit, R_risk_final_ij(y) = R_smooth_corrected_ij(y). This R_risk_final_ij(y) is the final soil and water loss risk level parameter output by the risk classification output layer for each spatiotemporal unit.

[0095] Step S150: Extract the key water and soil conservation monitoring area of ​​the target watershed based on the water and soil loss risk level parameters and the boundary sequence of water and soil loss sensitive areas in the spatiotemporal coupling field of water and soil loss risk, and generate intelligent water and soil conservation monitoring results including the boundary polygon of the monitoring area and the dynamic monitoring frequency of the monitoring area.

[0096] Step S151: Extract the final value of the soil erosion risk level parameter of all spatial units on each time section from the spatiotemporal coupling field of soil erosion risk. Mark the spatial unit on each time section whose final value of the soil erosion risk level parameter belongs to the third risk level as the third risk unit of that time section, and mark the spatial unit on each time section whose final value of the soil erosion risk level parameter belongs to the second risk level as the second risk unit of that time section.

[0097] From the spatiotemporal coupled field of soil erosion risk generated in step S145, extract the final value R_risk_final_ij(y) of the soil erosion risk level parameter for all spatial units (i, j) on each time segment y (y=1 to Y, Y=10). For each time segment y, traverse all spatial units. When R_risk_final_ij(y) is greater than or equal to the preset second risk level threshold H_risk2=0.6, mark the spatial unit as a third risk unit and set the flag bit Flag_3rd_ij(y)=1. When R_risk_final_ij(y) is greater than or equal to the preset first risk level threshold H_risk1=0.3 and less than H_risk2=0.6, mark the spatial unit as a second risk unit and set the flag bit Flag_2nd_ij(y)=1. Spatial units with a risk level lower than the second risk level are not marked.

[0098] Step S152: Perform spatial connectivity analysis on the third risk unit on each time segment, divide spatially adjacent third risk units into the same third risk connected region, calculate the total number of third risk units contained in each third risk connected region as the region size parameter of the third risk connected region, and retain the third risk connected regions whose region size parameter exceeds the preset first size threshold as candidate third risk regulatory regions.

[0099] For each time segment y, all spatial units with Flag_3rd_ij(y)=1 are extracted. Spatial connectivity analysis is performed using the eight-connectivity criterion: two spatial units are considered spatially adjacent when the absolute difference between their row indices and column indices is less than or equal to 1. A breadth-first search algorithm is used to traverse all third-risk units, aggregating spatially adjacent third-risk units into a single third-risk connected region. For each aggregated third-risk connected region, the total number of third-risk units contained within it (Count_3rd_region) is counted. A preset first scale threshold S_th1=10 is set, and third-risk connected regions with Count_3rd_region greater than S_th1 are retained as candidate third-risk regulatory regions, and the spatial coordinate sequence of all spatial units contained within these regions is recorded.

[0100] Step S153: Perform spatial connectivity analysis on the second risk unit on each time segment, divide the spatially adjacent second risk units into the same second risk connectivity region, calculate the total number of second risk units contained in each second risk connectivity region as the region size parameter of the second risk connectivity region, and retain the second risk connectivity region whose region size parameter exceeds the preset second size threshold as the candidate second risk regulatory region.

[0101] For each time segment y, extract all spatial units with Flag_2nd_ij(y)=1. Using the same eight-connectivity criterion and breadth-first search algorithm as in step S152, aggregate spatially adjacent second-risk units into the same second-risk connected region. For each aggregated second-risk connected region, count the total number of second-risk units contained in the region, Count_2nd_region. Set a preset second scale threshold S_th2=20, retain second-risk connected regions with Count_2nd_region greater than S_th2 as candidate second-risk monitoring regions, and record the spatial coordinate sequence of all spatial units contained in the region.

[0102] Step S154: Extract the spatial regions that are marked as candidate third risk regulatory regions on multiple consecutive time segments as stable third risk regions; extract the spatial regions that change from candidate second risk regulatory regions to candidate third risk regulatory regions on multiple consecutive time segments as risk level transition regions; and extract the spatial regions that change from candidate third risk regulatory regions to candidate second risk regulatory regions on multiple consecutive time segments as risk level decline regions.

[0103] For each spatial location (i, j), iterate through time segments y=1 to Y, recording the regulatory zone type to which the location belongs on each time segment. Set the continuous judgment window length L_cont=3. When a spatial location belongs to a candidate third-risk regulatory zone on three consecutive time segments, mark the spatial location as a member of a stable third-risk zone. Aggregate the spatial connectivity of all spatial locations marked as members of stable third-risk zones to obtain several stable third-risk zones. When a spatial location belongs to a candidate second-risk regulatory zone on time segment y, and belongs to a candidate third-risk regulatory zone on both time segments y+1 and y+2, mark the spatial location as a member of a risk level transition zone. Aggregate the spatial connectivity of all spatial locations marked as members of a risk level transition zone to obtain several risk level transition zones. When a spatial location belongs to a candidate third-risk regulatory zone on time segment y, and belongs to a candidate second-risk regulatory zone on both time segments y+1 and y+2, mark the spatial location as a member of a risk level decline zone. Spatial connectivity aggregation is performed on the spatial locations of all members marked as risk level decline areas to obtain several risk level decline areas.

[0104] Step S155: For each stable third risk area, extract the spatial unit coordinate sequence through which the spatial boundary of the stable third risk area passes on the latest time section, convert the spatial unit coordinate sequence into a regulatory area boundary polygon composed of multiple boundary inflection point coordinates, and use the regulatory area boundary polygon as the area boundary polygon of the stable third risk area.

[0105] For each stable third risk region extracted in step S154, obtain the spatial coordinates (X_ij, Y_ij) of all member spatial units of that region on the latest time segment y=Y. The convex hull algorithm is used to calculate the convex polygon boundary of these spatial coordinates. The specific process of the convex hull algorithm is as follows: First, find the point with the smallest spatial row coordinate among all points as the starting point. Then, sort the other points according to their polar angles. Next, use the Graham scan method to traverse the sorted points, maintaining a stack structure. When the direction formed by the new point and the top point of the stack is clockwise, pop the top element of the stack; otherwise, push the new point onto the stack. The remaining points in the stack are the vertices of the convex hull. Connect the spatial coordinates of the convex hull vertices in counter-clockwise order to form a closed polygon, which is the boundary polygon of the monitored area. Record the vertex coordinate sequence of this polygon.

[0106] Step S156: For each risk level transition region, extract the number of time segments from the first appearance of the second risk level to the first appearance of the third risk level as the level transition speed parameter of the risk level transition region. Risk level transition regions with a level transition speed parameter less than a preset speed threshold are marked as fast transition regions, and the regulatory area boundary polygon of the fast transition region is extracted.

[0107] For each risk level transition region extracted in step S154, obtain the record of regulatory area type changes for all member spatial locations within that region on the time cross-section. For each spatial location within that region, find the time cross-section y_first_2nd where the location was first marked as a candidate second-risk regulatory area, and the time cross-section y_first_3rd where it was first marked as a candidate third-risk regulatory area. Calculate the transition time difference Δy = y_first_3rd - y_first_2nd for that location. Take the average of the transition time differences for all spatial locations within that region, Δy_avg, as the risk level transition speed parameter for that risk level transition region. Set a preset speed threshold V_th = 2. When Δy_avg is less than V_th, the risk level transition region is marked as a fast transition region. For each fast transition region, use the same convex hull algorithm as in step S155 to extract the regulatory area boundary polygon of that region on the latest time cross-section y = Y.

[0108] Step S157: Determine the dynamic monitoring frequency of each regulatory area based on the area size of the boundary polygon of each regulatory area and the final value of the soil erosion risk level parameter corresponding to that regulatory area. Set the dynamic monitoring frequency of regulatory areas whose area exceeds the preset first area threshold and whose final value of the soil erosion risk level parameter exceeds the preset third risk level threshold as the first monitoring frequency. Set the dynamic monitoring frequency of regulatory areas whose area does not exceed the preset first area threshold but whose final value of the soil erosion risk level parameter exceeds the preset third risk level threshold as the second monitoring frequency. Set the dynamic monitoring frequency of regulatory areas whose final value of the soil erosion risk level parameter belongs to the second risk level as the third monitoring frequency.

[0109] For each stable third risk region generated in step S155 and each fast transition region generated in step S156, calculate the area Area_region enclosed by the boundary polygon of its regulatory area. The area calculation adopts the shoelace formula: for the polygon vertex coordinate sequence (X_1, Y_1), (X_2, Y_2), ..., (X_n, Y_n), the area Area_region = 0.5 × |Σ_{p=1}^{n}(X_p × Y_{p+1} - X_{p+1} × Y_p)|, where when p=n, X_{n+1}=X_1 and Y_{n+1}=Y_1. For each regulatory region, calculate the average value R_avg_region of R_risk_final_ij(Y) of all spatial units within the region at the latest time segment y=Y. Set the preset first area threshold A_th1 = 1,000,000 square meters and the preset third risk level threshold R_th3 = 0.6. When Area_region is greater than A_th1 and R_avg_region is greater than R_th3, the dynamic monitoring frequency of this regulatory area is set to the first monitoring frequency F1 = 1 time / day. When Area_region is less than or equal to A_th1 but R_avg_region is greater than R_th3, the dynamic monitoring frequency of this regulatory area is set to the second monitoring frequency F2 = 1 time / 3 days. For the portion of the candidate second risk regulatory area generated in step S153 that is not marked as a risk level transition area, the average value of R_risk_final_ij(Y) in this region, R_avg_region_2nd, is calculated. When R_avg_region_2nd belongs to the second risk level range (0.3 to 0.6), the dynamic monitoring frequency of this regulatory area is set to the third monitoring frequency F3 = 1 time / week.

[0110] Step S158: Associate and store all the boundary polygons of the monitored areas with their corresponding dynamic monitoring frequencies to generate intelligent water and soil conservation monitoring results containing the sequence of boundary polygons of the monitored areas and the dynamic monitoring frequency corresponding to each monitored area.

[0111] The boundary polygons of the monitored areas of the stable third risk zone generated in step S155 are associated with the corresponding dynamic monitoring frequencies determined in step S157 to generate a first type of data record. Each record includes an area identifier, a sequence of vertex coordinates of the boundary polygons, and the dynamic monitoring frequency. The boundary polygons of the monitored areas of the rapidly transitioning areas generated in step S156 are associated with the corresponding dynamic monitoring frequencies determined in step S157 to generate a second type of data record. The boundary polygons of the portion of the candidate second risk monitored areas not marked as risk level transition areas generated in step S153 are associated with the third monitoring frequencies determined in step S157 to generate a third type of data record. All three types of data records are organized into the aforementioned intelligent monitoring results for soil and water conservation.

[0112] For example, the method may further include: step S210, collecting vegetation multispectral remote sensing reflectance data and measured runoff sediment concentration data from hydrological stations in the target watershed, wherein the vegetation multispectral remote sensing reflectance data includes near-infrared band reflectance sequences and red band reflectance sequences for different vegetation types, and the measured runoff sediment concentration data from hydrological stations includes runoff sediment concentration and runoff velocity values ​​from different hydrological stations.

[0113] In this embodiment, the same target watershed is used as the scenario, and multispectral remote sensing reflectance data of vegetation in the target watershed are collected. Remote sensing reflectance products of the target watershed during the growing season (May to September each year) are acquired using a multispectral satellite sensor. These products include a near-infrared reflectance sequence (NIR(t)) and a red light reflectance sequence (RED(t)), with a spatial resolution of 30 meters by 30 meters and a temporal resolution of 16 days, acquiring data for five consecutive growing seasons, each containing 10 temporal phases. L hydrological stations are set up within the target watershed, L=15, each equipped with an automatic water level gauge and a suspended sediment sampler. At each hydrological station, whenever a rainfall-runoff event occurs, the automatic water level gauge records the runoff depth value H_l(t) at a frequency of once per hour, and the runoff velocity value V_l(t) is calculated by combining this data with the water level-discharge curve of that station. Meanwhile, the suspended sediment sampler collects water samples every hour. After filtration, drying, and weighing, the runoff sediment concentration C_l(t) is obtained, in grams per liter. All data collection processes above involve blurring the precise location information of the hydrological stations, retaining only the relative positions within the watershed.

[0114] Step S220: Perform vegetation canopy physiological state inversion processing on the vegetation multispectral remote sensing reflectance data. Calculate the normalized vegetation index time series for each pixel location using the ratio of the near-infrared band reflectance sequence to the red band reflectance sequence. Calculate the vegetation cover time series based on the normalized vegetation index time series, and generate vegetation cover dynamic change characteristic parameters for each pixel location. The vegetation cover dynamic change characteristic parameters include the interannual variation rate of vegetation cover and the maximum vegetation cover during the growing season.

[0115] The vegetation multispectral remote sensing reflectance data of the target watershed is divided into the same grid system as in step S110. Each grid cell corresponds to multiple pixel locations, and the average value of all pixels within each grid cell is taken as the reflectance value of that grid cell. For each grid cell (i, j) and each time phase τ (τ = 1 to Tau, Tau = 50, i.e., 5 growing seasons multiplied by 10 time phases per growing season), NIR_ij(τ) is extracted from the near-infrared reflectance sequence, and RED_ij(τ) is extracted from the red reflectance sequence. The normalized vegetation index NDVI_ij(τ) = (NIR_ij(τ) - RED_ij(τ)) / (NIR_ij(τ) + RED_ij(τ)) is calculated. The NDVI_ij(τ) of all time phases is arranged in chronological order to generate the normalized vegetation index time series of that grid cell. For each growing season g (g=1 to 5), extract 10 NDVI values ​​within that growing season and calculate the maximum NDVI value NDVI_max_ij(g) within that growing season. Then calculate the vegetation cover FVC_ij(g) = (NDVI_max_ij(g) - NDVI_soil) / (NDVI_veg - NDVI_soil), where NDVI_soil is the reference value of the Normalized Difference Vegetation Index for bare soil (0.1), and NDVI_veg is the reference value of the Normalized Difference Vegetation Index for dense vegetation (0.8). Divide the difference in vegetation cover between two consecutive growing seasons by the time interval (1 year) to obtain the interannual rate of change of vegetation cover ΔFVC_ij = (FVC_ij(g+1) - FVC_ij(g)) / 1. The maximum vegetation cover FVC_max_ij(g) within each growing season is taken as the peak vegetation cover of that growing season. Repeat the above calculations for each grid cell to generate vegetation cover dynamic change characteristic parameters for each grid cell, including the interannual vegetation cover change rate ΔFVC_ij and the maximum vegetation cover during the growing season FVC_max_ij(g).

[0116] Step S230: Perform runoff and sediment yield relationship analysis on the measured runoff sediment concentration data of the hydrological station. Multiply the runoff sediment concentration value by the corresponding runoff velocity value as the instantaneous unit sediment yield. Integrate the unit sediment yield of all instantaneous periods of the same hydrological station during a preset rainfall event to generate the total sediment yield of the rainfall event. Calculate the ratio of the total sediment yield to the total runoff of the rainfall event to generate the average sediment concentration of the hydrological station during a sub-rainfall event. Calculate the average value of the average sediment concentration of all sub-rainfall events within a preset time period to generate the annual average sediment yield load intensity parameter of the hydrological station.

[0117] For each hydrological station l (l=1 to L, L=15), the runoff sediment load C_l(t) and runoff velocity V_l(t) for that station are obtained from step S210, where t is the hourly time point during the rainfall event. For a complete rainfall event, the start time is T_s_event and the end time is T_e_event. The total sediment load of the rainfall event S_load_l = Σ_{t=T_s_event}^{T_e_event}(C_l(t)×V_l(t)×Δt), where Δt = 1 hour. The total runoff of the rainfall event Q_total_l = Σ_{t=T_s_event}^{T_e_event}(V_l(t)×A_l×Δt), where A_l is the catchment area of ​​the hydrological station. Calculate the average sediment concentration of this rainfall event: C_avg_event_l = S_load_l / Q_total_l. Statistically calculate the average sediment concentration of all rainfall events over the past five consecutive hydrological years, and then calculate their arithmetic mean to obtain the annual average sediment load intensity parameter for this hydrological station: Y_load_l = (Σ_{event=1}^{Event_total}C_avg_event_l) / Event_total.

[0118] Step S240: Construct a spatial response relationship network between vegetation cover dynamic change characteristic parameters and annual average sediment load intensity parameters. Use spatial kriging interpolation to interpolate the annual average sediment load intensity parameters of discrete hydrological stations to generate a spatial distribution map of annual average sediment load intensity covering the entire spatial range of the target watershed. Calculate the correlation coefficient between the interannual variation rate of vegetation cover at each pixel location and the interpolated annual average sediment load intensity parameter corresponding to that pixel location. Use this correlation coefficient as the response coefficient of vegetation change to sediment load intensity.

[0119] Step S220 yields the ΔFVC_ij for each grid cell, and step S230 yields the Y_load_l and its spatial coordinates (X_l, Y_l) for L hydrological stations. Ordinary Kriging interpolation is used to interpolate the discrete Y_load_l to generate a spatial distribution map of the annual average sediment load intensity Y_load_ij covering the entire target watershed. The ordinary Kriging interpolation process is as follows: First, calculate the semi-variogram value γ(h) = 0.5 × [Y_load_l - Y_load_m]^2 between all pairs of hydrological stations, where h is the distance between the two stations. Fit a semi-variogram model (e.g., a spherical model). Then, for each grid node (i, j), calculate the weight coefficient λ_l of the surrounding hydrological stations based on the semi-variogram model to minimize the interpolation variance. Finally, calculate Y_load_ij = Σλ_l × Y_load_l. For each grid cell (i, j), calculate the Pearson correlation coefficient R_corr_ij between ΔFVC_ij and Y_load_ij, where R_corr_ij = Cov(ΔFVC_ij, Y_load_ij) / (σ_ΔFVC × σ_Yload), and Cov is the covariance, σ is the standard deviation. R_corr_ij is used as the response coefficient β_ij to the influence of vegetation change on sediment yield intensity for that grid cell.

[0120] Step S250: Based on the response coefficient of vegetation change to sediment yield intensity, correct the soil and water loss risk level parameters in the spatiotemporal coupling field of soil and water loss risk. Extract the soil and water loss risk level parameters of each spatial unit from the spatiotemporal coupling field of soil and water loss risk. Multiply the soil and water loss risk level parameters with the response coefficient of vegetation change to sediment yield intensity of the spatial unit to generate the corrected soil and water loss risk level parameters. Remap the corrected soil and water loss risk level parameters of all spatial units to the corresponding spatial unit locations to generate the corrected spatiotemporal distribution field of soil and water loss risk.

[0121] Extract the final value of the soil erosion risk level parameter R_risk_final_ij(y) for each spatiotemporal unit from the spatiotemporal coupled field of soil erosion risk generated in step S150. Obtain the response coefficient β_ij for each spatial unit (i, j) from step S240. For each spatiotemporal unit, calculate the corrected soil erosion risk level parameter R_corrected_ij(y) = R_risk_final_ij(y) × β_ij. Reorganize the R_corrected_ij(y) of all spatiotemporal units according to the time dimension y and the spatial dimension (i, j) to generate a three-dimensional tensor R_corrected, where the dimension of R_corrected is Y × R_num × C_num. This three-dimensional tensor is the corrected spatiotemporal distribution field of soil erosion risk.

[0122] Step S260: Extract spatial units with corrected soil erosion risk level parameters belonging to the third risk level from the corrected spatiotemporal distribution field of soil erosion risk. Aggregate spatial units that are spatially adjacent and whose corrected soil erosion risk level parameters belong to the third risk level into third risk patches. Calculate the geometric center coordinates and area of ​​each third risk patch. Mark third risk patches with an area greater than a preset patch area threshold as candidate areas for soil and water conservation engineering intervention.

[0123] Extract the corrected risk level parameter R_corrected_ij(Y) for the latest time segment y=Y (i.e., the 10th hydrological year) from R_corrected generated in step S250. Set the third risk level threshold R_th3=0.6, and mark spatial units with R_corrected_ij(Y)≥R_th3 as high-risk units. Use the four-connectivity criterion to perform spatial clustering on high-risk units, aggregating spatially adjacent high-risk units into the same third-risk patch. For each third-risk patch, calculate its geometric center coordinates (X_center, Y_center), where X_center is equal to the average of the X_ij coordinates of all spatial units within the patch, and Y_center is equal to the average of the Y_ij coordinates of all spatial units within the patch. Calculate the patch area Area_patch, which is equal to the number of spatial units within the patch multiplied by the area D^2 of each unit. Set a preset patch area threshold A_patch=1,000,000 square meters, and mark the third risk patch with Area_patch>A_patch as a candidate area for soil and water conservation engineering intervention.

[0124] Step S270: Extract the interannual variation rate of vegetation cover from the original vegetation cover dynamic change characteristic parameters corresponding to each candidate area for soil and water conservation engineering intervention. Determine the vegetation cover change trend of the candidate area for soil and water conservation engineering intervention based on the positive or negative sign of the interannual variation rate of vegetation cover. If the interannual variation rate of vegetation cover is negative, it is determined to be a vegetation degradation trend. If the interannual variation rate of vegetation cover is positive, it is determined to be a vegetation recovery trend.

[0125] For each candidate area for soil and water conservation engineering intervention marked in step S260, the average value ΔFVC_avg is extracted from ΔFVC_ij of all grid cells within the candidate area stored in step S220. The sign of ΔFVC_avg is determined: when ΔFVC_avg < 0, the vegetation cover change trend of the candidate area is determined to be a vegetation degradation trend; when ΔFVC_avg > 0, the vegetation cover change trend of the candidate area is determined to be a vegetation recovery trend; when ΔFVC_avg = 0, it is determined to be no change trend.

[0126] Step S280: Search for the nearest suitable vegetation restoration project location within a preset distance range from the geometric center coordinates of the third risk patch in the candidate area for soil and water conservation engineering intervention with vegetation degradation trend. Connect the coordinates of the searched suitable vegetation restoration project location with the geometric center coordinates of the third risk patch in the candidate area to generate a vegetation restoration project implementation path.

[0127] For each candidate area for soil and water conservation engineering intervention identified as having a vegetation degradation trend, its geometric center coordinates (X_center, Y_center) are obtained. A search radius R_search_veg = 5000 meters is set, and suitable locations for vegetation restoration projects are searched within this radius. Suitable locations for vegetation restoration projects are defined as locations with a slope less than 15 degrees, soil texture of sandy loam or loam, and less than 1000 meters from existing roads. From all suitable locations found, the location (X_restore, Y_restore) closest to (X_center, Y_center) is selected. The straight line connecting (X_center, Y_center) and (X_restore, Y_restore) is used as the implementation path for the vegetation restoration project, and the coordinate sequence along this path is recorded.

[0128] Step S290: Search for the nearest existing soil and water conservation engineering facility location within a preset distance range from the geometric center coordinates of the third risk patch in the candidate area of ​​the soil and water conservation engineering intervention area with the searched existing soil and water conservation engineering facility location coordinates to generate an engineering facility maintenance and inspection path by connecting the geometric center coordinates of the third risk patch in the candidate area.

[0129] For each candidate area for soil and water conservation engineering intervention identified as having a vegetation restoration trend, its geometric center coordinates (X_center, Y_center) are obtained. A search radius R_search_maintain = 8000 meters is set, and the locations of existing soil and water conservation engineering facilities are searched within this radius. Existing soil and water conservation engineering facilities include horizontal terraces, silt-retaining dams, and vegetation filter belts, and their location coordinates are obtained from the soil and water conservation engineering database of the target watershed. From all the existing facilities found, the location (X_facility, Y_facility) closest to (X_center, Y_center) is selected. The straight line connecting (X_center, Y_center) and (X_facility, Y_facility) is used as the maintenance and inspection path for the engineering facilities, and the coordinate sequence along this path is recorded.

[0130] Step S2100: Overlay the vegetation restoration project implementation path and the engineering facility maintenance and inspection path with the modified spatiotemporal distribution field of soil and water loss risk to obtain a target management scheme that includes the geometric center coordinates of the third risk patch, the coordinate sequence of the vegetation restoration project implementation path, the coordinate sequence of the engineering facility maintenance and inspection path, and the vegetation cover change trend of the corresponding candidate area. Then, add the target management scheme to the intelligent monitoring results of soil and water conservation.

[0131] The coordinate sequences of vegetation restoration project implementation paths for each candidate area of ​​vegetation degradation trend generated in step S280, and the coordinate sequences of engineering facility maintenance and inspection paths for each candidate area of ​​vegetation restoration trend generated in step S290, are spatially overlaid with the corrected spatiotemporal distribution field R_corrected of soil and water loss risk generated in step S250. For each candidate area, a data record is created, containing the geometric center coordinates (X_center, Y_center) of the third risk patch in the candidate area, the vegetation cover change trend indicator (degradation or restoration), and the corresponding engineering path coordinate sequence. The data records of all candidate areas are organized into a list structure, which is the target management plan. This target management plan is appended to the intelligent soil and water conservation monitoring results generated in step S150 to form intelligent soil and water conservation monitoring results containing information on the monitored area and engineering management information.

[0132] For example, the method may further include: step S310, acquiring long-term land use change data and precipitation data from meteorological observation stations for the target watershed, wherein the long-term land use change data includes land use classification maps for multiple historical periods, and the precipitation data from meteorological observation stations includes annual precipitation sequences and annual precipitation erosivity sequences from multiple meteorological stations.

[0133] In this embodiment, the same target watershed is used as the scene to acquire long-term land use change data for the target watershed. Land use classification maps for the past 20 years are extracted from historical remote sensing image classification products, with each map defined in 5-year periods, for a total of 4 periods, corresponding to years t1, t2, t3, and t4. Each land use classification map has a spatial resolution of 30 meters by 30 meters and includes types such as forest land, grassland, cultivated land, construction land, water area, and unused land. Simultaneously, precipitation observation data from M meteorological stations within and around the target watershed are acquired, M=25. Each meteorological station contains the annual precipitation sequence P_annual_m(yr) and the annual precipitation erosivity sequence R_annual_m(yr) for the past 20 years, yr=1 to 20.

[0134] Step S320: Perform land cover transfer matrix analysis on the long-term land use change data, extract the land use type code for each pixel location in the land use classification map of two adjacent historical periods, combine the land use type code of the previous period with the land use type code of the next period to generate the land use transfer type code for each pixel location, perform statistical frequency calculation on the land use transfer type codes of all pixel locations, and generate a land use transfer matrix. The land use transfer matrix includes the transfer area parameters from forest land to cultivated land, the transfer area parameters from grassland to construction land, and the transfer area parameters from cultivated land to forest land.

[0135] The land use classification map is resampled to the same grid system as in step S110, and the land use type of each grid cell is the dominant type within that cell. For two adjacent periods t_p and t_{p+1}, where p = 1 to 3, the land use type code L_code_p_ij for each grid cell (i, j) in period t_p and the land use type code L_code_{p+1}_ij in period t_{p+1} are extracted. The two codes are combined to generate the transition type code T_code = L_code_p_ij × 10 + L_code_{p+1}_ij. The frequency of each T_code in all grid cells is counted, and multiplied by the area D^2 of each grid cell to obtain the area of ​​each transition type. The statistical results are organized into a U×U transition matrix, where U is the total number of land use types, U = 6, and the element in the a-th row and b-th column of the matrix represents the area transitioning from land type a to land type b. Three key parameters are extracted from this transfer matrix: the transfer area parameter A_forest_to_crop from forest land to cultivated land, the transfer area parameter A_grass_to_urban from grassland to construction land, and the transfer area parameter A_crop_to_forest from cultivated land to forest land.

[0136] Step S330: Perform spatiotemporal variation characteristic analysis on the precipitation data of the meteorological observation stations, calculate the multi-year average value and coefficient of variation of the annual precipitation sequence for each meteorological station, calculate the multi-year average value and coefficient of variation of the annual precipitation erosivity sequence for each meteorological station, multiply the multi-year average value and coefficient of variation of the annual precipitation sequence to generate the annual precipitation variability index, and multiply the multi-year average value and coefficient of variation of the annual precipitation erosivity sequence to generate the precipitation erosivity variability index.

[0137] For each meteorological station m, m=1 to M, M=25, extract its annual precipitation sequence P_annual_m(yr), yr=1 to 20. Calculate the multi-year average P_mean_m=(Σ_{yr=1}^{20}P_annual_m(yr)) / 20. Calculate the standard deviation of the annual precipitation sequence P_std_m=[(Σ_{yr=1}^{20}(P_annual_m(yr)-P_mean_m)^2) / 20]^{1 / 2}. Calculate the coefficient of variation P_cv_m=P_std_m / P_mean_m. Calculate the annual precipitation variability index P_var_m=P_mean_m×P_cv_m. Extract the annual precipitation erosivity sequence R_annual_m(yr), yr=1 to 20 for that station. Calculate the multi-year average value R_mean_m = (Σ_{yr=1}^{20}R_annual_m(yr)) / 20. Calculate the standard deviation of the precipitation erosivity series R_std_m = [(Σ_{yr=1}^{20}(R_annual_m(yr)-R_mean_m)^2) / 20]^{1 / 2}. Calculate the coefficient of variation R_cv_m = R_std_m / R_mean_m. Calculate the precipitation erosivity variability index R_var_m = R_mean_m × R_cv_m.

[0138] Step S340: Extract the transfer area parameter from forest land to cultivated land from the land use transfer matrix, extract the latest period's boundary sequence of soil and water loss sensitive areas from the spatiotemporal coupling field of soil and water loss risk, calculate the spatial overlap area ratio parameter between the spatial range involved in the transfer area parameter from forest land to cultivated land and the spatial range surrounded by the boundary sequence of soil and water loss sensitive areas, and use the spatial overlap area ratio parameter as the spatial impact intensity coefficient of land use change on soil and water loss risk.

[0139] Extract A_forest_to_crop from the transition matrix generated in step S320. This parameter corresponds to the total area converted from forest land to cultivated land. Extract the boundary sequence of soil erosion sensitive areas at the latest time section y=Y (the 10th hydrological year) from the spatiotemporal coupled field of soil erosion risk generated in step S150. The spatial range enclosed by this boundary sequence is denoted as region Z_sensitive. Calculate the spatial range involved in A_forest_to_crop, i.e., the set of grid cells where forest land to cultivated land conversion occurs, denoted as region Z_transition. Calculate the spatial overlap area A_overlap = Z_transition ∩ Z_sensitive. Calculate the spatial overlap area ratio parameter R_overlap = A_overlap / A_sensitive, where A_sensitive is the total area of ​​Z_sensitive. Use R_overlap as the spatial influence intensity coefficient γ_landuse of land use change on soil erosion risk.

[0140] Step S350: Extract the precipitation erosivity variability index of all meteorological stations from the precipitation data of the meteorological observation stations. Use spatial interpolation method to interpolate the precipitation erosivity variability index of discrete meteorological stations to generate a spatial distribution map of precipitation erosivity variability index covering the entire spatial range of the target watershed. Extract the average value of precipitation erosivity variability index of all pixel positions within the spatial range corresponding to the boundary sequence of each soil erosion sensitive area from the spatial distribution map of precipitation erosivity variability index. Use this average value as the dynamic intensity coefficient of precipitation erosivity of the soil erosion sensitive area.

[0141] From step S330, the precipitation erosivity variability index R_var_m and its spatial coordinates (X_m, Y_m) for each meteorological station m are obtained. Using the same inverse distance weighted interpolation method as in step S135, the discrete R_var_m is interpolated to generate a spatial distribution map R_var_ij of the precipitation erosivity variability index covering the entire target watershed. For each soil erosion-sensitive area generated in step S150, the R_var_ij values ​​of all grid cells within the spatial range enclosed by the boundary sequence of that area are extracted, and their arithmetic mean R_var_avg is calculated. This average value is used as the precipitation erosivity dynamic intensity coefficient δ_precip for that soil erosion-sensitive area.

[0142] Step S360: Construct a comprehensive stress assessment model for the spatial impact intensity coefficient of land use change on soil erosion risk and the dynamic intensity coefficient of precipitation erosion. Multiply the spatial impact intensity coefficient and the dynamic intensity coefficient of precipitation erosion to generate a comprehensive stress index for each soil erosion sensitive area. Normalize the comprehensive stress index of all soil erosion sensitive areas so that the value of the comprehensive stress index is within a preset range.

[0143] For each soil erosion-sensitive area generated in step S150, obtain its corresponding spatial influence intensity coefficient γ_landuse (all sensitive areas share the same γ_landuse value) and precipitation erosion dynamic intensity coefficient δ_precip. Calculate the comprehensive stress index U_stress = γ_landuse × δ_precip. Normalize the U_stress of all soil erosion-sensitive areas, find the maximum value U_max and minimum value U_min of U_stress in all areas, and calculate the normalized comprehensive stress index U_norm = (U_stress - U_min) / (U_max - U_min), where U_norm ranges from 0 to 1.

[0144] Step S370: Dynamically adjust the soil and water loss risk level parameters of each soil and water loss sensitive area based on the comprehensive stress index. Extract the average soil and water loss risk level parameters of each soil and water loss sensitive area from the spatiotemporal coupling field of soil and water loss risk. Convert the comprehensive stress index of each soil and water loss sensitive area into a dimensionless adjustment coefficient. Perform a weighted operation on the average soil and water loss risk level parameters and the adjustment coefficient to generate the adjusted soil and water loss risk level parameters. Reassign the adjusted soil and water loss risk level parameters to all spatial units within the soil and water loss sensitive area to generate the stress-corrected spatiotemporal coupling field of soil and water loss risk.

[0145] From the spatiotemporal coupled field of soil erosion risk generated in step S150, extract the average value R_avg_s of R_risk_final_ij(Y) for all spatiotemporal units (latest time section y=Y) within each soil erosion sensitive area s. Use the U_norm generated in step S360 for this area as an adjustment coefficient. Calculate the adjusted soil erosion risk level parameter R_adjusted_s = R_avg_s × (1 + U_norm). Assign R_adjusted_s to all spatial units within the sensitive area to obtain the stress-corrected risk level parameter R_corrected_final_ij for each spatial unit. Organize the R_corrected_final_ij of all spatial units into a two-dimensional matrix according to spatial coordinates to generate the stress-corrected spatiotemporal coupled field of soil erosion risk.

[0146] Step S380: Extract spatial units whose adjusted soil erosion risk level parameters are greater than the preset third risk level threshold from the spatiotemporal coupled field of soil erosion risk after stress correction. Aggregate spatially adjacent spatial units whose adjusted soil erosion risk level parameters are greater than the preset third risk level threshold to generate target stress risk patches. Calculate the area parameter and the average value of the adjusted soil erosion risk level parameters for each target stress risk patch. Mark the target stress risk patches whose area parameters are greater than the preset area threshold and whose average value of the adjusted soil erosion risk level parameters is greater than the preset level threshold as priority treatment areas.

[0147] From the stress-corrected spatiotemporal coupled field of soil erosion risk generated in step S370, spatial cells with R_corrected_final_ij > R_th3 (R_th3 = 0.6) are extracted. The four-connectivity criterion is used to perform spatial clustering on these spatial cells, aggregating spatially adjacent cells into the same target stress risk patch. For each target stress risk patch, its area Area_riskpatch = number of grid cells within the patch × D^2 is calculated, and the average value R_avg_patch of all R_corrected_final_ij values ​​within the patch is calculated. A preset area threshold A_risk = 2,000,000 square meters and a preset level threshold R_risk_th = 0.7 are set, and target stress risk patches with Area_riskpatch > A_risk and R_avg_patch > R_risk_th are marked as priority treatment areas.

[0148] Step S390: For each priority governance area, extract the land use transfer type coding sequence of the priority governance area in the past preset time period from the long-term land use change data, count the number of transfer events from forest land to cultivated land in the land use transfer type coding sequence, and use the number of transfer events as the human activity disturbance intensity parameter of the priority governance area.

[0149] For each priority governance area marked in step S380, the land use transfer type coding sequence of all grid cells within the spatial range of the priority governance area between adjacent periods is extracted from the four periods of land use classification data stored in step S320. The total number of occurrences, F_disturb, of the transfer events from forest land to cultivated land (i.e., L_code_p_ij represents forest land and L_code_{p+1}_ij represents cultivated land in the transfer type coding) in this sequence is counted. F_disturb is used as the human activity disturbance intensity parameter for the priority governance area.

[0150] Step S3100: Determine the comprehensive governance priority of each priority governance area based on the human activity interference intensity parameter and the comprehensive stress index of that priority governance area. Multiply the human activity interference intensity parameter and the comprehensive stress index to generate a comprehensive governance priority score. Sort the comprehensive governance priority scores of all priority governance areas in descending order to generate a comprehensive governance priority ranking list of priority governance areas.

[0151] For each priority governance area marked in step S380, obtain its human activity disturbance intensity parameter F_disturb (from step S390) and normalized comprehensive stress index U_norm (from step S360). Calculate the comprehensive governance priority score Q_priority = F_disturb × U_norm. Sort all priority governance areas in descending order of Q_priority to generate a comprehensive governance priority ranking list of priority governance areas.

[0152] Step S3110: Integrate the spatial boundary coordinates of the priority treatment area, the comprehensive treatment priority score, and the comprehensive treatment priority ranking list to generate a comprehensive treatment plan for soil and water loss risk under stress, and add the comprehensive treatment plan to the intelligent monitoring results of soil and water conservation.

[0153] For each priority treatment area marked in step S380, the spatial unit coordinate sequence traversed by its spatial boundary is extracted, and this sequence is converted into a priority treatment area boundary polygon composed of multiple boundary inflection point coordinates. The boundary polygon of each priority treatment area, its comprehensive treatment priority score Q_priority, and its sequence number in the comprehensive treatment priority ranking list are associated and stored to generate a data record. All priority treatment area data records are organized into a list structure, which represents the comprehensive treatment plan for soil erosion risk under stress. This comprehensive treatment plan is appended to the soil and water conservation intelligent monitoring results generated in step S2100 to form the final soil and water conservation intelligent monitoring results.

[0154] In one exemplary embodiment, a smart water and soil conservation monitoring system based on multi-source data fusion is provided. This system can be a terminal, server, etc., and its internal structure diagram can be as follows: Figure 2As shown, the system may include a processor, memory, input / output interface, communication interface, display unit, and input device. The processor, memory, and input / output interface are connected via a system bus, and the communication interface, display unit, and input device are also connected to the system bus via the input / output interface. The processor provides computing and control capabilities. The memory includes a non-volatile storage medium and internal memory. The non-volatile storage medium stores the operating system and computer programs. The internal memory provides an environment for the operation of the operating system and computer programs in the non-volatile storage medium. The input / output interface is used for exchanging information between the processor and external devices. The communication interface is used for wired or wireless communication with external terminals; wireless communication can be achieved through Wi-Fi, mobile cellular networks, near-field communication, or other technologies. When the computer program is executed by the processor, it implements a smart water and soil conservation monitoring method based on multi-source data fusion. The display unit is used to form a visually visible image and can be a display screen, projection device, or virtual reality imaging device. The display screen can be an LCD screen or an e-ink screen. The input device can be a touch layer covering the display screen, or a button, trackball, or touchpad set on the shell of the intelligent water and soil conservation monitoring system based on multi-source data fusion, or an external keyboard, touchpad, or mouse, etc.

[0155] 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 smart monitoring method for soil and water conservation based on multi-source data fusion, characterized in that, The method includes: The remote sensing inversion data of rainfall and the spatial distribution data of soil properties of the target watershed are collected. The remote sensing inversion data of rainfall includes the time series of rainfall intensity and the spatial distribution map of rainfall accumulation in multiple grid cells. The spatial distribution data of soil properties includes soil texture type parameters and soil organic matter content values ​​of multiple soil sampling points. The rainfall intensity time series and rainfall accumulation spatial distribution map in the rainfall remote sensing inversion data are subjected to spatiotemporal evolution analysis of rainfall erosivity to generate a spatiotemporal dynamic field of rainfall erosivity in the target watershed. The spatiotemporal dynamic field of rainfall erosivity includes the rainfall erosivity status identifier and rainfall erosivity peak migration trajectory for each spatiotemporal unit. The soil texture type parameter and soil organic matter content value in the spatial distribution data of soil properties are subjected to spatial differentiation analysis of soil erodibility to generate a spatial heterogeneous field of soil erodibility for the target watershed. The spatial heterogeneous field of soil erodibility includes the soil erodibility status identifier and soil structure stability index of each spatial unit. The spatiotemporal dynamic field of rainfall erosion force and the spatial heterogeneous field of soil erodibility are input into a pre-constructed water and soil loss coupled response model for spatiotemporal superposition and stress coupling processing to generate a spatiotemporal coupled field of water and soil loss risk for the target watershed. The spatiotemporal coupled field of water and soil loss risk includes water and soil loss risk level parameters and water and soil loss sensitive area boundary sequence for each spatiotemporal unit. Based on the soil and water loss risk level parameters and the boundary sequence of soil and water loss sensitive areas in the spatiotemporal coupling field of the soil and water loss risk, the key monitoring areas for soil and water conservation in the target watershed are extracted, and intelligent monitoring results for soil and water conservation, including the boundary polygons of the monitoring areas and the dynamic monitoring frequency of the monitoring areas, are generated.

2. The intelligent monitoring method for soil and water conservation based on multi-source data fusion according to claim 1, characterized in that, The step of performing spatiotemporal evolution analysis of rainfall erosivity on the rainfall intensity time series and rainfall accumulation spatial distribution map in the rainfall remote sensing inversion data to generate the spatiotemporal dynamic field of rainfall erosivity in the target watershed includes: The spatial coordinates of each grid cell, the rainfall intensity time series corresponding to the grid cell, and the rainfall accumulation value in the spatial distribution map of the rainfall accumulation corresponding to the grid cell are extracted from the rainfall remote sensing inversion data. The spatial coordinates of each grid cell are associated with the rainfall intensity time series and rainfall accumulation value of the grid cell and stored to generate a gridded set of rainfall parameters containing the spatial coordinates of the grid cells and the corresponding rainfall parameter series. The rainfall intensity time series of each grid cell in the rainfall parameter gridded set is processed by rainfall event segmentation. The time period in the rainfall intensity time series where the continuous rainfall intensity value exceeds the preset rainfall intensity threshold is divided into an independent rainfall event. The start time, end time and peak rainfall intensity value of each independent rainfall event are recorded to generate a rainfall event feature set for each grid cell. Extract the rainfall intensity time series of each independent rainfall event from the rainfall event feature set of each grid cell, calculate the rainfall erosivity index of each independent rainfall event based on the rainfall intensity time series, and sum the rainfall erosivity indices of all independent rainfall events in the same grid cell to generate the time-period cumulative rainfall erosivity index of the grid cell. Spatial distribution features of the time-period cumulative rainfall erosivity index of all grid units in the gridded set of rainfall parameters are extracted. The absolute value of the difference between the time-period cumulative rainfall erosivity index of adjacent grid units is calculated. The boundary between adjacent grid units whose absolute difference exceeds a preset difference threshold is marked as the spatial abrupt boundary of rainfall erosivity. The spatial location coordinate sequence of all spatial abrupt boundary of rainfall erosivity is recorded to generate a spatial differentiation boundary map of rainfall erosivity. For each grid cell, the cumulative rainfall erosivity index of all neighboring grid cells within a preset neighborhood range is extracted. The statistical average of these cumulative rainfall erosivity indices is calculated as the neighborhood average rainfall erosivity index of the grid cell. The cumulative rainfall erosivity index of each grid cell is compared with a preset erosivity classification threshold set. A classification code is assigned to the grid cell according to the threshold range it is in, and the rainfall erosivity status identifier of the grid cell is generated. The calculation of the cumulative rainfall erosivity index and the allocation of rainfall erosivity status identifiers are repeatedly performed on the gridded set of rainfall parameters obtained at different time sampling points. A spatial distribution map of rainfall erosivity status identifiers corresponding to each time sampling point is generated. The spatial distribution maps of rainfall erosivity status identifiers of adjacent time sampling points are superimposed and compared. The grid cells in which the rainfall erosivity status identifiers undergo category transformation are extracted as rainfall erosivity status transformation units. The spatial coordinate sequence of these rainfall erosivity status transformation units is recorded. The spatial coordinate sequence of the rainfall erosivity state transition unit corresponding to each time sampling point is connected in chronological order to generate the rainfall erosivity peak migration trajectory, which reflects the spatial transition trajectory of the rainfall erosivity state identifier. The rainfall erosivity state identifier and the rainfall erosivity peak migration trajectory of each spatiotemporal unit are organized according to spatiotemporal coordinates to generate a spatiotemporal dynamic field of rainfall erosivity containing the rainfall erosivity state identifier and the rainfall erosivity peak migration trajectory of each spatiotemporal unit.

3. The intelligent monitoring method for soil and water conservation based on multi-source data fusion according to claim 1, characterized in that, The step of performing spatial differentiation analysis on soil erodibility parameters and soil organic matter content values ​​in the spatial distribution data of soil properties to generate a spatial heterogeneous field of soil erodibility for the target watershed includes: The spatial coordinates of each soil sampling point, the soil texture type parameter corresponding to the soil sampling point, and the soil organic matter content value corresponding to the soil sampling point are extracted from the spatial distribution data of soil properties. The spatial coordinates of each soil sampling point are associated with the soil texture type parameter and the soil organic matter content value of the soil sampling point and stored in a related manner to generate a set of discrete points of soil properties containing the spatial coordinates of the soil sampling points and the corresponding soil property parameters. The soil texture type parameters in the set of discrete soil attribute points are converted into type codes. Sandy loam is encoded as the first texture type code, loam is encoded as the second texture type code, clay loam is encoded as the third texture type code, and clay is encoded as the fourth texture type code, generating the soil texture type code value for each soil sampling point. Soil organic matter content values ​​are extracted from the set of discrete soil attribute points. The soil organic matter content values ​​are compared with a preset organic matter content grading threshold sequence. The grading code corresponding to the grading interval into which the soil organic matter content value falls is used as the soil organic matter grading code value for that soil sampling point. The preset organic matter content grading threshold sequence includes a first grading threshold, a second grading threshold, and a third grading threshold. Sampling points with soil organic matter content values ​​less than the first grading threshold are coded as the first organic matter grading code. Sampling points with soil organic matter content values ​​between the first and second grading thresholds are coded as the second organic matter grading code. Sampling points with soil organic matter content values ​​between the second and third grading thresholds are coded as the third organic matter grading code. Sampling points with soil organic matter content values ​​greater than the third grading threshold are coded as the fourth organic matter grading code. The soil texture type code value and soil organic matter classification code value of each soil sampling point are combined and encoded. The soil texture type code value and soil organic matter classification code value are concatenated into strings to generate the comprehensive soil erodibility code value of each soil sampling point. A mapping relationship table between the comprehensive soil erodibility code value and the preset soil erodibility status identifier is established. The soil erodibility status identifier corresponding to the comprehensive soil erodibility code value of each soil sampling point is queried according to the mapping relationship table. Spatial interpolation is performed on the soil erodibility status markers of all soil sampling points. The inverse distance weighted interpolation method is used to calculate the interpolation results of the soil erodibility status markers at each grid node in the target watershed, generating a spatial distribution map of soil erodibility status markers covering the entire spatial range of the target watershed. The spatial distribution map of soil erodibility status markers includes the spatial coordinates of each grid node and the interpolation results of the soil erodibility status markers corresponding to that grid node. For each grid node in the spatial distribution map of soil erodibility status identifiers, the interpolation results of soil erodibility status identifiers of all grid nodes within the preset neighborhood window of that grid node are extracted. The statistical variance of these soil erodibility status identifier interpolation results is calculated as the soil erodibility spatial variation parameter of that grid node. The soil erodibility spatial variation parameters of all grid nodes are then organized according to spatial coordinates to generate a soil erodibility spatial variation distribution map. The interpolation result of the soil erodibility status identifier of each grid node is associated with the soil erodibility spatial variability parameter of that grid node and stored to generate the soil erodibility status identifier and soil erodibility spatial variability parameter of each grid node. The soil erodibility status identifier of each grid node is used as the soil erodibility status identifier of that grid node. The soil erodibility spatial variability parameter of each grid node is normalized after taking the reciprocal to generate the soil structure stability index of that grid node. The soil erodibility status identifier and soil structure stability index of each spatial unit are organized according to spatial coordinates to generate a soil erodibility spatial heterogeneous field containing the soil erodibility status identifier and soil structure stability index of each spatial unit.

4. The intelligent monitoring method for soil and water conservation based on multi-source data fusion according to claim 1, characterized in that, The process of inputting the spatiotemporal dynamic field of rainfall erosion and the spatial heterogeneous field of soil erodibility into a pre-constructed water and soil erosion coupled response model for spatiotemporal superposition and stress coupling to generate the spatiotemporal coupled field of water and soil erosion risk for the target watershed includes: The rainfall erosion force status identifier and rainfall erosion force peak migration trajectory of each spatiotemporal unit are extracted from the spatiotemporal dynamic field of rainfall erosion force. The soil erodibility status identifier and soil structure stability index of the same spatial unit are extracted from the spatial heterogeneous field of soil erodibility. The rainfall erosion force status identifier, rainfall erosion force peak migration trajectory, soil erodibility status identifier and soil structure stability index of the same spatiotemporal unit are horizontally spliced ​​to generate a multidimensional coupled input vector for each spatiotemporal unit. The multidimensional coupled input vector includes spatiotemporal coordinate identifier and corresponding feature components. The multidimensional coupled input vector is input into the input layer of the soil erosion coupled response model. The soil erosion coupled response model includes an input layer, a spatiotemporal attention mapping layer, a stress coupled response layer, and a risk classification output layer. The input layer performs numerical range compression processing on each feature component in the multidimensional coupled input vector of each spatiotemporal unit to generate a compressed value for each feature component. The compressed values ​​of all spatiotemporal units are arranged in spatiotemporal coordinate order to generate a multi-channel spatiotemporal feature tensor. The multi-channel spatiotemporal feature tensor includes a rainfall erosivity state identifier channel, a rainfall erosivity peak migration trajectory channel, a soil erosibility state identifier channel, and a soil structure stability index channel. The multi-channel spatiotemporal feature tensor is passed to the spatiotemporal attention mapping layer. The spatiotemporal attention mapping layer calculates the trend similarity of the rainfall erosivity state identifier between each spatiotemporal unit and its temporally neighboring units, and simultaneously calculates the spatial distribution similarity of the soil erodibility state identifier between each spatiotemporal unit and its spatially neighboring units. The trend similarity and spatial distribution similarity are weighted and summed to generate the spatiotemporal neighborhood influence coefficient of each spatiotemporal unit. The spatiotemporal neighborhood influence coefficient of each spatiotemporal unit is multiplied with the multidimensional coupling input vector of that spatiotemporal unit to generate the enhanced coupling vector of each spatiotemporal unit. The enhanced coupling vectors of all spatiotemporal units are rearranged according to spatiotemporal coordinates to generate the spatiotemporal attention enhanced feature tensor. The spatiotemporal attention-enhanced feature tensor is passed to the stress coupling response layer. The stress coupling response layer extracts the stress coupling relationship between different channels in the spatiotemporal attention-enhanced feature tensor. The state matching coefficient between the rainfall erosivity state identifier channel and the soil erodibility state identifier channel is calculated. The directional response coefficient between the rainfall erosivity peak migration trajectory channel and the soil structure stability index channel is calculated. The state matching coefficient and the directional response coefficient are multiplied to generate the initial value of the soil and water loss risk level parameter for each spatiotemporal unit. The initial value of the soil erosion risk level parameter of each spatiotemporal unit is passed to the risk classification output layer. The risk classification output layer performs spatiotemporal smoothing filtering on the initial value of the soil erosion risk level parameter of each spatiotemporal unit. The spatiotemporal Gaussian filter kernel is used to perform spatiotemporal neighborhood weighted average calculation on the initial value of the soil erosion risk level parameter to generate the final value of the soil erosion risk level parameter of each spatiotemporal unit. The final value of the soil erosion risk level parameter is divided into three categories: first risk level, second risk level and third risk level according to the preset risk level classification threshold. Spatial connectivity region extraction is performed on the final values ​​of soil erosion risk level parameters of all spatial units at each time segment. Spatial units that are spatially adjacent and whose final values ​​of soil erosion risk level parameters belong to the second or third risk level are divided into the same soil erosion sensitive area. The spatial unit coordinate sequence through which the boundary of each soil erosion sensitive area passes is calculated as the soil erosion sensitive area boundary sequence of that soil erosion sensitive area. The final values ​​of soil erosion risk level parameters of each spatiotemporal unit and the soil erosion sensitive area boundary sequences of each soil erosion sensitive area are organized according to spatiotemporal coordinates to generate a spatiotemporal coupled field of soil erosion risk containing the final values ​​of soil erosion risk level parameters of each spatiotemporal unit and the soil erosion sensitive area boundary sequences of each time segment.

5. The intelligent monitoring method for soil and water conservation based on multi-source data fusion according to claim 4, characterized in that, The process involves passing the multi-channel spatiotemporal feature tensor to the spatiotemporal attention mapping layer, calculating the trend similarity of rainfall erosivity state indicators between each spatiotemporal unit and its temporally neighboring units, and simultaneously calculating the spatial distribution similarity of soil erodibility state indicators between each spatiotemporal unit and its spatially neighboring units. A weighted sum of the trend similarity and spatial distribution similarity is then performed to generate the spatiotemporal neighborhood influence coefficient for each spatiotemporal unit, including: The rainfall erosion force status identifiers of all spatiotemporal units of the rainfall erosion force status identifier channel are extracted from the multi-channel spatiotemporal feature tensor. The rainfall erosion force status identifiers of different time points at the same spatial location are arranged into a rainfall erosion force status time sequence of that spatial location in chronological order. The difference between the rainfall erosion force status identifiers of adjacent time points in the rainfall erosion force status time sequence is used as the state change direction parameter of that spatial location at the corresponding time point. For each spatiotemporal unit, the state change direction parameter of the spatial location of the spatiotemporal unit is extracted, and the state change direction parameter of the previous time point adjacent to the time of the spatiotemporal unit is extracted. The cosine similarity value of the two state change direction parameters is calculated as the similarity of the change trend of the rainfall erosivity state identifier between the spatiotemporal unit and its time neighboring units. Soil erodibility status identifiers of all spatiotemporal units of the soil erodibility status identifier channel are extracted from the multi-channel spatiotemporal feature tensor. Soil erodibility status identifiers of different spatial locations at the same time point are arranged into a soil erodibility status spatial distribution matrix at that time point according to spatial coordinates. The difference between soil erodibility status identifiers of adjacent spatial locations in the soil erodibility status spatial distribution matrix is ​​used as the state space gradient vector at the corresponding spatial location at that time point. For each spatiotemporal unit, the state space gradient vector at the time point of the spatiotemporal unit is extracted, and the state space gradient vectors of the eight adjacent spatiotemporal units in the eight directions adjacent to the spatiotemporal unit are extracted. The dot product of the state space gradient vector of the spatiotemporal unit and the state space gradient vector of each adjacent spatiotemporal unit is calculated. All dot product values ​​are summed and divided by the total number of adjacent spatiotemporal units to generate the spatial distribution similarity of the soil erodibility status identifier between the spatiotemporal unit and its spatial neighbors. The temporal attention component is generated by multiplying the similarity of the change trend of each spatiotemporal unit by a preset time weight coefficient, and the spatial attention component is generated by multiplying the spatial distribution similarity of the spatiotemporal unit by a preset spatial weight coefficient. The original attention score of the spatiotemporal unit is generated by summing the temporal attention component and the spatial attention component. The original attention scores of all spatiotemporal units are normalized. The original attention score of each spatiotemporal unit is divided by the sum of the original attention scores of all spatiotemporal units to generate the spatiotemporal neighborhood influence coefficient of each spatiotemporal unit. The spatiotemporal neighborhood influence coefficient of each spatiotemporal unit is then multiplied element-wise with the multidimensional coupling input vector of that spatiotemporal unit to generate the weighted coupling vector of each spatiotemporal unit. The weighted coupling vectors of all spatiotemporal units are arranged in the original order of time and space coordinates to generate a spatiotemporal attention enhancement feature tensor. This spatiotemporal attention enhancement feature tensor is then passed as the output of the spatiotemporal attention mapping layer to the stress coupling response layer.

6. The intelligent monitoring method for soil and water conservation based on multi-source data fusion according to claim 4, characterized in that, The process of passing the spatiotemporal attention-enhanced feature tensor to the stress coupling response layer, extracting the stress coupling relationship between different channels in the spatiotemporal attention-enhanced feature tensor through the stress coupling response layer, calculating the state matching coefficient between the rainfall erosivity state identifier channel and the soil erodibility state identifier channel, and calculating the directional response coefficient between the rainfall erosivity peak migration trajectory channel and the soil structure stability index channel includes: The rainfall erosion force status identifiers of all spatiotemporal units of the rainfall erosion force status identifier channel are extracted from the spatiotemporal attention enhancement feature tensor. The rainfall erosion force status identifiers of each spatiotemporal unit are converted into corresponding status values. The first type of erosion force status identifiers are converted into first status values, the second type of erosion force status identifiers are converted into second status values, and the third type of erosion force status identifiers are converted into third status values, thereby generating the rainfall erosion force status value of each spatiotemporal unit. Soil erodibility status identifiers of all spatiotemporal units of the soil erodibility status identifier channel are extracted from the spatiotemporal attention-enhanced feature tensor. The soil erodibility status identifier of each spatiotemporal unit is converted into a corresponding state value. The soil erodibility status identifiers are converted into an increasing sequence of state values ​​in order from the first category to the fourth category to generate the soil erodibility status value of each spatiotemporal unit. For each spatiotemporal unit, the rainfall erosivity state value and the soil erodibility state value of the spatiotemporal unit are multiplied to generate a first product value. The sum of the squares of the rainfall erosivity state value and the soil erodibility state value of the spatiotemporal unit is added to generate a first sum value. The first product value is divided by the first sum value to generate the original state matching coefficient of the spatiotemporal unit. The original state matching coefficients of all spatiotemporal units are subjected to nonlinear mapping processing. The original state matching coefficients are used as input parameters and substituted into the S-shaped growth curve function for calculation to generate state matching coefficients with values ​​within a preset range. These state matching coefficients are then used as the first output parameter of the stress coupling relationship between the rainfall erosivity state identification channel and the soil erodibility state identification channel. The direction vector of the peak migration trajectory of rainfall erosivity in each spatiotemporal unit of the peak migration trajectory of rainfall erosivity is extracted from the spatiotemporal attention enhancement feature tensor. At the same time, the spatial gradient direction vector of the soil structure stability index of each spatiotemporal unit of the soil structure stability index is extracted. For each spatiotemporal unit, the cosine value of the angle between the direction vector of the peak migration trajectory of rainfall erosivity and the spatial gradient direction vector of the soil structure stability index is calculated as the direction consistency coefficient of the spatiotemporal unit. For each spatiotemporal unit, the directional consistency coefficient of the spatiotemporal unit is multiplied by the value of the soil structure stability index of the spatiotemporal unit to generate a second product value. The value of the soil structure stability index of the spatiotemporal unit is added to a preset constant parameter and the reciprocal is taken to generate a reciprocal coefficient. The second product value and the reciprocal coefficient are multiplied to generate the original directional response coefficient of the spatiotemporal unit. The original directional response coefficients of all spatiotemporal units are normalized. The original directional response coefficient of each spatiotemporal unit is divided by the maximum value of the original directional response coefficients of all spatiotemporal units to generate a directional response coefficient with a value range within a preset interval. This directional response coefficient is used as the second output parameter of the stress coupling relationship between the peak migration trajectory channel of rainfall erosion force and the soil structure stability index channel.

7. The intelligent monitoring method for soil and water conservation based on multi-source data fusion according to claim 4, characterized in that, The process of passing the initial value of the soil erosion risk level parameter for each spatiotemporal unit to the risk classification output layer, performing spatiotemporal smoothing filtering on the initial value of the soil erosion risk level parameter for each spatiotemporal unit through the risk classification output layer, and using a spatiotemporal Gaussian filter kernel to perform a spatiotemporal neighborhood weighted average calculation on the initial value of the soil erosion risk level parameter to generate the final value of the soil erosion risk level parameter for each spatiotemporal unit includes: The initial value of the soil erosion risk level parameter of each spatiotemporal unit is received from the stress coupling response layer, and the initial value of the soil erosion risk level parameter of all spatiotemporal units is organized into a three-dimensional risk parameter tensor according to the time dimension and the spatial dimension. The first dimension of the three-dimensional risk parameter tensor is the spatial row coordinate, the second dimension is the spatial column coordinate, and the third dimension is the time coordinate. A spatiotemporal Gaussian filter kernel function is constructed, which includes a spatial Gaussian component and a temporal Gaussian component. The standard deviation parameter of the spatial Gaussian component is used to control the smoothing range of the spatial neighborhood, and the standard deviation parameter of the temporal Gaussian component is used to control the smoothing range of the temporal neighborhood. The spatial Gaussian component and the temporal Gaussian component are multiplied to generate the weight coefficient matrix of the spatiotemporal Gaussian filter kernel. Each spatiotemporal unit in the three-dimensional risk parameter tensor is taken as the current processing center point. The initial values ​​of the soil and water loss risk level parameters corresponding to all spatial units and all time points in the temporal neighborhood window are extracted with the current processing center point as the center, forming the spatiotemporal neighborhood parameter set of the current processing center point. The weight coefficient matrix of the spatiotemporal Gaussian filter kernel is multiplied at corresponding positions with the initial value of each soil and water loss risk level parameter in the spatiotemporal neighborhood parameter set of the current processing center point to generate the weighted risk parameter value of each spatiotemporal neighborhood position. The weighted risk parameter values ​​of all spatiotemporal neighborhood positions are summed to generate a weighted sum value. All weight coefficients in the weight coefficient matrix of the spatiotemporal Gaussian filter kernel are summed to generate a total weighted sum value. Divide the weighted sum by the total weight to generate the spatiotemporally smoothed risk parameter value of the current processing center point. Use this spatiotemporally smoothed risk parameter value as the intermediate value of the soil and water loss risk level parameter of the current processing center point. Then, repeat the steps of extracting the spatiotemporal neighborhood parameter set, weighted multiplication and accumulation, and weighted average operation for each spatiotemporal unit in the three-dimensional risk parameter tensor to generate the intermediate value of the soil and water loss risk level parameter of all spatiotemporal units. Boundary effect correction processing is applied to the intermediate values ​​of soil erosion risk level parameters of all spatiotemporal units. The intermediate values ​​of soil erosion risk level parameters of spatiotemporal units located on the spatial boundary of the target watershed are replaced with one-sided neighborhood smoothing results. The intermediate values ​​of soil erosion risk level parameters of spatiotemporal units located at the start and end points of the time series are replaced with forward or backward smoothing results to generate corrected intermediate values ​​of soil erosion risk level parameters. The intermediate values ​​of the corrected soil erosion risk level parameters are reorganized according to the original spatiotemporal coordinates to generate the final value of the soil erosion risk level parameters for each spatiotemporal unit, and the final value of the soil erosion risk level parameters is used as the output result of the risk classification output layer.

8. The intelligent monitoring method for soil and water conservation based on multi-source data fusion according to claim 1, characterized in that, The process involves extracting key water and soil conservation monitoring areas for the target watershed based on the water and soil erosion risk level parameters and the boundary sequence of water and soil erosion sensitive areas in the spatiotemporal coupled field of water and soil erosion risk, and generating intelligent water and soil conservation monitoring results that include the boundary polygons of the monitoring areas and the dynamic monitoring frequency of the monitoring areas, including: Extract the final value of the soil erosion risk level parameter of all spatial units on each time section from the spatiotemporal coupling field of soil erosion risk. Mark the spatial unit whose final value of the soil erosion risk level parameter belongs to the third risk level on each time section as the third risk unit of that time section, and mark the spatial unit whose final value of the soil erosion risk level parameter belongs to the second risk level on each time section as the second risk unit of that time section. Spatial connectivity analysis is performed on the third risk units at each time segment. Spatially adjacent third risk units are divided into the same third risk connected region. The total number of third risk units contained in each third risk connected region is calculated as the region size parameter of the third risk connected region. Third risk connected regions whose region size parameter exceeds a preset first size threshold are retained as candidate third risk regulatory regions. Spatial connectivity analysis is performed on the second risk unit at each time segment. Spatially adjacent second risk units are divided into the same second risk connectivity region. The total number of second risk units contained in each second risk connectivity region is calculated as the region size parameter of the second risk connectivity region. Second risk connectivity regions whose region size parameter exceeds the preset second size threshold are retained as candidate second risk regulatory regions. Spatial regions marked as candidate third-risk regulatory areas across multiple consecutive time segments are extracted as stable third-risk areas; spatial regions that transition from candidate second-risk regulatory areas to candidate third-risk regulatory areas across multiple consecutive time segments are extracted as risk level transition areas; and spatial regions that transition from candidate third-risk regulatory areas to candidate second-risk regulatory areas across multiple consecutive time segments are extracted as risk level decline areas. For each stable third risk region, extract the spatial unit coordinate sequence through which the spatial boundary of the stable third risk region passes on the latest time section, convert the spatial unit coordinate sequence into a regulatory area boundary polygon composed of multiple boundary inflection point coordinates, and use the regulatory area boundary polygon as the region boundary polygon of the stable third risk region. For each risk level transition region, the number of time segments from the first appearance of the second risk level to the first appearance of the third risk level is extracted as the risk level transition speed parameter of the risk level transition region. Risk level transition regions with a risk level transition speed parameter less than a preset speed threshold are marked as fast transition regions, and the regulatory area boundary polygon of the fast transition region is extracted. The dynamic monitoring frequency of each regulatory area is determined based on the area size of the boundary polygon of each regulatory area and the final value of the soil and water loss risk level parameter corresponding to that regulatory area. The dynamic monitoring frequency of the regulatory area whose area exceeds the preset first area threshold and whose final value of the soil and water loss risk level parameter exceeds the preset third risk level threshold is set as the first monitoring frequency. The dynamic monitoring frequency of the regulatory area whose area does not exceed the preset first area threshold but whose final value of the soil and water loss risk level parameter exceeds the preset third risk level threshold is set as the second monitoring frequency. The dynamic monitoring frequency of the regulatory area whose final value of the soil and water loss risk level parameter belongs to the second risk level is set as the third monitoring frequency. All the boundary polygons of the monitored areas are associated with and stored with their corresponding dynamic monitoring frequencies, generating intelligent water and soil conservation monitoring results that include the sequence of boundary polygons of the monitored areas and the dynamic monitoring frequency corresponding to each monitored area.

9. A smart monitoring system for soil and water conservation based on multi-source data fusion, characterized in that, include: processor; A machine-readable storage medium for storing machine-executable instructions of the processor; The processor is configured to execute the intelligent water and soil conservation monitoring method based on multi-source data fusion as described in any one of claims 1 to 8 by executing the machine-executable instructions.

10. A computer program product, characterized in that, The computer program product includes machine-executable instructions stored in a computer-readable storage medium. The processor of the water and soil conservation intelligent monitoring system based on multi-source data fusion reads the machine-executable instructions from the computer-readable storage medium and executes the machine-executable instructions, causing the water and soil conservation intelligent monitoring system based on multi-source data fusion to perform the water and soil conservation intelligent monitoring method based on multi-source data fusion as described in any one of claims 1 to 8.

Images