A method for dynamic monitoring and evaluation of tributary sediment yield based on multi-source remote sensing data

Abstract

The application discloses a kind of tributary sediment yield dynamic monitoring and evaluation method based on multi-source remote sensing data fusion, the method is by obtaining the C band VV polarization SAR image and optical multispectral image of the same phase of target basin, cloud pollution pixel is identified using scene classification layer SCL;With the training sample of synchronous optical and SAR data in historical clear sky period, pixel-by-pixel matching least square regression establishes basin special SAR-SSC inversion equation;With DEM, sub-basin is demarcated, cloud pollution pixel is obtained using SAR-SSC equation Suspended sediment concentration, with clear sky optical SSC unified into sub-basin node feature matrix;Directed hydrological topological graph and space-time graph neural network are constructed, and NSE is used as core loss function training, to realize the space-time inversion of sediment yield;The application realizes the dynamic monitoring of sediment yield under all-weather, all-cloud and rain conditions, and fundamentally eliminates the observation blind area during heavy rain period.

Description

Technical Field

[0001] This invention belongs to the field of remote sensing hydrology technology, specifically relating to a method for dynamic monitoring and assessment of tributary sediment yield based on multi-source remote sensing data. Background Technology

[0002] Watershed sediment yield is a core parameter for soil and water conservation, reservoir siltation prediction, and river ecological management. The rainy season is the period of highest concentration and intensity of sediment yield, often contributing over 70% of the annual yield. However, the exact moment of the rainstorm coincides with the thickest cloud cover and the heaviest rainfall, which is also when optical satellites become completely ineffective. This observational blind spot highly overlaps with the peak sediment yield in time, resulting in a systematic lack of crucial sediment yield data in existing optical remote sensing methods. This leads to a severe underestimation of sediment yield, rendering the monitoring results almost useless in constraining rainstorm and flood events.

[0003] Currently, optical satellite remote sensing (such as Sentinel-2 and Landsat) has developed a relatively mature technology system for retrieving suspended sediment concentration (SSC) from water reflectance. Algorithms such as the red-green band ratio method can achieve high-precision SSC spatial mapping under clear-sky conditions. However, optical satellites have a fundamental physical limitation: they are completely opaque to clouds. Synthetic Aperture Radar (SAR) operates in the microwave band, and its signals can penetrate clouds and moderate rainfall, providing all-weather imaging capabilities. Studies have shown that there is a quantifiable physical correlation between the VV polarization backscattering coefficient of SAR and the suspended sediment concentration in water bodies, providing a theoretical basis for SSC retrieval. However, there is currently a lack of systematic methods to fully utilize synchronized optical and SAR data during clear-sky periods, establish reliable, basin-specific quantitative SAR-SSC retrieval models at the watershed scale, and put them into practical use during periods of heavy rain and cloud cover. Although the surface hydrological station network can provide continuous observation of cross-sectional flow and sediment transport rate, the station density is seriously insufficient, making it impossible to obtain spatial distribution information of sediment production processes. Moreover, instrument damage and measurement stoppages frequently occur under extreme flood conditions.

[0004] Most existing multi-source remote sensing fusion methods propose a general idea of ​​multi-source complementarity, but they lack a targeted system solution for the core contradiction that the peak of rainstorm sediment production is precisely when optical signals are completely ineffective. Specifically, there is a lack of operable methods to directly obtain high-quality SSCs under cloud and rain cover conditions, as well as a technical path to establish a watershed-specific SAR quantitative inversion model using limited clear-day data.

[0005] Therefore, there is an urgent need for a systematic method that can continuously acquire SSC information of the entire watershed under rainstorm and cloud cover conditions, and thereby realize dynamic monitoring and assessment of sediment yield. Summary of the Invention

[0006] The purpose of this invention is to address the core problem that existing optical remote sensing becomes completely ineffective during peak rainstorm sediment production periods due to cloud cover, resulting in a systematic data gap in sediment production monitoring. This invention provides a method for spatiotemporal inversion of suspended sediment concentration and dynamic monitoring of sediment production in watersheds based on the fusion of SAR and multi-source optical remote sensing.

[0007] To achieve the above objectives, this invention provides a method for dynamic monitoring and assessment of tributary sediment yield based on multi-source remote sensing data, the method comprising the following steps: Step S1: Acquire the first satellite C-band VV polarimetric synthetic aperture radar SAR image and the second satellite optical multispectral image covering the same time phase of the target watershed, as well as the measured flow and suspended sediment concentration SSC observation sequences of various surface hydrological stations in the watershed; using the scene classification layer SCL band of the second satellite, mark pixels with SCL values ​​of 8, 9 or 10 as cloud pollution pixels, and mark the remaining pixels falling into the water body range as effective optical pixels.

[0008] Step S2: Using the first satellite SAR image and the second satellite optical image acquired during the same period of historical clear skies in the basin as training samples; the suspended sediment concentration (SSC) in the water body is calculated using the red-green band ratio method for the effective optical pixels; the VV polarization backscattering coefficient of the SAR image at the same spatial location is paired with the suspended sediment concentration (SSC) pixel by pixel; and the SAR-SSC inversion equation of the basin is established by least squares regression using no less than 30 pairs of effective observation data.

