A method and system for monitoring soil erosion based on remote sensing technology

Abstract

The application discloses a kind of water and soil loss dynamic monitoring method and system based on remote sensing technology, the method comprises: extracting the monthly dynamic characteristics and the static characteristics in year of monitoring basin, constructs dynamic characteristic tensor and static characteristic tensor;Based on D8 single flow direction algorithm, the basin is gridded into node, with water flow direction as directed edge, constructs the directed acyclic graph of downstream;Dynamic and static characteristic tensor is input to the physical information graph neural network with the directed acyclic graph as topological constraint, and the annual soil erosion modulus of the whole basin is output per grid.The application constructs directed acyclic graph by water flow direction, solves the problem that traditional model ignores hydrological connectivity, and injects topographic and rainfall physics prior into graph attention mechanism in logarithmic bias form.At the same time, by using the compound training strategy of ordering consistency loss, measured correction loss and topographic adaptive spatial smoothing loss, high-precision prediction can still be obtained under the condition of sparse measured site.

Description

Technical Field

[0001] This application relates to the field of remote sensing intelligent monitoring technology for soil and water conservation, and relates to, but is not limited to, a method and system for dynamic monitoring of soil erosion based on remote sensing technology. Background Technology

[0002] Accurate monitoring of soil erosion is crucial for watershed management and ecological restoration. Existing methods have the following shortcomings: (1) Empirical models represented by RUSLE estimate the erosion modulus by independently multiplying factors such as rainfall erosivity, soil erodibility, topography, and vegetation cover. They ignore the nonlinear interaction between factors and the hydrological connectivity effect between upstream and downstream areas, resulting in a relatively reasonable spatial distribution trend but a large absolute numerical deviation.

[0003] (2) Traditional machine learning and CNN methods treat each pixel as an independent sample or only consider regular rectangular neighborhoods, which cannot model the directional spatial relationships formed along the water flow path in the watershed.

[0004] (3) Although standard graph neural networks (such as the GAT attention network) can process graph structure data, their design is geared towards undirected graphs. The attention weights are purely driven by data learning, which cannot express the directional characteristics of water flow from upstream to downstream, nor can it incorporate known hydrophysical priors such as topography and rainfall.

[0005] (4) The ground measurement stations are sparse, making it difficult to provide dense monitoring labels for the entire basin; while relying entirely on physical models cannot give full play to the advantages of remote sensing big data. Summary of the Invention

[0006] To address the problems identified in the prior art, this application provides a method for dynamic monitoring of soil erosion based on remote sensing technology.

[0007] The technical solution of this application embodiment is implemented as follows: Firstly, this application embodiment provides a method for dynamic monitoring of soil erosion based on remote sensing technology, comprising: For multi-source remote sensing data, digital elevation model data, meteorological data, soil property data, and ground-measured data collected in the monitoring basin for the year to be monitored, extract the dynamic features that change monthly and the static features that remain unchanged throughout the year; construct a dynamic feature tensor based on the dynamic features; and construct a static feature tensor based on the static features. Based on the digital elevation model, the D8 single-flow direction algorithm is used to construct a directed acyclic graph in the direction of flow. In this graph, multiple grids of the monitored watershed are used as multiple nodes, the direction of water flow is used as directed edges, and grids in the watershed whose slope to all neighboring directions is not greater than zero are marked as sinks. The dynamic feature tensor and the static feature tensor are input into the trained deep learning model for dynamic monitoring of soil erosion to obtain the annual soil erosion modulus for each grid in the entire watershed. The deep learning model for dynamic monitoring of soil erosion is based on a physical information graph neural network constructed with the directed acyclic graph as topological constraint, and is obtained through composite supervised training using a combination of sorting consistency loss, measured correction loss of actual measurement stations, and terrain adaptive spatial smoothing loss.

[0008] Preferably, the physical information graph neural network includes an input embedding layer, a physical injection graph attention layer, a temporal aggregation layer, and an erosion modulus output layer; The pre-construction process of the deep learning model for dynamic monitoring of soil erosion includes: The historical dynamic feature tensor and historical static feature tensor, as well as the corresponding historical annual soil erosion modulus for each grid in the whole watershed, are obtained and the corresponding directed acyclic graph is constructed. The physical weights of each directed edge, which are fused with static terrain weights and time-varying rainfall driving factors, are calculated and normalized to obtain a time-varying normalized adjacency matrix sequence. The historical dynamic feature tensor and historical static feature tensor are fused and time-encoded through the input embedding layer of the physical information graph neural network. Message passing and feature update processing are performed through the physical injection graph attention layer of the physical information graph neural network. During message passing, the logarithm of the physical weights in the time-varying normalized adjacency matrix sequence is used as a bias term and superimposed on the attention score calculation formula of the graph attention network to apply physical prior constraints to upstream neighbor edges. The temporal aggregation layer of the physical information graph neural network automatically weights the monthly data of 12 months under each node into annual data. The erosion modulus output layer of the physical information graph neural network outputs the historical annual soil erosion modulus for the entire watershed, grid by grid. Based on a composite loss function constructed using sorting consistency loss, measured correction loss, and terrain-adaptive spatial smoothing loss, the physical information graph neural network is trained using the historical dynamic feature tensor, the historical static feature tensor, and the historical annual soil erosion modulus to obtain the deep learning model for dynamic monitoring of soil and water loss.

[0009] Preferably, the construction of a directed acyclic graph in the direction of flow using the D8 single-flow algorithm based on the digital elevation model includes: dividing the monitored watershed into multiple grids as multiple nodes; Depression filling was performed on the digital elevation model of the monitored watershed to obtain a hydrologically consistent elevation surface; Based on the hydrologically consistent elevation surface, the D8 single-flow direction algorithm is used to determine the flow direction of each node: for each node, the slope from the node to its 8 neighboring nodes is calculated, and the direction with the largest and positive slope is determined as the flow direction. If the slope of a node to all neighborhood directions is not greater than zero, then mark it as a sink. Based on the sink, the direction of water flow, and the slope of each node to its 8 neighboring nodes, the cumulative flow of each node is recursively calculated in a topological order from upstream to downstream. A directed acyclic graph is constructed based on each node, the direction of water flow corresponding to the node as the directed edge of the node, the sink and the cumulative flow.

[0010] Preferably, the step of calculating the physical weight of each directed edge by fusing static terrain weight and time-varying rainfall driving factor, and obtaining a time-varying normalized adjacency matrix sequence after normalization, includes: For the slope from a node to a neighboring node in the directed acyclic graph, a static terrain weight calculation model is used to calculate the fused static terrain weight. For the nodes in the directed acyclic graph, the time-varying rainfall driving factor calculation model is used to calculate the time-varying rainfall driving factor; The product of the fused static terrain weight and the time-varying rainfall driving factor is determined as the physical weight of the directed edge from the node to the neighboring nodes. Normalize the in-degree of the target node for the physical weight of each directed edge, and output a time-varying normalized adjacency matrix sequence. The static terrain weight calculation model satisfies the following formula:

[0011] In the above formula, To integrate static terrain weights, For nodes To neighboring nodes The slope, The average gradient of the watershed; For nodes The cumulative amount of the confluence, This represents the maximum cumulative runoff in the basin. The calculation model for the time-varying rainfall driving factor satisfies the following formula:

[0012] In the above formula, As a driving factor for time-varying rainfall, For nodes In monthly time steps The erosive power of monthly rainfall, For nodes The maximum monthly rainfall erosion during the training period.

[0013] Preferably, the step of normalizing the in-degree of the target node for the physical weight of each directed edge and outputting a time-varying normalized adjacency matrix sequence includes: For the physical weight of each directed edge, the sum of the physical weights of all incoming edges of the downstream target node of the directed edge is used as the normalized denominator. The physical weight is divided by the denominator to obtain the normalized weight value of the directed edge. After performing the above normalization operation on all directed edges, we obtain the normalized weight values ​​of all directed edges. The normalized weights of all directed edges are combined to obtain a time-varying normalized adjacency matrix sequence.

[0014] Preferably, the step of fusing the historical dynamic feature tensor and the historical static feature tensor through the input embedding layer of the physical information graph neural network and adding time encoding includes: For each node in the input embedding layer of the physical information graph neural network at each time step, dynamic feature projection and static feature projection are performed to obtain the projected dynamic feature vector and the projected static feature vector. The projected dynamic feature vector is concatenated with the projected static feature vector, and then subjected to a nonlinear transformation through the fully connected layer of the physical information graph neural network. Time encoding is superimposed on the results of the nonlinear transformation.