[0009] Step S3: Delineate the sub-basin boundaries using the target watershed digital elevation model (DEM) and divide the target watershed into N sub-basin units. For all water body pixels within each sub-basin, calculate the optical suspended sediment concentration using the red-green band ratio method for effective optical pixels, and calculate the SAR suspended sediment concentration using the established SAR-SSC inversion equation for cloud-polluted pixels. Use the area-weighted average of all water body pixels within each sub-basin as the suspended sediment concentration feature value for that sub-basin node. Simultaneously, extract the Normalized Differential Vegetation Index (NDVI), Topographic Humidity Index (TWI), and the 24-hour cumulative rainfall from the nearest meteorological station for each sub-basin as auxiliary node features, which together constitute the node feature matrix.

[0010] Step S4: Using each sub-basin as a graph node, directed edges are generated based on the water flow direction extracted from the DEM. The function value of the channel length and average slope between adjacent sub-basins is used as the weight of the directed edges to construct a directed hydrological topology graph representing the upstream and downstream water and sediment transport relationship. A spatiotemporal graph neural network is constructed, consisting of a two-layer directed graph convolution module and a single-head temporal attention module connected in series. The measured sediment yield of each surface hydrological station is used as the supervision signal. The training set, validation set, and test set are divided in a 7:1.5:1.5 ratio according to time sequence. The spatiotemporal graph neural network is trained using a loss function with the negative Nash-Sutcliffe efficiency coefficient (NSE) as the core.

[0011] Step S5: The trained spatiotemporal graph neural network is used to perform real-time sediment yield inversion for each sub-basin of the target watershed. During the inference phase, the Dropout layer is kept active, and 50 random forward propagations are performed. The mean of the 50 predicted values ​​is used as the estimated sediment yield, and the coefficient of variation of the 50 predicted values ​​is used as the uncertainty index. Sub-basins with a coefficient of variation exceeding 0.3 are automatically marked as high uncertainty areas and a supplementary sampling request is triggered. The node features of the corresponding sub-basin are updated with the latest SAR or optical imagery from the supplementary sampling, and then local re-inference is performed to output the updated estimated sediment yield and confidence interval.

[0012] Furthermore, optically suspended sediment concentration The calculation formula is: ;in, The surface reflectance in the red band of the second satellite. denoted as the surface reflectance in the green band of the second satellite, and a and b are watershed-specific coefficients determined by logarithmic linear regression using ground-measured suspended sediment concentration data.

[0013] The SAR–SSC inversion equation is in the form of: ;in, It is the VV polarization backscattering coefficient (linear dimension). The local angle of incidence at the water surface pixel, in degrees. The least squares regression coefficient is used; the selection criterion for water body pixels is that the normalized difference water index (NDWI) is greater than 0.3.

[0014] Furthermore, in step S2, before establishing the SAR–SSC inversion equation, the following preprocessing is performed sequentially on the first satellite SAR image: Using the calibration vector in the metadata of the first satellite product, the original digital quantization (DN) value is converted into sigma-naught backscattering coefficients; Using a 30 m resolution as a reference, a range-Doppler terrain correction algorithm is used to eliminate geometric distortions caused by terrain undulations and orthorectify the SAR image to the same UTM projection coordinate system as the second satellite image. A 5×5 window Lee filter was applied to the radiometrically calibrated backscattering coefficient image to suppress speckle noise; The filtered backscattering coefficient image histogram is segmented using Otsu thresholding. Pixels below the threshold are identified as water surface pixels, generating a water surface mask to limit the spatial range of subsequent SAR-SSC paired samples.

[0015] Furthermore, in step S3, the terrain humidity index TWI is: ;in, The unit is the catchment area, in m² / m. The local slope is extracted from the DEM; the TWI of each sub-basin is the median of the TWI values ​​of all raster cells within that sub-basin.

[0016] For a sub-basin with a catchment area of ​​less than 1 km², if there are fewer than 5 effective water body pixels within it, the suspended sediment concentration characteristic value of the sub-basin is replaced by the area-weighted average of the suspended sediment concentration characteristic values ​​of its directly upstream adjacent sub-basins. The weight is the ratio of the confluence contribution area of ​​the corresponding upstream sub-basin to the sub-basin.

[0017] Furthermore, in step S4, the directed hydrological topology map is composed of upstream sub-basins. Pointing to downstream sub-basins Calculation of directed edge weights: ;in, The length of the river channel from the center point of the upstream sub-basin u to the center point of the downstream sub-basin v along the confluence path extracted by the DEM, in meters. The average gradient of the confluence path is expressed in m / m. For source sub-basin nodes without upstream adjacent sub-basins, their graph convolutional aggregation only uses their own node features and does not participate in cross-node message passing.

[0018] Furthermore, in step S4, the spatial graph convolutional module of the spatiotemporal graph neural network adopts an in-degree normalization aggregation strategy. Layer nodes The feature update formula is: ;in, Let v be the set of upstream neighboring sub-basin nodes of node v in a directed hydrological topology graph. The directed edge weights are normalized to the interval [0,1]. and The first The trainable weight matrix and bias vector of layer graph convolution, where CONCAT represents the vector concatenation operation.

[0019] The temporal attention module takes a 7-day node feature time series as input, performs a single-layer linear transformation, and then uses a softmax function to calculate the attention weights at each time point. The weighted sum is then input into the fully connected layer, and the predicted sediment yield of each sub-basin is output.

[0020] Furthermore, in step S4, the spatiotemporal graph neural network is trained using the following loss function: ;in, The predicted sediment yield values ​​for each sub-basin at each observation station at time t are: To correspond to the measured sediment yield, The training set is the time-averaged measured sediment yield; Training was performed using the Adam optimizer with an initial learning rate of 1×10⁻⁶. - ³, using a cosine annealing strategy to decay the learning rate to 1×10⁻⁶ over 200 training epochs. -5 After each graph convolution, a dropout rate of p=0.1 is applied to achieve both regularization and Monte Carlo uncertainty estimation during the inference stage.

[0021] Furthermore, in step S5, the estimated sediment yield and uncertainty index are calculated as follows: Perform 50 random forward propagations while maintaining Dropout activation. The prediction value for sub-basin S after the i-th propagation is denoted as... Forecast mean With the predicted standard deviation calculate: ; ;by Output as the 95% approximate confidence interval for sediment yield; coefficient of variation The calculation formula is: ,right If a sub-basin exceeds 0.3, triggering a supplementary imagery request, the latest first satellite imagery covering the sub-basin and within 3 days of the trigger time will be selected first. If no new imagery is acquired within 72 hours, the original predicted mean will be retained and marked "Insufficient data, doubtful estimation" in the results.

[0022] Furthermore, in step S4, for sub-basins lacking ground observation stations, training labels are constructed as follows: Measured sediment yield in downstream sub-basins with observation stations Based on this, the training label sediment yield of the station-free sub-basin is estimated using the following formula. : ;in, and These represent the catchment areas of the sub-basin without monitoring stations and the downstream sub-basin with monitoring stations, respectively, in km². and These are the surface erosion intensity indices for the corresponding sub-basins.

[0023] The estimated sediment yield is only used for auxiliary supervision during the training phase. It is assigned a weight coefficient of 0.3 in the loss function, while the measured sediment yield of the sub-basin with observation stations has a corresponding weight coefficient of 1.0, in order to distinguish the confidence difference between the measured data and the estimated data.

[0024] Furthermore, the method also includes an online rolling update step for the model, specifically: using a 30-day rolling time window, the root mean square error (RMSE) of the spatiotemporal graph neural network for predicting sediment yield in sub-basins with observation stations is statistically analyzed; when the RMSE exceeds 20% of the average measured sediment yield within the same window, a model update is triggered.

[0025] During model updates, the training set is the observation data from the 365 days prior to the current time. All weight parameters of the first layer graph convolution of the spatiotemporal graph neural network are frozen, and only the parameters of the second layer graph convolution and the temporal attention module are fine-tuned. The fine-tuning uses an initial learning rate of 1×10⁻⁶. -4 The Adam optimizer was trained for 50 epochs.

[0026] After fine-tuning, if the updated model's NSE on the validation set is not lower than that of the original model, the updated model will replace the original model and be put into use; otherwise, the original model will be retained and an anomaly log will be recorded for manual review.

[0027] Compared with the prior art, the beneficial effects of the present invention are: This invention utilizes the physical property that SAR radar signals can penetrate clouds and rainfall. Using SAR and optical image pairs acquired synchronously during historical clear-sky periods as training samples, a basin-specific SAR-SSC quantitative inversion equation is established. This allows for the acquisition of spatially continuous suspended sediment concentration information even during peak periods of heavy rainfall and sediment yield when optical satellites are completely ineffective. Furthermore, by combining a directed hydrological topology map with a spatiotemporal neural network, real-time dynamic estimation and uncertainty quantification of sediment yield across the entire basin can be achieved. Attached Figure Description

[0028] Figure 1 This is a remote sensing image of the same area as SAR and optical images in this invention; Figure 2 This is a schematic diagram of the DEM and directed hydrological topology map of the present invention. Detailed Implementation

[0029] To make the objectives, technical solutions, and advantages of this invention clearer, the technical solutions of this invention are described clearly and completely below. Obviously, the described embodiments are only a part of the embodiments of this invention, not all of them. Based on the embodiments of this invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of this invention.

[0030] A dynamic monitoring and assessment method for tributary sediment yield based on multi-source remote sensing data is presented, with a specific application example of a river basin. This river is a first-order tributary on the right bank of the middle reaches of the Yellow River, with a catchment area of ​​approximately 8706 km², a main channel length of 242 km, and a multi-year average sediment transport modulus of approximately 18200 t / (km²·a). It is one of the tributaries with the most severe soil erosion on the Loess Plateau.

[0031] Step S1: Acquire the first satellite C-band VV polarimetric synthetic aperture radar SAR image and the second satellite optical multispectral image covering the same time phase of the target watershed, as well as the measured flow and suspended sediment concentration SSC observation sequences of various surface hydrological stations in the watershed; using the scene classification layer SCL band of the second satellite, mark pixels with SCL values ​​of 8, 9 or 10 as cloud pollution pixels, and mark the remaining pixels falling into the water body range as effective optical pixels.