[0015] Preferably, the message passing and feature update processing through the physical injection graph attention layer of the physical information graph neural network includes: For each downstream node of the physical injection graph attention layer, the union of the upstream neighbor set of the downstream node and the self-loop is defined as the message aggregation neighborhood. The message aggregation neighborhood is used to limit the range of upstream nodes that the downstream node considers when updating its own state. The logarithmic value of the physical weight in the time-varying normalized adjacency matrix sequence is used as a logarithmic bias term. This is superimposed on the upstream neighbor edges into the attention score calculation formula of the graph attention network of the physical information graph neural network. It is added to the data-driven scaling dot product attention term of the graph attention network to form the attention score of each upstream node in the message aggregation neighborhood corresponding to the downstream node, thus forming the basis for updating the downstream node. The attention scores of each upstream node in the message aggregation neighborhood of the updated downstream node are normalized using Softmax to obtain normalized attention scores; the normalized attention scores are then converted into probability distributions to obtain the attention ratio of each upstream node to the downstream node in the message aggregation neighborhood. Using the attention ratio as the weight, the value vectors of each node in the message aggregation neighborhood are weighted and summed to obtain the aggregated message received by the downstream node from each attention head; The aggregated messages from multiple attention heads are concatenated and multi-view information is fused using a linear projection matrix to obtain the attention output of the physical injection graph attention layer. In the attention layer of the physical injection graph, a Post-Norm residual structure is adopted. The attention output is randomly deactivated and then added to the input residual. After layer normalization, the intermediate features are obtained. In the attention layer of the physical injection graph, the intermediate features are nonlinearly transformed by a feedforward network. The feedforward network is a two-layer fully connected network. The first layer expands the feature dimension and then uses the GELU activation function. The second layer restores the feature dimension to the original dimension, thus obtaining the spatiotemporal node feature representation of the physical injection graph attention layer output by fusing upstream multi-order neighborhood information and physical constraints. The spatiotemporal node feature representation output by the physical injection graph attention layer is used as the input of the next physical injection graph attention layer. By stacking multiple physical injection graph attention layers, each downstream node aggregates the feature information of its upstream multi-order neighborhood.

[0016] Preferably, the temporal aggregation layer of the physical information graph neural network automatically weights the monthly data for each node into annual data for the 12 months; the erosion modulus output layer of the physical information graph neural network outputs the historical annual soil erosion modulus for the entire watershed, grid by grid, including: The time-series aggregation layer employs a learnable time attention mechanism, which sums the time weights of each node at each monthly time step to compress monthly features into annual comprehensive features. The annual comprehensive features of the erosion modulus output layer are mapped sequentially through two fully connected networks, and a smooth non-negative activation function is used to constrain the output to be a positive value for the historical annual soil erosion modulus of the entire watershed grid by grid.

[0017] Preferably, the composite loss function satisfies the following formula:

[0018] In the above formula, For compound loss, For sorting consistency loss, To measure and correct the loss, For terrain-adaptive spatial smoothing loss, Weights are lost for consistency in sorting. To measure and correct the loss weights, The terrain-adaptive spatial smoothing loss weights are used; The order consistency loss satisfies the following formula:

[0019] In the above formula, For consistent sampling of sample pairs, For each batch, two nodes are randomly selected from N grids across the entire watershed to form an ordered pair. For sorting interval, This is the model's first predicted value. This is the model's second predicted value; The measured correction loss satisfies the following formula:

[0020] In the above formula, M represents the total number of measured stations. For the first Grid index of each measured site, For the model to the first Predicted values ​​for each measured site For the first Measured erosion modulus of each actual measurement site; The terrain-adaptive spatial smoothing loss satisfies the following formula:

[0021] In the above formula, The set of edges between adjacent grids in a graph structure. Let the total number of edges be . Indicates two adjacent grid nodes. For nodes The slope, For neighboring nodes The slope, The attenuation coefficient is... For terrain similarity weights, For nodes The erosion modulus, For nodes The erosion modulus.

[0022] Secondly, embodiments of this application provide a dynamic monitoring system for soil erosion based on remote sensing technology, comprising: the system using the method described herein to dynamically monitor soil erosion in the monitored watershed.

[0023] The beneficial effects of the technical solutions provided in this application include at least the following: 1. To address the problem that "empirical models ignore nonlinear interactions of factors and hydrological connectivity effects": This application constructs a deep learning model for dynamic monitoring of soil erosion based on a physical information graph neural network. It uses a directed acyclic graph constructed with the direction of water flow as a topological constraint, enabling the model to transmit information along the upstream and downstream hydrological paths. This breaks through the limitation of independent multiplication of factors in the traditional RUSLE model and achieves effective modeling of nonlinear interactions and watershed hydrological connectivity effects.

[0024] 2. To address the issue that "traditional machine learning and CNN cannot model the spatial correlation of water flow direction": This application divides the monitored watershed into multiple grids as nodes, uses the D8 single-flow direction algorithm to determine the water flow direction as directed edges, and constructs a downstream directed acyclic graph, enabling the model to pass messages from upstream to downstream water flow direction, rather than relying on regular rectangular neighborhoods or spatial distances, thereby accurately depicting the watershed's unique directional spatial correlation.

[0025] 3. To address the issue that "standard graph neural networks cannot express directionality and incorporate hydrophysical priors": In the attention mechanism of the physical information graph neural network, this application uses the logarithmic value of the physical weights calculated from the static topographic weights and the time-varying rainfall driving factor as a logarithmic bias term, which is superimposed on the data-driven attention score. This allows the model to be simultaneously controlled by the collaborative control of data-driven learning and physical prior constraints, thus expressing the directional characteristics of water flow and incorporating hydrophysical prior knowledge such as topography and rainfall.

[0026] 4. To address the problem of "sparse ground-based measurement stations and insufficient supervision signals": This application constructs a composite loss function that includes ordination consistency loss, measured correction loss, and terrain-adaptive spatial smoothing loss. It utilizes the RUSLE model covering the entire watershed to provide relative ordination relationships as weak supervision signals, and combines the absolute numerical calibration provided by sparse measurement stations as strong supervision signals. Even under the condition of sparse measured data, a high-precision grid-by-grid soil erosion modulus prediction model for the entire watershed can still be trained.

[0027] It should be understood that the above general description and the following detailed description are exemplary and explanatory only, and are not intended to limit the technical solutions provided in the embodiments of the present invention. Attached Figure Description

[0028] To more clearly illustrate the technical solutions in the embodiments of this application, the accompanying drawings used in the description of the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of this application. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort, wherein: Figure 1 A flowchart illustrating a method for dynamic monitoring of soil erosion based on remote sensing technology, provided in an embodiment of this application; Figure 2 This is a schematic diagram of the network structure of a deep learning model for dynamic monitoring of soil erosion based on remote sensing technology, provided in an embodiment of this application. Detailed Implementation

[0029] To make the objectives, technical solutions, and advantages of the embodiments of this application clearer, the technical solutions of the embodiments of this application will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of this application, not all embodiments. The following embodiments are used to illustrate this application, but are not intended to limit the scope of this application. Based on the embodiments in this application, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of this application.

[0030] In the following description, references are made to “some embodiments,” which describe a subset of all possible embodiments. However, it is understood that “some embodiments” may be the same subset or different subsets of all possible embodiments and may be combined with each other without conflict.

[0031] It should be noted that the terms "first, second, and third" used in the embodiments of this application are merely to distinguish similar objects and do not represent a specific ordering of objects. It is understood that "first, second, and third" can be interchanged in a specific order or sequence where permitted, so that the embodiments of this application described herein can be implemented in an order other than that illustrated or described herein.

[0032] It will be understood by those skilled in the art that, unless otherwise defined, all terms used herein (including technical and scientific terms) have the same meaning as commonly understood by one of ordinary skill in the art to which the embodiments of this application pertain. It should also be understood that terms such as those defined in general dictionaries should be understood to have a meaning consistent with their meaning in the context of the prior art, and should not be interpreted in an idealized or overly formal sense unless specifically defined as herein.

[0033] like Figure 1 As shown in the figure, this application provides a method for dynamic monitoring of soil erosion based on remote sensing technology, which may include the following steps: Step 1: Extract the monthly dynamic features and the year-invariant static features from the multi-source remote sensing data, digital elevation model data, meteorological data, soil property data, and ground-measured data collected within the monitoring basin for the year to be monitored; construct a dynamic feature tensor based on the dynamic features; construct a static feature tensor based on the static features. Step 2: Based on the digital elevation model, the D8 (eight-flow direction) single-flow direction algorithm is used to construct a directed acyclic graph in the direction of flow. In this algorithm, multiple grids of the monitoring watershed are used as multiple nodes, the direction of water flow is used as directed edges, and grids in the watershed whose slope to all neighboring directions is not greater than zero are marked as sinks. Step 3: Input the dynamic feature tensor and the static feature tensor into the trained deep learning model for dynamic monitoring of soil erosion to obtain the annual soil erosion modulus for each grid across the entire watershed. The deep learning model for dynamic monitoring of soil erosion is based on a physical information graph neural network constructed using the directed acyclic graph (DAG) as topological constraints, and is trained under composite supervision using a combination of ranking consistency loss, measured station correction loss, and terrain-adaptive spatial smoothing loss. In the attention mechanism of the physical information graph neural network, the logarithm of the time-varying physical weights of the DAG is added as a bias term to the attention score, coordinating with the data-driven attention term to control the information transmission intensity.