[0032] like Figure 1 The diagram shows remote sensing images of the same area where SAR and optical images are obtained. A directed flowchart illustrates the processing chain from multi-source remote sensing data acquisition to real-time sediment yield inversion. SAR and optical images acquired synchronously during clear-sky periods are preprocessed to jointly establish the SAR-SSC inversion equation. Imagery at the target time undergoes cloud detection and diversion processing; cloud-polluted pixels are supplemented with suspended sediment concentration using the SAR-SSC equation and incorporated into the sub-basin node feature matrix along with the clear-sky optical SSC. The feature matrix is ​​then processed by a directed hydrological topology map and a spatiotemporal neural network to output estimated sediment yield values ​​and 95% confidence intervals for each sub-basin. When prediction uncertainty exceeds a threshold, the system automatically triggers supplementary data collection and performs local re-inference, ultimately achieving a three-layer complementary all-weather sediment yield monitoring architecture of optical-SAR-neural network.

[0033] The target monitoring time was July 14, 2023. At 15:12 (Beijing time) on that day, Sentinel-1A passed over the target river basin in an ascending orbit and acquired C-band VV polarized GRD products with a ground distance resolution of 10 m. At 15:47, Sentinel-2B passed over and acquired L2A atmospheric corrected multispectral images. The passing times of the two satellites were about 35 minutes apart, which met the requirements for simultaneous phase registration.

[0034] Three national basic hydrological stations within the basin provided synchronously measured sequences of daily average flow and suspended sediment concentration (SSC) from 2016 to 2023, with a cumulative total of 2,372 valid station records.

[0035] Cloud detection was performed on the SCL band of Sentinel-2B: the area of ​​cloud-polluted pixels with SCL values ​​of 8, 9, or 10 accounted for 87.3% of the total water body pixels in the entire basin, leaving only 12.7% of the effective optical pixels, mainly concentrated in the high-altitude area of ​​the source on the southwest side of the basin. This detection result quantitatively confirms the serious limitations of optical remote sensing during rainstorm periods, and also indicates that this embodiment mainly relies on Sentinel-1A SAR imagery to obtain information on suspended sediment in the entire basin.

[0036] Step S2: Using the first satellite SAR image and the second satellite optical image acquired during the same period of historical clear skies in the basin as training samples; the suspended sediment concentration (SSC) in the water body is calculated using the red-green band ratio method for the effective optical pixels; the VV polarization backscattering coefficient of the SAR image at the same spatial location is paired with the suspended sediment concentration (SSC) pixel by pixel; and the SAR-SSC inversion equation of the basin is established by least squares regression using no less than 30 pairs of effective observation data.

[0037] Before establishing the SAR–SSC inversion equation, the following preprocessing steps are performed on the first satellite SAR image: Using the calibration vector in the metadata of the first satellite product, the original digital quantization (DN) values ​​are converted into sigma-naught backscattering coefficients. With a resolution of 30 m as a reference, a range-Doppler terrain correction algorithm is used to eliminate geometric distortions caused by terrain undulations, and the SAR image is orthorectified to the same UTM projection coordinate system as the second satellite image. A 5×5 window Lee filter is applied to the radiometrically calibrated backscattering coefficient image to suppress speckle noise. The histogram of the filtered backscattering coefficient image is segmented using Otsu thresholding, and pixels below the threshold are identified as water surface pixels, generating a water surface mask to limit the spatial range of subsequent SAR-SSC paired samples.

[0038] Optical suspended sediment concentration The calculation formula is: ;in, The surface reflectance in the red band of the second satellite. denoted as the surface reflectance in the green band of the second satellite, and a and b are watershed-specific coefficients determined by logarithmic linear regression using ground-measured suspended sediment concentration data.

[0039] The SAR–SSC inversion equation is in the form of: ;in, It is the VV polarization backscattering coefficient (linear dimension). The local angle of incidence at the water surface pixel, in degrees. The least squares regression coefficient is used; the selection criterion for water body pixels is that the normalized difference water index (NDWI) is greater than 0.3.

[0040] From March 2021 to November 2022, during non-rainy, clear-sky periods, 41 transit scenes with Sentinel-2B cloud cover below 5% were selected. Each scene was matched with Sentinel-1A orbital imagery within a time difference of no more than 6 hours. Radiometric calibration, range-Doppler topographic correction (referencing SRTM 30 m DEM, orthorectified to EPSG:32649 coordinate system), and 5×5 window Lee filtering were performed on the SAR images. Water surface masks were extracted using Otsu thresholding. Valid water pixels were selected based on NDWI > 0.3. A total of 1847 valid pixel-by-pixel pairing samples were extracted from the 41 paired images.

[0041] The optical SSC (calibration coefficient) of each paired point was calculated using the Sentinel-2B red-green band ratio method. = 84.3, =1.27), with this as the dependent variable, and the VV polarization backscattering coefficient of SAR as the dependent variable. with local angle of incidence Using the least squares regression rate as the independent variable, the coefficients of the watershed-specific SAR–SSC inversion equation are determined: = 12.47, = 0.83, = 0.021; Coefficient of determination of the equation on the independent validation set (20% random sampling) = 0.79, RMSE = 37.2 mg / L, which meets the accuracy requirements for engineering applications and can be used for full-basin SSC estimation during cloud and rain cover periods.

[0042] Step S3: Delineate the sub-basin boundaries using the target watershed digital elevation model (DEM) and divide the target watershed into N sub-basin units. For all water body pixels within each sub-basin, calculate the optical suspended sediment concentration using the red-green band ratio method for effective optical pixels, and calculate the SAR suspended sediment concentration using the established SAR-SSC inversion equation for cloud-polluted pixels. Use the area-weighted average of all water body pixels within each sub-basin as the suspended sediment concentration feature value for that sub-basin node. Simultaneously, extract the Normalized Differential Vegetation Index (NDVI), Topographic Humidity Index (TWI), and the 24-hour cumulative rainfall from the nearest meteorological station for each sub-basin as auxiliary node features, which together constitute the node feature matrix.

[0043] In step S3, the terrain humidity index TWI is: ;in, The unit is the catchment area, in m² / m. The local slope is extracted from the DEM; the TWI of each sub-basin is the median of the TWI values ​​of all raster cells within that sub-basin.

[0044] For a sub-basin with a catchment area of ​​less than 1 km², if there are fewer than 5 effective water body pixels within it, the suspended sediment concentration characteristic value of the sub-basin is replaced by the area-weighted average of the suspended sediment concentration characteristic values ​​of its directly upstream adjacent sub-basins. The weight is the ratio of the confluence contribution area of ​​the corresponding upstream sub-basin to the sub-basin.

[0045] Using SRTM 30 m DEM as the topographic data source, the D8 single-flow direction algorithm was used to extract the confluence direction grid. The sub-basin boundaries were delineated with a catchment area threshold of 1 km². The Kuye River Basin was divided into N = 47 sub-basin units, with the catchment area of ​​each sub-basin ranging from 0.82 km² to 1024 km², and the median area being 68.5 km².

[0046] For remote sensing data from July 14, 2023, the calibrated red-green band ratio method was directly applied to the sub-basins of the effective optical area. = 84.3 × ( / )^1.27) Optical SSC was calculated; the SSC of the sub-basins in the cloud-polluted area was calculated using the SAR-SSC inversion equation; among them, the SSC feature values ​​of the four source sub-basins with a catchment area of ​​less than 1 km² were replaced by the area-weighted average of the SSC of the upstream adjacent sub-basins because there were less than 5 effective water body pixels. The mean NDVI and median TWI were also extracted for each sub-basin. The 24-hour cumulative rainfall from the nearest weather station, together with the SSC feature values, constitutes a 47×4 node feature matrix, which serves as the input to the spatiotemporal graph neural network.

[0047] Figure 2 The diagram illustrates the structure of the DEM and the directed hydrological topology; sub-basin nodes are distributed according to geographical location, and directed edges point from upstream to downstream, with edge thickness determined by normalized edge weights. The size of the plot shows that the weight of the main channel is significantly greater than that of the tributary branches, which intuitively reflects the spatial differentiation of the intensity of water and sediment transfer between upstream and downstream.

[0048] Step S4: Using each sub-basin as a graph node, directed edges are generated based on the water flow direction extracted from the DEM. The function value of the channel length and average slope between adjacent sub-basins is used as the weight of the directed edges to construct a directed hydrological topology graph representing the upstream and downstream water and sediment transport relationship. A spatiotemporal graph neural network is constructed, consisting of a two-layer directed graph convolution module and a single-head temporal attention module connected in series. The measured sediment yield of each surface hydrological station is used as the supervision signal. The training set, validation set, and test set are divided in a 7:1.5:1.5 ratio according to time sequence. The spatiotemporal graph neural network is trained using a loss function with the negative Nash-Sutcliffe efficiency coefficient (NSE) as the core.

[0049] In step S4, the directed hydrological topology map consists of upstream sub-basins. Pointing to downstream sub-basins Calculation of directed edge weights: ;in, The length of the river channel from the center point of the upstream sub-basin u to the center point of the downstream sub-basin v along the confluence path extracted by the DEM, in meters. The average gradient of the confluence path is expressed in m / m. For source sub-basin nodes without upstream adjacent sub-basins, their graph convolutional aggregation only uses their own node features and does not participate in cross-node message passing.

[0050] The spatial graph convolutional module of the spatiotemporal graph neural network adopts an in-degree normalization aggregation strategy. Layer nodes The feature update formula is: ;in, Let v be the set of upstream neighboring sub-basin nodes of node v in a directed hydrological topology graph. The directed edge weights are normalized to the interval [0,1]. and The first The trainable weight matrix and bias vector of layer graph convolution, where CONCAT represents the vector concatenation operation.

[0051] The temporal attention module takes a 7-day node feature time series as input, performs a single-layer linear transformation, and then uses a softmax function to calculate the attention weights at each time point. The weighted sum is then input into the fully connected layer, and the predicted sediment yield of each sub-basin is output.

[0052] The spatiotemporal graph neural network is trained using the following loss function: ;in, The predicted sediment yield values ​​for each sub-basin at each observation station at time t are: To correspond to the measured sediment yield, The time mean of measured sediment yield for the training set.

[0053] Training was performed using the Adam optimizer with an initial learning rate of 1×10⁻⁶. - ³, using a cosine annealing strategy to decay the learning rate to 1×10⁻⁶ over 200 training epochs. -5 After each graph convolution, a dropout rate of p=0.1 is applied to achieve both regularization and Monte Carlo uncertainty estimation during the inference stage.