[0034] Furthermore, step 1 specifically includes the following steps: Step 1.1: Obtain five types of data from the monitored watershed: The first category consists of multi-source remote sensing data, including: Landsat-8 / 9 or Sentinel-2 multispectral imagery, used to acquire surface reflectance of the watershed and calculate vegetation indices and cover; Sentinel-1 Synthetic Aperture Radar (SAR) imagery, used to acquire VV (Vertical-Vertical polarization) and VH (Vertical-Horizontal polarization) backscattering coefficients of the watershed; 16-day vegetation index products from the Moderate Resolution Imaging Spectroradiometer (MODIS), used to assist in time series interpolation; MODIS surface temperature products; and annual land use classification data.

[0035] The second category is digital elevation model data, including data collected from watersheds by the Shuttle Radar Topography Mission 30-meter resolution (SRTM30m) or Advanced Land Observing Satellite Phased Array type L-band Synthetic Aperture Radar 12.5-meter resolution (ALOS PALSAR 12.5m) products.

[0036] The third category is meteorological data, including daily precipitation and average daily temperature provided by the ground meteorological station network.

[0037] The fourth category is soil property data, including the percentage of sand, silt, and clay by mass and the organic matter content in the watershed, which are derived from the watershed regional soil survey database.

[0038] The fifth category is ground-based measured data, including annual measured sediment transport of the watershed and annual measured soil erosion modulus of runoff plots using hydrological stations. This type of data can be used for numerical calibration in model training.

[0039] Step 1.2 Data preprocessing and feature construction includes the following steps: S1.2.1 Image and meteorological data preprocessing: Multispectral imagery underwent radiometric calibration, atmospheric correction, cloud and cloud shadow mask removal based on quality-marked bands, and topographic correction of Digital Elevation Model (DEM) data. Synthetic Aperture Radar (SAR) imagery underwent precise orbit file application, thermal noise removal, radiometric calibration to backscattering coefficient conversion, 7×7 window Lee Sigma filtering for speckle reduction, and range-Doppler topographic correction followed by conversion to decibel values. SAR imagery from multiple periods was synthesized with a monthly time step. The Normalized Difference Vegetation Index (NDVI) product was calculated using the maximum value from each period to reduce the residual impact of cloud pollution, while the median value was used for other bands. Surface meteorological station observations were interpolated to a grid scale using co-kriging with DEM elevation as a covariate. All data were uniformly resampled to a 30m spatial resolution and registered to the same Universal Transverse Mercator (UTM). In the Mercator projection coordinate system.

[0040] S1.2.2 Spatial Discretization of the Watershed Based on the spatial resolution of the digital elevation model (DEM), the monitored watershed is divided into regular square grid cells, with 30m×30m as the default resolution, consistent with the resolution of Landsat-8 / 9 or Sentinel-2 multispectral imagery. Watershed analysis is performed based on the DEM, extracting the catchment area boundaries, removing grid cells outside the boundaries, and excluding data-free areas, permanent water bodies, and construction land that do not participate in the soil erosion process. The remaining grid cells form a set of nodes. Each grid cell serves as a node in the subsequent graph network and has unique spatial coordinates.

[0041] S1.2.3 Feature Extraction Based on the preprocessed data, two types of variables, dynamic features and static features, are extracted for each grid cell; The dynamic features are variables that change monthly, with a total of 12 monthly time steps (T=12 months). Specifically, they include: Normalized Difference Vegetation Index (NDVI), Enhanced Vegetation Index (EVI), Bare Soil Index (BSI), and Fractional Vegetation Cover (FVC) obtained from multispectral image band calculations; VV and VH polarization backscattering coefficients obtained from synthetic aperture radar data; land surface temperature obtained from medium resolution imaging spectrometer products; and monthly cumulative precipitation, monthly rainfall erosivity, and monthly average temperature calculated from meteorological interpolation data. The dynamic feature dimension is... The monthly rainfall erosivity was calculated daily using a daily rainfall erosivity model and then accumulated over the month. Static characteristics remain unchanged throughout the year, specifically including: elevation, slope, sine and cosine components of aspect, plane curvature, profile curvature, and topographic wetness index (TWI) calculated by the digital elevation model; soil erodibility K-factor calculated from soil data using the Erosion-Productivity Impact Calculator (EPIC) formula; and unique thermal codes for land use types, including the number of categories. =6, specifically including arable land, forest land, grassland, water areas, construction land, and unused land. Therefore, the static feature dimension is... 14.

[0042] S1.2.4 Standardization and Feature Tensor Construction All numerical features in both dynamic and static features are z-score standardized. The mean and standard deviation used for standardization are obtained statistically from the training set, and the same set of parameters is used in the validation, testing, and inference phases. Dynamic features are arranged into a three-dimensional tensor. Arrange static features into a two-dimensional tensor Where N is the total number of effective grid nodes in the watershed, and T=12 is the total number of monthly time steps. For dynamic feature dimension, is the static feature dimension.

[0043] Furthermore, step 2, which involves constructing a downstream directed acyclic graph (DAG) using the D8 single-flow algorithm based on the digital elevation model, may include the following steps: Step 2.1: Divide the monitoring basin into multiple grids, which will serve as multiple nodes; Step 2.2: Perform depression filling on the digital elevation model of the monitored watershed to obtain a hydrologically consistent elevation surface; Step 2.3: Based on the hydrologically consistent elevation surface, the D8 single-flow direction algorithm is used to determine the flow direction of each node: For each node, the slope from the node to its 8 neighboring nodes is calculated, and the direction with the largest and positive slope is determined as the flow direction; Step 2.4: If the slope of a node to all neighborhood directions is not greater than zero, then mark it as a sink point; Step 2.5: Based on the sink, the direction of water flow, and the slope of each node to its 8 neighboring nodes, recursively calculate the cumulative flow of each node in the topological order from upstream to downstream. Step 2.6: Construct a directed acyclic graph based on each node, the direction of water flow corresponding to the node as the directed edge of the node, the sink and the cumulative flow.

[0044] In one specific implementation of this embodiment, steps 2.1-2.6 are as follows: First, depression filling is performed on the digital elevation model. The Planchon-Darboux algorithm is used to eliminate pseudo-depressions caused by data noise, resulting in a hydrologically consistent elevation surface. Then, the D8 single-flow direction algorithm is used to determine the flow direction at each node. Calculate its directions to its eight neighboring nodes respectively. The slope, then the node To neighboring nodes The slope is: ,in Distance between node centers For nodes The slope, For neighboring nodes The slope is calculated, and the direction with the largest and positive slope among all nodes is determined as the flow direction; if the slope of a node is not greater than zero in all neighboring directions, it is marked as a sink point; then, based on the D8 flow direction results, the flow direction of each node is recursively calculated in topological order from upstream to downstream. Convergence accumulation The cumulative value of the confluence is equal to the sum of the cumulative values ​​of the node itself (denoted as 1) and all the cumulative values ​​of the upstream nodes flowing into it. Finally, each valid node in the node set is considered as a node; in the edge set, if the flow direction of a node points to another node, a directed edge is established. The topological characteristics of this downstream directed acyclic graph are as follows: each non-sink node has exactly one outgoing edge pointing to the neighboring node with the largest slope among its eight neighborhoods; the number of incoming edges for each node equals the number of upstream neighbors that choose that node as a downstream target, with values ​​ranging from 0 to 8; source nodes (at the watershed boundary or local high points) have no incoming edges, and sink nodes (watershed outlets) have no outgoing edges; since the water flows from high to low according to gravity, a downstream directed acyclic graph is formed. ,in, Represents a set of nodes. Represents the set of edges.

[0045] Furthermore, the physical information graph neural network described in step 3 above includes an input embedding layer, a physical injection graph attention layer, a temporal aggregation layer, and an erosion modulus output layer; The pre-construction process of the deep learning model for dynamic monitoring of soil erosion includes: Step 3.1a: Obtain the historical dynamic feature tensor and the historical static feature tensor, as well as the corresponding historical annual soil erosion modulus for each grid in the whole watershed and construct the corresponding directed acyclic graph. Calculate the fusion of static terrain weight and time-varying rainfall driving factor for each directed edge, and then normalize it to obtain the time-varying normalized adjacency matrix sequence. Step 3.2a: The historical dynamic feature tensor and historical static feature tensor are fused through the input embedding layer of the physical information graph neural network and time encoding is added; message passing and feature update processing are performed through the physical injection graph attention layer of the physical information graph neural network, wherein during message passing, the logarithm of the physical weights in the time-varying normalized adjacency matrix sequence is used as a bias term and superimposed on the attention score calculation formula of the graph attention network to apply physical prior constraints to the upstream neighbor edges; the temporal aggregation layer of the physical information graph neural network automatically weights the monthly data of 12 months under each node into annual data; the erosion modulus output layer of the physical information graph neural network outputs the historical annual soil erosion modulus of the entire watershed grid by grid. Step 3.3a: Based on the composite loss function constructed using sorting consistency loss, measured correction loss and terrain adaptive spatial smoothing loss, the physical information graph neural network is trained using the historical dynamic feature tensor, the historical static feature tensor and the historical annual soil erosion modulus to obtain the deep learning model for dynamic monitoring of soil and water loss.