[0054] For sub-basins lacking ground-based observation stations, training labels are constructed as follows: Measured sediment yield in downstream sub-basins with observation stations Based on this, the training label sediment yield of the station-free sub-basin is estimated using the following formula. : ;in, and These represent the catchment areas of the sub-basin without monitoring stations and the downstream sub-basin with monitoring stations, respectively, in km². and These are the surface erosion intensity indices for the corresponding sub-basins; The estimated sediment yield is only used for auxiliary supervision during the training phase. It is assigned a weight coefficient of 0.3 in the loss function, while the measured sediment yield of the sub-basin with observation stations has a corresponding weight coefficient of 1.0, in order to distinguish the confidence difference between the measured data and the estimated data.

[0055] The training set consisted of 1825 days of data from 2016 to 2021, the validation set consisted of 365 days of data from 2022, and the test set consisted of 365 days of data from 2023. The training labels for the remaining 44 sub-basins without monitoring stations were estimated from measured values ​​of downstream sub-basins with monitoring stations, based on the ratio of catchment area to surface erosion intensity index, with a weighting coefficient of 0.3.

[0056] The spatiotemporal graph neural network employs two layers of directed graph convolutions (64 hidden dimensions per layer) and single-head temporal attention (sequence length 7 days), with a total of approximately 52,000 parameters; it uses the Adam optimizer and an initial learning rate of 1×10⁻⁶. With a batch size of 16 and 200 training rounds, cosine annealing reduces the learning rate to 1×10⁻⁶. The convolutional layer dropout rate was p = 0.1, and negative NSE was used as the core loss function. Training took approximately 4.2 hours on a single NVIDIA V100 GPU.

[0057] On the 2023 test set, the mean NSE of sediment yield prediction for the three sub-basins with stations was 0.81, and the mean RMSE was 128.3 t / km², which is 0.19 higher than the NSE of the single-station interpolation benchmark method, verifying the effectiveness of the spatiotemporal graph neural network in integrating upstream information and temporal features.

[0058] Step S5: The trained spatiotemporal graph neural network is used to perform real-time sediment yield inversion for each sub-basin of the target watershed. During the inference phase, the Dropout layer is kept active, and 50 random forward propagations are performed. The mean of the 50 predicted values ​​is used as the estimated sediment yield, and the coefficient of variation of the 50 predicted values ​​is used as the uncertainty index. Sub-basins with a coefficient of variation exceeding 0.3 are automatically marked as high uncertainty areas and a supplementary sampling request is triggered. The node features of the corresponding sub-basin are updated with the latest SAR or optical imagery from the supplementary sampling, and then local re-inference is performed to output the updated estimated sediment yield and confidence interval.

[0059] The estimated sediment yield and uncertainty indicators are calculated as follows: Perform 50 random forward propagations while maintaining Dropout activation. The prediction value for sub-basin S after the i-th propagation is denoted as... Forecast mean With the predicted standard deviation calculate: ; ;by Output as the 95% approximate confidence interval for sediment yield; coefficient of variation The calculation formula is: ,right If a sub-basin exceeds 0.3, triggering a supplementary imagery request, the latest first satellite imagery covering the sub-basin and within 3 days of the trigger time will be selected first. If no new imagery is acquired within 72 hours, the original predicted mean will be retained and marked "Insufficient data, doubtful estimation" in the results.

[0060] For the heavy rainfall event from July 14 to 16, 2023, a trained spatiotemporal graph neural network was used to perform 50 random forward propagations daily (with Dropout kept active) on 47 sub-basins. The initial inference results on July 14 showed that the estimated total daily sediment yield for the entire basin was 1.23 × 10⁻⁶. t, the 95% confidence interval is [0.98 × 10 t, 1.48×10 t]; where the coefficient of variation (CV) of the five upstream source sub-basins is (= / If the CV exceeds 0.3 (maximum CV = 0.47), the system will automatically trigger a supplementary acquisition request, prioritizing the latest Sentinel-1A transit imagery with a time difference of no more than 3 days.

[0061] After acquiring new imagery at 09:23 on July 15th, local re-inference was completed, and the CV of the five high-uncertainty sub-basins all decreased to below 0.25; the estimated sediment yield of the corresponding sub-basin at Wenjiachuan station was updated to 4.87 × 10⁻⁶. t, compared with the measured sediment yield of 4.65 × 10 at Wenjiachuan Station on that day. The relative error of t is only 4.7%, indicating that the supplementary sampling and re-inference mechanism effectively suppresses the uncertainty caused by rapid changes in ground conditions during the rainstorm period.

[0062] The method also includes an online rolling update step, specifically: using a 30-day rolling time window, the root mean square error (RMSE) of the spatiotemporal graph neural network for predicting sediment yield in sub-basins with observation stations is statistically analyzed; when the RMSE exceeds 20% of the average measured sediment yield within the same window, a model update is triggered.

[0063] During model updates, the training set is the observation data from the 365 days prior to the current time. All weight parameters of the first layer graph convolution of the spatiotemporal graph neural network are frozen, and only the parameters of the second layer graph convolution and the temporal attention module are fine-tuned. The fine-tuning uses an initial learning rate of 1×10⁻⁶. -4 The Adam optimizer was trained for 50 epochs.

[0064] After fine-tuning, if the updated model's NSE on the validation set is not lower than that of the original model, the updated model will replace the original model and be put into use; otherwise, the original model will be retained and an anomaly log will be recorded for manual review.

[0065] Taking the 30-day rolling window in August 2023 as an example, the RMSE of the spatiotemporal graph neural network for predicting sediment yield at the three stations were 102.3 t / km², 87.6 t / km², and 115.4 t / km², respectively. The corresponding measured average sediment yields within the window were 430.2 t / km², 305.7 t / km², and 512.1 t / km², respectively. The ratios of RMSE to the average were 23.8%, 28.7%, and 22.5%, respectively. All three stations exceeded the 20% trigger threshold, and the online rolling update of the model was initiated on August 31, 2023.

[0066] Using 365 days of observation data from September 2022 to August 2023 as the training set, all weights of the first-layer graph convolution were frozen, and only the parameters of the second-layer graph convolution and the temporal attention module were fine-tuned. The Adam optimizer was used (initial learning rate 1×10⁻⁶). The training process took 50 epochs and approximately 23 minutes. After fine-tuning, the NSE of the updated model on the September 2023 validation set improved from 0.78 to 0.83, meeting the admission criteria of no less than that of the model before the update. The new model then replaced the original model and was put into real-time use, achieving adaptive online iteration without human intervention.

[0067] The specific embodiments described above further illustrate the purpose, technical solution, and beneficial effects of the present invention. It should be understood that the above description is only a specific embodiment of the present invention and is not intended to limit the scope of protection of the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. A method for dynamic monitoring and assessment of tributary sediment yield based on multi-source remote sensing data, characterized in that, The method includes the following steps: Step S1: Acquire the first satellite C-band VV polarimetric synthetic aperture radar SAR image and the second satellite optical multispectral image covering the same time phase of the target watershed, as well as the measured flow and suspended sediment concentration SSC observation sequences of various surface hydrological stations in the target watershed; using the scene classification layer SCL band of the second satellite, mark pixels with SCL values ​​of 8, 9 or 10 as cloud pollution pixels, and mark the remaining pixels falling into the water body range as effective optical pixels; Step S2: Using the first satellite SAR image and the second satellite optical image acquired during the same period of historical clear sky in the target watershed as training samples; the suspended sediment concentration (SSC) in the water body is calculated using the red-green band ratio method for the effective optical pixels; the VV polarization backscattering coefficient of the SAR image at the same spatial location is paired with the suspended sediment concentration (SSC) pixel by pixel; and the SAR-SSC inversion equation for the target watershed is established by least squares regression using no less than 30 pairs of effective observation data. Step S3: Delineate the sub-basin boundaries using the target watershed digital elevation model (DEM) and divide the target watershed into N sub-basin units. For all water body pixels within each sub-basin, calculate the optical suspended sediment concentration using the red-green band ratio method for effective optical pixels, and calculate the SAR suspended sediment concentration using the established SAR-SSC inversion equation for cloud-polluted pixels. Use the area-weighted average of all water body pixels within each sub-basin as the suspended sediment concentration feature value of the sub-basin node. Simultaneously, extract the Normalized Differential Vegetation Index (NDVI), Topographic Humidity Index (TWI), and the 24-hour cumulative rainfall from the nearest meteorological station for each sub-basin as auxiliary node features, which together constitute the node feature matrix. Step S4: Using each sub-basin as a graph node, directed edges are generated based on the water flow direction extracted from the DEM. The function value of the channel length and average slope between adjacent sub-basins is used as the weight of the directed edges to construct a directed hydrological topology graph representing the upstream and downstream water and sediment transfer relationship. A spatiotemporal graph neural network is constructed, consisting of a two-layer directed graph convolution module and a single-head temporal attention module connected in series. The measured sediment yield of each surface hydrological station is used as the supervision signal. The training set, validation set and test set are divided in a 7:1.5:1.5 ratio according to the time sequence. The spatiotemporal graph neural network is trained using a loss function with the negative Nash-Sutcliffe efficiency coefficient (NSE) as the core. Step S5: The trained spatiotemporal graph neural network is used to perform real-time sediment yield inversion for each sub-basin of the target watershed. During the inference phase, the Dropout layer is kept active, and 50 random forward propagations are performed. The mean of the 50 predicted values ​​is used as the estimated sediment yield, and the coefficient of variation of the 50 predicted values ​​is used as the uncertainty index. Sub-basins with a coefficient of variation exceeding 0.3 are automatically marked as high uncertainty areas and a supplementary sampling request is triggered. The node features of the corresponding sub-basin are updated with the latest SAR or optical imagery from the supplementary sampling, and then local re-inference is performed to output the updated estimated sediment yield and confidence interval.

2. The method according to claim 1, characterized in that, Optical suspended sediment concentration The calculation formula is: ;in, The surface reflectance in the red band of the second satellite. For the surface reflectance in the green band of the second satellite, and The watershed-specific coefficients were determined by log-linear regression using ground-measured suspended sediment concentration data. The SAR–SSC inversion equation is in the form of: ;in, The VV polarization backscattering coefficient is... The local angle of incidence at the water surface pixel, in degrees. The least squares regression coefficient is used; the selection criterion for water body pixels is that the normalized difference water index (NDWI) is greater than 0.

3.

3. The method according to claim 2, characterized in that, In step S2, before establishing the SAR–SSC inversion equation, the following preprocessing is performed on the first satellite SAR image in sequence: Using the calibration vector in the metadata of the first satellite product, the original digital quantization (DN) value is converted into sigma-naught backscattering coefficients; Using a 30 m resolution as a reference, a range-Doppler terrain correction algorithm is used to eliminate geometric distortions caused by terrain undulations and orthorectify the SAR image to the same UTM projection coordinate system as the second satellite image. A 5×5 window Lee filter was applied to the radiometrically calibrated backscattering coefficient image to suppress speckle noise; The filtered backscattering coefficient image histogram is segmented using Otsu thresholding. Pixels below the threshold are identified as water surface pixels, generating a water surface mask to limit the spatial range of subsequent SAR-SSC paired samples.

4. The method according to claim 3, characterized in that, In step S3, the terrain humidity index TWI is: ;in, The unit is the catchment area, in m² / m. The local slope is extracted from the DEM; the TWI of each sub-basin is the median of the TWI values ​​of all raster cells within that sub-basin; For a sub-basin with a catchment area of ​​less than 1 km², if there are fewer than 5 effective water body pixels within it, the suspended sediment concentration characteristic value of the sub-basin is replaced by the area-weighted average of the suspended sediment concentration characteristic values ​​of its directly upstream adjacent sub-basins. The weight is the ratio of the confluence contribution area of ​​the corresponding upstream sub-basin to the sub-basin.

5. The method according to claim 4, characterized in that, In step S4, the directed hydrological topology map consists of upstream sub-basins. Pointing to downstream sub-basins Calculation of directed edge weights: ;in, The length of the river channel from the center point of the upstream sub-basin u to the center point of the downstream sub-basin v along the confluence path extracted by the DEM, in meters. The average gradient of the confluence path is expressed in m / m. For source sub-basin nodes without upstream adjacent sub-basins, their graph convolutional aggregation only uses their own node features and does not participate in cross-node message passing.

6. The method according to claim 5, characterized in that, In step S4, the spatial graph convolutional module of the spatiotemporal graph neural network adopts an in-degree normalization aggregation strategy. Layer nodes The feature update formula is: ;in, Let v be the set of upstream neighboring sub-basin nodes of node v in a directed hydrological topology graph. The directed edge weights are normalized to the interval [0,1]. and The first The trainable weight matrix and bias vector of layer graph convolution, where CONCAT represents the vector concatenation operation; The temporal attention module takes a 7-day node feature time series as input, performs a single-layer linear transformation, and then uses a softmax function to calculate the attention weights at each time point. The weighted sum is then input into the fully connected layer, and the predicted sediment yield of each sub-basin is output.

7. The method according to claim 6, characterized in that, In step S4, the spatiotemporal graph neural network is trained using the following loss function: ;in, The predicted sediment yield values ​​for each sub-basin at each observation station at time t are: To correspond to the measured sediment yield, The training set is the time-averaged measured sediment yield; Training was performed using the Adam optimizer with an initial learning rate of 1×10⁻⁶. - ³, using a cosine annealing strategy to decay the learning rate to 1×10⁻⁶ over 200 training epochs. -5 After each graph convolution, a dropout rate of p=0.1 is applied to achieve both regularization and Monte Carlo uncertainty estimation during the inference stage.

8. The method according to claim 7, characterized in that, In step S5, the estimated sediment yield and uncertainty index are calculated as follows: Perform 50 random forward propagations while maintaining Dropout activation, the first... The predicted value of the second propagation for the sub-basin S is denoted as... Forecast mean With the predicted standard deviation calculate: ; ;by Output as the 95% approximate confidence interval for sediment yield; coefficient of variation The calculation formula is: ,right If a sub-basin triggers a supplementary imagery request with a value exceeding 0.3, the latest first satellite imagery covering that sub-basin and with a time difference of no more than 3 days from the trigger time will be selected first. If no new images are acquired within 72 hours, the original predicted mean is retained and marked "Insufficient data, estimation questionable" in the results.

9. The method according to claim 8, characterized in that, In step S4, for sub-basins lacking ground observation stations, training labels are constructed as follows: Measured sediment yield in downstream sub-basins with observation stations Based on this, the training label sediment yield of the station-free sub-basin is estimated using the following formula. : ;in, and These represent the catchment areas of the sub-basin without monitoring stations and the downstream sub-basin with monitoring stations, respectively, in km². and These are the surface erosion intensity indices for the corresponding sub-basins; The estimated sediment yield is only used for auxiliary supervision during the training phase. It is assigned a weight coefficient of 0.3 in the loss function, while the measured sediment yield of the sub-basin with observation stations has a corresponding weight coefficient of 1.0, in order to distinguish the confidence difference between the measured data and the estimated data.

10. The method according to claim 9, characterized in that, The method also includes an online rolling update step for the model, specifically: using a 30-day rolling time window, the root mean square error (RMSE) of the spatiotemporal graph neural network for predicting sediment yield in sub-basins with observation stations is statistically analyzed; when the RMSE exceeds 20% of the average measured sediment yield within the same window, the model is updated. During model updates, the training set is the observation data from the 365 days prior to the current time. All weight parameters of the first layer graph convolution of the spatiotemporal graph neural network are frozen, and only the parameters of the second layer graph convolution and the temporal attention module are fine-tuned. The fine-tuning uses an initial learning rate of 1×10⁻⁶. -4 The Adam optimizer was trained for 50 epochs; After fine-tuning, if the updated model's NSE on the validation set is not lower than that of the original model, the updated model will replace the original model and be put into use; otherwise, the original model will be retained and an anomaly log will be recorded for manual review.

Images