[0046] Here, in the attention mechanism of the physical information graph neural network, the logarithm of the physical weights in the time-varying normalized adjacency matrix sequence is superimposed as a bias term to the attention score, and together with the data-driven scaling dot product attention term, the information transmission intensity is controlled to obtain the annual soil erosion modulus of the whole watershed grid by grid.

[0047] Further, step 3.1a, which involves calculating the physical weights of each directed edge by fusing the static terrain weights and the time-varying rainfall driving factors, and then normalizing the weights to obtain a time-varying normalized adjacency matrix sequence, may include the following steps: Step 3.1.1a: For the slope from a node to a neighboring node in the directed acyclic graph, use the static terrain weight calculation model to calculate the fused static terrain weight; Step 3.1.1b: For the nodes in the directed acyclic graph, use the time-varying rainfall driving factor calculation model to calculate the time-varying rainfall driving factor; Step 3.1.1c: The product of the fused static terrain weight and the time-varying rainfall driving factor is determined as the physical weight of the directed edge from the node to the neighboring nodes; Step 3.1.1d: Normalize the in-degree of the target node for the physical weight of each directed edge, and output the time-varying normalized adjacency matrix sequence; The static terrain weight calculation model described in step 3.1.1a satisfies the following formula:

[0048] In the above formula, To integrate static terrain weights, For nodes To neighboring nodes The slope, The average gradient of the watershed; For nodes The cumulative amount of the confluence, The maximum runoff accumulation in the basin is represented by the static topographic weight calculation model. In the design, the exponents of the slope ratio (0.5) and runoff accumulation ratio (0.4) are determined by comparative experiments within the range of [0.3, 0.7] with a step size of 0.1. This can control the degree of nonlinear influence of slope and runoff area on water and sediment transport capacity, respectively.

[0049] The calculation model for time-varying rainfall driving factors described in step 3.1.1b satisfies the following formula:

[0050] In the above formula, As a driving factor for time-varying rainfall, For nodes In monthly time steps The erosive power of monthly rainfall, For nodes The maximum monthly rainfall erosivity during the training period. In the design of the time-varying rainfall driving factor calculation model, the exponent was set to 0.6, which allows the physical weights to exhibit non-linear saturation characteristics in response to rainfall.

[0051] In step 3.1.1c, the aim is to assign a physical weight reflecting the hydrological transport capacity to each directed edge. physical weight By fusing static terrain weights With time-varying rainfall driving factors The product is determined as follows: .

[0052] In step 3.1.1d, the step of normalizing the in-degree of the target node for the physical weight of each directed edge and outputting a time-varying normalized adjacency matrix sequence includes: Step 3.1.1.1: For the physical weight of each directed edge, use the sum of the physical weights of all incoming edges of the downstream target node of the directed edge as the normalized denominator, divide the physical weight by the denominator to obtain the normalized weight value of the directed edge. Step 3.1.1.2: After performing the above normalization operation on all directed edges, the normalized weight values ​​of all directed edges are obtained; Step 3.1.1.3: Combine the normalized weights of all directed edges to obtain a time-varying normalized adjacency matrix sequence.

[0053] In one specific implementation of this embodiment, in steps 3.1.1.1-3.1.1.3, the obtained physical weights are normalized by the in-degree of the target node to obtain time-varying normalized weights. : ,in, This is the sum of the original weights of all incoming edges. To prevent division by zero, the numerical stability constant is normalized. This normalization ensures that the sum of the weights of all incoming edges to each node is 1, making the contributions of each upstream source comparable during information aggregation. The final output is a time-varying normalized adjacency matrix sequence. ,in Let N be the sparse matrix of effective grid nodes within the watershed, and let N be the total number of effective grid nodes within the watershed.

[0054] Further, step 3.2a, which involves fusing the historical dynamic feature tensor and the historical static feature tensor through the input embedding layer of the physical information graph neural network and adding time encoding, includes: Step 3.2.1a: For each node in the input embedding layer of the physical information graph neural network at each time step, perform dynamic feature projection and static feature projection to obtain the projected dynamic feature vector and the projected static feature vector; Step 3.2.2a: Concatenate the projected dynamic feature vector with the projected static feature vector, and then perform a nonlinear transformation through the fully connected layer of the physical information graph neural network; Step 3.2.3a: Superimpose time encoding on the results after nonlinear transformation.

[0055] Steps 3.2.1a-3.2.3a mainly implement the separation and linear projection of historical dynamic features and historical static features by the input embedding layer. After splicing, fusion and nonlinear transformation, time coding is superimposed and layer normalization is performed to generate an initial node representation that integrates heterogeneous information and encodes temporal position, providing a spatiotemporally unified feature input for subsequent physical injection graph convolution.

[0056] Example for each node At each monthly time step, perform the following operations: (1) Dynamic feature projection in, For nodes In monthly time steps The dynamic characteristics, These are learnable, dynamic feature vectors.

[0057] (2) Static feature projection ,in, For nodes In monthly time steps static characteristics, These are learnable static feature vectors.

[0058] (3) Concatenation: The projection results are concatenated into a d-dimensional vector, which is then subjected to a nonlinear transformation through a fully connected layer, time-encoded, and the initial node representation is obtained through layer normalization. : ,in This represents the concatenated vector. This represents the feature fusion weight matrix. This is for time encoding. d represents the embedding dimension. In comparative experiments with d∈{32,64,128,256}, this value provides the optimal balance between prediction accuracy and computational efficiency: when d=32, insufficient model capacity leads to a decrease in NSE (Nash-Sutcliffe Efficiency coefficient) of approximately 0.05; when d=128 and d=256, the accuracy improvement is less than 0.01, but the training time increases by approximately 1.8 times and 3.5 times, respectively. In practical applications, d can be selected within the range of 32 to 256 based on the total feature dimension and computational resources. In this embodiment... Time coding Monthly time step number Calculated using sine and cosine functions: ,in Indexed by feature dimensions, This represents the 2k-th even-numbered dimension in the position encoding vector of the t-th month's time step. This represents the encoded value of the (2k+1)th odd-numbered dimension in the position encoding vector of the t-th month's time step.

[0059] Further, the message passing and feature update processing via the physical injection graph attention layer of the physical information graph neural network described in step 3.2a includes the following steps: Step 3.2.1b: For each downstream node of the physical injection graph attention layer, the union of the upstream neighbor set and the self-loop of the downstream node is defined as the message aggregation neighborhood. The message aggregation neighborhood is used to limit the range of upstream nodes considered by the downstream node when updating its own state. Step 3.2.2b: The logarithmic value of the physical weight in the time-varying normalized adjacency matrix sequence is used as a logarithmic bias term. This is then added to the attention score calculation formula of the graph attention network of the physical information graph neural network, along with the data-driven scaling dot product attention term of the graph attention network. Together, these terms constitute the attention scores of each upstream node in the message aggregation neighborhood corresponding to the downstream node, and the downstream node is updated. Step 3.2.3b: Perform Softmax normalization on the attention scores of each upstream node in the message aggregation neighborhood of the updated downstream node to obtain a normalized attention score; convert the normalized attention score into a probability distribution to obtain the attention ratio of each upstream node to the downstream node in the message aggregation neighborhood. Step 3.2.4b: Using the attention ratio as the weight, perform a weighted summation of the value vectors of each node in the message aggregation neighborhood to obtain the aggregated message received by the downstream node from each attention head; Step 3.2.5b: Concatenate the aggregated messages from multiple attention heads, fuse the multi-view information using a linear projection matrix, and obtain the attention output of the physical injection graph attention layer; Step 3.2.6b: In the attention layer of the physical injection graph, a Post-Norm residual structure is adopted. After random deactivation of the attention output, it is added to the input residual and then normalized by the layer to obtain the intermediate features. Step 3.2.7b: In the attention layer of the physical injection graph, the intermediate features are nonlinearly transformed by a feedforward network. The feedforward network is a two-layer fully connected network. The first layer expands the feature dimension and uses the GELU (Gaussian Error Linear Unit) activation function. The second layer restores the feature dimension to the original dimension, thus obtaining the spatiotemporal node feature representation of the physical injection graph attention layer output by fusing upstream multi-order neighborhood information and physical constraints. Step 3.2.8b: Use the spatiotemporal node feature representation output by the physical injection graph attention layer as the input of the next physical injection graph attention layer. By stacking multiple physical injection graph attention layers, each downstream node aggregates the feature information of its upstream multi-order neighborhood.

[0060] Steps 3.2.1b-3.2.8b, specifically the improvement to the physical injection graph attention layer, are the core of this invention and the key difference from standard graph attention networks. (Co-stacking) Layered Physical Injection Graph Attention Layer The value range is 1 to 4. Layered graph convolutions are pre-selected based on the monitoring basin's resolution and scale. Each layer is processed independently at each monthly time step, with the time-varying normalized adjacency matrix corresponding to that monthly time step providing physical weight constraints. This ensures that the node feature aggregation process balances data-driven feature similarity learning with constraints from the physical laws of hydrological connectivity. Each layer employs a multi-head attention mechanism to capture multi-subspace feature associations in parallel. After residual connections and layer normalization, a feedforward network performs nonlinear feature transformation, outputting a spatiotemporal node feature representation that integrates upstream multi-order neighborhood information and physical constraints.

[0061] In one specific embodiment of this example, in steps 3.2.1b-3.2.8b, the first... layer( Taking the t-th month's time step as an example, the input is the node features output from the previous layer. The processing flow is as follows: Neighborhood definition For each node (The target node to be updated) has a message aggregation neighborhood that is the union of its upstream neighbor set and its self-loop: ,in For nodes In a hydrological topology graph, the upstream neighbor set is the set of all nodes to which water flows. .node It is contained in the neighborhood. This allows the attention mechanism to retain the original characteristics of nodes while aggregating information from upstream neighbors. Adding self-loops ensures that nodes retain their own characteristics when aggregating upstream information; this operation is performed only during message passing computation and does not alter the constructed directed acyclic graph.

[0062] Multi-head physical injection attention mechanism: Employing a multi-head attention mechanism, the number of heads In this embodiment, a value of 4 is chosen (this value performs best in the comparative experiment with H∈{1,2,4,8}: H=1 (single head) cannot distinguish multi-scale neighborhood associations, resulting in a decrease in NSE of approximately 0.04; H=2 has slightly lower accuracy than H=4; when H=8, the dimension per head decreases to 8 dimensions, limiting expressive power and providing no further improvement in accuracy). The dimension per head... For the first Size ( A downstream node relates to its upstream neighbor x. The attention score is composed of the sum of the data-driven term and the physical bias term:

[0063] in, The feature vector of the downstream node. The feature vector of the upstream node; Let h be the query projection matrix and key projection matrix of the h-th head; the first term in the formula is the scaled dot product attention, which is calculated from the node features through learnable parameters; For calculating the edges Normalized physical weights at monthly time step t, indicator function exist The value is 1 if the condition is met, and 0 otherwise, ensuring that physical bias is applied only to upstream neighbor edges, and no bias is applied to the node itself. It is a very small positive number (used for numerical stability to prevent the logarithm from being undefined).

[0064] Then, for the neighborhood of the downstream node The attention scores of each upstream node are normalized using Softmax, and gradient clipping is used to control the magnitude of the backpropagation gradient, thus converting the attention scores into a probability distribution. , For the first The proportion of attention paid by upstream nodes to downstream nodes in the head, satisfying .

[0065] Then, the messages are aggregated, and the value vectors of each node in the neighborhood are weighted and summed to obtain the value of each node from the first node. The aggregated message received by the head Physical meaning: The new state of a node is a weighted average of its neighborhood information according to its "importance" (attention weight), where, For the first The projection matrix of the head value, =1...H, where H is the total number of heads.

[0066] Multi-head splicing and output projection Each attention head captures node association patterns from different subspaces, The aggregated message vector received by the header is concatenated into a complete representation as follows: via linear projection matrix By fusing information from multiple perspectives, the attention output of this layer is obtained: It is automatically optimized through backpropagation during model training.

[0067] Residual connectivity and feedforward network The feature update is performed in two steps using a Post-Norm residual structure (meaning that normalization is placed after the residual connection): Step 1: Injecting Attention Information After randomly deactivating the attention output, it is added to the input residual, and then normalized to obtain intermediate features. : , This indicates the output of multi-head splicing and output projection. This indicates that the attention output is randomly deactivated to suppress overfitting; residual connections allow nodes to superimpose new information from neighborhood aggregation while retaining their original features.

[0068] Step 2: Nonlinear Feature Transformation The feed-forward network (FFN) of the attention layer of the physical injection graph undergoes deep feature processing to satisfy the following formula: The feedforward network FFN is a two-layer fully connected network. The first layer expands the features from d-dimensional to 4d-dimensional and uses the GELU activation function. The second layer reduces the features from 4d-dimensional back to d-dimensional.

[0069] Further, in step 3.2a, the automatic weighting of monthly data for the 12 months under each node into annual data through the temporal aggregation layer of the physical information graph neural network; and the output of the historical annual soil erosion modulus for the entire watershed grid-by-grid through the erosion modulus output layer of the physical information graph neural network, including: Step 3.2.1c: A learnable temporal attention mechanism is adopted for the temporal aggregation layer to perform weighted summation of the temporal weights of each node at each monthly time step, compressing the monthly features into annual comprehensive features; Step 3.2.2c: The annual comprehensive features of the erosion modulus output layer are mapped sequentially through two fully connected networks, and a smooth non-negative activation function is used to constrain the output to be a positive value for the historical annual soil erosion modulus of the entire watershed grid by grid.

[0070] In step 3.2.1c above, a learnable temporal attention mechanism is used for the temporal aggregation layer, with the goal of aggregating the feature vectors of each node at 12 monthly time steps into a single annual comprehensive feature: Let the total number of graph convolution layers be L. After L layers of physical information graph convolution processing, each node has a monthly feature vector of T=12. It satisfies the following formula:

[0071] in, Let t be the total number of convolutional layers in the graph, t be the monthly time step, and T=12 be the total number of monthly time steps.

[0072] The temporal aggregation layer employs a learnable temporal attention mechanism to complete the aggregation, as follows: Calculation of time weights Node In the The features of each monthly time step are projected, activated, and then dot-producted with the attention vector to obtain an attention score. This score is then normalized using Softmax to obtain the time weights. It satisfies the following formula:

[0073] in, The time projection matrix is ​​randomly initialized using Xavier and its element values ​​are updated during training via backpropagation. Linear projection of dimensional node features to 3D attention hidden space, For learnable attention vectors, For attention hidden layer dimensions.

[0074] The time weights of each node at each monthly time step are weighted and summed. By time weight coefficient, for The annual comprehensive characteristics are obtained by weighting and summing the total number of monthly time steps according to their time weights. ,in, Indicates the process After layer graph convolution processing, nodes In the The eigenvectors of the month.

[0075] The erosion modulus output layer mentioned in step 3.2.2c above is specifically implemented as follows: The erosion modulus output layer maps the annual comprehensive features output by the time-series aggregation layer to predicted soil erosion modulus values. This layer extracts high-order nonlinear correlations through a layer-by-layer dimensionality-reduced fully connected network and constrains the output to be strictly positive by a smooth non-negative activation function, ultimately outputting positive historical annual soil erosion modulus values ​​for the entire watershed gridwise.

[0076] The annual comprehensive features are mapped to the erosion modulus via a two-layer fully connected network. :

[0077] in, , , , All are learnable parameters. Softplus function. Ensure output It conforms to the physical constraints of soil erosion modulus, and the output unit is t / (km²). 2 ·a), It is a natural constant. For nodes The annual comprehensive feature vector output by a multi-scale graph neural network.

[0078] Step 3.3a describes the use of a composite loss function constructed with ranking consistency loss, measured correction loss, and terrain-adaptive spatial smoothing loss. This is followed by training the physical information graph neural network using the historical dynamic feature tensor, the historical static feature tensor, and the historical annual soil erosion modulus to obtain the deep learning model for dynamic monitoring of soil erosion. Specifically, this model includes: S4.1 This step constructs two complementary types of supervision signals to support the training of the physical information graph neural network under sparse data conditions: Weakly supervised signal: derived from the Revised Universal Soil Loss Equation (RUSLE) model calculations, covering N grids across the entire watershed. Its core purpose is to provide the model with a relative ranking of erosion intensity. Due to the systematic bias in RUSLE, this signal only retains the spatial distribution trend and is not used as an absolute numerical label. The strong monitoring signal originates from M measured stations (including runoff plots and hydrological stations), covering only... The historical dynamic feature tensor, historical static feature tensor, and historical annual soil erosion modulus of the total number of effective grid nodes N within the watershed are used as training data. Their core purpose is to provide absolute numerical calibration anchor points for the physical information graph neural network. This signal is numerically reliable but has limited spatial coverage, complementing the weakly supervised signal and jointly supporting model training under sparse data conditions.

[0079] S4.2 Construct the composite loss function, satisfying the following formula:

[0080] In the above formula, For compound loss, For sorting consistency loss, To measure and correct the loss, For terrain-adaptive spatial smoothing loss, Weights are lost for consistency in sorting. To measure and correct the loss weights, The terrain-adaptive spatial smoothing loss weights are used; (The loss weights are determined through a 3D grid search, with the search range being...) ∈{0.1,0.3,0.5,1.0}、 ∈{1.0,5.0,10.0,20.0}、 There are 64 combinations in total, ∈{0.01,0.05,0.1,0.5}, and the final result is... The combination achieves the highest NSE on the validation set and the lowest relative error at the actual test sites. During training, the ratio of the norm of each loss gradient is monitored; if the ratio is greater than 20, the weights of the measured corrective loss are temporarily reduced. Upgraded to 5.0 to prevent overfitting of measured data.

[0081] The order consistency loss satisfies the following formula:

[0082] In the above formula, For consistent sampling of sample pairs, For each batch, two nodes are randomly selected from N grids across the entire watershed to form an ordered pair. For sorting interval, This is the model's first predicted value. This is the model's second prediction; each batch randomly selects two nodes from N grids across the entire watershed to form an ordered pair. Pairs with significant differences in RUSLE erosion modulus are preferentially retained (i.e., those satisfying the condition of...). ,in Take 0.5 times the Rusle standard deviation (τ is experimentally verified to eliminate about 38% of fuzzy pairings). and They are nodes and nodes RUSLE erosion modulus), total sampling right( The value of 500 was determined through comparative experiments on {100, 200, 500, 1000, 2000}. Through The comparative experiment on {0.01, 0.05, 0.1, 0.2, 0.5} determined the value to be 0.1. When the value is 0.01, the prediction results have only 78% consistency with RUSLE's spatial ranking. When the value is 0.5, the measured MAE of the site increases by approximately 15%. and The model's predicted values ​​are used. This loss only constrains the relative magnitudes of the model's predicted values ​​to maintain the spatial ranking trend of RUSLE, without constraining the absolute values, thus effectively avoiding the systematic bias inherited from RUSLE.

[0083] The measured correction loss satisfies the following formula:

[0084] In the above formula, M represents the total number of measured stations. For the first Grid index of each measured site, For the model to the first Predicted values ​​for each measured site For the first The measured erosion modulus of each actual site is used; the +1 in the logarithmic transformation is used to avoid taking the logarithm of zero values. The absolute values ​​after the logarithmic transformation are rigorously fitted on sparse measured sites to provide numerical calibration anchors for the model.

[0085] The terrain-adaptive spatial smoothing loss satisfies the following formula:

[0086] In the above formula, The set of edges between adjacent grids in a graph structure. Let the total number of edges be . Indicates two adjacent grid nodes. For nodes The slope, For neighboring nodes The slope, The attenuation coefficient (this value is within) (Determined in a comparative experiment ∈{0.5, 1.0, 2.0, 5.0}), controlling the decay rate of slope difference on the smoothing constraint. For terrain similarity weights, For nodes The erosion modulus, For nodes The erosion modulus. The smoothing constraint between adjacent grids with large slope differences is automatically weakened, allowing reasonable discontinuities at abrupt topographic changes; while the spatial continuity of the prediction results is maintained in areas with similar slopes.

[0087] S4.3 Training Strategy (1) Dataset partitioning The data used for training can be divided into training set, validation set and test set by year (e.g., training in the first 3 years, validation in the 4th year, and testing in the 5th year) to ensure that the data for all 12 months of the same year belong to the same set, avoid time leakage, and ensure the time independence of the model's generalization ability.

[0088] (2) Subgraph sampling training For large watersheds (>50,000 nodes), a DEM-based sub-watershed partitioning strategy is adopted, with each sub-watershed containing 1,000-5,000 core nodes. Four to eight subgraphs are randomly sampled per batch (a subgraph is a subset of nodes corresponding to the sub-watershed and its neighborhood). Each subgraph is expanded with two layers of neighborhood buffers to provide complete neighbor information for edge nodes, ensuring correct message aggregation. The composite loss is calculated on the core nodes of the subgraphs. By using sub-watershed samples as partial training data, full watershed coverage training is achieved, improving GPU memory efficiency.

[0089] (3) Optimizer and learning rate scheduling An adaptive moment estimation optimizer is selected, with an initial learning rate of [missing information]. The cosine annealing strategy is used to decay the temperature to... Weight decay coefficient To balance convergence speed and generalization performance.

[0090] (4) Regularization strategy The dropout rate was set to 0.1 (determined through comparative experiments), and the maximum norm of gradient clipping was set to 1.0 to prevent gradient explosion.

[0091] (5) Early termination and model selection Using the validation set NSE (Nash efficiency coefficient) as the monitoring metric, training is terminated when the NSE does not improve for 20 consecutive rounds, and the model parameters with the best performance on the validation set are saved.

[0092] Step S5: Model Reasoning and Result Generation (1) Data preparation and feature construction For the year to be monitored, remote sensing images and meteorological data are acquired according to S1.1-S1.2, and the S1.2 preprocessing workflow is executed to construct dynamic feature tensors. ; Synchronously update static feature tensor .

[0093] (2) Time-varying physical weight update Based on monthly rainfall erosivity in the year to be monitored Update the time-varying rainfall driving factors according to S2.2. This process then updates the time-varying normalized physical weight matrix. This step is particularly crucial because the annual variation in rainfall erosivity directly drives the dynamic adjustment of the spatial relationships of the graph structure.

[0094] (3) Model forward propagation Load the optimal model parameters saved during the training phase and perform forward propagation: input embedding layer → physical injection graph attention layer (L layer) → temporal aggregation layer → erosion modulus output layer, outputting grid-by-grid annual soil erosion modulus prediction values. The unit is .

[0095] (4) Classification of erosion intensity According to the "Classification and Grading Standard for Soil Erosion" (SL 190-2007), the erosion modulus is divided into six levels:

[0096] (5) Output of results Generate thematic maps of soil erosion intensity in GeoTIFF (geotagged image file) format, with the projection coordinate system consistent with the input DEM; statistically analyze the area proportion and spatial distribution characteristics of each erosion level; output a complete monitoring report including data description, method overview, accuracy verification, and results analysis, providing data support for dynamic assessment and governance decisions of watershed soil erosion.

[0097] like Figure 2 As shown, the present invention constructs a deep learning model for dynamic monitoring of soil erosion. Based on the physical characteristics of the watershed hydrological processes, a graph neural network architecture integrating physical priors is designed. Its core ideas include: (1) Topological constraints: defining node associations based on hydrological connectivity rather than spatial proximity; (2) Physical guidance: determining information transmission intensity based on the physical laws of topography and rainfall rather than pure data-driven approaches; (3) Spatiotemporal integration: simultaneously modeling monthly dynamic changes and annual cumulative effects. The model consists of four parts: an input embedding layer, a physical injection graph attention layer, a temporal aggregation layer, and an erosion modulus output layer. The input embedding layer integrates the dynamic and static features obtained in step 1 and adds time encoding; the physical injection graph attention layer adds the topographic and rainfall physical weights generated based on the D8 algorithm to the attention score calculation formula of the graph attention network (GAT) in the form of a logarithmic bias term; the temporal aggregation layer compresses the monthly data of 12 months into annual comprehensive features; and the output layer finally outputs the annual soil erosion modulus for the entire watershed, grid-by-grid. In the model training phase (step 4), the relative ranking relationship of the RUSLE model is learned using ranking consistency loss, and combined with the measured correction loss of the actual measurement stations and the terrain adaptive spatial smoothing loss for composite supervised training. In the model inference and result generation phase (step 5), forward propagation is performed on the annual data to be monitored, and the annual soil erosion modulus of the entire watershed is output grid by grid, and the thematic map of erosion intensity is generated according to the water conservancy industry standards.

[0098] To verify the effectiveness of the deep learning model for dynamic monitoring of soil erosion described in this invention, two sets of comparative experiments were designed: Experiment 1 verifies the overall performance advantages of this invention compared to existing technologies; Experiment 2 verifies the independent contributions of each innovative module through an ablation experiment. The experiments were conducted in a typical sub-basin of the Yellow River Basin (approximately 5000 km²). 2 The grid consisted of N=12,500 grids, including 12 hydrological stations and 8 runoff plots. Data from 2015 to 2019 was used (training in the first 3 years, validation in the 4th year, and testing in the 5th year). Rusle calculations were used as weakly supervised labels, and the measured station data were used as strongly supervised labels. All experiments were independently repeated 5 times under the same hardware and software environment (Python 3.8, PyTorch 1.12, NVIDIA RTX 3090), and the average performance metrics are reported.

[0099] Experiment 1: Performance Comparison with Existing Technologies Experimental Design: The present invention was compared with three existing methods as shown in Table 1 below: (1) Traditional physical model RUSLE (Modified General Soil Loss Equation); (2) Traditional machine learning methods (Random Forest RF, Support Vector Machine SVM); (3) Pure data-driven deep learning methods (Standard GAT, Spatiotemporal Graph Convolutional Network STGCN). All comparison methods were run under the same training / test partition, and the input features were uniformly the extracted dynamic and static feature tensors to ensure fair comparison.

[0100] Table 1 Comparison of soil erosion modulus prediction performance of different models

[0101] Note: ↑ indicates higher is better, ↓ indicates lower is better; NSE is the Nash efficiency coefficient, R 2 The coefficient of determination is RMSE, and the mean square error is MAE. The measured match rate refers to the percentage of stations where the relative error between the predicted and measured values ​​is less than 20%. All indicators are the average of 5 independent runs, with a standard deviation of less than 0.03.

[0102] Experiment 2: Ablation Experiment Verifies the Contribution of Each Module Experimental Design: Based on the deep learning model for dynamic monitoring of soil erosion of the present invention, four ablation versions were constructed by removing key modules one by one: (1) No-Physics: attention calculation only retains the data-driven term and removes the logarithmic physical bias; (2) No-Temporal: the temporal aggregation layer is changed to simple averaging and the learnable attention weights are removed; (3) No-RankLoss: only the measured correction loss and smoothing loss are used for training; (4) No-SmoothLoss: only the ranking loss and measured correction loss are used for training. All ablation versions maintain the same network depth (L=2) and hyperparameters, and differ from the deep learning model for dynamic monitoring of soil erosion of the present invention only in specified modules.

[0103] Table 2 Comparison of Ablation Test Performance

[0104] Note: No-Physics and No-Temporal are structural ablation methods, while No-RankLoss and No-SmoothLoss are loss function ablation methods; all metrics are the average of 5 independent runs.

[0105] In summary, the solution of this invention can solve the challenges of soil erosion monitoring involving both spatial topology modeling and sparse supervised learning. It injects hydrological and physical laws (flow direction, topography, rainfall) as "prior knowledge" into the deep learning model. By constructing a time-varying physical directed graph constrained by the DEM, it solves the spatial topology modeling problem and ensures the physical rationality of hydrological connectivity. At the algorithm level, the logarithm of the physical weights is added as a bias term to the attention score, achieving physics-data synergy. During model training, the reliable spatial distribution trend (relative ranking) of RUSLE is extracted as a weak supervision signal, combined with strong supervision from sparse measured data for composite training, solving the problem of sparse supervised learning and enabling the training of a high-precision model even with insufficient data. The core improvements of this invention compared to existing graph neural networks are threefold: Improvement 1: Construction of time-varying physical directed graph constrained by DEM Existing graph neural networks typically employ undirected graphs or adjacency relationships based on spatial distance, failing to represent the directional characteristics of water flow in a watershed. This invention constructs a strictly downstream directed graph based on D8 flow direction analysis, with edge weights jointly determined by static topographic weights and time-varying rainfall driving factors. This graph structure encodes the topological constraints and seasonal connectivity variations of watershed hydrological transmission, providing a physically sound topological basis for message passing. The core lies in combining the strictly downstream DAG constructed by the D8 algorithm with time-varying rainfall-driven physical weights, specifically designed for graph neural network message passing in soil erosion monitoring.

[0106] Improvement 2: Attention Injection Mechanism Based on Log-Bias of Physical Weights Existing graph attention networks rely entirely on data-driven learning for their attention weights, failing to incorporate known hydrophysical laws. This invention adds the logarithmic value of the time-varying physical weights calculated in step S2 as a bias term to the attention score. This design allows the attention weights to be synergistically controlled by both data-driven learning (query-key dot product) and physical prior constraints (logarithmic bias): edges with strong physical connectivity receive attention enhancement, while edges with weak connectivity are automatically masked, simultaneously preserving the model's ability to learn complex erosion patterns from remote sensing big data.

[0107] Improvement 3: Weakly supervised training strategy for maintaining order Unlike directly using RUSLE values ​​as labels (inheriting their systematic bias) or relying solely on sparse measured data (insufficient supervision signals), this invention employs ranking consistency loss to extract reliable spatial distribution knowledge (relative ranking relationships) from RUSLE, and combines it with numerical calibration loss from measured stations and terrain-adaptive spatial smoothing loss to construct a composite supervision framework suitable for sparse ground data conditions.

[0108] The synergistic effect of the three improvements—improvement one provides a physically reasonable graph structure, improvement two enables physical-data collaborative message passing on this structure, and improvement three establishes a reasonable supervisory balance between physical priors and measured data—can solve the core challenges of "spatial topology modeling" and "sparse supervised learning" in soil erosion monitoring.

[0109] Example 2:

[0110] This application provides a dynamic monitoring system for soil erosion based on remote sensing technology, comprising: the system using the method described above to dynamically monitor soil erosion in the monitored watershed.

[0111] It should be understood that the phrase "one embodiment" or "an embodiment" throughout the specification means that a specific feature, structure, or characteristic related to the embodiment is included in at least one embodiment of this application. Therefore, "in one embodiment" or "in an embodiment" appearing throughout the specification does not necessarily refer to the same embodiment. Furthermore, these specific features, structures, or characteristics can be combined in any suitable manner in one or more embodiments. It should be understood that in the various embodiments of this application, the sequence numbers of the above-described processes do not imply a sequential order of execution; the execution order of each process should be determined by its function and internal logic, and should not constitute any limitation on the implementation process of the embodiments of this application. The sequence numbers of the above-described embodiments are merely descriptive and do not represent the superiority or inferiority of the embodiments.

[0112] It should be noted that, in this document, the terms "comprising," "including," or any other variations thereof are intended to cover non-exclusive inclusion, such that a process, method, article, or apparatus that comprises a list of elements includes not only those elements but also other elements not expressly listed, or elements inherent to such a process, method, article, or apparatus. Unless otherwise specified, an element defined by the phrase "comprising one..." does not exclude the presence of other identical elements in the process, method, article, or apparatus that includes that element.

[0113] In the several embodiments provided in this application, it should be understood that the disclosed devices and methods can be implemented in other ways. The device embodiments described above are merely illustrative. For example, the division of units is only a logical functional division, and in actual implementation, there may be other division methods, such as: multiple units or components can be combined, or integrated into another system, or some features can be ignored or not executed. In addition, the coupling, direct coupling, or communication connection between the various components shown or discussed can be through some interfaces, and the indirect coupling or communication connection between devices or units can be electrical, mechanical, or other forms.

[0114] The units described above as separate components may or may not be physically separate. The components shown as units may or may not be physical units. They may be located in one place or distributed across multiple network units. Some or all of the units may be selected to achieve the purpose of the embodiments of this application, depending on actual needs.

[0115] In addition, each functional unit in the various embodiments of this application can be integrated into one processing unit, or each unit can be a separate unit, or two or more units can be integrated into one unit; the integrated unit can be implemented in hardware or in the form of hardware plus software functional units.

[0116] Alternatively, if the integrated units described above are implemented as software functional modules and sold or used as independent products, they can also be stored in a computer-readable storage medium. Based on this understanding, the technical solutions of the embodiments of this application, or the parts that contribute to related technologies, can be embodied in the form of a software product. This computer software product is stored in a storage medium and includes several instructions to cause the device automatic test line to execute all or part of the methods described in the various embodiments of this application. The aforementioned storage medium includes various media capable of storing program code, such as mobile storage devices, ROMs, magnetic disks, or optical disks.

[0117] The methods disclosed in the several method embodiments provided in this application can be arbitrarily combined without conflict to obtain new method embodiments.

[0118] The features disclosed in the several method or device embodiments provided in this application can be arbitrarily combined without conflict to obtain new method or device embodiments.

[0119] The above description is merely an embodiment of this application, but the scope of protection of this application is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in this application should be included within the scope of protection of this application. Therefore, the scope of protection of this application should be determined by the scope of the claims.

Claims

1. A method for dynamic monitoring of water and soil erosion based on remote sensing technology, characterized in that, include: Extract the dynamic features that change monthly and the static features that remain unchanged throughout the year from the multi-source remote sensing data, digital elevation model data, meteorological data, soil property data and ground measured data collected in the monitoring basin for the year to be monitored. Based on the dynamic features, a dynamic feature tensor is constructed; based on the static features, a static feature tensor is constructed. Based on the digital elevation model, the D8 single-flow direction algorithm is used to construct a directed acyclic graph in the direction of flow. In this graph, multiple grids of the monitored watershed are used as multiple nodes, the direction of water flow is used as directed edges, and grids in the watershed whose slope to all neighboring directions is not greater than zero are marked as sinks. The dynamic feature tensor and the static feature tensor are input into the trained deep learning model for dynamic monitoring of soil erosion to obtain the annual soil erosion modulus for each grid in the entire watershed. The deep learning model for dynamic monitoring of soil erosion is based on a physical information graph neural network constructed with the directed acyclic graph as topological constraint, and is obtained through composite supervised training using a combination of sorting consistency loss, measured correction loss of actual measurement stations, and terrain adaptive spatial smoothing loss.

2. The method of claim 1, wherein, The physical information graph neural network includes an input embedding layer, a physical injection graph attention layer, a temporal aggregation layer, and an erosion modulus output layer; The pre-construction process of the deep learning model for dynamic monitoring of soil erosion includes: The historical dynamic feature tensor and historical static feature tensor, as well as the corresponding historical annual soil erosion modulus for each grid in the whole watershed, are obtained and the corresponding directed acyclic graph is constructed. For each directed edge, the static terrain weight and time-varying rainfall driving factor are calculated and then normalized to obtain the time-varying normalized adjacency matrix sequence. The historical dynamic feature tensor and historical static feature tensor are fused through the input embedding layer of the physical information graph neural network and time encoding is added. Message passing and feature update processing are performed through the physical injection graph attention layer of the physical information graph neural network. During message passing, the logarithm of the physical weights in the time-varying normalized adjacency matrix sequence is used as a bias term and superimposed on the attention score calculation formula of the graph attention network to apply physical prior constraints to upstream neighbor edges. The temporal aggregation layer of the physical information graph neural network automatically weights the monthly data of 12 months under each node into annual data. The erosion modulus output layer of the physical information graph neural network outputs the historical annual soil erosion modulus for the entire watershed, grid-by-grid. Based on a composite loss function constructed using sorting consistency loss, measured correction loss, and terrain-adaptive spatial smoothing loss, the physical information graph neural network is trained using the historical dynamic feature tensor, the historical static feature tensor, and the historical annual soil erosion modulus to obtain the deep learning model for dynamic monitoring of soil and water loss.

3. The method of claim 1, wherein, The method for constructing a directed acyclic graph in the direction of flow using the D8 unidirectional flow algorithm based on the digital elevation model includes: dividing the monitored watershed into multiple grids, which serve as multiple nodes; Depression filling was performed on the digital elevation model of the monitored watershed to obtain a hydrologically consistent elevation surface; Based on the hydrologically consistent elevation surface, the D8 single-flow direction algorithm is used to determine the flow direction of each node: for each node, the slope from the node to its 8 neighboring nodes is calculated, and the direction with the largest and positive slope is determined as the flow direction. If the slope of a node to all neighborhood directions is not greater than zero, then mark it as a sink. Based on the sink, the direction of water flow, and the slope of each node to its 8 neighboring nodes, the cumulative flow of each node is recursively calculated in a topological order from upstream to downstream. A directed acyclic graph is constructed based on each node, the direction of water flow corresponding to the node as the directed edge of the node, the sink and the cumulative flow.

4. The method according to claim 2, characterized in that, The process involves calculating the physical weights for each directed edge, fusing static terrain weights and time-varying rainfall driving factors. After normalization, a time-varying normalized adjacency matrix sequence is obtained, including: For the slope from a node to a neighboring node in the directed acyclic graph, a static terrain weight calculation model is used to calculate the fused static terrain weight. For the nodes in the directed acyclic graph, the time-varying rainfall driving factor calculation model is used to calculate the time-varying rainfall driving factor; The product of the fused static terrain weight and the time-varying rainfall driving factor is determined as the physical weight of the directed edge from the node to the neighboring nodes. Normalize the in-degree of the target node for the physical weight of each directed edge, and output a time-varying normalized adjacency matrix sequence. The static terrain weight calculation model satisfies the following formula: In the above formula, To integrate static terrain weights, For nodes To neighboring nodes The slope, The average gradient of the watershed; For nodes The cumulative amount of the confluence, This represents the maximum cumulative runoff in the basin. The calculation model for the time-varying rainfall driving factor satisfies the following formula: In the above formula, As a driving factor for time-varying rainfall, For nodes In monthly time steps The erosive power of monthly rainfall, For nodes The maximum monthly rainfall erosion during the training period.

5. The method according to claim 3, characterized in that, The step of normalizing the in-degree of the target node for the physical weight of each directed edge and outputting a time-varying normalized adjacency matrix sequence includes: For the physical weight of each directed edge, the sum of the physical weights of all incoming edges of the downstream target node of the directed edge is used as the normalized denominator. The physical weight is divided by the denominator to obtain the normalized weight value of the directed edge. After performing the above normalization operation on all directed edges, we obtain the normalized weight values ​​of all directed edges. The normalized weights of all directed edges are combined to obtain a time-varying normalized adjacency matrix sequence.

6. The method according to claim 2, characterized in that, The process of fusing the historical dynamic feature tensor and the historical static feature tensor through the input embedding layer of the physical information graph neural network and adding time encoding includes: For each node in the input embedding layer of the physical information graph neural network at each time step, dynamic feature projection and static feature projection are performed to obtain the projected dynamic feature vector and the projected static feature vector. The projected dynamic feature vector is concatenated with the projected static feature vector, and then subjected to a nonlinear transformation through the fully connected layer of the physical information graph neural network. Time encoding is superimposed on the results of the nonlinear transformation.

7. The method according to claim 2, characterized in that, The message passing and feature update processing through the physical injection graph attention layer of the physical information graph neural network includes: For each downstream node of the physical injection graph attention layer, the union of the upstream neighbor set of the downstream node and the self-loop is defined as the message aggregation neighborhood. The message aggregation neighborhood is used to limit the range of upstream nodes that the downstream node considers when updating its own state. The logarithmic value of the physical weight in the time-varying normalized adjacency matrix sequence is used as a logarithmic bias term. This is then superimposed on the upstream neighbor edges into the attention score calculation formula of the graph attention network of the physical information graph neural network. This is added to the data-driven scaling dot product attention term of the graph attention network to form the attention score of each upstream node in the message aggregation neighborhood corresponding to the downstream node, thus updating the downstream node. The attention scores of each upstream node in the message aggregation neighborhood of the updated downstream node are normalized using Softmax to obtain normalized attention scores; the normalized attention scores are then converted into probability distributions to obtain the attention ratio of each upstream node to the downstream node in the message aggregation neighborhood. Using the attention ratio as the weight, the value vectors of each node in the message aggregation neighborhood are weighted and summed to obtain the aggregated message received by the downstream node from each attention head; The aggregated messages from multiple attention heads are concatenated and multi-view information is fused using a linear projection matrix to obtain the attention output of the physical injection graph attention layer. In the attention layer of the physical injection graph, a Post-Norm residual structure is adopted. The attention output is randomly deactivated and then added to the input residual. After layer normalization, the intermediate features are obtained. In the attention layer of the physical injection graph, the intermediate features are nonlinearly transformed by a feedforward network. The feedforward network is a two-layer fully connected network. The first layer expands the feature dimension and then uses the GELU activation function. The second layer restores the feature dimension to the original dimension, thus obtaining the spatiotemporal node feature representation of the physical injection graph attention layer output by fusing upstream multi-order neighborhood information and physical constraints. The spatiotemporal node feature representation output by the physical injection graph attention layer is used as the input of the next physical injection graph attention layer. By stacking multiple physical injection graph attention layers, each downstream node aggregates the feature information of its upstream multi-order neighborhood.

8. The method according to claim 2, characterized in that, The temporal aggregation layer of the physical information graph neural network automatically weights the monthly data for each node into annual data for the 12 months below; the erosion modulus output layer of the physical information graph neural network outputs the historical annual soil erosion modulus for the entire watershed, grid by grid, including: The time-series aggregation layer employs a learnable time attention mechanism, which sums the time weights of each node at each monthly time step to compress monthly features into annual comprehensive features. The annual comprehensive features of the erosion modulus output layer are mapped sequentially through two fully connected networks, and a smooth non-negative activation function is used to constrain the output to be a positive value for the historical annual soil erosion modulus of the entire watershed grid by grid.

9. The method according to claim 2, characterized in that, The composite loss function satisfies the following formula: In the above formula, For compound loss, For sorting consistency loss, To measure and correct the loss, For terrain-adaptive spatial smoothing loss, Weights are lost for consistency in sorting. To measure and correct the loss weights, The terrain-adaptive spatial smoothing loss weights are used; The order consistency loss satisfies the following formula: In the above formula, For consistent sampling of sample pairs, For each batch, two nodes are randomly selected from N grids across the entire watershed to form an ordered pair. For sorting interval, This is the model's first predicted value. This is the model's second predicted value; The measured correction loss satisfies the following formula: In the above formula, M represents the total number of measured stations. For the first Grid index of each measured site, For the model to the first Predicted values ​​for each measured site For the first Measured erosion modulus of each actual measurement site; The terrain-adaptive spatial smoothing loss satisfies the following formula: In the above formula, The set of edges between adjacent grids in a graph structure. Let the total number of edges be . Indicates two adjacent grid nodes. For nodes The slope, For neighboring nodes The slope, The attenuation coefficient is... For terrain similarity weights, For nodes The erosion modulus, For nodes The erosion modulus.

10. A dynamic monitoring system for soil erosion based on remote sensing technology, characterized in that, include: The system uses the method described in any one of claims 1 to 9 to perform dynamic monitoring of soil and water loss in the monitored watershed.

Images